add option for evaluation at fixed heights, fix bug in poisson version

R/sedrate_gen_helpers.R

@@ -46,51 +46,77 @@ crppp = function(x_min, x_max, rate = 1){
   return(points)
 }
 
-sed_rate_from_matrix = function(height, sedrate, matrix, rate = 1, expand_domain = TRUE){
+sed_rate_from_matrix = function(height, sedrate, matrix, mode = "deterministic", rate = 1, expand_domain = TRUE, transform = identity){
   #' @export
   #' @title make sed rate gen from matrix
   #' 
   #' @param height vector of heights
   #' @param sedrate vector of sed. rates x values
-  #' @param matrix matrix of sed rate y values
+  #' @param matrix matrix of sed rate y values. Must have as many columns as length(height) and as many rows as length(sedrate).
+  #' @param mode character, "deterministic" or "poisson". Determines at which stratigraphic heights the sed rate is determined. If "deterministic" this will be the heights in `height`, if "poisson" the heights where the sed rate is determined follows a poisson point process with rate specified by `rate`
   #' @param rate numeric, rate of the Poisson point process determining frequency of sedimentation rate changes.
   #' @param expand_domain should sedimentation rates be defined below/above the highest/lowest height in the section? If TRUE, the sed rate values are the values at the closest interpolated point, if FALSE it will be NA
+  #' @param transform a function, the identity function by default. How should the values of the (pseudo)pdf defined by the entries of `matrix` be transformed? Using this function allows to (nonlinearly) rescale the values in matrix to put more emphasis on higher/lower values
   #' 
   #' @returns a function factory for usage with `sedrate_to_multiadm`
   #' 
   #' @seealso [sedrate_to_multiadm()] for estimating sedimentation rates based on the outputs, [get_data_from_eTimeOpt()] for extracting data from the `eTimeOpt` function of the astrochron package.
   #' 
   #' @description
-    #' at height `height[i]`, the sedimentation rate is specified by the pdf `approxfun(sedrate, matrix[i,]`)
-    #' 
+    #' Construct a sedimentation rate generator (function factory) from a matrix, e.g. one returned from `get_data_from_eTimeOpt`. This generator can be passed on to `sedrate_to_multiadm` to estimate age-depth models from it. 
+    #' If mode is "deterministic", the generator evaluates the sedimentation rates at heights specified by `height`, if the mode is "poisson" it is evaluated at heights that are determined based on a poisson point process. At these heights, the value of the sedimentation rate is determined based on the (pseudo) pdf that is determined by the matrix values.
+  if (!all(dim(matrix) == c(length(sedrate), length(height)))){
+    stop("dimension mismatch. \"matrix\" must have length(height) columns and length(sedrate) rows")
+  }
 
-  # x_min = -2
-  # x_max = 3
-  # height = seq(x_min, x_max, by = 0.25)
-  # sedrate = seq(0.1, 10, by = 0.1)
-  # matrix = matrix(data = runif(n = length(height) * length(sedrate)), nrow = length(height), ncol = length(sedrate))
   rule = 1
   if (expand_domain == TRUE){
     rule = 2
   }
-  f = function(){
-  x_max = max(height)
-  x_min = min(height)
-  interp_points = sort(c(x_min, x_max, crppp(x_min, x_max, rate)))
-  interp_heights = rep(NA, length(interp_points))
-  interp_vals = rep(NA, length(interp_points))
-  se = rep(NA, length(interp_points))
-  for (i in seq_along(interp_points)){
-    interp_index = which.min(abs(interp_points[i] - height))
-    sed_rate_vals = matrix[,interp_index]
-    sed_rate_val = rej_sampling(sedrate, sed_rate_vals)
-    interp_heights[i] = height[interp_index]
-    se[i] = sed_rate_val
-
+  if (mode == "poisson"){
+    f = function(){
+      x_max = max(height)
+      x_min = min(height)
+      interp_points = sort(c(x_min, x_max, crppp(x_min, x_max, rate)))
+      interp_heights = rep(NA, length(interp_points))
+      interp_vals = rep(NA, length(interp_points))
+      se = rep(NA, length(interp_points))
+      for (i in seq_along(interp_points)){
+        low_int = height[height <= interp_points[i]]   #heights below poi
+        high_int = height[height >= interp_points[i]]  # heights above poi
+        h_low = max(low_int)  # highest point below poi
+        h_high = min(high_int) # lowest point above poi
+        low_ind = which(height ==h_low)
+        high_ind = which(height == h_high)
+        if (h_high == h_low){
+          lam = 0
+        } else {
+          lam = (interp_points[i] - h_low)/(h_high - h_low) #relative position of the poi between h_low and h_high
+        }
+        ppdf = lam * matrix[, low_ind] + (1-lam)* matrix[, high_ind]
+        sed_rate_vals = transform(ppdf)
+        sed_rate_val = rej_sampling(sedrate, sed_rate_vals)
+        se[i] = sed_rate_val
+
+      }
+      return(stats::approxfun(x = interp_points, y = se, ties = function(x) sample(x, 1), rule = rule)) # for ties, randomly select one sample
+    }
+    return(f)
   }
-  return(stats::approxfun(x = interp_heights, y = se, ties = function(x) sample(x, 1), rule = rule)) # for ties, randomly select one sample
+  if (mode == "deterministic"){
+    f = function(){
+      interp_points = sort(height)
+      se = rep(NA, length(interp_points))
+      for (i in seq_along(interp_points)){
+        sed_rate_vals = transform(matrix[, i])
+        se[i] = rej_sampling(sedrate, sed_rate_vals)
+      }
+      return(stats::approxfun(x = interp_points, y = se, rule = rule))
+    }
+    return(f)
   }
-  return(f)
+  stop("unknown mode, use either \"poisson\" or \"deterministic\"")
+
 }
 
 sed_rate_gen_from_bounds = function(h_l, s_l, h_u, s_u, rate = 1){

tests/testthat/test_sedrate_gen_helpers.R

@@ -0,0 +1,27 @@
+## sedrate_from_matrix is a higher order function. The code below can be used to interactively debug it
+
+
+# ex=cycles(freqs=c(1/405.6795,1/130.719,1/123.839,1/98.86307,1/94.87666,1/23.62069, 
+#                   1/22.31868,1/19.06768,1/18.91979),end=4000,dt=5)
+# 
+# # convert to meters with a linearly increasing sedimentation rate from 0.01 m/kyr to 0.03 m/kyr
+# ex=sedRamp(ex,srstart=0.01,srend=0.03)
+# 
+# # interpolate to median sampling interval
+# ex=linterp(ex)
+# 
+# # evaluate precession & eccentricity power, and precession modulations
+# res=eTimeOpt(ex,win=20,step=1,fit=1,output=1)
+# 
+# aa = get_data_from_eTimeOpt(res, index = 3)
+# 
+# f = sed_rate_from_matrix(aa$heights, aa$sed_rate, aa$results, rate = 0.1, mode = "poisson")
+# plot(aa$heights, f()(aa$heights), type = "l")
+# 
+# 
+# height = aa$heights
+# sedrate = aa$sed_rate
+# matrix = aa$results
+# rate = 1
+# i = 1
+# transform = identity