# this script is to be used with the software 'R' (R Development Core Team, 2010. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org.)
# citation for this script is 'Stanford, J. D., Hemingway, R., Rohling, E. J., Challenor, P. G., Medina-Elizalde, M., Lester, A. J., In press. Sea-level probabilities for the last deglaciation: a statistical analysis of far-field records, Global and Planetary Change.'
  
# load the data.  Replace ... with the pathway of the directory in which the .csv files have been saved 

barbados_new <- read.csv(".../Barbados.csv")
sunda2m_new <- read.csv(".../Sunda_2m.csv")
sunda7m_new <- read.csv(".../Sunda_7m.csv")
bonaparte_bm3_new <- read.csv(".../Bonaparte_bm3.csv")
bonaparte_mm5_new <- read.csv(".../Bonaparte_mm5.csv")
tahiti_new <- read.csv(".../Tahiti_2010_1996.csv")
huon_new <- read.csv(".../Huon.csv")
florida_cor_new <- read.csv(".../Florida_corals.csv")
tosP21_new <- read.csv(".../Tos_peat_21.csv")
tosP50_new <- read.csv(".../Tos_peat_50.csv")
tosP66_new <- read.csv(".../Tos_peat_66.csv")


# n_barbados and n_sunda are the number of data points in each data set

n_barbados <- dim(barbados_new)[1]
n_sunda2m <- dim(sunda2m_new)[1]
n_sunda7m <- dim(sunda7m_new)[1]
n_bonaparte_bm3 <- dim(bonaparte_bm3_new)[1]
n_bonaparte_mm5 <- dim(bonaparte_mm5_new)[1]
n_tahiti <- dim(tahiti_new)[1]
n_huon <-dim(huon_new)[1]
n_florida_cor <-dim(florida_cor_new)[1]
n_tosP21 <-dim(tosP21_new)[1]
n_tosP50 <-dim(tosP50_new)[1]
n_tosP66 <-dim(tosP66_new)[1]

# number of simulations
n_simul <- 1000

# time values
# startt is the start time
# lent is the length of the period
# n_points is the number of points
startt <- 0
lent <- 22
n_points <- 220

# setting up the smoothing spline

smooth_x <- c(1:n_points)/n_points
smooth_x <- lent*smooth_x+startt

# setting the depth uncertainty for the lognormal depth parameters

coralsd <- 6
coralvar <- coralsd^2
bonaparte_mm5var <- 5^2

# setup the lognormal parameters
barbados_new$logmean <- log(-barbados_new$sealevel)-0.5*log(1+coralvar/(barbados_new$sealevel^2))
barbados_new$logsd <- sqrt(log(1+coralvar/(barbados_new$sealevel^2)))
bonaparte_mm5_new$logmean <- log(-bonaparte_mm5_new$sealevel)-0.5*log(1+bonaparte_mm5var/(bonaparte_mm5_new$sealevel^2))
bonaparte_mm5_new$logsd <- sqrt(log(1+bonaparte_mm5var/(bonaparte_mm5_new$sealevel^2)))
tahiti_new$logmean <- log(-tahiti_new$sealevel)-0.5*log(1+coralvar/(tahiti_new$sealevel^2))
tahiti_new$logsd <- sqrt(log(1+coralvar/(tahiti_new$sealevel^2)))
huon_new$logmean <- log(-huon_new$sealevel)-0.5*log(1+coralvar/(huon_new$sealevel^2))
huon_new$logsd <- sqrt(log(1+coralvar/(huon_new$sealevel^2)))
florida_cor_new$logmean <- log(-florida_cor_new$sealevel)-0.5*log(1+coralvar/(florida_cor_new$sealevel^2))
florida_cor_new$logsd <- sqrt(log(1+coralvar/(florida_cor_new$sealevel^2)))



# Simulate n_simul points from each data set

agesb <- as.data.frame(matrix(data=rnorm(n_barbados*n_simul/2,mean=barbados_new$mean_age,sd=sqrt(barbados_new$var)),nrow=n_barbados,ncol=n_simul))
agess2 <- as.data.frame(matrix(data=rnorm(n_sunda2m*n_simul/2,mean=sunda2m_new$mean_age,sd=sqrt(sunda2m_new$var)),nrow=n_sunda2m,ncol=n_simul))
agess7 <- as.data.frame(matrix(data=rnorm(n_sunda7m*n_simul/2,mean=sunda7m_new$mean_age,sd=sqrt(sunda7m_new$var)),nrow=n_sunda7m,ncol=n_simul))
agesbon3 <- as.data.frame(matrix(data=rnorm(n_bonaparte_bm3*n_simul/2,mean=bonaparte_bm3_new$mean_age,sd=sqrt(bonaparte_bm3_new$var)),nrow=n_bonaparte_bm3,ncol=n_simul))
agesbon5 <- as.data.frame(matrix(data=rnorm(n_bonaparte_mm5*n_simul/2,mean=bonaparte_mm5_new$mean_age,sd=sqrt(bonaparte_mm5_new$var)),nrow=n_bonaparte_mm5,ncol=n_simul))
agest <- as.data.frame(matrix(data=rnorm(n_tahiti*n_simul/2,mean=tahiti_new$mean_age,sd=sqrt(tahiti_new$var)),nrow=n_tahiti,ncol=n_simul))
agesh <- as.data.frame(matrix(data=rnorm(n_huon*n_simul/2,mean=huon_new$mean_age,sd=sqrt(huon_new$var)),nrow=n_huon,ncol=n_simul))
agesflcor <- as.data.frame(matrix(data=rnorm(n_florida_cor*n_simul/2,mean=florida_cor_new$mean_age,sd=sqrt(florida_cor_new$var)),nrow=n_florida_cor,ncol=n_simul))
agestP21 <- as.data.frame(matrix(data=rnorm(n_tosP21*n_simul/2,mean=tosP21_new$mean_age,sd=sqrt(tosP21_new$var)),nrow=n_tosP21,ncol=n_simul))
agestP50 <- as.data.frame(matrix(data=rnorm(n_tosP50*n_simul/2,mean=tosP50_new$mean_age,sd=sqrt(tosP50_new$var)),nrow=n_tosP50,ncol=n_simul))
agestP66 <- as.data.frame(matrix(data=rnorm(n_tosP66*n_simul/2,mean=tosP66_new$mean_age,sd=sqrt(tosP66_new$var)),nrow=n_tosP66,ncol=n_simul))
sealevb <- as.data.frame(matrix(data=-rlnorm(n_barbados*n_simul/2,meanlog=barbados_new$logmean,sdlog=barbados_new$logsd),nrow=n_barbados,ncol=n_simul))
sealevs2 <- as.data.frame(matrix(data=runif(n_sunda2m*n_simul/2,min=sunda2m_new$sealevel-2,max=sunda2m_new$sealevel+2),nrow=n_sunda2m,ncol=n_simul))
sealevs7 <- as.data.frame(matrix(data=runif(n_sunda7m*n_simul/2,min=sunda7m_new$sealevel-7,max=sunda7m_new$sealevel+7),nrow=n_sunda7m,ncol=n_simul))
sealevbon3 <- as.data.frame(matrix(data=runif(n_bonaparte_bm3*n_simul/2,min=bonaparte_bm3_new$sealevel-3,max=bonaparte_bm3_new$sealevel+3),nrow=n_bonaparte_bm3,ncol=n_simul))
sealevbon5 <- as.data.frame(matrix(data=-rlnorm(n_bonaparte_mm5*n_simul/2,meanlog=bonaparte_mm5_new$logmean,sdlog=bonaparte_mm5_new$logsd),nrow=n_bonaparte_mm5,ncol=n_simul))
sealevt <- as.data.frame(matrix(data=-rlnorm(n_tahiti*n_simul/2,meanlog=tahiti_new$logmean,sdlog=tahiti_new$logsd),nrow=n_tahiti,ncol=n_simul))
sealevh <- as.data.frame(matrix(data=-rlnorm(n_huon*n_simul/2,meanlog=huon_new$logmean,sdlog=huon_new$logsd),nrow=n_huon,ncol=n_simul))
sealevflcor <- as.data.frame(matrix(data=-rlnorm(n_florida_cor*n_simul/2,meanlog=florida_cor_new$logmean,sdlog=florida_cor_new$logsd),nrow=n_florida_cor,ncol=n_simul))
sealevtP21 <- as.data.frame(matrix(data=runif(n_tosP21*n_simul/2,min=tosP21_new$sealevel-0.21,max=tosP21_new$sealevel+0.21),nrow=n_tosP21,ncol=n_simul))
sealevtP50 <- as.data.frame(matrix(data=runif(n_tosP50*n_simul/2,min=tosP50_new$sealevel-0.5,max=tosP50_new$sealevel+0.5),nrow=n_tosP50,ncol=n_simul))
sealevtP66 <- as.data.frame(matrix(data=runif(n_tosP66*n_simul/2,min=tosP66_new$sealevel-0.66,max=tosP66_new$sealevel+0.66),nrow=n_tosP66,ncol=n_simul))


age <- rbind(agesb,agess2,agess7,agesbon3,agesbon5,agest,agesh,agesflcor,agestP21,agestP50,agestP66)
sealev <- rbind(sealevb,sealevs2,sealevs7,sealevbon3,sealevbon5,sealevt,sealevh,sealevflcor,sealevtP21,sealevtP50,sealevtP66)


#  ordering the data and put one simulation in each row of age and sealev

for (i in 1:n_simul) {
 perm <- order(age[,i])
 age[,i] <- age[perm,i]
 sealev[,i] <- sealev[perm,i]
 }

# fit a smoothing spline to all the data

smoothed <- matrix(nrow = 220,ncol=n_simul)
firstderiv <- matrix(nrow = 220,ncol=n_simul)
for (i in 1:n_simul){
tmp <- smooth.spline(age[,i],sealev[,i])
smoothed[,i] <- predict(tmp,smooth_x)$y
firstderiv[,i] <- predict(tmp,smooth_x,deriv=1)$y
}

percent5 <- matrix(ncol=220,nrow=2)
percent5d <- matrix(ncol=220,nrow=2)
for (i in 1:220){
	percent5[,i] <-quantile(smoothed[i,],c(0.025,.975),na.rm=T)
    percent5d[,i] <-quantile(firstderiv[i,],c(0.025,.975),na.rm=T)

}

percent1 <- matrix(ncol=220,nrow=2)
percent1d <- matrix(ncol=220,nrow=2)
for (i in 1:220){
	percent1[,i] <-quantile(smoothed[i,],c(0.01,.99),na.rm=T)
	percent1d[,i] <-quantile(firstderiv[i,],c(0.01,.99),na.rm=T)

}

percent33 <- matrix(ncol=220,nrow=2)
percent33d <- matrix(ncol=220,nrow=2)
for (i in 1:220){
	percent33[,i] <-quantile(smoothed[i,],c(0.165,.835),na.rm=T)
	percent33d[,i] <-quantile(firstderiv[i,],c(0.165,.835),na.rm=T)

}



plot(age$V1,sealev$V1,ylim=c(-130,0),xlim=c(0,22),type='n',xlab='ka BP',ylab='Sea Level')

for (i in 1:n_simul){
lines(smooth_x,smoothed[,i],col='pink')}

lines(smooth_x,percent5[1,],col='blue')
lines(smooth_x,percent5[2,],col='blue')
lines(smooth_x,percent1[1,],col='red')
lines(smooth_x,percent1[2,],col='red')
lines(smooth_x,percent33[1,],col='black')
lines(smooth_x,percent33[2,],col='black')


points(barbados_new$mean_age, barbados_new$sealevel,col='blue')
points(sunda2m_new$mean_age, sunda2m_new$sealevel,col='red')
points(sunda7m_new$mean_age, sunda7m_new$sealevel,col='yellow')
points(bonaparte_bm3_new$mean_age, bonaparte_bm3_new$sealevel,col='green')
points(bonaparte_mm5_new$mean_age, bonaparte_mm5_new$sealevel,col='green')
points(tahiti_new$mean_age, tahiti_new$sealevel,col='grey')
points(huon_new$mean_age, huon_new$sealevel,col='darkgreen')
points(florida_cor_new$mean_age, florida_cor_new$sealevel,col='aquamarine')
points(tosP21_new$mean_age, tosP21_new$sealevel,col='violet')
points(tosP50_new$mean_age, tosP50_new$sealevel,col='violet')
points(tosP66_new$mean_age, tosP66_new$sealevel,col='violet')


#  plotting the derivatives
quartz()

plot(age$V1,sealev$V1,ylim=c(-20,40),xlim=c(0,22),type='n',xlab='ka BP',ylab='Derivative of Sea Level')

for (i in 1:n_simul){
lines(smooth_x,-firstderiv[,i],col='pink')}

lines(smooth_x,-percent5d[1,],col='blue')
lines(smooth_x,-percent5d[2,],col='blue')
lines(smooth_x,-percent1d[1,],col='red')
lines(smooth_x,-percent1d[2,],col='red')
lines(smooth_x,-percent33d[1,],col='black')
lines(smooth_x,-percent33d[2,],col='black')

# plotting the histograms

quartz()
par(mfcol=c(2,2))

hist(smoothed[smooth_x==5],main='5 ka BP',xlab ='Sea Level')
hist(smoothed[smooth_x==10],main='10 ka BP',xlab ='Sea Level')
hist(smoothed[150,], main='15 ka BP',xlab ='Sea Level')
hist(smoothed[smooth_x==20], main='20 ka BP',xlab ='Sea Level')

par(mfcol=c(1,1))