From c65486f256df869e3932a62e85daad9df428df06 Mon Sep 17 00:00:00 2001 From: Ed Ionides Date: Sat, 20 Apr 2024 10:26:58 -0400 Subject: [PATCH] more unit tests for spatPomp() and ibpf() --- DESCRIPTION | 4 +- R/ibpf.R | 18 +++---- R/spatPomp.R | 23 +++++---- tests/bm.R | 25 ++++++++-- tests/bm.Rout.save | 41 ++++++++++++---- tests/he10.R | 77 ++++++++++++++++++++++++++---- tests/he10.Rout.save | 109 ++++++++++++++++++++++++++++++++++++------- 7 files changed, 237 insertions(+), 60 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 2c00426..7ba81da 100644 --- a/DESCRIPTION +++ b/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 = "kidusasfaw1990@gmail.com", role = c("aut")), person("Edward", "Ionides", email = "ionides@umich.edu",role = c("cre","aut")), diff --git a/R/ibpf.R b/R/ibpf.R index 852a3b6..aea798b 100644 --- a/R/ibpf.R +++ b/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,7 +167,7 @@ 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)) @@ -174,7 +175,7 @@ setMethod( 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) diff --git a/R/spatPomp.R b/R/spatPomp.R index 3b1531e..0d11a3f 100644 --- a/R/spatPomp.R +++ b/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') } } diff --git a/tests/bm.R b/tests/bm.R index a8d0417..12c8c7b 100644 --- a/tests/bm.R +++ b/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 ## _________________________________________________________________ diff --git a/tests/bm.Rout.save b/tests/bm.Rout.save index da218d4..3e8256c 100644 --- a/tests/bm.Rout.save +++ b/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 > ## _________________________________________________________________ diff --git a/tests/he10.R b/tests/he10.R index 590583e..480e263 100644 --- a/tests/he10.R +++ b/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,16 +88,48 @@ 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)) @@ -102,20 +137,44 @@ try(ibpf(h_model,Nbpf=NULL,block_list=block_list,Np=10,rw.sd=rw_sd(mu1=0.1),shar 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)) diff --git a/tests/he10.Rout.save b/tests/he10.Rout.save index 7716e1d..1bde922 100644 --- a/tests/he10.Rout.save +++ b/tests/he10.Rout.save @@ -60,7 +60,7 @@ Loading required package: pomp > # 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) @@ -73,8 +73,10 @@ amplitude1 amplitude2 sigmaSE1 sigmaSE2 rho1 rho2 S_01 S_02 E_01 E_02 I_01 I_02 3.2e-02 5.0e-05 5.0e-05 4.0e-05 4.0e-05 > -> paste("bpfilter logLik for he10 model:",logLik(bpfilter(h_model,Np=10,block_size=1))) -[1] "bpfilter logLik for he10 model: -3827.56515459983" +> h_bpfilter <- bpfilter(h_model,Np=10,block_size=1) +> +> paste("bpfilter logLik for he10 model:",logLik(h_bpfilter)) +[1] "bpfilter logLik for he10 model: -2110.11342383345" > > > h_U <- length(unit_names(h_model)) @@ -102,6 +104,7 @@ amplitude1 amplitude2 sigmaSE1 sigmaSE2 rho1 rho2 S_01 > 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, @@ -114,32 +117,71 @@ amplitude1 amplitude2 sigmaSE1 sigmaSE2 rho1 rho2 S_01 + block_list=block_list + ) > -> paste("ibpf logLik for he10 model:",logLik(bpfilter(h_ibpf,Np=10,block_size=1))) -[1] "ibpf logLik for he10 model: -6589.70922701501" +> h_bpfilter <- bpfilter(h_ibpf,Np=10,block_size=1) +> +> paste("ibpf logLik for he10 model:",logLik(h_bpfilter)) +[1] "ibpf logLik for he10 model: -2814.71744861074" > +> # 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) alpha1 iota1 cohort1 mu1 R01 R02 -1.000000e+00 0.000000e+00 0.000000e+00 2.000000e-02 3.194959e+01 3.194959e+01 +1.000000e+00 0.000000e+00 0.000000e+00 2.000000e-02 2.998682e+01 2.998682e+01 psi1 psi2 g1 g2 sigma1 sigma2 -1.489545e-01 1.489545e-01 4.104272e+02 4.104272e+02 4.821791e+01 4.821791e+01 +1.478142e-01 1.478142e-01 3.890605e+02 3.890605e+02 5.146403e+01 5.146403e+01 gamma1 gamma2 amplitude1 amplitude2 sigmaSE1 sigmaSE2 -4.722066e+01 4.722066e+01 4.736288e-01 4.736288e-01 1.461788e-01 1.461788e-01 +5.263311e+01 5.263311e+01 5.011599e-01 5.011599e-01 1.417224e-01 1.417224e-01 rho1 rho2 S_01 S_02 E_01 E_02 -5.423109e-01 5.009793e-01 3.257674e-02 3.188634e-02 4.835306e-05 4.617649e-05 +4.847137e-01 5.354589e-01 3.025523e-02 3.154374e-02 4.923364e-05 5.325866e-05 I_01 I_02 -3.749471e-05 4.049693e-05 +4.123391e-05 4.097914e-05 > -> # test various errors +> # test errors for ibpf on class 'missing' or character > try(ibpf()) Error : in ‘ibpf’: ‘data’ is a required argument. > try(ibpf("h_model")) Error : ‘ibpf’ is undefined for ‘data’ of class ‘character’. +> +> # test errors for ibpf on class spatPomp +> try(ibpf(h_model)) +Error : in ‘ibpf’: Nbpf is required +> try(ibpf(h_model,Nbpf=2)) +Error : in ‘ibpf’: rw.sd is required +> try(ibpf(h_model,Nbpf=2,rw.sd=rw_sd(mu1=0.1))) +Error : in ‘ibpf’: Np is required > try(ibpf(h_model,Nbpf=NA,Np=10)) Error : in ‘ibpf’: rw.sd is required > try(ibpf(h_model,Nbpf=NA,Np=10,block_size=1)) Error : in ‘ibpf’: rw.sd is required > try(ibpf(h_model,Nbpf=NA,Np=10,block_size=1,sharedParNames=NULL)) Error : in ‘ibpf’: rw.sd is required +> try(ibpf(h_model,Nbpf=2,rw.sd=rw_sd(mu1=0.1),Np=10,sharedParNames=sharedParNames,unitParNames=unitParNames)) +Error : in ‘ibpf’: ‘block_list’ or ‘block_size’ must be specified to the call +> 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)) +Error : in ‘ibpf’: Exactly one of ‘block_size’ and ‘block_list’ should be provided, but not both. +> try(ibpf(h_model,Nbpf=2,rw.sd=rw_sd(mu1=0.1),Np=10,sharedParNames=sharedParNames,unitParNames=unitParNames,block_list=block_list)) +Error : in ‘ibpf’: ‘spat_regression’ should be provided when there are shared parameters +> 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)) +Error : in ‘ibpf’: ‘block_size’ cannot be greater than the number of spatial units +> > try(ibpf(h_model,Nbpf=NULL,block_list=block_list,Np=10,rw.sd=rw_sd(mu1=0.1))) Error : in ‘ibpf’: sharedParNames is required > try(ibpf(h_model,Nbpf=NULL,block_list=block_list,Np=10,rw.sd=rw_sd(mu1=0.1),sharedParNames=NULL)) @@ -147,24 +189,57 @@ Error : in ‘ibpf’: unitParNames is required > try(ibpf(h_model,Nbpf=NULL,block_list=block_list,Np=10,rw.sd=rw_sd(mu1=0.1),sharedParNames=NULL,unitParNames=NULL)) Error : in ‘ibpf’: cooling.fraction.50 is required > 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)) -Error : in ‘ibpf’: ‘Nbpf’ must be a positive integer. +Error : in ‘ibpf’: ‘spat_regression’ should be provided when there are shared parameters > 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)) Error : in ‘ibpf’: ‘Nbpf’ must be a positive integer. > 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)) Error : in ‘ibpf’: ‘cooling.fraction.50’ must be in (0,1]. > 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)) Error : in ‘ibpf’: ‘Nbpf’ must be a positive integer. +> +> # 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)) Error : in ‘ibpf’: ‘Np’ must be specified. > 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)) Error : in ‘ibpf’: time-varying ‘Np’ is currently unsupported In addition: Warning message: in ‘ibpf’: Np[k] ignored for k > ‘length(time(object))’. +> 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)) +Error : in ‘ibpf’: ‘Np’ must be a number, a vector of numbers, or a function. +> 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)) +Error : in ‘ibpf’: ‘Np’ must have length 1 or ‘length(time(object))’. +> > -> ## 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 ibpfd_spatPomp +> ibpf(h_ibpf,sharedParNames=sharedParNames,unitParNames=unitParNames) + +> try(ibpf(h_ibpf,sharedParNames=sharedParNames,unitParNames=unitParNames,block_size=1,Nbpf <- 0.1)) +Error : in ‘ibpf’: ‘Nbpf’ should be a positive integer +> try(ibpf(h_ibpf,sharedParNames=sharedParNames,unitParNames=unitParNames,block_size=3)) +Error : in ‘ibpf’: ‘block_size’ cannot be greater than the number of spatial units > -> # 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 errors on class bpfilterd_spatPomp +> try(ibpf(h_bpfilter,block_list=block_list,block_size=1)) +Error : in ‘in ‘ibpf’: ’: Exactly one of ‘block_size’ and ‘block_list’ can be provided, but not both. +> try(ibpf(h_bpfilter,block_size=23)) +Error : in ‘ibpf’: ‘block_size’ cannot be greater than the number of spatial units +> try(ibpf(h_bpfilter)) +Error : in ‘argument "Nbpf" is missing, with no default’: ibpf +> +> # 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 ++ )) +Error : in ‘ibpf’: ‘rprocess’, ‘dunit_measure’ are needed basic components. > > ## test error message when munit_measure is undefined > try(girf(h_model,kind="moment", @@ -173,7 +248,7 @@ Error in .local(object, ...) : girf with kind = 'moment' requires munit_measure > > # 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, @@ -190,7 +265,7 @@ Error in .local(object, ...) : [1] TRUE > > # Throws error because they have different chain lengths -> try(c(h_ibpf, h_ibpf2)) +> try(c(h_ibpf, h_ibpf3)) Error in validObject(.Object) : invalid class “ibpfList” object: error in ‘c’: to be combined, ‘ibpfd_spatPomp’ objects must have chains of equal length >