diff options
| -rw-r--r-- | R/distrib.R | 8 | ||||
| -rw-r--r-- | src/lossdistrib.c | 38 |
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) { |