#################################################################################
# These functions can be used to investigate the effect of lipid extraction or  #
# the presence of lipids on the solution of a mixing analysis carried out       #
# with Siar (Parnell et al. 2010). A stepwise increase in d13C (corresponding   #
# to lipid exraction) or decrease in d13C (corresponding to the presence of     #
# lipids) is simulated for one or several sources, as well as for the consumer, #
# and Siar is run for all considered combinations of modified data.             # 
#-------------------------------------------------------------------------------#																																						
# Predator data and source data should be entered as for Siar.                  #
# Values for carbon should be first and values for Nitrogen second.    	  			#
#-------------------------------------------------------------------------------#
# The script is linked to the following paper:                                  #
# Tarroux, A., D. Ehrich, N. Lecomte, T. Jardine, J. Bety and                   #
# D. Berteaux (2010) Sensitivity of stable isotope mixing models to variation   #
# in isotopic ratios: evaluating consequences of lipid extraction. Methods in   #
# Ecology and Evolution, accepted 13-04-2010                                    #
#################################################################################  


# Main simulation function
# The result of the function is a list with the outputs of all siar runs
# In addition, the function writes a file with a summary of the simulations
# called fatsim-output.txt

Fatsim <- function(pred, sources, maxmod, incrC=0.5, correc=0, conc = 0,  
    incrN=0, consumer=0, incrNcons = 0, shift=0.1, cct = "mode", int = 95,
    burn=50000, iter=200000, thin=15) {

  nbsources <- dim(sources)[1]
  if (sum(incrN)==0) incrN <- c(rep(0, nbsources))
  if (dim(pred)[1] == 1) solo <- TRUE else solo <- FALSE

  modwhich <- which(maxmod != 0)
  nbmod <- length(modwhich)
  nbsteps <- maxmod[modwhich] / incrC 
  steps <- 1/incrC
  nbvar <- prod(abs(nbsteps) + 1)
  
  if (nbmod > 0) {
  incrCs <- vector(mode="list", nbmod)
  for (i in 1:nbmod) {
    incrCs[[i]] <- c(0:nbsteps[i])/steps
    }     # list of all steps considered for each source to modify

  mods <- matrix(0, nbvar, nbmod)
  mods [ ,1] <- incrCs[[1]]
  a <- 1
  vect <- incrCs[[1]]
  while (a < nbmod) {
    b <- length (vect)
    a <- a + 1
    vect <- sort(c(rep(incrCs[[a]], b)))
    mods[ , a] <- vect 
  }  # matrix whith all combination of modifications for the sources to modify
  if (dim(mods)[2] > 1){
    stemp <- apply(mods, 1, sum)
    mods <- mods[stemp!=0, ]
    mods <- rbind(c(rep(0, nbmod)), mods)
    rm(stemp)   
    } else {
    mods1 <- mods[mods[ ,1]!=0, ]
    mods <- as.matrix(c(0, mods1), dim=c(length(mods1+1), 1))
    rm(mods1)
    }
  
  variants <- matrix(0, nbvar, nbsources)
  for (i in 1:nbmod) {
    variants[ , modwhich[i]] <- mods[ , i]
  } # matrix whith all combination of modifications for the sources to modify,
    # including sources which will not be modified
   } else {
   variants <- matrix(0, nbvar, nbsources)   }

if (solo == FALSE){
if (consumer==0) {  
  mixes <- vector("list", (nbvar + 6))

  for (i in 1:nbvar) {
    sources.mod <- sources
    sources.mod[ , 4] <- sources[ , 4] + incrN
    sources.mod[ , 2] <- sources[ , 2] + variants[i, ]
    mixes[[i]] <-  siarmcmcdirichletv4(pred, sources.mod, corrections=correc, 
          concdep=conc, iterations=iter, burnin=burn, thinby=thin)[[15]]
    } 
  allvars <- variants
  }
  
if (consumer != 0){
  nbconssteps <- consumer / incrC
  conssteps <- c(0:nbconssteps)/steps
  mixes <- vector("list", (nbvar * (abs(nbconssteps) + 1) + 6))
  allvars <- cbind(c(rep(conssteps[1], nbvar)), variants)
      
  for (f in 1:(abs(nbconssteps)+1)) {
    pred.mod <- pred
    pred.mod[ , 2] <- pred [ ,2] + incrNcons
    pred.mod[ , 1] <- pred [ ,1] + conssteps[f]  
    for (i in 1:nbvar) {
      sources.mod <- sources
      sources.mod[ , 4] <- sources[ , 4] + incrN
      sources.mod[ , 2] <- sources[ , 2] + variants[i, ]
      mixes[[(f-1)*nbvar + i]] <-  siarmcmcdirichletv4(pred.mod, sources.mod, 
            corrections=correc, concdep=conc, iterations=iter, burnin=burn, thinby=thin)[[15]]
      } 
      if(f != 1) allvars <- rbind(allvars, cbind(c(rep(conssteps[f], nbvar)), variants))
      } }
  }

if (solo == TRUE){
if (consumer==0) {  
  mixes <- vector("list", (nbvar + 6))                       

  for (i in 1:nbvar) {
    sources.mod <- sources
    sources.mod[ , 4] <- sources[ , 4] + incrN
    sources.mod[ , 2] <- sources[ , 2] + variants[i, ]
    mixes[[i]] <-  siarsolomcmcv4(pred, sources.mod, corrections=correc, 
          concdep=conc, iterations=iter, burnin=burn, thinby=thin)[[15]]
    } 
  allvars <- variants
  }
  
if (consumer != 0){
  nbconssteps <- consumer / incrC
  conssteps <- c(0:nbconssteps)/steps
  mixes <- vector("list", (nbvar * (abs(nbconssteps) + 1) + 6))
  allvars <- cbind(c(rep(conssteps[1], nbvar)), variants)
      
  for (f in 1:(abs(nbconssteps)+1)) {
    pred.mod <- pred
    pred.mod[ , 2] <- pred [ ,2] + incrNcons
    pred.mod[ , 1] <- pred [ ,1] + conssteps[f]  
    for (i in 1:nbvar) {
      sources.mod <- sources
      sources.mod[ , 4] <- sources[ , 4] + incrN
      sources.mod[ , 2] <- sources[ , 2] + variants[i, ]
      mixes[[(f-1)*nbvar + i]] <-  siarsolomcmcv4(pred.mod, sources.mod, 
            corrections=correc, concdep=conc, iterations=iter, burnin=burn, thinby=thin)[[15]]
      } 
      if(f != 1) allvars <- rbind(allvars, cbind(c(rep(conssteps[f], nbvar)), variants))
      } }
  }

# prepare the output

nrepl <- dim(allvars)[1]

if (cct == "median") {
  obs <- apply(mixes[[1]], 2, median)[1:nbsources]}
if (cct == "mean") {
  obs <- apply(mixes[[1]], 2, mean)[1:nbsources]}
if (cct == "mode") {
  obs <- vector("numeric", nbsources)
  probs = c(95, 75, 50)
  for (b in 1:nbsources){  
    obs[b] <- hdr(mixes[[1]][, b], probs, h = bw.nrd0(mixes[[1]][, b]))$mode
  }
  rm(probs)}

a <- length (mixes)
mixes[[a-5]] <- incrN
mixes[[a-4]] <- incrNcons
mixes[[a-3]] <- allvars
mixes[[a-2]] <- pred
mixes[[a-1]] <- sources
mixes[[a]] <- correc

#write the output

cat("Output of the R function Fatsim:", "\n", "\n", file="fatsim-output.txt", sep="")
  if (consumer == 0) cat("For this run the consumer has not been modified.", "\n", 
  file="fatsim-output.txt", sep="", append=TRUE)
if (solo==T) cat("For this run a single consumer individual was used.", "\n", sep="\t", 
  file="fatsim-output.txt", append=TRUE) 
if (solo==F) cat("For this run several consumer individuals were used.", "\n", sep="\t", 
  file="fatsim-output.txt", append=TRUE) 
cat("For the combinations marked with DIFFERENT, the proportion of at least one of",
  "\n", "the sources in the estimated mixture changed with more than ", 
  shift, "\n", "compared to the mixture resulting from the original data.", "\n", 
  file="fatsim-output.txt", sep="", append=TRUE) 
if(cct == "median") cat("Comparisons and proportion estimates are based on the median of the posterior distributions.", "\n", 
  file="fatsim-output.txt", sep="", append=TRUE) 
if(cct == "mean") cat("Comparisons and proportion estimates are based on the mean of the posterior distributions.", "\n", 
  file="fatsim-output.txt", sep="", append=TRUE) 
if(cct == "mode") cat("Comparisons and proportion estimates are based on the mode of the posterior distributions.", "\n", 
  file="fatsim-output.txt", sep="", append=TRUE)
cat("Line 0 in the table below refers to the original data.", "\n", 
  file="fatsim-output.txt", sep="", append=TRUE) 
cat("\n", "\n", file="fatsim-output.txt", sep="", append=TRUE)
temp <- vector("character", nbsources*3)
    b <- 1
    c <- 1
    while (b < nbsources*3){
      temp[b] <- paste(sources[c, 1], "_est", sep="")
      temp[b+1] <- paste(sources[c, 1], "_lo", sep="")
      temp[b+2] <- paste(sources[c, 1], "_up", sep="")  
      b <- b + 3
      c <- c + 1
      }
if (consumer ==0 ) cat("line", as.vector(sources[ , 1]), "maximal difference", "\t", temp, "\n", sep="\t", 
  file="fatsim-output.txt", append=TRUE) 
if (consumer !=0 ) cat("line", "consumer", as.vector(sources[ , 1]), "maximal difference", "\t", temp, "\n", sep="\t", 
  file="fatsim-output.txt", append=TRUE) 
  
for (i in 1:nrepl){
  cat(i-1, allvars[i, ], sep="\t", file="fatsim-output.txt", append=TRUE) 
  if (cct == "median") {
    modres <- apply(mixes[[i]], 2, median)[1:nbsources] }
  if (cct == "mean") {
    modres <- apply(mixes[[i]], 2, mean)[1:nbsources] }
  if (cct == "mode") {
    modres <- vector("numeric", nbsources)
    probs = c(95, 75, 50)
    for (b in 1:nbsources){  
      modres[b] <- hdr(mixes[[i]][, b], probs, h = bw.nrd0(mixes[[i]][, b]))$mode
      }
    rm(probs)}
 
  cat("\t", max(abs(obs-modres)), file="fatsim-output.txt", sep="", append=TRUE)
  if (max(abs(obs-modres)) >= shift) {
  cat("\t", "DIFFERENT", file="fatsim-output.txt", sep="", append=TRUE)
   } else cat("\t", file="fatsim-output.txt", sep="", append=TRUE)

# check for problems with the data. Code taken from the function siarhdrs 
  hdrsummary <- matrix(0, ncol = 4, nrow = ncol(mixes[[i]]))
  for (j in 1:ncol(mixes[[i]])) {
        if (all(mixes[[i]][, j] == 0)) {
            hdrsummary[j, ] <- 0
            }
        else {
            temp <- hdr(mixes[[i]][, j], h = bw.nrd0(mixes[[i]][, j]))
            hdrsummary[j, 1] <- max(0, temp$hdr[2, 1])
            hdrsummary[j, 2] <- temp$hdr[2, 2]
            hdrsummary[j, 3] <- temp$mode
            hdrsummary[j, 4] <- mean(mixes[[i]][, j])
        }    }
  if (any(hdrsummary > 5)) {
    cat("\t", "problem with the data!", "\t", file="fatsim-output.txt", sep="", append=TRUE)
    } else {cat("\t", "\t", file="fatsim-output.txt", sep="", append=TRUE) }
  
# add credibility intervals
    ints <- vector("numeric", (nbsources*3))
    probs = c(int, 75, 50)
    b <- 1
    c <- 1
    while (b < nbsources*3){  
      ints[(b+1)] <- max(0, hdr(mixes[[i]][, c], probs, h = bw.nrd0(mixes[[i]][, c]))$hdr[1, 1])
      ints[(b+2)] <- hdr(mixes[[i]][, c], probs, h = bw.nrd0(mixes[[i]][, c]))$hdr[1, 2]
      if (cct == "median") {
        ints[b] <- median(mixes[[i]][, c]) }
      if (cct == "mean") {
        ints[b] <- mean(mixes[[i]][, c]) }
      if (cct == "mode") {
        ints[b] <- hdr(mixes[[i]][, c], probs, h = bw.nrd0(mixes[[i]][, c]))$mode }      
      b <- b + 3
      c <- c + 1
      }
    rm(probs)
    cat(ints, "\n", file="fatsim-output.txt", sep="\t", append=TRUE)
  } 
 
return(mixes)                             
} #end of function  

# plotting function

densplot <- function (resdata, mtitle, probs = c(95, 75, 50), ct="median") {
    clr = gray((9:1)/10)
    scl = 1
    xspc = 0.5
    nbsources <- dim(resdata)[2]-2
    a <- length(resdata)
    
    groupseq <- seq(1, nbsources, by = 1)
    shift <- nbsources + 2
    usepars <- resdata[ , 1:nbsources]

        plot(1, 1, xlab = "Source", ylab = "Proportion", main = mtitle, xlim = c(min(groupseq) - 
            xspc, max(groupseq) + xspc), ylim = c(0, 1), type = "n", 
            xaxt = "n")
        axis(side = 1, at = min(groupseq):max(groupseq), 
          labels = colnames(resdata)[1:nbsources])

    clrs <- rep(clr, 5)
    for (j in 1:ncol(usepars)) {
        temp <- hdr(usepars[, j], probs, h = bw.nrd0(usepars[, j]))$hdr
        mod <- hdr(usepars[, j], probs, h = bw.nrd0(usepars[, j]))$mode
        line_widths <- seq(2, 20, by = 4) * scl
        bwd <- c(0.1, 0.15, 0.2, 0.25, 0.3) * scl
         for (k in 1:length(probs)) {
            temp2 <- temp[k, ]
            polygon(c(groupseq[j] - bwd[k], groupseq[j] - 
            bwd[k], groupseq[j] + bwd[k], groupseq[j] + 
            bwd[k]), c(max(min(temp2[!is.na(temp2)]), 0), 
            min(max(temp2[!is.na(temp2)]), 1), min(max(temp2[!is.na(temp2)]), 
            1), max(min(temp2[!is.na(temp2)]), 0)), col = clrs[k])
                
            if (ct == "mode") {points(j, mod, pch=16, col="white", cex=1)}
            if (ct == "mean") {points(j, mean(usepars[,j]), pch=16, col="white", cex=1)}
            if (ct == "median") {points(j, median(usepars[,j]), pch=16, col="white", cex=1)}

            }
        } }

# Function to plot the results for a certain modification

Fatsim.resplot <- function(resdat, nb, cct = "median", pair = TRUE){
a <- length(resdat)
mods <- resdat[[a-3]][nb + 1, ]
nbsources <- length(resdat[[a-1]][ ,2])
if (length(mods) > nbsources) consumer <- TRUE
if (length(mods) == nbsources) consumer <- FALSE

if (consumer==F){
  textwhich <- paste("sources modified by: ") 
  for (j in 1:(nbsources-1)) {
    textwhich <- paste(textwhich, colnames(resdat[[1]][j]), " ", mods[j], ", ", sep="")}
  textwhich <- paste(textwhich, colnames(resdat[[1]][nbsources]), " ", mods[nbsources], sep="")
    }  
  if (consumer==T & sum(mods)==0){
    textwhich <- paste("consumer modified by ", mods[1]) }
  if (consumer==T & sum(mods)!=0){
      textwhich <- paste("consumer modified by ", mods[1], ", ", "sources modified by: ", sep="") 
      for (j in 1:(nbsources-1)) {
        textwhich <- paste(textwhich, colnames(resdat[[1]][j]), " ", mods[j+1], ", ", sep="")}
      textwhich <- paste(textwhich, colnames(resdat[[1]][nbsources]), " ", mods[nbsources], sep="")
      }

if (pair == TRUE) {
  par(mfrow=c(1, 2))
  densplot(resdat[[1]], "Original data", ct=cct)
  densplot(resdat[[nb+1]], "Modified data", ct=cct)
  mtext(textwhich, side=4, cex = 0.8)
  } else {
  densplot(resdat[[nb+1]], "Modified data", ct=cct)
  mtext(textwhich, side=4, cex = 0.8) }
}

# Function to plot the data for a certain modification

Fatsim.datplot <- function(resdat, nb, pair = TRUE) {
  a <- length(resdat)
  xpred <- resdat[[a-2]][ ,1]
  ypred <- resdat[[a-2]][ ,2]
  x <- resdat[[a-1]][ ,2] 
  y <- resdat[[a-1]][ ,4]
  if (is.null(dim(resdat[[a]]))==F) {
    x <- x + resdat[[a]][ , 2]
    y <- y + resdat[[a]][ , 4] }
  xsd <- resdat[[a-1]][ ,3] 
  ysd <- resdat[[a-1]][ ,5]
  modn <- resdat[[a-5]]
  modnp <- resdat[[a-4]]
  mods <- resdat[[a-3]][nb + 1, ]
  nbsources <- length(x)

  if (length(mods) > nbsources) {
    consumer <- TRUE  # the cosumer has been modified
    xpred.mod <- xpred + mods[1]
    ypred.mod <- ypred + modnp
    x.mod <- x + mods[2:length(mods)]
    y.mod <- y + modn
    } else { # the consumer has not been modified
    consumer <- FALSE
    xpred.mod <- xpred
    ypred.mod <- ypred
    x.mod <- x + mods
    y.mod <- y + modn
    }

    library(graphics)
    
    xmax <- max(c(x, x.mod, xpred, xpred.mod)) + max(xsd) + 1
    xmin <- min(c(x, x.mod, xpred, xpred.mod)) - max(xsd) - 1
    ymax <- max(c(y, y.mod, ypred, ypred.mod)) + max(ysd) + 1
    ymin <- min(c(y, y.mod, ypred, ypred.mod)) - max(ysd) - 1

    symbcodes <- c(0:2, 4:7, 21:25, 8:18)
    symbs <- symbcodes[1:nbsources]

    if (consumer==F){
      textwhich <- paste("sources modified by: ") 
      for (j in 1:(nbsources-1)) {
        textwhich <- paste(textwhich, colnames(resdat[[1]][j]), " ", mods[j], ", ", sep="")}
      textwhich <- paste(textwhich, colnames(resdat[[1]][nbsources]), " ", mods[nbsources], sep="")}
    if (consumer==T & sum(mods)==0){
      textwhich <- paste("consumer modified by ", mods[1]) }
    if (consumer==T & sum(mods)!=0){
      textwhich <- paste("consumer modified by ", mods[1], ", ", "sources modified by: ", sep="") 
      for (j in 1:(nbsources-1)) {
        textwhich <- paste(textwhich, colnames(resdat[[1]][j]), " ", mods[j+1], ", ", sep="")}
      textwhich <- paste(textwhich, colnames(resdat[[1]][nbsources]), " ", mods[nbsources], sep="")
      }
      
  if (pair==TRUE) {
  par(mfrow=c(1, 2))
    # plot original data
    plot (xpred, ypred, las=1, mgp=c(2,0.3,0), tcl=-0.2, pch=16, col="black",
    xlab=expression(paste(delta^13,"C (â)")), ylab=expression(paste(delta^15,"N (â)")), 
    xlim=c(xmin, xmax),ylim=c(ymin, ymax), 
    main=c("Original data" ))  
    points(x, y, pch=symbs)
    arrows(x, y-ysd, x, y+ysd, code=3, angle=90, length=0.1)
    arrows(x-xsd, y, x+xsd, y, code=3, angle=90, length=0.1) 

    #plot modified data
    plot (xpred.mod, ypred.mod, las=1, mgp=c(2,0.3,0), tcl=-0.2, pch=16, col="black",
    xlab=expression(paste(delta^13,"C (â)")), ylab=expression(paste(delta^15,"N (â)")), 
    xlim=c(xmin, xmax),ylim=c(ymin, ymax), 
    main=c("Modified data" ))  
    points(x.mod, y.mod, pch=symbs)
    arrows(x.mod, y.mod-ysd, x.mod, y.mod+ysd, code=3, angle=90, length=0.1)
    arrows(x.mod-xsd, y.mod, x.mod+xsd, y.mod, code=3, angle=90, length=0.1) 
    mtext(textwhich, side = 4, cex =0.8)

    legend(locator(1), c("consumer", colnames(resdat[[1]][1:nbsources])), pch=c(16,symbs),cex=0.8)
  }  

  if (pair==FALSE) {
    #plot modified data
    plot (xpred.mod, ypred.mod, las=1, mgp=c(2,0.3,0), tcl=-0.2, pch=16, col="black",
    xlab=expression(paste(delta^13,"C (â)")), ylab=expression(paste(delta^15,"N (â)")), 
    xlim=c(xmin, xmax),ylim=c(ymin, ymax), 
    main=c("Modified data" ))  
    points(x.mod, y.mod, pch=symbs)
    arrows(x.mod, y.mod-ysd, x.mod, y.mod+ysd, code=3, angle=90, length=0.1)
    arrows(x.mod-xsd, y.mod, x.mod+xsd, y.mod, code=3, angle=90, length=0.1) 

    legend(locator(1), c("consumer", colnames(resdat[[1]][1:nbsources])), pch=c(16,symbs),cex=0.8)
    mtext(textwhich, side = 4, cex =0.8)
  } } 

# example datasets
# consumer:

x <- c(-20.01, 10.01)
y <- c(-20.51, 10.51)
z <- c(-20.51, 9.81)
w <- c(-19.61, 10.31)
qu <- c(-19.81, 10.11)
r <- c(-19.91, 9.81)
sa <- c(-20.21, 9.91)
ta <- c(-19.91, 9.71)
predator <- rbind(x, y, z, w, qu, r, sa, ta)
colnames(predator) <- c("dC13", "dN15")
rownames(predator) <- NULL
rm(x, y, z, w, qu, r, sa, ta)

# Sources
a <- c(-21.0, 0.5, 11.0, 0.5)
b <- c(-15.0, 0.5, 15.0, 0.5)
cc <- c(-25.0, 0.5, 5.0, 0.5)
d <- c(-19.0, 0.5, 9.0, 0.5)
sour <- as.data.frame(rbind(a, b, cc, d))
rm(a, b, cc, d)
so <- c("bread", "cheese", "chocolate", "potato")
sour <- as.data.frame(cbind(so, sour))
names(sour) <- c("source", "d13Cmean", "d13Csd", "d15Nmean", "d15Nsd") 
rm(so)

# References
# Parnell, A.C., Inger, R., Bearhop, S. & Jackson, A.L. (2010) Source partitioning using stable isotopes: coping with too much variation. PLoS ONE, 5, e9672.