Merge pull request #1 from MindTheGap-ERC/dev

Dev

.Rbuildignore

@@ -1,3 +1,4 @@
 ^StratPal\.Rproj$
 ^\.Rproj\.user$
 ^LICENSE$
+^\.github$

.github/.gitignore

@@ -0,0 +1 @@
+*.html

.github/workflows/R-CMD-check.yaml

@@ -0,0 +1,52 @@
+# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples
+# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
+on:
+  push:
+    branches: [main, master]
+  pull_request:
+    branches: [main, master]
+
+name: R-CMD-check
+
+permissions: read-all
+
+jobs:
+  R-CMD-check:
+    runs-on: ${{ matrix.config.os }}
+
+    name: ${{ matrix.config.os }} (${{ matrix.config.r }})
+
+    strategy:
+      fail-fast: false
+      matrix:
+        config:
+          - {os: macos-latest,   r: 'release'}
+          - {os: windows-latest, r: 'release'}
+          - {os: ubuntu-latest,   r: 'devel', http-user-agent: 'release'}
+          - {os: ubuntu-latest,   r: 'release'}
+          - {os: ubuntu-latest,   r: 'oldrel-1'}
+
+    env:
+      GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+      R_KEEP_PKG_SOURCE: yes
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - uses: r-lib/actions/setup-pandoc@v2
+
+      - uses: r-lib/actions/setup-r@v2
+        with:
+          r-version: ${{ matrix.config.r }}
+          http-user-agent: ${{ matrix.config.http-user-agent }}
+          use-public-rspm: true
+
+      - uses: r-lib/actions/setup-r-dependencies@v2
+        with:
+          extra-packages: any::rcmdcheck
+          needs: check
+
+      - uses: r-lib/actions/check-r-package@v2
+        with:
+          upload-snapshots: true
+          build_args: 'c("--no-manual","--compact-vignettes=gs+qpdf")'

DESCRIPTION

@@ -10,7 +10,7 @@ Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.3.2
 Imports: 
-    admtools
+    admtools (>= 0.3.0)
 Suggests: 
     knitr,
     rmarkdown

R/apply_niche_pref.R

@@ -1,13 +1,44 @@
-apply_niche_pref = function(x, niche, temp_change){
+apply_niche_pref = function(x, niche_def, gc){
   #' @export
   #'
-  #' @title apply niche preference to events
+  #' @title apply niche preference
   #'
-  #' @param x events
-  #' @param niche function, defining a niche along a gradient
-  #' @param temp_change function, change in the niche with time
+  #' @param x events, e.g. fossil occurrences
+  #' @param niche_def function, specifying the niche along a gradient. should return 0 when taxon is outside of niche, and 1 when fully inside niche. Values between 0 and 1 are interpreted as probabilities.
+  #' @param gc function, stands for "gradient change". Specifies how the gradient changes, e.g. with time
   #'
-  niche_with_time = function(t) niche(temp_change(t))
+  #' @description
+    #' models niche preferences by removing events (fossil occurrences) when they are outside of their preferred niche using the function `thin`
+    #' Combines the functions niche_def and gc ("gradient change") to determine how the taxons preferred environment changes with time. This is done by composing `niche_def` and `gc`. The result is then used as a thinning.
+  #'
+  #' @examples
+    #' \dontrun{
+    #' ## setup
+    #' # using water depth as gradient
+    #'  t = scenarioA$t_myr
+    #'  wd = scenarioA$wd_m[,"8km"]
+    #'  gc = approxfun(t, wd)
+    #'  plot(t, gc(t), type = "l", xlab = "Time", ylab = "water depth [m]",
+    #'   main = "gradient change with time")
+    #'  # define niche
+    #'  # preferred wd 10 m, tolerant to intermediate wd changes (standard deviation 10 m), non-terrestrial
+    #'  niche_def = cnd(mean = 10, sd = 10, inc  = 40, cut_neg = TRUE)
+    #'  plot(seq(-1, 50, by = 0.5), niche_def(seq(-1, 50, by = 0.5)), type = "l",
+    #'  xlab = "water depth", ylab = "preference", main = "Niche def")
+    #'  # niche pref with time
+    #'  plot(t, niche_def(gc(t)), type = "l", xlab = "time", ylab = "preference", main = "pref with time")
+    #'
+    #'  ## simulate fossil occurrences
+    #'  foss_occ = p3(rate = 100, from = 0, to = max(t))
+    #'  # foss occ without niche pref
+    #'  hist(foss_occ, xlab = "time")
+    #'  foss_occ_niche = apply_niche_pref(foss_occ, niche_def, gc)
+    #'  # fossil occurrences with niche preference
+    #'  hist(foss_occ_niche, xlab = "time")
+    #'
+    #' }
+    #'
+  niche_with_time = function(t) niche_def(gc(t))
   r = thin(x, niche_with_time)
   return(r)
   }

R/cnd.R

@@ -11,6 +11,19 @@ cnd = function(mean, sd, inc = 1, cut_neg = TRUE){
   #'
   #' @returns a function
   #'
+  #' @examples
+    #' \dontrun{
+    #' # using water depht as niche
+    #' wd = seq(-3, 40, by = 0.5)
+    #' f = cnd(mean = 10, sd = 5, inc = 15, cut_neg = FALSE)
+    #' # 1 indicates high preference, 0 indicates low preference
+    #' plot(wd, f(wd), xlab = "Water depth", ylab = "Env. preference")
+    #' # set value at neg wd to 0 - non-terrestrial species.
+    #' f = cnd(mean = 10, sd = 5, inc = 15, cut_neg = TRUE)
+    #' plot(wd, f(wd), xlab = "Water depth", ylab = "Env. preference")
+    #' }
+    #'
+  #'
   #' @description
     #' returns a function that defines a niche based on a capped normal distribution, i.e. a pdf of a normal distribution where all values above 1 are capped. Mathematically, this is f(x) = min( inc * pdf(x), 1)
 

R/ornstein_uhlenbeck.R

@@ -13,6 +13,17 @@ ornstein_uhlenbeck = function(t, mu = 0, theta = 1, sigma = 1, y0 = 0){
     #' simulates ornstein-uhlenbeck using the Euler-Maruyame method. The process is simulated on a scale of 0.25 * min(diff(t)) and then interpolated to the values of `t`.
     #'
   #' @returns a list with two elements: `t` and `y`. `t` is a duplicate of the input `t`, `y` are the values of the OU process at these times. Outputs are of class `timelist` and can thus be plotted directly using `plot`, see `?plot.timelist`
+  #'
+  #' @examples
+    #' \dontrun{
+    #' library("admtools") # required for plotting of results
+    #' t = seq(0, 3, by = 0.01)
+    #' l = ornstein_uhlenbeck(t, y0 = 3) # start away from optimum (mu)
+    #' plot(l, type = "l")
+    #' l2 = ornstein_uhlenbeck(t, y0 = 0) # start in optimum
+    #' lines(l2$t, l2$y, col = "red")
+    #' }
+    #'
 
   if (y0 == "stationary"){
     y0 = stats::rnorm(1, mean = mu, sd = sigma / sqrt(2 * theta))

R/p3.R

@@ -2,16 +2,28 @@ p3 = function(rate, from, to, n = NULL){
   #'
   #' @export
   #'
-  #' @title  simulate fossil occurrences
+  #' @title  simulate Poisson point process
   #'
-  #' @param rate rate of fossil occurrences (avg. no of fossils per unit)
-  #' @param from lowest possible fossil occurrence
-  #' @param to highest possible fossil occurrence
-  #' @param n integer of NULL (default). Number of fossil occurrences to return. If null, the number is random and determined by the rate parameter
+  #' @param rate strictly positive scalar, rate of events (avg events per unit)
+  #' @param from lowest boundary of observed interval
+  #' @param to upper boundary of observed interval
+  #' @param n integer of NULL (default). Number of events to return. If NULL, the number is random and determined by the rate parameter
   #'
   #' @description
-    #' Simulates fossil occurrences in the interval from, to based on a Poisson point process with rate `rate`. If the parameter `n` is used, the number of fossils is conditioned to be `n`
+    #' Simulates events in the interval `from` to `to` based on a Poisson point process with rate `rate`. If the parameter `n` is used, the number of fossils is conditioned to be `n`
+    #' In the context of paleontology, these events can be interpreted as fossil occurrences or first/last occurrences of species. In this case, the rate is the average number of fossil occurrences (resp first/last occurrences) per unit
+  #' @examples
+    #' \dontrun{
+    #' # for fossil occ.
+    #' x = p3(rate = 5, from = 0, to = 1) # 5 fossil occurrences per myr on avg.
+    #' hist(x, xlab = "Time (Myr"), ylab = "Fossil Occurrences" )
     #'
+    #' x = p3(rate = 3, from = 0, to = 4)
+    #' hist(x, main = paste0(length(x), " samples")) # no of events is random
+    #'
+    #' x = p3(rate = 3, from = 0, to = 4, n = 10)
+    #' hist(x, main = paste0(length(x), " samples")) # no of events is fixed to n
+    #' }
     #'
 
   if (rate <= 0){ stop("Need strictly positive rate")}

R/p3_var_rate.R

@@ -13,14 +13,16 @@ p3_var_rate = function(x, y = NULL, from = 0, to = 1, f_max = 1, n = NULL){
   #'
   #' @description
     #' simulates events based on a variable rate Poisson point process. Rates can be either specified by a function passed to `x`, or by providing two vectors `x` and `y`. In this case the rate is specified by approxfun(x, y, rule = 2), i.e. by linear interpolation between the values of x (abscissa) and y (ordinate)
-    #'
+    #' In the context of paleontology, these events can be interpreted as fossil occurrences or first/last occurrences of species. In this case, the rate is the average number of fossil occurrences (resp first/last occurrences) per unit
   #' @examples
     #' \dontrun{
+    #' # assuming events are fossil occurrences
+    #' # then rate is the avg rate of fossil occ. per unit
     #' linear decrease in rate from 50 at x = 0 to 0 at x = 1
     #' x = c(0, 1)
     #' y = c(50, 0)
     #' s = p3_var_rate(x, y, f_max = 50)
-    #' hist(s)
+    #' hist(s, xlab = "Time (myr)", main = "Fossil Occurrences")
     #' # conditoned to return 100 samples
     #' s = p3_var_rate(x, y, f_max = 50, n = 100)
     #' # hand over function

R/random_walk.R

@@ -13,6 +13,16 @@ random_walk = function(t, sigma = 1, mu = 0, y0 = 0){
     #'
   #' @returns a list with elements `t` and `y`. `t` is a duplicate of the input parameter and is the times at which the random walk is evaluated. `y` are the values of the random walk at said times
   #'
+  #' @examples
+    #' \dontrun{
+    #' library("admtools") # required for plotting of results
+    #' t = seq(0, 1, by = 0.01)
+    #' l = random_walk(t, sigma = 3) # high variability, no direction
+    #' plot(l, type = "l")
+    #' l2 = random_walk(t, mu = 1) # low variabliity, increasing trend
+    #' lines(l2$t, l2$y, col = "red")
+    #' }
+    #'
   increments = diff(t)
   acc = cumsum(c(0, stats::rnorm(n = increments, mean = 0, sd = sqrt(increments))))
   y = sigma * acc + mu * (t-min(t)) + y0

R/rej_samp.R

@@ -4,9 +4,9 @@ rej_samp = function(f, x_min, x_max, n = 1L, f_max = 1){
   #' @title rejection sampling
   #'
   #' @description
-    #' rejection sampling from the (pseudo) pdf `f` in the interval between `x_min` and `x_max`. Returns `n` samples. Note that values of `f` below 0 are rounded up to zero
+    #' rejection sampling from the (pseudo) pdf `f` in the interval between `x_min` and `x_max`. Returns `n` samples. Note that values of `f` below 0 are capped to zero
     #'
-  #' @param f function. (pseude) pdf from which the sample is drawn
+  #' @param f function. (pseudo) pdf from which the sample is drawn
   #' @param x_min scalar. lower limit of the examined interval
   #' @param x_max scalar. upper limit of the examined interval
   #' @param n integer. number of samples drawn

R/stasis.R

@@ -10,6 +10,16 @@ stasis = function(t, mean = 0, sd = 1){
   #' @description
     #' simulates stasis as independent, normally distributed random variables with mean `mean` and standard deviatin `sd`
     #'
+  #' @examples
+    #' \dontrun{
+    #'  library("admtools") # required for plotting of results
+    #'  t = seq(0, 1, by = 0.01)
+    #'  l = stasis(t)
+    #'  plot(l, type = "l") # plot lineage
+    #'  l2 = stasis(t, mean = 0.5, sd = 0.3) # simulate second lineage
+    #'  lines(l2$t, l2$y, col = "red") # plot second lineage
+    #' }
+    #'
   l = list(t = t,
            y = stats::rnorm(n = length(t), mean = mean, sd = sd))
   class(l) = c("timelist", "list")

README.md

@@ -1,5 +1,10 @@
 # StratPal
 
+<!-- badges: start -->
+[![R-CMD-check](https://github.com/MindTheGap-ERC/StratPal/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/MindTheGap-ERC/StratPal/actions/workflows/R-CMD-check.yaml)
+<!-- badges: end -->
+
+
 R package for Stratigraphic Paleobiology. IN DEVELOPMENT
 
 ## Authors