#Set working Directory
setwd("")

####################################################################################
# Access libraries of needed functions
####################################################################################
library("rgdal")
library("raster")
library("maps")
library("maptools")
library("rasterVis")
library("mvnmle")
library("mixtools")
library("fossil")
library(sp)
library(maptools)
library(rgdal)
library(classInt)
library(assignR)
library(colorRamps)

####################################################################################
# Define new functions
####################################################################################

## function returns raster of posterior probabilities for univariate normal data
calcCellProb <- function(x,isoscape,std){
  m <- getValues(isoscape)
  s <- getValues(std) # use this is you have raster of variances
  m <- dnorm(x,mean=m,sd=s)
  cell.dens <- setValues(isoscape,m)
  return(cell.dens)
}

## function returns raster of posterior probabilities for bivariate normal data
## x is the unknown tissue of interest, will have two values, one for each isotope
## m is a 2-D vector, all the values in the raster for each isotope
## v is the same as m, but for variances
## r is a single number - the covariance. Can be vector if estimated as non-stationary
## ras is a raster that will serve as a template for the final product
calcCellProb2D <- function(x,m,v,r,ras) {
  pd <- 1/(2*pi*sqrt(v[,1])*sqrt(v[,2])*sqrt(1-r^2))*exp(-(1/(2*(1-r^2)))*
                                                           ((x[1]-m[,1])^2/v[,1]+(x[2]-m[,2])^2/v[,2]-(2*r*(x[1]-m[,1])*
                                                                                                         (x[2]-m[,2]))/(sqrt(v[,1])*sqrt(v[,2]))))
  pdras <- setValues(ras,pd)
  return(pdras)
}

## function returns raster of posterior probability distribution
calcPostProb <- function(x){
  pp <- x/cellStats(x,sum)
  return(pp)
}


## function to get normalized cell probabilites
calcNormProb <- function(x){
  np <- x/cellStats(x,max)
  return(np)
}


####################################################################################
## Data input: read in tissue isotope data, GIS raster models
####################################################################################
iso.data <- read.csv("iso.data_2020.csv") # tissue isotope data
iso.data$d18OVSMOW<-iso.data$d18O*1.03092 +30.92
iso.data$d18Ow<-iso.data$d18OVSMOW*1.59-48.634###Using equation from Dotsika et al. 1.59*iso.data$d18OVSMOW-48.634

iso.data$d18Op<-iso.data$d18OVSMOW/1.015-8.79/1.015###Using equation from Iacumin et al. 1996 ??18OC=1.015(±0.043)* ??18OP+8.79(±0.79)
iso.data$d18Ow<-(iso.data$d18Op-21.28)/0.68 ###Using equation from Hoppe et al. 2006 : ??18OP=21.28(±0.51)+[0.68(±0.04)*??18Ow].
data(wrld_simpl)

### Download d18O isoscape rasters at RCWIP http://www-naweb.iaea.org/napc/ih/IHS_resources_rcwip.html
### input the annual amount-weighted mean and sd from RCWIP
rain_d18O_aaw <- raster("O18_H2_Ann_GS\\RCWIP_grid_5_14_100_20130502.tif")  # this is baseline isoscape
rain_d18O_err_aaw <- raster("O18_H2_Ann_GS\\RCWIP_grid_err_5_14_100_20130502.tif")  # this is baseline isoscape

### Download 87Sr/86Sr isoscape rasters and uncertainty at https://drive.google.com/drive/folders/1g9rCGo3Kd3hz2o5JKkSbgNsGJclvsuQm?usp=sharing
### See Bataille et al. 2020 for details
rf_sr<-raster("rf_plantsoilmammal1.tif")
rf_sr_err <-raster("srse.tif")

#Clip data to study area Europe
europe<-as.vector(c(-500000,1000000,5000000,6500000))

rangemap<-crop(rf_sr/rf_sr, europe, snap='near')
writeRaster(rangemap, "range.tif",overwrite=TRUE)
sr<-rf_sr*rangemap
sr.se<-rf_sr_err*rangemap

###Reproject d18O isoscape to Sr projection and resolution
d18O <- projectRaster(rain_d18O_aaw,sr, method="bilinear")
d18O.se<-projectRaster(rain_d18O_err_aaw,sr, method="bilinear")+1###Add 1 per mile uncertainty due to regression equation

###Download d34S isoscape generated from R_Script_S1 from this manuscript
d34S<-raster("rf_d34S.tif")
d34S<-d34S*rangemap
###Define uncertainty for d34S at 2.8 per mile constant over the entire area
d34S.se<-d34S/d34S*2.8

###Plot points on d34S isoscape
d34S_proj<-project(as.matrix(iso.data[,c("Lat","Long")]), "+proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs")
iso.data$X<-d34S_proj[,1]
iso.data$Y<-d34S_proj[,2]


###########################d18O single isotope assignments#########################################################
####################################################################################
distance.list <- vector("list", length(iso.data[,1]))
origins <- stack()
for (i in seq(along=iso.data[,1])){
  pp <- calcCellProb(iso.data$d18Ow[i],d18O,d18O.se) # compute probs
  np <- calcNormProb(pp)
  origins <- raster::stack(origins,np)
}


## summary map plots
pdf("maps_d18O_annual.pdf")
opar<-par()
par(mfrow=c(2,2),mar=c(2,3,2,2)+0.01)
xl<-length(iso.data[,1])
for (i in seq(along=iso.data[,1])){
  plot(origins[[i]],axes = FALSE)
  #plot(wrld_simpl,add=T)
  # plot known origin and individual label info
  #points(s.data$Long,s.data$Lat,pch=21, col="black", bg=colcode, lwd=0.4, cex=symbol.size)
  points(iso.data$X[i],iso.data$Y[i],col="black",cex=0.75, pch=4)
  #text(-5,47,paste("??34S=",iso.data$d34S[i],"???"),cex=0.70)
  text(0,6.4e6,paste("ID",iso.data$ID[i]),cex=0.70)
  #legend(-13, 65, legend=names(attr(colcode, "table")),
  #      fill=attr(colcode, "palette"), cex=0.6, bty="n")
  #hist(distance.list[[i]],xlim=c(0,xl),xlab="",ylab="",main="")
}
par(opar)
dev.off()
rm(i)

###Store posterior probability maps in raster stack
O<-origins

###Duchal individuals
brk <- c(0, 0.25,0.5,0.75,1)
par(mfrow=c(2,1))
duchal_d18O=O[[7]]+O[[8]]+O[[9]]+O[[10]]
plot(duchal_d18O/4, breaks=brk, col=rev(heat.colors(4)))

###French army individuals
french_d18O=O[[1]]+O[[2]]+O[[3]]+O[[4]]+O[[5]]+O[[6]]
plot(french_d18O/6, breaks=brk, col=rev(heat.colors(4)))

###Local individuals
O_local<-O[[11]]+O[[12]]+O[[13]]+O[[14]]+O[[15]]+O[[16]]+O[[17]]
O_local <- O_local/cellStats(O_local,sum)

qprob<-assignR::qtlRaster(O_local, 0.333,thresholdType="prob",genplot=TRUE, outDir = NULL)
plot(qprob)
cellStats(qprob,sum)

###########################SR single isotope assignments#########################################################

distance.list <- vector("list", length(iso.data[,1]))
origins <- stack()
for (i in seq(along=iso.data[,1])){
  pp <- calcCellProb(iso.data$Sr[i],sr,sr.se) # compute probs
  np <- calcNormProb(pp)
  origins <- raster::stack(origins,np)
  
}

## summary map plots
pdf("maps_sr.pdf")
opar<-par()
par(mfrow=c(2,2),mar=c(2,3,2,2)+0.01)
xl<-length(iso.data[,1])
for (i in seq(along=iso.data[,1])){
  plot(origins[[i]],axes = FALSE)
  #plot(wrld_simpl,add=T)
  # plot known origin and individual label info
  #points(s.data$Long,s.data$Lat,pch=21, col="black", bg=colcode, lwd=0.4, cex=symbol.size)
  points(iso.data$X[i],iso.data$Y[i],col="black",cex=0.75, pch=4)
  #text(-5,47,paste("??34S=",iso.data$d34S[i],"???"),cex=0.70)
  text(0,6.4e6,paste("ID",iso.data$ID[i]),cex=0.70)
  #legend(-13, 65, legend=names(attr(colcode, "table")),
  #      fill=attr(colcode, "palette"), cex=0.6, bty="n")
  #hist(distance.list[[i]],xlim=c(0,xl),xlab="",ylab="",main="")
}
par(opar)
dev.off()
rm(i)

###Store posterior probability maps in raster stack
Sr87_86Sr<-origins

###Store posterior probability maps in raster stacks to generate figures
###Duchal individuals
brk <- c(0, 0.25,0.5,0.75,1)
par(mfrow=c(2,1))
duchal_Sr87_86Sr=Sr87_86Sr[[7]]+Sr87_86Sr[[8]]+Sr87_86Sr[[9]]+Sr87_86Sr[[10]]
plot(duchal_Sr87_86Sr/4, breaks=brk, col=rev(heat.colors(4)))

###French army individuals
french_Sr87_86Sr=Sr87_86Sr[[1]]+Sr87_86Sr[[2]]+Sr87_86Sr[[3]]+Sr87_86Sr[[4]]+Sr87_86Sr[[5]]+Sr87_86Sr[[6]]
plot(french_Sr87_86Sr/6, breaks=brk, col=rev(heat.colors(4)))

###Local individuals
Sr_local<-Sr87_86Sr[[11]]+Sr87_86Sr[[12]]+Sr87_86Sr[[13]]+Sr87_86Sr[[14]]+Sr87_86Sr[[15]]+Sr87_86Sr[[16]]+Sr87_86Sr[[17]]
Sr_local <- Sr_local/cellStats(Sr_local,sum)

qprob<-assignR::qtlRaster(Sr_local, 0.333,thresholdType="prob",genplot=TRUE, outDir = NULL)
plot(qprob)
cellStats(qprob,sum)


###########################d34S single assignments#########################################################
distance.list <- vector("list", length(iso.data[,1]))
origins <- stack()
for (i in seq(along=iso.data[,1])){
  pp <- calcCellProb(iso.data$d34S[i],d34S,d34S.se) # compute probs
  np <- calcNormProb(pp)
  origins <- raster::stack(origins,np)
  # rem the next 5 lines if assigning unknown origin samples
}

###Summary map plots
pdf("maps_d34S.pdf")
opar<-par()
par(mfrow=c(2,2),mar=c(2,3,2,2)+0.01)
xl<-length(iso.data[,1])
for (i in seq(along=iso.data[,1])){
  plot(origins[[i]],axes = FALSE)
  #plot(wrld_simpl,add=T)
  # plot known origin and individual label info
  #points(s.data$Long,s.data$Lat,pch=21, col="black", bg=colcode, lwd=0.4, cex=symbol.size)
  points(iso.data$X[i],iso.data$Y[i],col="black",cex=0.75, pch=4)
  #text(-5,47,paste("??34S=",iso.data$d34S[i],"???"),cex=0.70)
  text(0,6.4e6,paste("ID",iso.data$ID[i]),cex=0.70)
  #legend(-13, 65, legend=names(attr(colcode, "table")),
  #      fill=attr(colcode, "palette"), cex=0.6, bty="n")
  #hist(distance.list[[i]],xlim=c(0,xl),xlab="",ylab="",main="")
}
par(opar)
dev.off()
rm(i)

###Store posterior probablity raster stack
S34<-origins

###Duchal individuals
brk <- c(0, 0.25,0.5,0.75,1)
par(mfrow=c(2,1))
duchal_S=S34[[7]]+S34[[8]]+S34[[9]]+S34[[10]]
plot(duchal_S/4, breaks=brk, col=rev(heat.colors(4)))

###French army individuals
french_S=S34[[2]]+S34[[3]]+S34[[4]]+S34[[5]]
plot(french_S/4, breaks=brk, col=rev(heat.colors(4)))

###Local individuals
S_local<-S34[[11]]+S34[[12]]+S34[[13]]+S34[[15]]+S34[[16]]+S34[[17]]
S_local <- S_local/cellStats(S_local,sum)

qprob<-assignR::qtlRaster(S_local, 0.333,thresholdType="prob",genplot=TRUE, outDir = NULL)
cellStats(qprob,sum)
plot(qprob)



###########################d18O and 87Sr/86Sr dual assignments#########################################################
###Combine d18O and 87Sr/86Sr assuming independence
OSr<-O*Sr87_86Sr


pdf("maps_OSr.pdf")
opar<-par()
par(mfrow=c(2,2),mar=c(2,3,2,2)+0.01)
xl<-length(iso.data[,1])
for (i in seq(along=iso.data[,1])){
  plot(OSr[[i]],axes = FALSE)
  #plot(wrld_simpl,add=T)
  # plot known origin and individual label info
  #points(s.data$Long,s.data$Lat,pch=21, col="black", bg=colcode, lwd=0.4, cex=symbol.size)
  points(iso.data$X[i],iso.data$Y[i],col="black",cex=0.75, pch=4)
  #text(-5,47,paste("??34S=",iso.data$d34S[i],"???"),cex=0.70)
  text(0,6.4e6,paste("ID",iso.data$ID[i]),cex=0.70)
  #legend(-13, 65, legend=names(attr(colcode, "table")),
  #      fill=attr(colcode, "palette"), cex=0.6, bty="n")
  #hist(distance.list[[i]],xlim=c(0,xl),xlab="",ylab="",main="")
}
par(opar)
dev.off()
rm(i)

###Duchal individuals
brk <- c(0, 0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1)
par(mfrow=c(2,1))
duchal_OSr=OSr[[7]]+OSr[[8]]+OSr[[9]]+OSr[[10]]
plot(duchal_OSr/4, breaks=brk, col=rev(heat.colors(10)))

###French Army individuals
french_OSr=OSr[[1]]+OSr[[2]]+OSr[[3]]+OSr[[4]]+OSr[[5]]+OSr[[6]]
plot(french_OSr/6, breaks=brk, col=rev(heat.colors(10)))

###Summary plot individuals by group french army vs. Duchal army
brk <- c(0, 0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1)
pal.1 <- colorRampPalette(c("gray92","gray1"), bias=1)
pdf("OSr_group.pdf",width=8, height=4)
par(mfrow=c(1,2))
par(mar=c(1,1,1,1))
plot(french_OSr/4, breaks=brk, col=pal.1(10),axes = FALSE)
plot(duchal_OSr/4, breaks=brk, col=pal.1(10),axes = FALSE)
dev.off()

###Local individuals
OSr_local<-OSr[[11]]+OSr[[12]]+OSr[[13]]+OSr[[14]]+OSr[[15]]+OSr[[16]]+OSr[[17]]
OSr_local <- OSr_local/cellStats(OSr_local,sum)

qprob<-assignR::qtlRaster(OSr_local, 0.333,thresholdType="prob",genplot=TRUE, outDir = NULL)
cellStats(qprob,sum)


###########################d34S and 87Sr/86Sr dual assignments#########################################################
###Combine d34S and 87Sr/86Sr probability maps assuming independence
SSr<-S34*Sr87_86Sr

pdf("maps_SSr.pdf")
opar<-par()
par(mfrow=c(2,2),mar=c(2,3,2,2)+0.01)
xl<-length(iso.data[,1])
for (i in seq(along=iso.data[,1])){
  plot(SSr[[i]],axes = FALSE)
  #plot(wrld_simpl,add=T)
  # plot known origin and individual label info
  #points(s.data$Long,s.data$Lat,pch=21, col="black", bg=colcode, lwd=0.4, cex=symbol.size)
  points(iso.data$X[i],iso.data$Y[i],col="black",cex=0.75, pch=4)
  #text(-5,47,paste("??34S=",iso.data$d34S[i],"???"),cex=0.70)
  text(0,6.4e6,paste("ID",iso.data$ID[i]),cex=0.70)
  #legend(-13, 65, legend=names(attr(colcode, "table")),
  #      fill=attr(colcode, "palette"), cex=0.6, bty="n")
  #hist(distance.list[[i]],xlim=c(0,xl),xlab="",ylab="",main="")
}
par(opar)
dev.off()
rm(i)

###Duchal individuals
brk <- c(0, 0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1)
par(mfrow=c(2,1))
duchal_SSr=SSr[[7]]+SSr[[8]]+SSr[[9]]+SSr[[10]]
plot(duchal_SSr/4, breaks=brk, col=rev(heat.colors(10)))

###French army individuals
french_SSr=SSr[[2]]+SSr[[3]]+SSr[[4]]+SSr[[5]]
plot(french_SSr/4, breaks=brk, col=rev(heat.colors(10)))

###Summary plot individuals by group french army vs. Duchal army
brk <- c(0, 0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1)
pal.1 <- colorRampPalette(c("gray92","gray1"), bias=1)
pdf("SSr_group.pdf",width=8, height=4)
par(mfrow=c(1,2))
par(mar=c(1,1,1,1))
plot(french_SSr/4, breaks=brk, col=pal.1(10),axes = FALSE)
plot(duchal_SSr/4, breaks=brk, col=pal.1(10),axes = FALSE)
dev.off()

###Local individuals
SSr_local<-SSr[[11]]+SSr[[12]]+SSr[[13]]+SSr[[15]]+SSr[[16]]+SSr[[17]]
SSr_local <- SSr_local/cellStats(SSr_local,sum)

qprob<-assignR::qtlRaster(SSr_local, 0.333,thresholdType="prob",genplot=TRUE, outDir = NULL)
cellStats(qprob,sum)



###########################d18O and d34S dual assignments#########################################################
###Combine d18O and d34S probability maps assuming independence
OS<-O*S34

pdf("maps_SO.pdf")
opar<-par()
par(mfrow=c(2,2),mar=c(2,3,2,2)+0.01)
xl<-length(iso.data[,1])
for (i in seq(along=iso.data[,1])){
  plot(OS[[i]],axes = FALSE)
  #plot(wrld_simpl,add=T)
  # plot known origin and individual label info
  #points(s.data$Long,s.data$Lat,pch=21, col="black", bg=colcode, lwd=0.4, cex=symbol.size)
  points(iso.data$X[i],iso.data$Y[i],col="black",cex=0.75, pch=4)
  #text(-5,47,paste("??34S=",iso.data$d34S[i],"???"),cex=0.70)
  text(0,6.4e6,paste("ID",iso.data$ID[i]),cex=0.70)
  #legend(-13, 65, legend=names(attr(colcode, "table")),
  #      fill=attr(colcode, "palette"), cex=0.6, bty="n")
  #hist(distance.list[[i]],xlim=c(0,xl),xlab="",ylab="",main="")
}
par(opar)
dev.off()
rm(i)

###Duchal individuals
duchal_OS=OS[[7]]+OS[[8]]+OS[[9]]+OS[[10]]
plot(duchal_OS/4, breaks=brk, col=rev(heat.colors(10)))

###French army individuals
french_OS=OS[[2]]+OS[[3]]+OS[[4]]+OS[[5]]
plot(french_OS/4, breaks=brk, col=rev(heat.colors(10)))

###Summary plot individuals by group french army vs. Duchal army
brk <- c(0, 0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1)
pal.1 <- colorRampPalette(c("gray92","gray1"), bias=1)
pdf("OS_group.pdf",width=8, height=4)
par(mfrow=c(1,2))
par(mar=c(1,1,1,1))
plot(french_OS/4, breaks=brk, col=pal.1(10),axes = FALSE)
plot(duchal_OS/4, breaks=brk, col=pal.1(10),axes = FALSE)
dev.off()

###Local individuals
OS_local<-OS[[11]]+OS[[12]]+OS[[13]]+OS[[15]]+OS[[16]]+OS[[17]]
OS_local <- OS_local/cellStats(OS_local,sum)

qprob<-assignR::qtlRaster(OS_local, 0.333,thresholdType="prob",genplot=TRUE, outDir = NULL)
cellStats(qprob,sum)



###########################Triple isotope assignments#########################################################
###Combine triple isotope probability maps assuming independence
OSSr<-O*Sr87_86Sr*S34

pdf("maps_OSSr.pdf")
opar<-par()
par(mfrow=c(2,2),mar=c(2,3,2,2)+0.01)
xl<-length(iso.data[,1])
for (i in seq(along=iso.data[,1])){
  plot(OSSr[[i]],axes = FALSE)
  #plot(wrld_simpl,add=T)
  # plot known origin and individual label info
  #points(s.data$Long,s.data$Lat,pch=21, col="black", bg=colcode, lwd=0.4, cex=symbol.size)
  points(iso.data$X[i],iso.data$Y[i],col="black",cex=0.75, pch=4)
  #text(-5,47,paste("??34S=",iso.data$d34S[i],"???"),cex=0.70)
  text(0,6.4e6,paste("ID",iso.data$ID[i]),cex=0.70)
  #legend(-13, 65, legend=names(attr(colcode, "table")),
  #      fill=attr(colcode, "palette"), cex=0.6, bty="n")
  #hist(distance.list[[i]],xlim=c(0,xl),xlab="",ylab="",main="")
}
par(opar)
dev.off()
rm(i)

###Duchal individuals
duchal_OSSr=OSSr[[7]]+OSSr[[8]]+OSSr[[9]]+OSSr[[10]]
plot(duchal_OSSr/4, breaks=brk, col=rev(heat.colors(10)))

###French army individuals
french_OSSr=OSSr[[2]]+OSSr[[3]]+OSSr[[4]]+OSSr[[5]]
plot(french_OSSr/4, breaks=brk, col=rev(heat.colors(10)))

french_cav_OSSr=OSSr[[4]]+OSSr[[5]]

french_nocav_OSSr=OSSr[[2]]+OSSr[[3]]

###Summary plot individuals by group french army vs. Duchal army
brk <- c(0, 0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1)
pal.1 <- colorRampPalette(c("gray92","gray1"), bias=1)
pdf("OSSr_group_cav.pdf",width=8, height=4)
par(mfrow=c(1,3))
par(mar=c(1,1,1,1))
plot((duchal_OSSr/4*1.2), breaks=brk, col=matlab.like2(10),axes = FALSE)
#plot(french_OSSr/4, breaks=brk, col=matlab.like2(10),axes = FALSE)
plot(french_cav_OSSr/2*1.2, breaks=brk, col=matlab.like2(10),axes = FALSE)
plot(french_nocav_OSSr/2*1.2, breaks=brk, col=matlab.like2(10),axes = FALSE)
dev.off()

###Local individuals
OSSr_local<-OSSr[[11]]+OSSr[[12]]+OSSr[[13]]+OSSr[[15]]+OSSr[[16]]+OSSr[[17]]
OSSr_local <- OSSr_local/cellStats(OSSr_local,sum)

qprob<-assignR::qtlRaster(OSSr_local, 0.333,thresholdType="prob",genplot=TRUE, outDir = NULL)
cellStats(qprob,sum)