From 7124dd6d792d12e56e5491f85998f2d1617d6a8a Mon Sep 17 00:00:00 2001 From: Ed Ionides Date: Fri, 19 Apr 2024 15:26:24 -0400 Subject: [PATCH] some more ibpf tests --- DESCRIPTION | 4 +- R/ibpf.R | 88 ++++++++++++++++++++++---------------------- tests/he10.R | 21 +++++++++++ tests/he10.Rout.save | 43 ++++++++++++++++++++-- 4 files changed, 107 insertions(+), 49 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 76f2b7f..2c00426 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: spatPomp Type: Package Title: Inference for Spatiotemporal Partially Observed Markov Processes -Version: 0.34.2 -Date: 2024-02-17 +Version: 0.34.3 +Date: 2024-04-19 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 0bbe6e8..852a3b6 100644 --- a/R/ibpf.R +++ b/R/ibpf.R @@ -103,44 +103,47 @@ setMethod( ..., verbose = getOption("verbose", FALSE) ){ ep <- paste0("in ",sQuote("ibpf"),": ") + if (missing(Nbpf)) pStop_(ep, "Nbpf is required") + if (missing(rw.sd)) pStop_(ep, "rw.sd is required") + 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)) { - stop(ep,sQuote("block_list"), " or ", sQuote("block_size"), - " must be specified to the call",call.=FALSE) + pStop_(ep,sQuote("block_list"), " or ", sQuote("block_size"), + " must be specified to the call") } if (!missing(block_list) & !missing(block_size)){ - stop(ep,"Exactly one of ",sQuote("block_size"), " and ", - sQuote("block_list"), " should be provided, but not both.",call.=FALSE) + 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)) - stop(ep, sQuote("spat_regression"), - " should be provided when there are shared parameters",call.=FALSE) - if (missing(Nbpf)) stop(ep, "Nbpf is required") - if (missing(rw.sd)) stop(ep, "rw.sd is required") - if (missing(Np)) stop(ep, "Np is required") - if (missing(sharedParNames)) stop(ep, "sharedParNames is required") - if (missing(unitParNames)) stop(ep, "unitParNames is required") - if (missing(spat_regression)) spat_regression <- numeric() - if (missing(block_list)){ + pStop_(ep, sQuote("spat_regression"), + " should be provided when there are shared parameters") + if (missing(block_list)){ 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) block_list = split(all_units, sort(all_units %% nblocks)) } block_list <- lapply(block_list, as.integer) + if (missing(cooling.fraction.50)) pStop_(ep, "cooling.fraction.50 is required") tryCatch( - ibpf_internal(data,Nbpf=Nbpf,spat_regression=spat_regression, + ibpf_internal(data, + Nbpf=Nbpf, + spat_regression=spat_regression, Np=Np,rw.sd=rw.sd, sharedParNames=sharedParNames, unitParNames=unitParNames, - cooling.type=cooling.type,cooling.fraction.50=cooling.fraction.50, + cooling.type=cooling.type, + cooling.fraction.50=cooling.fraction.50, block_list=block_list, ...,verbose=verbose ), - error = "error in ibpf_internal" + error = function(e) pStop("ibpf",conditionMessage(e)) ) } ) @@ -181,6 +184,7 @@ 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 (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 @@ -190,10 +194,11 @@ setMethod( Np=Np,rw.sd=rw.sd, sharedParNames=sharedParNames, unitParNames=unitParNames, - cooling.type=cooling.type,cooling.fraction.50=cooling.fraction.50, + cooling.type=cooling.type, + cooling.fraction.50=cooling.fraction.50, block_list=block_list, ...,verbose=verbose), - error = "error in ibpf_internal" + error = function(e) pStop("ibpf",conditionMessage(e)) ) } ) @@ -235,20 +240,17 @@ setMethod( Np=Np,rw.sd=rw.sd, sharedParNames=sharedParNames, unitParNames=unitParNames, - cooling.type=cooling.type,cooling.fraction.50=cooling.fraction.50, + cooling.type=cooling.type, + cooling.fraction.50=cooling.fraction.50, block_list=block_list, .paramMatrix=NULL, ...,verbose=verbose), - error = "error in ibpf_internal" + error = function (e) pStop(who="ibpf",conditionMessage(e)) ) } ) - - - - ibpf_internal <- function (object,Nbpf,spat_regression, Np, rw.sd,cooling.type,cooling.fraction.50, sharedParNames, @@ -308,7 +310,7 @@ ibpf_internal <- function (object,Nbpf,spat_regression, Np <- rep(Np,times=ntimes) } else if (length(Np) > ntimes) { if (Np[1L] != Np[ntimes+1] || length(Np) > ntimes+1) { - pWarn("igirf","Np[k] ignored for k > ",sQuote("length(time(object))"),".") + pWarn("ibpf","Np[k] ignored for k > ",sQuote("length(time(object))"),".") } Np <- head(Np,ntimes) } else if (length(Np) < ntimes) { @@ -413,7 +415,6 @@ ibpf_bpfilter <- function (object,block_list,params, verbose <- as.logical(verbose) nbpf <- as.integer(nbpf) Np <- as.integer(Np) - ep <- paste0("in ",sQuote("ibpf"),": ") do_ta <- length(.indices)>0L if (do_ta && length(.indices)!=Np[1L]) @@ -425,22 +426,22 @@ ibpf_bpfilter <- function (object,block_list,params, nblocks <- length(block_list) loglik <- rep(NA,ntimes) - if (length(Np)==1) + if (length(Np)==1) { Np <- rep(Np,times=ntimes+1) - else if (length(Np)!=(ntimes+1)) - stop(ep,sQuote("Np")," must have length 1 or length ",ntimes+1,call.=FALSE) + } else if (length(Np)!=(ntimes+1)) { + pStop_(sQuote("Np")," must have length 1 or length ",ntimes+1) + } else if (!all(Np==Np[1])) + pStop_("time-varying ", sQuote("Np")," is currently unsupported") if (any(Np<=0)) - stop(ep,"number of particles, ",sQuote("Np"), - ", must always be positive",call.=FALSE) + pStop_("number of particles, ",sQuote("Np"),", must always be positive") if (!is.numeric(Np)) - stop(ep,sQuote("Np"), - " must be a number, a vector of numbers, or a function",call.=FALSE) + pStop_(sQuote("Np"), + " must be a number, a vector of numbers, or a function") Np <- as.integer(Np) if (is.matrix(params)) { if (!all(Np==ncol(params))) - stop(ep,"when ",sQuote("params"), - " is provided as a matrix, do not specify ", - sQuote("Np"),"!",call.=FALSE) + pStop_("when ",sQuote("params"), + " is provided as a matrix, do not specify ", sQuote("Np"),"!") } if (NCOL(params)==1) { coef(object) <- params @@ -448,7 +449,7 @@ ibpf_bpfilter <- function (object,block_list,params, } paramnames <- rownames(params) if (is.null(paramnames)) - stop(ep,sQuote("params")," must have rownames",call.=FALSE) + pStop_(sQuote("params")," must have rownames") # create array to store weights per particle per block_list weights <- array(data = numeric(0), dim=c(nblocks,Np[1L])) @@ -489,8 +490,7 @@ ibpf_bpfilter <- function (object,block_list,params, .gnsi=gnsi ), error = function (e) { - stop(ep,"process simulation error: ", - conditionMessage(e),call.=FALSE) + pStop_("process simulation error: ", conditionMessage(e)) } ) @@ -512,8 +512,8 @@ ibpf_bpfilter <- function (object,block_list,params, .gnsi=gnsi ), error = function (e) { - stop(ep,"error in calculation of weights: ", - conditionMessage(e),call.=FALSE) + pStop_("error in calculation of weights: ", + conditionMessage(e)) } ) log_d <- apply(log_vd[,,1,drop=FALSE], 2, function(x) sum(x)) @@ -549,7 +549,7 @@ ibpf_bpfilter <- function (object,block_list,params, weights=weights[i,] ), error = function (e) { - stop(ep,conditionMessage(e),call.=FALSE) # nocov + pStop_(conditionMessage(e)) } ) diff --git a/tests/he10.R b/tests/he10.R index fdc781e..590583e 100644 --- a/tests/he10.R +++ b/tests/he10.R @@ -89,6 +89,27 @@ paste("ibpf logLik for he10 model:",logLik(bpfilter(h_ibpf,Np=10,block_size=1))) coef(h_ibpf) +# test various errors +try(ibpf()) +try(ibpf("h_model")) +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=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)) +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)) + +## 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)) + +# 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 error message when munit_measure is undefined try(girf(h_model,kind="moment", Np=10,Ninter=2,Nguide=10,lookahead=1,tol=1e-5)) diff --git a/tests/he10.Rout.save b/tests/he10.Rout.save index 814be77..7716e1d 100644 --- a/tests/he10.Rout.save +++ b/tests/he10.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. @@ -129,6 +129,43 @@ amplitude1 amplitude2 sigmaSE1 sigmaSE2 rho1 rho2 S_01 I_01 I_02 3.749471e-05 4.049693e-05 > +> # test various errors +> try(ibpf()) +Error : in ‘ibpf’: ‘data’ is a required argument. +> try(ibpf("h_model")) +Error : ‘ibpf’ is undefined for ‘data’ of class ‘character’. +> 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=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)) +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. +> 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. +> 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))’. +> +> ## 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)) +> +> # 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 error message when munit_measure is undefined > try(girf(h_model,kind="moment", + Np=10,Ninter=2,Nguide=10,lookahead=1,tol=1e-5))