[Pomp-commits] r121 - in pkg: inst/doc inst/include src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu May 7 22:41:46 CEST 2009
Author: kingaa
Date: 2009-05-07 22:41:44 +0200 (Thu, 07 May 2009)
New Revision: 121
Modified:
pkg/inst/doc/compiled_code_in_pomp.pdf
pkg/inst/doc/intro_to_pomp.pdf
pkg/inst/include/pomp.h
pkg/src/eulermultinom.c
pkg/src/pomp.h
pkg/tests/sir.R
pkg/tests/sir.Rout.save
Log:
streamline "reulermultinom" by moving GetRNGstate() and PutRNGstate() outside of a level of looping. If you want to use 'reulermultinom' at the C level, you must take care that it is wrapped inside of the RNGstate functions. An approx. 10-fold increase in speed for the euler.simulate function is realized by this simple change.
Modified: pkg/inst/doc/compiled_code_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/include/pomp.h
===================================================================
--- pkg/inst/include/pomp.h 2009-05-05 22:45:04 UTC (rev 120)
+++ pkg/inst/include/pomp.h 2009-05-07 20:41:44 UTC (rev 121)
@@ -8,22 +8,13 @@
#include <Rdefines.h>
// prototypes for C-level access to Euler-multinomial distribution functions
+// NB: 'reulermultinom' does not call GetRNGstate() and PutRNGstate() internally, so the user must do so
void reulermultinom (int ntrans, double size, double *rate, double dt, double *trans);
double deulermultinom (int ntrans, double size, double *rate, double dt, double *trans, int give_log);
// facility for dotting a vector of parameters ('coef') against a vector of basis-function values ('basis')
double dot_product (int dim, const double *basis, const double *coef);
-// Prototype for stochastic simulation algorithm reaction-rate function.
-typedef double pomp_ssa_rate_fn(int j, double t, const double *x, const double *p);
-// Description:
-// on input:
-// j = integer specifying the number of the reaction whose rate is desired
-// t = time at which the rates are to be evaluated
-// x = vector of state variables
-// p = vector of parameters
-// returns the rate of the j-th reaction
-
// Prototype for one-step simulator, as used by "euler.simulate" and "onestep.simulate":
typedef void pomp_onestep_sim(double *x, const double *p,
const int *stateindex, const int *parindex, const int *covindex,
Modified: pkg/src/eulermultinom.c
===================================================================
--- pkg/src/eulermultinom.c 2009-05-05 22:45:04 UTC (rev 120)
+++ pkg/src/eulermultinom.c 2009-05-07 20:41:44 UTC (rev 121)
@@ -21,7 +21,6 @@
}
p += rate[k]; // total event rate
}
- GetRNGstate();
size = rbinom(size,1-exp(-p*dt)); // total number of events
if (!(R_FINITE(size)))
warning("reulermultinom: result of binomial draw is not finite");
@@ -35,7 +34,6 @@
p -= rate[k];
}
trans[m] = size;
- PutRNGstate();
}
// probability density of Euler-multinomial transitions
@@ -107,9 +105,11 @@
PROTECT(nn = AS_INTEGER(n)); nprotect++;
dim[0] = ntrans;
dim[1] = INTEGER(nn)[0];
- PROTECT(x = makearray(2,dim)); nprotect++;
+ PROTECT(x = makearray(2,dim)); nprotect++;
setrownames(x,GET_NAMES(rate),2);
+ GetRNGstate();
reulermultinom_multi(INTEGER(nn),&ntrans,REAL(size),REAL(rate),REAL(dt),REAL(x));
+ PutRNGstate();
UNPROTECT(nprotect);
return x;
}
@@ -127,7 +127,7 @@
warning("deulermultinom: only the first element of 'size' is meaningful");
if (length(dt)>1)
warning("deulermultinom: only the first element of 'dt' is meaningful");
- PROTECT(f = NEW_NUMERIC(n)); nprotect++;
+ PROTECT(f = NEW_NUMERIC(n)); nprotect++;
deulermultinom_multi(&n,&ntrans,REAL(size),REAL(rate),REAL(dt),REAL(x),INTEGER(log),REAL(f));
UNPROTECT(nprotect);
return f;
Modified: pkg/src/pomp.h
===================================================================
--- pkg/src/pomp.h 2009-05-05 22:45:04 UTC (rev 120)
+++ pkg/src/pomp.h 2009-05-07 20:41:44 UTC (rev 121)
@@ -8,22 +8,13 @@
#include <Rdefines.h>
// prototypes for C-level access to Euler-multinomial distribution functions
+// NB: 'reulermultinom' does not call GetRNGstate() and PutRNGstate() internally, so the user must do so
void reulermultinom (int ntrans, double size, double *rate, double dt, double *trans);
double deulermultinom (int ntrans, double size, double *rate, double dt, double *trans, int give_log);
// facility for dotting a vector of parameters ('coef') against a vector of basis-function values ('basis')
double dot_product (int dim, const double *basis, const double *coef);
-// Prototype for stochastic simulation algorithm reaction-rate function.
-typedef double pomp_ssa_rate_fn(int j, double t, const double *x, const double *p);
-// Description:
-// on input:
-// j = integer specifying the number of the reaction whose rate is desired
-// t = time at which the rates are to be evaluated
-// x = vector of state variables
-// p = vector of parameters
-// returns the rate of the j-th reaction
-
// Prototype for one-step simulator, as used by "euler.simulate" and "onestep.simulate":
typedef void pomp_onestep_sim(double *x, const double *p,
const int *stateindex, const int *parindex, const int *covindex,
Modified: pkg/tests/sir.R
===================================================================
--- pkg/tests/sir.R 2009-05-05 22:45:04 UTC (rev 120)
+++ pkg/tests/sir.R 2009-05-07 20:41:44 UTC (rev 121)
@@ -189,7 +189,7 @@
set.seed(3049953)
## simulate from the model
tic <- Sys.time()
-x <- simulate(euler.sir,nsim=3)
+x <- simulate(euler.sir,nsim=100)
toc <- Sys.time()
print(toc-tic)
plot(x[[1]],variables=c("S","I","R","cases","W"))
Modified: pkg/tests/sir.Rout.save
===================================================================
--- pkg/tests/sir.Rout.save 2009-05-05 22:45:04 UTC (rev 120)
+++ pkg/tests/sir.Rout.save 2009-05-07 20:41:44 UTC (rev 121)
@@ -148,7 +148,7 @@
> x <- simulate(po,params=log(params),nsim=3)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 2.599395 secs
+Time difference of 4.062482 secs
>
> pdf(file='sir.pdf')
>
@@ -165,7 +165,7 @@
> X3 <- trajectory(po,params=log(params),times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 5.949128 secs
+Time difference of 9.65457 secs
> plot(t3,X3['I',1,],type='l')
>
> f1 <- dprocess(
@@ -220,10 +220,10 @@
> set.seed(3049953)
> ## simulate from the model
> tic <- Sys.time()
-> x <- simulate(euler.sir,nsim=3)
+> x <- simulate(euler.sir,nsim=100)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 0.1701131 secs
+Time difference of 2.457672 secs
> plot(x[[1]],variables=c("S","I","R","cases","W"))
>
> t3 <- seq(0,20,by=1/52)
@@ -231,7 +231,7 @@
> X4 <- trajectory(euler.sir,times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 4.644418 secs
+Time difference of 7.418138 secs
> plot(t3,X4['I',1,],type='l')
>
> f2 <- dprocess(
More information about the pomp-commits
mailing list