R Under development (unstable) (2024-11-13 r87330 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
> set.seed(1)
> 
> ## a default measles model with fixed shared initial value parameters (IVPs)
> ## compiled into the object
> 
> m_U <- 2 
> m_model <- measles(U=m_U)
> 
> range(time(m_model))
[1] 1950 1965
> 
> m <- window(m_model,1950,1955)
> 
> m_params <- c(R0=56.8,amplitude=0.554,gamma=30.4,sigma=28.9,mu=0.02,sigmaSE=0.02,rho=0.488,psi=0.116,g=100,alpha=1,iota=0,cohort=0.1)
> 
> m_pf <- pfilter(m_model,Np=20,params=m_params)
> logLik(m_pf)
[1] -85104.57
> 
> ## setting IVPs
> 
> m2_model <- measles(U=m_U,fixed_ivps=FALSE)
> m2_model <- window(m2_model,1950,1955)
> 
> m2_shared_params <- c(R0=56.8,amplitude=0.554,gamma=30.4,sigma=28.9,mu=0.02,sigmaSE=0.02,rho=0.488,psi=0.116,g=100,alpha=1,iota=0,cohort=0.1)
> m2_unit_ivp_names <- c('S','E','I')
> m2_statenames <- paste0(rep(m2_unit_ivp_names,each=m_U),1:m_U)
> m2_ivp_names <- paste0(m2_statenames,"_0")
> m2_ivps <- rep(c(0.032,0.00005,0.00004),each=m_U)
> names(m2_ivps) <- m2_ivp_names
> m2_params <- c(m2_shared_params,m2_ivps)
> 
> set.seed(1)
> m2_pf <- pfilter(m2_model,Np=20,params=m2_params)
> logLik(m2_pf)
[1] -46672.26
> 
> ##
> ## Note: the measles skeleton is correct only when there is no cohort effect,
> ## i.e., cohort=0.
> ##
> 
> ## 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. 
> ## It also tests igirf (method=moment) with covariates
> 
> m3_params <- m_params
> m3_params["cohort"] <- 0
> m3_igirf_lookahead <- 1
> m3_igirf_ninter <- 2
> m3_igirf_np <- 5
> m3_igirf_nguide <- 5
> m3_igirf_ngirf <- 2
> 
> suppressWarnings(
+   m3_igirf_out <- igirf(m_model, Ngirf = m3_igirf_ngirf,
+     params=m3_params,
+     rw.sd=rw_sd(g=0.02),
+     cooling.type = "geometric",
+     cooling.fraction.50 = 0.5,
+     Np=m3_igirf_np,
+     Ninter = m3_igirf_ninter,
+     lookahead = m3_igirf_lookahead,
+     Nguide = m3_igirf_nguide,
+     kind = 'moment',
+     verbose = FALSE
+   )
+ )
> logLik(m3_igirf_out)
[1] -86380.55
> 
> ## another call to igirf tests igirf (method=bootstrap) with covariates
> suppressWarnings(m3_igirf_bootstrap_out <- igirf(m3_igirf_out,kind='bootstrap'))
> logLik(m3_igirf_bootstrap_out)
[1] -86123.26
> 
> ## test error message
> try(measles(U=1000))
Error in measles(U = 1000) : 
  Require U<=40 since data are only available for 40 cities
> 
> 
> ## --------------------------------------------------------------------
> ## using measles to test spatPomp() covariate replacement functionality
> ## ____________________________________________________________________
> 
> m_covrep <- spatPomp(m,
+   covar=spatPomp::measlesUK[
+     measlesUK$city %in% m@unit_names,
+     c("year","city","pop","births")]
+ )
> for(slt in slotNames(m)) if(!identical(slot(m,slt),slot(m_covrep,slt))) print(slt)
[1] "unit_covarnames"
[1] "covar"
> ## covariates are not the same since this replacement does not deal with
> ## birth delay
> 
> proc.time()
   user  system elapsed 
   5.59    0.35   13.20