fix lossdistrib_truncated

2 files changed, 23 insertions, 23 deletions

diff --git a/R/distrib.R b/R/distrib.R
index 8b10394..277b324 100644
--- a/R/distrib.R
+++ b/R/distrib.R
@@ -294,10 +294,12 @@ lossdistCZ <- function(p, w, S, N, defaultflag=FALSE, rho, Z){
        q = matrix(0, N, length(Z)))$q

 }

 

-lossdistC.truncated <- function(p, w, S, N, T=N){

-    ## C version of lossdistrib2, roughly 50 times faster

+lossdistC.truncated <- function(p, w, S, N, T=N, defaultflag=FALSE){

+    ## truncated version of lossdistrib

+    ## q[i] is 0 for i>=T

     .C("lossdistrib_truncated", as.double(p), as.integer(length(p)),

-       as.double(w), as.double(S), as.integer(N), as.integer(T), q = double(T))$q

+       as.double(w), as.double(S), as.integer(N), as.integer(T), as.logical(defaultflag),

+       q = double(N))$q

 }

 

 recovdistC <- function(dp, pp, w, S, N){

diff --git a/src/lossdistrib.c b/src/lossdistrib.c
index 61f3c84..d28b1b2 100644
--- a/src/lossdistrib.c
+++ b/src/lossdistrib.c
@@ -145,7 +145,7 @@ void lossdistrib_Z(double *p, int *np, double *w, double *S, int *N, int *defaul
     int i, j;
     double* pshocked = malloc(sizeof(double) * (*np) * (*nZ));
 
-#pragma omp parallel for private(j)
+    #pragma omp parallel for private(j)
     for(i = 0; i < *nZ; i++){
         for(j = 0; j < *np; j++){
             pshocked[j + (*np) * i] = shockprob(p[j], rho[j], Z[i], 0);
@@ -172,34 +172,32 @@ void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N,
 
     int i, j, d1, d2, M;
     double lu, d, p1, p2;
-    double *q1, *q2;
+    double *qtemp;
+    int one = 1;
+    double pbar;
+    int bound;
 
     lu = 1./(*N-1);
-    q1 = calloc( *T, sizeof(double));
-    q2 = calloc( *T, sizeof(double));
+    qtemp = malloc( *T * sizeof(double));
     q[0] = 1;
     M = 1;
     for(i=0; i<(*np); i++){
         d = (*defaultflag)? w[i] / lu : S[i] * w[i] / lu;
-        d1 = floor(d);
-        d2 = ceil(d);
-        p1 = p[i] * (d2-d);
+        d1 = (int)floor(d);
+        d2 = d1 + 1;
+        p1 = p[i] * ((double)d2-d);
         p2 = p[i] - p1;
-        for(j=0; j < MIN(M, *T); j++){
-            q1[j] = p1 * q[j];
-            q2[j] = p2 * q[j];
-            q[j] = (1-p[i]) * q[j];
-        }
-        for(j=0; j < MIN(M, *T-d1); j++){
-            q[d1+j] += q1[j];
-        };
-        for(j=0; j < MIN(M, *T-d2); j++){
-            q[d2+j] += q2[j];
-        };
+        memcpy(qtemp, q, MIN(M, *T) * sizeof(double));
+        pbar = 1-p[i];
+        bound = MIN(M, *T);
+        dscal_(&bound, &pbar, q, &one);
+        bound = MIN(M, *T-d1);
+        daxpy_(&bound, &p1, qtemp, &one, q+d1, &one);
+        bound = MIN(M, *T-d2);
+        daxpy_(&bound, &p2, qtemp, &one, q+d2, &one);
         M += d2;
     }
-    free(q1);
-    free(q2);
+    free(qtemp);
 }
 
 void lossdistrib_joint(double *p, int *np, double *w, double *S, int *N, int *defaultflag, double *q) {