[Pomp-commits] r625 - in pkg: . R data inst inst/doc inst/include man src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Mar 9 01:41:25 CET 2012
Author: kingaa
Date: 2012-03-09 01:41:25 +0100 (Fri, 09 Mar 2012)
New Revision: 625
Added:
pkg/src/partrans.c
pkg/tests/dimchecks.R
pkg/tests/dimchecks.Rout.save
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/init-state-pomp.R
pkg/R/pomp-fun.R
pkg/R/pomp-methods.R
pkg/R/pomp.R
pkg/data/blowflies.rda
pkg/data/dacca.rda
pkg/data/euler.sir.rda
pkg/data/gillespie.sir.rda
pkg/data/gompertz.rda
pkg/data/ou2.rda
pkg/data/ricker.rda
pkg/data/rw2.rda
pkg/data/verhulst.rda
pkg/inst/NEWS
pkg/inst/doc/advanced_topics_in_pomp.Rnw
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/doc/gompertz-multi-mif.rda
pkg/inst/doc/gompertz-trajmatch.rda
pkg/inst/doc/intro_to_pomp.pdf
pkg/inst/doc/nlf-block-boot.rda
pkg/inst/doc/nlf-boot.rda
pkg/inst/doc/nlf-fit-from-truth.rda
pkg/inst/doc/nlf-fits.rda
pkg/inst/doc/nlf-lag-tests.rda
pkg/inst/doc/nlf-multi-short.rda
pkg/inst/doc/ricker-mif.rda
pkg/inst/doc/ricker-probe-match.rda
pkg/inst/include/pomp.h
pkg/man/pomp-class.Rd
pkg/src/dmeasure.c
pkg/src/dprocess.c
pkg/src/initstate.c
pkg/src/pomp.h
pkg/src/pomp_internal.h
pkg/src/rmeasure.c
pkg/src/rprocess.c
pkg/src/skeleton.c
pkg/tests/logistic.R
pkg/tests/logistic.Rout.save
pkg/tests/ou2-procmeas.R
pkg/tests/ou2-procmeas.Rout.save
pkg/tests/rw2.Rout.save
pkg/tests/sir.Rout.save
Log:
- parameter transformations are handled through 'partrans'
- all the horsemen now recycle parameters or states as appropriate
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2012-03-09 00:22:08 UTC (rev 624)
+++ pkg/DESCRIPTION 2012-03-09 00:41:25 UTC (rev 625)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.40-7
-Date: 2012-03-07
+Version: 0.40-8
+Date: 2012-03-08
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>
Property changes on: pkg/DESCRIPTION
___________________________________________________________________
Added: keywords
+ Rev
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2012-03-09 00:22:08 UTC (rev 624)
+++ pkg/NAMESPACE 2012-03-09 00:41:25 UTC (rev 625)
@@ -18,6 +18,7 @@
probe_acf,probe_ccf,
probe_nlar,
synth_loglik,
+ do_partrans,
do_rprocess,
do_dprocess,
do_rmeasure,
@@ -46,7 +47,7 @@
pomp,
plot,show,print,coerce,summary,logLik,window,"$",
dprocess,rprocess,rmeasure,dmeasure,init.state,skeleton,
- data.array,obs,coef,"coef<-",time,"time<-",timezero,"timezero<-",
+ data.array,obs,partrans,coef,"coef<-",time,"time<-",timezero,"timezero<-",
simulate,pfilter,
particles,mif,continue,states,trajectory,
pred.mean,pred.var,filter.mean,conv.rec,
Modified: pkg/R/init-state-pomp.R
===================================================================
--- pkg/R/init-state-pomp.R 2012-03-09 00:22:08 UTC (rev 624)
+++ pkg/R/init-state-pomp.R 2012-03-09 00:41:25 UTC (rev 625)
@@ -7,12 +7,6 @@
function (object, params, t0, ...) { # the package algorithms will only use these arguments
if (missing(t0)) t0 <- object at t0
if (missing(params)) params <- coef(object)
- x <- try(
- .Call(do_init_state,object,as.matrix(params),t0),
- silent=FALSE
- )
- if (inherits(x,'try-error'))
- stop("init.state error: error in user ",sQuote("initializer"),call.=FALSE)
- x
+ .Call(do_init_state,object,params,t0)
}
)
Modified: pkg/R/pomp-fun.R
===================================================================
--- pkg/R/pomp-fun.R 2012-03-09 00:22:08 UTC (rev 624)
+++ pkg/R/pomp-fun.R 2012-03-09 00:41:25 UTC (rev 625)
@@ -45,7 +45,7 @@
R.fun=function(...)stop(sQuote(fname)," not specified"),
native.fun=character(0),
PACKAGE=PACKAGE,
- use=1L
+ use=-1L
)
}
retval
@@ -55,13 +55,16 @@
'show',
'pomp.fun',
function (object) {
- if (object at use==1) {
+ use <- object at use
+ if (use==1L) {
show(object at R.fun)
- } else {
+ } else if (use==2L) {
cat("native function ",sQuote(object at native.fun),sep="")
if (length(object at PACKAGE)>0)
cat(", dynamically loaded from ",sQuote(object at PACKAGE),sep="")
cat ("\n")
+ } else {
+ cat("function not specified\n")
}
}
)
@@ -73,13 +76,13 @@
)
get.pomp.fun <- function (object) {
- use <- as.integer(object at use-1)
- if (use==0L) {
+ use <- object at use
+ if (use==1L) {
f <- object at R.fun
- } else if (use==1L) {
+ } else if (use==2L) {
f <- getNativeSymbolInfo(name=object at native.fun,PACKAGE=object at PACKAGE)$address
} else {
- stop("invalid ",sQuote("use")," value")
+ stop("function not specified")
}
- list(f,use)
+ list(f,as.integer(use-1))
}
Modified: pkg/R/pomp-methods.R
===================================================================
--- pkg/R/pomp-methods.R 2012-03-09 00:22:08 UTC (rev 624)
+++ pkg/R/pomp-methods.R 2012-03-09 00:41:25 UTC (rev 625)
@@ -2,18 +2,13 @@
## functions to extract or call the components of a "pomp" object
setGeneric("data.array",function(object,...)standardGeneric("data.array"))
-
setGeneric("obs",function(object,...)standardGeneric("obs"))
-
setGeneric("time<-",function(object,...,value)standardGeneric("time<-"))
-
setGeneric("coef<-",function(object,...,value)standardGeneric("coef<-"))
-
setGeneric("states",function(object,...)standardGeneric("states"))
-
setGeneric("timezero",function(object,...)standardGeneric("timezero"))
-
setGeneric("timezero<-",function(object,...,value)standardGeneric("timezero<-"))
+setGeneric("partrans",function(object,params,dir,...)standardGeneric("partrans"))
## 'coerce' method: allows for coercion of a "pomp" object to a data-frame
setAs(
@@ -39,6 +34,20 @@
as.data.frame.pomp <- function (x, row.names, optional, ...) as(x,"data.frame")
+## parameter transformations
+pomp.transform <- function (object, params, dir = c("forward","inverse"), ...) {
+ if (!object at has.trans) return(params)
+ dir <- match.arg(dir)
+ tfunc <- switch(
+ dir,
+ forward=get.pomp.fun(object at par.trans),
+ inverse=get.pomp.fun(object at par.untrans)
+ )
+ .Call(do_partrans,object,params,tfunc)
+}
+
+setMethod("partrans","pomp",pomp.transform)
+
## a simple method to extract the data array
setMethod(
"data.array",
@@ -167,37 +176,6 @@
}
)
-pomp.transform <- function (object, params, dir = c("forward","inverse")) {
- dir <- match.arg(dir)
- r <- length(dim(params))
- nm <- if (r>0) rownames(params) else names(params)
- tfunc <- switch(
- dir,
- forward=function (x) do.call(object at par.trans,c(list(x),object at userdata)),
- inverse=function (x) do.call(object at par.untrans,c(list(x),object at userdata))
- )
- if (r > 1) {
- retval <- apply(params,seq.int(from=2,to=r),tfunc)
- no.names <- is.null(rownames(retval))
- } else {
- retval <- tfunc(params)
- no.names <- is.null(names(retval))
- }
- if (no.names)
- switch(
- dir,
- forward=stop(
- "invalid ",sQuote("pomp")," object: ",
- sQuote("parameter.transform")," must return a named numeric vector"
- ),
- inverse=stop(
- "invalid ",sQuote("pomp")," object: ",
- sQuote("parameter.inv.transform")," must return a named numeric vector"
- )
- )
- retval
-}
-
## extract the coefficients
setMethod(
"coef",
@@ -205,7 +183,7 @@
function (object, pars, transform = FALSE, ...) {
if (length(object at params)>0) {
if (transform)
- params <- pomp.transform(object,params=object at params,dir="inverse")
+ params <- partrans(object,params=object at params,dir="inverse")
else
params <- object at params
if (missing(pars))
@@ -235,7 +213,7 @@
if (missing(pars)) { ## replace the whole params slot with 'value'
if (length(value)>0) {
if (transform)
- value <- pomp.transform(object,params=value,dir="forward")
+ value <- partrans(object,params=value,dir="forward")
pars <- names(value)
if (is.null(pars)) {
if (transform)
@@ -252,15 +230,12 @@
" names of ",sQuote("value")," are being discarded",
call.=FALSE
)
-### comment these lines out because 'coef(obj,c("a","b")) <- 3' should be legal
-### if (length(pars)!=length(value))
-### stop(sQuote("pars")," and ",sQuote("value")," must be of equal length")
if (length(object at params)==0) { ## no pre-existing 'params' slot
val <- numeric(length(pars))
names(val) <- pars
val[] <- value
if (transform)
- value <- pomp.transform(object,params=val,dir="forward")
+ value <- partrans(object,params=val,dir="forward")
object at params <- value
} else { ## pre-existing params slot
params <- coef(object,transform=transform)
@@ -280,7 +255,7 @@
}
params[pars] <- val
if (transform)
- params <- pomp.transform(object,params=params,dir="forward")
+ params <- partrans(object,params=params,dir="forward")
object at params <- params
}
}
Modified: pkg/R/pomp.R
===================================================================
--- pkg/R/pomp.R 2012-03-09 00:22:08 UTC (rev 624)
+++ pkg/R/pomp.R 2012-03-09 00:41:25 UTC (rev 625)
@@ -22,8 +22,9 @@
paramnames = 'character',
covarnames = 'character',
zeronames = 'character',
- par.trans = 'function',
- par.untrans = 'function',
+ has.trans = 'logical',
+ par.trans = 'pomp.fun',
+ par.untrans = 'pomp.fun',
PACKAGE = 'character',
userdata = 'list'
)
@@ -100,11 +101,6 @@
dmeasure <- mm$dmeasure
}
- if (missing(rmeasure))
- rmeasure <- function(x,t,params,covars,...)stop(sQuote("rmeasure")," not specified")
- if (missing(dmeasure))
- dmeasure <- function(y,x,t,params,log,covars,...)stop(sQuote("dmeasure")," not specified")
-
skeleton.type <- match.arg(skeleton.type)
skelmap.delta.t <- as.numeric(skelmap.delta.t)
if (skelmap.delta.t <= 0)
@@ -131,8 +127,15 @@
call.=TRUE
)
- rmeasure <- pomp.fun(f=rmeasure,PACKAGE=PACKAGE,proto=quote(rmeasure(x,t,params,...)))
- dmeasure <- pomp.fun(f=dmeasure,PACKAGE=PACKAGE,proto=quote(dmeasure(y,x,t,params,log,...)))
+ if (missing(rmeasure))
+ rmeasure <- pomp.fun()
+ else
+ rmeasure <- pomp.fun(f=rmeasure,PACKAGE=PACKAGE,proto=quote(rmeasure(x,t,params,...)))
+
+ if (missing(dmeasure))
+ dmeasure <- pomp.fun()
+ else
+ dmeasure <- pomp.fun(f=dmeasure,PACKAGE=PACKAGE,proto=quote(dmeasure(y,x,t,params,log,...)))
if (!is.function(initializer))
stop(
@@ -244,24 +247,26 @@
}
}
if (has.trans) {
- par.trans <- match.fun(parameter.transform)
- par.untrans <- match.fun(parameter.inv.transform)
- if (!all(c('params','...')%in%names(formals(par.trans))))
- stop(
- "pomp error: ",sQuote("parameter.transform")," must be a function of prototype ",
- sQuote("parameter.transform(params,...)"),
- call.=TRUE
- )
- if (!all(c('params','...')%in%names(formals(par.untrans))))
- stop(
- "pomp error: ",sQuote("parameter.inv.transform")," must be a function of prototype ",
- sQuote("parameter.inv.transform(params,...)"),
- call.=TRUE
- )
+ par.trans <- pomp.fun(
+ f=parameter.transform,
+ PACKAGE=PACKAGE,
+ proto=quote(par.trans(params,...))
+ )
+ par.untrans <- pomp.fun(
+ f=parameter.inv.transform,
+ PACKAGE=PACKAGE,
+ proto=quote(par.untrans(params,...))
+ )
} else {
- par.untrans <- par.trans <- function(params, ...) params
+ par.trans <- pomp.fun()
+ par.untrans <- pomp.fun()
}
-
+ if (
+ has.trans &&
+ par.trans at use<0 &&
+ par.untrans at use<0
+ )
+ has.trans <- FALSE
new(
'pomp',
@@ -283,6 +288,7 @@
paramnames = paramnames,
covarnames = covarnames,
zeronames = zeronames,
+ has.trans = has.trans,
par.trans = par.trans,
par.untrans = par.untrans,
PACKAGE = PACKAGE,
@@ -541,8 +547,8 @@
stop("pomp error: if ",sQuote("parameter.transform")," is supplied, then " ,
sQuote("parameter.inv.transform")," must also be supplied")
} else {
- par.trans <- match.fun(parameter.transform)
- par.untrans <- match.fun(parameter.inv.transform)
+ par.trans <- parameter.transform
+ par.untrans <- parameter.inv.transform
}
}
Modified: pkg/data/blowflies.rda
===================================================================
(Binary files differ)
Modified: pkg/data/dacca.rda
===================================================================
(Binary files differ)
Modified: pkg/data/euler.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/gillespie.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/gompertz.rda
===================================================================
(Binary files differ)
Modified: pkg/data/ou2.rda
===================================================================
(Binary files differ)
Modified: pkg/data/ricker.rda
===================================================================
(Binary files differ)
Modified: pkg/data/rw2.rda
===================================================================
(Binary files differ)
Modified: pkg/data/verhulst.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/NEWS
===================================================================
--- pkg/inst/NEWS 2012-03-09 00:22:08 UTC (rev 624)
+++ pkg/inst/NEWS 2012-03-09 00:41:25 UTC (rev 625)
@@ -1,4 +1,10 @@
NEWS
+0.40-8
+ o A new method 'partrans' allows transformation of vectors or matrices of parameters.
+ The parameter transformations have been pushed into C for speed.
+ It is possible to specify native parameter transformation routines in addition to C functions, but this last facility has not yet been extensively tested.
+ A new slot 'has.trans' has been introduced into the 'pomp' class, and the types of slots 'par.trans' and 'par.untrans' have changed.
+
0.40-7
o 'parmat' now gracefully handles the case when 'params' is already a matrix.
Modified: pkg/inst/doc/advanced_topics_in_pomp.Rnw
===================================================================
--- pkg/inst/doc/advanced_topics_in_pomp.Rnw 2012-03-09 00:22:08 UTC (rev 624)
+++ pkg/inst/doc/advanced_topics_in_pomp.Rnw 2012-03-09 00:41:25 UTC (rev 625)
@@ -518,7 +518,7 @@
The initial states and parameters must be matrices, and they are checked for commensurability.
The method returns a rank-3 array containing simulated state trajectories, sampled at the times specified.
<<>>=
-x <- rprocess(ou2,xstart=x0,times=time(ou2,t0=T),params=as.matrix(true.p))
+x <- rprocess(ou2,xstart=x0,times=time(ou2,t0=T),params=true.p)
dim(x)
x[,,1:5]
@
@@ -530,7 +530,7 @@
The \code{rmeasure} method gives access to the measurement model simulator:
<<>>=
x <- x[,,-1,drop=F]
-y <- rmeasure(ou2,x=x,times=time(ou2),params=as.matrix(true.p))
+y <- rmeasure(ou2,x=x,times=time(ou2),params=true.p)
dim(y)
y[,,1:5]
@
@@ -539,12 +539,12 @@
The \code{dmeasure} and \code{dprocess} methods give access to the measurement and process model densities, respectively.
<<>>=
-fp <- dprocess(ou2,x=x,times=time(ou2),params=as.matrix(true.p))
+fp <- dprocess(ou2,x=x,times=time(ou2),params=true.p)
dim(fp)
fp[,36:40]
@
<<>>=
-fm <- dmeasure(ou2,y=y[,1,],x=x,times=time(ou2),params=as.matrix(true.p))
+fm <- dmeasure(ou2,y=y[,1,],x=x,times=time(ou2),params=true.p)
dim(fm)
fm[,36:40]
@
@@ -554,7 +554,7 @@
There are a number of example \code{pomp} objects included with the package.
These can be found by running
-<<all-data-loaable,eval=F>>=
+<<all-data-loadable,eval=F>>=
data(package="pomp")
@
The \R\ scripts that generated these are included in the \code{data-R} directory of the installed package.
Modified: pkg/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/gompertz-multi-mif.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/gompertz-trajmatch.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/nlf-block-boot.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/nlf-boot.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/nlf-fit-from-truth.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/nlf-fits.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/nlf-lag-tests.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/nlf-multi-short.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/ricker-mif.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/ricker-probe-match.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/include/pomp.h
===================================================================
--- pkg/inst/include/pomp.h 2012-03-09 00:22:08 UTC (rev 624)
+++ pkg/inst/include/pomp.h 2012-03-09 00:41:25 UTC (rev 625)
@@ -37,6 +37,9 @@
// facility for computing evaluating a basis of periodic bsplines
void periodic_bspline_basis_eval (double x, double period, int degree, int nbasis, double *y);
+// Prototype for parameter transformation function.
+typedef void pomp_transform_fn (double *pt, double *p, int *parindex);
+
// Prototype for stochastic simulation algorithm reaction-rate function, as used by "gillespie.sim":
typedef double pomp_ssa_rate_fn(int j, double t, const double *x, const double *p,
int *stateindex, int *parindex, int *covindex,
Modified: pkg/man/pomp-class.Rd
===================================================================
--- pkg/man/pomp-class.Rd 2012-03-09 00:22:08 UTC (rev 624)
+++ pkg/man/pomp-class.Rd 2012-03-09 00:41:25 UTC (rev 625)
@@ -62,8 +62,9 @@
\item{zeronames}{
Names of accumulator variables.
}
- \item{par.trans, par.untrans}{
- The forward and inverse parameter transformations.
+ \item{has.trans, par.trans, par.untrans}{
+ \code{has.trans=TRUE} if there is a non-trivial parameter transformation.
+ In this case, the forward and inverse parameter transformations are stored in \code{par.trans} and \code{par.untrans}.
}
\item{PACKAGE}{
Character variable giving the name of the dynamically loadable library containing native routines.
Modified: pkg/src/dmeasure.c
===================================================================
--- pkg/src/dmeasure.c 2012-03-09 00:22:08 UTC (rev 624)
+++ pkg/src/dmeasure.c 2012-03-09 00:41:25 UTC (rev 625)
@@ -17,16 +17,20 @@
double t, *lp, *xp, *pp, *yp;
int nvar = ndim[0];
int npar = ndim[1];
- int nrep = ndim[2];
- int ntimes = ndim[3];
- int covlen = ndim[4];
- int covdim = ndim[5];
- int nobs = ndim[6];
+ int nrepp = ndim[2];
+ int nrepx = ndim[3];
+ int ntimes = ndim[4];
+ int covlen = ndim[5];
+ int covdim = ndim[6];
+ int nobs = ndim[7];
double covar_fn[covdim];
+ int nrep;
int k, p;
// set up the covariate table
struct lookup_table covariate_table = {covlen, covdim, 0, time_table, covar_table};
+
+ nrep = (nrepp > nrepx) ? nrepp : nrepx;
for (k = 0; k < ntimes; k++) { // loop over times
@@ -42,8 +46,8 @@
for (p = 0; p < nrep; p++) { // loop over replicates
lp = &lik[p+nrep*k];
- xp = &x[nvar*(p+nrep*k)];
- pp = ¶ms[npar*p];
+ xp = &x[nvar*((p%nrepx)+nrepx*k)];
+ pp = ¶ms[npar*(p%nrepp)];
(*f)(lp,yp,xp,pp,*give_log,obsindex,stateindex,parindex,covindex,covdim,covar_fn,t);
@@ -107,10 +111,9 @@
SEXP do_dmeasure (SEXP object, SEXP y, SEXP x, SEXP times, SEXP params, SEXP log, SEXP fun)
{
int nprotect = 0;
- int *dim, nvars, npars, nreps, ntimes, covlen, covdim, nobs;
+ int *dim, nvars, npars, nrepsp, nrepsx, nreps, ntimes, covlen, covdim, nobs;
int ndim[7];
SEXP F, fn;
- SEXP dimP, dimX, dimY;
SEXP tcovar, covar;
SEXP statenames, paramnames, covarnames, obsnames;
SEXP sindex, pindex, cindex, oindex;
@@ -125,27 +128,29 @@
if (ntimes < 1)
error("dmeasure error: no work to do");
- PROTECT(dimY = GET_DIM(y)); nprotect++;
- if ((isNull(dimY)) || (length(dimY)!=2))
- error("dmeasure error: 'y' must be a rank-2 array");
- dim = INTEGER(dimY); nobs = dim[0];
+ PROTECT(y = as_matrix(y)); nprotect++;
+ dim = INTEGER(GET_DIM(y));
+ nobs = dim[0];
+
if (ntimes != dim[1])
error("dmeasure error: length of 'times' and 2nd dimension of 'y' do not agree");
- PROTECT(dimX = GET_DIM(x)); nprotect++;
- if ((isNull(dimX)) || (length(dimX)!=3))
- error("dmeasure error: 'x' must be a rank-3 array");
- dim = INTEGER(dimX); nvars = dim[0]; nreps = dim[1];
+ PROTECT(x = as_state_array(x)); nprotect++;
+ dim = INTEGER(GET_DIM(x));
+ nvars = dim[0]; nrepsx = dim[1];
+
if (ntimes != dim[2])
error("dmeasure error: length of 'times' and 3rd dimension of 'x' do not agree");
- PROTECT(dimP = GET_DIM(params)); nprotect++;
- if ((isNull(dimP)) || (length(dimP)!=2))
- error("dmeasure error: 'params' must be a rank-2 array");
- dim = INTEGER(dimP); npars = dim[0];
- if (nreps != dim[1])
- error("dmeasure error: 2nd dimensions of 'params' and 'x' do not agree");
+ PROTECT(params = as_matrix(params)); nprotect++;
+ dim = INTEGER(GET_DIM(params));
+ npars = dim[0]; nrepsp = dim[1];
+ nreps = (nrepsp > nrepsx) ? nrepsp : nrepsx;
+
+ if ((nreps % nrepsp != 0) || (nreps % nrepsx != 0))
+ error("dmeasure error: larger number of replicates is not a multiple of smaller");
+
PROTECT(tcovar = GET_SLOT(object,install("tcovar"))); nprotect++;
PROTECT(covar = GET_SLOT(object,install("covar"))); nprotect++;
PROTECT(obsnames = GET_SLOT(object,install("obsnames"))); nprotect++;
@@ -237,8 +242,8 @@
cidx = 0;
}
- ndim[0] = nvars; ndim[1] = npars; ndim[2] = nreps; ndim[3] = ntimes;
- ndim[4] = covlen; ndim[5] = covdim; ndim[6] = nobs;
+ ndim[0] = nvars; ndim[1] = npars; ndim[2] = nrepsp; ndim[3] = nrepsx; ndim[4] = ntimes;
+ ndim[5] = covlen; ndim[6] = covdim; ndim[7] = nobs;
dens_meas(ff,REAL(F),REAL(y),REAL(x),REAL(times),REAL(params),INTEGER(log),ndim,
oidx,sidx,pidx,cidx,REAL(tcovar),REAL(covar));
Modified: pkg/src/dprocess.c
===================================================================
--- pkg/src/dprocess.c 2012-03-09 00:22:08 UTC (rev 624)
+++ pkg/src/dprocess.c 2012-03-09 00:41:25 UTC (rev 625)
@@ -5,37 +5,74 @@
#include <Rdefines.h>
#include <Rinternals.h>
+#include "pomp_internal.h"
+
SEXP do_dprocess (SEXP object, SEXP x, SEXP times, SEXP params, SEXP log)
{
int nprotect = 0;
- int *xdim, nreps, ntimes;
+ int *xdim, npars, nvars, nreps, nrepsx, ntimes;
SEXP X, fn, fcall, rho;
- SEXP dimP, dimX, dimF;
+ SEXP dimP, dimF;
+
ntimes = length(times);
- if (ntimes < 2) {
- UNPROTECT(nprotect);
- error("dprocess error: no transitions, no work to do");
+ if (ntimes < 2)
+ error("dprocess error: length(times)==0: no transitions, no work to do");
+
+ PROTECT(x = as_state_array(x)); nprotect++;
+ xdim = INTEGER(GET_DIM(x));
+ nvars = xdim[0]; nrepsx = xdim[1];
+ if (ntimes != xdim[2])
+ error("dprocess error: length of 'times' and 3rd dimension of 'x' do not agree");
+
+ PROTECT(params = as_matrix(params)); nprotect++;
+ xdim = INTEGER(GET_DIM(params));
+ npars = xdim[0]; nreps = xdim[1];
+
+ if (nrepsx > nreps) { // more states than parameters
+ if (nrepsx % nreps != 0) {
+ error("dprocess error: larger number of replicates is not a multiple of smaller");
+ } else {
+ SEXP copy;
+ double *src, *tgt;
+ int dims[2];
+ int j, k;
+ dims[0] = npars; dims[1] = nrepsx;
+ PROTECT(copy = duplicate(params)); nprotect++;
+ PROTECT(params = makearray(2,dims)); nprotect++;
+ setrownames(params,GET_ROWNAMES(GET_DIMNAMES(copy)),2);
+ src = REAL(copy);
+ tgt = REAL(params);
+ for (j = 0; j < nrepsx; j++) {
+ for (k = 0; k < npars; k++, tgt++) {
+ *tgt = src[k+npars*(j%nreps)];
+ }
+ }
+ }
+ nreps = nrepsx;
+ } else if (nrepsx < nreps) { // more parameters than states
+ if (nreps % nrepsx != 0) {
+ error("dprocess error: larger number of replicates is not a multiple of smaller");
+ } else {
+ SEXP copy;
+ double *src, *tgt;
+ int dims[3];
+ int i, j, k;
+ dims[0] = nvars; dims[1] = nreps; dims[2] = ntimes;
+ PROTECT(copy = duplicate(x)); nprotect++;
+ PROTECT(x = makearray(3,dims)); nprotect++;
+ setrownames(x,GET_ROWNAMES(GET_DIMNAMES(copy)),3);
+ src = REAL(copy);
+ tgt = REAL(x);
+ for (i = 0; i < ntimes; i++) {
+ for (j = 0; j < nreps; j++) {
+ for (k = 0; k < nvars; k++, tgt++) {
+ *tgt = src[k+nvars*((j%nrepsx)+nrepsx*i)];
+ }
+ }
+ }
+ }
}
- PROTECT(dimX = GET_DIM(x)); nprotect++;
- if ((isNull(dimX)) || (length(dimX)!=3)) {
- UNPROTECT(nprotect);
- error("dprocess error: 'x' must be a rank-3 array");
- }
- PROTECT(dimP = GET_DIM(params)); nprotect++;
- if ((isNull(dimP)) || (length(dimP)!=2)) {
- UNPROTECT(nprotect);
- error("dprocess error: 'params' must be a rank-2 array");
- }
- xdim = INTEGER(dimX);
- nreps = xdim[1];
- if (nreps != INTEGER(dimP)[1]) {
- UNPROTECT(nprotect);
- error("dprocess error: number of columns of 'params' and 'x' do not agree");
- }
- if (ntimes != INTEGER(dimX)[2]) {
- UNPROTECT(nprotect);
- error("rprocess error: length of 'times' and 3rd dimension of 'x' do not agree");
- }
+
// extract the process function
PROTECT(fn = GET_SLOT(object,install("dprocess"))); nprotect++;
// construct the call
@@ -61,8 +98,11 @@
PROTECT(fcall = LCONS(x,fcall)); nprotect++;
SET_TAG(fcall,install("x"));
PROTECT(fcall = LCONS(fn,fcall)); nprotect++;
+
PROTECT(rho = (CLOENV(fn))); nprotect++; // environment of the function
+
PROTECT(X = eval(fcall,rho)); nprotect++; // do the call
+
PROTECT(dimF = GET_DIM(X)); nprotect++;
if ((isNull(dimF)) || (length(dimF) != 2)) {
UNPROTECT(nprotect);
@@ -73,6 +113,7 @@
UNPROTECT(nprotect);
error("dprocess error: user 'dprocess' must return a %d x %d array",nreps,ntimes-1);
}
+
UNPROTECT(nprotect);
return X;
}
Modified: pkg/src/initstate.c
===================================================================
--- pkg/src/initstate.c 2012-03-09 00:22:08 UTC (rev 624)
+++ pkg/src/initstate.c 2012-03-09 00:41:25 UTC (rev 625)
@@ -19,6 +19,7 @@
int xdim[2], j, k;
double *p, *pp, *xp, *xpp;
+ PROTECT(params = as_matrix(params)); nprotect++;
dim = INTEGER(GET_DIM(params));
npar = dim[0]; nrep = dim[1];
p = REAL(params);
Added: pkg/src/partrans.c
===================================================================
--- pkg/src/partrans.c (rev 0)
+++ pkg/src/partrans.c 2012-03-09 00:41:25 UTC (rev 625)
@@ -0,0 +1,136 @@
+// dear emacs, please treat this as -*- C++ -*-
+
+#include <R.h>
+#include <Rmath.h>
+#include <Rdefines.h>
+#include <Rinternals.h>
+
+#include "pomp_internal.h"
+
+static SEXP pomp_R_transform_params (SEXP object, SEXP params, SEXP fun) {
+ int nprotect = 0;
+ SEXP Dim, nm;
+ SEXP params_transformed, par, fcall, rho, ans;
+ double *p = 0, *pp = 0, *pt = 0, *ps = 0;
+ int nreps, npar1, npar2;
+ int k;
+
+ PROTECT(Dim = GET_DIM(params)); nprotect++;
+ if (isNull(Dim)) { // a single vector
+ nreps = 1;
+ } else { // a parameter matrix
+ int *dim = INTEGER(Dim);
+ npar1 = dim[0]; nreps = dim[1];
+ }
+
+ if (nreps > 1) { // matrix case
+ PROTECT(par = NEW_NUMERIC(npar1)); nprotect++;
+ SET_NAMES(par,GET_ROWNAMES(GET_DIMNAMES(params)));
+ }
+
+ // set up the function call
+ PROTECT(rho = (CLOENV(fun))); nprotect++;
+ PROTECT(fcall = VectorToPairList(GET_SLOT(object,install("userdata")))); nprotect++;
+ if (nreps > 1) { // matrix case
+ PROTECT(fcall = LCONS(par,fcall)); nprotect++;
+ } else { // vector case
+ PROTECT(fcall = LCONS(params,fcall)); nprotect++;
+ }
+ SET_TAG(fcall,install("params"));
+ PROTECT(fcall = LCONS(fun,fcall)); nprotect++;
+
+ if (nreps > 1) {
+ p = REAL(params);
+ pp = REAL(par);
+
+ memcpy(pp,p,npar1*sizeof(double));
+
+ PROTECT(ans = eval(fcall,rho)); nprotect++;
+
+ PROTECT(nm = GET_NAMES(ans)); nprotect++;
+ if (isNull(nm))
+ error("user transformation functions must return a named numeric vector");
+
+ {
+ int ndim[2];
+ npar2 = LENGTH(ans);
+ ndim[0] = npar2; ndim[1] = nreps;
+ PROTECT(params_transformed = makearray(2,ndim)); nprotect++;
+ setrownames(params_transformed,nm,2);
+ }
+
+ ps = REAL(AS_NUMERIC(ans));
+ pt = REAL(params_transformed);
+ memcpy(pt,ps,npar2*sizeof(double));
+
+ p += npar1;
+ pt += npar2;
+ for (k = 1; k < nreps; k++, p += npar1, pt += npar2) {
+ memcpy(pp,p,npar1*sizeof(double));
+ // PROTECT(ans = AS_NUMERIC(eval(fcall,rho)));
+ // ps = REAL(ans);
+ ps = REAL(AS_NUMERIC(eval(fcall,rho)));
+ memcpy(pt,ps,npar2*sizeof(double));
+ // UNPROTECT(1);
+ }
+
+ } else {
+
+ PROTECT(params_transformed = eval(fcall,rho)); nprotect++;
+ if (isNull(GET_NAMES(params_transformed)))
+ error("user transformation functions must return a named numeric vector");
+
+ }
+
+ UNPROTECT(nprotect);
+ return params_transformed;
+}
+
+SEXP do_partrans (SEXP object, SEXP params, SEXP fun)
+{
+ int nprotect = 0;
+ SEXP ptrans, fn;
+ int use, drop;
+
+ PROTECT(fn = VECTOR_ELT(fun,0)); nprotect++;
+ use = INTEGER(VECTOR_ELT(fun,1))[0];
+
+ switch (use) {
+ case 0: // use user-supplied R function
+ PROTECT(ptrans = pomp_R_transform_params(object,params,fn)); nprotect++;
+ break;
+ case 1: // use native routine
+ {
+ pomp_transform_fn *ff = (pomp_transform_fn *) R_ExternalPtrAddr(fn);
+ SEXP paramnames, pindex;
+ int *idx, npar, nrep;
+ double *ps, *pt;
+ int k;
+
+ idx = INTEGER(GET_DIM(params));
+ npar = idx[0]; nrep = idx[1];
+
+ PROTECT(paramnames = GET_SLOT(object,install("paramnames"))); nprotect++;
+ if (LENGTH(paramnames) > 0) {
+ PROTECT(pindex = matchnames(GET_ROWNAMES(GET_DIMNAMES(params)),paramnames)); nprotect++;
+ idx = INTEGER(pindex);
+ } else {
+ idx = 0;
+ }
+
+ PROTECT(ptrans = duplicate(params)); nprotect++;
+
+ for (k = 0, ps = REAL(params), pt = REAL(ptrans); k < nrep; k++, ps += npar, pt += npar) {
+ R_CheckUserInterrupt();
+ (*ff)(pt,ps,idx);
+ }
+
+ }
+ break;
+ default:
+ error("unrecognized 'use' slot in 'partrans'");
+ }
+
+ UNPROTECT(nprotect);
+ return ptrans;
+}
Modified: pkg/src/pomp.h
===================================================================
--- pkg/src/pomp.h 2012-03-09 00:22:08 UTC (rev 624)
+++ pkg/src/pomp.h 2012-03-09 00:41:25 UTC (rev 625)
@@ -37,6 +37,9 @@
// facility for computing evaluating a basis of periodic bsplines
void periodic_bspline_basis_eval (double x, double period, int degree, int nbasis, double *y);
+// Prototype for parameter transformation function.
+typedef void pomp_transform_fn (double *pt, double *p, int *parindex);
+
// Prototype for stochastic simulation algorithm reaction-rate function, as used by "gillespie.sim":
typedef double pomp_ssa_rate_fn(int j, double t, const double *x, const double *p,
int *stateindex, int *parindex, int *covindex,
Modified: pkg/src/pomp_internal.h
===================================================================
--- pkg/src/pomp_internal.h 2012-03-09 00:22:08 UTC (rev 624)
+++ pkg/src/pomp_internal.h 2012-03-09 00:41:25 UTC (rev 625)
@@ -122,6 +122,95 @@
UNPROTECT(nprotect);
}
+static R_INLINE SEXP as_matrix (SEXP x) {
+ int nprotect = 0;
+ SEXP dim, names;
+ int *xdim, nrow, ncol;
+ PROTECT(dim = GET_DIM(x)); nprotect++;
+ if (isNull(dim)) {
+ PROTECT(x = duplicate(x)); nprotect++;
+ PROTECT(names = GET_NAMES(x)); nprotect++;
+ dim = NEW_INTEGER(2);
+ xdim = INTEGER(dim);
+ xdim[0] = LENGTH(x); xdim[1] = 1;
+ SET_DIM(x,dim);
+ SET_NAMES(x,R_NilValue);
+ setrownames(x,names,2);
+ } else if (LENGTH(dim) == 1) {
+ PROTECT(x = duplicate(x)); nprotect++;
+ PROTECT(names = GET_ROWNAMES(GET_DIMNAMES(x))); nprotect++;
+ dim = NEW_INTEGER(2);
+ xdim = INTEGER(dim);
+ xdim[0] = LENGTH(x); xdim[1] = 1;
+ SET_DIM(x,dim);
+ SET_NAMES(x,R_NilValue);
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 625
More information about the pomp-commits
mailing list