[Pomp-commits] r573 - in pkg: R src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jan 5 20:33:01 CET 2012
Author: kingaa
Date: 2012-01-05 20:33:00 +0100 (Thu, 05 Jan 2012)
New Revision: 573
Modified:
pkg/R/trajectory-pomp.R
pkg/src/euler.c
pkg/src/pomp_internal.h
pkg/src/sir.c
pkg/src/trajectory.c
pkg/tests/sir.R
pkg/tests/sir.Rout.save
Log:
- put calculation of Euler steps into a separate routine, now used both by the Euler plugin and by the deterministic map iterator
- 'trajectory' can now make use of 'zeronames' a la the Euler plugin
- 'trajectory' for maps can now handle a general timestep
Modified: pkg/R/trajectory-pomp.R
===================================================================
--- pkg/R/trajectory-pomp.R 2012-01-05 13:42:50 UTC (rev 572)
+++ pkg/R/trajectory-pomp.R 2012-01-05 19:33:00 UTC (rev 573)
@@ -56,12 +56,22 @@
type <- object at skeleton.type # map or vectorfield?
+ if ("zeronames"%in%names(object at userdata))
+ znames <- object at userdata$zeronames
+ else
+ znames <- character(0)
+
if (is.na(type))
stop("trajectory error: no skeleton specified",call.=FALSE)
if (type=="map") {
- x <- .Call(iterate_map,object,times,t0,x0,params)
+ if ("skelmap.delta.t"%in%names(object at userdata))
+ dt <- as.numeric(object at userdata$skelmap.delta.t)
+ else
+ dt <- 1
+
+ x <- .Call(iterate_map,object,times,t0,x0,params,dt,znames)
} else if (type=="vectorfield") {
@@ -86,6 +96,9 @@
)
}
+ if (length(znames)>0)
+ x0[znames,,] <- 0
+
X <- try(
ode(
y=x0,
@@ -103,6 +116,9 @@
x <- array(data=t(X[-1,-1]),dim=c(nvar,nrep,ntimes),dimnames=list(statenames,NULL,NULL))
+ if (length(znames)>0)
+ x[znames,,-1] <- apply(x[znames,,,drop=FALSE],c(1,2),diff)
+
} else {
stop("deterministic skeleton not specified")
Modified: pkg/src/euler.c
===================================================================
--- pkg/src/euler.c 2012-01-05 13:42:50 UTC (rev 572)
+++ pkg/src/euler.c 2012-01-05 19:33:00 UTC (rev 573)
@@ -3,6 +3,31 @@
#include "pomp_internal.h"
#include <R_ext/Constants.h>
+int num_euler_steps (double t1, double t2, double *dt) {
+ double tol = sqrt(DOUBLE_EPS);
+ int nstep;
+ // nstep will be the number of Euler steps to take in going from t1 to t2.
+ // note also that the stepsize changes.
+ // this choice is meant to be conservative
+ // (i.e., so that the actual dt does not exceed the specified dt
+ // by more than the relative tolerance 'tol')
+ // and to counteract roundoff error.
+ // It seems to work well, but is not guaranteed:
+ // suggestions would be appreciated.
+
+ if (t1 >= t2) {
+ *dt = 0.0;
+ nstep = 0;
+ } else if (t1+*dt >= t2) {
+ *dt = t2-t1;
+ nstep = 1;
+ } else {
+ nstep = (int) ceil((t2-t1)/(*dt)/(1+tol));
+ *dt = (t2-t1)/((double) nstep);
+ }
+ return nstep;
+}
+
// take nstep Euler-Poisson steps from t1 to t2
static void euler_simulator (pomp_onestep_sim *estep,
double *x, double *xstart, double *times, double *params,
@@ -24,8 +49,6 @@
struct lookup_table covariate_table = {covlen, covdim, 0, time_table, covar_table};
- tol = sqrt(DOUBLE_EPS); // relative tolerance in choosing Euler stepsize
-
// copy the start values into the result array
for (p = 0; p < nrep; p++)
for (k = 0; k < nvar; k++)
@@ -39,29 +62,12 @@
t = times[step-1];
dt = *deltat;
- // neuler is the number of Euler steps to take in going from
- // times[step-1] to times[step].
- // this choice is meant to be conservative
- // (i.e., so that the actual dt does not exceed the specified dt
- // by more than the relative tolerance 'tol')
- // and to counteract roundoff error.
- // It seems to work well, but is not guaranteed:
- // suggestions would be appreciated.
-
if (t > times[step]) {
error("'times' is not an increasing sequence");
}
- else if (t == times[step]) {
- dt = 0.0;
- neuler = 0;
- } else if (t+dt >= times[step]) {
- dt = times[step] - t;
- neuler = 1;
- } else {
- neuler = (int) ceil((times[step]-t)/dt/(1+tol));
- dt = (times[step]-t)/((double) neuler);
- }
+ neuler = num_euler_steps(t,times[step],&dt);
+
for (p = 0; p < nrep; p++) {
xp = &x[nvar*(p+nrep*step)];
// copy in the previous values of the state variables
Modified: pkg/src/pomp_internal.h
===================================================================
--- pkg/src/pomp_internal.h 2012-01-05 13:42:50 UTC (rev 572)
+++ pkg/src/pomp_internal.h 2012-01-05 19:33:00 UTC (rev 573)
@@ -22,6 +22,10 @@
double *y;
};
+// routine to compute number of Euler steps to take.
+// used by plugins in 'euler.c' and map iterator in 'trajectory.c'
+int num_euler_steps (double t1, double t2, double *dt);
+
// simple linear interpolation of the lookup table (with derivative if desired)
// setting dydt = 0 in the call to 'table_lookup' will bypass computation of the derivative
void table_lookup (struct lookup_table *tab, double x, double *y, double *dydt);
Modified: pkg/src/sir.c
===================================================================
--- pkg/src/sir.c 2012-01-05 13:42:50 UTC (rev 572)
+++ pkg/src/sir.c 2012-01-05 19:33:00 UTC (rev 573)
@@ -202,7 +202,7 @@
DSDT = term[0]-term[1]-term[2];
DIDT = term[1]-term[3]-term[4];
DRDT = term[3]-term[5];
- DCDT = DIDT*gamma/52; // roughly the weekly number of new cases
+ DCDT = term[3]; // accumulate the new I->R transitions
}
Modified: pkg/src/trajectory.c
===================================================================
--- pkg/src/trajectory.c 2012-01-05 13:42:50 UTC (rev 572)
+++ pkg/src/trajectory.c 2012-01-05 19:33:00 UTC (rev 573)
@@ -1,6 +1,7 @@
// -*- C++ -*-
#include <Rdefines.h>
+#include <R_ext/Constants.h>
#include "pomp_internal.h"
@@ -29,15 +30,19 @@
return ans;
}
-SEXP iterate_map (SEXP object, SEXP times, SEXP t0, SEXP x0, SEXP params)
+SEXP iterate_map (SEXP object, SEXP times, SEXP t0, SEXP x0, SEXP params,
+ SEXP dt, SEXP zeronames)
{
int nprotect = 0;
SEXP ans;
- SEXP x, f, skel, time;
+ SEXP x, f, skel, time, zindex;
int nvars, nreps, ntimes;
- int i, j, k;
+ int nzeros, nsteps;
+ int h, i, j, k;
int ndim[3];
double *tm, *tp, *xp, *fp, *ap;
+ int *zidx;
+ double deltat;
PROTECT(t0 = AS_NUMERIC(t0)); nprotect++;
PROTECT(x = duplicate(AS_NUMERIC(x0))); nprotect++;
@@ -63,24 +68,57 @@
PROTECT(skel = get_pomp_fun(GET_SLOT(object,install("skeleton")))); nprotect++;
+ PROTECT(dt = AS_NUMERIC(dt)); nprotect++;
+ if (REAL(dt)[0]<=0)
+ error("timestep 'dt' must be positive!");
+
+ PROTECT(zeronames = AS_CHARACTER(zeronames)); nprotect++;
+ nzeros = LENGTH(zeronames);
+
+ if (nzeros>0) {
+ PROTECT(zindex = MATCHROWNAMES(x0,zeronames)); nprotect++;
+ zidx = INTEGER(zindex);
+ } else {
+ zidx = 0;
+ }
+
+ if (nzeros>0) {
+ for (j = 0; j < nreps; j++) {
+ for (i = 0; i < nzeros; i++) {
+ xp[zidx[i]+nvars*j] = 0.0;
+ }
+ }
+ }
+
for (k = 0; k < ntimes; k++, tp++) {
- if (*tm < *tp) {
- while (*tm < *tp) {
- PROTECT(f = do_skeleton(object,x,time,params,skel));
- fp = REAL(f);
- for (j = 0; j < nreps; j++) {
- for (i = 0; i < nvars; i++) {
- xp[i+nvars*j] = fp[i+nvars*j];
- }
+
+ if (k > 0) *tm = *(tp-1);
+ // compute number of steps to take
+ deltat = REAL(dt)[0];
+ nsteps = num_euler_steps(*tm,*tp,&deltat);
+
+ for (h = 0; h < nsteps; h++) {
+ PROTECT(f = do_skeleton(object,x,time,params,skel));
+ fp = REAL(f);
+ for (j = 0; j < nreps; j++) {
+ for (i = 0; i < nvars; i++) {
+ xp[i+nvars*j] = fp[i+nvars*j];
}
- *tm += 1.0;
- UNPROTECT(1);
}
+ UNPROTECT(1);
+ *tm += deltat;
+ if (h == nsteps-2) { // penultimate step
+ deltat = *(tp+1)-*tm;
+ *tm = *(tp+1)-deltat;
+ }
}
for (j = 0; j < nreps; j++) {
for (i = 0; i < nvars; i++) {
ap[i+nvars*(j+nreps*k)] = xp[i+nvars*j];
}
+ for (i = 0; i < nzeros; i++) {
+ xp[zidx[i]+nvars*j] = 0.0;
+ }
}
}
Modified: pkg/tests/sir.R
===================================================================
--- pkg/tests/sir.R 2012-01-05 13:42:50 UTC (rev 572)
+++ pkg/tests/sir.R 2012-01-05 19:33:00 UTC (rev 573)
@@ -102,7 +102,7 @@
terms[1]-terms[2]-terms[3],
terms[2]-terms[4]-terms[5],
terms[4]-terms[6],
- gamma/52*(terms[2]-terms[4]-terms[5])
+ terms[4]
)
xdot
}
Modified: pkg/tests/sir.Rout.save
===================================================================
--- pkg/tests/sir.Rout.save 2012-01-05 13:42:50 UTC (rev 572)
+++ pkg/tests/sir.Rout.save 2012-01-05 19:33:00 UTC (rev 573)
@@ -1,5 +1,5 @@
-R version 2.13.1 (2011-07-08)
+R version 2.14.1 (2011-12-22)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-unknown-linux-gnu (64-bit)
@@ -123,7 +123,7 @@
+ terms[1]-terms[2]-terms[3],
+ terms[2]-terms[4]-terms[5],
+ terms[4]-terms[6],
-+ gamma/52*(terms[2]-terms[4]-terms[5])
++ terms[4]
+ )
+ xdot
+ }
@@ -394,7 +394,7 @@
zeronames = zeronames, tcovar = tcovar, covar = covar,
args = pairlist(...))
}
-<environment: 0x14e87b8>
+<environment: 0x21be458>
process model density, dprocess =
function (x, times, params, ..., statenames = character(0), paramnames = character(0),
covarnames = character(0), tcovar, covar, log = FALSE)
@@ -404,7 +404,7 @@
covarnames = covarnames, tcovar = tcovar, covar = covar,
log = log, args = pairlist(...))
}
-<environment: 0x193a340>
+<environment: 0x2132e10>
measurement model simulator, rmeasure =
function (x, t, params, covars, ...)
{
@@ -442,8 +442,7 @@
terms <- c(mu * pop, foi * S, mu * S, gamma * I, mu *
I, mu * R)
xdot[1:4] <- c(terms[1] - terms[2] - terms[3], terms[2] -
- terms[4] - terms[5], terms[4] - terms[6], gamma/52 *
- (terms[2] - terms[4] - terms[5]))
+ terms[4] - terms[5], terms[4] - terms[6], terms[4])
xdot
})
}
@@ -463,11 +462,11 @@
parameter transform function =
function (params, ...)
params
-<environment: 0x15beeb0>
+<environment: 0x237d2c8>
parameter inverse transform function =
function (params, ...)
params
-<environment: 0x15beeb0>
+<environment: 0x237d2c8>
userdata =
$zeronames
[1] "cases"
@@ -479,7 +478,7 @@
> x <- simulate(po,params=log(params),nsim=3)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 2.027774 secs
+Time difference of 1.124028 secs
>
> pdf(file='sir.pdf')
>
@@ -496,7 +495,7 @@
> X3 <- trajectory(po,params=log(params),times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 2.722191 secs
+Time difference of 1.465972 secs
> plot(t3,X3['I',1,],type='l')
>
> f1 <- dprocess(
@@ -556,7 +555,7 @@
> x <- simulate(po,nsim=100)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 2.520875 secs
+Time difference of 1.205045 secs
> plot(x[[1]],variables=c("S","I","R","cases","W"))
>
> t3 <- seq(0,20,by=1/52)
@@ -564,7 +563,7 @@
> X4 <- trajectory(po,times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 0.1533842 secs
+Time difference of 0.07390642 secs
> plot(t3,X4['I',1,],type='l')
>
> g2 <- dmeasure(
More information about the pomp-commits
mailing list