Running 'abf.R' [0s]
  Running 'abfir.R' [0s]
  Running 'bm.R' [14s]
  Running 'bm2.R' [8s]
  Comparing 'bm2.Rout' to 'bm2.Rout.save' ... OK
  Running 'bpfilter.R' [0s]
  Running 'enkf.R' [0s]
  Running 'gbm.R' [7s]
  Comparing 'gbm.Rout' to 'gbm.Rout.save' ... OK
  Running 'girf.R' [0s]
  Running 'he10.R' [9s]
  Comparing 'he10.Rout' to 'he10.Rout.save' ... OK
  Running 'ibpf.R' [0s]
  Running 'ienkf.R' [0s]
  Running 'igirf.R' [0s]
  Running 'iter_filter.R' [3s]
  Comparing 'iter_filter.Rout' to 'iter_filter.Rout.save' ... OK
  Running 'iubf.R' [0s]
  Running 'lorenz.R' [7s]
  Comparing 'lorenz.Rout' to 'lorenz.Rout.save' ... OK
  Running 'measles.R' [17s]
  Comparing 'measles.Rout' to 'measles.Rout.save' ... OK
  Running 'pfilter.R' [0s]
  Running 'plot.R' [0s]
  Running 'simulate.R' [0s]
  Running 'spatPomp.R' [0s]
  Running 'xtests.R' [3s]
Running the tests in 'tests/bm.R' failed.
Complete output:
  > png(filename="bm-%02d.png",res=100)
  > library(spatPomp)
  Loading required package: pomp
  
  Welcome to pomp!
  
  As of version 4.6, no user-visible pomp function has a name that
  includes a dot ('.'). Function names have been changed to replace the
  dot with an underscore ('_'). For more information, see the pomp blog:
  https://kingaa.github.io/pomp/blog.html.
  
  > 
  > 
  > # For CRAN tests, need to limit to two cores
  > # https://stackoverflow.com/questions/50571325/r-cran-check-fail-when-using-parallel-functions
  > chk <- Sys.getenv("_R_CHECK_LIMIT_CORES_", "")
  > 
  > if (nzchar(chk) && chk == "TRUE") {
  +   # use 2 cores for CRAN
  +   num_workers <- 2L
  + } else {
  +   # use all cores when testing
  +   num_workers <- parallel::detectCores()
  + }
  > num_workers <- 2L
  > 
  > # if(.Platform$OS.type != "windows")
  > #   doParallel::registerDoParallel(num_workers)
  > 
  > # For covr, needs to be single core (https://github.com/r-lib/covr/issues/227)
  > 
  > # CRAN win-builder test fails in foreach for iubf when using a single
  > # core registered with 
  > ###  doParallel::registerDoParallel(1)
  > # so run without registering parallel backend at all
  > # this generates an R warning
  > # Warning message:
  > # executing %dopar% sequentially: no parallel backend registered 
  > # but that is not a major problem
  > 
  > set.seed(2)
  > ## doRNG::registerDoRNG(2)
  > ## using doRNG with 1 core leads to warnings: it seems to make
  > ## foreach confused about whether it is running in parallel or not.
  > 
  > b_model <- bm(U=2,N=2) 
  > 
  > ## ------------------------------------------------------------
  > ## The bm model provides a simple example to test other methods.
  > ## First, we test the filtering methods
  > ## ____________________________________________________________
  > 
  > ##
  > ## exact likelihood via the Kalman filter
  > ##
  > 
  > paste("bm kalman filter loglik: ",round(bm_kalman_logLik(b_model),10))
  [1] "bm kalman filter loglik:  -10.5025878626"
  > 
  > ##
  > ## pfilter tested on bm
  > ##
  > 
  > b_pf <- pfilter(b_model,Np=10)
  > paste("bm pfilter loglik: ",round(logLik(b_pf),10))
  [1] "bm pfilter loglik:  -10.1282113761"
  > 
  > ##
  > ## abf tested on bm. abf uses parallelization, so we also test that
  > ##
  > 
  > b_bag_nbhd <- function(object, time, unit) {
  +   nbhd_list <- list()
  +   if(unit>1) nbhd_list <- c(nbhd_list, list(c(unit-1, time)))
  +   return(nbhd_list)
  + }
  > 
  > set.seed(7)
  > b_abf <- abf(b_model,Nrep=3,Np=10, nbhd = b_bag_nbhd)
  Warning message:
  executing %dopar% sequentially: no parallel backend registered 
  > paste("bm abf loglik: ",round(logLik(b_abf),10))
  [1] "bm abf loglik:  -9.6260083071"
  > 
  > set.seed(7)
  > b_abf_repeat <- abf(b_abf)
  > paste("check abf on abfd_spatPomp: ",
  +   logLik(b_abf_repeat)==logLik(b_abf))
  [1] "check abf on abfd_spatPomp:  TRUE"
  > 
  > ##
  > ## abfir tested on bm
  > ##
  > 
  > b_abfir <- abfir(b_model, Nrep = 3, Np = 10, nbhd = b_bag_nbhd)
  > paste("bm abfir loglik: ",round(logLik(b_abfir),10))
  [1] "bm abfir loglik:  -8.7556872332"
  > 
  > ##
  > ## bpfilter tested on bm
  > ##
  > 
  > set.seed(5)
  > b_bpfilter <- bpfilter(b_model, Np = 10, block_size = 1)
  > paste("bm bpfilter loglik: ",round(logLik(b_bpfilter),10))
  [1] "bm bpfilter loglik:  -11.2089477605"
  > set.seed(5)
  > b_bpfilter_repeat <- bpfilter(b_bpfilter)
  > paste("check bpfilter on bpfilterd_spatPomp: ",
  +   logLik(b_bpfilter)==logLik(b_bpfilter_repeat))
  [1] "check bpfilter on bpfilterd_spatPomp:  TRUE"
  > 
  > ##
  > ## enkf tested on bm
  > ##
  > 
  > b_enkf <- enkf(b_model, Np = 10)
  > paste("bm enkf loglik: ",round(logLik(b_enkf),10))
  [1] "bm enkf loglik:  -11.038124567"
  > 
  > ##
  > ## girf tested on bm, both moment and bootstrap methods
  > ##
  > 
  > b_girf_mom <- girf(b_model,Np = 5,lookahead = 1,Nguide = 5,
  +   kind = 'moment',Ninter=2)
  > paste("bm girf loglik, moment guide: ",round(logLik(b_girf_mom),10))
  [1] "bm girf loglik, moment guide:  -13.2689915662"
  > 
  > set.seed(0)
  > b_girf_boot <- girf(b_model,Np = 5,lookahead = 1,Nguide = 5,
  +   kind = 'bootstrap',Ninter=2)
  > paste("bm girf loglik, bootstrap guide: ",round(logLik(b_girf_boot),10))
  [1] "bm girf loglik, bootstrap guide:  -14.3228315101"
  > 
  > set.seed(0)
  > b_girf_boot_repeat <- girf(b_girf_boot)
  > paste("check girf on girfd_spatPomp: ",
  +   logLik(b_girf_boot)==logLik(b_girf_boot_repeat))
  [1] "check girf on girfd_spatPomp:  TRUE"
  > 
  > ## ------------------------------------------------------------
  > ## Now, we test the inference methods
  > ## ____________________________________________________________
  > 
  > b_rw.sd <- rw_sd(rho=0.02,X1_0=ivp(0.02))
  > 
  > ##
  > ## igirf on bm
  > ##
  > ## A call to igirf using the moment-based guide function can test compiled
  > ## code for eunit_measure, munit_measure, vunit_measure, dunit_measure,
  > ## runit_measure, rprocess, skeleton, rinit and partrans. 
  > ##
  > ## we test both geometric and hyperbolic cooling
  > 
  > set.seed(1)
  > b_igirf_geom <- igirf(b_model,
  +   Ngirf = 2,
  +   rw.sd = b_rw.sd,
  +   cooling.type = "geometric",
  +   cooling.fraction.50 = 0.5,
  +   Np=10,
  +   Ninter = 2,
  +   lookahead = 1,
  +   Nguide = 5,
  +   kind = 'moment',
  +   verbose = FALSE
  + )
  > paste("bm igirf loglik, geometric cooling, verbose=F: ",round(logLik(b_igirf_geom),10))
  [1] "bm igirf loglik, geometric cooling, verbose=F:  -9.9887892835"
  > 
  > set.seed(1)
  > b_igirf_geom_repeat <- igirf(b_igirf_geom,params=coef(b_model))
  > paste("check igirf on igirfd_spatPomp: ",
  +   logLik(b_igirf_geom)==logLik(b_igirf_geom_repeat))
  [1] "check igirf on igirfd_spatPomp:  FALSE"
  > 
  > 
  > b_igirf_hyp <- igirf(b_model,
  +   Ngirf = 2,
  +   rw.sd = b_rw.sd,
  +   cooling.type = "hyperbolic",
  +   cooling.fraction.50 = 0.5,
  +   Np=10,
  +   Ninter = 2,
  +   lookahead = 1,
  +   Nguide = 5,
  +   kind = 'moment',
  +   verbose = TRUE
  + )
  igirf iteration 1 of 2 completed with likelihood  -8.864864 
           rho        sigma          tau         X1_0         X2_0 
  0.4042940956 1.0000000000 1.0000000000 0.0004004701 0.0000000000 
  igirf iteration 2 of 2 completed with likelihood  -10.7886 
          rho       sigma         tau        X1_0        X2_0 
   0.40021397  1.00000000  1.00000000 -0.03993169  0.00000000 
  > paste("bm igirf loglik, hyperbolic cooling, verbose=T: ",round(logLik(b_igirf_hyp),10))
  [1] "bm igirf loglik, hyperbolic cooling, verbose=T:  -10.7886016055"
  > 
  > plot(b_igirf_geom) -> b_igirf_plot
  > head(b_igirf_plot$data)
    iteration  param       value
  2         2 loglik -10.8072342
  3         3 loglik  -9.9887893
  4         1    rho   0.4000000
  5         2    rho   0.4014437
  6         3    rho   0.3984085
  7         1  sigma   1.0000000
  > 
  > ##
  > ## ienkf on bm, with geometric and hyperbolic cooling
  > ##
  > 
  > b_ienkf_geom <- ienkf(b_model,
  +   Nenkf=2,
  +   Np = 10,
  +   rw.sd=b_rw.sd,
  +   cooling.type="geometric",
  +   cooling.fraction.50 = 0.5,
  +   verbose=FALSE
  + )
  > paste("bm ienkf loglik, geometric cooling, verbose=F: ",round(logLik(b_ienkf_geom),10))
  [1] "bm ienkf loglik, geometric cooling, verbose=F:  -10.7262815137"
  > 
  > b_ienkf_hyp <- ienkf(b_model,
  +   Nenkf=2,
  +   Np = 10,
  +   rw.sd=b_rw.sd,
  +   cooling.type="hyperbolic",
  +   cooling.fraction.50 = 0.5,
  +   verbose=TRUE
  + )
  ienkf iteration 1 of 2 completed
  ienkf iteration 2 of 2 completed
  > paste("bm ienkf loglik, hyperbolic cooling, verbose=T: ",round(logLik(b_ienkf_hyp),10))
  [1] "bm ienkf loglik, hyperbolic cooling, verbose=T:  -10.4670683154"
  > 
  > ##
  > ## iubf on bm, with geometric and hyperbolic cooling
  > ##
  > 
  > b_iubf_geom <- iubf(b_model,
  +   Nubf = 2,
  +   Nrep_per_param = 5,
  +   Nparam = 5,
  +   nbhd = b_bag_nbhd,
  +   prop = 0.8,
  +   rw.sd =b_rw.sd,
  +   cooling.type = "geometric",
  +   cooling.fraction.50 = 0.5,
  +   verbose=FALSE
  + )
  > paste("bm iubf loglik, geometric cooling, verbose=F: ",round(logLik(b_iubf_geom),10))
  [1] "bm iubf loglik, geometric cooling, verbose=F:  -9.1342209153"
  > 
  > b_iubf_hyp <- iubf(b_model,
  +   Nubf = 2,
  +   Nrep_per_param = 5,
  +   Nparam = 5,
  +   nbhd = b_bag_nbhd,
  +   prop = 0.8,
  +   rw.sd =b_rw.sd,
  + #  cooling.type = "hyperbolic",
  +   cooling.type = "geometric",
  +   cooling.fraction.50 = 0.5,
  +   verbose=TRUE
  + )
  working on observation times  1  in iteration  1 
  working on observation times  2  in iteration  1 
  iubf iteration 1 of 2 completed with log-likelihood -10.75135 
         rho      sigma        tau       X1_0       X2_0 
  0.39859037 1.00000000 1.00000000 0.02810935 0.00000000 
  working on observation times  1  in iteration  2 
  working on observation times  2  in iteration  2 
  iubf iteration 2 of 2 completed with log-likelihood -12.00495 
         rho      sigma        tau       X1_0       X2_0 
  0.39755252 1.00000000 1.00000000 0.01749735 0.00000000 
  > paste("bm ienkf loglik, hyperbolic cooling, verbose=T: ",round(logLik(b_iubf_hyp),10))
  [1] "bm ienkf loglik, hyperbolic cooling, verbose=T:  -12.0049482564"
  > 
  > ## --------------------------------------------
  > ## using bm to test simulate and plot
  > ## ____________________________________________
  > 
  > b_sim1 <- simulate(b_model,nsim=2,format='data.frame')
  > head(b_sim1,10)
  # A tibble: 8 x 5
     time .id   unitname       X      Y
    <dbl> <ord> <chr>      <dbl>  <dbl>
  1     1 1     U1        0.358   1.13 
  2     1 2     U1       -0.0475 -0.546
  3     1 1     U2        1.17    0.815
  4     1 2     U2        0.700   1.98 
  5     2 1     U1        1.17   -0.211
  6     2 2     U1       -0.411  -0.363
  7     2 1     U2        1.60    0.596
  8     2 2     U2        2.51    2.31 
  > b_sim2 <- simulate(b_model,nsim=2,format='data.frame',include.data=TRUE)
  > head(b_sim2,10)
  # A tibble: 10 x 5
      time .id   unitname      X      Y
     <dbl> <ord> <chr>     <dbl>  <dbl>
   1     1 data  U1        2.98   3.61 
   2     1 1     U1        0.425  0.282
   3     1 2     U1        1.15   0.411
   4     1 data  U2        0.824  2.81 
   5     1 1     U2        0.286 -1.66 
   6     1 2     U2        0.224 -0.285
   7     2 data  U1        2.03   1.73 
   8     2 1     U1        3.13   2.95 
   9     2 2     U1        0.518 -0.312
  10     2 data  U2       -0.759 -0.849
  > b_sim3 <- simulate(b_model,nsim=2,format='spatPomps')
  > 
  > plot(b_model,type="l",log=FALSE)
  > b_sim3v2 <- b_sim3[[1]]
  > b_sim3v2@data <- exp(b_sim3v2@data)
  > plot(b_sim3v2,type="l",log=TRUE)
  > plot(b_sim3[[2]],type="h")
  > 
  > ## --------------------------------------------
  > ## using bm to test spatPomp workhorse functions, extending pomp:
  > ## vunit_measure, eunit_measure, munit_measure, dunit_measure
  > ##
  > ## these are tested implicitly in the methods, but here is
  > ## a more direct test
  > ## ____________________________________________
  > 
  > 
  > b_s <- states(b_model)[,1,drop=FALSE]
  > dim(b_s) <- c(dim(b_s),1)
  > dimnames(b_s) <- list(variable=dimnames(states(b_model))[[1]], rep=NULL)
  > b_p <- coef(b_model)
  > dim(b_p) <- c(length(b_p),1)
  > dimnames(b_p) <- list(param=names(coef(b_model)))
  > b_y <- obs(b_model)[,1,drop=FALSE]
  > 
  > vunit_measure(b_model, x=b_s, unit=2, time=1, params=b_p)
  , , 1
  
        rep
  unit   [,1]
    [1,]    1
  
  > 
  > eunit_measure(b_model, x=b_s, unit=2, time=1, params=b_p)
  , , 1
  
        rep
  unit        [,1]
    [1,] 0.8235357
  
  > 
  > b_array.params <- array(b_p,
  +   dim = c(length(b_p),length(unit_names(b_model)), 1, 1),
  +   dimnames = list(params = rownames(b_p)))
  > b_vc <- c(4, 9) # this should have length equal to the number of units
  > dim(b_vc) <- c(length(b_vc), 1, 1)
  > 
  > munit_measure(b_model, x=b_s, vc=b_vc, Np=1, unit = 1, time=1,
  +   params=b_array.params)
  , , 1, 1
  
         
  params  [,1]
    rho    0.4
    sigma  1.0
    tau    2.0
    X1_0   0.0
    X2_0   0.0
  
  > 
  > dunit_measure(b_model, y=b_y,
  +   x=b_s, unit=1, time=1, params=b_p)
        time
  rep         [,1]
    [1,] -1.112688
  > 
  > runit_measure(b_model, x=b_s, unit=2, time=1, params=b_p)
  , , 1
  
          rep
  variable      [,1]
        Y1 0.9852548
  
  > 
  > vec_rmeasure(b_model,x=b_s,time=1, params=b_p)
  , , 1
  
            [,1]
  [1,] 1.7469874
  [2,] 0.7308696
  
  > 
  > ## --------------------------------------------
  > ## using bm to test edge cases and utility functions
  > ## perhaps only of technical interest
  > ## ____________________________________________
  > 
  > print(b_model)
  <object of class 'spatPomp'>
  > 
  > # check how u is treated by dunit_measure, runit_measure, eunit_measure,
  > # vunit_measure and munit_measure. this should output unit-1 to
  > # be consistent with Csnippet indexing.
  > 
  > b_u <- spatPomp(b_model,
  +   dunit_measure=spatPomp_Csnippet("lik=u;"),
  +   eunit_measure=spatPomp_Csnippet("ey=u;"),
  +   munit_measure=spatPomp_Csnippet("M_tau=u;"),
  +   vunit_measure=spatPomp_Csnippet("vc=u;"),
  +   runit_measure=spatPomp_Csnippet("Y=u;")  
  + )
  > 
  > vunit_measure(b_u, x=b_s, unit=2, time=1, params=b_p)
  , , 1
  
        rep
  unit   [,1]
    [1,]    1
  
  > eunit_measure(b_u, x=b_s, unit=2, time=1, params=b_p)
  , , 1
  
        rep
  unit   [,1]
    [1,]    1
  
  > munit_measure(b_u, x=b_s, vc=b_vc, Np=1, unit = 2, time=1,params=b_array.params)
  , , 1, 1
  
         
  params  [,1]
    rho    0.4
    sigma  1.0
    tau    1.0
    X1_0   0.0
    X2_0   0.0
  
  > dunit_measure(b_u, y=b_y,x=b_s, unit=2, time=1, params=b_p)
        time
  rep    [,1]
    [1,]    1
  > runit_measure(b_u, x=b_s, unit=2, time=1, params=b_p)
  , , 1
  
          rep
  variable [,1]
        Y1    1
  
  > 
  > dev.off()
  null device 
            1 
  > 
  > ## test spatPomp_Csnippet variable construction
  > spatPomp_Csnippet("lik=u;",unit_statenames="A",unit_obsnames=c("B","C"), unit_covarnames="D",
  +   unit_ivpnames="E",unit_paramnames="F",unit_vfnames="G")
  An object of class "Csnippet"
  Slot "text":
  [1] "double *A = &A1;\nconst double *B = &B1;\nconst double *C = &C1;\nconst double *D = &D1;\nconst double *E_0 = &E1_0;\nconst double *F = &F1;\ndouble *DG = &DG1;\nlik=u;"
  
  > 
  > ## --------------------------------------------
  > ## using bm to test spatPomp() replacement functionality
  > ## ____________________________________________
  > 
  > b_rep1 <- spatPomp(b_model,params=coef(b_model))
  > for(slt in slotNames(b_model)) if(!identical(slot(b_model,slt),slot(b_rep1,slt))) print(slt)
  [1] "states"
  > 
  > try(spatPomp(data=as.data.frame(b_model),units=NULL),outFile=stdout())
  Error in value[[3L]](cond) : 
    'times' should be a single name identifying the column of data that represents the observation times. 'units' should be likewise for column that represents the observation units.
  > 
  >   
  > 
  > 
  > 
  > proc.time()
     user  system elapsed 
     7.83    0.49   13.93