[Pomp-commits] r401 - in pkg: src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Oct 24 17:14:06 CEST 2010
Author: kingaa
Date: 2010-10-24 17:14:06 +0200 (Sun, 24 Oct 2010)
New Revision: 401
Modified:
pkg/src/pomp_mat.h
pkg/src/probe_marginal.c
pkg/src/probe_nlar.c
pkg/src/synth_lik.c
pkg/tests/ou2-probe.Rout.save
pkg/tests/ricker-probe.Rout.save
Log:
- get rid of all pivoting in synthetic likelihood calculation
- make the inline routines in 'pomp_mat.h' a bit more flexible
Modified: pkg/src/pomp_mat.h
===================================================================
--- pkg/src/pomp_mat.h 2010-10-23 22:02:50 UTC (rev 400)
+++ pkg/src/pomp_mat.h 2010-10-24 15:14:06 UTC (rev 401)
@@ -4,10 +4,17 @@
#include "pomp_internal.h"
#include <R_ext/Lapack.h>
-static R_INLINE void pomp_backsolve (double *a, int m, int n, double *x, int incx,
+static R_INLINE void pomp_backsolve (double *a, int lda, int n, double *x, int incx,
char *uplo, char *transpose, char *unit) {
// Level 2 BLAS triangular-matrix solver DTRSV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
- F77_NAME(dtrsv)(uplo,transpose,unit,&n,a,&m,x,&incx);
+ // DTRSV: x <- A ^{-1} x -or- x <- A ^{-T} x, A triangular
+ // N is the order of A
+ // LDA is A's leading dimension
+ // INCX is the increment between successive X locations
+ // UPLO is "U" or "L" depending on whether A is upper or lower triangular
+ // TRANSPOSE is "T" or "N" depending on whether the transpose is desired
+ // DIAG is "U" or "N" depending on whether A is unit triangular or not
+ F77_NAME(dtrsv)(uplo,transpose,unit,&n,a,&lda,x,&incx);
}
static R_INLINE void pomp_qr (double *a, int m, int n, int *pivot, double *tau) {
@@ -25,22 +32,17 @@
for (j = 0; j < n; j++) pivot[j]--;
}
-static R_INLINE void pomp_qrqy (double *c, double *a, double *tau, int m, int n, int k, int left, int tp) {
- int lda, info, lwork = -1;
- char *side, *trans;
+static R_INLINE void pomp_qrqy (double *c, double *a, int lda, double *tau, int m, int n, int k, char *side, char *transpose) {
+ int info, lwork = -1;
double work1;
- side = (left) ? "left" : "right";
- lda = (left) ? m : n;
- trans = (tp) ? "transpose" : "no transpose";
-
// workspace query
// DORMQR(SIDE,TRANS,M,N,K,A,LDA,TAU,C,LDC,WORK,LWORK,INFO)
- F77_NAME(dormqr)(side,trans,&m,&n,&k,a,&lda,tau,c,&m,&work1,&lwork,&info);
+ F77_NAME(dormqr)(side,transpose,&m,&n,&k,a,&lda,tau,c,&m,&work1,&lwork,&info);
lwork = (int) ceil(work1);
{ // actual call
double work[lwork];
- F77_NAME(dormqr)(side,trans,&m,&n,&k,a,&lda,tau,c,&m,work,&lwork,&info);
+ F77_NAME(dormqr)(side,transpose,&m,&n,&k,a,&lda,tau,c,&m,work,&lwork,&info);
}
}
Modified: pkg/src/probe_marginal.c
===================================================================
--- pkg/src/probe_marginal.c 2010-10-23 22:02:50 UTC (rev 400)
+++ pkg/src/probe_marginal.c 2010-10-24 15:14:06 UTC (rev 401)
@@ -135,7 +135,7 @@
R_qsort(x,1,nx);
// solve R b = Q'x for b
- pomp_qrqy(x,mm,tau,nx,1,np,1,1); // y <- Q'y
+ pomp_qrqy(x,mm,nx,tau,nx,1,np,"left","transpose"); // y <- Q'y
pomp_backsolve(mm,nx,np,x,1,"Upper","No transpose","Non-unit"); // y <- R^{-1} y
// unpivot and store the coefficients in beta
Modified: pkg/src/probe_nlar.c
===================================================================
--- pkg/src/probe_nlar.c 2010-10-23 22:02:50 UTC (rev 400)
+++ pkg/src/probe_nlar.c 2010-10-24 15:14:06 UTC (rev 401)
@@ -120,7 +120,7 @@
// first QR decompose the model matrix
pomp_qr(X,nx,nterms,pivot,tau);
// then solve R b = Q'y for b
- pomp_qrqy(yp,X,tau,nx,1,nterms,1,1); // y <- Q'y
+ pomp_qrqy(yp,X,nx,tau,nx,1,nterms,"left","transpose"); // y <- Q'y
pomp_backsolve(X,nx,nterms,yp,1,"Upper","No transpose","Non-unit"); // y <- R^{-1} y
// unpivot and store coefficients in beta
Modified: pkg/src/synth_lik.c
===================================================================
--- pkg/src/synth_lik.c 2010-10-23 22:02:50 UTC (rev 400)
+++ pkg/src/synth_lik.c 2010-10-24 15:14:06 UTC (rev 401)
@@ -17,7 +17,6 @@
int ncol = dim[1];
double alpha = 2.0, beta = 1.25;
double w[nrow], tau[ncol], work[ncol];
- int pivot[ncol];
int info = 0, ione = 1;
double one = 1.0;
double *y1, *y2, *yp;
@@ -32,22 +31,20 @@
// compute column means, center each column, precondition
memcpy(y1,y,nrow*ncol*sizeof(double));
for (yp = y1, j = 0; j < ncol; j++, yp += nrow) {
- for (xx = 0, i = 0; i < nrow; i++) xx += yp[i];
- xx /= nrow;
- for (i = 0; i < nrow; i++) yp[i] -= xx; // center the column
- for (xx = 0, i = 0; i < nrow; i++) xx += yp[i]*yp[i];
- x = sqrt(xx/(nrow-1)); // column SD
- for (i = 0; i < nrow; i++) yp[i] /= x; // precondition
+ for (x = 0, i = 0; i < nrow; i++) x += yp[i];
+ x /= nrow;
+ for (i = 0; i < nrow; i++) yp[i] -= x; // center the column
+ for (x = 0, i = 0; i < nrow; i++) x += yp[i]*yp[i];
+ d = sqrt(x/(nrow-1)); // column SD
+ for (i = 0; i < nrow; i++) yp[i] /= d; // precondition
}
// do first QR decomposition & backsolve
memcpy(y2,y1,nrow*ncol*sizeof(double));
- {
- // LAPACK QR decomposition without pivoting DGEQR2(M,N,A,LDA,TAU,WORK,INFO)
- F77_NAME(dgeqr2)(&nrow,&ncol,y2,&nrow,tau,work,&info);
- // Level-3 BLAS triangular matrix solver DTRSM(SIDE,UPLO,TRANS,DIAG,M,N,ALPHA,A,LDA,B,LDB)
- F77_NAME(dtrsm)("Right","Upper","No transpose","Non-unit",&nrow,&ncol,&one,y2,&nrow,y1,&nrow);
- }
+ // LAPACK QR decomposition without pivoting DGEQR2(M,N,A,LDA,TAU,WORK,INFO)
+ F77_NAME(dgeqr2)(&nrow,&ncol,y2,&nrow,tau,work,&info);
+ // Level-3 BLAS triangular matrix solver DTRSM(SIDE,UPLO,TRANS,DIAG,M,N,ALPHA,A,LDA,B,LDB)
+ F77_NAME(dtrsm)("right","upper","no transpose","non-unit",&nrow,&ncol,&one,y2,&nrow,y1,&nrow);
// create Campbell weight vector
d0 = sqrt(ncol)+alpha/sqrt(2.0);
@@ -71,9 +68,9 @@
// compute weighted column means, center each column, precondition
memcpy(y1,y,nrow*ncol*sizeof(double));
for (yp = y1, j = 0; j < ncol; j++, yp += nrow) {
- for (xx = 0, x = 0, i = 0; i < nrow; i++) {
- xx += w[i]*yp[i];
+ for (x = 0, xx = 0, i = 0; i < nrow; i++) {
x += w[i];
+ xx += w[i]*yp[i];
}
xx /= x; // column mean
for (i = 0; i < nrow; i++) yp[i] -= xx; // center the column
@@ -89,13 +86,13 @@
}
// do second QR decomposition & backsolve
- pomp_qr(y1,nrow,ncol,pivot,tau); // Q*R = diag(w)*Y*inv(diag(d))
- for (j = 0; j < ncol; j++) tau[j] = ydat[pivot[j]]; // unpivot
- pomp_backsolve(y1,nrow,ncol,tau,1,"Upper","Transpose","Non-unit");
+ // LAPACK QR decomposition without pivoting DGEQR2(M,N,A,LDA,TAU,WORK,INFO)
+ F77_NAME(dgeqr2)(&nrow,&ncol,y1,&nrow,tau,work,&info);
+ pomp_backsolve(y1,nrow,ncol,ydat,1,"upper","transpose","non-unit");
// compute residual sum of squares and add up logs of diag(R)
for (yp = y1, rss = 0, i = nrow+1, j = 0; j < ncol; j++, yp += i) { // yp marches along the diagonal of R
- x = tau[j];
+ x = ydat[j];
rss += x*x;
half_log_det += log(fabs(*yp)); // log(diag(R))
}
Modified: pkg/tests/ou2-probe.Rout.save
===================================================================
--- pkg/tests/ou2-probe.Rout.save 2010-10-23 22:02:50 UTC (rev 400)
+++ pkg/tests/ou2-probe.Rout.save 2010-10-24 15:14:06 UTC (rev 401)
@@ -219,7 +219,7 @@
0.3016983 0.2957043 0.3216783 0.3016983 0.2957043 0.3216783
$synth.loglik
-[1] 92.59923
+[1] 86.8972
>
> pb <- probe(
Modified: pkg/tests/ricker-probe.Rout.save
===================================================================
--- pkg/tests/ricker-probe.Rout.save 2010-10-23 22:02:50 UTC (rev 400)
+++ pkg/tests/ricker-probe.Rout.save 2010-10-24 15:14:06 UTC (rev 401)
@@ -266,7 +266,7 @@
+ )
+ )
user system elapsed
- 12.803 0.030 12.836
+ 12.860 0.031 12.893
> plot(pm)
>
> cbind(truth=coef(ricker),est=coef(pm),guess=coef(po))
@@ -325,7 +325,7 @@
+ )
+ )
user system elapsed
- 9.545 0.013 9.559
+ 9.608 0.009 9.618
> plot(pm)
> plot(as(pm,"pomp"),variables="y")
> plot(simulate(pm),variables="y")
@@ -484,7 +484,7 @@
0.8171828 0.8171828 0.1358641 0.6313686 0.2117882 0.2677323
$synth.loglik
-[1] 29.79205
+[1] 29.83407
> plot(pb)
>
More information about the pomp-commits
mailing list