more tests

DESCRIPTION

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

R/bm2.R

@@ -191,9 +191,8 @@ bm2 <- function(U=5,N=100,delta_t=0.1,
 
 log_unit_names <- c("sigma", "tau")
 logit_unit_names <- c("rho")
-log_names <- unlist(lapply(log_unit_names,
-  function(x,y,U){if(x%in%y)paste0(x,"1") else paste0(x,1:U)},
-  y=shared_names,U=U))
+log_name_func <- function(x,y,U) if(x%in%y)paste0(x,"1") else paste0(x,1:U)
+log_names <- unlist(lapply(log_unit_names,log_name_func,y=shared_names,U=U))
 logit_names <- unlist(lapply(logit_unit_names,
   function(x,y,U){if(x%in%y)paste0(x,"1") else paste0(x,1:U)},
   y=shared_names,U=U))

R/ienkf.R

@@ -293,19 +293,17 @@ ienkf.filter <- function (object, params, Np, enkfiter, rw.sd, cooling.fn,
     R <- diag(rowMeans(meas_var))
     sqrtR <- tryCatch(
       t(chol(R)),                     # t(sqrtR)%*%sqrtR == R
-      error = function (e) {
-        pStop_("degenerate ",sQuote("R"), "at time ", sQuote(nt), ": ",conditionMessage(e))
-      }
+      error = function (e) pStop_("degenerate ",sQuote("R"), "at time ", sQuote(nt), ": ",conditionMessage(e))
     )
 
-                                        # expand the state space
+    # expand the state space
     XT <- rbind(X[,,1],params)
     pm <- rowMeans(XT) # prediction mean
 
-                                        # forecast mean
+    # forecast mean
     ym <- rowMeans(Y)
 
-                                        # center prediction and forecast ensembles
+    # center prediction and forecast ensembles
     XT <- XT-pm
     Y <- Y-ym
 
@@ -320,8 +318,8 @@ ienkf.filter <- function (object, params, Np, enkfiter, rw.sd, cooling.fn,
     XT <- XT+pm+crossprod(Kt,resid-Y+Ek)
     params <- XT[pnames,,drop = FALSE]
     X <- XT[xnames,,drop = FALSE]
-    loglik[nt] <- sum(dnorm(x=crossprod(svdS$u,resid),mean=0,sd=sqrt(svdS$d),log=TRUE))
-                                        # print(rowMeans(partrans(object,params,dir="fromEst",.gnsi=gnsi)))
+    loglik[nt] <- sum(dnorm(x=crossprod(svdS$u,resid),
+      mean=0,sd=sqrt(svdS$d),log=TRUE))
 
     ## compute mean at last timestep
     if (nt == ntimes) {

R/iubf.R

@@ -307,10 +307,12 @@ iubf_internal <- function (object, Nrep_per_param, Nparam, nbhd, Nubf, prop, rw.
   if (is.null(.paramMatrix)) {
     rep_param_init <- coef(object)
     start <- coef(object)
-  } else {
-    rep_param_init <- .paramMatrix
-    start <- apply(.paramMatrix,1L,mean)
   }
+  ## continuation using .paramMatrix not yet implemented
+  ## else {
+  ##   rep_param_init <- .paramMatrix
+  ##   start <- apply(.paramMatrix,1L,mean)
+  ## }
 
   rw.sd <- perturbn.kernel.sd(rw.sd,time=time(object),paramnames=names(start))
 

R/undefined.R

@@ -51,15 +51,6 @@ setMethod(
     object@mode == pompfunmode$undef
   }
 )
-##' @rdname undefined
-setMethod(
-  "undefined",
-  signature=signature(object="partransPlugin"),
-  definition=function (object, ...) {
-    undefined(object@to) || undefined(object@from)
-  }
-)
-
 ##' @rdname undefined
 setMethod(
   "undefined",
@@ -68,20 +59,3 @@ setMethod(
     undefined(object@step.fn) && undefined(object@rate.fn)
   }
 )
-##' @rdname undefined
-setMethod(
-  "undefined",
-  signature=signature(object="covartable"),
-  definition=function (object) {
-    nrow(object@table) == 0L
-  }
-)
-
-##' @rdname undefined
-setMethod(
-  "undefined",
-  signature=signature(object="skelPlugin"),
-  definition=function (object, ...) {
-    object@type==skeletontype$undef || undefined(object@skel.fn)
-  }
-)

src/abf.c

@@ -8,16 +8,17 @@ SEXP abf_computations (SEXP x, SEXP params, SEXP Np,
 			   SEXP resamp_weights)
 {
   int nprotect = 0;
-  SEXP anc = R_NilValue;
   SEXP newstates = R_NilValue;
   SEXP retval, retvalnames;
   const char *dimnm[2] = {"variable","rep"};
   double *xx = 0;
-  int *xanc = 0;
   SEXP dimX, dimP, Xnames, Pnames;
   int *dim, np;
   int nvars, nreps;
-  int do_ta;
+  // ancestor tracking currently not implemented
+  // int do_ta;
+  // SEXP anc = R_NilValue;
+  // int *xanc = 0;
   int j, k, l;
 
   PROTECT(dimX = GET_DIM(x)); nprotect++;
@@ -34,12 +35,12 @@ SEXP abf_computations (SEXP x, SEXP params, SEXP Np,
 
   np = 1;       // number of particles to resample for abf is just one
 
-  do_ta = *(LOGICAL(AS_LOGICAL(trackancestry))); // track ancestry?
-
-  if (do_ta) {
-    PROTECT(anc = NEW_INTEGER(np)); nprotect++;
-    xanc = INTEGER(anc);
-  }
+  // ancestry tracking not implemented
+  // do_ta = *(LOGICAL(AS_LOGICAL(trackancestry))); // track ancestry?
+  // if (do_ta) {
+  //   PROTECT(anc = NEW_INTEGER(np)); nprotect++;
+  //   xanc = INTEGER(anc);
+  // }
   GetRNGstate();
 
   int xdim[2];
@@ -64,7 +65,7 @@ SEXP abf_computations (SEXP x, SEXP params, SEXP Np,
 
   for (k = 0; k < copies; k++) { // copy the particles
     for (j = 0, xx = ss+nvars*sample[k]; j < nvars; j++, st++, xx++) *st = *xx;
-    if (do_ta) xanc[k] = sample[k]+1;
+  //  if (do_ta) xanc[k] = sample[k]+1;
   }
 
   PutRNGstate();
@@ -78,9 +79,9 @@ SEXP abf_computations (SEXP x, SEXP params, SEXP Np,
   SET_ELEMENT(retval,0,newstates);
   SET_ELEMENT(retval,1,params);
 
-  if (do_ta) {
-    SET_ELEMENT(retval,2,anc);
-  }
+  // if (do_ta) {
+  //   SET_ELEMENT(retval,2,anc);
+  // }
 
   UNPROTECT(nprotect);
   return(retval);

src/resample.c

@@ -3,6 +3,25 @@
 #include "spatPomp_defines.h"
 #include <Rdefines.h>
 
+SEXP spatPomp_systematic_resampling (SEXP weights, SEXP np);
+void nosort_resamp (int nw, double *w, int np, int *p, int offset);
+
+SEXP spatPomp_systematic_resampling (SEXP weights, SEXP np)  
+{
+  int m, n;
+  SEXP perm;
+
+  m = *(INTEGER(AS_INTEGER(np)));
+  n = LENGTH(weights);
+  PROTECT(perm = NEW_INTEGER(m));
+  PROTECT(weights = AS_NUMERIC(weights));
+  GetRNGstate();
+  nosort_resamp(n,REAL(weights),m,INTEGER(perm),1);
+  PutRNGstate();
+  UNPROTECT(2);
+  return(perm);
+}
+
 void nosort_resamp (int nw, double *w, int np, int *p, int offset)
 {
   int i, j;
@@ -23,7 +42,7 @@ void nosort_resamp (int nw, double *w, int np, int *p, int offset)
     while ((u > w[i]) && (i < nw-1)) i++;
     p[j] = i;
   }
-  if (offset)			// add offset if needed
+  if (offset)     // add offset if needed
     for (j = 0; j < np; j++) p[j] += offset;
 
 }

tests/bm.R

@@ -158,7 +158,7 @@ b_bpfilter_repeat <- bpfilter(b_bpfilter)
 paste("check bpfilter on bpfilterd_spatPomp: ",
   logLik(b_bpfilter)==logLik(b_bpfilter_repeat))
 
-bpfilter(b_model_t0_equal_t1,Np = Np, block_size = 1)
+bpfilter(b_model_t0_equal_t1,Np = Np, block_size = 1,filter_traj=TRUE)
 
 set.seed(5)
 b_bpfilter_filter_traj <- bpfilter(b_bpfilter,filter_traj=TRUE)
@@ -442,7 +442,7 @@ b_ienkf_hyp <- ienkf(b_model,
 paste("bm ienkf loglik, hyperbolic cooling, verbose=T: ",round(logLik(b_ienkf_hyp),10))
 
 ## test error messages for ienkf
-try(ienkf(b_model_no_rproc))
+try(ienkf(b_model_no_rprocess))
 try(ienkf(b_model, Nenkf="JUNK"))
 ienkf(b_model,Nenkf=2,Np = 3,rw.sd=b_rw.sd,cooling.type="geometric",
   cooling.fraction.50 = 0.5,
@@ -611,22 +611,25 @@ try(dunit_measure(b_model, y=b_y,
   x=b_s, unit=1, time=1:10, params=b_p))  
 b_s2 <- 1:6
 dim(b_s2) <- c(2,3,1)
+b_s3 <- 1:6
+dim(b_s3) <- c(1,3,2)
 dimnames(b_s2) <- list(variable=dimnames(states(b_model))[[1]], rep=NULL)
 b_p2 <- c(coef(b_model),coef(b_model))
 dim(b_p2) <- c(length(coef(b_model)),2)
 dimnames(b_p2) <- list(param=names(coef(b_model)))
 try(dunit_measure(b_model, y=b_y, x=b_s2, unit=1, time=1, params=b_p2))
+try(dunit_measure(b_model, y=b_y, x=b_s3, unit=1, time=1, params=b_p2))
 
 ## trigger error messages in runit_measure.c
 runit_measure(b_model_no_runit_measure, x=b_s, unit=2, time=1, params=b_p)
 try(runit_measure(b_model_no_runit_measure, x=b_s, unit=2, time=numeric(0), params=b_p))
 try(runit_measure(b_model_no_runit_measure, x=b_s, unit=2, time=1:10, params=b_p))
-try(runit_measure(b_model_no_runit_measure, x=b_s2, unit=2, time=1:10, params=b_p2))
+try(runit_measure(b_model_no_runit_measure, x=b_s2, unit=2, time=1,params=b_p2))
 
 ## trigger error messages in vunit_measure.c
 vunit_measure(b_model_no_vunit_measure, x=b_s, unit=2, time=1, params=b_p)
 try(vunit_measure(b_model, x=b_s, unit=2, time=1:10, params=b_p))
-try(vunit_measure(b_model, x=b_s2, unit=2, time=1, params=b_p2))
+try(vunit_measure(b_model, x=b_s3, unit=2, time=1, params=b_p2))
 
 ## trigger error messages in eunit_measure.c
 eunit_measure(b_model_no_eunit_measure, x=b_s, unit=2, time=1, params=b_p)
@@ -727,3 +730,14 @@ try(spatPomp(data=model_unit_covar,covar=b_unit_covar,
   shared_covarnames="Z"))
 
 
+## --------------------------------------------
+## test utility functions
+## ____________________________________________
+
+spatPomp:::undefined()
+spatPomp:::undefined(NULL)
+spatPomp:::undefined("JUNK")
+
+try(.Call("spatPomp_systematic_resampling",c(-0.1,-0.2),2))
+
+

tests/bm.Rout.save

@@ -216,7 +216,7 @@ Error in spatPomp:::abfir_internal(object = object, Np = Np, nbhd = nbhd,  :
 +   logLik(b_bpfilter)==logLik(b_bpfilter_repeat))
 [1] "check bpfilter on bpfilterd_spatPomp:  TRUE"
 > 
-> bpfilter(b_model_t0_equal_t1,Np = Np, block_size = 1)
+> bpfilter(b_model_t0_equal_t1,Np = Np, block_size = 1,filter_traj=TRUE)
 <object of class ‘bpfilterd_spatPomp’>
 > 
 > set.seed(5)
@@ -635,9 +635,8 @@ ienkf iteration 2 of 2 completed
 [1] "bm ienkf loglik, hyperbolic cooling, verbose=T:  -12.8788856955"
 > 
 > ## test error messages for ienkf
-> try(ienkf(b_model_no_rproc))
-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_no_rprocess))
+Error : in ‘ienkf’: ‘rprocess’, ‘eunit_measure’, ‘vunit_measure’ are needed basic components.
 > 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",
@@ -958,12 +957,16 @@ Warning message:
 Error : length of 'times' and 2nd dimension of 'y' do not agree.
 > b_s2 <- 1:6
 > dim(b_s2) <- c(2,3,1)
+> b_s3 <- 1:6
+> dim(b_s3) <- c(1,3,2)
 > dimnames(b_s2) <- list(variable=dimnames(states(b_model))[[1]], rep=NULL)
 > b_p2 <- c(coef(b_model),coef(b_model))
 > dim(b_p2) <- c(length(coef(b_model)),2)
 > dimnames(b_p2) <- list(param=names(coef(b_model)))
 > try(dunit_measure(b_model, y=b_y, x=b_s2, unit=1, time=1, params=b_p2))
 Error : larger number of replicates is not a multiple of smaller.
+> try(dunit_measure(b_model, y=b_y, x=b_s3, unit=1, time=1, params=b_p2))
+Error : length of 'times' and 3rd dimension of 'x' do not agree.
 > 
 > ## trigger error messages in runit_measure.c
 > runit_measure(b_model_no_runit_measure, x=b_s, unit=2, time=1, params=b_p)
@@ -979,8 +982,8 @@ Warning message:
 Error : length('times') = 0, no work to do.
 > try(runit_measure(b_model_no_runit_measure, x=b_s, unit=2, time=1:10, params=b_p))
 Error : length of 'times' and 3rd dimension of 'x' do not agree.
-> try(runit_measure(b_model_no_runit_measure, x=b_s2, unit=2, time=1:10, params=b_p2))
-Error : length of 'times' and 3rd dimension of 'x' do not agree.
+> try(runit_measure(b_model_no_runit_measure, x=b_s2, unit=2, time=1,params=b_p2))
+Error : larger number of replicates is not a multiple of smaller.
 > 
 > ## trigger error messages in vunit_measure.c
 > vunit_measure(b_model_no_vunit_measure, x=b_s, unit=2, time=1, params=b_p)
@@ -994,8 +997,8 @@ Warning message:
 'vunit_measure' unspecified. 
 > try(vunit_measure(b_model, x=b_s, unit=2, time=1:10, params=b_p))
 Error : length of 'times' and 3rd dimension of 'x' do not agree.
-> try(vunit_measure(b_model, x=b_s2, unit=2, time=1, params=b_p2))
-Error : larger number of replicates is not a multiple of smaller.
+> try(vunit_measure(b_model, x=b_s3, unit=2, time=1, params=b_p2))
+Error : length of 'times' and 3rd dimension of 'x' do not agree.
 > 
 > ## trigger error messages in eunit_measure.c
 > eunit_measure(b_model_no_eunit_measure, x=b_s, unit=2, time=1, params=b_p)
@@ -1144,4 +1147,19 @@ Error : in ‘spatPomp’: for unit-specific covariates, there should be a colum
 Error : in ‘spatPomp’: ‘shared_covarnames’ currently supported only when there are no unit-specific covariates
 > 
 > 
+> ## --------------------------------------------
+> ## test utility functions
+> ## ____________________________________________
+> 
+> spatPomp:::undefined()
+[1] NA
+> spatPomp:::undefined(NULL)
+[1] TRUE
+> spatPomp:::undefined("JUNK")
+[1] NA
+> 
+> try(.Call("spatPomp_systematic_resampling",c(-0.1,-0.2),2))
+Error : in 'systematic_resampling': non-positive sum of weights
+> 
+> 
 > 

tests/measles2.R

@@ -31,7 +31,8 @@ par1[paste0('iota',1:U)] <- 10
 # all shared parameters
 m2 <- measles2(U=U,N=N,expandedParNames=NULL,
   contractedParNames=c("R0", "c", "A", "muIR",
-    "muEI", "sigmaSE", "rho", "psi", "g", "S_0", "E_0", "I_0")
+    "muEI", "sigmaSE", "rho", "psi", "g", "S_0", "E_0", "I_0"),
+  simulated=TRUE   ## to test this flag for measles2()
 )
 
 set.seed(1)

tests/measles2.Rout.save

@@ -59,7 +59,8 @@ name     [,1] [,2]
 > # all shared parameters
 > m2 <- measles2(U=U,N=N,expandedParNames=NULL,
 +   contractedParNames=c("R0", "c", "A", "muIR",
-+     "muEI", "sigmaSE", "rho", "psi", "g", "S_0", "E_0", "I_0")
++     "muEI", "sigmaSE", "rho", "psi", "g", "S_0", "E_0", "I_0"),
++   simulated=TRUE   ## to test this flag for measles2()
 + )
 > 
 > set.seed(1)