###########################################################################
#FUNCTION FOR ONE-PHASE EXPONENIAL MODEL OF LOSS FROM THE MIXED LAYER
#WITH ABRUPT INCREAE AND TERMINATION IN PRODUCTION 
###########################################################################
OneStateConstantTwiceTruncatedExp=function(theta){
lambda1=exp(theta)
N=length(ages)
first01=0
first02=(lambda1*(exp(-lambda1*TIME)-exp(-lambda1*TIMEMAX)))/lambda1
first=(1/(first01+first02))
second=lambda1*exp(-lambda1*ages)
prsecond=sum(log(second)) 
loglik=N*log(first)+prsecond
-(loglik)
}

###########################################################################
#FUNCTION FOR TWO-PHASE EXPONENIAL MODEL OF LOSS FROM THE MIXED LAYER
#WITH ABRUPT INCREAE AND TERMINATION IN PRODUCTION 
###########################################################################
OneStateConstantTwiceTruncatedExpMix=function(theta){
lambda1=exp(theta[1])
lambda2=exp(theta[2])
tau=exp(theta[3])
N=length(ages)
first01=(tau*(exp(-lambda2*TIME)-exp(-lambda2*TIMEMAX)))/lambda2
first02=((lambda1-lambda2)*(exp(-(tau+lambda1)*TIME)-exp(-(tau+lambda1)*TIMEMAX)))/(tau+lambda1)
#first01=ifelse(first01>0,first01,10e-50)
#first02=ifelse(first02>0,first02,10e-50)
first=(1/(first01+first02))
second=(tau*exp(-lambda2*ages) + (lambda1-lambda2)*exp(-(tau+lambda1)*ages))
#second=ifelse(second>0,second,10e-50)
prsecond=sum(log(second)) 
loglik=N*log(first)+prsecond
if (lambda2>=(lambda1)) {loglik=-100000}
-(loglik)
}


#####################################################################################
#REQUIRED PACKAGES
#####################################################################################
require(msm)
require(TruncatedDistributions)
require(gridBase)
require(modeest)
require(grid)

#####################################################################################
#OTHER FUNCTIONS FOR BOOTSTRAP AND HISTOGRAM
bootcimode<-function(temp) {
x.resample<-numeric(); 
xmode.resample<-numeric(); 
for (i in 1:1000) {
for (j in 1:length(temp)) {x.resample[j]<-sample(temp,1,replace=T)}
xmode.resample[i]<-mlv(x.resample, method = "hsm", bw=0.1)$M
}
out=list(resamplemode=xmode.resample)
out
}

horiz.hist <- function(Data, breaks="Sturges", col="transparent", las=1, 
ylim=range(HBreaks), labelat=pretty(ylim), labels=labelat, border=par("fg"), ... )
  {a <- hist(Data, plot=FALSE, breaks=breaks)
  HBreaks <- a$breaks
  HBreak1 <- a$breaks[1]
  hpos <<- function(Pos) (Pos-HBreak1)*(length(HBreaks)-1)/ diff(range(HBreaks))
  barplot(a$counts, space=0, horiz=T, ylim=hpos(ylim), col=col, border=border,...)      
  axis(2, at=hpos(labelat), labels=labels, las=las, ...) 
  print("use hpos() to address y-coordinates") 
}

######################################################################################################################################################
#VECTOR WITH POSTMORTEM AGES FROM AAR CALIBRATED BY RADIOCARBON
ageLaqueusSCB=c(108,164,4623,182,380,269,211,471,405,269,4708,168,211,943,4068,339,2969,4832,217,1465,136,
250,128,1417,147,228,461,155,302,159,222,118,1015,222,155,4014,1846,1568,222,3385,256,
1214,2759,723,1300,461,738,461,6179,2236,5800,1417,233,979,461,797,177,405,143,263,239,
155,168,196,317,115,121,132,186,738,355,105,177,115,132,244,168,201,206,211,1072,
481,288,159,217,182,159,710,355,324,461,943,275,191,233,217,196,220,155,147,452,
186,172,2488,1974,1137,228,1227,196,132,397,151,524,414,228,159,347,2925,324,7089,7622,
2137,442,172,581,1111,1417,617,1489,164,1214,2341,1192,683,828,546,1595,405,1675,630,2005,
3901,5874,182,2718,339,442,683,4122,6416,5800,3688,388,324,1489,452,535,3846,843,1621,2527,
1758,960,1130,1893,1788,1846,2203,1072,1940,2236,2203,8423,875,1702,8482,943,8472,2137,1130,1541,
1393,960,1111,1256,1465,1072,3194,2104,9431)

#SITE ASSIGNMENTS OF SHELLS: 4134 - WESTERN PALOS VERDES SHELF (2003),
#OC - WESTERN PALOS VERDES SHELF (2008), 10C - EASTERN PALOS VERDES SHELF (2008 AND 2009), 24205 - SANTA MONICA BAY (1975)
site=c("4134","4134","4134","4134","4134","4134","4134","4134","4134","4134","4134","4134","4134", 
"4134","4134","4134","4134","4134","4134","4134","4134","4134","4134","4134","4134","4134", 
"4134","4134","4134","4134","4134","4134","4134","4134","4134","4134","4134","4134","4134", 
"10C", "10C", "10C", "10C", "10C", "10C", "10C", "10C", "10C", "10C", "10C", "10C", "10C",
"10C", "10C", "10C", "10C", "0C","0C","0C","0C","0C","0C","0C","0C","0C", 
"0C","0C","0C","0C","0C","0C","0C","0C","0C","0C","0C","0C","0C", 
"0C","0C","0C","0C","0C","0C","0C","0C","0C","0C","0C","0C","0C",   
"0C","0C","0C","0C","0C","0C","0C","0C","0C","0C","0C","0C","0C",   
"0C","0C","0C","0C","0C","0C","0C","0C","0C","0C","0C","0C","10C",
"10C","10C","10C","10C","10C","10C","10C","10C","10C","10C","10C","10C","10C", 
"10C","10C","10C","10C","10C","10C","10C","10C","10C","10C","10C","10C","10C",  
"10C","10C","10C","10C","10C","10C","10C","10C","10C","10C","10C","10C","10C",  
"10C","10C","10C","24205","24205","24205","24205","24205","24205","24205","24205","24205","24205",
"24205","24205","24205","24205","24205","24205","24205","24205","24205","24205","24205","24205","24205",
"24205","24205","24205","24205","24205","24205","24205","24205")

#logged standard deviation from calibration (for lognormal ucertainty) 
ln.d=-2.291966

######################################################################################################################################################
#REGIONAL LAQUEUS HISTOGRAM FIGURE
######################################################################################################################################################
BREAKS=100
mode.est=quantile(bootcimode(ageLaqueusSCB)$resamplemode, c(0.025,0.5,0.975))
mlv(ageLaqueusSCB, method="hsm", bw=0.1)

par(mfrow=c(2,2))
par(mar=c(10,4,2,1))
out1=hist(ageLaqueusSCB, breaks=c(0,seq(13,10000,by=BREAKS)), col="gray", ylim=c(0,50), main="", freq=T,
xlab="Shell age (years before 2013)", ylab="Number of specimens", cex.axis=1)
mode.est
range(ageLaqueusSCB)
INSET=TRUE
if (INSET==TRUE)  {
vp <- baseViewports()
pushViewport(vp$inner,vp$figure,vp$plot)
pushViewport(viewport(x=0.25,y=0.99,width=0.7,height=0.7,just=c("left","top")))
par(plt=gridPLT(),new=T) 
grid.rect(gp=gpar(lwd=1,col="black"))
breaks2=c(0,seq(13,1000,by=25)); smallage<-ageLaqueusSCB[ageLaqueusSCB<max(breaks2)]; 
hdata<-hist(smallage,breaks=breaks2,plot=F)
yaxis=c(0,max(hdata$counts)+5)
hdata<-hist(smallage,breaks=breaks2,xlab=NULL,ylab=NULL, main=NULL, freq=T,col="gray", ylim=yaxis,cex.axis=1, cex.lab=1, cex.axis=1)
times=hdata$mids
abline(v=mode.est[2], lty=2)
popViewport(4)
par(mar=c(4,4,2,1))
}

IQR(ageLaqueusSCB)
######################################################################################################################################################
#CONFIDENCE INTERVALS ON MINIMUM AGE, MEDIAN AGE, MODE, ACCOUNTING FOR CALIBRATION ERRORS
######################################################################################################################################################
REPS=1000
Full.Laqueus.bootout.250=array(NA, dim=c(REPS,40))
Full.Laqueus.bootout.50=array(NA, dim=c(REPS,200))
Full.Laqueus.bootout.10=array(NA, dim=c(REPS,1000))
Full.Laqueus.bootout.25=array(NA, dim=c(REPS,400))

boot.mode=numeric()
boot.min=numeric()
boot.max=numeric()
boot.IQR=numeric()
boot.range95=numeric()
for (j in 1:REPS) {
temp=ageLaqueusSCB
sample.age=sample(temp,length(ageLaqueusSCB),replace=T)
bootages=rtlnorm(length(ageLaqueusSCB),log(temp),sdlog=sqrt(exp(ln.d)),a=0,b=10000)
out1=hist(bootages[bootages<10000], breaks=seq(0,10000,by=250)) 
out2=hist(bootages[bootages<10000], breaks=seq(0,10000,by=50)) 
out3=hist(bootages[bootages<10000], breaks=seq(0,10000,by=10)) 
out4=hist(bootages[bootages<10000], breaks=seq(0,10000,by=25)) 
Full.Laqueus.bootout.250[j,]=out1$counts
Full.Laqueus.bootout.50[j,]=out2$counts
Full.Laqueus.bootout.10[j,]=out3$counts
Full.Laqueus.bootout.25[j,]=out4$counts

boot.mode[j]=unlist(mlv(bootages, method = "hsm", bw=0.1)[1])
boot.min[j]=min(bootages)
boot.max[j]=max(bootages)
boot.min[j]=quantile(bootages, 0.025)
boot.max[j]=quantile(bootages, 0.975)
boot.IQR[j]=IQR(bootages)
boot.range95[j]=diff(quantile(bootages, c(0.025,0.975)))
}

quantile(boot.IQR, c(0.025,0.5,0.975))
2013-quantile(boot.mode, c(0.025,0.5,0.975))
2013-quantile(boot.min, c(0.025,0.5,0.975))
2013-quantile(boot.max, c(0.025,0.5,0.975))

#######################################################################
#ACCOUNTING FOR CALIBRATION ERRORS - TIME AVERAGING CORRECTION 
#######################################################################
mean.IQR.shell=numeric()
mean.range95.shell=numeric()
mean.025.shell=numeric()
mean.975.shell=numeric()
for (i in 1:length(ageLaqueusSCB)) {
temp=ageLaqueusSCB[i]
sample.age=rtlnorm(1000,log(temp),sdlog=sqrt(exp(ln.d)),a=0,b=10000)
mean.IQR.shell[i]=IQR(sample.age)
mean.range95.shell[i]=diff(quantile(sample.age, c(0.025,0.975)))
mean.025.shell[i]=quantile(sample.age, c(0.025))
mean.975.shell[i]=quantile(sample.age, c(0.975))
}

mean(mean.IQR.shell)
mean(mean.range95.shell)
IQR(ageLaqueusSCB)
diff(quantile(ageLaqueusSCB, c(0.025,0.975)))
IQR(ageLaqueusSCB)-mean(mean.IQR.shell)
diff(quantile(ageLaqueusSCB, c(0.025,0.975)))-mean(mean.range95.shell)

##############################################################################################
#NUMBER OF INDIVIDUALS IN 250y, 50y, 10y, 25y COHORTS ON THE BASIS OF CALIBRATION RESAMPLING
##############################################################################################
Laqueus.250.means=apply(Full.Laqueus.bootout.250,MARGIN=2,mean)
Laqueus.250.025=apply(Full.Laqueus.bootout.250,MARGIN=2,quantile, 0.025)
Laqueus.250.975=apply(Full.Laqueus.bootout.250,MARGIN=2,quantile, 0.975)
names(Laqueus.250.means)=seq(125,9875,by=250)
scale.250=seq(125,9875,by=250)

Laqueus.50.means=apply(Full.Laqueus.bootout.50,MARGIN=2,mean)
Laqueus.50.025=apply(Full.Laqueus.bootout.50,MARGIN=2,quantile, 0.025)
Laqueus.50.975=apply(Full.Laqueus.bootout.50,MARGIN=2,quantile, 0.975)
names(Laqueus.50.means)=seq(25,9975,by=50)
scale.50=seq(25,9975,by=50)

Laqueus.10.means=apply(Full.Laqueus.bootout.10,MARGIN=2,mean)
Laqueus.10.025=apply(Full.Laqueus.bootout.10,MARGIN=2,quantile, 0.025)
Laqueus.10.975=apply(Full.Laqueus.bootout.10,MARGIN=2,quantile, 0.975)
names(Laqueus.10.means)=seq(5,9995,by=10)
scale.10=seq(5,9995,by=10)

Laqueus.25.means=apply(Full.Laqueus.bootout.25,MARGIN=2,mean)
Laqueus.25.025=apply(Full.Laqueus.bootout.25,MARGIN=2,quantile, 0.025)
Laqueus.25.975=apply(Full.Laqueus.bootout.25,MARGIN=2,quantile, 0.975)
names(Laqueus.25.means)=seq(12.5,9987.5,by=25)
scale.25=seq(12.5,9987.5,by=25)
Laqueus.25.mode=apply(Full.Laqueus.bootout.25,MARGIN=1,which.max)

#################################################
#MODE
#################################################
par(mfrow=c(2,3))
hist(2013-boot.mode[boot.mode<1700], breaks=rev(c(seq(1600,2000,by=10),2013)), col="gray", freq=T,
xlab="Date of mode of AFD (A.D.)", main="Modes expected under AAR age uncertainty", cex.main=0.9, xlim=c(1600,2020))
abline(v=2013-quantile(boot.mode,c(0.025,0.5,0.975)),lty=c(2,1,2),lwd=c(1,2,1))
2013-quantile(boot.mode,c(0.025,0.5,0.975))

#################################################
#MISSING PRODUCTION - TRUNCATED RECENT PRODUCTON#
#################################################
age=ageLaqueusSCB
ages=ageLaqueusSCB
TIMEMAX=10000
TRUNCATEDTIME=min(ages)
TIME=TRUNCATEDTIME
TIMEMIN=TIME


out1=optim(par=log(0.0001), OneStateConstantTwiceTruncatedExp)
fit1phaseUNDERtruncation=exp(out1$par)
lambda=exp(out1$par)
lambda
LIK=out1$value
k=3
AIC=k*1-(-2*LIK)
AICc=AIC+((2*k*(k+1))/(length(age)-k-1))
AICc

start1=log(c(0.002,0.0001,0.00001))
out=optim(par=start1, OneStateConstantTwiceTruncatedExpMix)
fit2phasetauUNDERtruncation=exp(out$par[3])
fit2phaselambda2UNDERtruncation=exp(out$par[2])
fit2phaselambda1UNDERtruncation=exp(out$par[1])
out$value
LIK=out$value
k=5
AIC=k*1-(-2*LIK)
AICc=AIC+((2*k*(k+1))/(length(age)-k-1))
AICc


lambda1=fit2phaselambda1UNDERtruncation
lambda2=fit2phaselambda2UNDERtruncation
tau=fit2phasetauUNDERtruncation
randomtwophaserate1<-lambda2
randomtwophaserate2<-tau+lambda1
rate1<-lambda2
rate2<-tau+lambda1
randomtwophasealpha<-tau/(tau+lambda1-lambda2)
alpha<-tau/(tau+lambda1-lambda2)
randomtwophasebeta<-((1/rate1)*alpha)/((1/rate1)*alpha+(1/rate2)*(1-alpha))

lambda1
tau
lambda2
VECTOR=seq(5,9995,by=10)

#Tomasovych et al. 2016 Paleobiology equation 14
term1=(tau*(exp(-lambda2*TIMEMIN)-exp(-lambda2*TIMEMAX)))/lambda2
term2=(lambda1-lambda2)*(exp(-(tau+lambda1)*TIMEMIN)-exp(-(tau+lambda1)*TIMEMAX))
term3=tau+lambda1
CONSTANT=1/(term1+term2/term3)
temp1=tau*exp(-lambda2*VECTOR)+(lambda1-lambda2)*exp(-(tau+lambda1)*VECTOR)
temp2=randomtwophasealpha*exp(-randomtwophaserate1*seq(0,9990,by=10))+(1-randomtwophasealpha)*exp(-randomtwophaserate2*seq(0,9990,by=10))
temp.lambda1=lambda1
temp.lambda2=lambda2
temp.tau=tau
temp.alpha=temp.tau/(temp.tau+temp.lambda1-temp.lambda2)

VECTOR.10=seq(5,9995,by=10)
invsurv.10=(temp.alpha*exp(-temp.lambda2*VECTOR.10)+(1-temp.alpha)*exp(-(temp.tau+temp.lambda1)*VECTOR.10))
VECTOR.25=seq(12.5,9987.5,by=25)
invsurv.25=(temp.alpha*exp(-temp.lambda2*VECTOR.25)+(1-temp.alpha)*exp(-(temp.tau+temp.lambda1)*VECTOR.25))
VECTOR.50=seq(25,9975,by=50)
invsurv.50=(temp.alpha*exp(-temp.lambda2*VECTOR.50)+(1-temp.alpha)*exp(-(temp.tau+temp.lambda1)*VECTOR.50))



#################################################
#LAST PEAK OF PRODUCTION
#################################################
out=t(Full.Laqueus.bootout.10[,1:24])/invsurv.10[1:24]
invsurv.mode.10=apply(out,MARGIN=2,which.max)
hist(2013-invsurv.mode.10*10, breaks=seq(1350,2013,by=10), xlim=c(1600,2020), col="gray", cex.main=0.9,
main="Date of last peak in production",xlab="Date of last peak in production (A.D.)")
abline(v=2013-quantile(invsurv.mode.10*10,c(0.025)),lty=c(2),lwd=c(1,2,1))
2013-quantile(invsurv.mode.10*10,c(0.025,0.5,0.975))
#which.max(Laqueus.10.means[1:40]/invsurv.10[1:40])

########################################################################################################
#RECONSTRUCTING THE TRAJECTORY, ASSUMING STEADY-STATE PRODUCTION PRIOR TO THE DECLINE
########################################################################################################
plot(VECTOR.25[1:40],Laqueus.25.means[1:40],type="l", ylim=c(0,80))
lines(VECTOR.25[1:40],Laqueus.25.means[1:40]/invsurv.25[1:40],lwd=2)
lines(VECTOR.25[1:40],Laqueus.25.025[1:40]/invsurv.25[1:40],lty=2)
lines(VECTOR.25[1:40],Laqueus.25.975[1:40]/invsurv.25[1:40],lty=2)


##############################################################################################################################
#THE PULL OF THE RECENT EFFECT ON RANDOM DATA
##############################################################################################################################
par(mfrow=c(2,2))
RATE=0.005
time=seq(1,1000,by=1)
input=rep(100,1000)
input=dnorm(1:1000,mean=500,sd=200)
plot(time, input, ylab="Density", xlab="Time", type="l", col="blue", lwd=2, ylim=c(0,0.0025))
preserved=input*exp(-RATE*time)
lines(time, preserved, col="red", lwd=2)
legend(x="topright",col=c("blue","red"), lwd=2, legend=c("Production","Preserved"))

##############################################################################################
#FIGURE HISTOGRAM VS RECONSTRUCTED
##############################################################################################
par(mfrow=c(1,3))
AD.ageLaqueusSCB=2013-ageLaqueusSCB
par(mar=c(4,4,2,1))
horiz.hist(AD.ageLaqueusSCB[AD.ageLaqueusSCB>1000], breaks=seq(1000,2020,by=25), col="gray", ylab="", xlab="Number of specimens", xlim=c(0,15))
mode.est
range(ageLaqueusSCB)

plot(Laqueus.25.means[scale.25<1000], -scale.25[scale.25<1000], xlab="Abundance (25 year cohorts)", type="n", 
frame.plot=FALSE,ylab="Time (years)",pch=16, log="", xlim=c(0,40), ylim=c(-1000,0))
lines(Laqueus.25.means[scale.25<1000], -scale.25[scale.25<1000], lwd=3)
lines(Laqueus.25.025[scale.25<1000],-scale.25[scale.25<1000], lty=2)
lines(Laqueus.25.975[scale.25<1000],-scale.25[scale.25<1000], lty=2)
scale.25[which.max(Laqueus.25.means)]

VECTOR=seq(12.5,9987.5,by=25)
invsurv=Laqueus.25.means/exp(-0.002*seq(12.5,9987.5,by=25))
invsurv=Laqueus.25.means/(temp.alpha*exp(-temp.lambda2*VECTOR)+(1-temp.alpha)*exp(-(temp.tau+temp.lambda1)*VECTOR))
seg.UCI=Laqueus.25.025/(temp.alpha*exp(-temp.lambda2*VECTOR)+(1-temp.alpha)*exp(-(temp.tau+temp.lambda1)*VECTOR))
seg.LCI=Laqueus.25.975/(temp.alpha*exp(-temp.lambda2*VECTOR)+(1-temp.alpha)*exp(-(temp.tau+temp.lambda1)*VECTOR))
lines(invsurv[scale.25<1000], -scale.25[scale.25<1000], col="gray", lwd=3)
lines(seg.LCI[scale.25<1000],-scale.25[scale.25<1000], lty=2, col="gray")
lines(seg.UCI[scale.25<1000],-scale.25[scale.25<1000], lty=2, col="gray")
legend(x="topright",lwd=c(3,3), col=c("black","gray"),legend=c("Preserved","Reconstructed"), bty="n")


##############################################################################################
#FIGURE HISTOGRAM VS RECONSTRUCTED IN ONE FIGURE AND LIFESPAN TO STANDING DENSITY
##############################################################################################
AD.ageLaqueusSCB=2013-ageLaqueusSCB
par(mfrow=c(1,4))
par(mar=c(25,4,4,0))
out=horiz.hist(AD.ageLaqueusSCB[AD.ageLaqueusSCB>1400], breaks=seq(1400,2020,by=25), col="gray", ylab="", xlab="Number of specimens", xlim=c(0,40), ylim=c(1400,2020))
out=hist(AD.ageLaqueusSCB[AD.ageLaqueusSCB>1400], breaks=seq(1400,2020,by=25), plot=F)

par(new=T)
plot(Laqueus.25.means[scale.25<600], 2013-scale.25[scale.25<600], xlab="", type="l", ylab="", frame.plot=FALSE, xlim=c(0,40), ylim=c(1400,2020), yaxt="n", xaxt="n")
lines(Laqueus.25.means[scale.25<1000], 2013-scale.25[scale.25<1000], lwd=3, col="gray41")
lines(Laqueus.25.025[scale.25<1000],2013-scale.25[scale.25<1000], lty=2, col="gray41")
lines(Laqueus.25.975[scale.25<1000],2013-scale.25[scale.25<1000], lty=2, col="gray41")
scale.25[which.max(Laqueus.25.means)]

VECTOR=seq(12.5,9987.5,by=25)
invsurv=Laqueus.25.means/exp(-0.002*seq(12.5,9987.5,by=25))
invsurv=Laqueus.25.means/(temp.alpha*exp(-temp.lambda2*VECTOR)+(1-temp.alpha)*exp(-(temp.tau+temp.lambda1)*VECTOR))
seg.UCI=Laqueus.25.025/(temp.alpha*exp(-temp.lambda2*VECTOR)+(1-temp.alpha)*exp(-(temp.tau+temp.lambda1)*VECTOR))
seg.LCI=Laqueus.25.975/(temp.alpha*exp(-temp.lambda2*VECTOR)+(1-temp.alpha)*exp(-(temp.tau+temp.lambda1)*VECTOR))
lines(invsurv[scale.25<1000], 2013-scale.25[scale.25<1000], col="black", lwd=3)
lines(seg.LCI[scale.25<1000],2013-scale.25[scale.25<1000], lty=2, col="black")
lines(seg.UCI[scale.25<1000],2013-scale.25[scale.25<1000], lty=2, col="black")
legend(x="topright",lwd=c(3,3), col=c("black","gray"),legend=c("Preserved","Reconstructed"), bty="n")
#LIFESPAN
#total area = 0.5 m2, i.e., five grabs, therefore, multiply by 2 to get individuals/m2
#lifespan of Laqueus = 12 years, abundance is per 25 year cohort, therefore 25/12 = 2.08, divide by 2.08 to get abundance per year
#finally, if all pedicle and brachial valves would be counted, the number should be divided by two because each individual consist of two valves
#we use 1.3 rather than 2 because in some samples, just one type of valves was dated
axis(3,at=c(0,20,40),labels= round(c(0,20,40)*2/2.08/1.3))
(invsurv[scale.25<600]*2)/2.08/1.3