R/calc.R

Defines functions fromslump toslump agemodel.it Bacon.rng Bacon.hist hiatus.slopes Bacon.d.Age Bacon.Age.d

Documented in agemodel.it Bacon.Age.d Bacon.d.Age Bacon.hist

#' @name Bacon.Age.d
#' @title Output all ages for a single depth.
#' @description Output all MCMC-derived age estimates for a given depth.
#' @details Obtaining an age-depth model is often only a step towards a goal, e.g., plotting a core's
#' fossil series ('proxies') against calendar time. Bacon.Age.d can be used to list all MCMC-derived age estimates for a given (single) depth, for example to calculate mean ages for a depth. See also Bacon.d.Age which calculates the depths of a single age estimate.
#' @param d The depth of which Bacon age estimates are to be returned. Has to be a single depth.
#' @param set Detailed information of the current run, stored within this session's memory as variable info.
#' @param its The set of MCMC iterations to be used. Defaults to the entire MCMC output, \code{its=set$output}.
#' @param BCAD The calendar scale of graphs and age output-files is in \code{cal BP} by default, but can be changed to BC/AD using \code{BCAD=TRUE}.
#' @param na.rm Whether or not to remove NA values (ages within slumps)
#' @author Maarten Blaauw, J. Andres Christen
#' @return Outputs all MCMC-derived ages for a given depth.
#' @examples
#' \dontrun{
#'   Bacon(run=FALSE, coredir=tempfile())
#'   agedepth(age.res=50, d.res=50, d.by=10)
#'   ages.d20 = Bacon.Age.d(20)
#'   mean(ages.d20)
#' }
#' @export
Bacon.Age.d <- function(d, set=get('info'), its=set$output, BCAD=set$BCAD, na.rm=FALSE) {
  if(length(d) > 1)
    stop("Bacon.Age.d can handle one depth at a time only", call.=FALSE)
  if(length(its) == 0)
    stop("core's data not yet in memory. Please run agedepth first\n", call.=FALSE)

  hiatus.depths <- set$hiatus.depths
  if(length(set$slump) > 0) {
   d <- toslump(d, set$slump, remove=na.rm)
   if(!is.na(hiatus.depths[1]))
     hiatus.depths <- set$slumphiatus
  }

  ages <- numeric(nrow(its)) # to sort R-base problem with c() and loops
  if(!is.na(d))
    if(d >= min(set$elbows)) { # we cannot calculate ages of depths above the top depth; was > d.min before Feb 2021
      topages <- as.vector(its[,1]) # ages for the core top
      maxd <- max(which(set$elbows <= d)) # find the relevant sections
      accs <- as.matrix(its[,1+(1:maxd)]) # the accumulation rates xi for each relevant section
      cumaccs <- cbind(0, t(apply(accs, 1, cumsum))) # cumulative accumulation
      ages <- topages + (set$thick * cumaccs[,maxd]) + # topages + xi * dC + ...
	   ((d-set$elbows[maxd]) * accs[,maxd]) # ... remaining bit of lowest section

      # now using a uniform jump, not gamma
      if(!is.na(hiatus.depths[1]))
        for(i in 1:length(hiatus.depths)) {
          above <- max(which(set$elbows < hiatus.depths[i]), 1)[1]
          below <- above + 1
          if(d > set$elbows[above] && d <= set$elbows[below]) { # adapt ages for sections with hiatus
            if(d > hiatus.depths[i]) # April 2020, changed snippets below from [[i]]] to [,i]
              ages <- set$elbow.below[,i] - (set$slope.below[,i] * (set$elbows[below] - d)) else
                ages <- set$elbow.above[,i] + (set$slope.above[,i] * (d - set$elbows[above]))
            }
        }
    } else
         ages <- NA # Feb 2021

    if(BCAD)
      ages <- 1950 - ages
    return(c(ages))
}



#' @name Bacon.d.Age
#' @title Output all depths for a single age.
#' @description Output all depths of a single given MCMC-derived age estimate.
#' @details Obtaining an age-depth model is often only a step towards a goal, e.g., plotting a core's
#' fossil series ('proxies') against calendar time. Bacon.d.Age can be used to list all MCMC-derived depths belonging to a given (single) age, for example to calculate mean depths belonging to a modelled depth. 
#' This function was kindly written and provided by Timon Netzel (Bonn University). See also Bacon.Age.d, which calculates the ages for a single depth.
#' @param age The age estimate for which depths are to be returned. Has to be a single age.
#' @param set Detailed information of the current run, stored within this session's memory as variable info.
#' @param its The set of MCMC iterations to be used. Defaults to the entire MCMC output, \code{its=set$output}.
#' @param BCAD The calendar scale of graphs and age output-files is in \code{cal BP} by default, but can be changed to BC/AD using \code{BCAD=TRUE}.
#' @param na.rm Whether or not to remove NA values (ages within slumps)
#' @author Maarten Blaauw, J. Andres Christen
#' @return Outputs all MCMC-derived ages for a given depth.
#' @examples
#' \dontrun{
#'   Bacon(run=FALSE, coredir=tempfile())
#'   agedepth(age.res=50, d.res=50, d.by=10)
#'   ages.d20 = Bacon.Age.d(20)
#'   mean(ages.d20)
#' }
#' @export
Bacon.d.Age <- function(age, set=get("info"), BCAD=set$BCAD, its=set$output, na.rm=FALSE ){

  if(BCAD) 
    age <- 1950 - age # if the customer decides to enter the age in BCAD
  if(length(age) > 1) 
    stop("Bacon.d.Age can handle one age at a time only", call. = FALSE)
  if(length(its) == 0 ) 
    stop("core's data not yet in memory. Please run agedepth first\n", call. = FALSE)

  # create elbow ages
  elbows <- set$elbows
  topages <- as.vector(its[,1]) # ages for the core top
  accs <- as.matrix(its[,1+(1:set$K)]) # the accumulation rates xi for each section
  cumaccs <- set$thick * cbind(0, t(apply(accs, 1, cumsum))) # cumulative accumulation
  elbow.ages <- topages + cumaccs

  if ( age <= min(elbow.ages) || age > max(elbow.ages) ) 
    message(" Warning, age outside the core's age range!\n")

  # prepare slump
  hiatus.depths <- set$hiatus.depths
  if(length(set$slump) > 0) {
    diff.slumps <- apply(set$slump,1,diff,na.rm=T)  # vector containing the difference of each slump
    if ( !is.na(hiatus.depths[1]) ) 
      hiatus.depths <- set$slumphiatus
  }
    
    # prepare hiatus
    hiatus.check <- rep(0,ncol(elbow.ages)) # vector containing the information for each elbow whether there is a hiatus (1) or no hiatus (0)
    if ( !is.na(hiatus.depths[1]) ){ 
        # if there is no slump
        if ( length(set$slump) == 0 ) 
            hiatus.depths <- set$hiatus.depths
        # where are the hiatuses
        for ( i in 1:length(hiatus.depths) ) {
            hiatus.check[ min(which(elbows > hiatus.depths[i])) ] <- 1
        }
    }

    # intersections which are below the last elbow
    these.bottom <- which( elbow.ages[, ncol(elbow.ages)] < age ) 
    
    # summary of the depths corresponding to the chosen age
    depths <- rep(NA,nrow(its))
    for ( i in 2:ncol(elbow.ages) ){
        # consideration of hiatuses
        if ( hiatus.check[i] == 1 ){ 
            for ( j in 1:length(hiatus.depths) ) {
                above <- max(which(elbows < hiatus.depths[j]), 1)[1]
                below <- above + 1
                # two age ranges for one hiatus (below/above) with the adapted slopes (could be included in the info list)
                ages.above <- set$elbow.above[,j] + (set$slope.above[,j] * (hiatus.depths[j] - elbows[above]))
                ages.below <- set$elbow.below[,j] - (set$slope.below[,j] * (elbows[below] - hiatus.depths[j]))
                # intersections below and above
                these.above <- (elbow.ages[, above] < age) * (ages.above > age) 
                these.below <- (ages.below < age) * (elbow.ages[, below] > age) 
                # find the depths which correspond to the intersections below and above the hiatus
                if ( sum(these.above, na.rm=T) > 0 ){
                    # slopes above hiatus
                    accs.above <- (1/set$slope.above[,j])[which(these.above==1)]
                    # age difference with the last age range above the hiatus
                    age.diff.above <- age - elbow.ages[which(these.above==1), above]
                    # depths: last depth above hiatus depth + slopes * age difference
                    depths[which(these.above==1)] <- elbows[above] + accs.above * age.diff.above
                }
                if ( sum(these.below, na.rm=T) > 0 ){
                    # slopes below hiatus
                    accs.below <- (1/set$slope.below[,j])[which(these.below==1)]
                    # age difference with the lower hiatus age range 
                    age.diff.below  <- age - ages.below[which(these.below==1)]
                    # depths: hiatus depth + slopes * age difference 
                    depths[which(these.below==1)] <- hiatus.depths[j] + accs.below * age.diff.below
                }
            }
        } else {
            # consideration of the cases without hiatuses
            
            # find the correct intersections
            these <- (elbow.ages[, i - 1] < age) * (elbow.ages[, i] > age)

            # find the depths which correspond to these intersections
            if ( sum(these, na.rm=T) > 0 ){ 
                # which slopes
                accs <- 1/set$output[which(these == 1), i]
                # difference between the age and the elbow ages found
                age.diff <- age - elbow.ages[which(these == 1), i - 1]
                # depths: elbow depth + slopes *  age difference
                depths[which(these==1)] <- elbows[i-1] + accs * age.diff
            }
            
            # find these intersections which are below the last elbow
            if( length(these.bottom)  > 0 & i == ncol(elbow.ages) ){
                # slopes below last elbow
                accs.bottom <- 1/set$output[these.bottom, i]
                # difference between the age and the last elbow ages found
                age.diff.bottom <-  age - elbow.ages[these.bottom, i ]
                # depth: last elbow + slopes *  age difference 
                depths[these.bottom] <- elbows[i-1] + accs.bottom * age.diff.bottom
            }
        }  
    }
    
    # Query: remove NA entries if some entries are not NA (happens for some hiatuses).
    # In these cases the depth vector does not have the same length for different age.
    # In cases where the age entered is younger than min(elbow.ages), each entry in depth is NA and is returned as such.
    if( !all(is.na(depths)) & na.rm ) depths <- depths[!is.na(depths)]
 
    # consideration of slumps
    if ( length(set$slump) > 0 ) {
        for( j in 1:nrow(set$slump) ){
            # which depths are below the slump
            these.below.slump <- (depths >= set$slump[j,1]) 
            below.slump <- sum(these.below.slump, na.rm = T)
            if ( below.slump > 0 ){
                # shift the found depths with the slump depths above
                these.shift <- which(these.below.slump==1)
                depths[these.shift] <- depths[these.shift] + diff.slumps[j]  
            }
        }
    }
 
    # output
    return(depths)
}



# First get ages and slopes of c's just above and below hiatus.
# then calculate slope.below and slope.above as extrapolations from the sections below resp. above (option 1, default), no adapting of slopes (option 0),
# or as w-weighted mix of prior and accrates, resp. prior only (option 2).
# then check if these slopes work (no reversals). Those with reversals revert to the original slopes.
# check approach for boundaries. check that elbows in correct order, then set slopes to either ages.below or ages.above???
hiatus.slopes <- function(set=get('info'), hiatus.option=1) {
  elbows <- set$elbows
  hiatus.depths <- set$hiatus.depths
  if(length(set$slump) > 0)
    hiatus.depths <- set$slumphiatus

  its <- cbind(set$output)
  w <- its[,ncol(its)]^(1/set$thick)
  fillvals <- array(NA, dim=c(nrow(its), length(hiatus.depths))) # 17 Apr 2020
#  set$slope.below <- numeric(nrow(its))  # 17 Apr 2020
#  set$slope.above <- numeric(nrow(its))  # 17 Apr 2020
  set$slope.below <- fillvals  # 17 Apr 2020
  set$slope.above <- fillvals  # 17 Apr 2020
  set$elbow.below <- fillvals # 17 Apr 2020
  set$elbow.above <- fillvals # 17 Apr 2020
  set$above <- numeric(1)
  topages <- as.vector(its[,1]) # ages for the core top
  accs <- as.matrix(its[,1+(1:set$K)]) # the accumulation rates xi for each section
  cumaccs <- set$thick * cbind(0, t(apply(accs, 1, cumsum))) # cumulative accumulation
  elbow.ages <- topages + cumaccs

  for(i in 1:length(hiatus.depths)) {
    above <- max(which(elbows < hiatus.depths[i]), 1) ### is this the correct one?
      elbow.below <- elbow.ages[,above+1]
      elbow.above <- elbow.ages[,above]
      orig.slope <- accs[,above]

    if(hiatus.option == 0) { # then do nothing
      slope.below <- orig.slope
      slope.above <- orig.slope
    }
    if(hiatus.option == 1) { # then extrapolate slopes above/below section w hiatus
      slope.below <- its[,above+2]
      slope.above <- its[,above]
    }
    if(hiatus.option == 2) { # then w-weighted for below, and prior-only for above
      slope.below <- w*set$output[,above+2] + (1-w)*set$acc.mean[i+1]
      slope.above <- rep(set$acc.mean[i], nrow(its))
    }

    # now calculate the ages at the hiatus/boundary, coming from below and from above
    ages.above <- elbow.above + (slope.above * (hiatus.depths[i] - elbows[above]))
    ages.below <- elbow.below - (slope.below * (elbows[above+1] - hiatus.depths[i]))

    if(!is.na(set$boundary[1])) { # then set the boundary's elbow at ages.below for both sections
      if(length(set$slumpboundary) > 0)
        boundary <- set$slumpboundary else
          boundary <- set$boundary
      ages.boundary <- elbow.below - (slope.below * (elbows[above+1] - set$boundary[i]))
      slope.above <- (ages.boundary - elbow.above) / (set$boundary[i] - elbows[above])
      slope.below <- (ages.below - ages.boundary) / (elbows[above+1] - set$boundary[i])
    }

    # for sections with reversals, use the original slopes
    reversed <- c(which(elbow.above > ages.above),
      which(ages.above > ages.below),
      which(ages.below > elbow.below),
      which(slope.below < 0), which(slope.above < 0))
   if(length(reversed) > 0) {
     slope.above[reversed] <- orig.slope[reversed]
     slope.below[reversed] <- orig.slope[reversed]
   }

    # store the updated information
    # set$elbow.below[[i]] <- elbow.below # commented Apr 2020
    # set$elbow.above[[i]] <- elbow.above
    # set$slope.below[[i]] <- slope.below
    # set$slope.above[[i]] <- slope.above
    set$elbow.below[,i] <- elbow.below # new April 2020
    set$elbow.above[,i] <- elbow.above # new 
    set$slope.below[,i] <- slope.below # new
    set$slope.above[,i] <- slope.above # new
    
    set$above <- above
  }
  return(set)
}



#' @name Bacon.hist
#' @title Calculate age distributions of depths.
#' @description Calculate the distribution of age estimates of single or multiple depths.
#' @details Age estimates of specific depths can also be plotted.
#' @param d The depth or depths for which a histogram and age ranges should be provided. If multiple depths are given, then just the age ranges, median and means (no graphs) are provided for each depth.
#' @param set Detailed information of the current run, stored within this session's memory as variable info.
#' @param BCAD The calendar scale of graphs and age output-files is in \code{cal BP} by default, but can be changed to BC/AD using \code{BCAD=TRUE}.
#' @param age.lab The labels for the calendar axis (default \code{age.lab="cal BP"} or \code{"BC/AD"} if \code{BCAD=TRUE}).
#' @param age.lim Minimum and maximum calendar age ranges, calculated automatically by default (\code{age.lim=c()}).
#' @param hist.lab The y-axis is labelled \code{ylab="Frequency"} by default.
#' @param calc.range Calculate ranges? Takes time so can be left out
#' @param hist.lim Limits of the y-axis.
#' @param draw Draw a plot or not. Defaults to \code{draw=TRUE}, however no plots are made if more than one depth \code{d} is provided.
#'  If \code{draw=FALSE}, then the age ranges, median and mean are given for each depth (as four columns).
#' @param prob Age ranges are given as quantiles, e.g., 2.5\% and 97.5\% for the default of 95\% confidence limits (\code{prob=0.95})).
#' @param hist.col Colour of the histogram. Default grey, \code{hist.col=grey(0.5)}.
#' @param hist.border Colour of the histogram's outline. Default dark grey, \code{hist.border=grey(0.2)}.
#' @param range.col Colour of confidence ranges. Defaults to \code{range.col="blue"}.
#' @param med.col Colour of the median. Defaults to \code{med.col="green"}.
#' @param mean.col Colour of the mean. Defaults to \code{mn.col="red"}.
#' @param verbose Provide feedback on what is happening (default \code{verbose=TRUE}).
#' @author Maarten Blaauw, J. Andres Christen
#' @return A plot with the histogram and the age ranges, median and mean, or just the age ranges, medians and means if more than one depth \code{d} is given.
#' @examples
#' \dontrun{
#'   Bacon(run=FALSE, coredir=tempfile())
#'   agedepth(age.res=50, d.res=50, d.by=10)
#'   Bacon.hist(20)
#'   Bacon.hist(20:30)
#' }
#' @export
Bacon.hist <- function(d, set=get('info'), BCAD=set$BCAD, age.lab=c(), age.lim=c(), hist.lab="Frequency", calc.range=TRUE, hist.lim=c(), draw=TRUE, prob=set$prob, hist.col=grey(0.5), hist.border=grey(.2), range.col="blue", med.col="green", mean.col="red", verbose=TRUE) {
  outfile <- paste(set$prefix, ".out", sep="")
  if(length(set$output) == 0 || length(set$Tr) == 0) {
    set <- Bacon.AnaOut(outfile, set)
    assign_to_global("set", set)
  }
  hist3 <- function(d, BCAD) {
    hsts <- list(); maxhist <- 0; minhist <- 1
    pb <- txtProgressBar(min=0, max=max(1,length(d)-1), style = 3)
    for(i in 1:length(d)) {
      if(length(d) > 1)
        setTxtProgressBar(pb, i)
      ages <- Bacon.Age.d(d[i], set, BCAD=BCAD)
      if(length(ages[!is.na(ages)]) > 0) { # added !is.na(ages) 21 April 21
        hst <- density(ages)
        th0 <- min(hst$x)
        th1 <- max(hst$x)
        maxhist <- max(maxhist, hst$y)
        minhist <- min(minhist, max(hst$y))
        n <- length(hst$x)
        counts <- hst$y
        ds <- d[i]
        hsts <- append(hsts, pairlist(list(d=ds, th0=th0, th1=th1, n=n, counts=counts, max=maxhist, min=minhist)))
      } else hsts$d[[i]] <- d[i]
    }
    return(hsts)
  }
  hists <- hist3(d, BCAD)
  assign_to_global("hists", hists)

  rng <- array(NA, dim=c(length(d), 4)) # R > 4.0 does not like to fill c() using loops
  if(calc.range)
    rng <- Bacon.rng(d, set, BCAD, prob)

  if(length(d)==1)
    if(draw==TRUE) {
      hst <- hists[[1]]
      if(length(age.lab) == 0)
        age.lab <- ifelse(BCAD, "BC/AD", "cal yr BP")
      if(length(age.lim) == 0)
        age.lim <- c(hst$th0, hst$th1)
      if(BCAD)
        age.lim <- rev(age.lim)
      if(length(hist.lim) == 0)
        hist.lim <- c(0, 1.1*hst$max)

      pol <- cbind(c(hst$th0, seq(hst$th0, hst$th1, length=hst$n), hst$th1), c(0, hst$counts, 0))
      plot(0, type="n", xlim=age.lim, ylim=hist.lim, xlab=age.lab, ylab=hist.lab, yaxs="i")
      polygon(pol, col=hist.col, border=hist.border)
      segments(rng[,1], 0, rng[,2], 0, col=range.col, lwd=3)
      points(rng[,3], 0, col=med.col, pch=20)
      points(rng[,4], 0, col=mean.col, pch=20)

      if(verbose) {
        message("mean (", mean.col, "): ", round(rng[4],1), " ", age.lab,
          ", median (", med.col, "): ",  round(rng[3],1), " ", age.lab, "\n")
        message(100*prob, "% range (", range.col, "): ", round(rng[1],1), " to ", round(rng[2],1), " ", age.lab, "\n")
      }
    }
  invisible(rng)
}


# to calculate age ranges
Bacon.rng <- function(d, set=get('info'), BCAD=set$BCAD, prob=set$prob) {
  outfile <- paste(set$prefix, ".out", sep="")
  if(length(set$output) == 0 || length(set$Tr) == 0) {
    set <- Bacon.AnaOut(outfile, set)
    assign_to_global("set", set)
  }

  if(length(d) > 1)
    pb <- txtProgressBar(min=0, max=max(1, length(d)-1), style=3)
  rng <- array(NA, dim=c(length(d), 4))

  for(i in 1:length(d)) {
    ages <- Bacon.Age.d(d[i], set, BCAD=BCAD)
    ages <- ages[!is.na(ages)]
    if(length(ages) > 0) {
      rng[i,1:3] <- quantile(ages, c(((1-prob)/2), 1-((1-prob)/2), .5), na.rm=TRUE)
      rng[i,4] <- mean(ages)
    }

  if(length(d) > 1)
    setTxtProgressBar(pb, i)
  }
  #if(length(d) > 1)
  #  close(pb)
  return(rng)
}

#' @name agemodel.it
#' @title Extract one age-model iteration
#' @description For one MCMC iteration (it), extract the corresponding age-depth model.
#' @param it The MCMC iteration of which the age-model should be calculated.
#' @param set Detailed information of the current run, stored within this session's memory as variable info.
#' @param BCAD The calendar scale of graphs and age output-files is in \code{cal BP} by default, but can be changed to BC/AD using \code{BCAD=TRUE}.
#' @author Maarten Blaauw, J. Andres Christen
#' @return A variable with two columns - depth and the age-depth model of a single iteration.
#' @examples
#' \dontrun{
#'   Bacon(run=FALSE, coredir=tempfile())
#'   agedepth(age.res=50, d.res=50, d.by=10)
#'   lines(agemodel.it(5), col="red")
#' }
#' @export
agemodel.it <- function(it, set=get('info'), BCAD=set$BCAD) {
  outfile <- paste(set$prefix, ".out", sep="")
  if(length(set$output) == 0 || length(set$Tr) == 0) {
    set <- Bacon.AnaOut(outfile, set)
    assign_to_global("set", set)
  }
  # does this function work in cores with slumps? also doesn't take into account hiatuses.
  if(length(set$hiatus.depths) > 0)
    age <- sort(c(set$d, set$hiatus.depths+.001, set$hiatus.depths))
  age <- c()
  for(i in 1:length(set$elbows))
    age[i] <- Bacon.Age.d(set$elbows[i], set, BCAD=BCAD)[it]
  cbind(set$elbows, age)
}


# calculate slumpfree depths
toslump <- function(d, slump, remove=FALSE) {
  d <- sort(d)
  slump <- matrix(sort(slump), ncol=2, byrow=TRUE)
  slices <- c(0, slump[,2] - slump[,1])
  dfree <- d
  for(i in 1:nrow(slump)) {

    inside <- which(d <= slump[i,2]) # find depths within slump, step 1
    inside <- which(d[inside] >= slump[i,1]) # step 2
    below <- which(d >= slump[i,2]) # adapt depths below slumps

    if(length(below) > 0) # depths below slump
      dfree[below] <- dfree[below] - slices[i+1]

    if(length(inside) > 0) # depths within slump
      if(min(d) < max(slump[i,]))
        if(remove)
          dfree[inside] <- NA else
            dfree[inside] <- slump[i,1] - sum(slices[1:i])
  }
  return(dfree)
}



# calculate original depths. Needed?
fromslump <- function(d, slump) {
  slump <- matrix(sort(slump), ncol=2, byrow=TRUE)
  slices <- slump[,2] - slump[,1]
  dorig <- d # original depths
  for(i in 1:nrow(slump)) {
    below <- which(d > min(slump[i,]))
  if(length(below) > 0)
    dorig[below] <- dorig[below] - slices[i]
  }
  return(dorig)
}

Try the rbacon package in your browser

Any scripts or data that you put into this service are public.

rbacon documentation built on March 29, 2022, 5:05 p.m.