onestep takes one step

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: pomp
 Type: Package
 Title: Statistical Inference for Partially Observed Markov Processes
-Version: 5.9.0.0
-Date: 2024-06-02
+Version: 5.9.1.0
+Date: 2024-06-07
 Authors@R: c(person(given=c("Aaron","A."),family="King",role=c("aut","cre"),email="[email protected]",comment=c(ORCID="0000-0001-6159-3207")),
 	  person(given=c("Edward","L."),family="Ionides",role="aut",comment=c(ORCID="0000-0002-4190-0174"))	,
 	  person(given="Carles",family="Bretó",role="aut",comment=c(ORCID="0000-0003-4695-4902")),

R/rprocess_spec.R

@@ -52,7 +52,7 @@
 ##' They differ as to how this is done.
 ##' Specifically,
 ##' \enumerate{
-##' \item \code{onestep} takes a single step to go from any given time \code{t1} to any later time \code{t2} (\code{t1 < t2}).
+##' \item \code{onestep} takes a single step to go from any given time \code{t1} to any later time \code{t2} (\code{t1 <= t2}).
 ##' Thus, this plug-in is designed for use in situations where a closed-form solution to the process exists.
 ##' \item To go from \code{t1} to \code{t2}, \code{euler} takes \code{n} steps of equal size, where
 ##' \preformatted{

inst/NEWS

@@ -1,5 +1,16 @@
 _N_e_w_s _f_o_r _p_a_c_k_a_g_e '_p_o_m_p'
 
+_C_h_a_n_g_e_s _i_n '_p_o_m_p' _v_e_r_s_i_o_n _5._9._1:
+
+        • The rprocess plugin ‘onestep’ now always takes exactly one
+          step to get from one observation time to the next, even when
+          the interval between them is zero.  That is, if ‘P’ is a
+          ‘pomp’ object with an rprocess built using
+          ‘rprocess=onestep(f)’, ‘t0=timezero(P)’, and ‘t=time(P)’,
+          then ‘f’ will be called exactly once for the interval
+          (‘t0’,‘t[1]’) and once for each of the intervals
+          (‘t[i-1]’,‘t[i]’), ‘i=2,...,length(t)’.
+
 _C_h_a_n_g_e_s _i_n '_p_o_m_p' _v_e_r_s_i_o_n _5._8._4:
 
         • To add additional elements for use by the basic model

inst/NEWS.Rd

@@ -1,5 +1,11 @@
 \name{NEWS}
 \title{News for package `pomp'}
+\section{Changes in \pkg{pomp} version 5.9.1}{
+  \itemize{
+    \item The rprocess plugin \code{onestep} now always takes exactly one step to get from one observation time to the next, even when the interval between them is zero.
+    That is, if \code{P} is a \sQuote{pomp} object with an rprocess built using \code{rprocess=onestep(f)}, \code{t0=timezero(P)}, and \code{t=time(P)}, then \code{f} will be called exactly once for the interval (\code{t0},\code{t[1]}) and once for each of the intervals (\code{t[i-1]},\code{t[i]}), \code{i=2,...,length(t)}.
+  }
+}
 \section{Changes in \pkg{pomp} version 5.8.4}{
   \itemize{
     \item To add additional elements for use by the basic model components (i.e., \dQuote{userdata}), the \code{userdata} argument should now be used.

src/euler.c

@@ -208,7 +208,7 @@ SEXP euler_simulator
     switch (method) {
     case onestep: default:      // one step
       dt = time[step]-t;
-      nstep = (dt > 0) ? 1 : 0;
+      nstep = 1;
       break;
     case discrete:              // fixed step
       dt = deltat;

src/rprocess.c

@@ -111,7 +111,7 @@ SEXP do_rprocess (SEXP object, SEXP xstart, SEXP tstart, SEXP times, SEXP params
       double deltat = 1.0;
       PROTECT(fn = GET_SLOT(rproc,install("step.fn")));
       PROTECT(X = euler_simulator(fn,xstart,tstart,times,params,deltat,type,
-                                        accumvars,covar,args,gnsi));
+                                  accumvars,covar,args,gnsi));
       nprotect += 2;
     }
     break;
@@ -122,7 +122,7 @@ SEXP do_rprocess (SEXP object, SEXP xstart, SEXP tstart, SEXP times, SEXP params
       PROTECT(fn = GET_SLOT(rproc,install("step.fn")));
       deltat = *(REAL(AS_NUMERIC(GET_SLOT(rproc,install("delta.t")))));
       PROTECT(X = euler_simulator(fn,xstart,tstart,times,params,deltat,type,
-                                        accumvars,covar,args,gnsi));
+                                  accumvars,covar,args,gnsi));
       nprotect += 2;
     }
     break;

tests/steps.R

@@ -1,9 +1,10 @@
 options(digits=3)
 
-library(pomp)
 suppressPackageStartupMessages({
+  library(tidyr)
   library(dplyr)
 })
+library(pomp)
 
 tm <- freeze(sort(runif(n=20,max=50)),seed=1930502785L)
 simulate(
@@ -98,3 +99,45 @@ stopifnot(
     as.numeric(states(sm,"steps"))
   )
 )
+
+simulate(
+  t0=0,
+  times=c(0,0,0,1,3,3,7),
+  rinit=\(t0,...)c(x=t0,n=0),
+  rprocess=onestep(\(t,x,n,delta.t,...) c(x=t+delta.t,n=n+1)),
+  rmeasure=\(x,n,...) c(y1=x,y2=n),
+  dmeasure=\(...,log) if (log) 0 else 1
+) -> s1
+
+simulate(
+  t0=0,
+  times=c(0,0,0,1,3,3,7),
+  rinit=Csnippet("x=t; n=0;"),
+  rprocess=onestep(Csnippet("x=t+dt; n=n+1;")),
+  rmeasure=Csnippet("y1 = x; y2 = n;"),
+  dmeasure=Csnippet("lik = (give_log) ? 0 : 1;"),
+  statenames=c("x","n"),
+  obsnames=c("y1","y2")
+) -> s2
+
+s1 |>
+  as.data.frame() |>
+  with(
+    stopifnot(
+      time==x,
+      n==seq_along(time)
+    )
+  )
+
+s2 |>
+  pfilter(Np=1,save.states=TRUE) |>
+  saved_states(format="list") |>
+  melt() |>
+  pivot_wider() |>
+  mutate(time=time(s2)[.L1]) |>
+  with(
+    stopifnot(
+      time==x,
+      n==seq_along(time)
+    )
+  )

tests/steps.Rout.save

@@ -19,10 +19,11 @@ Type 'q()' to quit R.
 
 > options(digits=3)
 > 
-> library(pomp)
 > suppressPackageStartupMessages({
++   library(tidyr)
 +   library(dplyr)
 + })
+> library(pomp)
 > 
 > tm <- freeze(sort(runif(n=20,max=50)),seed=1930502785L)
 > simulate(
@@ -120,3 +121,45 @@ Warning message:
 +   )
 + )
 > 
+> simulate(
++   t0=0,
++   times=c(0,0,0,1,3,3,7),
++   rinit=\(t0,...)c(x=t0,n=0),
++   rprocess=onestep(\(t,x,n,delta.t,...) c(x=t+delta.t,n=n+1)),
++   rmeasure=\(x,n,...) c(y1=x,y2=n),
++   dmeasure=\(...,log) if (log) 0 else 1
++ ) -> s1
+> 
+> simulate(
++   t0=0,
++   times=c(0,0,0,1,3,3,7),
++   rinit=Csnippet("x=t; n=0;"),
++   rprocess=onestep(Csnippet("x=t+dt; n=n+1;")),
++   rmeasure=Csnippet("y1 = x; y2 = n;"),
++   dmeasure=Csnippet("lik = (give_log) ? 0 : 1;"),
++   statenames=c("x","n"),
++   obsnames=c("y1","y2")
++ ) -> s2
+> 
+> s1 |>
++   as.data.frame() |>
++   with(
++     stopifnot(
++       time==x,
++       n==seq_along(time)
++     )
++   )
+> 
+> s2 |>
++   pfilter(Np=1,save.states=TRUE) |>
++   saved_states(format="list") |>
++   melt() |>
++   pivot_wider() |>
++   mutate(time=time(s2)[.L1]) |>
++   with(
++     stopifnot(
++       time==x,
++       n==seq_along(time)
++     )
++   )
+>