fix resample bug (thanks to Kidus and Ed)

kingaa · Jun 14, 2018 · 15ad754 · 15ad754
1 parent 47629be
commit 15ad754
Show file tree

Hide file tree

Showing 16 changed files with 59 additions and 75 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: pomp
 Type: Package
 Title: Statistical Inference for Partially Observed Markov Processes
-Version: 1.16.2.3
-Date: 2018-05-16
+Version: 1.16.2.4
+Date: 2018-06-14
 Authors@R: c(person(given=c("Aaron","A."),family="King",
 		role=c("aut","cre"),email="[email protected]"),
 	  person(given=c("Edward","L."),family="Ionides",role=c("aut")),

diff --git a/inst/NEWS b/inst/NEWS
@@ -2,6 +2,8 @@ _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 _1._1_6._2:
 
+        • Fix bug in systematic resampler when ‘np != nw’.
+
         • Fix bug in ‘conv.rec’ when only one parameter is requested.
 
         • Correct and improve documentation for ‘conv.rec’.

diff --git a/inst/NEWS.Rd b/inst/NEWS.Rd
@@ -2,6 +2,7 @@
 \title{News for package `pomp'}
 \section{Changes in \pkg{pomp} version 1.16.2}{
   \itemize{
+    \item Fix bug in systematic resampler when \code{np != nw}.
     \item Fix bug in \code{conv.rec} when only one parameter is requested.
     \item Correct and improve documentation for \code{conv.rec}.
     \item \code{coef<-} can now take a list of parameters: it simply applies \code{unlist} to turn the list into a numeric vector.

diff --git a/src/resample.c b/src/resample.c
@@ -20,7 +20,7 @@ void nosort_resamp (int nw, double *w, int np, int *p, int offset)
     u += du;
     // In the following line, the second test is needed to correct 
     // the infamous Bug of St. Patrick, 2017-03-17.
-    while ((u > w[i]) && (i < np-1)) i++;
+    while ((u > w[i]) && (i < nw-1)) i++;
     p[j] = i;
   }
   if (offset)			// add offset if needed

diff --git a/tests/bake.Rout.save b/tests/bake.Rout.save
@@ -29,7 +29,7 @@ attr(,"seed")
 [1] 499586
 attr(,"system.time")
    user  system elapsed 
-  0.000   0.004   0.000 
+      0       0       0 
 > w2 <- runif(2)
 > freeze(runif(5),seed=499586)
 [1] 0.6719240 0.5092283 0.0643100 0.3672672 0.7148970

diff --git a/tests/blowflies-02.png b/tests/blowflies-02.png
diff --git a/tests/demos.Rout b/tests/demos.Rout
@@ -67,7 +67,7 @@ Type 'q()' to quit R.
 > toc <- Sys.time()
 
 > print(toc - tic)
-Time difference of 0.004065514 secs
+Time difference of 0.006081343 secs
 
 > X <- reshape(X, dir = "wide", v.names = "X", timevar = "traj", 
 +     idvar = "time")
@@ -82,7 +82,7 @@ Time difference of 0.004065514 secs
 > toc <- Sys.time()
 
 > print(toc - tic)
-Time difference of 0.003564119 secs
+Time difference of 0.004382133 secs
 
 > x <- reshape(x, dir = "wide", v.names = c("Y", "X"), 
 +     timevar = "sim", idvar = "time")
@@ -105,7 +105,7 @@ Time difference of 0.003564119 secs
 > toc <- Sys.time()
 
 > print(toc - tic)
-Time difference of 0.5819807 secs
+Time difference of 0.6539347 secs
 
 > print(logmeanexp(sapply(pf, logLik), se = TRUE))
                    se 
@@ -194,7 +194,7 @@ Time difference of 0.5819807 secs
 > toc <- Sys.time()
 
 > print(toc - tic)
-Time difference of 0.003418207 secs
+Time difference of 0.004687071 secs
 
 > plot(incid ~ time, data = X, type = "l")
 
@@ -205,7 +205,7 @@ Time difference of 0.003418207 secs
 > toc <- Sys.time()
 
 > print(toc - tic)
-Time difference of 0.03128076 secs
+Time difference of 0.03365993 secs
 
 > plot(incid ~ time, data = x, col = as.factor(x$sim), 
 +     pch = 16)
@@ -221,7 +221,7 @@ Time difference of 0.03128076 secs
 > toc <- Sys.time()
 
 > print(toc - tic)
-Time difference of 7.518277 secs
+Time difference of 7.390307 secs
 Warning messages:
 1: in 'table_lookup': extrapolating at 1.934011e+03 
 2: in 'table_lookup': extrapolating at 1.934011e+03 
@@ -231,5 +231,5 @@ null device
           1 
 > 
 > print(packageVersion("pomp"))
-[1] ‘1.16.2.1’
+[1] ‘1.16.2.4’
 > 
diff --git a/tests/gompertz.Rout.save b/tests/gompertz.Rout.save
@@ -50,12 +50,12 @@ process model simulator, rprocess = function (xstart, times, params, ..., zerona
             stop(ep, conditionMessage(e), call. = FALSE)
         })
 }
-<bytecode: 0x39702d8>
-<environment: 0x3332bc8>
+<bytecode: 0x23c14d0>
+<environment: 0x1d84270>
 process model density, dprocess = function (x, times, params, log = FALSE, ...) 
 stop(sQuote("dprocess"), " not specified", call. = FALSE)
-<bytecode: 0x4354938>
-<environment: 0x3228d20>
+<bytecode: 0x2da9cd0>
+<environment: 0x1c79e18>
 measurement model simulator, rmeasure = native function ‘_gompertz_normal_rmeasure’
 measurement model density, dmeasure = native function ‘_gompertz_normal_dmeasure’
 prior simulator, rprior = not specified
@@ -66,12 +66,12 @@ parameter transformation (to estimation scale) = function (params, ...)
 {
     log(params)
 }
-<environment: 0x32354a0>
+<environment: 0x1c86598>
 parameter transformation (from estimation scale) = function (params, ...) 
 {
     exp(params)
 }
-<environment: 0x32354a0>
+<environment: 0x1c86598>
 parameter(s):
     K     r sigma   tau   X.0 
   1.0   0.1   0.1   0.1   1.0 
@@ -91,8 +91,8 @@ function (xstart, times, params, ..., zeronames = character(0),
             stop(ep, conditionMessage(e), call. = FALSE)
         })
 }
-<bytecode: 0x39702d8>
-<environment: 0x3332bc8>
+<bytecode: 0x23c14d0>
+<environment: 0x1d84270>
 > 
 > po <- gompertz
 > coef(po)

diff --git a/tests/mif2-07.png b/tests/mif2-07.png
diff --git a/tests/mif2-08.png b/tests/mif2-08.png
diff --git a/tests/mif2.Rout.save b/tests/mif2.Rout.save
@@ -138,11 +138,10 @@ truth   -482
 +                alpha.3=0.2*(time<10)))
 > 
 > m4 <- mif2(m4,Nmif=10,Np=c(rep(200,20),rep(100,80),200))
-Warning message:
-in ‘mif2.pfilter’: 3 filtering failures occurred. 
+Warning messages:
+1: in ‘mif2.pfilter’: 17 filtering failures occurred. 
+2: in ‘mif2.pfilter’: 1 filtering failure occurred. 
 > m4 <- continue(m4,Nmif=10,cooling.fraction=0.1)
-Warning message:
-in ‘mif2.pfilter’: 25 filtering failures occurred. 
 > try(continue(m4,Np=function(k)if(k<10) "B" else 500))
 Error : in ‘mif2’: if ‘Np’ is a function, it must return a single positive integer
 > try(continue(m4,Np=function(k)if(k<10) -30 else 500))
@@ -160,31 +159,13 @@ Error : in ‘mif2’: number of particles, ‘Np’, must always be positive
 +                    alpha.2=ifelse(time==1,0.2,0.1),
 +                    alpha.3=0.2*(time<10)))
 + ) -> out
-filtering failure at time t = 33
-filtering failure at time t = 34
-filtering failure at time t = 35
-filtering failure at time t = 36
-filtering failure at time t = 37
-filtering failure at time t = 38
-filtering failure at time t = 39
-filtering failure at time t = 40
-filtering failure at time t = 41
-filtering failure at time t = 42
-filtering failure at time t = 43
-filtering failure at time t = 44
-filtering failure at time t = 45
-filtering failure at time t = 46
-filtering failure at time t = 47
-filtering failure at time t = 48
-filtering failure at time t = 49
-filtering failure at time t = 50
-filtering failure at time t = 51
-filtering failure at time t = 74
-filtering failure at time t = 75
 filtering failure at time t = 76
 filtering failure at time t = 77
+filtering failure at time t = 78
+filtering failure at time t = 79
+filtering failure at time t = 80
 Warning message:
-in ‘mif2.pfilter’: 23 filtering failures occurred. 
+in ‘mif2.pfilter’: 5 filtering failures occurred. 
 > stopifnot(length(out)==210)
 > stopifnot(sum(grepl("mif2 iteration",out))==10)
 > stopifnot(sum(grepl("mif2 pfilter",out))==200)
@@ -194,22 +175,22 @@ in ‘mif2.pfilter’: 23 filtering failures occurred.
 > try(mif2(m5,Np=100,Nmif=2))
 Error : in ‘mif2’: in ‘mif2.pfilter’: ‘dmeasure’ returns non-finite value.
 likelihood, data, states, and parameters are:
-   time:            3
+   time:            2
     lik:          NaN
-     y1:      3.76759
-     y2:      2.97521
+     y1:      0.24217
+     y2:      7.03942
      x1:          Inf
      x2:          NaN
 alpha.1:          0.8
 alpha.2:         -Inf
-alpha.3:     0.229848
+alpha.3:      0.85781
 alpha.4:          0.9
 sigma.1:            3
 sigma.2:         -0.5
 sigma.3:            2
     tau:            1
-   x1.0:     -8.67937
-   x2.0:      2.73245
+   x1.0:     0.121014
+   x2.0:      1.74011
 > 
 > library(ggplot2)
 > library(reshape2)
@@ -241,8 +222,8 @@ likelihood, data, states, and parameters are:
   lik:          NaN
     Y:      1.11862
     X:          NaN
-    K:     -0.95075
-    r:      0.14582
+    K:     -1.01089
+    r:      0.17549
 sigma:          0.1
   tau:          0.1
   X.0:            1
@@ -264,7 +245,7 @@ likelihood, data, states, and parameters are:
   lik:          NaN
     Y:      1.11862
     X:          NaN
-    K:      1.70288
+    K:      6.43134
     r:          NaN
 sigma:          0.1
   tau:          0.1

diff --git a/tests/pfilter.Rout.save b/tests/pfilter.Rout.save
@@ -44,7 +44,7 @@ Error : in ‘pfilter’: if ‘Np’ is a function, it must return a single pos
 > ll.se <- sd(exp(pf.ll-mean(pf.ll)))/exp(ll.est-mean(pf.ll))/sqrt(length(pf))
 > print(round(c(loglik=ll.est,loglik.se=ll.se),digits=2))
    loglik loglik.se 
-  -480.47      0.90 
+  -479.79      0.41 
 > 
 > pompExample(euler.sir)
 In ‘pomp’: the following unrecognized argument(s) will be stored for use by user-defined functions: ‘nbasis’,‘degree’,‘period’

diff --git a/tests/ricker.Rout.save b/tests/ricker.Rout.save
@@ -134,12 +134,12 @@ process model simulator, rprocess = function (xstart, times, params, ..., zerona
             stop(ep, conditionMessage(e), call. = FALSE)
         })
 }
-<bytecode: 0x31dd858>
-<environment: 0x2f4c528>
+<bytecode: 0x38de4d0>
+<environment: 0x364c4f0>
 process model density, dprocess = function (x, times, params, log = FALSE, ...) 
 stop(sQuote("dprocess"), " not specified", call. = FALSE)
-<bytecode: 0x463a728>
-<environment: 0x46407d0>
+<bytecode: 0x4d3a368>
+<environment: 0x4d40410>
 measurement model simulator, rmeasure = native function ‘_ricker_poisson_rmeasure’, dynamically loaded from ‘pomp’
 measurement model density, dmeasure = native function ‘_ricker_poisson_dmeasure’, dynamically loaded from ‘pomp’
 prior simulator, rprior = not specified
@@ -152,14 +152,14 @@ parameter transformation (to estimation scale) = function (params, ...)
         "sigma", "phi", "c", "N.0")])
     params
 }
-<environment: 0x4407a88>
+<environment: 0x4b076c8>
 parameter transformation (from estimation scale) = function (params, ...) 
 {
     params[c("r", "sigma", "phi", "c", "N.0")] <- exp(params[c("r", 
         "sigma", "phi", "c", "N.0")])
     params
 }
-<environment: 0x4407a88>
+<environment: 0x4b076c8>
 parameter(s):
        r    sigma      phi        c      N.0      e.0 
 44.70118  0.30000 10.00000  1.00000  7.00000  0.00000 

diff --git a/tests/rw2.Rout.save b/tests/rw2.Rout.save
@@ -182,8 +182,8 @@ measurement model simulator, rmeasure = function (x, t, params, covars, ...)
     }
     y
 }
-<bytecode: 0x480f4b8>
-<environment: 0x2498338>
+<bytecode: 0x3a730f8>
+<environment: 0x16fc2c8>
 measurement model density, dmeasure = function (y, x, t, params, log, covars, ...) 
 {
     f <- 0
@@ -195,8 +195,8 @@ measurement model density, dmeasure = function (y, x, t, params, log, covars, ..
         f
     else exp(f)
 }
-<bytecode: 0x480a338>
-<environment: 0x2498338>
+<bytecode: 0x3a6df78>
+<environment: 0x16fc2c8>
 prior simulator, rprior = not specified
 prior density, dprior = native function ‘_pomp_default_dprior’, dynamically loaded from ‘pomp’
 skeleton = not specified

diff --git a/tests/sir.Rout.save b/tests/sir.Rout.save
@@ -193,8 +193,8 @@ process model simulator, rprocess = function (xstart, times, params, ..., zerona
             stop(ep, conditionMessage(e), call. = FALSE)
         })
 }
-<bytecode: 0x2ee39a0>
-<environment: 0x255aa98>
+<bytecode: 0x36d7678>
+<environment: 0x2d59a28>
 process model density, dprocess = function (x, times, params, ..., tcovar, covar, log = FALSE, 
     .getnativesymbolinfo = TRUE) 
 {
@@ -204,8 +204,8 @@ process model density, dprocess = function (x, times, params, ..., tcovar, covar
         stop(ep, conditionMessage(e), call. = FALSE)
     })
 }
-<bytecode: 0x2e2ba08>
-<environment: 0x23ea800>
+<bytecode: 0x36264d8>
+<environment: 0x2be9758>
 measurement model simulator, rmeasure = native function ‘__pomp_rmeasure’, defined by a Csnippet
 measurement model density, dmeasure = function (y, x, t, params, log, covars, ...) 
 {
@@ -236,7 +236,7 @@ skeleton (vectorfield) = function (x, t, params, covars, ...)
         xdot
     })
 }
-<environment: 0x4446bd0>
+<environment: 0x4c45810>
 initializer = function (params, t0, ...) 
 {
     with(as.list(params), {

diff --git a/tests/verbosity.Rout b/tests/verbosity.Rout
@@ -101,16 +101,16 @@ filtering failure at time t = 3
 filtering failure at time t = 20
 filtering failure at time t = 20
 make[2]: Entering directory `/userdata/kingaa/projects/Rpkg/pomp/tests'
-gcc -std=gnu99 -I"/usr/local/apps/R/R-3.5.0/lib64/R/include" -DNDEBUG   -I/usr/local/include  -I/userdata/kingaa/projects/Rpkg/library/pomp/include -fpic  -g -O2 -Wall -pedantic -c /tmp/RtmpMe0KBl/57083/pomp_75e4ce3619203951f1d91ddcea2d584a.c -o /tmp/RtmpMe0KBl/57083/pomp_75e4ce3619203951f1d91ddcea2d584a.o
-/tmp/RtmpMe0KBl/57083/pomp_75e4ce3619203951f1d91ddcea2d584a.c: In function ‘__pomp_rmeasure’:
-/tmp/RtmpMe0KBl/57083/pomp_75e4ce3619203951f1d91ddcea2d584a.c:16:1: error: expected ‘;’ before ‘}’ token
+gcc -std=gnu99 -I"/usr/local/apps/R/R-3.5.0/lib64/R/include" -DNDEBUG -I/userdata/kingaa/projects/Rpkg/library/pomp/include  -I/usr/local/include   -fpic  -g -O2 -Wall -pedantic -c /tmp/Rtmp9WaSWM/105624/pomp_1645b0f13d33cd7dfcde6d7eee46f917.c -o /tmp/Rtmp9WaSWM/105624/pomp_1645b0f13d33cd7dfcde6d7eee46f917.o
+/tmp/Rtmp9WaSWM/105624/pomp_1645b0f13d33cd7dfcde6d7eee46f917.c: In function ‘__pomp_rmeasure’:
+/tmp/Rtmp9WaSWM/105624/pomp_1645b0f13d33cd7dfcde6d7eee46f917.c:16:1: error: expected ‘;’ before ‘}’ token
  }
  ^
-make[2]: *** [/tmp/RtmpMe0KBl/57083/pomp_75e4ce3619203951f1d91ddcea2d584a.o] Error 1
+make[2]: *** [/tmp/Rtmp9WaSWM/105624/pomp_1645b0f13d33cd7dfcde6d7eee46f917.o] Error 1
 make[2]: Leaving directory `/userdata/kingaa/projects/Rpkg/pomp/tests'
 make[2]: Entering directory `/userdata/kingaa/projects/Rpkg/pomp/tests'
 make[2]: Leaving directory `/userdata/kingaa/projects/Rpkg/pomp/tests'
-Error : in ‘pomp’: error in building shared-object library from C snippets: in ‘pompCBuilder’: compilation error: cannot compile shared-object library ‘/tmp/RtmpMe0KBl/57083/pomp_75e4ce3619203951f1d91ddcea2d584a.so’: status = 1
+Error : in ‘pomp’: error in building shared-object library from C snippets: in ‘pompCBuilder’: compilation error: cannot compile shared-object library ‘/tmp/Rtmp9WaSWM/105624/pomp_1645b0f13d33cd7dfcde6d7eee46f917.so’: status = 1
 In addition: Warning messages:
 1: in ‘bsmc2’: 2 filtering failures occurred. 
 2: in ‘pfilter’: 1 filtering failure occurred.