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
> 
> ## test error messages
> try(measles2(U=1000,N=5))
Error in measles2(U = 1000, N = 5) : U <= 40
> try(measles2(U=2,N=1000))
Error in measles2(U = 2, N = 1000) : N <= 391
> 
> i <- 1
> DEBUG=FALSE
> U <- switch(i,2,10,40)
> N <- switch(i,2,50,391)
> Np <- switch(i,5,100,1000)
> 
> m1 <- measles2(U=U,N=N)
> if(DEBUG){
+   plot(simulate(m1),log=T)
+ }
> 
> # test for all parameters expanded, by default
> set.seed(1)
> s1 <- simulate(m1)
> head(obs(s1))
        year
name     [,1] [,2]
  cases1   95  123
  cases2   58   17
> 
> if(DEBUG){
+ par1 <- coef(m1)
+ par1[paste0('iota',1:U)] <- 10
+  plot(simulate(m1,params=par1),ty="l",log=T)
+ }
> 
> # test for all parameters contracted, e.g., for mif2 with
> # all shared parameters
> m2 <- measles2(U=U,N=N,expandedParNames=NULL,
+   contractedParNames=c("R0", "c", "A", "muIR",
+     "muEI", "sigmaSE", "rho", "psi", "g", "S_0", "E_0", "I_0"),
+   simulated=TRUE   ## to test this flag for measles2()
+ )
> 
> set.seed(1)
> s2 <- simulate(m2)
> 
> # test for all parameters contracted, e.g., for mif2 with
> # all shared parameters
> m2 <- measles2(U=U,N=N,expandedParNames=NULL,
+   contractedParNames=c("R0", "c", "A", "muIR",
+     "muEI", "sigmaSE", "rho", "psi", "g", "S_0", "E_0", "I_0")
+ )
> 
> # test for both expanded and contracted parameters
> m3 <- measles2(U=U,N=N,expandedParNames=c("A","muIR"),
+   contractedParNames=c("R0", "c", "S_0")
+ )
> 
> set.seed(1)
> s3 <- simulate(m3)
> 
> if(any(obs(s2)!=obs(s1)))stop("s1 and s2 should be identical")
> if(any(obs(s3)!=obs(s1)))stop("s1 and s3 should be identical")
> 
> partrans(m1,coef(m1),dir="toEst")
        R01         R02          c1          c2          A1          A2 
  3.4011974   3.4011974  -0.8472979  -0.8472979   0.0000000   0.0000000 
      muIR1       muIR2       muEI1       muEI2    sigmaSE1    sigmaSE2 
  3.9512437   3.9512437   3.9512437   3.9512437  -1.8971200  -1.8971200 
       rho1        rho2        psi1        psi2          g1          g2 
  0.0000000   0.0000000  -1.8971200  -1.8971200   5.9914645   5.9914645 
       S_01        S_02        E_01        E_02        I_01        I_02 
 -3.4094962  -3.4094962  -9.9034376  -9.9034376 -10.1265911 -10.1265911 
     alpha1       iota1        muD1 
  0.9800000   0.1000000   0.0200000 
> partrans(m2,coef(m2),dir="toEst")
        R01          c1          A1       muIR1       muEI1    sigmaSE1 
  3.4011974  -0.8472979   0.0000000   3.9512437   3.9512437  -1.8971200 
       rho1        psi1          g1        S_01        E_01        I_01 
  0.0000000  -1.8971200   5.9914645  -3.4094962  -9.9034376 -10.1265911 
     alpha1       iota1        muD1 
  0.9800000   0.1000000   0.0200000 
> partrans(m3,coef(m3),dir="toEst")
         A1          A2       muIR1       muIR2         R01          c1 
  0.0000000   0.0000000   3.9512437   3.9512437   3.4011974  -0.8472979 
       S_01      alpha1       iota1        psi1       muEI1    sigmaSE1 
 -3.4094962   0.9800000   0.1000000   0.1500000  52.0000000   0.1500000 
       muD1        rho1          g1        E_01        I_01 
  0.0200000   0.5000000 400.0000000   0.0000500   0.0000400 
> 
> set.seed(2)
> pf1 <- pfilter(m1,Np=Np)
> logLik(pf1)
[1] -38.565
> set.seed(2)
> pf2 <- pfilter(m2,Np=Np)
> logLik(pf2)
[1] -38.565
> set.seed(2)
> pf3 <- pfilter(m3,Np=Np)
> logLik(pf3)
[1] -38.565
> 
> e1 <- enkf(m1,Np=Np)
> logLik(e1)
[1] -98.95214
> 
> if(DEBUG){
+   # compare to measles()
+   n1 <- measles(U=40)
+   coef(n1) <- c(
+     alpha = 1,
+     iota = 0,  
+     R0 = 30,
+     cohort = 0,
+     amplitude = 0.5,
+     gamma = 52,
+     sigma = 52,
+     mu = 0.02,
+     sigmaSE = 0.15, 
+     rho = 0.5,
+     psi = 0.15,
+     g = 400,
+     S_0 = 0.032, 
+     E_0 = 0.00005, 
+     I_0 = 0.00004
+   )
+   t1 <- simulate(n1)
+   plot(simulate(n1),ty="l",log=T)
+ }
> 
> 
> proc.time()
   user  system elapsed 
   1.62    0.23   17.10