[Pomp-commits] r360 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Oct 4 21:53:57 CEST 2010
Author: kingaa
Date: 2010-10-04 21:53:57 +0200 (Mon, 04 Oct 2010)
New Revision: 360
Modified:
pkg/R/trajectory-pomp.R
Log:
- put guts of trajectory calculation into the named function 'trajectory.internal'
Modified: pkg/R/trajectory-pomp.R
===================================================================
--- pkg/R/trajectory-pomp.R 2010-10-04 16:10:27 UTC (rev 359)
+++ pkg/R/trajectory-pomp.R 2010-10-04 19:53:57 UTC (rev 360)
@@ -1,97 +1,95 @@
-trajectory <- function (object, params, times, ...)
+trajectory <- function (object, params, times, t0, ...)
stop("function ",sQuote("trajectory")," is undefined for objects of class ",sQuote(class(object)))
setGeneric('trajectory')
-setMethod(
- "trajectory",
- "pomp",
- function (object, params, times, ...) {
- if (missing(params)) {
- params <- coef(object)
- if (length(params)==0) {
- stop("trajectory error: ",sQuote("params")," must be supplied",call.=FALSE)
- }
- }
- nrep <- NCOL(params)
- if (is.null(dim(params))) {
- params <- matrix(
- params,
- nrow=length(params),
- ncol=nrep,
- dimnames=list(
- names(params),
- NULL
- )
- )
- }
- paramnames <- rownames(params)
- if (is.null(paramnames))
- stop("pfilter error: ",sQuote("params")," must have rownames",call.=FALSE)
- params <- as.matrix(params)
+trajectory.internal <- function (object, params, times, t0, ...) {
+ if (missing(params)) {
+ params <- coef(object)
+ if (length(params)==0) {
+ stop("trajectory error: ",sQuote("params")," must be supplied",call.=FALSE)
+ }
+ }
+ nrep <- NCOL(params)
+ if (is.null(dim(params))) {
+ params <- matrix(
+ params,
+ nrow=length(params),
+ ncol=nrep,
+ dimnames=list(
+ names(params),
+ NULL
+ )
+ )
+ }
+ paramnames <- rownames(params)
+ if (is.null(paramnames))
+ stop("pfilter error: ",sQuote("params")," must have rownames",call.=FALSE)
+ params <- as.matrix(params)
- if (missing(times))
- times <- time(object,t0=TRUE)
- else
- times <- as.numeric(times)
+ if (missing(times))
+ times <- time(object,t0=TRUE)
+ else
+ times <- as.numeric(times)
- tm <- times[1]
- x0 <- init.state(object,params=params,t0=tm)
- nm <- rownames(x0)
- dim(x0) <- c(nrow(x0),nrep,1)
- dimnames(x0) <- list(nm,NULL,NULL)
-
- x <- array(
- dim=c(nrow(x0),nrep,length(times)),
- dimnames=list(rownames(x0),NULL,NULL)
- )
- switch(
- object at skeleton.type,
- map={ # iterate the map
- for (k in seq_along(times)) {
- if (tm < times[k]) {
- while (tm < times[k]) {
- x0[,,] <- skeleton(object,x=x0,t=tm,params=params)
- tm <- tm+1
- }
- }
- x[,,k] <- x0
- }
- },
- vectorfield={ # integrate the vectorfield
- for (j in seq_len(nrep)) {
- X <- try(
- deSolve::lsoda(
- y=x0[,j,1],
- times=times,
- func=function(t,y,parms){
- list(
- skeleton(
- object,
- x=array(
- data=y,
- dim=c(length(y),1,1),
- dimnames=list(names(y),NULL,NULL)
- ),
- t=t,
- params=as.matrix(parms)
- ),
- NULL
- )
- },
- parms=params[,j],
- ...
- ),
- silent=FALSE
- )
- if (inherits(X,'try-error'))
- stop("trajectory error: error in ",sQuote("lsoda"),call.=FALSE)
- if (attr(X,'istate')[[1]]!=2)
- warning("abnormal exit from ",sQuote("lsoda"),", istate = ",attr(X,'istate'),call.=FALSE)
- x[,j,] <- t(X[,-1])
- }
- },
- unspecified=stop("deterministic skeleton not specified")
- )
- x
- }
- )
+ tm <- times[1]
+ x0 <- init.state(object,params=params,t0=tm)
+ nm <- rownames(x0)
+ dim(x0) <- c(nrow(x0),nrep,1)
+ dimnames(x0) <- list(nm,NULL,NULL)
+
+ x <- array(
+ dim=c(nrow(x0),nrep,length(times)),
+ dimnames=list(rownames(x0),NULL,NULL)
+ )
+ switch(
+ object at skeleton.type,
+ map={ # iterate the map
+ for (k in seq_along(times)) {
+ if (tm < times[k]) {
+ while (tm < times[k]) {
+ x0[,,] <- skeleton(object,x=x0,t=tm,params=params)
+ tm <- tm+1
+ }
+ }
+ x[,,k] <- x0
+ }
+ },
+ vectorfield={ # integrate the vectorfield
+ for (j in seq_len(nrep)) {
+ X <- try(
+ deSolve::lsoda(
+ y=x0[,j,1],
+ times=times,
+ func=function(t,y,parms){
+ list(
+ skeleton(
+ object,
+ x=array(
+ data=y,
+ dim=c(length(y),1,1),
+ dimnames=list(names(y),NULL,NULL)
+ ),
+ t=t,
+ params=as.matrix(parms)
+ ),
+ NULL
+ )
+ },
+ parms=params[,j],
+ ...
+ ),
+ silent=FALSE
+ )
+ if (inherits(X,'try-error'))
+ stop("trajectory error: error in ",sQuote("lsoda"),call.=FALSE)
+ if (attr(X,'istate')[[1]]!=2)
+ warning("abnormal exit from ",sQuote("lsoda"),", istate = ",attr(X,'istate'),call.=FALSE)
+ x[,j,] <- t(X[,-1])
+ }
+ },
+ unspecified=stop("deterministic skeleton not specified")
+ )
+ x
+}
+
+setMethod("trajectory",signature=signature(object="pomp"),definition=trajectory.internal)
More information about the pomp-commits
mailing list