From d524fed8ff4e6f7487d6873a8afb21c8b9246ed0 Mon Sep 17 00:00:00 2001 From: Ed Ionides Date: Mon, 29 Apr 2024 10:26:59 -0400 Subject: [PATCH] tidying igirf --- DESCRIPTION | 4 +- R/bpfilter.R | 2 - R/igirf.R | 140 +++++++++++++++++++++------------------------ tests/bm.R | 50 ++++++++++++---- tests/bm.Rout.save | 89 ++++++++++++++++++++++------ 5 files changed, 180 insertions(+), 105 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 7ba81da..d7191e2 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: spatPomp Type: Package Title: Inference for Spatiotemporal Partially Observed Markov Processes -Version: 0.34.4 -Date: 2024-04-21 +Version: 0.34.5 +Date: 2024-04-29 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/bpfilter.R b/R/bpfilter.R index cf9fd44..e811e75 100644 --- a/R/bpfilter.R +++ b/R/bpfilter.R @@ -120,7 +120,6 @@ setMethod( } ) - ##' @name bpfilter-spatPomp ##' @aliases bpfilter,spatPomp-method ##' @rdname bpfilter @@ -245,7 +244,6 @@ bpfilter.internal <- function (object, Np, block_list, save_states, filter_traj, pStop_(ep,sQuote("Np")," must be a number, a vector of numbers, or a function") Np <- as.integer(Np) params <- as.matrix(coef(object)) - paramnames <- rownames(params) ## returns an nvars by nsim matrix init.x <- rinit(object,params=params,nsim=Np[1L],.gnsi=gnsi) diff --git a/R/igirf.R b/R/igirf.R index 5621ad4..ab7aa49 100644 --- a/R/igirf.R +++ b/R/igirf.R @@ -84,7 +84,7 @@ setMethod( "igirf", signature=signature(data="ANY"), definition=function (data, ...) { - pStop_("igirf is undefined for ", sQuote(data), "of class ", sQuote(class(data)), ".") + pStop_("igirf is undefined for ", sQuote("data"), "of class ", sQuote(class(data)), ".") } ) @@ -99,14 +99,16 @@ setMethod( Ninter,lookahead=1,Nguide,kind=c('bootstrap', 'moment'), tol = 1e-300, ..., verbose = getOption("verbose", FALSE)) { + ep <- paste0("in ", sQuote("igirf") , " : ") if(missing(Ninter)) Ninter <- length(unit_names(data)) if(missing(Ngirf)) pStop_(ep, sQuote("Ngirf"), " is a required argument") if (missing(rw.sd)) pStop_(ep, sQuote("rw.sd")," must be specified.") if (missing(cooling.fraction.50)) pStop_(ep, sQuote("cooling.fraction.50")," is a required argument.") + if (missing(cooling.type)) + pStop_(ep, sQuote("cooling.type")," is a required argument.") if (missing(Np)) pStop_(ep, sQuote("Np")," must be specified.") - if (missing(Ninter)) pStop_(ep, sQuote("Ninter")," must be specified.") if (missing(Nguide)) pStop_(ep, sQuote("Nguide")," must be specified.") kind = match.arg(kind) tryCatch( @@ -179,39 +181,15 @@ if (undefined(object@rprocess) || undefined(object@dmeasure)) start <- apply(.paramMatrix,1L,mean) } - ntimes <- length(time(object)) - - if (is.null(Np)) { + if (missing(Np) || is.null(Np)) { pStop_(sQuote("Np")," must be specified.") - } else if (is.function(Np)) { - Np <- tryCatch( - vapply(seq_len(ntimes),Np,numeric(1)), - error = function (e) { - pStop_("if ",sQuote("Np"), - " is a function, it must return a single positive integer.") - } - ) - } else if (!is.numeric(Np)) { - pStop_(sQuote("Np")," must be a number, a vector of numbers, or a function.") - } - - if (length(Np) == 1) { - Np <- rep(Np,times=ntimes) - } else if (length(Np) > ntimes) { - if (Np[1L] != Np[ntimes+1] || length(Np) > ntimes+1) { - warning("in igirf ","Np[k] ignored for k > ",sQuote("length(time(object))"),".") - } - Np <- head(Np,ntimes) - } else if (length(Np) < ntimes) { - pStop_(sQuote("Np")," must have length 1 or ", - sQuote("length(time(object))"),".") + } else if (length(Np)>1) { + pStop_(sQuote("Np")," must be a single positive integer") + } else if (!is.numeric(Np)|| !is.finite(Np) || (Np < 0)) { + pStop_(sQuote("Np")," must be a single positive integer") } - if (!all(is.finite(Np)) || any(Np <= 0)) - pStop_(sQuote("Np")," must be a positive integer.") - Np <- as.integer(Np) - Np <- c(Np,Np[1L]) ##rw.sd <- perturbn.kernel.sd(rw.sd,time=time(object),paramnames=names(start)) rw.sd <- perturbn.kernel.sd(rw.sd, @@ -230,7 +208,7 @@ if (undefined(object@rprocess) || undefined(object@dmeasure)) ) if (is.null(.paramMatrix)) { - paramMatrix <- array(data=start,dim=c(length(start),Np[1L]), + paramMatrix <- array(data=start,dim=c(length(start),Np), dimnames=list(variable=names(start),rep=NULL)) } else { paramMatrix <- .paramMatrix @@ -298,14 +276,13 @@ igirf.momgirf <- function (object, params, Ninter, lookahead, Nguide, gnsi <- as.logical(.gnsi) verbose <- as.logical(verbose) girfiter <- as.integer(girfiter) - Np <- as.integer(Np) ep <- paste0("in ",sQuote("igirf")," : ") if (length(tol) != 1 || !is.finite(tol) || tol < 0) pStop_(ep, sQuote("tol")," should be a small positive number.") do_ta <- length(.indices)>0L - if (do_ta && length(.indices)!=Np[1L]) + if (do_ta && length(.indices)!=Np) pStop_(ep, sQuote(".indices")," has improper length.") times <- time(object,t0=TRUE) @@ -317,7 +294,7 @@ igirf.momgirf <- function (object, params, Ninter, lookahead, Nguide, znames <- object@accumvars # Initialize filter guide function - log_filter_guide_fun <- array(0, dim = Np[1]) + log_filter_guide_fun <- array(0, dim = Np) for (nt in 0:(ntimes-1)) { pmag <- cooling.fn(nt,girfiter)$alpha*rw.sd[,nt] params <- .Call(randwalk_perturbation_spatPomp,params,pmag) @@ -331,7 +308,7 @@ igirf.momgirf <- function (object, params, Ninter, lookahead, Nguide, lookahead_steps <- min(lookahead, ntimes-nt) ## Four-dimensional array: nvars by nguide by ntimes by nreps Xg <- array(0, - dim=c(length(statenames), Nguide, lookahead_steps, Np[1]), + dim=c(length(statenames), Nguide, lookahead_steps, Np), dimnames = list( nvars = statenames, ng = NULL, @@ -341,16 +318,17 @@ igirf.momgirf <- function (object, params, Ninter, lookahead, Nguide, ) # For each particle get K guide particles, and fill in sample variance # over K for each (lookahead value - unit - particle) combination - fcst_samp_var <- array(0, dim = c(length(unit_names(object)), lookahead_steps, Np[1])) - x_with_guides <- x[,rep(1:Np[1], rep(Nguide, Np[1]))] - tp_with_guides <- tparams[,rep(1:Np[1], rep(Nguide, Np[1]))] - Xg <- rprocess(object, x0=x_with_guides, t0=times[nt+1], times=times[(nt+2):(nt+1+lookahead_steps)], + fcst_samp_var <- array(0, dim = c(length(unit_names(object)), lookahead_steps, Np)) + x_with_guides <- x[,rep(1:Np, rep(Nguide, Np))] + tp_with_guides <- tparams[,rep(1:Np, rep(Nguide, Np))] + Xg <- rprocess(object, x0=x_with_guides, t0=times[nt+1], + times=times[(nt+2):(nt+1+lookahead_steps)], params=tp_with_guides,.gnsi=gnsi) xx <- tryCatch( .Call(do_fcst_samp_var, object=object, X=Xg, - Np = as.integer(Np[1]), + Np = as.integer(Np), times=times[(nt+2):(nt+1+lookahead_steps)], params=tp_with_guides, gnsi=TRUE), @@ -359,7 +337,7 @@ igirf.momgirf <- function (object, params, Ninter, lookahead, Nguide, } ) fcst_samp_var <- xx - dim(fcst_samp_var) <- c(length(unit_names(object)), lookahead_steps, Np[1]) + dim(fcst_samp_var) <- c(length(unit_names(object)), lookahead_steps, Np) for (s in 1:Ninter){ tparams <- partrans(object,params,dir="fromEst",.gnsi=gnsi) @@ -390,7 +368,7 @@ igirf.momgirf <- function (object, params, Ninter, lookahead, Nguide, .Call(do_vunit_measure, object=object, X=skel, - Np = as.integer(Np[1]), + Np = as.integer(Np), times=times[(nt+2):(nt+1+lookahead_steps)], params=tparams, gnsi=TRUE), @@ -398,27 +376,27 @@ igirf.momgirf <- function (object, params, Ninter, lookahead, Nguide, pStop_("in ", sQuote("igirf"), " : ",conditionMessage(e)) } ) - dim(meas_var_skel) <- c(length(unit_names(object)), lookahead_steps, Np[1]) - fcst_var_upd <- array(0, dim = c(length(unit_names(object)), lookahead_steps, Np[1])) + 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)) for(l in 1:lookahead_steps) fcst_var_upd[,l,] <- fcst_samp_var[,l,]*(times[nt+1+l] - tt[s+1])/(times[nt+1+l] - times[nt+1]) inflated_var <- meas_var_skel + fcst_var_upd - array.tparams <- array(NA, dim = c(dim(tparams)[1], length(unit_names(object)), Np[1], lookahead_steps), dimnames = list(tparams = rownames(tparams))) + array.tparams <- array(NA, dim = c(dim(tparams)[1], length(unit_names(object)), Np, lookahead_steps), dimnames = list(tparams = rownames(tparams))) for(i in 1:length(unit_names(object))) array.tparams[,i,,1:lookahead_steps] <- tparams mmp <- tryCatch( .Call(do_munit_measure, object=object, X=skel, vc=inflated_var, - Np = as.integer(Np[1]), + Np = as.integer(Np), times=times[(nt+2):(nt+1+lookahead_steps)], params=array.tparams, gnsi=TRUE), error = function (e) pStop_("in ", sQuote("igirf"), " : ", conditionMessage(e)) ) mom_match_param <- mmp - dim(mom_match_param) <- c(dim(tparams)[1], length(unit_names(object)), lookahead_steps, Np[1]) + dim(mom_match_param) <- c(dim(tparams)[1], length(unit_names(object)), lookahead_steps, Np) dimnames(mom_match_param) <- list(tparam = rownames(tparams)) - log_guide_fun <- vector(mode = "numeric", length = Np[1]) + 1 + log_guide_fun <- vector(mode = "numeric", length = Np) + 1 for(l in 1:lookahead_steps){ if(nt+1+l-lookahead <= 0) discount_denom_init <- object@t0 @@ -493,7 +471,7 @@ igirf.momgirf <- function (object, params, Ninter, lookahead, Nguide, .Call(girf_computations, x=X, params=params, - Np=Np[nt+1], + Np=Np, trackancestry=FALSE, doparRS=TRUE, weights=weights, @@ -524,7 +502,7 @@ igirf.momgirf <- function (object, params, Ninter, lookahead, Nguide, Nguide=Nguide, lookahead=lookahead, paramMatrix=params, - Np=Np[1], + Np=Np, tol=tol, loglik=sum(cond.loglik) ) @@ -545,7 +523,7 @@ igirf.bootgirf <- function (object, params, Ninter, lookahead, Nguide, stop(ep,sQuote("tol")," should be a small positive number.",call.=FALSE) do_ta <- length(.indices)>0L - if (do_ta && length(.indices)!=Np[1L]) + if (do_ta && length(.indices)!=Np) stop(ep,sQuote(".indices")," has improper length.",call.=FALSE) times <- time(object,t0=TRUE) @@ -556,7 +534,7 @@ igirf.bootgirf <- function (object, params, Ninter, lookahead, Nguide, znames <- object@accumvars - log_filter_guide_fun <- array(0, dim = Np[1]) + log_filter_guide_fun <- array(0, dim = Np) for (nt in 0:(ntimes-1)) { pmag <- cooling.fn(nt,girfiter)$alpha*rw.sd[,nt] params <- .Call(randwalk_perturbation_spatPomp,params,pmag) @@ -569,7 +547,7 @@ igirf.bootgirf <- function (object, params, Ninter, lookahead, Nguide, lookahead_steps <- min(lookahead, ntimes-nt) ## Four-dimensional array: nvars by nguide by ntimes by nreps Xg <- array(0, - dim=c(length(statenames), Nguide, lookahead_steps, Np[1]), + dim=c(length(statenames), Nguide, lookahead_steps, Np), dimnames = list( nvars = statenames, ng = NULL, @@ -577,7 +555,7 @@ igirf.bootgirf <- function (object, params, Ninter, lookahead, Nguide, nreps = NULL ) ) - guidesim_index <- seq_len(Np[1]) ## the index for guide simulations (to be updated each time resampling occurs) + guidesim_index <- seq_len(Np) ## the index for guide simulations (to be updated each time resampling occurs) x_with_guides <- x[,rep(guidesim_index, each = Nguide)] tp_with_guides <- tparams[,rep(guidesim_index, each = Nguide)] Xg <- rprocess(object, x0=x_with_guides, t0=times[nt+1], times=times[(nt+2):(nt+1+lookahead_steps)], @@ -588,11 +566,11 @@ igirf.bootgirf <- function (object, params, Ninter, lookahead, Nguide, params=tparams, times = times[(nt+2):(nt+1+lookahead_steps)], method = 'adams') - resids <- Xg - Xskel[,rep(1:Np[1], each=Nguide),,drop=FALSE] # residuals + resids <- Xg - Xskel[,rep(1:Np, each=Nguide),,drop=FALSE] # residuals rm(Xg, Xskel, x_with_guides) for (s in 1:Ninter){ tparams <- partrans(object,params,dir="fromEst",.gnsi=gnsi) - tp_with_guides <- tparams[,rep(1:Np[1], rep(Nguide, Np[1]))] + tp_with_guides <- tparams[,rep(1:Np, rep(Nguide, Np))] X <- rprocess(object,x0=x, t0 = tt[s], times= tt[s+1], params=tparams,.gnsi=gnsi) if(s>1 && length(znames)>0){ @@ -617,18 +595,26 @@ igirf.bootgirf <- function (object, params, Ninter, lookahead, Nguide, } else { skel <- X } - log_guide_fun <- vector(mode = "numeric", length = Np[1]) + 1 + log_guide_fun <- vector(mode = "numeric", length = Np) + 1 for(l in 1:lookahead_steps){ if(nt+1+l-lookahead <= 0) discount_denom_init <- object@t0 else discount_denom_init <- times[nt+1+l - lookahead] - discount_factor <- 1 - (times[nt+1+l] - tt[s+1])/(times[nt+1+l] - discount_denom_init)/ifelse(lookahead==1,2,1) ## to ensure that the discount factor does not become too small for L=1 and small s (which can lead to very uninformative guide function), increase the discount factor to at least 1/2 when L=1. - - # construct pseudo-simulations by adding simulated noise terms (residuals) to the skeletons - pseudosims <- skel[,rep(1:Np[1], each=Nguide),l,drop=FALSE] + - resids[,rep(guidesim_index-1, each=Nguide)*Nguide+rep(1:Nguide, Np[1]),l,drop=FALSE] - - resids[,rep(guidesim_index-1, each=Nguide)*Nguide+rep(1:Nguide, Np[1]),1,drop=FALSE] + - resids[,rep(guidesim_index-1, each=Nguide)*Nguide+rep(1:Nguide, Np[1]),1,drop=FALSE] * sqrt((times[nt+2]-tt[s+1])/(times[nt+2]-times[nt+1])) + discount_factor <- 1 - (times[nt+1+l] - tt[s+1])/(times[nt+1+l] - discount_denom_init)/ifelse(lookahead==1,2,1) + ## to ensure that the discount factor does not become too small for L=1 + ## and small s (which can lead to very uninformative guide function), increase + ## the discount factor to at least 1/2 when L=1. + + ## construct pseudo-simulations by adding simulated noise terms (residuals) + ## to the skeletons + pseudosims <- skel[,rep(1:Np, each=Nguide),l,drop=FALSE] + + resids[,rep(guidesim_index-1, each=Nguide)* + Nguide+rep(1:Nguide, Np),l,drop=FALSE] - + resids[,rep(guidesim_index-1, each=Nguide)* + Nguide+rep(1:Nguide, Np),1,drop=FALSE] + + resids[,rep(guidesim_index-1, each=Nguide)* + Nguide+rep(1:Nguide, Np),1,drop=FALSE] * + sqrt((times[nt+2]-tt[s+1])/(times[nt+2]-times[nt+1])) rm(skel) log_dmeas_weights <- tryCatch( (vec_dmeasure( @@ -658,9 +644,14 @@ igirf.bootgirf <- function (object, params, Ninter, lookahead, Nguide, ## params=tp_with_guides[,na_ix_col] ## ) ## } - - ldw <- array(log_dmeas_weights, c(nunits,Nguide,Np[1])) # log_dmeas_weights is an array with dim U*(Np*Nguide)*1. Reorder it as U*Nguide*Np - log_fcst_lik <- colSums(log(apply(exp(ldw),c(1,3),sum)/Nguide)) # average dmeas (natural scale) over Nguide sims, then take log, and then sum over 1:U (for each particle) + + ## log_dmeas_weights is an array with dim U*(Np*Nguide)*1. + ## Reorder it as U*Nguide*Np + ldw <- array(log_dmeas_weights, c(nunits,Nguide,Np)) + + ## average dmeas (natural scale) over Nguide sims, then take log, + ## and then sum over 1:U (for each particle) + log_fcst_lik <- colSums(log(apply(exp(ldw),c(1,3),sum)/Nguide)) log_resamp_weights <- log_fcst_lik*discount_factor log_guide_fun = log_guide_fun + log_resamp_weights } @@ -707,19 +698,20 @@ igirf.bootgirf <- function (object, params, Ninter, lookahead, Nguide, .Call(girf_computations, x=X, params=params, - Np=Np[nt+1], + Np=Np, trackancestry=TRUE, doparRS=TRUE, weights=weights, lgps=log_guide_fun, - fsv=array(0,dim=c(length(unit_names(object)), lookahead_steps, Np[1])), # bootgirf2 doesn't use fsv, set to an arbitrary val. + ## bootgirf2 doesn't use fsv, set to an arbitrary value + fsv=array(0,dim=c(length(unit_names(object)), lookahead_steps, Np)), tol=tol ), error = function (e) { - stop(ep,conditionMessage(e),call.=FALSE) # nocov + stop(ep,conditionMessage(e),call.=FALSE) } ) - guidesim_index = guidesim_index[xx$ancestry] # update guidesim index + guidesim_index = guidesim_index[xx$ancestry] cond.loglik[nt+1, s] <- xx$loglik + max_log_weights x <- xx$states log_filter_guide_fun <- xx$logfilterguides @@ -739,7 +731,7 @@ igirf.bootgirf <- function (object, params, Ninter, lookahead, Nguide, Nguide=Nguide, lookahead=lookahead, paramMatrix=params, - Np=Np[1], + Np=Np, tol=tol, loglik=sum(cond.loglik) ) @@ -760,7 +752,7 @@ igirf.bootgirf <- function (object, params, Ninter, lookahead, Nguide, igirf_rw_sd_times <- function(times, Ninter){ rw_sd_times <- c() for(i in seq(length(times)-1)) rw_sd_times <- c( - rw_sd_times, seq(from=times[i],to=times[i+1],length.out=Ninter+1)) + rw_sd_times, seq(from=times[i],to=times[i+1],length.out=Ninter+1)) return(rw_sd_times) } diff --git a/tests/bm.R b/tests/bm.R index b246350..c94c922 100644 --- a/tests/bm.R +++ b/tests/bm.R @@ -1,4 +1,3 @@ -png(filename="bm-%02d.png",res=100) library(spatPomp) i <- 1 U <- switch(i,2,10) @@ -38,7 +37,7 @@ set.seed(2) ## foreach confused about whether it is running in parallel or not. b_model <- bm(U=U,N=N) -b_model_no_rproc <- spatPomp(b_model,rprocess=NULL) +b_model_no_rprocess <- spatPomp(b_model,rprocess=NULL) b_model_no_eunit_measure <- spatPomp(b_model,eunit_measure=NULL) b_model_no_vunit_measure <- spatPomp(b_model,vunit_measure=NULL) b_model_t0_equal_t1 <- spatPomp(b_model,t0=1) @@ -135,7 +134,7 @@ try(bpfilter(b_bpfilter,params=test_params_matrix)) ## test error messages try(enkf()) try(enkf("JUNK")) -try(enkf(b_model_no_rproc)) +try(enkf(b_model_no_rprocess)) try(enkf(b_model_no_eunit_measure)) try(enkf(b_model_no_vunit_measure)) @@ -172,11 +171,11 @@ try(girf(b_girf_boot,Np=c(Inf))) try(girf(b_girf_boot,Np=seq(from=10,length=N+1,by=2))) try(girf(b_model_no_eunit_measure,kind='moment')) try(girf(b_model_no_vunit_measure,kind='moment')) -try(girf(b_model_no_rproc,kind='moment')) +try(girf(b_model_no_rprocess,kind='moment')) try(girf(b_model,kind='moment')) try(girf(b_model,kind='moment',Np=5)) try(girf(b_model,kind='moment',Np=5,Nguide=3,tol=1:1000)) -try(girf(b_model_no_rproc,kind='boot')) +try(girf(b_model_no_rprocess,kind='boot')) try(girf(b_model,kind='boot')) try(girf(b_model,kind='boot',Np=5)) try(girf(b_model,kind='boot',Np=5,Nguide=3,tol=1:1000)) @@ -268,18 +267,47 @@ print("The following deliver an error message, to test it") try(igirf()) try(igirf(data="JUNK")) try(igirf(b_igirf_boot_geom,Np=c(Inf))) -try(igirf(b_igirf_boot_geom,Np=3)) +igirf(b_igirf_boot_geom,Np=3) +try(igirf(b_igirf_boot_geom,Np=NULL)) try(igirf(b_model_no_eunit_measure,kind='moment', Ngirf = 2, Nguide=2, rw.sd = b_rw.sd, cooling.type = "hyperbolic", cooling.fraction.50 = 0.5, Np=floor(Np/2), Ninter = 2)) +try(igirf(b_model_no_vunit_measure,kind='moment', Ngirf = 2, Nguide=2, + rw.sd = b_rw.sd, cooling.type = "hyperbolic", cooling.fraction.50 = 0.5, + Np=floor(Np/2), Ninter = 2)) +try(igirf(b_model_no_rprocess,kind='moment', Ngirf = 2, Nguide=2, + rw.sd = b_rw.sd, cooling.type = "hyperbolic", cooling.fraction.50 = 0.5, + Np=floor(Np/2), Ninter = 2)) + +try(igirf(b_model)) +try(igirf(b_model,Ngirf=2)) +try(igirf(b_model,Ngirf=2,rw.sd=b_rw.sd)) +try(igirf(b_model,Ngirf=2,rw.sd=b_rw.sd,cooling.fraction.50=0.5)) +try(igirf(b_model,Ngirf=2,rw.sd=b_rw.sd,cooling.fraction.50=0.5, + cooling.type="geometric")) +try(igirf(b_model,Ngirf=2,rw.sd=b_rw.sd,cooling.fraction.50=0.5, + cooling.type="geometric",Np=4)) +try(igirf(b_model,Ngirf=2,rw.sd=b_rw.sd,cooling.fraction.50=0.5, + cooling.type="geometric",Np="JUNK",Nguide=4)) + +igirf(b_igirf_boot_geom,Np=3, + .paramMatrix=cbind(coef(b_model),coef(b_model),coef(b_model))) +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_model,kind='moment', Ngirf = 2, Nguide=2, + rw.sd = b_rw.sd, cooling.type = "hyperbolic", cooling.fraction.50 = 0.5, + Np=floor(Np/2), Ninter = 2)) # try(igirf(b_model_no_vunit_measure,kind='moment')) -# try(igirf(b_model_no_rproc,kind='moment')) +# try(igirf(b_model_no_rprocess,kind='moment')) # try(igirf(b_model,kind='moment')) # try(igirf(b_model,kind='moment',Np=5)) # try(igirf(b_model,kind='moment',Np=5,Nguide=3,tol=1:1000)) -# try(igirf(b_model_no_rproc,kind='boot')) +# try(igirf(b_model_no_rprocess,kind='boot')) # try(igirf(b_model,kind='boot')) # try(igirf(b_model,kind='boot',Np=5)) # try(igirf(b_model,kind='boot',Np=5,Nguide=3,tol=1:1000)) @@ -289,6 +317,8 @@ try(igirf(b_model_no_eunit_measure,kind='moment', Ngirf = 2, Nguide=2, ## ienkf on bm, with geometric and hyperbolic cooling ## +set.seed(55) + b_ienkf_geom <- ienkf(b_model, Nenkf=2, Np = Np, @@ -372,11 +402,13 @@ b_sim2 <- simulate(b_model,nsim=2,format='data.frame',include.data=TRUE) head(b_sim2,10) b_sim3 <- simulate(b_model,nsim=2,format='spatPomps') +png(filename="bm-%02d.png",res=100) plot(b_model,type="l",log=FALSE) b_sim3v2 <- b_sim3[[1]] b_sim3v2@data <- exp(b_sim3v2@data) plot(b_sim3v2,type="l",log=TRUE) plot(b_sim3[[2]],type="h") +dev.off() ## -------------------------------------------- ## using bm to test spatPomp workhorse functions, extending pomp: @@ -445,8 +477,6 @@ munit_measure(b_u, x=b_s, vc=b_vc, Np=1, unit = 2, time=1, dunit_measure(b_u, y=b_y,x=b_s, unit=2, time=1, params=b_p) runit_measure(b_u, x=b_s, unit=2, time=1, params=b_p) -dev.off() - ## test spatPomp_Csnippet variable construction spatPomp_Csnippet("lik=u;",unit_statenames="A",unit_obsnames=c("B","C"), unit_covarnames="D", diff --git a/tests/bm.Rout.save b/tests/bm.Rout.save index 324c159..c787531 100644 --- a/tests/bm.Rout.save +++ b/tests/bm.Rout.save @@ -17,7 +17,6 @@ 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. -> png(filename="bm-%02d.png",res=100) > library(spatPomp) Loading required package: pomp > i <- 1 @@ -58,7 +57,7 @@ Loading required package: pomp > ## foreach confused about whether it is running in parallel or not. > > b_model <- bm(U=U,N=N) -> b_model_no_rproc <- spatPomp(b_model,rprocess=NULL) +> b_model_no_rprocess <- spatPomp(b_model,rprocess=NULL) > b_model_no_eunit_measure <- spatPomp(b_model,eunit_measure=NULL) > b_model_no_vunit_measure <- spatPomp(b_model,vunit_measure=NULL) > b_model_t0_equal_t1 <- spatPomp(b_model,t0=1) @@ -189,7 +188,7 @@ Error : ‘params’ must be a named numeric vector. Error : in ‘enkf’ : ‘data’ is a required argument. > try(enkf("JUNK")) Error : ‘enkf’ is undefined for ‘data’ of class ‘character’. -> try(enkf(b_model_no_rproc)) +> try(enkf(b_model_no_rprocess)) Error : in ‘enkf’ : ‘rprocess’ is a needed basic component. > try(enkf(b_model_no_eunit_measure)) Error : in ‘enkf’ : ‘eunit_measure’ is a needed basic component. @@ -240,7 +239,7 @@ Error : in ‘girf’: ‘Np’ must be a single positive integer Error : girf with kind = 'moment' requires eunit_measure > try(girf(b_model_no_vunit_measure,kind='moment')) Error : girf with kind = 'moment' requires vunit_measure -> try(girf(b_model_no_rproc,kind='moment')) +> try(girf(b_model_no_rprocess,kind='moment')) Error : in ‘girf’: ‘rprocess’, ‘dmeasure’ are needed basic components. > try(girf(b_model,kind='moment')) Error : in ‘girf’: ‘Np’ must be specified. @@ -248,7 +247,7 @@ Error : in ‘girf’: ‘Np’ must be specified. Error : in ‘girf’: argument "Nguide" is missing, with no default > try(girf(b_model,kind='moment',Np=5,Nguide=3,tol=1:1000)) Error : in ‘girf’: ‘tol’ should be a small positive number. -> try(girf(b_model_no_rproc,kind='boot')) +> try(girf(b_model_no_rprocess,kind='boot')) Error : in ‘girf’: ‘rprocess’, ‘dmeasure’ are needed basic components. > try(girf(b_model,kind='boot')) Error : in ‘girf’: ‘Np’ must be specified. @@ -369,12 +368,15 @@ igirf iteration 2 of 2 completed with likelihood -11.11327 > try(igirf()) Error : igirf: data is a required argument. > try(igirf(data="JUNK")) -Error : igirf is undefined for ‘JUNK’of class ‘character’. +Error : igirf is undefined for ‘data’of class ‘character’. > try(igirf(b_igirf_boot_geom,Np=c(Inf))) -Error : in ‘igirf’ : Error: ‘Np’ must be a positive integer. +Error : in ‘igirf’ : Error: ‘Np’ must be a single positive integer -> try(igirf(b_igirf_boot_geom,Np=3)) +> igirf(b_igirf_boot_geom,Np=3) +> try(igirf(b_igirf_boot_geom,Np=NULL)) +Error : in ‘igirf’ : Error: ‘Np’ must be specified. + > try(igirf(b_model_no_eunit_measure,kind='moment', Ngirf = 2, Nguide=2, + rw.sd = b_rw.sd, cooling.type = "hyperbolic", cooling.fraction.50 = 0.5, + Np=floor(Np/2), Ninter = 2)) @@ -382,14 +384,64 @@ Error : in ‘igirf’ : Error in if (max_log_weights > -Inf) {: missing value w In addition: Warning message: 'eunit_measure' unspecified. +> try(igirf(b_model_no_vunit_measure,kind='moment', Ngirf = 2, Nguide=2, ++ rw.sd = b_rw.sd, cooling.type = "hyperbolic", cooling.fraction.50 = 0.5, ++ Np=floor(Np/2), Ninter = 2)) +Error : in ‘igirf’ : Error in if (max_log_weights > -Inf) {: missing value where TRUE/FALSE needed + +In addition: Warning message: +'vunit_measure' unspecified. +> try(igirf(b_model_no_rprocess,kind='moment', Ngirf = 2, Nguide=2, ++ rw.sd = b_rw.sd, cooling.type = "hyperbolic", cooling.fraction.50 = 0.5, ++ Np=floor(Np/2), Ninter = 2)) +Error : in ‘igirf’ : Error: ‘rprocess’‘dmeasure’ are needed basic components. + +> +> try(igirf(b_model)) +Error : in ‘igirf’ : ‘Ngirf’ is a required argument +> try(igirf(b_model,Ngirf=2)) +Error : in ‘igirf’ : ‘rw.sd’ must be specified. +> try(igirf(b_model,Ngirf=2,rw.sd=b_rw.sd)) +Error : in ‘igirf’ : ‘cooling.fraction.50’ is a required argument. +> try(igirf(b_model,Ngirf=2,rw.sd=b_rw.sd,cooling.fraction.50=0.5)) +Error : in ‘igirf’ : ‘cooling.type’ is a required argument. +> try(igirf(b_model,Ngirf=2,rw.sd=b_rw.sd,cooling.fraction.50=0.5, ++ cooling.type="geometric")) +Error : in ‘igirf’ : ‘Np’ must be specified. +> try(igirf(b_model,Ngirf=2,rw.sd=b_rw.sd,cooling.fraction.50=0.5, ++ cooling.type="geometric",Np=4)) +Error : in ‘igirf’ : ‘Nguide’ must be specified. +> try(igirf(b_model,Ngirf=2,rw.sd=b_rw.sd,cooling.fraction.50=0.5, ++ cooling.type="geometric",Np="JUNK",Nguide=4)) +Error : in ‘igirf’ : Error: ‘Np’ must be a single positive integer + +> +> igirf(b_igirf_boot_geom,Np=3, ++ .paramMatrix=cbind(coef(b_model),coef(b_model),coef(b_model))) + +> try(igirf(b_igirf_boot_geom,Np=function(x) 4)) +Error : in ‘igirf’ : Error: ‘Np’ must be a single positive integer + +> try(igirf(b_igirf_boot_geom,Np=function(x) "JUNK")) +Error : in ‘igirf’ : Error: ‘Np’ must be a single positive integer + +> try(igirf(b_igirf_boot_geom,Np=5:15)) +Error : in ‘igirf’ : Error: ‘Np’ must be a single positive integer + +> +> +> try(igirf(b_model,kind='moment', Ngirf = 2, Nguide=2, ++ rw.sd = b_rw.sd, cooling.type = "hyperbolic", cooling.fraction.50 = 0.5, ++ Np=floor(Np/2), Ninter = 2)) + > > > # try(igirf(b_model_no_vunit_measure,kind='moment')) -> # try(igirf(b_model_no_rproc,kind='moment')) +> # try(igirf(b_model_no_rprocess,kind='moment')) > # try(igirf(b_model,kind='moment')) > # try(igirf(b_model,kind='moment',Np=5)) > # try(igirf(b_model,kind='moment',Np=5,Nguide=3,tol=1:1000)) -> # try(igirf(b_model_no_rproc,kind='boot')) +> # try(igirf(b_model_no_rprocess,kind='boot')) > # try(igirf(b_model,kind='boot')) > # try(igirf(b_model,kind='boot',Np=5)) > # try(igirf(b_model,kind='boot',Np=5,Nguide=3,tol=1:1000)) @@ -399,6 +451,8 @@ In addition: Warning message: > ## ienkf on bm, with geometric and hyperbolic cooling > ## > +> set.seed(55) +> > b_ienkf_geom <- ienkf(b_model, + Nenkf=2, + Np = Np, @@ -408,7 +462,7 @@ In addition: Warning message: + verbose=FALSE + ) > paste("bm ienkf loglik, geometric cooling, verbose=F: ",round(logLik(b_ienkf_geom),10)) -[1] "bm ienkf loglik, geometric cooling, verbose=F: -11.6995986214" +[1] "bm ienkf loglik, geometric cooling, verbose=F: -12.7660513791" > > b_ienkf_hyp <- ienkf(b_model, + Nenkf=2, @@ -422,11 +476,12 @@ ienkf iteration 1 of 2 completed ienkf iteration 2 of 2 completed > > paste("bm ienkf loglik, hyperbolic cooling, verbose=T: ",round(logLik(b_ienkf_hyp),10)) -[1] "bm ienkf loglik, hyperbolic cooling, verbose=T: -9.7365825263" +[1] "bm ienkf loglik, hyperbolic cooling, verbose=T: -10.2194988409" > > ## test error messages for ienkf > try(ienkf(b_model_no_rproc)) -Error : in ‘ienkf’: ‘rprocess’, ‘eunit_measure’, ‘vunit_measure’ are needed basic components. +Error in h(simpleError(msg, call)) : + error in evaluating the argument 'data' in selecting a method for function 'ienkf': object 'b_model_no_rproc' not found > try(ienkf(b_model, Nenkf="JUNK")) Error : in ‘ienkf’: ‘Nenkf’ must be a positive integer. > ienkf(b_model,Nenkf=2,Np = 3,rw.sd=b_rw.sd,cooling.type="geometric", @@ -532,11 +587,15 @@ iubf iteration 2 of 2 completed with log-likelihood -10.16506 10 2 data U2 -0.759 -0.849 > b_sim3 <- simulate(b_model,nsim=2,format='spatPomps') > +> png(filename="bm-%02d.png",res=100) > plot(b_model,type="l",log=FALSE) > b_sim3v2 <- b_sim3[[1]] > b_sim3v2@data <- exp(b_sim3v2@data) > plot(b_sim3v2,type="l",log=TRUE) > plot(b_sim3[[2]],type="h") +> dev.off() +null device + 1 > > ## -------------------------------------------- > ## using bm to test spatPomp workhorse functions, extending pomp: @@ -680,10 +739,6 @@ variable [,1] Y1 1 > -> dev.off() -null device - 1 -> > ## test spatPomp_Csnippet variable construction > spatPomp_Csnippet("lik=u;",unit_statenames="A",unit_obsnames=c("B","C"), + unit_covarnames="D",