[Pomp-commits] r322 - in pkg: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Sep 24 15:55:35 CEST 2010
Author: kingaa
Date: 2010-09-24 15:55:35 +0200 (Fri, 24 Sep 2010)
New Revision: 322
Added:
pkg/src/probe.c
Modified:
pkg/R/probe-match.R
pkg/R/probe.R
Log:
- put the most computationally-intensive bits of the probe calculations into C
Modified: pkg/R/probe-match.R
===================================================================
--- pkg/R/probe-match.R 2010-09-24 10:31:20 UTC (rev 321)
+++ pkg/R/probe-match.R 2010-09-24 13:55:35 UTC (rev 322)
@@ -38,17 +38,29 @@
params[est] <- par
- simval <- apply.probe.sim(object,probes=probes,params=params,nsim=nsim,seed=seed) # apply probes to model simulations
+ ## apply probes to model simulations
+ simval <- .Call(
+ apply_probe_sim,
+ object=object,
+ nsim=nsim,
+ params=params,
+ seed=seed,
+ probes=probes,
+ pnames=names(datval)
+ )
## compute a measure of the discrepancies between simulations and data
sim.means <- colMeans(simval)
simval <- sweep(simval,2,sim.means)
discrep <- ((datval-sim.means)^2)/colMeans(simval^2)
-
+ if ((length(weights)>1) && (length(weights)!=length(discrep)))
+ stop(length(discrep)," probes have been computed, but ",length(weights)," have been supplied")
if (!all(is.finite(discrep))) {
mismatch <- fail.value
+ } else if (length(weights)>1) {
+ mismatch <- sum(discrep*weights)/sum(weights)
} else {
- mismatch <- sum(discrep*weights)/sum(weights)
+ mismatch <- sum(discrep)
}
mismatch
@@ -82,15 +94,17 @@
if (!is.list(probes)) probes <- list(probes)
if (!all(sapply(probes,is.function)))
stop(sQuote("probes")," must be a function or a list of functions")
+ if (!all(sapply(probes,function(f)length(formals(f))==1)))
+ stop("each probe must be a function of a single argument")
- if (missing(weights)) weights <- rep(1,length(probes))
+ if (missing(weights)) weights <- 1
method <- match.arg(method)
params <- start
guess <- params[par.index]
- datval <- apply.probe.data(object,probes=probes) # apply probes to data
+ datval <- .Call(apply_probe_data,object,probes) # apply probes to data
if (eval.only) {
val <- probe.mismatch(
Modified: pkg/R/probe.R
===================================================================
--- pkg/R/probe.R 2010-09-24 10:31:20 UTC (rev 321)
+++ pkg/R/probe.R 2010-09-24 13:55:35 UTC (rev 322)
@@ -26,57 +26,6 @@
setGeneric("probe",function(object,probes,...)standardGeneric("probe"))
-apply.probe.data <- function (object, probes) {
- val <- vector(mode="list",length=length(probes))
- names(val) <- names(probes)
- yy <- data.array(object)
- for (p in seq_along(probes)) {
- f <- probes[[p]]
- if (length(formals(f))>1)
- stop("each element of ",sQuote("probes")," must be a function of a single argument")
- vv <- f(yy)
- val[[p]] <- vv
- }
- do.call(c,val)
-}
-
-apply.probe.sim <- function (object, probes, params, nsim, seed) {
- if (!is.numeric(nsim)||(length(nsim)>1)||(nsim<=0))
- stop(sQuote("nsim")," must be a positive integer")
- nsim <- as.integer(nsim)
-
- y <- simulate(
- object,
- nsim=nsim,
- seed=seed,
- params=params,
- times=time(object,t0=TRUE),
- obs=TRUE
- )[,,-1,drop=FALSE]
- nmy <- dimnames(y)
- dy <- dim(y)
- yy <- array(dim=dy[c(1,3)],dimnames=nmy[c(1,3)])
-
- val <- vector(mode="list",length=length(probes))
- names(val) <- names(probes)
- for (s in seq_len(nsim)) {
- yy[,] <- y[,s,]
- for (p in seq_along(probes)) {
- f <- probes[[p]]
- vv <- f(yy)
- if (s==1) {
- val[[p]] <- array(dim=c(nsim,length(vv)))
- } else {
- if (length(vv)!=ncol(val[[p]]))
- if (length(vv)!=ncol(val))
- stop("the size of the vector returned by a probe must not vary")
- }
- val[[p]][s,] <- vv
- }
- }
- do.call(cbind,val)
-}
-
setMethod(
"probe",
signature(object="pomp"),
@@ -93,16 +42,10 @@
}
## apply probes to data
- datval <- apply.probe.data(object,probes=probes)
+ datval <- .Call(apply_probe_data,object,probes)
## apply probes to model simulations
- simval <- apply.probe.sim(
- object,
- probes=probes,
- params=coef(object),
- nsim=nsim,
- seed=seed
- )
-
+ simval <- .Call(apply_probe_sim,object,nsim,coef(object),seed,probes,names(datval))
+
nprobes <- length(datval)
pvals <- numeric(nprobes)
names(pvals) <- names(datval)
Added: pkg/src/probe.c
===================================================================
--- pkg/src/probe.c (rev 0)
+++ pkg/src/probe.c 2010-09-24 13:55:35 UTC (rev 322)
@@ -0,0 +1,119 @@
+// -*- mode: C++; -*-
+#include "pomp_internal.h"
+
+SEXP apply_probe_data (SEXP object, SEXP probes) {
+ int nprotect = 0;
+ SEXP retval, data, vals;
+ int nprobe;
+ int i;
+
+ nprobe = LENGTH(probes);
+ PROTECT(data = GET_SLOT(object,install("data"))); nprotect++;
+ PROTECT(vals = NEW_LIST(nprobe)); nprotect++;
+ SET_NAMES(vals,GET_NAMES(probes));
+
+ for (i = 0; i < nprobe; i++) {
+ SET_ELEMENT(vals,i,eval(lang2(VECTOR_ELT(probes,i),data),CLOENV(VECTOR_ELT(probes,i))));
+ if (!IS_NUMERIC(VECTOR_ELT(vals,i))) {
+ UNPROTECT(nprotect);
+ error("probe %ld returns a non-numeric result",i);
+ }
+ }
+
+ PROTECT(retval = eval(LCONS(install("c"),VectorToPairList(vals)),R_BaseEnv)); nprotect++;
+
+ UNPROTECT(nprotect);
+ return retval;
+}
+
+SEXP apply_probe_sim (SEXP object, SEXP nsim, SEXP params, SEXP seed, SEXP probes, SEXP names) {
+ int nprotect = 0;
+ SEXP y, obs, call;
+ SEXP retval, val, valnames, x;
+ int nprobe, nsims, nvars, ntimes, nvals;
+ int xdim[2];
+ double *xp, *yp;
+ int p, s, i, j, k, len0, len = 0;
+
+ PROTECT(nsim = AS_INTEGER(nsim)); nprotect++;
+ if ((LENGTH(nsim)>1) || (INTEGER(nsim)[0]<=0))
+ error("'nsim' must be a positive integer");
+
+ // 'names' holds the names of the probe values
+ // we get these from a previous call to 'apply_probe_data'
+ nprobe = LENGTH(probes);
+ nvals = LENGTH(names);
+ PROTECT(names = AS_CHARACTER(names)); nprotect++;
+
+ // call 'simulate' to get simulated data sets
+ PROTECT(obs = NEW_LOGICAL(1)); nprotect++;
+ LOGICAL(obs)[0] = 1; // we set obs=TRUE
+ PROTECT(call = LCONS(obs,R_NilValue)); nprotect++;
+ SET_TAG(call,install("obs"));
+ PROTECT(call = LCONS(params,call)); nprotect++;
+ SET_TAG(call,install("params"));
+ PROTECT(call = LCONS(seed,call)); nprotect++;
+ SET_TAG(call,install("seed"));
+ PROTECT(call = LCONS(nsim,call)); nprotect++;
+ SET_TAG(call,install("nsim"));
+ PROTECT(call = LCONS(object,call)); nprotect++;
+ SET_TAG(call,install("object"));
+ PROTECT(call = LCONS(install("simulate"),call)); nprotect++;
+ PROTECT(y = eval(call,R_GlobalEnv)); nprotect++;
+
+ nvars = INTEGER(GET_DIM(y))[0];
+ nsims = INTEGER(GET_DIM(y))[1];
+ ntimes = INTEGER(GET_DIM(y))[2]; // recall that 'simulate' returns a value for time zero
+
+ // set up temporary storage
+ xdim[0] = nvars; xdim[1] = ntimes-1;
+ PROTECT(x = makearray(2,xdim)); nprotect++;
+ setrownames(x,GET_ROWNAMES(GET_DIMNAMES(y)),2);
+
+ // set up matrix to hold results
+ xdim[0] = nsims; xdim[1] = nvals;
+ PROTECT(retval = makearray(2,xdim)); nprotect++;
+ PROTECT(valnames = NEW_LIST(2)); nprotect++;
+ SET_ELEMENT(valnames,1,names); // set column names
+ SET_DIMNAMES(retval,valnames);
+
+ for (p = 0, k = 0; p < nprobe; p++, k += len) { // loop over probes
+
+ for (s = 0; s < nsims; s++) { // loop over simulations
+
+ // copy the data from y[,s,-1] to x[,]
+ xp = REAL(x);
+ yp = REAL(y)+nvars*(s+nsims);
+ for (j = 1; j < ntimes; j++, yp += nvars*nsims) {
+ for (i = 0; i < nvars; i++, xp++) *xp = yp[i];
+ }
+
+ // evaluate the probe on the simulated data
+ PROTECT(val = eval(lang2(VECTOR_ELT(probes,p),x),CLOENV(VECTOR_ELT(probes,p)))); nprotect++;
+ if (!IS_NUMERIC(val)) {
+ UNPROTECT(nprotect);
+ error("probe %ld returns a non-numeric result",p);
+ }
+
+ len = LENGTH(val);
+ if (s == 0)
+ len0 = len;
+ else if (len != len0) {
+ UNPROTECT(nprotect);
+ error("variable-sized results returned by probe %ld",p);
+ }
+ if (k+len > nvals)
+ error("probes return different number of values on different datasets");
+
+ xp = REAL(retval); yp = REAL(val);
+ for (i = 0; i < len; i++) xp[s+nsims*(i+k)] = yp[i];
+
+ }
+
+ }
+ if (k != nvals)
+ error("probes return different number of values on different datasets");
+
+ UNPROTECT(nprotect);
+ return retval;
+}
Property changes on: pkg/src/probe.c
___________________________________________________________________
Added: svn:eol-style
+ native
More information about the pomp-commits
mailing list