[Pomp-commits] r302 - pkg/R

[noreply at r-forge.r-project.org]

Author: kingaa
Date: 2010-08-24 20:06:08 +0200 (Tue, 24 Aug 2010)
New Revision: 302

Modified:
   pkg/R/trajectory-pomp.R
Log:
- modify the algorithm in the case of a discrete-time map so that the assumption that 'times' are consecutive is no longer implicitly made


Modified: pkg/R/trajectory-pomp.R
===================================================================
--- pkg/R/trajectory-pomp.R	2010-08-24 18:04:54 UTC (rev 301)
+++ pkg/R/trajectory-pomp.R	2010-08-24 18:06:08 UTC (rev 302)
@@ -27,14 +27,19 @@
             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)
 
-            params <- as.matrix(params)
-            if (missing(times))
-              times <- time(object,t0=TRUE)
-            x0 <- init.state(object,params=params,t0=times[1])
+            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)
@@ -42,21 +47,21 @@
             switch(
                    object at skeleton.type,
                    map={                # iterate the map
-                     x[,,1] <- x0
-                     for (k in seq(from=2,to=length(times),by=1)) {
-                       x[,,k] <- skeleton(
-                                          object,
-                                          x=x[,,k-1,drop=FALSE],
-                                          t=times[k-1],
-                                          params=params
-                                          )
+                     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],
+                                               y=x0[,j,1],
                                                times=times,
                                                func=function(t,y,parms){
                                                  list(