## R functions for computing age-depth model following Colombaroli et al. (2018)
## Example usage of these functions is in code in file silt-model-workflow-R.txt

## Colombaroli, D., Gavin, D.G., Morey, A.E. & Thorndycraft, V.R. 
## High resolution lake sediment record reveals self-organized criticality in erosion 
## processes regulated by internal feedbacks. EARTH SURFACE PROCESSES AND LANDFORMS 43, 
## 2181–2192 (2018).

require(compiler)

## Function silt.model  This function is called by the "optimize" function, to compute 
## optimal value of crit.s. See code in silt-model-workflow-R.txt

# crit.s	value to be optimized to minimize residuals of straight-line fit between age control points
#			and effective depths
# sval		vector of proxy values for silt content of sediment. Could be density, magnetic 
#			susceptibility, or grayscale color. Ordered from core top to base at equal depth intervals.
# min.s		single value of minimum sval corresponding to 0% silt content
# date.dep	vector of depths of age control points, as indices of the sval vector.
# date.age	vector of ages of age control points

silt.model <- function(crit.s,sval,min.s,date.dep,date.age){
	len<-length(sval)
	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])
		}
	fit<-lm(date.age~ed[date.dep])
	return(sum(fit$residuals^2))
	}

silt.model.c<-cmpfun(silt.model)


#### Monotonic spline functions. Unpublished algorithm by CJC Kruger.
# For given sequence of control points (depths (depths) and ages (ag)), 
# provide the ages at a sequence of depths (depthseq) using constrained cubic spline
# interpolation.

dxx <- function(x1,x0){
	dx <- x1-x0
	if(dx == 0) dx <- 10^30
	return(dx)
	}

# monspline2 = fit one depth
monspline2<-function(x,ag,depths) 
	{
    gxx<-array(NA,2); ggxx<-array(NA,2)
	Nmax <- length(depths)
	#(1a) Find LineNumber or segment. Linear extrapolate if outside range.
	Num <- 1
	if(x<depths[1] || x>depths[Nmax]){
		if(x<depths[1]) Num<-2 else Num <- Nmax
		B <- (ag[Num]-ag[(Num-1)])/dxx(depths[Num],depths[(Num-1)])
      	A <- ag[Num]-B*depths[Num]
      	return(A + B*x)
      	
      	#(1b) Find LineNumber or segment. Linear extrapolate if outside range.
      	} else { 
      	for(z in 2:Nmax){
      		if(x <= depths[z]){
      			Num <- z
      			break
      			}
      		}
      	}
    #(2) Calc first derivative (slope) for intermediate points
    for(j in 1:2){  #Two points around line
		z<-Num-2+j
		if(z==1 || z==Nmax) gxx[j] <- 10^30 else #Set very large slope at ends
		#Only check for 0 dy. dx assumed NEVER equals 0 !
		if(ag[(z+1)]-ag[z]==0 & ag[z]-ag[(z-1)]==0) gxx[j]<-0 else 
		#Pos PLUS neg slope is 0. Prevent div by zero.
		if(((depths[z+1] - depths[z])) / (ag[(z+1)] - ag[z]) + (depths[z] -depths[(z-1)]) / (ag[z]-ag[(z-1)]) == 0) gxx[j]<-0 else 
		#Pos AND neg slope, assume slope = 0 to prevent overshoot
		if((ag[(z+1)] - ag[z]) * (ag[z] - ag[(z-1)]) < 0) gxx[j]<-0 else 
		#Calculate an average slope for point based on connecting lines
		gxx[j] <- 2 / (dxx(depths[(z+1)], depths[z]) / (ag[(z+1)] - ag[z]) + dxx(depths[z], depths[(z-1)]) / (ag[z]-ag[(z-1)])) 
		} #next j

	#(3) Reset first derivative (slope) at first and last point
	if(Num==2){ #First point has 0 2nd derivative
		gxx[1] = 3 / 2 * (ag[Num] - ag[(Num-1)]) / dxx(depths[Num], depths[(Num-1)]) - gxx[2] / 2
		}
	if(Num==Nmax){ #Last point has 0 2nd derivative
		gxx[2] <- 3 / 2 * (ag[Num] - ag[(Num-1)]) / dxx(depths[Num], depths[(Num-1)]) - gxx[1] / 2
		}
	#(4) Calc second derivative at points
	ggxx[1] <- -2 * (gxx[2] + 2 * gxx[1]) / dxx(depths[Num], depths[(Num-1)]) + 6 * (ag[Num] - ag[(Num-1)]) / dxx(depths[Num], depths[(Num-1)]) ^ 2
	ggxx[2] <- 2 * (2 * gxx[2] + gxx[1]) / dxx(depths[Num], depths[(Num-1)]) - 6 * (ag[Num] - ag[(Num-1)]) / dxx(depths[Num], depths[(Num-1)]) ^ 2

	#(5) Calc constants for cubic
	D <- 1 / 6 * (ggxx[2] - ggxx[1]) / dxx(depths[Num], depths[(Num-1)])
	C <- 1 / 2 * (depths[Num] * ggxx[1] - depths[(Num-1)] * ggxx[2]) / dxx(depths[Num], depths[(Num-1)])
	B <- (ag[Num] - ag[(Num-1)] - C * (depths[Num]^2 - depths[(Num-1)]^2) - D * (depths[Num]^3 - depths[(Num-1)]^3)) / dxx(depths[Num], depths[(Num-1)])
	A <- ag[(Num-1)] - B * depths[(Num-1)] - C * depths[(Num-1)] ^ 2 - D * depths[(Num-1)]^3      		
	fitsp<-A + B * x + C * x ^ 2 + D * x ^ 3
	return(fitsp)
	} 
monspline2.c<-cmpfun(monspline2)


monspline1<-function(ag,depths,depthseq){
	#loop through all depths where an age is needed to be calculated
	temp<-array(NA,length(depthseq))
	for(x in 1:length(depthseq)){ 
    	temp[x]<-monspline2.c(depthseq[x],ag,depths)
    	}
    return(temp)
    }
monspline1.c<-cmpfun(monspline1)