Merge pull request #41 from MindTheGap-ERC/dev

Dev

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: admtools
 Title: Estimate and Manipulate Age-Depth Models
-Version: 0.3.0
+Version: 0.3.1
 Authors@R: 
     person("Niklas", "Hohmann", , "[email protected]", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0003-1559-1838"))

NAMESPACE

@@ -2,9 +2,13 @@
 
 S3method(condensation,adm)
 S3method(condensation,multiadm)
+S3method(get_L_tp,adm)
+S3method(get_L_tp,sac)
 S3method(get_L_unit,adm)
 S3method(get_L_unit,multiadm)
 S3method(get_L_unit,sac)
+S3method(get_T_tp,adm)
+S3method(get_T_tp,sac)
 S3method(get_T_unit,adm)
 S3method(get_T_unit,multiadm)
 S3method(get_T_unit,sac)
@@ -68,7 +72,9 @@ export(condensation_fun)
 export(flux_const)
 export(flux_linear)
 export(flux_quad)
+export(get_L_tp)
 export(get_L_unit)
+export(get_T_tp)
 export(get_T_unit)
 export(get_completeness)
 export(get_data_from_eTimeOpt)

NEWS.md

@@ -1,4 +1,8 @@
-# admtools (development version)
+# admtools 0.3.1
+
+* more options for sedimentation rate generator `sed_rate_from_matrix`
+
+* abstracted exctraction of tie points
 
 # admtools 0.3.0
 

R/get_L_tp.R

@@ -0,0 +1,29 @@
+get_L_tp = function(x, ...){
+  #' @export
+  #' 
+  #' @title get height/length tie point
+  #' 
+  #' @param x age-depth model (adm) or sediment accumulation curve (sac)
+  #' @param ... other options, currently not used
+  #' 
+  #' @description
+    #' extracts the height/length time points from an age-depth model or sediment accumulation curve
+    #' 
+  #' @returns numeric vector of the time/length tie points
+  #' 
+  #' @seealso [get_T_tp()] to extract time tie points
+  #' 
+  UseMethod("get_L_tp")
+}
+
+get_L_tp.adm = function(x, ...){
+  #' @export
+  #' 
+  return(x$h)
+}
+
+get_L_tp.sac = function(x, ...){
+  #' @export
+  #' 
+  return(x$h)
+}

R/get_T_tp.R

@@ -0,0 +1,27 @@
+get_T_tp = function(x, ...){
+  #' @export
+  #' 
+  #' @title extract time tie points
+  #' 
+  #' @param x age-depth model (adm) or sediment accumulation curve (sac)
+  #' @param ... other options, currently unused
+  #' 
+  #' @description
+    #' Extracts the time tie points from an age-depth model or sediment accumulation curve
+    #' 
+  #' @returns a vector, containing the time tie points
+  #' 
+  #' @seealso [get_L_tp()] to extract length/height tie points
+  #' 
+  UseMethod("get_T_tp")
+}
+
+get_T_tp.adm = function(x, ...){
+  #' @export
+  return(x$t)
+}
+
+get_T_tp.sac = function(x, ...){
+  #' @export
+  return(x$t)
+}

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){

R/tp_to_adm.R

@@ -22,7 +22,7 @@ tp_to_adm = function(t, h, T_unit = NULL, L_unit = NULL){
   #' 
   #' @returns object of class `adm`
   #' 
-  #' @seealso [is_adm()] to check validity of `adm` objects
+  #' @seealso [is_adm()] to check validity of `adm` objects, [get_T_tp()] and [get_L_tp()] to extract time and height/length tie points
   #'
   #'
   #' @examples

R/tp_to_sac.R

@@ -12,6 +12,8 @@ tp_to_sac = function(t, h, T_unit = NULL, L_unit = NULL){
   #' @param L_unit length unit
   #' 
   #' @returns a _sac_ object reflecting a sediment accumulation curve
+  #' 
+  #' @seealso [sac_to_adm()] to transform sediment accumulation curves into age-depth models,  [get_T_tp()] and [get_L_tp()] to extract time and height/length tie points
 
   li = list(t = t,
             h = h,

tests/testthat/test_get_L_tp.R

@@ -0,0 +1,10 @@
+test_that("correct heights are returned",{
+  h = 1:4
+  t = 1:4
+  adm = tp_to_adm(t, h)
+  expect_equal(get_L_tp(adm), h)
+
+  h = rev(h)
+  sac = tp_to_sac(t, h)
+  expect_equal(get_L_tp(sac), h)
+})

tests/testthat/test_get_T_tp.R

@@ -0,0 +1,10 @@
+test_that("correct times are returned",{
+  h = 1:4
+  t = 1:4
+  adm = tp_to_adm(t, h)
+  expect_equal(get_T_tp(adm), t)
+
+  h = rev(h)
+  sac = tp_to_sac(t, h)
+  expect_equal(get_T_tp(sac), t)
+})

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

vignettes/admtools.Rmd

@@ -141,7 +141,7 @@ If you want to inspect the insides of the object, use `str`:
 str(my_adm)
 ```
 
-You can manually manipulate the fields of the `adm` object by treating it like a list. I do not recommend doing so, as it might result in unexpected downstream behavior.
+You can manually manipulate the fields of the `adm` object by treating it like a list. I do not recommend doing so, as it might result in unexpected downstream behavior. If you want to extract tie points use `get_L_tp` and `get_T_tp`.
 
 You can plot `adm` objects via the standard `plot` function. Here, I use the option to highlight hiatuses in red, and increase the linw width of the conservative ( = non-destructive) intervals.