better exception trapping in eulermultinom

kingaa · kingaa · commit 387daf3d79ab · 2023-06-10T07:33:49.000-04:00
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: pomp
 Type: Package
 Title: Statistical Inference for Partially Observed Markov Processes
-Version: 5.2.0.0
-Date: 2023-06-06
+Version: 5.2.1.0
+Date: 2023-06-09
 Authors@R: c(person(given=c("Aaron","A."),family="King",role=c("aut","cre"),email="kingaa@umich.edu",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/inst/NEWS b/inst/NEWS
@@ -1,5 +1,11 @@
 _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._2._1:
+
+        • ‘reulermultinom’ now returns ‘NA’ rather than ‘NaN’, in
+          keeping with the behavior of ‘rbinom’.  Thanks to John Drake
+          for calling attention to this issue.
+
 _C_h_a_n_g_e_s _i_n '_p_o_m_p' _v_e_r_s_i_o_n _5._1._1:
 
         • A bug in ‘stew’ resulting in namespace conflicts has been
diff --git a/inst/NEWS.Rd b/inst/NEWS.Rd
@@ -1,5 +1,11 @@
 \name{NEWS}
 \title{News for package `pomp'}
+\section{Changes in \pkg{pomp} version 5.2.1}{
+  \itemize{
+    \item \code{reulermultinom} now returns \code{NA} rather than \code{NaN}, in keeping with the behavior of \code{rbinom}.
+    Thanks to John Drake for calling attention to this issue.
+  }
+}
 \section{Changes in \pkg{pomp} version 5.1.1}{
   \itemize{
     \item A bug in \code{stew} resulting in namespace conflicts has been repaired.
diff --git a/inst/include/pomp.h b/inst/include/pomp.h
@@ -44,27 +44,26 @@ static R_INLINE void reulermultinom (int m, double size, const double *rate,
                                      double dt, double *trans) {
   double p = 0.0;
   int j, k;
-  if ((size < 0.0) || (dt < 0.0) || (floor(size+0.5) != size)) {
-    for (k = 0; k < m; k++) trans[k] = R_NaN;
+  if ( !R_FINITE(size) || size < 0.0 || floor(size+0.5) != size ||
+       !R_FINITE(dt) || dt < 0.0) {
+    for (k = 0; k < m; k++) trans[k] = R_NaReal;
+    warningcall(R_NilValue,"in 'reulermultinom': NAs produced.");
     return;
   }
   for (k = 0; k < m; k++) {
-    if (rate[k] < 0.0) {
-      for (j = 0; j < m; j++) trans[j] = R_NaN;
+    if (!R_FINITE(rate[k]) || rate[k] < 0.0) {
+      for (j = 0; j < m; j++) trans[j] = R_NaReal;
+      warningcall(R_NilValue,"in 'reulermultinom': NAs produced.");
       return;
     }
     p += rate[k]; // total event rate
   }
   if (p > 0.0) {
     size = rbinom(size,1-exp(-p*dt)); // total number of events
-    if (!(R_FINITE(size)))
-      warningcall(R_NilValue,"in 'reulermultinom': result of binomial draw is not finite.");
     m -= 1;
     for (k = 0; k < m; k++) {
       if (rate[k] > p) p = rate[k];
       trans[k] = ((size > 0) && (p > 0)) ? rbinom(size,rate[k]/p) : 0;
-      if (!(R_FINITE(size)&&R_FINITE(p)&&R_FINITE(rate[k])&&R_FINITE(trans[k])))
-        warningcall(R_NilValue,"in 'reulermultinom': result of binomial draw is not finite.");
       size -= trans[k];
       p -= rate[k];
     }
diff --git a/src/distributions.c b/src/distributions.c
@@ -4,42 +4,45 @@
 
 #include "pomp_internal.h"
 
-static void reulermultinom_multi (int *n, int *ntrans, double *size, double *rate, double *deltat, double *trans) {
+static void reulermultinom_multi (int m, int n, double *size, double *rate, double *deltat, double *trans) {
   int k;
-  int m = *ntrans;
-  for (k = 0; k < *n; k++) {
-    reulermultinom(*ntrans,*size,rate,*deltat,trans);
+  for (k = 0; k < n; k++) {
+    reulermultinom(m,*size,rate,*deltat,trans);
     trans += m;
   }
 }
 
-static void deulermultinom_multi (int *n, int *ntrans, double *size, double *rate, double *deltat, double *trans, int *give_log, double *f) {
+static void deulermultinom_multi (int m, int n, double *size, double *rate, double *deltat, double *trans, int *give_log, double *f) {
   int k;
-  int m = *ntrans;
-  for (k = 0; k < *n; k++) {
-    f[k] = deulermultinom(*ntrans,*size,rate,*deltat,trans,*give_log);
+  for (k = 0; k < n; k++) {
+    f[k] = deulermultinom(m,*size,rate,*deltat,trans,*give_log);
     trans += m;
   }
 }
 
 SEXP R_Euler_Multinom (SEXP n, SEXP size, SEXP rate, SEXP deltat) {
   int ntrans = length(rate);
   int dim[2];
-  SEXP nn, x, nm;
+  SEXP x, nm;
   if (length(size)>1)
     warn("in 'reulermultinom': only the first element of 'size' is meaningful");
   if (length(deltat)>1)
     warn("in 'reulermultinom': only the first element of 'dt' is meaningful");
-  PROTECT(nn = AS_INTEGER(n));
+  PROTECT(n = AS_INTEGER(n));
+  PROTECT(size = AS_NUMERIC(size));
+  PROTECT(rate = AS_NUMERIC(rate));
+  PROTECT(deltat = AS_NUMERIC(deltat));
   dim[0] = ntrans;
-  dim[1] = INTEGER(nn)[0];
+  dim[1] = *INTEGER(n);
+  if (dim[1] == NA_INTEGER || dim[1] < 0)
+    err("in 'reulermultinom': 'n' must be a non-negative integer.");
   PROTECT(x = makearray(2,dim));
   PROTECT(nm = GET_NAMES(rate));
   setrownames(x,nm,2);
   GetRNGstate();
-  reulermultinom_multi(INTEGER(nn),&ntrans,REAL(size),REAL(rate),REAL(deltat),REAL(x));
+  reulermultinom_multi(dim[0],dim[1],REAL(size),REAL(rate),REAL(deltat),REAL(x));
   PutRNGstate();
-  UNPROTECT(3);
+  UNPROTECT(6);
   return x;
 }
 
@@ -56,8 +59,12 @@ SEXP D_Euler_Multinom (SEXP x, SEXP size, SEXP rate, SEXP deltat, SEXP log) {
   if (length(deltat)>1)
     warn("in 'deulermultinom': only the first element of 'dt' is meaningful");
   PROTECT(f = NEW_NUMERIC(n));
-  deulermultinom_multi(&n,&ntrans,REAL(size),REAL(rate),REAL(deltat),REAL(x),INTEGER(log),REAL(f));
-  UNPROTECT(1);
+  PROTECT(size = AS_NUMERIC(size));
+  PROTECT(rate = AS_NUMERIC(rate));
+  PROTECT(deltat = AS_NUMERIC(deltat));
+  PROTECT(log = AS_LOGICAL(log));
+  deulermultinom_multi(ntrans,n,REAL(size),REAL(rate),REAL(deltat),REAL(x),INTEGER(log),REAL(f));
+  UNPROTECT(5);
   return f;
 }
 
diff --git a/src/pomp.h b/src/pomp.h
@@ -44,27 +44,26 @@ static R_INLINE void reulermultinom (int m, double size, const double *rate,
                                      double dt, double *trans) {
   double p = 0.0;
   int j, k;
-  if ((size < 0.0) || (dt < 0.0) || (floor(size+0.5) != size)) {
-    for (k = 0; k < m; k++) trans[k] = R_NaN;
+  if ( !R_FINITE(size) || size < 0.0 || floor(size+0.5) != size ||
+       !R_FINITE(dt) || dt < 0.0) {
+    for (k = 0; k < m; k++) trans[k] = R_NaReal;
+    warningcall(R_NilValue,"in 'reulermultinom': NAs produced.");
     return;
   }
   for (k = 0; k < m; k++) {
-    if (rate[k] < 0.0) {
-      for (j = 0; j < m; j++) trans[j] = R_NaN;
+    if (!R_FINITE(rate[k]) || rate[k] < 0.0) {
+      for (j = 0; j < m; j++) trans[j] = R_NaReal;
+      warningcall(R_NilValue,"in 'reulermultinom': NAs produced.");
       return;
     }
     p += rate[k]; // total event rate
   }
   if (p > 0.0) {
     size = rbinom(size,1-exp(-p*dt)); // total number of events
-    if (!(R_FINITE(size)))
-      warningcall(R_NilValue,"in 'reulermultinom': result of binomial draw is not finite.");
     m -= 1;
     for (k = 0; k < m; k++) {
       if (rate[k] > p) p = rate[k];
       trans[k] = ((size > 0) && (p > 0)) ? rbinom(size,rate[k]/p) : 0;
-      if (!(R_FINITE(size)&&R_FINITE(p)&&R_FINITE(rate[k])&&R_FINITE(trans[k])))
-        warningcall(R_NilValue,"in 'reulermultinom': result of binomial draw is not finite.");
       size -= trans[k];
       p -= rate[k];
     }
diff --git a/tests/R_v_C.Rout b/tests/R_v_C.Rout
@@ -117,16 +117,16 @@ https://kingaa.github.io/pomp/blog.html.
 > ## ----comparison----------------------------------------------------------
 > system.time(simulate(gompertz,nsim=10000,format="arrays"))
    user  system elapsed 
-  7.079   0.038   7.118 
+  7.012   0.069   7.082 
 > system.time(simulate(Gompertz,nsim=10000,format="arrays"))
    user  system elapsed 
-  0.148   0.011   0.160 
+  0.156   0.004   0.160 
 > system.time(pfilter(gompertz,Np=1000))
    user  system elapsed 
-  0.653   0.031   0.683 
+  0.649   0.000   0.649 
 > system.time(pfilter(Gompertz,Np=1000))
    user  system elapsed 
-  0.028   0.000   0.028 
+  0.029   0.000   0.028 
 > 
 > dev.off()
 null device 
diff --git a/tests/eulermultinom.R b/tests/eulermultinom.R
@@ -15,6 +15,12 @@ reulermultinom(n=1,size=c(100,200,300),rate=c(1,2,3),dt=0.2)
 reulermultinom(n=1,size=0,rate=c(1,2,3),dt=0.2)
 reulermultinom(n=1,size=10,rate=c(1,Inf,1),0.1)
 reulermultinom(n=1,size=Inf,rate=c(1,100,1),0.1)
+try(reulermultinom(n=NA,size=100,rate=c(1,2,3),dt=1))
+reulermultinom(n=5,size=NA,rate=c(1,2,3),dt=1)
+reulermultinom(n=5,size=100,rate=c(1,NA,3),dt=1)
+reulermultinom(n=5,size=100,rate=c(1,2,3),dt=NA)
+reulermultinom(n=5,size=100,rate=c(0,0,0,0),dt=1)
+
 x <- reulermultinom(n=3,size=100,rate=c(3,2,1),dt=0.1)
 try(deulermultinom(rbind(x,c(0,1,0)),size=100,rate=c(3,2,1),dt=0.1))
 deulermultinom(x,size=c(100,NA),rate=c(3,2,1),dt=0.1)
@@ -24,6 +30,7 @@ deulermultinom(x=c(3,4,0),size=10,rate=c(1,1,1),dt=0-.1)
 deulermultinom(x=c(3,4,0),size=10,rate=c(1,1,-1),dt=0.1)
 deulermultinom(x=c(3,-4,0),size=10,rate=c(1,1,1),dt=0.1)
 deulermultinom(x=c(3,6,3),size=10,rate=c(1,1,0),dt=0.1,log=TRUE)
+
 rgammawn(n=5,sigma=2,dt=0.1)
 rgammawn(n=5,sigma=1:5,dt=0.1)
 rgammawn(n=5,sigma=1,dt=rep(1,5))
diff --git a/tests/eulermultinom.Rout.save b/tests/eulermultinom.Rout.save
@@ -36,39 +36,62 @@ https://kingaa.github.io/pomp/blog.html.
 [3,]   20   22   21   23   20
 > reulermultinom(n=5,size=-3,rate=c(1,2,3),dt=0.1)
      [,1] [,2] [,3] [,4] [,5]
-[1,]  NaN  NaN  NaN  NaN  NaN
-[2,]  NaN  NaN  NaN  NaN  NaN
-[3,]  NaN  NaN  NaN  NaN  NaN
+[1,]   NA   NA   NA   NA   NA
+[2,]   NA   NA   NA   NA   NA
+[3,]   NA   NA   NA   NA   NA
+Warning messages:
+1: in 'reulermultinom': NAs produced. 
+2: in 'reulermultinom': NAs produced. 
+3: in 'reulermultinom': NAs produced. 
+4: in 'reulermultinom': NAs produced. 
+5: in 'reulermultinom': NAs produced. 
 > reulermultinom(n=5,size=100,rate=c(1,-2,3),dt=0.1)
      [,1] [,2] [,3] [,4] [,5]
-[1,]  NaN  NaN  NaN  NaN  NaN
-[2,]  NaN  NaN  NaN  NaN  NaN
-[3,]  NaN  NaN  NaN  NaN  NaN
+[1,]   NA   NA   NA   NA   NA
+[2,]   NA   NA   NA   NA   NA
+[3,]   NA   NA   NA   NA   NA
+Warning messages:
+1: in 'reulermultinom': NAs produced. 
+2: in 'reulermultinom': NAs produced. 
+3: in 'reulermultinom': NAs produced. 
+4: in 'reulermultinom': NAs produced. 
+5: in 'reulermultinom': NAs produced. 
 > reulermultinom(n=5,size=100,rate=c(1,NA,3),dt=0.1)
      [,1] [,2] [,3] [,4] [,5]
-[1,]    0    0    0    0    0
-[2,]    0    0    0    0    0
-[3,]    0    0    0    0    0
+[1,]   NA   NA   NA   NA   NA
+[2,]   NA   NA   NA   NA   NA
+[3,]   NA   NA   NA   NA   NA
+Warning messages:
+1: in 'reulermultinom': NAs produced. 
+2: in 'reulermultinom': NAs produced. 
+3: in 'reulermultinom': NAs produced. 
+4: in 'reulermultinom': NAs produced. 
+5: in 'reulermultinom': NAs produced. 
 > reulermultinom(n=5,size=100.3,rate=c(1,2,3),dt=0.1)
      [,1] [,2] [,3] [,4] [,5]
-[1,]  NaN  NaN  NaN  NaN  NaN
-[2,]  NaN  NaN  NaN  NaN  NaN
-[3,]  NaN  NaN  NaN  NaN  NaN
+[1,]   NA   NA   NA   NA   NA
+[2,]   NA   NA   NA   NA   NA
+[3,]   NA   NA   NA   NA   NA
+Warning messages:
+1: in 'reulermultinom': NAs produced. 
+2: in 'reulermultinom': NAs produced. 
+3: in 'reulermultinom': NAs produced. 
+4: in 'reulermultinom': NAs produced. 
+5: in 'reulermultinom': NAs produced. 
 > reulermultinom(n=0,size=100,rate=c(1,2,3),dt=0.1)
     
 [1,]
 [2,]
 [3,]
 > try(reulermultinom(n=-2,size=100,rate=c(1,2,3),dt=0.1))
-Error : in 'reulermultinom': negative length vectors are not allowed
+Error : in 'reulermultinom': in 'reulermultinom': 'n' must be a non-negative integer.
 > reulermultinom(n=1,size=100,rate=c(1,2e400,3),dt=0.1)
      [,1]
-[1,]    0
-[2,]  NaN
-[3,]  NaN
-Warning messages:
-1: in 'reulermultinom': result of binomial draw is not finite. 
-2: in 'reulermultinom': result of binomial draw is not finite. 
+[1,]   NA
+[2,]   NA
+[3,]   NA
+Warning message:
+in 'reulermultinom': NAs produced. 
 > reulermultinom(n=1,size=100,rate=c(1,2,3),dt=c(0.1,0.2,0.3,Inf))
      [,1]
 [1,]   11
@@ -90,21 +113,60 @@ in 'reulermultinom': only the first element of 'size' is meaningful
 [3,]    0
 > reulermultinom(n=1,size=10,rate=c(1,Inf,1),0.1)
      [,1]
-[1,]    0
-[2,]  NaN
-[3,]  NaN
-Warning messages:
-1: in 'reulermultinom': result of binomial draw is not finite. 
-2: in 'reulermultinom': result of binomial draw is not finite. 
+[1,]   NA
+[2,]   NA
+[3,]   NA
+Warning message:
+in 'reulermultinom': NAs produced. 
 > reulermultinom(n=1,size=Inf,rate=c(1,100,1),0.1)
      [,1]
-[1,]    0
-[2,]    0
-[3,]  NaN
+[1,]   NA
+[2,]   NA
+[3,]   NA
+Warning message:
+in 'reulermultinom': NAs produced. 
+> try(reulermultinom(n=NA,size=100,rate=c(1,2,3),dt=1))
+Error : in 'reulermultinom': in 'reulermultinom': 'n' must be a non-negative integer.
+> reulermultinom(n=5,size=NA,rate=c(1,2,3),dt=1)
+     [,1] [,2] [,3] [,4] [,5]
+[1,]   NA   NA   NA   NA   NA
+[2,]   NA   NA   NA   NA   NA
+[3,]   NA   NA   NA   NA   NA
+Warning messages:
+1: in 'reulermultinom': NAs produced. 
+2: in 'reulermultinom': NAs produced. 
+3: in 'reulermultinom': NAs produced. 
+4: in 'reulermultinom': NAs produced. 
+5: in 'reulermultinom': NAs produced. 
+> reulermultinom(n=5,size=100,rate=c(1,NA,3),dt=1)
+     [,1] [,2] [,3] [,4] [,5]
+[1,]   NA   NA   NA   NA   NA
+[2,]   NA   NA   NA   NA   NA
+[3,]   NA   NA   NA   NA   NA
+Warning messages:
+1: in 'reulermultinom': NAs produced. 
+2: in 'reulermultinom': NAs produced. 
+3: in 'reulermultinom': NAs produced. 
+4: in 'reulermultinom': NAs produced. 
+5: in 'reulermultinom': NAs produced. 
+> reulermultinom(n=5,size=100,rate=c(1,2,3),dt=NA)
+     [,1] [,2] [,3] [,4] [,5]
+[1,]   NA   NA   NA   NA   NA
+[2,]   NA   NA   NA   NA   NA
+[3,]   NA   NA   NA   NA   NA
 Warning messages:
-1: in 'reulermultinom': result of binomial draw is not finite. 
-2: in 'reulermultinom': result of binomial draw is not finite. 
-3: in 'reulermultinom': result of binomial draw is not finite. 
+1: in 'reulermultinom': NAs produced. 
+2: in 'reulermultinom': NAs produced. 
+3: in 'reulermultinom': NAs produced. 
+4: in 'reulermultinom': NAs produced. 
+5: in 'reulermultinom': NAs produced. 
+> reulermultinom(n=5,size=100,rate=c(0,0,0,0),dt=1)
+     [,1] [,2] [,3] [,4] [,5]
+[1,]    0    0    0    0    0
+[2,]    0    0    0    0    0
+[3,]    0    0    0    0    0
+[4,]    0    0    0    0    0
+> 
 > x <- reulermultinom(n=3,size=100,rate=c(3,2,1),dt=0.1)
 > try(deulermultinom(rbind(x,c(0,1,0)),size=100,rate=c(3,2,1),dt=0.1))
 Error : in 'deulermultinom': NROW('x') should match length of 'rate'
@@ -130,6 +192,7 @@ in 'deulermultinom': NaNs produced.
 [1] 0
 > deulermultinom(x=c(3,6,3),size=10,rate=c(1,1,0),dt=0.1,log=TRUE)
 [1] -Inf
+> 
 > rgammawn(n=5,sigma=2,dt=0.1)
 [1] 2.593438e-18 2.417425e-05 5.489565e-25 1.344180e-18 3.775061e-03
 > rgammawn(n=5,sigma=1:5,dt=0.1)
diff --git a/tests/spy.Rout b/tests/spy.Rout

-Original file line number
+Diff line change
@@ @@ -1,5 +1,11 @@ @@
 _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._2._1:
++
 +        • ‘reulermultinom’ now returns ‘NA’ rather than ‘NaN’, in
 +          keeping with the behavior of ‘rbinom’.  Thanks to John Drake
 +          for calling attention to this issue.
++
 _C_h_a_n_g_e_s _i_n '_p_o_m_p' _v_e_r_s_i_o_n _5._1._1:
         • A bug in ‘stew’ resulting in namespace conflicts has been