[Pomp-commits] r647 - in pkg: . R inst src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Apr 7 22:47:07 CEST 2012
Author: kingaa
Date: 2012-04-07 22:47:07 +0200 (Sat, 07 Apr 2012)
New Revision: 647
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/trajectory-pomp.R
pkg/inst/TODO
pkg/src/pomp_internal.h
pkg/src/skeleton.c
Log:
- put vectorfield evaluation function (for ODE integration of deterministic skeleton) into C for speed
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2012-04-06 12:41:58 UTC (rev 646)
+++ pkg/DESCRIPTION 2012-04-07 20:47:07 UTC (rev 647)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.41-2
-Date: 2012-04-06
+Version: 0.41-3
+Date: 2012-04-07
Revision: $Rev$
Author: Aaron A. King, Edward L. Ionides, Carles Breto, Steve Ellner, Bruce Kendall, Helen Wearing, Matthew J. Ferrari, Michael Lavine, Daniel C. Reuman
Maintainer: Aaron A. King <kingaa at umich.edu>
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2012-04-06 12:41:58 UTC (rev 646)
+++ pkg/NAMESPACE 2012-04-07 20:47:07 UTC (rev 647)
@@ -18,6 +18,8 @@
probe_acf,probe_ccf,
probe_nlar,
synth_loglik,
+ pomp_desolve_init,
+ pomp_vf_eval,
do_partrans,
do_rprocess,
do_dprocess,
Modified: pkg/R/trajectory-pomp.R
===================================================================
--- pkg/R/trajectory-pomp.R 2012-04-06 12:41:58 UTC (rev 646)
+++ pkg/R/trajectory-pomp.R 2012-04-07 20:47:07 UTC (rev 647)
@@ -56,6 +56,9 @@
dim(x0) <- c(nvar,nrep,1)
dimnames(x0) <- list(statenames,NULL,NULL)
+ znames <- object at zeronames
+ if (length(znames)>0) x0[znames,,] <- 0
+
type <- object at skeleton.type # map or vectorfield?
if (is.na(type))
@@ -68,37 +71,17 @@
} else if (type=="vectorfield") {
skel <- get.pomp.fun(object at skeleton)
+ .Call(pomp_desolve_init,object,params,skel,statenames,nvar,nrep);
- ## vectorfield function (RHS of ODE) in 'deSolve' format
- vf.eval <- function (t, y, ...) {
- list(
- .Call(
- do_skeleton,
- object,
- x=array(
- data=y,
- dim=c(nvar,nrep,1),
- dimnames=list(statenames,NULL,NULL)
- ),
- t=t,
- params=params,
- skel=skel
- ),
- NULL
- )
- }
-
- znames <- object at zeronames
-
- if (length(znames)>0)
- x0[znames,,] <- 0
-
X <- try(
ode(
y=x0,
times=c(t0,times),
- func=vf.eval,
method="lsoda",
+ func="pomp_vf_eval",
+ dllname="pomp",
+ initfunc=NULL,
+ parms=NULL,
...
),
silent=FALSE
@@ -110,15 +93,15 @@
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")
}
+ if (length(znames)>0)
+ x[znames,,-1] <- apply(x[znames,,,drop=FALSE],c(1,2),diff)
+
if (as.data.frame) {
x <- lapply(
seq_len(ncol(x)),
Modified: pkg/inst/TODO
===================================================================
--- pkg/inst/TODO 2012-04-06 12:41:58 UTC (rev 646)
+++ pkg/inst/TODO 2012-04-07 20:47:07 UTC (rev 647)
@@ -4,8 +4,6 @@
* sorting vs no-sorting systematic resampling algorithm?
-* native C routine for vectorfield evaluation inside lsoda
-
* parameter transformations: put 'transform' option into each estimation routine
(spect.match)
Modified: pkg/src/pomp_internal.h
===================================================================
--- pkg/src/pomp_internal.h 2012-04-06 12:41:58 UTC (rev 646)
+++ pkg/src/pomp_internal.h 2012-04-07 20:47:07 UTC (rev 647)
@@ -84,7 +84,7 @@
PROTECT(index = match(x,names,0)); nprotect++;
idx = INTEGER(index);
for (k = 0; k < n; k++) {
- if (idx[k]==0) error("variable %s not found",CHARACTER_DATA(STRING_ELT(nm,k)));
+ if (idx[k]==0) error("variable '%s' not found",CHARACTER_DATA(STRING_ELT(nm,k)));
idx[k] -= 1;
}
UNPROTECT(nprotect);
Modified: pkg/src/skeleton.c
===================================================================
--- pkg/src/skeleton.c 2012-04-06 12:41:58 UTC (rev 646)
+++ pkg/src/skeleton.c 2012-04-07 20:47:07 UTC (rev 647)
@@ -5,6 +5,7 @@
#include <Rmath.h>
#include <Rdefines.h>
#include <Rinternals.h>
+#include <R_ext/Rdynload.h>
static void eval_skel (pomp_skeleton *vf,
double *f,
@@ -265,3 +266,48 @@
#undef RHO
#undef FCALL
#undef VNAMES
+
+
+static SEXP *_pomp_vf_eval_object;
+static SEXP *_pomp_vf_eval_params;
+static SEXP *_pomp_vf_eval_skelfun;
+static SEXP *_pomp_vf_eval_xnames;
+static int _pomp_vf_eval_xdim[3];
+#define OBJECT (_pomp_vf_eval_object)
+#define PARAMS (_pomp_vf_eval_params)
+#define SKELFUN (_pomp_vf_eval_skelfun)
+#define XNAMES (_pomp_vf_eval_xnames)
+#define XDIM (_pomp_vf_eval_xdim)
+
+void pomp_desolve_init (SEXP object, SEXP params, SEXP fun, SEXP statenames, SEXP nvar, SEXP nrep) {
+ OBJECT = &object;
+ PARAMS = ¶ms;
+ SKELFUN = &fun;
+ XNAMES = &statenames;
+ XDIM[0] = INTEGER(AS_INTEGER(nvar))[0];
+ XDIM[1] = INTEGER(AS_INTEGER(nrep))[0];
+ XDIM[2] = 1;
+}
+
+
+void pomp_vf_eval (int *neq, double *t, double *y, double *ydot, double *yout, int *ip)
+{
+ SEXP T, X, dXdt;
+ int dim[3];
+
+ PROTECT(T = NEW_NUMERIC(1));
+ PROTECT(X = makearray(3,XDIM));
+ setrownames(X,*(XNAMES),3);
+ REAL(T)[0] = *t;
+ memcpy(REAL(X),y,(*neq)*sizeof(double));
+ PROTECT(dXdt = do_skeleton(*(OBJECT),X,T,*(PARAMS),*(SKELFUN)));
+ memcpy(ydot,REAL(dXdt),(*neq)*sizeof(double));
+
+ UNPROTECT(3);
+}
+
+#undef XDIM
+#undef XNAMES
+#undef SKELFUN
+#undef PARAMS
+#undef OBJECT
More information about the pomp-commits
mailing list