more unit tests for spatPomp() and ibpf()

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: spatPomp
 Type: Package
 Title: Inference for Spatiotemporal Partially Observed Markov Processes
-Version: 0.34.3
-Date: 2024-04-19
+Version: 0.34.4
+Date: 2024-04-21
 Authors@R: c(
   person("Kidus", "Asfaw", email = "[email protected]", role = c("aut")),
   person("Edward", "Ionides", email = "[email protected]",role = c("cre","aut")),

R/ibpf.R

@@ -108,7 +108,6 @@ setMethod(
     if (missing(Np)) pStop_(ep, "Np is required")
     if (missing(sharedParNames)) pStop_(ep, "sharedParNames is required")
     if (missing(unitParNames)) pStop_(ep, "unitParNames is required")
-    if (missing(spat_regression)) spat_regression <- numeric()
     if(missing(block_list) && missing(block_size)) {
       pStop_(ep,sQuote("block_list"), " or ", sQuote("block_size"),
         " must be specified to the call")
@@ -117,10 +116,12 @@ setMethod(
       pStop_(ep,"Exactly one of ",sQuote("block_size"), " and ",
         sQuote("block_list"), " should be provided, but not both.")
     }
-    if (missing(spat_regression) && !is.null(sharedParNames)) 
+    if (missing(spat_regression) && !is.null(sharedParNames)){
       pStop_(ep, sQuote("spat_regression"),
        " should be provided when there are shared parameters")
-      if (missing(block_list)){
+    }
+    if (missing(spat_regression)) spat_regression <- numeric()
+    if (missing(block_list)){
       if(block_size > length(unit_names(data))){
         pStop_(ep,sQuote("block_size"), " cannot be greater than the number of spatial units")
       }
@@ -166,15 +167,15 @@ setMethod(
   ){
     ep <- paste0("in ",sQuote("ibpf"),": ")
     if (!missing(block_list) & !missing(block_size)){
-      stop(ep,"Exactly one of ",sQuote("block_size"), " and ", sQuote("block_list"), " can be provided, but not both.",call.=FALSE)
+      pStop_(ep,"Exactly one of ",sQuote("block_size"), " and ", sQuote("block_list"), " can be provided, but not both.")
     }
 
     if(missing(block_list) && missing(block_size))
       block_list <- data@block_list
 
     if (!missing(block_size)){
       if(block_size > length(unit_names(data))){
-        stop(ep,sQuote("block_size"), " cannot be greater than the number of spatial units",call.=FALSE)
+        pStop_(ep,sQuote("block_size"), " cannot be greater than the number of spatial units")
       }
       all_units = seq_len(length(unit_names(data)))
       nblocks = round(length(all_units)/block_size)
@@ -184,10 +185,11 @@ setMethod(
 
     if (missing(Np)) Np <- data@Np
     if (missing(Nbpf)) Nbpf <- data@Nbpf
-    if (!is.numeric(Nbpf)|| Nbpf <1) stop(ep, sQuote("Nbpf"), " should be a positive integer")
+    if (!is.numeric(Nbpf)|| Nbpf <1) pStop_(ep, sQuote("Nbpf"), " should be a positive integer")
     if (missing(rw.sd)) rw.sd <- data@rw.sd
     if (missing(cooling.type)) cooling.type <- data@cooling.type
     if (missing(cooling.fraction.50)) cooling.fraction.50 <- data@cooling.fraction.50
+    if (missing(spat_regression)) spat_regression <- data@spat_regression
 
     tryCatch(
       ibpf_internal(data,Nbpf=Nbpf,spat_regression=spat_regression,
@@ -220,12 +222,12 @@ setMethod(
   ) {
     ep <- paste0("in ",sQuote("ibpf"),": ")
     if (!missing(block_list) & !missing(block_size)){
-      stop(ep,"Exactly one of ",sQuote("block_size"), " and ", sQuote("block_list"), " can be provided, but not both.",call.=FALSE)
+      pStop(ep,"Exactly one of ",sQuote("block_size"), " and ", sQuote("block_list"), " can be provided, but not both.")
     }
     if(missing(block_list) && missing(block_size)) block_list <- data@block_list
     if (!missing(block_size)){
       if(block_size > length(unit_names(data))){
-        stop(ep,sQuote("block_size"), " cannot be greater than the number of spatial units",call.=FALSE)
+        pStop_(ep,sQuote("block_size"), " cannot be greater than the number of spatial units")
       }
       all_units = seq_len(length(unit_names(data)))
       nblocks = round(length(all_units)/block_size)

R/spatPomp.R

@@ -98,10 +98,10 @@ spatPomp <- function (data, units, times, covar, t0, ...,
   ep <- paste0("in ",sQuote("spatPomp"),": ")
 
   if (missing(data))
-    stop(ep,sQuote("data")," is a required argument",call.=FALSE)
+    pStop_(ep,sQuote("data")," is a required argument")
 
   if (!inherits(data,what=c("data.frame","spatPomp")))
-    pStop("spatPomp",sQuote("data")," must be a data frame or an object of ",
+    pStop_(ep, sQuote("data")," must be a data frame or an object of ",
       "class ",sQuote("spatPomp"),".")
 
   ## return as quickly as possible if no work is to be done
@@ -133,7 +133,7 @@ spatPomp <- function (data, units, times, covar, t0, ...,
       globals=globals,cdir=cdir,cfile=cfile,shlib.args=shlib.args,
       compile=compile, verbose=verbose
     ),
-    error = function (e) stop(conditionMessage(e))
+    error = function (e) pStop_(conditionMessage(e))
   )
 }
 
@@ -142,7 +142,7 @@ setMethod(
   "construct_spatPomp",
   signature=signature(data="ANY", times="ANY", units="ANY"),
   definition = function (data, times, t0, ...) {
-    stop(sQuote("times")," should be a single name identifying the column of data that represents",
+    pStop_("in spatPomp : ", sQuote("times")," should be a single name identifying the column of data that represents",
       " the observation times. ", sQuote("units"), " should be likewise for column that represents",
       " the observation units.")
   }
@@ -158,19 +158,21 @@ setMethod(
     paramnames, shared_covarnames, PACKAGE, globals,
     cdir, cfile, shlib.args, compile, verbose) {
 
-    if (anyDuplicated(names(data)))
-      stop("names of data variables must be unique.")
+   ep <- paste0("in ",sQuote("spatPomp"),": ")
 
-    if (missing(t0)) reqd_arg(NULL,"t0")
+   if (anyDuplicated(names(data)))
+      pStop_(ep, "names of data variables must be unique.")
+
+    if (missing(t0)) pStop_(ep, sQuote("t0"), " is a required argument")
 
     tpos <- match(times,names(data),nomatch=0L)
     upos <- match(units,names(data),nomatch=0L)
 
     if (length(times) != 1 || tpos == 0L)
-      stop(sQuote("times")," does not identify a single column of ",
+      pStop_(ep, sQuote("times")," does not identify a single column of ",
         sQuote("data")," by name.")
     if (length(units) != 1 || upos == 0L)
-      stop(sQuote("units")," does not identify a single column of ",
+      pStop_(ep, sQuote("units")," does not identify a single column of ",
         sQuote("data")," by name.")
 
     timename <- times
@@ -229,7 +231,8 @@ setMethod(
     if(!missing(covar)){
       if(timename %in% names(covar)) tcovar <- timename
       else{
-        stop(sQuote("covariate"), ' data.frame should have a time column with the same name as the ',
+        pStop_(ep, sQuote("covariate"),
+	  ' data.frame should have a time column with the same name as the ',
           'time column of the observation data.frame')
       }
     }

tests/bm.R

@@ -341,14 +341,31 @@ spatPomp_Csnippet("lik=u;",unit_statenames="A",unit_obsnames=c("B","C"), unit_co
 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)
 
-# test an error message
-print("The following delivers an error message, to test it")
-try(spatPomp(data=as.data.frame(b_model),units=NULL),outFile=stdout())
-
 # 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')
 
+# test do-nothing behavior
+b_rep3 <- spatPomp(b_model)
+
+## --------------------------------------------
+## using bm to test spatPomp() warning messages
+## ____________________________________________
+
+print("The following deliver error messages, to test them")
+try(spatPomp(data=as.data.frame(b_model),units=NULL),outFile=stdout())
+try(spatPomp("test on type character"))
+
+b_data <- as.data.frame(b_model)
+try(spatPomp(data=b_data,times="time",units="unit"))
+try(spatPomp(data=b_data,times="NONSENSE",units="unit",t0=0))
+try(spatPomp(data=b_data,times="time",units="NONSENSE",t0=0))
+b_data2 <- b_data
+names(b_data2) <- c("time","unit","X","X")
+try(spatPomp(data=b_data2,times="time",units="unit"))
+b_data_only_model <- spatPomp(data=b_data,times="time",units="unit",
+  t0=0)
+
 ## -----------------------------------------------------------------
 ## using bm to test behavior of inference methods when logLik = -Inf
 ## _________________________________________________________________

tests/bm.Rout.save

@@ -1,7 +1,7 @@
 
-R version 4.2.3 (2023-03-15) -- "Shortstop Beagle"
-Copyright (C) 2023 The R Foundation for Statistical Computing
-Platform: x86_64-apple-darwin17.0 (64-bit)
+R version 4.3.3 (2024-02-29) -- "Angel Food Cake"
+Copyright (C) 2024 The R Foundation for Statistical Computing
+Platform: x86_64-apple-darwin20 (64-bit)
 
 R is free software and comes with ABSOLUTELY NO WARRANTY.
 You are welcome to redistribute it under certain conditions.
@@ -522,17 +522,38 @@ Slot "text":
 > 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')
 > 
+> # test do-nothing behavior
+> b_rep3 <- spatPomp(b_model)
+> 
+> ## --------------------------------------------
+> ## using bm to test spatPomp() warning messages
+> ## ____________________________________________
+> 
+> print("The following deliver error messages, to test them")
+[1] "The following deliver error messages, to test them"
+> try(spatPomp(data=as.data.frame(b_model),units=NULL),outFile=stdout())
+Error : in spatPomp : ‘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.
+> try(spatPomp("test on type character"))
+Error : in ‘spatPomp’: ‘data’ must be a data frame or an object of class ‘spatPomp’.
+> 
+> b_data <- as.data.frame(b_model)
+> try(spatPomp(data=b_data,times="time",units="unit"))
+Error : in ‘spatPomp’: ‘t0’ is a required argument
+> try(spatPomp(data=b_data,times="NONSENSE",units="unit",t0=0))
+Error : in ‘spatPomp’: ‘times’ does not identify a single column of ‘data’ by name.
+> try(spatPomp(data=b_data,times="time",units="NONSENSE",t0=0))
+Error : in ‘spatPomp’: ‘units’ does not identify a single column of ‘data’ by name.
+> b_data2 <- b_data
+> names(b_data2) <- c("time","unit","X","X")
+> try(spatPomp(data=b_data2,times="time",units="unit"))
+Error : in ‘spatPomp’: names of data variables must be unique.
+> b_data_only_model <- spatPomp(data=b_data,times="time",units="unit",
++   t0=0)
+> 
 > ## -----------------------------------------------------------------
 > ## using bm to test behavior of inference methods when logLik = -Inf
 > ## _________________________________________________________________

tests/he10.R

@@ -40,12 +40,14 @@ if(model_type == "mostly fixed"){
 # Note: here we assume that there are no unestimated unit-specific
 # parameters. That could readily be accommodated if needed.
 
-h_model <- he10(U=2,dt=4/365,Tmax=1950.5,
+h_model <- he10(U=2,dt=4/365,Tmax=1950.1,
   expandedParNames=estParNames)
 
 coef(h_model)
 
-paste("bpfilter logLik for he10 model:",logLik(bpfilter(h_model,Np=10,block_size=1)))
+h_bpfilter <- bpfilter(h_model,Np=10,block_size=1)
+
+paste("bpfilter logLik for he10 model:",logLik(h_bpfilter))
 
 
 h_U <- length(unit_names(h_model))
@@ -73,6 +75,7 @@ nblocks = 2
 block_list = split(all_units, sort(all_units %% nblocks))
 block_list <- lapply(block_list, as.integer)
 
+set.seed(3)
 h_ibpf <- ibpf(h_model,
   params=coef(h_model),
   sharedParNames=sharedParNames,
@@ -85,37 +88,93 @@ h_ibpf <- ibpf(h_model,
   block_list=block_list
 )
 
-paste("ibpf logLik for he10 model:",logLik(bpfilter(h_ibpf,Np=10,block_size=1)))
+h_bpfilter <- bpfilter(h_ibpf,Np=10,block_size=1)
+
+paste("ibpf logLik for he10 model:",logLik(h_bpfilter))
+
+# test whether specifying Np as a function gives the same result
+set.seed(3)
+h_ibpf2 <- ibpf(
+h_model,
+  params=coef(h_model),
+  sharedParNames=sharedParNames,
+  unitParNames=unitParNames,
+  Nbpf=2,
+  spat_regression=0.1,
+  Np=function(k) 10,
+  rw.sd=h_rw.sd,
+  cooling.fraction.50=0.5,
+  block_list=block_list
+)
+
+h_bpfilter2 <- bpfilter(h_ibpf2,Np=10,block_size=1)
 
+if (logLik(h_bpfilter2)!=logLik(h_bpfilter))
+  stop("in ibpf: Np specified as a function gives a different result from Np as a scalar")
+
 coef(h_ibpf)
 
-# test various errors
+# test errors for ibpf on class 'missing' or character
 try(ibpf())
 try(ibpf("h_model"))
+
+# test errors for ibpf on class spatPomp
+try(ibpf(h_model))
+try(ibpf(h_model,Nbpf=2))
+try(ibpf(h_model,Nbpf=2,rw.sd=rw_sd(mu1=0.1)))
 try(ibpf(h_model,Nbpf=NA,Np=10))
 try(ibpf(h_model,Nbpf=NA,Np=10,block_size=1))
 try(ibpf(h_model,Nbpf=NA,Np=10,block_size=1,sharedParNames=NULL))
+try(ibpf(h_model,Nbpf=2,rw.sd=rw_sd(mu1=0.1),Np=10,sharedParNames=sharedParNames,unitParNames=unitParNames))
+try(ibpf(h_model,Nbpf=2,rw.sd=rw_sd(mu1=0.1),Np=10,sharedParNames=sharedParNames,unitParNames=unitParNames,block_list=block_list,block_size=1))
+try(ibpf(h_model,Nbpf=2,rw.sd=rw_sd(mu1=0.1),Np=10,sharedParNames=sharedParNames,unitParNames=unitParNames,block_list=block_list))
+try(ibpf(h_model,Nbpf=2,rw.sd=rw_sd(mu1=0.1),Np=5,sharedParNames=sharedParNames,unitParNames=unitParNames,spat_regression=0.5,block_size=10))
+
 try(ibpf(h_model,Nbpf=NULL,block_list=block_list,Np=10,rw.sd=rw_sd(mu1=0.1)))
 try(ibpf(h_model,Nbpf=NULL,block_list=block_list,Np=10,rw.sd=rw_sd(mu1=0.1),sharedParNames=NULL))
 try(ibpf(h_model,Nbpf=NULL,block_list=block_list,Np=10,rw.sd=rw_sd(mu1=0.1),sharedParNames=NULL,unitParNames=NULL))
 try(ibpf(h_model,Nbpf=NULL,block_list=block_list,Np=10,rw.sd=rw_sd(mu1=0.1),sharedParNames=sharedParNames,unitParNames=unitParNames,cooling.fraction.50=0.5))
 try(ibpf(h_model,Nbpf=NULL,block_list=block_list,Np=10,rw.sd=rw_sd(mu1=0.00001),sharedParNames=sharedParNames,unitParNames=unitParNames,cooling.fraction.50=0.5,spat_regression=0.5))
 try(ibpf(h_model,Nbpf=1,block_list=block_list,Np=10,rw.sd=rw_sd(mu1=0.00001),sharedParNames=sharedParNames,unitParNames=unitParNames,cooling.fraction.50=12,spat_regression=0.5))
 try(ibpf(h_model,Nbpf=-1,block_list=block_list,Np=10,rw.sd=rw_sd(mu1=0.00001),sharedParNames=sharedParNames,unitParNames=unitParNames,cooling.fraction.50=0.5,spat_regression=0.5))
+
+# test errors on Np specification
 try(ibpf(h_model,Nbpf=2,block_list=block_list,Np=NULL,rw.sd=rw_sd(mu1=0.00001),sharedParNames=sharedParNames,unitParNames=unitParNames,cooling.fraction.50=0.5,spat_regression=0.5))
 try(ibpf(h_model,Nbpf=2,block_list=block_list,Np=1:100,rw.sd=rw_sd(mu1=0.00001),sharedParNames=sharedParNames,unitParNames=unitParNames,cooling.fraction.50=0.5,spat_regression=0.5))
+try(ibpf(h_model,Nbpf=2,block_list=block_list,Np="a character vector",rw.sd=rw_sd(mu1=0.00001),sharedParNames=sharedParNames,unitParNames=unitParNames,cooling.fraction.50=0.5,spat_regression=0.5))
+try(ibpf(h_model,Nbpf=2,block_list=block_list,Np=c(10,10),rw.sd=rw_sd(mu1=0.00001),sharedParNames=sharedParNames,unitParNames=unitParNames,cooling.fraction.50=0.5,spat_regression=0.5))
+
+
+# test ibpf errors on class ibpfd_spatPomp
+ibpf(h_ibpf,sharedParNames=sharedParNames,unitParNames=unitParNames)
+try(ibpf(h_ibpf,sharedParNames=sharedParNames,unitParNames=unitParNames,block_size=1,Nbpf <- 0.1))
+try(ibpf(h_ibpf,sharedParNames=sharedParNames,unitParNames=unitParNames,block_size=3))
 
-## this works, but was tested previously
-## try(ibpf(h_model,Nbpf=1,block_list=block_list,Np=10,rw.sd=rw_sd(mu1=0.00001),sharedParNames=sharedParNames,unitParNames=unitParNames,cooling.fraction.50=0.5,spat_regression=0.5))
+# test ibpf errors on class bpfilterd_spatPomp
+try(ibpf(h_bpfilter,block_list=block_list,block_size=1))
+try(ibpf(h_bpfilter,block_size=23))
+try(ibpf(h_bpfilter))
 
-# tmp <- ibpf(spatPomp(h_model,dmeasure=function(log,...)log(-1)), Nbpf=1,block_list=block_list,Np=10,rw.sd=rw_sd(mu1=0.00001),sharedParNames=sharedParNames,unitParNames=unitParNames,cooling.fraction.50=0.5,spat_regression=0.5)
+# test ibpf with missing basic model component
+h_model2 <- spatPomp(h_model,rprocess=NULL)
+try(h_ibpf2 <- ibpf(h_model2,
+  params=coef(h_model),
+  sharedParNames=sharedParNames,
+  unitParNames=unitParNames,
+  Nbpf=2,
+  spat_regression=0.1,
+  Np=10,
+  rw.sd=h_rw.sd,
+  cooling.fraction.50=0.5,
+  block_list=block_list
+))
 
 ## test error message when munit_measure is undefined
 try(girf(h_model,kind="moment",
   Np=10,Ninter=2,Nguide=10,lookahead=1,tol=1e-5))
 
 # Create second ibpfd_spatPomp object with different chain length, to test error
-h_ibpf2 <- ibpf(h_model,
+h_ibpf3 <- ibpf(h_model,
                 params=coef(h_model),
                 sharedParNames=sharedParNames,
                 unitParNames=unitParNames,
@@ -131,5 +190,5 @@ h_ibpf2 <- ibpf(h_model,
 is(c(h_ibpf, h_ibpf), "ibpfList")
 
 # Throws error because they have different chain lengths
-try(c(h_ibpf, h_ibpf2))
+try(c(h_ibpf, h_ibpf3))