fix BCloss_recov_trunc

1 files changed, 60 insertions, 48 deletions

diff --git a/src/lossdistrib.c b/src/lossdistrib.c
index a3b3866..412437d 100644
--- a/src/lossdistrib.c
+++ b/src/lossdistrib.c
@@ -527,9 +527,9 @@ void BCloss_recov_dist(const double *defaultprob, const int *dim1, const int *di
         #pragma omp threadprivate(gshocked, Sshocked, Rshocked)
         #pragma omp parallel
         {
-            gshocked = malloc((*dim1) * sizeof(double));
-            Sshocked = malloc((*dim1) * sizeof(double));
-            Rshocked = malloc((*dim1) * sizeof(double));
+            gshocked = malloc(3 * (*dim1) * sizeof(double));
+            Sshocked = gshocked + *dim1;
+            Rshocked = Sshocked + *dim1;
         }
         #pragma omp parallel for
         for(int i = 0; i < *n; i++) {
@@ -547,8 +547,6 @@ void BCloss_recov_dist(const double *defaultprob, const int *dim1, const int *di
         #pragma omp parallel
         {
             free(gshocked);
-            free(Sshocked);
-            free(Rshocked);
         }
         cblas_dgemv(CblasColMajor, CblasNoTrans, *N, *n, alpha, Lw, *N, w, 1, beta, L + t * (*N), 1);
         cblas_dgemv(CblasColMajor, CblasNoTrans, *N, *n, alpha, Rw, *N, w, 1, beta, R + t * (*N), 1);
@@ -586,63 +584,68 @@ void BCloss_recov_trunc(const double *defaultprob, const int *dim1, const int *d
     double lu = 1./(*N+1);
     double* valL = malloc( T * sizeof(double));
     posK(T, *K, lu, valL);
+
     double* EL = malloc( (*n) * sizeof(double));
     double* ER = malloc( (*n) * sizeof(double));
     double* Lw = malloc(T * (*n) * sizeof(double));
     const double alpha = 1;
     const double beta = 0;
 
-    for(int t=0; t < (*dim2); t++) {
+    for(int t = 0; t < (*dim2); t++) {
         memset(Lw, 0, T * (*n) * sizeof(double));
+        static double* gshocked;
+        static double* gshockedbar;
+        static double* Sshocked;
+        static double* Rshocked;
+
+        #pragma omp threadprivate(gshocked, gshockedbar,Sshocked, Rshocked)
         #pragma omp parallel
         {
-            double g, Ktilde;
-            int Ttilde;
-            #pragma omp for
-            for(int i=0; i < *n; i++){
-                double* Rw = NULL;
-                double* Rshocked = malloc((*dim1) * sizeof(double));
-                double* Sshocked = malloc((*dim1) * sizeof(double));
-                double* gshocked = malloc((*dim1) * sizeof(double));
-                double* gshockedbar = malloc((*dim1) * sizeof(double));
-                double* valR = NULL;
-                for(int j=0; j < (*dim1); j++){
-                    g = defaultprob[j + (*dim1) * t];
-                    gshocked[j] = shockprob(g, rho[j], Z[i], 0);
-                    Sshocked[j] = shockseverity(1-recov[j], Z[i], rho[j], g);
-                    gshockedbar[j] = 1 - gshocked[j];
-                    Rshocked[j] = 1 - Sshocked[j];
-                }
+            gshocked = malloc(4 * (*dim1) * sizeof(double));
+            gshockedbar = gshocked + *dim1;
+            Rshocked = gshockedbar + *dim1;
+            Sshocked = Rshocked + *dim1;
+        }
+        #pragma omp for
+        for(int i = 0; i < *n; i++) {
 
-                lossdistrib(gshocked, dim1, issuerweights, Sshocked,
-                            N, &T, defaultflag, Lw + i * T);
-                ER[i] = 0;
-                Ktilde = *K - cblas_ddot(*dim1, issuerweights, 1, Sshocked, 1);
-                if(Ktilde > 0){
-                    Ttilde = (int) floor(Ktilde * (*N))+1;
-                    Rw = calloc(Ttilde, sizeof(double));
-                    lossdistrib(gshockedbar, dim1, issuerweights, Rshocked,
-                                N, &Ttilde, defaultflag, Rw);
-                    valR = malloc(Ttilde * sizeof(double));
-                    posK(Ttilde, Ktilde, lu, valR);
-                    ER[i] = cblas_ddot(Ttilde, Rw, 1, valR, 1);
-                }
-                if(Rw) {
-                    free(Rw);
-                }
-                if(valR) {
-                    free(valR);
-                }
-                free(Rshocked);
-                free(Sshocked);
-                free(gshocked);
-                free(gshockedbar);
+            double g;
+            for(int j = 0; j < (*dim1); j++){
+                g = defaultprob[j + (*dim1) * t];
+                gshocked[j] = shockprob(g, rho[j], Z[i], 0);
+                Sshocked[j] = shockseverity(1-recov[j], Z[i], rho[j], g);
+                gshockedbar[j] = 1 - gshocked[j];
+                Rshocked[j] = 1 - Sshocked[j];
+            }
+
+            lossdistrib(gshocked, dim1, issuerweights, Sshocked,
+                        N, &T, defaultflag, Lw + i * T);
+            ER[i] = 0;
+            double Ktilde = *K - cblas_ddot(*dim1, issuerweights, 1, Sshocked, 1);
+
+            if(Ktilde > 0) {
+                int Ttilde = (int) floor(Ktilde * (*N))+1;
+
+                double* Rw = calloc(2 * Ttilde, sizeof(double));
+                double* valR = Rw + Ttilde;
+                lossdistrib(gshockedbar, dim1, issuerweights, Rshocked,
+                            N, &Ttilde, defaultflag, Rw);
+
+                posK(Ttilde, Ktilde, lu, valR);
+                ER[i] = cblas_ddot(Ttilde, Rw, 1, valR, 1);
+                free(Rw);
             }
         }
+        #pragma omp parallel
+        {
+            free(gshocked);
+        }
         cblas_dgemv(CblasColMajor, CblasTrans, T, *n, alpha, Lw, T, valL, 1, beta, EL, 1);
+
         ELt[t] = cblas_ddot(*n, EL, 1, w, 1);
         ERt[t] = cblas_ddot(*n, ER, 1, w, 1);
     }
+
     free(Lw);
     free(EL);
     free(ER);
@@ -663,7 +666,7 @@ int main() {
     gsl_rng* gen = gsl_rng_alloc(gsl_rng_mt19937);
     for(int i = 0; i<n; i++) {
         weights[i] = 1./n;
-        rho[i] = 0.45;
+        rho[i] = 0.6;
         p[i] = gsl_rng_uniform(gen);
     }
     double* S = malloc(n * nZ * sizeof(double));
@@ -690,7 +693,16 @@ int main() {
     double* L = calloc(N * T, sizeof(double));
     double* R = calloc(N * T, sizeof(double));
 
-    BCloss_recov_dist(p, &n, &T, weights, S, Z, w, &nZ, rho, &N, &flag, L, R);
+    double K = 0.5;
+    //BCloss_recov_dist(p, &n, &T, weights, S, Z, w, &nZ, rho, &N, &flag, L, R);
+    double* ELt = malloc(sizeof(double) * T);
+    double* ERt = malloc(sizeof(double) * T);
+    BCloss_recov_trunc(p, &n, &T, weights, S, Z, w, &nZ, rho, &N, &K, &flag, ELt, ERt);
+    for(int t = 0; t < T; t++) {
+        printf("%f %f\n", ELt[t], ERt[t]);
+    }
+    free(ELt);
+    free(ERt);
     free(p);
     free(S);
     free(weights);