[Pomp-commits] r375 - in pkg: inst/doc src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Oct 11 00:07:06 CEST 2010
Author: kingaa
Date: 2010-10-11 00:07:05 +0200 (Mon, 11 Oct 2010)
New Revision: 375
Modified:
pkg/inst/doc/advanced_topics_in_pomp.Rnw
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/doc/intro_to_pomp.pdf
pkg/src/rprocess.c
pkg/src/simulate.c
Log:
- accelerate 'rprocess' and 'simulate' by using 'memcpy' where appropriate
- change the order of topics in 'advanced_topics_in_pomp' vignette
Modified: pkg/inst/doc/advanced_topics_in_pomp.Rnw
===================================================================
--- pkg/inst/doc/advanced_topics_in_pomp.Rnw 2010-10-07 13:25:40 UTC (rev 374)
+++ pkg/inst/doc/advanced_topics_in_pomp.Rnw 2010-10-10 22:07:05 UTC (rev 375)
@@ -61,8 +61,7 @@
set.seed(5384959)
@
-This document serves to give some examples of the use of native (C or FORTRAN) codes in \code{pomp}
-and to introduce the low-level interface to \code{pomp} objects.
+This document gives some examples of the use of native (C or FORTRAN) codes in \code{pomp} and introduces the low-level interface to \code{pomp} objects.
\section{Acceleration using native codes: using plug-ins with native code}
@@ -122,7 +121,7 @@
The source code for the native routines \verb+_sir_euler_simulator+, \verb+_sir_ODE+, \verb+_sir_binom_rmeasure+, and \verb+_sir_binom_dmeasure+ is provided with the package (in the \code{examples} directory).
To see the source code, do
<<view-sir-source,eval=F>>=
-edit(file=system.file("examples/sir.c",package="pomp"))
+file.show(file=system.file("examples/sir.c",package="pomp"))
@
Also in the \code{examples} directory is an \R\ script that shows how to compile \verb+sir.c+ into a shared-object library and link it with \R.
Note that the native routines for this model are included in the package, which is why we give the \verb+PACKAGE="pomp"+ argument to \pomp.
@@ -248,7 +247,7 @@
For convenience, the source codes are provided with the package in the \code{examples} directory.
Do
<<view-ou2-source,eval=F>>=
-edit(file=system.file("examples/ou2.c",package="pomp"))
+file.show(file=system.file("examples/ou2.c",package="pomp"))
@
to view the source code.
Modified: pkg/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/src/rprocess.c
===================================================================
--- pkg/src/rprocess.c 2010-10-07 13:25:40 UTC (rev 374)
+++ pkg/src/rprocess.c 2010-10-10 22:07:05 UTC (rev 375)
@@ -4,6 +4,7 @@
#include <Rmath.h>
#include <Rdefines.h>
#include <Rinternals.h>
+#include <string.h>
#include "pomp_internal.h"
@@ -70,14 +71,10 @@
error("rprocess error: user 'rprocess' must return an array with rownames");
}
if (off > 0) {
- int i, n;
- double *s, *t;
xdim[2] -= off;
PROTECT(Xoff = makearray(3,xdim)); nprotect++;
setrownames(Xoff,GET_ROWNAMES(GET_DIMNAMES(X)),3);
- s = REAL(X)+off*nvars*nreps;
- t = REAL(Xoff);
- for (i = 0, n = nvars*nreps*(ntimes-off); i < n; i++, s++, t++) *t = *s;
+ memcpy(REAL(Xoff),REAL(X)+off*nvars*nreps,(ntimes-off)*nvars*nreps*sizeof(double));
UNPROTECT(nprotect);
return Xoff;
} else {
Modified: pkg/src/simulate.c
===================================================================
--- pkg/src/simulate.c 2010-10-07 13:25:40 UTC (rev 374)
+++ pkg/src/simulate.c 2010-10-10 22:07:05 UTC (rev 375)
@@ -1,12 +1,15 @@
// -*- C++ -*-
-#include "pomp_internal.h"
#include <Rdefines.h>
+#include <string.h>
+#include "pomp_internal.h"
+
SEXP simulation_computations (SEXP object, SEXP params, SEXP times, SEXP t0, SEXP nsim, SEXP obs, SEXP states)
{
int nprotect = 0;
- SEXP xstart, x, xx, x0, y, p, alltimes, coef, yy, offset;
+ SEXP xstart, x, y, alltimes, coef, yy, offset;
+ SEXP p = R_NilValue, x0 = R_NilValue, xx = R_NilValue;
SEXP ans, ans_names;
SEXP po, popo;
SEXP statenames, paramnames, obsnames, statedim, obsdim;
@@ -14,10 +17,10 @@
int qobs, qstates;
int *dim, dims[3];
double *s, *t, *xs, *xt, *ys, *yt, *ps, *pt, tt;
- int i, j, k;
+ int i, j, k, np, nx;
PROTECT(offset = NEW_INTEGER(1)); nprotect++;
- INTEGER(offset)[0] = 1;
+ *(INTEGER(offset)) = 1;
nsims = INTEGER(AS_INTEGER(nsim))[0]; // number of simulations per parameter set
if (LENGTH(nsim)>1)
@@ -25,8 +28,8 @@
if (nsims < 1)
return R_NilValue; // no work to do
- qobs = LOGICAL(AS_LOGICAL(obs))[0]; // 'obs' flag set?
- qstates = LOGICAL(AS_LOGICAL(states))[0]; // 'states' flag set?
+ qobs = *(LOGICAL(AS_LOGICAL(obs))); // 'obs' flag set?
+ qstates = *(LOGICAL(AS_LOGICAL(states))); // 'states' flag set?
PROTECT(paramnames = GET_ROWNAMES(GET_DIMNAMES(params))); nprotect++;
dim = INTEGER(GET_DIM(params));
@@ -47,7 +50,7 @@
error("if 'times' is empty, there is no work to do");
PROTECT(alltimes = NEW_NUMERIC(ntimes+1)); nprotect++;
- tt = REAL(t0)[0];
+ tt = *(REAL(t0));
s = REAL(times);
t = REAL(alltimes);
@@ -64,12 +67,8 @@
}
// call 'rprocess' to simulate state process
- if (nsims == 1) { // only one simulation per parameter set
+ if (nsims > 1) { // multiple simulations per parameter set
- PROTECT(x = do_rprocess(object,xstart,alltimes,params,offset)); nprotect++;
-
- } else { // nsim > 1
-
dims[0] = npars; dims[1] = nreps;
PROTECT(p = makearray(2,dims)); nprotect++;
setrownames(p,paramnames,2);
@@ -79,15 +78,19 @@
setrownames(x0,statenames,2);
// make 'nsims' copies of the parameters and initial states
- for (k = 0, pt = REAL(p), xt = REAL(x0); k < nsims; k++) {
- for (j = 0, ps = REAL(params), xs = REAL(xstart); j < nparsets; j++) {
- for (i = 0; i < npars; i++, pt++, ps++) *pt = *ps;
- for (i = 0; i < nvars; i++, xt++, xs++) *xt = *xs;
- }
+ ps = REAL(params); np = LENGTH(params);
+ xs = REAL(xstart); nx = LENGTH(xstart);
+ for (k = 0, pt = REAL(p), xt = REAL(x0); k < nsims; k++, pt += np, xt += nx) {
+ memcpy(pt,ps,np*sizeof(double));
+ memcpy(xt,xs,nx*sizeof(double));
}
PROTECT(x = do_rprocess(object,x0,alltimes,p,offset)); nprotect++;
+ } else { // nsim == 1
+
+ PROTECT(x = do_rprocess(object,xstart,alltimes,params,offset)); nprotect++;
+
}
if (!qobs && qstates) { // obs=F,states=T: return states only
@@ -95,12 +98,12 @@
UNPROTECT(nprotect);
return x;
- } else {
+ } else { // we must do 'rmeasure'
- if (nsims == 1) {
+ if (nsims > 1) {
+ PROTECT(y = do_rmeasure(object,x,times,p)); nprotect++;
+ } else {
PROTECT(y = do_rmeasure(object,x,times,params)); nprotect++;
- } else {
- PROTECT(y = do_rmeasure(object,x,times,p)); nprotect++;
}
if (qobs) {
@@ -109,9 +112,9 @@
PROTECT(ans = NEW_LIST(2)); nprotect++;
PROTECT(ans_names = NEW_CHARACTER(2)); nprotect++;
- SET_NAMES(ans,ans_names);
SET_STRING_ELT(ans_names,0,mkChar("states"));
SET_STRING_ELT(ans_names,1,mkChar("obs"));
+ SET_NAMES(ans,ans_names);
SET_ELEMENT(ans,0,x);
SET_ELEMENT(ans,1,y);
UNPROTECT(nprotect);
@@ -157,7 +160,7 @@
ps = REAL(params);
pt = REAL(GET_SLOT(po,install("params")));
- for (i = 0; i < npars; i++, ps++, pt++) *pt = *ps;
+ memcpy(pt,ps,npars*sizeof(double));
UNPROTECT(nprotect);
return po;
@@ -168,32 +171,32 @@
PROTECT(ans = NEW_LIST(nreps)); nprotect++;
// create an array for the 'states' slot
- PROTECT(xx = makearray(2,INTEGER(statedim))); nprotect++;
+ PROTECT(xx = makearray(2,INTEGER(statedim)));
setrownames(xx,statenames,2);
SET_SLOT(po,install("states"),xx);
+ UNPROTECT(1);
// create an array for the 'data' slot
- PROTECT(yy = makearray(2,INTEGER(obsdim))); nprotect++;
+ PROTECT(yy = makearray(2,INTEGER(obsdim)));
setrownames(yy,obsnames,2);
SET_SLOT(po,install("data"),yy);
+ UNPROTECT(1);
xs = REAL(x);
ys = REAL(y);
- if (nsims == 1)
- ps = REAL(params);
- else
- ps = REAL(p);
+ ps = (nsims > 1) ? REAL(p) : REAL(params);
- for (k = 0; k < nreps; k++) { // loop over replicates
+ for (k = 0; k < nreps; k++, ps += npars) { // loop over replicates
PROTECT(popo = duplicate(po));
+ // copy parameters
+ pt = REAL(GET_SLOT(popo,install("params")));
+ memcpy(pt,ps,npars*sizeof(double));
+
// copy x[,k,] and y[,k,] into popo
xt = REAL(GET_SLOT(popo,install("states")));
yt = REAL(GET_SLOT(popo,install("data")));
- pt = REAL(GET_SLOT(popo,install("params")));
-
- for (i = 0; i < npars; i++, pt++, ps++) *pt = *ps;
for (j = 0; j < ntimes; j++) {
for (i = 0; i < nvars; i++, xt++) *xt = xs[i+nvars*(k+nreps*j)];
for (i = 0; i < nobs; i++, yt++) *yt = ys[i+nobs*(k+nreps*j)];
More information about the pomp-commits
mailing list