experimental lookup table repair

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: pomp
 Type: Package
 Title: Statistical Inference for Partially Observed Markov Processes
-Version: 5.8.1.1
+Version: 5.8.2.0
 Date: 2024-05-13
 Authors@R: c(person(given=c("Aaron","A."),family="King",role=c("aut","cre"),email="kingaa@umich.edu",comment=c(ORCID="0000-0001-6159-3207")),
 	  person(given=c("Edward","L."),family="Ionides",role="aut",comment=c(ORCID="0000-0002-4190-0174"))	,

NAMESPACE

@@ -62,6 +62,7 @@ export(probe_sd)
 export(probe_var)
 export(profile_design)
 export(rbetabinom)
+export(repair_lookup_table)
 export(reulermultinom)
 export(rgammawn)
 export(ricker)

R/covariate_table.R

@@ -27,6 +27,9 @@
 ##' @param \dots numeric vectors or data frames containing time-varying covariates.
 ##' It must be possible to bind these into a data frame.
 ##'
+##' @return
+##' \code{covariate_table} returns a lookup table suitable for inclusion of covariates in a \sQuote{pomp} object.
+##' Specifically, this is an object of class \sQuote{covartable}.
 ##' @section Extrapolation:
 ##' If \code{t} is outside the range of the lookup table, the values will be extrapolated, and a warning will be issued.
 ##' The type of extrapolation performed will be constant or linear according to the \code{order} flag used when creating the table.

R/lookup.R

@@ -17,3 +17,21 @@ lookup <- function (table, t) {
   d <- .Call(P_lookup_in_table,table,t)
   data.frame(t=t,t(d))
 }
+
+##' @rdname covariate_table
+##' @return
+##' \code{repair_lookup_table} returns a lookup table with entries at the provided values of \code{t}.
+##' @param table a \sQuote{covartable} object created by a call to \code{\link{covariate_table}}
+##' @param t numeric vector;
+##' times at which interpolated values of the covariates in \code{table} are required.
+##' @details
+##' \code{repair_lookup_table} applies \code{\link{lookup}} at the provided values of \code{t} and returns the resulting lookup table.
+##' \strong{\code{repair_lookup_table} should be considered experimental: its interface may change without notice}.
+##' @export
+repair_lookup_table <- function (table, t) {
+  covariate_table(
+    lookup(table,t=t),
+    order=if (table@order==0L) "constant" else "linear",
+    times="t"
+  )
+}

inst/NEWS

@@ -1,5 +1,13 @@
 _N_e_w_s _f_o_r _p_a_c_k_a_g_e '_p_o_m_p'
 
+_C_h_a_n_g_e_s _i_n '_p_o_m_p' _v_e_r_s_i_o_n _5._8._2:
+
+        • A bug in covariate-table extrapolation for the case
+          ‘order="constant"’ has been fixed.
+
+        • A new experimental function ‘repair_lookup_table’ is provided
+          to help eliminate unnecessary warnings about extrapolation.
+
 _C_h_a_n_g_e_s _i_n '_p_o_m_p' _v_e_r_s_i_o_n _5._8._1:
 
         • Functions with names that include dots (‘.’), which have been

inst/NEWS.Rd

@@ -1,5 +1,11 @@
 \name{NEWS}
 \title{News for package `pomp'}
+\section{Changes in \pkg{pomp} version 5.8.2}{
+  \itemize{
+    \item A bug in covariate-table extrapolation for the case \code{order="constant"} has been fixed.
+    \item A new experimental function \code{repair_lookup_table} is provided to help eliminate unnecessary warnings about extrapolation.
+  }
+}
 \section{Changes in \pkg{pomp} version 5.8.1}{
   \itemize{
     \item Functions with names that include dots (\code{.}), which have been defunct since version 4.7, have been expunged.

src/lookup_table.c

@@ -53,25 +53,26 @@ SEXP lookup_in_table (SEXP covar, SEXP t) {
 void table_lookup (lookup_table_t *tab, double x, double *y)
 {
   int flag = 0;
-  int j, k;
+  int j, k, n;
   double e;
   if ((tab == 0) || (tab->length < 1) || (tab->width < 1)) return;
-  tab->index = findInterval(tab->x,tab->length,x,TRUE,TRUE,tab->index,&flag);
+  tab->index = findInterval(tab->x,tab->length,x,1,1,tab->index,&flag);
   // warn only if we are *outside* the interval
-  if ((x < tab->x[0]) || (x > tab->x[(tab->length)-1]))
+  if ((x < tab->x[0]) || (x > tab->x[tab->length-1]))
     warn("in 'table_lookup': extrapolating at %le.", x);
   switch (tab->order) {
   case 1: default: // linear interpolation
     e = (x - tab->x[tab->index-1]) / (tab->x[tab->index] - tab->x[tab->index-1]);
-    for (j = 0; j < tab->width; j++) {
-      k = j+(tab->width)*(tab->index);
-      y[j] = e*(tab->y[k])+(1-e)*(tab->y[k-tab->width]);
+    for (j = 0, k = tab->index*tab->width, n = k-tab->width; j < tab->width; j++, k++, n++) {
+      y[j] = e*(tab->y[k])+(1-e)*(tab->y[n]);
     }
     break;
   case 0: // piecewise constant
-    for (j = 0; j < tab->width; j++) {
-      k = j+(tab->width)*(tab->index);
-      y[j] = tab->y[k-tab->width];
+    if (flag < 0) n = 0;
+    else if (flag > 0) n = tab->index;
+    else n = tab->index-1;
+    for (j = 0, k = n*tab->width; j < tab->width; j++, k++) {
+      y[j] = tab->y[k];
     }
     break;
   }

tests/lookup.R

@@ -52,4 +52,25 @@ covariate_table(
     labeller=labeller(variable=c(f1="derivative",f0="function")))+
   theme_bw()
 
+covariate_table(
+  x=seq(1,10,by=1),
+  y=seq(2,20,by=2),
+  times="x"
+) -> tab
+lookup(tab,c(0,2.5,3.6,10,20))
+
+repair_lookup_table(tab,t=c(seq(0,10),20)) -> tab2
+lookup(tab2,c(0,2.5,3.6,10,20))
+
+covariate_table(
+  x=seq(1,10,by=1),
+  y=seq(1,10,by=1),
+  order="const",
+  times="x"
+) -> tab
+lookup(tab,c(0,2.5,3.6,10,20))
+
+repair_lookup_table(tab,t=c(seq(0,10,by=1),20)) -> tab2
+lookup(tab2,c(0,2.5,3.6,10,20))
+
 dev.off()

tests/lookup.Rout.save

@@ -127,6 +127,63 @@ Error : 'covariate_table' is undefined for 'times' of class 'logical'.
 +     labeller=labeller(variable=c(f1="derivative",f0="function")))+
 +   theme_bw()
 > 
+> covariate_table(
++   x=seq(1,10,by=1),
++   y=seq(2,20,by=2),
++   times="x"
++ ) -> tab
+> lookup(tab,c(0,2.5,3.6,10,20))
+     t    y
+1  0.0  0.0
+2  2.5  5.0
+3  3.6  7.2
+4 10.0 20.0
+5 20.0 40.0
+Warning messages:
+1: in 'table_lookup': extrapolating at 0.000000e+00. 
+2: in 'table_lookup': extrapolating at 2.000000e+01. 
+> 
+> repair_lookup_table(tab,t=c(seq(0,10),20)) -> tab2
+Warning messages:
+1: in 'table_lookup': extrapolating at 0.000000e+00. 
+2: in 'table_lookup': extrapolating at 2.000000e+01. 
+> lookup(tab2,c(0,2.5,3.6,10,20))
+     t    y
+1  0.0  0.0
+2  2.5  5.0
+3  3.6  7.2
+4 10.0 20.0
+5 20.0 40.0
+> 
+> covariate_table(
++   x=seq(1,10,by=1),
++   y=seq(1,10,by=1),
++   order="const",
++   times="x"
++ ) -> tab
+> lookup(tab,c(0,2.5,3.6,10,20))
+     t  y
+1  0.0  1
+2  2.5  2
+3  3.6  3
+4 10.0 10
+5 20.0 10
+Warning messages:
+1: in 'table_lookup': extrapolating at 0.000000e+00. 
+2: in 'table_lookup': extrapolating at 2.000000e+01. 
+> 
+> repair_lookup_table(tab,t=c(seq(0,10,by=1),20)) -> tab2
+Warning messages:
+1: in 'table_lookup': extrapolating at 0.000000e+00. 
+2: in 'table_lookup': extrapolating at 2.000000e+01. 
+> lookup(tab2,c(0,2.5,3.6,10,20))
+     t  y
+1  0.0  1
+2  2.5  2
+3  3.6  3
+4 10.0 10
+5 20.0 10
+> 
 > dev.off()
 null device 
           1