Nothing
#' @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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.