Skip to content

Commit

Permalink
creeping toward full coverage
Browse files Browse the repository at this point in the history
  • Loading branch information
ionides committed May 13, 2024
1 parent 5b4558f commit bd81c81
Show file tree
Hide file tree
Showing 20 changed files with 104 additions and 60 deletions.
15 changes: 7 additions & 8 deletions R/abf.R
Original file line number Diff line number Diff line change
Expand Up @@ -108,7 +108,7 @@ setMethod(
"abf",
signature=signature(object="spatPomp"),
function (object, Nrep, Np, nbhd, params,
tol = 1e-300,
tol = 1.0e-100,
..., verbose=getOption("verbose",FALSE)) {

ep <- paste0("in ",sQuote("abf"),": ")
Expand Down Expand Up @@ -222,9 +222,6 @@ abf_internal <- function (object, Np, nbhd, tol, ..., verbose, .gnsi = TRUE) {
params <- coef(object)
pompLoad(object,verbose=FALSE)
gnsi <- as.logical(.gnsi)
if (length(params)==0) pStop_(ep,sQuote("params")," must be specified")
if (missing(tol)) pStop(ep,sQuote("tol")," must be specified")

times <- time(object,t0=TRUE)
ntimes <- length(times)-1
nunits <- length(unit_names(object))
Expand All @@ -246,10 +243,12 @@ abf_internal <- function (object, Np, nbhd, tol, ..., verbose, .gnsi = TRUE) {
if (!is.numeric(Np))
pStop(ep,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)))
pStop_(ep,"when ",sQuote("params")," is provided as a matrix, do not specify ",sQuote("Np"))
}
## matrix parameters not currently supported
## if (is.matrix(params)) {
## if (!all(Np==ncol(params)))
## pStop_(ep,"when ",sQuote("params"),
## " is provided as a matrix, do not specify ",sQuote("Np"))
## }

params <- as.matrix(params)

Expand Down
15 changes: 3 additions & 12 deletions R/abfir.R
Original file line number Diff line number Diff line change
Expand Up @@ -91,9 +91,8 @@ setMethod(
"abfir",
signature=signature(object="spatPomp"),
function (object, Np, Nrep, nbhd,
Ninter, tol = (1e-300), params, ...,
Ninter, tol = (1e-100), params, ...,
verbose=getOption("verbose",FALSE) ) {

ep <- paste0("in ",sQuote("abfir"),": ")
## declare global variable since foreach's u uses non-standard evaluation
i <- 1
Expand Down Expand Up @@ -205,15 +204,11 @@ abfir_internal <- function (object, Np, nbhd,
unit_obsnames = object@unit_obsnames,
unit_accumvars = object@unit_accumvars)
params <- coef(object)
verbose = FALSE
pompLoad(object,verbose=verbose)
on.exit(pompUnload(object,verbose=verbose))
gnsi <- as.logical(.gnsi)

if (length(params)==0) pStop_(ep,sQuote("params")," must be specified")
if (missing(tol)) pStop_(ep,sQuote("tol")," must be specified")

times <- time(object,t0=TRUE)
times <- time(object,t0=TRUE)
N <- length(times)-1
U <- length(unit_names(object))

Expand Down Expand Up @@ -278,10 +273,7 @@ abfir_internal <- function (object, Np, nbhd,
log=TRUE,
.gnsi=gnsi
),
error = function (e) {
pStop_("abfir error in calculation of weights: ",
conditionMessage(e))
}
error = function (e) pStop_("abfir error in calculation of weights: ", conditionMessage(e))
)

log_cond_densities[,,n] <- log_weights[,,1]
Expand Down Expand Up @@ -375,7 +367,6 @@ abfir_internal <- function (object, Np, nbhd,
}
## resample down to one particle, making Np copies of, say, particle #1.
xas <- xf[,1]

if (verbose) cat("abfir timestep",n,"of",N,"finished\n")

}
Expand Down
2 changes: 1 addition & 1 deletion R/girf.R
Original file line number Diff line number Diff line change
Expand Up @@ -129,7 +129,7 @@ setMethod(
...,
verbose = getOption("verbose", FALSE)) {

if (missing(tol)) tol <- 1e-300
if (missing(tol)) tol <- 1e-100
if (missing(Ninter)) Ninter <- length(unit_names(object))
kind = match.arg(kind)

Expand Down
26 changes: 13 additions & 13 deletions R/igirf.R
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,7 @@ setMethod(
signature=signature(data="spatPomp"),
definition=function (data,Ngirf,Np,rw.sd,cooling.type,cooling.fraction.50,
Ninter,lookahead=1,Nguide,kind=c('bootstrap', 'moment'),
tol = 1e-300,
tol = 1e-100,
..., verbose = getOption("verbose", FALSE)) {

ep <- paste0("in ", sQuote("igirf") , " : ")
Expand Down Expand Up @@ -276,7 +276,7 @@ igirf.momgirf <- function (object, params, Ninter, lookahead, Nguide,
gnsi <- as.logical(.gnsi)
verbose <- as.logical(verbose)
girfiter <- as.integer(girfiter)
ep <- paste0("in ",sQuote("igirf")," : ")
ep <- paste0("in ",sQuote("igirf"), " (method=moment) : ")

if (length(tol) != 1 || !is.finite(tol) || tol < 0)
pStop_(ep, sQuote("tol")," should be a small positive number.")
Expand Down Expand Up @@ -333,7 +333,7 @@ igirf.momgirf <- function (object, params, Ninter, lookahead, Nguide,
times=times[(nt+2):(nt+1+lookahead_steps)],
params=tp_with_guides,
gnsi=TRUE),
error = function (e) pStop_("in ", sQuote("igirf"), " : ", conditionMessage(e))
error = function (e) pStop_("in ", sQuote("igirf_moment"), " : ", conditionMessage(e))
)
fcst_samp_var <- xx
dim(fcst_samp_var) <- c(length(unit_names(object)), lookahead_steps, Np)
Expand Down Expand Up @@ -371,7 +371,7 @@ igirf.momgirf <- function (object, params, Ninter, lookahead, Nguide,
times=times[(nt+2):(nt+1+lookahead_steps)],
params=tparams,
gnsi=TRUE),
error = function (e) pStop_("in ", sQuote("igirf"), " : ",conditionMessage(e))
error = function (e) pStop_("in ", sQuote("igirf_moment"), " : ",conditionMessage(e))
)
dim(meas_var_skel) <- c(length(unit_names(object)), lookahead_steps, Np)
fcst_var_upd <- array(0, dim = c(length(unit_names(object)), lookahead_steps, Np))
Expand All @@ -388,7 +388,7 @@ igirf.momgirf <- function (object, params, Ninter, lookahead, Nguide,
times=times[(nt+2):(nt+1+lookahead_steps)],
params=array.tparams,
gnsi=TRUE),
error = function (e) pStop_("in ", sQuote("igirf"), " : ", conditionMessage(e))
error = function (e) pStop_("in ", sQuote("igirf_moment"), " : ", conditionMessage(e))
)
mom_match_param <- mmp
dim(mom_match_param) <- c(dim(tparams)[1], length(unit_names(object)), lookahead_steps, Np)
Expand All @@ -409,7 +409,7 @@ igirf.momgirf <- function (object, params, Ninter, lookahead, Nguide,
log=TRUE,
.gnsi=gnsi
)),
error = function (e) pStop_("in igirf : error in calculation of log_dmeas_weights: ", conditionMessage(e))
error = function (e) pStop_("in igirf_moment : error in calculation of log_dmeas_weights: ", conditionMessage(e))
)

## uncomment for debugging: tracking down particles with NA weight
Expand Down Expand Up @@ -446,7 +446,7 @@ igirf.momgirf <- function (object, params, Ninter, lookahead, Nguide,
log=TRUE,
.gnsi=gnsi
),
error = function (e) pStop_("in igirf : error in calculation of log_meas_weights: ", conditionMessage(e))
error = function (e) pStop_("in igirf_moment : error in calculation of log_meas_weights: ", conditionMessage(e))
)
gnsi <- FALSE
log_weights <- as.numeric(log_meas_weights) + log_s_not_1_weights
Expand All @@ -455,7 +455,7 @@ igirf.momgirf <- function (object, params, Ninter, lookahead, Nguide,
if (any(log_weights>-Inf)) {
coef(object,transform=TRUE) <- apply(params,1L,weighted.mean,w=exp(log_weights))
} else {
warning("igirf: ","filtering failure at last filter iteration; using ",
warning("igirf_moment: ","filtering failure at last filter iteration; using ",
"unweighted mean for point estimate.")
coef(object,transform=TRUE) <- apply(params,1L,mean)
}
Expand All @@ -476,7 +476,7 @@ igirf.momgirf <- function (object, params, Ninter, lookahead, Nguide,
fsv=fcst_samp_var,
tol=tol
),
error = function (e) pStop_("in igirf : ", conditionMessage(e))
error = function (e) pStop_("in igirf_moment : ", conditionMessage(e))
)
cond.loglik[nt+1, s] <- xx$loglik + max_log_weights
x <- xx$states
Expand Down Expand Up @@ -514,7 +514,7 @@ igirf.bootgirf <- function (object, params, Ninter, lookahead, Nguide,
verbose <- as.logical(verbose)
girfiter <- as.integer(girfiter)
Np <- as.integer(Np)
ep <- paste0("in ",sQuote("ibootgirf"),": ")
ep <- paste0("in ",sQuote("igirf")," (method=bootstrap) : ")

if (length(tol) != 1 || !is.finite(tol) || tol < 0)
stop(ep,sQuote("tol")," should be a small positive number.",call.=FALSE)
Expand Down Expand Up @@ -623,7 +623,7 @@ igirf.bootgirf <- function (object, params, Ninter, lookahead, Nguide,
log=TRUE,
.gnsi=gnsi
)),
error = function (e) pStop_("error in calculation of log_dmeas_weights: ", conditionMessage(e))
error = function (e) pStop_("error in igirf_bootstrap calculation of log_dmeas_weights: ", conditionMessage(e))
)

## uncomment for debugging: tracking down particles with NA weight
Expand Down Expand Up @@ -668,7 +668,7 @@ igirf.bootgirf <- function (object, params, Ninter, lookahead, Nguide,
log=TRUE,
.gnsi=gnsi
),
error = function (e) pStop_("error in calculation of log_meas_weights: ", conditionMessage(e))
error = function (e) pStop_("error in igirf_bootstrap calculation of log_meas_weights: ", conditionMessage(e))
)
gnsi <- FALSE
log_weights <- as.numeric(log_meas_weights) + log_s_not_1_weights
Expand Down Expand Up @@ -699,7 +699,7 @@ igirf.bootgirf <- function (object, params, Ninter, lookahead, Nguide,
fsv=array(0,dim=c(length(unit_names(object)), lookahead_steps, Np)),
tol=tol
),
error = function (e) pStop_("error in igirf computations : ", conditionMessage(e))
error = function (e) pStop_("error in igirf_bootsrap computations : ", conditionMessage(e))
)
guidesim_index = guidesim_index[xx$ancestry]
cond.loglik[nt+1, s] <- xx$loglik + max_log_weights
Expand Down
2 changes: 1 addition & 1 deletion man/abf.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

2 changes: 1 addition & 1 deletion man/abfir.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

2 changes: 1 addition & 1 deletion man/igirf.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

5 changes: 0 additions & 5 deletions src/fcstsampvar.c
Original file line number Diff line number Diff line change
Expand Up @@ -80,11 +80,6 @@ SEXP do_fcst_samp_var (SEXP object, SEXP X, SEXP Np, SEXP times, SEXP params, SE
PROTECT(F = ret_array(nunits, nreps, ntimes)); nprotect++;
switch (mode) {

case Rfun: {
}

break;

case native: case regNative: {
int *oidx, *sidx, *pidx, *cidx;
spatPomp_unit_measure_mean *ff = NULL;
Expand Down
21 changes: 20 additions & 1 deletion tests/bm.R
Original file line number Diff line number Diff line change
Expand Up @@ -137,7 +137,9 @@ abf(b_model_zero_dmeasure,Nrep=2,Np=Np,verbose=TRUE)
b_abfir <- abfir(b_model, Nrep = 2, Np = Np, nbhd = b_bag_nbhd)
paste("bm abfir loglik: ",round(logLik(b_abfir),10))

abfir(b_abfir,verbose=TRUE,accumvars="X1")
capture.output(
abfir(b_abfir,verbose=TRUE,accumvars="X1")
) -> b_abfir_out
abfir(b_model, Nrep = 2, Np = Np)
try(abfir(b_model))
try(abfir(b_model,Nrep=2))
Expand All @@ -146,6 +148,10 @@ try(abfir(b_model,Nrep=2,Np=1:10))
try(abfir(b_model,Nrep=2,Np=Np,params=unname(coef(b_model))))
try(abfir(b_model_zero_dmeasure,Nrep = 3, Np = Np))

# test abfir when all particles fail...
# make this happen by setting a high tol
abfir(b_abfir,tol=1000)

## --------------------------------------------------------------
## bpfilter tested on bm
## ______________________________________________________________
Expand Down Expand Up @@ -274,6 +280,11 @@ girf(b_model_with_accumvars,Np = 3,lookahead = 2, Nguide = 3,
girf(b_model_with_accumvars,Np = 3,lookahead = 2, Nguide = 3,
kind = 'boot',Ninter=2)

# test girf when all particles fail...
# make this happen by setting a high tol
girf(b_girf_mom,tol=1000)



## ------------------------------------------------------------
## Now, we test the inference methods
Expand Down Expand Up @@ -397,6 +408,7 @@ igirf(b_igirf_boot_geom,Np=3,
try(igirf(b_igirf_boot_geom,Np=function(x) 4))
try(igirf(b_igirf_boot_geom,Np=function(x) "JUNK"))
try(igirf(b_igirf_boot_geom,Np=5:15))
try(igirf(b_igirf_boot_geom,tol=-1))


igirf(b_model_with_accumvars,kind='moment', Ngirf = 2, Nguide=2,
Expand Down Expand Up @@ -529,6 +541,7 @@ b_sim3v2@data <- exp(b_sim3v2@data)
plot(b_sim3v2,type="l",log=TRUE)
plot(b_sim3[[2]],type="h",plot_unit_names=FALSE)
dev.off()
plot(b_sim3[[2]],type="h",plot_unit_names=TRUE) -> b_sim3_plot

print(b_model)

Expand Down Expand Up @@ -645,6 +658,11 @@ b_array.params2 <- array(c(b_p,b_p),
dimnames = list(params = rownames(b_p)))
try(munit_measure(b_model, x=b_s2, vc=b_vc, Np=3, unit = 2, time=1,params=b_array.params2))

## trigger error messages in fcstsampvar.c
try(.Call("do_fcst_samp_var",b_model,b_s,3,1:10,b_array.params,FALSE))
try(.Call("do_fcst_samp_var",b_model,b_s2,3,1,b_array.params2,FALSE))


## test spatPomp_Csnippet variable construction
spatPomp_Csnippet("lik=u;",unit_statenames="A",unit_obsnames=c("B","C"),
unit_covarnames="D",
Expand Down Expand Up @@ -699,6 +717,7 @@ try(spatPomp(data=b_data,times="time",units="unit",t0=0,
b_shared_covar <- data.frame(time=0:2,Z=3:5)
model_shared_covar <- spatPomp(data=b_data,times="time",units="unit",
t0=0,covar=b_shared_covar, shared_covarnames="Z")
as.data.frame(b_shared_covar)

b_unit_covar <- data.frame(time=c(0:2,0:2),unit=rep(c("U1","U2"),each=3),
Z=rep(3:5,times=2))
Expand Down
Loading

0 comments on commit bd81c81

Please sign in to comment.