add eeulermultinom

kingaa · Dec 7, 2024 · e25844e · e25844e
1 parent 5af2a3d
commit e25844e
Show file tree

Hide file tree

Showing 103 changed files with 335 additions and 116 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: pomp
 Type: Package
 Title: Statistical Inference for Partially Observed Markov Processes
-Version: 5.11.0.3
-Date: 2024-11-08
+Version: 5.11.1.0
+Date: 2024-12-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")),

diff --git a/TODO.md b/TODO.md
@@ -4,6 +4,7 @@
 
 ## For pomp:
 
+- Rosenzweig-MacArthur example
 - replace `melt` with `data.table::melt`?
 - YAML interface
 - "Getting Started" vignette:

diff --git a/inst/NEWS b/inst/NEWS
@@ -2,6 +2,10 @@ _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._1_1._1:
 
+        • A new addition to the C API allows one to easily compute the
+          expectation of the Euler-multinomial.  The function is
+          ‘eeulermultinom’.
+
         • The help-page discussion of accumulator variables has been
           expanded and clarified.
 

diff --git a/inst/NEWS.Rd b/inst/NEWS.Rd
@@ -2,6 +2,8 @@
 \title{News for package `pomp'}
 \section{Changes in \pkg{pomp} version 5.11.1}{
   \itemize{
+    \item A new addition to the C API allows one to easily compute the expectation of the Euler-multinomial.
+    The function is \code{eeulermultinom}.
     \item The help-page discussion of accumulator variables has been expanded and clarified.
   }
 }

diff --git a/inst/include/pomp.h b/inst/include/pomp.h
@@ -125,6 +125,34 @@ double rgammawn
   return (sigmasq > 0) ? rgamma(dt/sigmasq,sigmasq) : dt;
 }
 
+static R_INLINE
+void eeulermultinom
+(int m, double size, const double *rate, double dt, double *trans)
+{
+  double lambda = 0.0;
+  int j, k;
+  if ( !R_FINITE(size) || size < 0.0 || !R_FINITE(dt) || dt < 0.0) {
+    for (k = 0; k < m; k++) trans[k] = R_NaReal;
+    warn("in 'eeulermultinom': NAs produced.");
+    return;
+  }
+  for (k = 0; k < m; k++) {
+    if (!R_FINITE(rate[k]) || rate[k] < 0.0) {
+      for (j = 0; j < m; j++) trans[j] = R_NaReal;
+      warn("in 'eeulermultinom': NAs produced.");
+      return;
+    }
+    lambda += rate[k];
+  }
+  if (lambda > 0.0) {
+    size = size*(1-exp(-lambda*dt));
+    for (k = 0; k < m; k++)
+      trans[k] = size*rate[k]/lambda;
+  } else {
+    for (k = 0; k < m; k++) trans[k] = 0.0;
+  }
+}
+
 static R_INLINE
 void reulermultinom
 (int m, double size, const double *rate, double dt, double *trans)

diff --git a/src/pomp.h b/src/pomp.h
@@ -133,6 +133,34 @@ double rgammawn
   return (sigmasq > 0) ? rgamma(dt/sigmasq,sigmasq) : dt;
 }
 
+static R_INLINE
+void eeulermultinom
+(int m, double size, const double *rate, double dt, double *trans)
+{
+  double lambda = 0.0;
+  int j, k;
+  if ( !R_FINITE(size) || size < 0.0 || !R_FINITE(dt) || dt < 0.0) {
+    for (k = 0; k < m; k++) trans[k] = R_NaReal;
+    warn("in 'eeulermultinom': NAs produced.");
+    return;
+  }
+  for (k = 0; k < m; k++) {
+    if (!R_FINITE(rate[k]) || rate[k] < 0.0) {
+      for (j = 0; j < m; j++) trans[j] = R_NaReal;
+      warn("in 'eeulermultinom': NAs produced.");
+      return;
+    }
+    lambda += rate[k];
+  }
+  if (lambda > 0.0) {
+    size = size*(1-exp(-lambda*dt));
+    for (k = 0; k < m; k++)
+      trans[k] = size*rate[k]/lambda;
+  } else {
+    for (k = 0; k < m; k++) trans[k] = 0.0;
+  }
+}
+
 static R_INLINE
 void reulermultinom
 (int m, double size, const double *rate, double dt, double *trans)

diff --git a/tests/R_v_C.Rout b/tests/R_v_C.Rout
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 
@@ -109,16 +109,16 @@ Type 'q()' to quit R.
 > ## ----comparison----------------------------------------------------------
 > system.time(simulate(gompertz,nsim=10000,format="arrays"))
    user  system elapsed 
-  5.415   0.114   5.530 
+  5.480   0.064   5.544 
 > system.time(simulate(Gompertz,nsim=10000,format="arrays"))
    user  system elapsed 
-  0.155   0.004   0.160 
+  0.156   0.005   0.161 
 > system.time(pfilter(gompertz,Np=1000))
    user  system elapsed 
-  0.547   0.003   0.549 
+  0.554   0.000   0.554 
 > system.time(pfilter(Gompertz,Np=1000))
    user  system elapsed 
-  0.028   0.000   0.029 
+  0.027   0.001   0.028 
 > 
 > dev.off()
 null device 

diff --git a/tests/abc.Rout.save b/tests/abc.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/baddm.Rout.save b/tests/baddm.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/bake.Rout.save b/tests/bake.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/basic_probes.Rout.save b/tests/basic_probes.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/betabinom.Rout.save b/tests/betabinom.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/blowflies.Rout.save b/tests/blowflies.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/bsmc2.Rout.save b/tests/bsmc2.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/bsplines1.Rout.save b/tests/bsplines1.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/bsplines2.Rout.save b/tests/bsplines2.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/concat.Rout.save b/tests/concat.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/covmat.Rout.save b/tests/covmat.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/csnippet.Rout.save b/tests/csnippet.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/dacca-02.png b/tests/dacca-02.png
diff --git a/tests/dacca.Rout.save b/tests/dacca.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/dinit.Rout.save b/tests/dinit.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/dp-01.png b/tests/dp-01.png
diff --git a/tests/dp-02.png b/tests/dp-02.png
diff --git a/tests/dp-03.png b/tests/dp-03.png
diff --git a/tests/dp-04.png b/tests/dp-04.png
diff --git a/tests/dp-05.png b/tests/dp-05.png
diff --git a/tests/dp-06.png b/tests/dp-06.png
diff --git a/tests/dp.Rout.save b/tests/dp.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/dprocess.Rout.save b/tests/dprocess.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/ebola.Rout.save b/tests/ebola.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/emeasure.Rout.save b/tests/emeasure.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/eulermultinom.Rout.save b/tests/eulermultinom.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/flow.Rout.save b/tests/flow.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/gillespie.Rout.save b/tests/gillespie.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/gompertz.Rout.save b/tests/gompertz.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/helpers.Rout.save b/tests/helpers.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/hitch.Rout.save b/tests/hitch.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/issue109.Rout.save b/tests/issue109.Rout.save
@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

diff --git a/tests/issue222-01.png b/tests/issue222-01.png
diff --git a/tests/issue222.R b/tests/issue222.R
@@ -0,0 +1,67 @@
+library(pomp)
+
+png(filename="issue222-%02d.png",res=100)
+
+set.seed(974650257)
+
+simulate(
+  t0=-1/52,
+  times=seq(0,10,by=1/52),
+  params=c(
+    gamma = 26, mu = 0.02, iota = 0.1,
+    Beta = 400, Beta_sd = 0.01,
+    rho = 0.6, k = 10,
+    pop = 1e6,
+    S_0 = 25/400, I_0 = 0.001, R_0 = 375/400
+  ),
+  rprocess = euler(
+    step.fun=Csnippet(r"{
+      double dW = rgammawn(Beta_sd,dt);
+      double rate[6];
+      double trans[6];
+      rate[0] = mu*pop;
+      rate[1] = Beta*(I+iota)/pop*dW/dt;
+      rate[2] = mu;
+      rate[3] = gamma;
+      rate[4] = mu;
+      rate[5] = mu;
+      trans[0] = rate[0]*dt;
+      eeulermultinom(2,S,&rate[1],dt,&trans[1]);
+      eeulermultinom(2,I,&rate[3],dt,&trans[3]);
+      eeulermultinom(1,R,&rate[5],dt,&trans[5]);
+      // balance equations
+      S += trans[0] - trans[1] - trans[2];
+      I += trans[1] - trans[3] - trans[4];
+      R += trans[3] - trans[5];
+      C += trans[3];
+      W += (dW-dt)/Beta_sd;}"
+    ),
+    delta.t=1/365
+  ),
+  rmeasure = Csnippet(r"{
+      reports = rnbinom_mu(k,rho*C);}"
+  ),
+  rinit=initlz_pf <- Csnippet(r"{
+      double m = pop/(S_0+I_0+R_0);
+      S = m*S_0;
+      I = m*I_0;
+      R = m*R_0;
+      C = 0;
+      W = 0;}"
+  ),
+  partrans=parameter_trans(
+    logit=c("rho"),
+    log=c("gamma","mu","k","Beta","Beta_sd","iota"),
+    barycentric=c("S_0","I_0","R_0")
+  ),
+  obsnames = "reports",
+  accumvars = c("W","C"),
+  statenames = c("S","I","R","C","W"),
+  paramnames = c(
+    "pop","rho","k","gamma","mu","Beta","Beta_sd","iota",
+    "S_0","I_0","R_0"
+  )
+) |>
+  plot()
+
+dev.off()