# Workflow for obtaining age-depth model from Langlois Lake using silt as a proxy for 
# sedimentation rate. 

# Citation: Daniel G. Gavin, William T. Struble and Mark A. Fonstad. Holocene Lake Sediments 
# Reveal Alluvial Fan History with Links to Climate, Wildfire, and Earthquakes. 
# Journal of Geophysical Research: Earth Surfaces.


# This code should be applicable to other data sets.  Must specify several variables
# for other data sets: radiocarbon dates, slump depths (if any)

## Note that subsequent applications of this code will result in slightly varying results
## due to the randomizations underlying the age-depth model fitting 

source("silt_model_functions_R.txt")
require(Bchron)
require(svMisc)

# Input data

#### Read raw sediment color data: ls_composite_rgb_Lab.csv
# CIELab lightness values obtained from composite of core photos. Values are at 
# irregular intervals due to resolution of images.  Therefore we need to interpolate values. 
# This will interpolate to a constant depth interval (intv) between depths d.min and d.max.

ls<-read.csv("ls_composite_rgb_Lab.csv")
d.min <- round(min(ls$depth),1)
d.max <- round(max(ls$depth),1)
intv <- 0.05
f5<-floor(ls$depth*(1/intv))/(1/intv) #index for intv cm intervals
cie.l.5<-tapply(ls$CIE.L,f5,mean)  
prox <- data.frame(depths=as.numeric(rownames(cie.l.5)),gs=cie.l.5)
# ensure there are values every .05 cm: interpolate over gaps
prox2 <- approx(x=prox$depth,y=prox$gs,xout=seq(d.min,d.max,intv))

# remove slumps from depths, adjust depths to effective depths.  Add slumps back in 
# after computing age-depth model
# Slumps are 80-108 cm, and 370-375 cm
slump <- c(80,108.05,370,375.05)
slump <- (slump*(1/intv))-(d.min*(1/intv))  # as indices in the prox2 array

# sval=sediment lightness values at 0.05 cm resolution
sval <- prox2$y[c(1:slump[1],slump[2]:slump[3],slump[4]:length(prox2$x))]

min.s.in <- min(sval)  # minimum value of silt proxy
max.s.in <- max(sval)
# d.sval = depths of sval values, with no slumps
d.sval <- seq(d.min,(length(sval)/(1/intv))+(d.min-intv),intv) 

# number of randomizations of calibrated age pdfs to calculate age uncertainty
n.s<-1000

#### Calibrate radiocarbon dates, obtain 1000 draws from PDFs of calibrated ages
c14s <- BchronCalibrate(ages=c(999,1230,2997,3681,4000,5130),ageSds=c(23,50,23,23,60,60))
date.samples <- sampleAges(c14s,n.s)
# depth of 23.5 cm has Pb-210 date at 1867 AD +/- 3 (74 BP)
top.date <- rnorm(n=n.s,mean=83,sd=3)
date.samples <- cbind(top.date,date.samples)
depths<-c(24.1,66,126,180,330,383,522)
# remove slumps
ds<-c(24.1,66,(126-28),(180-28),(330-28),(383-28-5),(522-28-5))

# indexes in grayscale data where age controls exist
date.dep.in<-which(round(d.sval,2) %in% round(ds,1))
num.dates<-length(date.dep.in)

# test the fit to mean of 14C ages
date.age.in <- c(83,917,1149,3184,4032,4474,5869)  # date.samples[1,]
c.min <- optimize(silt.model,interval=c(20,100),sval=sval,min.s=min.s.in,date.dep=date.dep.in,date.age=date.age.in,maximum=FALSE)
# RMSE (root mean square error)
(c.min$objective/7)^.5   #187 years

#### Obtain a distribution of optimum crit.s based on 1000 resampling 14C dates.
#	Requires a long time (hours) to run 1000 iterations
c.min<-matrix(NA,n.s,2)
s.dates<-matrix(NA,n.s,num.dates)  #dates saved for each sim.

for(s in 1:n.s){	
	date.age.in <- date.samples[s,]
	out<-optimize(silt.model.c,interval=c(20,100),sval=sval,min.s=min.s.in,date.dep=date.dep.in,date.age=date.age.in,maximum=FALSE)
	c.min[s,1]<-out$minimum
	c.min[s,2]<-out$objective
	progress(100*s/n.s)
	} #next iteration

# Write file of s.crit and fit values for all iterations 
write.table(c.min,"cmin.csv",sep=",")
# median and range of $minimum (crit.s values)
quantile(c.min[,1],probs=c(0.05,0.5,0.95))
#       5%      50%      95% 
#  88.70763 91.50779 94.31825 

# median and range of ($objective/7)^.5 = RMSE = 223
quantile((c.min[,2]/7)^.5,probs=c(0.05,0.5,0.95))
#       5%      50%      95% 
# 155.4883 190.2176 230.8188 

#### Run spline fits on each resampling (using crit.s and resampled dates)
# ed=effective depths
min.s <- min.s.in
len <- length(sval)
chron <- matrix(NA,len,n.s)  #chron =  n.s age models for all depths in the input data

for(s in 1:n.s){
	crit.s<-c.min[s,1]
	theta.est<-array(NA,len)
	ed<-array(NA,len)
	for(i in 1:len){
		theta.est[i]<-1-min(1,((sval[i]-min.s)/(crit.s-min.s)))
		}
	for(i in 1:len){
		ed[i]<-sum(theta.est[1:i])*.05
		}

	# Fits ages to each actual depth, using the effective depths in the age spline.
	chron[,s] <- monspline1.c(date.samples[s,],ed[date.dep.in],ed)
	progress(100*s/n.s)
}

## In case you want to write out the (very large) table of all of the simulated fitted ages
## all depths * n.s simulations.  
write.table(chron,"chron.csv",sep=",")

#### Spaghetti plot of estimated depth vs age for all iterations (note: slow to calculate)
###  The array 'check' contains T/F of which iterations  have age control points that 
###  are not in stratigraphic order.  These iterations will be removed from age estimates.

## set up axes
plot(x=1:50,y=1:50,type="n",xlim=c(-50,6000),ylim=c(0,180),axes=F,xlab="",ylab="")
axis(1,at=seq(0,6000,1000))
axis(2,at=seq(0,180,20),las=1)
mtext("Year BP", side=1, las=1,line=2)
mtext("Effective depth (cm)",side=2, las=3,line=2.5)

check<-rep(TRUE,n.s)

# Calculate age-depth lines from the crit.s values. And, add lines to plot.
for(s in 1:n.s){
	crit.s<-c.min[s,1]
	theta.est<-array(NA,len)
	ed<-array(NA,len)
	for(i in 1:len){
		theta.est[i]<-1-min(1,((sval[i]-min.s)/(crit.s-min.s)))
		}
	for(i in 1:len){
		ed[i]<-sum(theta.est[1:i])*.05
		}

	#trim out age reversals in sampled ages
	check[s] <- TRUE
	for(a in 2:7){
		if(date.samples[s,a]<date.samples[s,(a-1)]){
			check[s] <- FALSE
			}
		}
	if(check[s]){
		lines(x=chron[,s],y=ed,lwd=1,col=rgb(0,0,0,0.05))
		}
	progress(100*s/n.s)
	}

# plot symbols for age control points
for(s in 1:n.s){  
	if(check[s]){
		crit.s<-c.min[s,1]
		theta.est<-array(NA,len)
		ed<-array(NA,len)
		for(i in 1:len){
			theta.est[i]<-1-min(1,((sval[i]-min.s)/(crit.s-min.s)))
			}
		for(i in 1:len){
			ed[i]<-sum(theta.est[1:i])*.05
			}
		points(date.samples[s,],ed[date.dep.in],pch=19,cex=0.4,col=rgb(.2,.2,1,0.05))
		}
	progress(100*s/n.s)
	}

# end plot commands

# Select median and the 90% prob intervals and write out. 
# Use only models with increasing ages ('check' array)
# these are the values used in Gavin et al. 

sum(check)  # Note: all 1000 iterations have ages increasing with depth

probs<-c(0.05,0.5,0.95)
p.out<-matrix(NA,len,length(probs))
for(i in 1:len){
	for(p in 1:length(probs)){
		p.out[i,p]<-quantile(chron[i,check],probs[p])
		}
	}

#### Calculate a chron using median crit.s (91.507) and fit to the median of the cal ages.
#### Note: this chron was not used in Gavin et al. 

for(i in 1:7){
print(median(date.samples[,i]))
}
#[1] 82.91128
#[1] 917
#[1] 1149
#[1] 3184.5
#[1] 4032.5
#[1] 4474.5
#[1] 5869

crit.s<-91.507
theta.est<-array(NA,len)
ed<-array(NA,len)
for(i in 1:len){
	theta.est[i]<-1-min(1,((sval[i]-min.s)/(crit.s-min.s)))
	}
for(i in 1:len){
	ed[i]<-sum(theta.est[1:i])*.05
	}
fit.mean <- monspline1.c(c(83,917,1149,3184,4032,4474,5869),ed[date.dep.in],ed)

# View the fitted line
plot(d.sval,fit.mean,type="l")
lines(d.sval,p.out[,1],type="l",col=2)
lines(d.sval,p.out[,3],type="l",col=2)

# Add in slumps to depths

depth2 <- seq(d.min,d.max,intv)
p05 <- array(NA,length(depth2))
p50 <- array(NA,length(depth2))
p95 <- array(NA,length(depth2))

p05[1:slump[1]] <- p.out[1:slump[1],1]
p05[(slump[2]+1):slump[3]] <- p.out[(slump[1]+1):(slump[1]+slump[3]-slump[2]),1]
p05[(slump[4]):length(depth2)] <- p.out[(slump[1]+slump[3]-slump[2]+2):length(p.out[,1]),1]

p50[1:slump[1]] <- p.out[1:slump[1],2]
p50[(slump[2]+1):slump[3]] <- p.out[(slump[1]+1):(slump[1]+slump[3]-slump[2]),2]
p50[(slump[4]):length(depth2)] <- p.out[(slump[1]+slump[3]-slump[2]+2):length(p.out[,1]),2]

p95[1:slump[1]] <- p.out[1:slump[1],3]
p95[(slump[2]+1):slump[3]] <- p.out[(slump[1]+1):(slump[1]+slump[3]-slump[2]),3]
p95[(slump[4]):length(depth2)] <- p.out[(slump[1]+slump[3]-slump[2]+2):length(p.out[,1]),3]

# Write out the age-depth file (values each 0.5 mm)
mm.ages<-data.frame(cm=depth2,p05=p05,p50=p50,p95=p95,gs=prox2$y)
write.table(mm.ages,"langlois_depth_age_model_output.csv",sep=",")

#### Bin theta to decadal values, corresponding to mm silt accumulation per decade
# To obtain a 90% confidence interval for silt accumulation by using a high and low age model
# and a high and low value of s.crit.  
# Low crit.s corresponds to high sedimentation rates, so using p.05 age model.

quantile(c.min[,1],probs=c(0.05,.5,.95))

# percentage silt (theta) for each crit.s
perc.silt<-matrix(NA,len,3)
for(i in 1:len){
	perc.silt[i,1]<-min(1,((sval[i]-min.s)/(88.7076-min.s)))
	perc.silt[i,2]<-min(1,((sval[i]-min.s)/(91.5078-min.s)))
	perc.silt[i,3]<-min(1,((sval[i]-min.s)/(94.3183-min.s)))
	}

# To integrate precisely to decadal bins, first interpolate perc.silt and ages to
# very fine intervals: 0.001 cm
perc.silt.hr.05 <- approx(x=d.sval,y=perc.silt[,1],xout=seq(min(d.sval),max(d.sval),0.001))
perc.silt.hr.50 <- approx(x=d.sval,y=perc.silt[,2],xout=seq(min(d.sval),max(d.sval),0.001))
perc.silt.hr.95 <- approx(x=d.sval,y=perc.silt[,3],xout=seq(min(d.sval),max(d.sval),0.001))

# Median age model. Interpolate ages to the same fine intervals
chron.hr <- approx(x=d.sval,y=p.out[,2],xout=seq(min(d.sval),max(d.sval),0.001))
perc.silt.decadal <- tapply(X=perc.silt.hr.50$y, INDEX=floor(chron.hr$y/10), FUN=mean)
thickness.decadal <- tapply(X=perc.silt.hr.50$y, INDEX=floor(chron.hr$y/10), FUN=length)
silt.decadal.p50 <- perc.silt.decadal*thickness.decadal/100
out.p50 <- cbind(10*as.numeric(rownames(perc.silt.decadal)),silt.decadal.p50)
colnames(out.p50) <- c("yrBP","silt.mm")
write.table(out.p50,file="Fig5a_silt_decadal_p50.csv",sep=",")

# for lo and hi age models: Not used.
chron.hr <- approx(x=d.sval,y=p.out[,1],xout=seq(min(d.sval),max(d.sval),0.001))
perc.silt.decadal <- tapply(X=perc.silt.hr.05$y, INDEX=floor(chron.hr$y/10), FUN=mean)
thickness.decadal <- tapply(X=perc.silt.hr.05$y, INDEX=floor(chron.hr$y/10), FUN=length)
silt.decadal.p05 <- perc.silt.decadal*thickness.decadal/100
write.table(cbind(10*as.numeric(rownames(perc.silt.decadal)),silt.decadal.p05),file="silt_decadal_p05.csv",sep=",")

chron.hr <- approx(x=d.sval,y=p.out[,3],xout=seq(min(d.sval),max(d.sval),0.001))
perc.silt.decadal <- tapply(X=perc.silt.hr.95$y, INDEX=floor(chron.hr$y/10), FUN=mean)
thickness.decadal <- tapply(X=perc.silt.hr.95$y, INDEX=floor(chron.hr$y/10), FUN=length)
silt.decadal.p95 <- perc.silt.decadal*thickness.decadal/100
write.table(cbind(10*as.numeric(rownames(perc.silt.decadal)),silt.decadal.p95),file="silt_decadal_p95.csv",sep=",")