R Under development (unstable) (2024-08-19 r87028 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.

> library(spatPomp)
Loading required package: pomp
> 
> ## test error message
> try(bm2(U=3,N=2,unit_specific_names="rho",shared_names="sigma"))
Error in bm2(U = 3, N = 2, unit_specific_names = "rho", shared_names = "sigma") : 
  both shared_names and unit_specific names cannot be given to bm2
> try(bm2_kalman_logLik(bm2(U=1,N=2)))
Error in bm2_kalman_logLik(bm2(U = 1, N = 2)) : bm2 designed for U>1
> 
> b2_U <- 4
> 
> set.seed(0)
> b2 <- bm2(U=b2_U,N=2,unit_specific_names="rho")
> 
> set.seed(0)
> b2alt <- bm2(U=b2_U,N=2,shared_names=c("tau","sigma","X_0"))
> if(any(obs(b2)!=obs(b2alt))) stop("two different ways to specify the same bm2 model should match")
> 
> b2_bpfilter <- bpfilter(b2,Np=5,block_size=1)
> paste("bpfilter logLik for bm2 model:",logLik(b2_bpfilter))
[1] "bpfilter logLik for bm2 model: -12.9121635525516"
> 
> # here there are no transformations so use small rw.sd. to avoid negatives
> b2_rw_list <- rep(list(0.001),times=b2_U) 
> names(b2_rw_list) <-paste0("rho",1:b2_U)
> b2_rw.sd <- do.call(rw_sd,b2_rw_list)
> 
> b2_units = seq_len(b2_U)
> b2_nblocks = b2_U/2
> b2_block_list = split(b2_units, sort(b2_units %% b2_nblocks))
> b2_block_list <- lapply(b2_block_list, as.integer)
> 
> b2_ibpf <- ibpf(b2,
+   params=coef(b2),
+   sharedParNames=NULL,
+   unitParNames="rho",
+   Nbpf=2,
+   spat_regression=0.1,
+   Np=5,
+   rw.sd=b2_rw.sd,
+   cooling.fraction.50=0.5,
+   block_list=b2_block_list
+ )
> 
> paste("ibpf logLik for b2 model:",logLik(b2_ibpf))
[1] "ibpf logLik for b2 model: -14.1817035202796"
> 
> paste("kf logLik for b2:",bm2_kalman_logLik(b2))
[1] "kf logLik for b2: -13.2725611269272"
> 
> #######################################################################
> # test ibpf with argument re-use and replacement
> # check that results match after re-use and replacement
> #######################################################################
> 
> set.seed(5)
> b2_ibpf <- ibpf(b2,Np=5,block_size=1,Nbpf=2,
+   rw.sd=b2_rw.sd,
+   cooling.frac=0.5, spat_regression=0.1,
+   unitParNames="rho",sharedParNames=NULL)
> 
> paste("bm2 ibpf loglik: ",round(logLik(b2_ibpf),10))
[1] "bm2 ibpf loglik:  -13.2406696135"
> 
> set.seed(5)
> b2_ibpf_repeat <- ibpf(b2_ibpf,params=coef(b2), unitParNames="rho",
+   sharedParNames=NULL,spat_regression=0.1)
> paste("check ibpf on ibpfd_spatPomp: ",
+   logLik(b2_ibpf_repeat)==logLik(b2_ibpf))
[1] "check ibpf on ibpfd_spatPomp:  TRUE"
> 
> set.seed(5)
> b2_ibpf_bpfilterd <- ibpf(b2_bpfilter,Np=5,block_size=1,Nbpf=2,
+   rw.sd=b2_rw.sd,
+   cooling.frac=0.5, spat_regression=0.1,
+   unitParNames="rho",sharedParNames=NULL)
> paste("check ibpf on bpfilterd_spatPomp: ",
+   logLik(b2_ibpf_bpfilterd)==logLik(b2_ibpf))
[1] "check ibpf on bpfilterd_spatPomp:  TRUE"
> 
> 
> 
> 
> proc.time()
   user  system elapsed 
   1.57    0.42   13.15