[Pomp-commits] r687 - pkg/pomp/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Apr 29 16:24:53 CEST 2012
Author: kingaa
Date: 2012-04-29 16:24:52 +0200 (Sun, 29 Apr 2012)
New Revision: 687
Modified:
pkg/pomp/src/dmeasure.c
pkg/pomp/src/skeleton.c
Log:
- more work on the guts
Modified: pkg/pomp/src/dmeasure.c
===================================================================
--- pkg/pomp/src/dmeasure.c 2012-04-28 18:01:22 UTC (rev 686)
+++ pkg/pomp/src/dmeasure.c 2012-04-29 14:24:52 UTC (rev 687)
@@ -11,7 +11,6 @@
SEXP do_dmeasure (SEXP object, SEXP y, SEXP x, SEXP times, SEXP params, SEXP log, SEXP fun)
{
int nprotect = 0;
- int first = 1;
int mode = -1;
int give_log;
int ntimes, nvars, npars, ncovars, nreps, nrepsx, nrepsp, nobs;
@@ -20,12 +19,8 @@
SEXP tvec, xvec, yvec, pvec, cvec;
SEXP fn, fcall, rho, ans;
SEXP F;
- double *ts, *xs, *ps, *ys, *fs;
- double *cp, *xp, *pp, *tp, *yp;
- double *ft;
- int *dim, ndim[2];
int *sidx, *pidx, *cidx, *oidx;
- int i, j, k;
+ int *dim;
struct lookup_table covariate_table;
pomp_measure_model_density *ff = NULL;
@@ -37,7 +32,6 @@
PROTECT(y = as_matrix(y)); nprotect++;
dim = INTEGER(GET_DIM(y));
nobs = dim[0];
- ys = REAL(y);
if (ntimes != dim[1])
error("dmeasure error: length of 'times' and 2nd dimension of 'y' do not agree");
@@ -45,7 +39,6 @@
PROTECT(x = as_state_array(x)); nprotect++;
dim = INTEGER(GET_DIM(x));
nvars = dim[0]; nrepsx = dim[1];
- xs = REAL(x);
if (ntimes != dim[2])
error("dmeasure error: length of 'times' and 3rd dimension of 'x' do not agree");
@@ -53,7 +46,6 @@
PROTECT(params = as_matrix(params)); nprotect++;
dim = INTEGER(GET_DIM(params));
npars = dim[0]; nrepsp = dim[1];
- ps = REAL(params);
nreps = (nrepsp > nrepsx) ? nrepsp : nrepsx;
@@ -73,11 +65,7 @@
// vector for interpolated covariates
PROTECT(cvec = NEW_NUMERIC(ncovars)); nprotect++;
SET_NAMES(cvec,Cnames);
- cp = REAL(cvec);
- ndim[0] = nreps; ndim[1] = ntimes;
- PROTECT(F = makearray(2,ndim)); nprotect++;
-
// extract the user-defined function
PROTECT(fn = unpack_pomp_fun(fun,&mode)); nprotect++;
@@ -90,19 +78,12 @@
case 0: // R function
PROTECT(tvec = NEW_NUMERIC(1)); nprotect++;
- tp = REAL(tvec);
-
PROTECT(xvec = NEW_NUMERIC(nvars)); nprotect++;
+ PROTECT(yvec = NEW_NUMERIC(nobs)); nprotect++;
+ PROTECT(pvec = NEW_NUMERIC(npars)); nprotect++;
SET_NAMES(xvec,Snames);
- xp = REAL(xvec);
-
- PROTECT(yvec = NEW_NUMERIC(nobs)); nprotect++;
SET_NAMES(yvec,Onames);
- yp = REAL(yvec);
-
- PROTECT(pvec = NEW_NUMERIC(npars)); nprotect++;
SET_NAMES(pvec,Pnames);
- pp = REAL(pvec);
// set up the function call
PROTECT(fcall = LCONS(cvec,fcall)); nprotect++;
@@ -144,45 +125,63 @@
}
+ // create array to store results
+ {
+ int dim[2] = {nreps, ntimes};
+ PROTECT(F = makearray(2,dim)); nprotect++;
+ }
+
// now do computations
switch (mode) {
case 0: // R function
- for (k = 0, ft = REAL(F), ts = REAL(times); k < ntimes; k++, ts++) { // loop over times
+ {
+ int first = 1;
+ double *ys = REAL(y);
+ double *xs = REAL(x);
+ double *ps = REAL(params);
+ double *cp = REAL(cvec);
+ double *tp = REAL(tvec);
+ double *xp = REAL(xvec);
+ double *yp = REAL(yvec);
+ double *pp = REAL(pvec);
+ double *ft = REAL(F);
+ double *time = REAL(times);
+ int j, k;
- R_CheckUserInterrupt(); // check for user interrupt
+ for (k = 0; k < ntimes; k++, time++, ys += nobs) { // loop over times
- *tp = *ts; // copy the time
- table_lookup(&covariate_table,*tp,cp,0); // interpolate the covariates
+ R_CheckUserInterrupt(); // check for user interrupt
- for (i = 0; i < nobs; i++) yp[i] = ys[i+nobs*k];
-
- for (j = 0; j < nreps; j++, ft++) { // loop over replicates
+ *tp = *time; // copy the time
+ table_lookup(&covariate_table,*time,cp,0); // interpolate the covariates
- // copy the states and parameters into place
- for (i = 0; i < nvars; i++) xp[i] = xs[i+nvars*((j%nrepsx)+nrepsx*k)];
- for (i = 0; i < npars; i++) pp[i] = ps[i+npars*(j%nrepsp)];
+ memcpy(yp,ys,nobs*sizeof(double));
+
+ for (j = 0; j < nreps; j++, ft++) { // loop over replicates
+
+ // copy the states and parameters into place
+ memcpy(xp,&xs[nvars*((j%nrepsx)+nrepsx*k)],nvars*sizeof(double));
+ memcpy(pp,&ps[npars*(j%nrepsp)],npars*sizeof(double));
- if (first) {
- // evaluate the call
- PROTECT(ans = eval(fcall,rho)); nprotect++;
+ if (first) {
+ // evaluate the call
+ PROTECT(ans = eval(fcall,rho)); nprotect++;
+ if (LENGTH(ans) != 1)
+ error("user 'dmeasure' return a vector of length %d when it should return a scalar",LENGTH(ans));
- if (LENGTH(ans) != 1)
- error("user 'dmeasure' return a vector of length %d when it should return a scalar",LENGTH(ans));
+ *ft = *(REAL(AS_NUMERIC(ans)));
- fs = REAL(AS_NUMERIC(ans));
+ first = 0;
- first = 0;
+ } else {
- } else {
+ *ft = *(REAL(AS_NUMERIC(eval(fcall,rho))));
- fs = REAL(AS_NUMERIC(eval(fcall,rho)));
+ }
}
-
- *ft = *fs;
-
}
}
@@ -192,22 +191,31 @@
set_pomp_userdata(fcall);
- for (k = 0, ft = REAL(F), ts = REAL(times); k < ntimes; k++, ts++) { // loop over times
+ {
+ double *yp = REAL(y);
+ double *xs = REAL(x);
+ double *ps = REAL(params);
+ double *cp = REAL(cvec);
+ double *ft = REAL(F);
+ double *time = REAL(times);
+ double *xp, *pp;
+ int j, k;
- R_CheckUserInterrupt(); // check for user interrupt
+ for (k = 0; k < ntimes; k++, time++, yp += nobs) { // loop over times
+
+ R_CheckUserInterrupt(); // check for user interrupt
- // interpolate the covar functions for the covariates
- table_lookup(&covariate_table,*ts,cp,0);
+ // interpolate the covar functions for the covariates
+ table_lookup(&covariate_table,*time,cp,0);
- yp = &ys[nobs*k];
-
- for (j = 0; j < nreps; j++, ft++) { // loop over replicates
+ for (j = 0; j < nreps; j++, ft++) { // loop over replicates
- xp = &xs[nvars*((j%nrepsx)+nrepsx*k)];
- pp = &ps[npars*(j%nrepsp)];
+ xp = &xs[nvars*((j%nrepsx)+nrepsx*k)];
+ pp = &ps[npars*(j%nrepsp)];
- (*ff)(ft,yp,xp,pp,give_log,oidx,sidx,pidx,cidx,ncovars,cp,*ts);
+ (*ff)(ft,yp,xp,pp,give_log,oidx,sidx,pidx,cidx,ncovars,cp,*time);
+ }
}
}
Modified: pkg/pomp/src/skeleton.c
===================================================================
--- pkg/pomp/src/skeleton.c 2012-04-28 18:01:22 UTC (rev 686)
+++ pkg/pomp/src/skeleton.c 2012-04-29 14:24:52 UTC (rev 687)
@@ -12,37 +12,29 @@
{
int nprotect = 0;
int nvars, npars, nrepp, nrepx, nreps, ntimes, ncovars;
- int mode;
- int use_names;
- int *dim, ndim[3], *op;
+ int mode = -1;
+ int *dim;
SEXP fn, fcall, rho, ans, nm;
SEXP tvec, xvec, pvec, cvec;
SEXP Snames, Cnames, Pnames;
SEXP statenames, paramnames, covarnames;
SEXP F;
- double *xs, *ts, *ps, *fs, *ft;
- double *tp, *xp, *pp, *cp;
int *sidx, *pidx, *cidx;
- int first = 1;
pomp_skeleton *ff = NULL;
struct lookup_table covariate_table;
- int i, j, k;
PROTECT(t = AS_NUMERIC(t)); nprotect++;
ntimes = LENGTH(t);
- ts = REAL(t);
PROTECT(x = as_state_array(x)); nprotect++;
dim = INTEGER(GET_DIM(x));
nvars = dim[0]; nrepx = dim[1];
if (ntimes != dim[2])
error("skeleton error: length of 't' and 3rd dimension of 'x' do not agree");
- xs = REAL(x);
PROTECT(params = as_matrix(params)); nprotect++;
dim = INTEGER(GET_DIM(params));
npars = dim[0]; nrepp = dim[1];
- ps = REAL(params);
// 2nd dimension of 'x' and 'params' need not agree
nreps = (nrepp > nrepx) ? nrepp : nrepx;
@@ -59,13 +51,7 @@
// vector for interpolated covariates
PROTECT(cvec = NEW_NUMERIC(ncovars)); nprotect++;
SET_NAMES(cvec,Cnames);
- cp = REAL(cvec);
- // set up the array to be returned
- ndim[0] = nvars; ndim[1] = nreps; ndim[2] = ntimes;
- PROTECT(F = makearray(3,ndim)); nprotect++;
- setrownames(F,Snames,3);
-
// extract the user-defined function
PROTECT(fn = unpack_pomp_fun(fun,&mode)); nprotect++;
@@ -77,16 +63,11 @@
case 0: // R skeleton
PROTECT(tvec = NEW_NUMERIC(1)); nprotect++;
- tp = REAL(tvec);
-
PROTECT(xvec = NEW_NUMERIC(nvars)); nprotect++;
+ PROTECT(pvec = NEW_NUMERIC(npars)); nprotect++;
SET_NAMES(xvec,Snames);
- xp = REAL(xvec);
+ SET_NAMES(pvec,Pnames);
- PROTECT(pvec = NEW_NUMERIC(npars)); nprotect++;
- SET_NAMES(pvec,GET_ROWNAMES(GET_DIMNAMES(params)));
- pp = REAL(pvec);
-
// set up the function call
PROTECT(fcall = LCONS(cvec,fcall)); nprotect++;
SET_TAG(fcall,install("covars"));
@@ -118,54 +99,77 @@
}
+ // set up the array to hold results
+ {
+ int dim[3] = {nvars, nreps, ntimes};
+ PROTECT(F = makearray(3,dim)); nprotect++;
+ setrownames(F,Snames,3);
+ }
+
// now do computations
switch (mode) {
case 0: // R skeleton
- for (k = 0, ft = REAL(F); k < ntimes; k++, ts++) { // loop over times
+ {
+ int first = 1;
+ int use_names;
+ double *time = REAL(t);
+ double *xs = REAL(x);
+ double *ps = REAL(params);
+ double *cp = REAL(cvec);
+ double *tp = REAL(tvec);
+ double *xp = REAL(xvec);
+ double *pp = REAL(pvec);
+ double *ft = REAL(F);
+ double *fs;
+ int *posn;
+ int i, j, k;
- R_CheckUserInterrupt(); // check for user interrupt
+ for (k = 0; k < ntimes; k++, time++) { // loop over times
+
+ R_CheckUserInterrupt(); // check for user interrupt
- *tp = *ts; // copy the time
+ *tp = *time; // copy the time
- // interpolate the covar functions for the covariates
- table_lookup(&covariate_table,*tp,cp,0);
+ // interpolate the covar functions for the covariates
+ table_lookup(&covariate_table,*time,cp,0);
- for (j = 0; j < nreps; j++, ft += nvars) { // loop over replicates
+ for (j = 0; j < nreps; j++, ft += nvars) { // loop over replicates
- for (i = 0; i < nvars; i++) xp[i] = xs[i+nvars*((j%nrepx)+nrepx*k)];
- for (i = 0; i < npars; i++) pp[i] = ps[i+npars*(j%nrepp)];
+ memcpy(xp,&xs[nvars*((j%nrepx)+nrepx*k)],nvars*sizeof(double));
+ memcpy(pp,&ps[npars*(j%nrepp)],npars*sizeof(double));
- if (first) {
+ if (first) {
- PROTECT(ans = eval(fcall,rho)); nprotect++;
- if (LENGTH(ans)!=nvars)
- error("user 'skeleton' returns a vector of %d state variables but %d are expected",LENGTH(ans),nvars);
+ PROTECT(ans = eval(fcall,rho)); nprotect++;
+ if (LENGTH(ans)!=nvars)
+ error("user 'skeleton' returns a vector of %d state variables but %d are expected",LENGTH(ans),nvars);
- // get name information to fix possible alignment problems
- PROTECT(nm = GET_NAMES(ans)); nprotect++;
- use_names = !isNull(nm);
- if (use_names) {
- op = INTEGER(PROTECT(matchnames(Snames,nm))); nprotect++;
- } else {
- op = 0;
- }
+ // get name information to fix possible alignment problems
+ PROTECT(nm = GET_NAMES(ans)); nprotect++;
+ use_names = !isNull(nm);
+ if (use_names) {
+ posn = INTEGER(PROTECT(matchnames(Snames,nm))); nprotect++;
+ } else {
+ posn = 0;
+ }
- fs = REAL(AS_NUMERIC(ans));
+ fs = REAL(AS_NUMERIC(ans));
- first = 0;
+ first = 0;
- } else {
+ } else {
- fs = REAL(AS_NUMERIC(eval(fcall,rho)));
+ fs = REAL(AS_NUMERIC(eval(fcall,rho)));
- }
+ }
- if (use_names)
- for (i = 0; i < nvars; i++) ft[op[i]] = fs[i];
- else
- for (i = 0; i < nvars; i++) ft[i] = fs[i];
+ if (use_names)
+ for (i = 0; i < nvars; i++) ft[posn[i]] = fs[i];
+ else
+ for (i = 0; i < nvars; i++) ft[i] = fs[i];
+ }
}
}
@@ -175,20 +179,30 @@
set_pomp_userdata(fcall);
- for (k = 0, ft = REAL(F); k < ntimes; k++, ts++) { // loop over times
+ {
+ double *time = REAL(t);
+ double *xs = REAL(x);
+ double *ps = REAL(params);
+ double *cp = REAL(cvec);
+ double *ft = REAL(F);
+ double *xp, *pp;
+ int j, k;
- R_CheckUserInterrupt(); // check for user interrupt
+ for (k = 0; k < ntimes; k++, time++) { // loop over times
+
+ R_CheckUserInterrupt(); // check for user interrupt
- // interpolate the covar functions for the covariates
- table_lookup(&covariate_table,*ts,cp,0);
+ // interpolate the covar functions for the covariates
+ table_lookup(&covariate_table,*time,cp,0);
- for (j = 0; j < nreps; j++, ft += nvars) { // loop over replicates
+ for (j = 0; j < nreps; j++, ft += nvars) { // loop over replicates
- xp = &xs[nvars*((j%nrepx)+nrepx*k)];
- pp = &ps[npars*(j%nrepp)];
+ xp = &xs[nvars*((j%nrepx)+nrepx*k)];
+ pp = &ps[npars*(j%nrepp)];
- (*ff)(ft,xp,pp,sidx,pidx,cidx,ncovars,cp,*ts);
+ (*ff)(ft,xp,pp,sidx,pidx,cidx,ncovars,cp,*time);
+ }
}
}
More information about the pomp-commits
mailing list