R Under development (unstable) (2024-02-17 r85935 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > png(filename="bm-%02d.png",res=100) > library(spatPomp) Loading required package: pomp > i <- 1 > U <- switch(i,2,10) > N <- switch(i,2,10) > Np <- switch(i,10,100) > > # 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=U,N=N) > > ## ------------------------------------------------------------ > ## 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=Np) > 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=2,Np=Np, 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.6241814039" > > 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 = 2, Np = Np, nbhd = b_bag_nbhd) > paste("bm abfir loglik: ",round(logLik(b_abfir),10)) [1] "bm abfir loglik: -9.6437796804" > > ## > ## bpfilter tested on bm > ## > > set.seed(5) > b_bpfilter <- bpfilter(b_model, Np = Np, 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" > > set.seed(5) > b_bpfilter_filter_traj <- bpfilter(b_bpfilter,filter_traj=TRUE) > paste("bpfilter filter trajectory final particle: ") [1] "bpfilter filter trajectory final particle: " > round(b_bpfilter_filter_traj@filter.traj[,1,],3) time name [,1] [,2] [,3] X1 0 1.373 2.349 X2 0 2.035 0.796 > > > set.seed(5) > b_bpfilter_save_states <- bpfilter(b_bpfilter,save_states=TRUE) > paste("bpfilter final particles: ") [1] "bpfilter final particles: " > round(b_bpfilter_save_states@saved.states[[N]],3) .id name [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] X1 2.611 2.236 2.186 2.186 2.349 1.077 1.077 1.155 1.442 1.063 X2 0.796 0.442 0.442 0.231 0.231 0.231 0.648 0.648 1.337 1.019 > > ## > ## enkf tested on bm > ## > > set.seed(5) > b_enkf <- enkf(b_model, Np = Np) > paste("bm enkf loglik: ",round(logLik(b_enkf),10)) [1] "bm enkf loglik: -10.8897941231" > > ## > ## girf on bm: moment and bootstrap methods, followed by error tests > ## > > set.seed(0) > b_girf_mom <- girf(b_model,Np = floor(Np/2),lookahead = 1, + Nguide = floor(Np/2), + kind = 'moment',Ninter=2) > paste("bm girf loglik, moment guide: ",round(logLik(b_girf_mom),10)) [1] "bm girf loglik, moment guide: -14.9828908185" > > ## for boostrap girf, we do not set Ninter, to test the default which is Ninter=U > set.seed(0) > b_girf_boot <- girf(b_model,Np = floor(Np/2),lookahead = 1, + Nguide = floor(Np/2), + kind = 'bootstrap') > 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" > > print("The following delivers an error message, to test it") [1] "The following delivers an error message, to test it" > try(girf()) Error in girf() : girf: data is a required argument. > try(girf(object="nonsense")) Error in girf(object = "nonsense") : girf is undefined for 'nonsense'of class 'character'. > try(girf(b_girf_boot,Np=c(Inf))) Error : in 'girf': 'Np' must be a single positive integer > try(girf(b_girf_boot,Np=seq(from=10,length=N+1,by=2))) Error : in 'girf': 'Np' must be a single positive integer > > ## ------------------------------------------------------------ > ## 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=Np, + 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=floor(Np/2), + Ninter = 2, + lookahead = 1, + Nguide = floor(Np/2), + kind = 'moment', + verbose = TRUE + ) igirf iteration 1 of 2 completed with likelihood -11.56221 rho sigma tau X1_0 X2_0 0.39957108 1.00000000 1.00000000 0.01379269 0.00000000 igirf iteration 2 of 2 completed with likelihood -10.59957 rho sigma tau X1_0 X2_0 0.40074887 1.00000000 1.00000000 0.04083778 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.599572921" > > 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 = Np, + 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: -11.2659371179" > > b_ienkf_hyp <- ienkf(b_model, + Nenkf=2, + Np = Np, + 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.4969604389" > > ## > ## iubf on bm, with geometric and hyperbolic cooling > ## > > b_iubf_geom <- iubf(b_model, + Nubf = 2, + Nrep_per_param = floor(Np/2), + Nparam = floor(Np/2), + 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: -11.4595005548" > > b_iubf_hyp <- iubf(b_model, + Nubf = 2, + Nrep_per_param = floor(Np/2), + Nparam = floor(Np/2), + 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 -12.04241 rho sigma tau X1_0 X2_0 0.40312565 1.00000000 1.00000000 -0.01326395 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.83948 rho sigma tau X1_0 X2_0 0.40523197 1.00000000 1.00000000 -0.00474667 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.8394821627" > > ## -------------------------------------------- > ## using bm to test simulate and plot > ## ____________________________________________ > > b_sim1 <- simulate(b_model,nsim=2,format='data.frame') > head(b_sim1,10) # A tibble: 8 × 5 time .id unitname X Y 1 1 1 U1 -1.12 -0.918 2 1 2 U1 -0.670 -0.945 3 1 1 U2 -1.07 -0.409 4 1 2 U2 -0.648 -1.24 5 2 1 U1 -2.81 -4.87 6 2 2 U1 0.850 0.589 7 2 1 U2 -1.26 -1.96 8 2 2 U2 -0.491 -0.0759 > b_sim2 <- simulate(b_model,nsim=2,format='data.frame',include.data=TRUE) > head(b_sim2,10) # A tibble: 10 × 5 time .id unitname X Y 1 1 data U1 2.98 3.61 2 1 1 U1 -0.827 -1.76 3 1 2 U1 0.417 3.17 4 1 data U2 0.824 2.81 5 1 1 U2 -0.228 -0.282 6 1 2 U2 -0.158 -1.84 7 2 data U1 2.03 1.73 8 2 1 U1 0.164 -0.282 9 2 2 U1 -0.270 0.730 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 1.181054 > > vec_rmeasure(b_model,x=b_s,time=1, params=b_p) , , 1 [,1] [1,] 0.8543340 [2,] -0.5666138 > > ## -------------------------------------------- > ## using bm to test edge cases and utility functions > ## perhaps only of technical interest > ## ____________________________________________ > > print(b_model) > > # 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" > > # test an error message > print("The following delivers an error message, to test it") [1] "The following delivers an error message, to test it" > 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. > > # test parameter replacement > b_rep2 <- spatPomp(b_model,params=coef(b_model)+1) > if(!identical(coef(b_rep2),coef(b_model)+1)) stop('problem with parameter replacement') > > ## ----------------------------------------------------------------- > ## using bm to test behavior of inference methods when logLik = -Inf > ## _________________________________________________________________ > > > b_inf <- spatPomp(b_model, + dmeasure = spatPomp_Csnippet( + method="dmeasure", + unit_statenames="X", + unit_obsnames="Y", + code = " + lik=0; + if(give_log) lik = log(lik); + " + ), + dunit_measure = spatPomp_Csnippet(" + lik = 0; + if(give_log) lik = log(lik); + ") + ) > > dmeasure(b_inf,x=states(b_model)) time .id [,1] [,2] [1,] 0 0 > > dmeasure(b_inf,x=states(b_model),log=T) time .id [,1] [,2] [1,] -Inf -Inf > > vec_dmeasure(b_inf,y=obs(b_inf),x=states(b_model),units=1:U,times=1:2,params=coef(b_inf),log=T)[,,1] [,1] [,2] [1,] -Inf -Inf [2,] -Inf -Inf > > b_girf_mom_inf <- girf(b_inf,Np = floor(Np/2),lookahead = 1, + Nguide = floor(Np/2), + kind = 'moment',Ninter=2) > paste("bm moment girf loglik, zero measurement: ",round(logLik(b_girf_mom_inf),10)) [1] "bm moment girf loglik, zero measurement: -Inf" > > set.seed(0) > b_girf_boot_inf <- girf(b_inf,Np = floor(Np/2),lookahead = 1, + Nguide = floor(Np/2), + kind = 'bootstrap',Ninter=2) > paste("bm bootstrap girf loglik, zero measurement: ",round(logLik(b_girf_boot_inf),10)) [1] "bm bootstrap girf loglik, zero measurement: -Inf" > > set.seed(5) > b_bpfilter_inf <- bpfilter(b_inf, Np = Np, block_size = 1) > paste("bm bpfilter loglik, zero measurement: ",round(logLik(b_bpfilter_inf),10)) [1] "bm bpfilter loglik, zero measurement: NaN" > > > > proc.time() user system elapsed 3.70 0.32 12.26