Age_OSLC14.R

#' Bayesian analysis for age estimation of OSL measerments and C-14 ages of various samples
#'
#' This function compute an age of OSL data of at least two samples and calibrate 14C ages of samples to get an age (in ka).\cr
#' Age of OSL data are computed according to the model given in Combes and Philippe (2017).
#' Single-grain or Multi-grain OSL measurements can be analysed simultaneouly (with output of \code{\link{Generate_DataFile}}
#' or \code{\link{Generate_DataFile_MG}} or both of them using \code{\link{combine_DataFiles}}).
#' Samples, for which data is avalilable in several BIN files, can be analysed.\cr
#' For C14 data, the user can choose one of the following radiocarbon calibration curve:
#' Northern or Sourthen Hemisphere or marine atmospheric.
#'
#' @param DATA list of objects: LT, sLT, ITimes, dLab, ddot_env, regDose, J, K, Nb_measurement,
#' provided by the function \code{\link{Generate_DataFile}} or \code{\link{Generate_DataFile_MG}} or \code{\link{combine_DataFiles}}.
#' \code{DATA} contains information for more than one sample.
#' If there is stratigraphic relations between samples, informations in DATA must be ordered by order of incresing ages.
#' See the details section to for more informations.
#' @param Data_C14Cal numeric vector: corresponding to 14C age estimate (in years, conversion in ka is automatically donne in the function).
#' If there is stratigraphic relations between samples, \code{Data_C14Cal} must be ordered by order of incresing ages.
#' @param Data_SigmaC14Cal numeric vector: correponding to the error of 14C age estimates.
#' @param Nb_sample integer: number of samples (OSL data and 14C age),
#' (\code{Nb_sample>3}, at least to sample of OSL data and one sample of 14C age).
#' @param SampleNames character vector: sample names for both OSL data and C14 data. The length of this vector is equal to \code{Nb_sample}.
#' If there is stratigephic relation, this vector must be ordered by increasing order (to mix OSL samples and 14C ages if it is needed).
#' @param SampleNature numeric matrix: states the nature of the sample.
#' Row number of \code{SampleNature} matrix is equal to \code{2} and column number is equal to \code{Nb_sample}.
#' First line of the matrix: \code{SampleNature[1,i]} states if sample whose number ID is equal to \code{i}, is an OSL sample \bold{1} or not \bold{0}.
#' Second line of the matrix: \code{SampleNature[2,i]} states if sample whose number ID is equal to \code{i}, is an 14C sample \bold{1} or not \bold{0}.
#' @param PriorAge numeric vector (with default): lower and upper bounds for age parameter of each sample \bold{(in ka)}.
#'  Note that, \code{length(PriorAge)=2*Nb_sample}
#'  and \code{PriorAge[2i-1,2i]} correponds to the lower and upper bounds of sample whose number ID is equal to \code{i}.
#' @param SavePdf boolean (with default): if TRUE save graphs in pdf file named \code{OutputFileName} in folder \code{OutputFilePath}.
#' @param OutputFileName character (with default): name of the pdf file that will be generated by the function if \code{SavePdf}=TRUE,
#' \code{length(OutputFileName)=3}, see \bold{PLOT OUTPUT} in \bold{Value} section for more informations.
#' @param OutputFilePath character (with default): path to the pdf file that will be generated by the function if \code{SavePdf}=TRUE.
#' If it is not equal to "", it must be terminated by "/".
#' @param SaveEstimates boolean (with default): if TRUE save Bayes estimates, credible interval at level 68% and 95%  and
#' the result of the gelman en Rubin test of convergency, in a csv table named \code{OutputFileName} in folder \code{OutputFilePath}.
#' @param OutputTableName character (with default): name of the table that will be generated by the function if \code{SaveEstimates}=TRUE.
#' @param OutputTablePath character (with default): path to the table that will be generated by the function if \code{SaveEstimates}=TRUE.
#' If it is not equal to "", it must be terminated by "/".
#' @param StratiConstraints numeric matrix or character(with default): input object for the statigraphic relation between samples.
#' If there is stratigraphic relation between samples see the details section for instructions regarding how to correctly fill \code{StratiConstraints},
#' the user can refer to a matrix (numeric matrix) or to a csv file (character).
#' Otherwise, default value is suitable.
#' @param sepSC character (with default): if \code{StratiConstraints} is character,
#' indicate column separator in \code{StratiConstraints} csv file.
#' @param BinPerSample integer vector (with default): vector with the number of BIN files per OSL sample.
#' The length of this vector is equal to the number of OSL samples.
#' `BinPerSample[i]` correponds to the number of BIN files for the sample whose number ID is equal to \code{i}.
#' For more information to fill this vector, we refer to detatils in \code{\link{Generate_DataFile}} or in \code{\link{Generate_DataFile_MG}}.
#' @param THETA numeric matrix or character (with default): input object for systematic and individual error for OSL samples.
#' If systematic errors are considered, see the details section for instructions regarding how to correctly fill \code{THETA};
#' the user can refer to a matrix (numeric matrix) or to a csv file (character).
#' Otherwise, default value is suitable, and only individual error is considered.
#' @param sepTHETA character (with default): if \code{THETA} is character,
#' indicate column separator in \code{THETA} csv file.
#' @param LIN_fit logical (with default): if TRUE (default) allows a linear component,
#' on top of the (default) saturating exponential curve, for the fitting of dose response curves, for OSL samples.
#' See details for more informations on the proposed dose response curves.
#' @param Origin_fit logical (with default): if TRUE, forces the dose response curves to pass through the origin.
#' See details for more informations on the proposed growth curves, for OSL samples.
#' @param distribution character (with default): type of distribution that defines
#' how individual equivalent dose values are distributed around the palaeodose, for OSL samples.
#' Allowed inputs are \bold{"cauchy"}, \bold{"gaussian"}, \bold{"lognormal_A"} and \bold{"lognormal_M"}, see details for more informations.
#' @param Model_C14 character (with default): if \bold{"full"}, error on estimate calibration curve is taken account, for 14C samples.
#' If \bold{"naive"} this error is not taken account in the age estimate.
#' @param CalibrationCurve character (with default): calibration curve choosen, for 14C samples. Allowed inputs are
#' \itemize{
#'   \item \bold{"AtmosphericNorth"} for Northern Hemisphere atmospheric radiocarbon calibration curve,
#'   \item \bold{"Marine"} for Marine radiocarbon calibration curve,
#'   \item \bold{"AtmosphericSouth"} for Southern Hemisphere atmospheric radiocarbon calibration curve,
#'   \item \bold{a csv file, with tree columns, the first column is dedicated to "Cal.BP", the second to "X14C.age", the third to "Error".
#'   The decimal of this file must be a dot, and the separator must be a comma. }
#' }
#' @param Iter integer (with default): number of iterations for the MCMC computation (for more information see \code{\link{jags.model}}).
#' @param t integer (with default): 1 every \code{t} iterations of the MCMC is considered for sampling the posterior distribution
#' (for more information see \code{\link{jags.model}}).
#' @param n.chains integer (with default): number of independent chains for the model (for more information see \code{\link{jags.model}}).
#'
#' @param quiet \code{\link{logical}} (with default): enables/disables \link{rjags} messages
#'
#' @details
#'
#' Note that there is tree type of arguments in the previous list.
#' There are arguments for informtations concerning only OSL samples: \code{DATA}, \code{BinPerSample}, \code{THETA},
#' \code{sepTHETA}, \code{LIN_fit}, \code{Origin_fit}, \code{distribution}.
#'
#' There are arguments for informtations concerning only C14 samples: \code{Data_C14Cal}, \code{Data_SigmaC14Cal},
#' \code{Model_C14}, \code{CalibrationCurve}.
#'
#' There are arguments for informtations concerning all the samples: \code{Nb_sample}, \code{SampleNames}, \code{SampleNature},
#' \code{PriorAge}, \code{SavePdf}, \code{OutputFileName}, \code{OutputFilePath}, \code{SaveEstimates}, \code{OutputTableName},
#' \code{OutputTablePath}, \code{StratiConstraints}, \code{sepSC}.\cr
#'
#' \bold{** How to fill} \code{StratiConstraints} \bold{? **}\cr
#'
#' If there is stratigraphic relations between samples, \bold{14C estimate age in \code{Data_C14Cal} must be ordered by order of increasing ages,
#' as informations in \code{DATA}}. Names in \code{SampleNames} must be ordered and correponds to the order in \code{Data_C14Cal} and in \code{DATA},
#' also if it is needed to mix names of OSL samples and 14C samples.
#'
#' The user can fill the \code{StratiConstraints} matrix as follow.
#' \enumerate{
#'  \item \bold{Size of the matrix}: row number of \code{StratiConstraints} matrix is equal to \code{Nb_sample+1},
#' and column number is equal to \code{Nb_sample}.
#'  \item \bold{First line of the matrix}:
#' for all \code{i in {1,...,Nb_Sample}}, \code{StratiConstraints[1,i]=1} that means the lower bound of the sample age (given in \code{PriorAge[2i-1]})
#' for the sample whose number ID is equal to \code{i}, is taken into account.
#'  \item \bold{Sample relations}: for all  \code{j in {2,...,Nb_Sample+1}} and all \code{i in {j,...,Nb_Sample}},
#' \code{StratiConstraints[j,i]=1} if sample age whose number ID is equal to \code{j-1} is lower than sample age whose number ID is equal to \code{i}.
#' Otherwise, \code{StratiConstraints[j,i]=0}.
#' }
#' Note that \code{StratiConstraints_{2:Nb_sample+1,1:Nb_sample}} is a upper triangular matrix.
#'
#' The user can also use \code{\link{SCMatrix}}  or \code{\link{SC_Ordered}} (if all samples are ordered) function to construc
#' the \code{StratiConstraints} matrix.
#'
#' The user can also refer to a csv file that containts the relation between samples as defined above.
#' The user must take care about the separator used in the csv file using the argument \code{sepSC}.\cr
#'
#'  \bold{** How to fill} \code{THETA} \bold{covariance matrix concerning common and individual error? **}\cr
#'
#' If systematic errors are considered, the user can fill the \code{THETA} matrix as follow.
#' \itemize{
#'  \item row number of \code{THETA} is equal the column number, equal to \code{Nb_sample}.
#'  \item For all \code{i in {1,...,Nb_sample}}, \code{THETA[i,i]} containts individual error
#'  plus systematic error of the sample whose number ID is equal to \code{i}.
#'  \item For all \code{i,j in {1,...,Nb_sample}} and \code{i} different from \code{j} ,
#' \code{THETA[i,j]} containts common error between samples whose number ID are equal to \code{i} and \code{j}.
#' }
#' Note that \code{THETA[i,j]} is a symetric matrix.
#'
#' The user can also refer to a .csv file that containts the errors as defined above.\cr
#'
#' \bold{** Option on growth curves **}\cr
#'
#' As for \code{\link{Age_Computation}} and \code{\link{Palaeodose_Computation}}, the user can choose from 4 dose response curves:
#' \itemize{
#'   \item \bold{Saturating exponential plus linear growth} (\code{AgesMultiCS2_EXPLIN}):
#'
#'   for all \code{x} in IR+, \code{f(x)=a(1-exp(-x/b))+cx+d}; select
#'   \itemize{
#'     \item \code{LIN_fit=TRUE}
#'     \item \code{Origin_fit=FALSE}
#'   }
#'   \item \bold{Saturating exponential growth} (\code{AgesMultiCS2_EXP}):
#'
#'   for all \code{x} in IR+, \code{f(x)=a(1-exp(-x/b))+d}; select
#'   \itemize{
#'     \item \code{LIN_fit=FALSE}
#'     \item \code{Origin_fit=FALSE}
#'   }
#'   \item \bold{Saturating exponential plus linear growth and fitting through the origin} (\code{AgesMultiCS2_EXPLINZO}):
#'
#'   for all \code{x} in IR+, \code{f(x)=a(1-exp(-x/b))+cx}; select
#'   \itemize{
#'     \item \code{LIN_fit=TRUE}
#'     \item \code{Origin_fit=TRUE}
#'   }
#'   \item \bold{Saturating exponential growth and fitting through the origin} (\code{AgesMultiCS2_EXPZO}):
#'
#'   for all \code{x} in IR+, \code{f(x)=a(1-exp(-x/b))}; select
#'   \itemize{
#'     \item \code{LIN_fit=FALSE}
#'     \item \code{Origin_fit=TRUE}
#'   }
#' }
#'
#' \bold{** Option on equivalent dose distribution around the palaeodose **}\cr
#'
#' The use can choose between :
#' \itemize{
#'   \item \code{cauchy}: a Cauchy distribution with location parameter equal to the palaeodose of the sample
#'   \item \code{gaussian}: a Gaussian distribution with mean equal to the palaeodose of the sample
#'   \item \code{lognormal_A}: a log-normal distribution with mean or \bold{A}verage equal to the palaeodose of the sample
#'   \item \code{lognormal_M}: a log-normal distribution with \bold{M}edian equal to the palaeodose of the sample
#' }
#'
#' \bold{** More precision on \code{Model} **}\cr
#'
#' We propose two models "full" or "naive". If \code{Model='full'} that means measurement error and error on calibration curve are taken account in
#' the Bayesian model; if \code{Model="naive"} that means only error on measurement are taken account in the mode.
#'
#' More precisely, the model considered here, as the one developped by Christen, JA (1994), assume multiplicative effect of errors to address the
#' problem of outliers. In addition, to not penalyse variables that are not outliers and damage theirs estimation,
#' we introduce a structure of mixture, that means only variable that are considered as outlier have in addition a multiplicative error.
#'
#' @return
#' \bold{NUMERICAL OUTPUT}\cr
#'
#' \enumerate{
#' \item \bold{A list containing the following objects:}
#'  \itemize{
#'   \item \bold{Sampling}: that corresponds to a sample of the posterior distributions of the age parameters (in ka for both C14 samples and OSL samples);
#'   \item \bold{PriorAge}: stating the priors used for the age parameter;
#'   \item \bold{StratiConstraints}: stating the stratigraphic relations between samples considered in the model;
#'   \item \bold{Model_OSL_GrowthCurve}: stating which dose response fitting option was chosen;
#'   \item \bold{Model_OSL_Distribution}: stating which distribution was chosen to model the dispersion of
#'  individual equivalent dose values around the palaeodose of the sample;
#'   \item \bold{Model_C14}: stating which model was chosen (\code{"full"} or \code{"naive"});
#'   \item \bold{CalibrationCurve}: stating which radiocarbon calibration curve was chosen;
#'   \item \bold{Outlier}: stating the names of samples that must be outliers.
#'  }
#'
#'   \item \bold{The Gelman and Rubin test of convergency}: prints the result of the Gelman and Rubin test of convergency for the age estimate for each sample.
#' A result close to one is expected.\cr
#' In addition, the user must visually assess the convergency of the trajectories by looking at the graph
#' generated by the function (see \bold{PLOT OUTPUT} for more informations).\cr
#' If both convergencies (Gelman and Rubin test and plot checking) are satisfactory,
#' the user can consider the estimates as valid.
#' Otherwise, the user may try increasing the number of MCMC interations (\code{Iter})
#' or be more precise on the \code{PriorAge} parameter to reach convergency.
#'   \item \bold{Credible intervals and Bayes estimates}: prints the Bayes esitmates, the credible intervals at 95% and 68% for
#' the age parameters for each sample.
#' }
#'
#' \bold{PLOT OUTPUT}
#'
#' \enumerate{
#'  \item\bold{MCMC trajectories}: A graph with the MCMC trajectories and posterior distributions of the age parameter is displayed. \cr
#' On each line, the plot on the left represents the MCMC trajectories, and the one on the right the posterior distribution of the parameter.
#'  \item \bold{Age estimate and HPD at 95% of 14C samples on calibration curve}: plot age estimate and HPD on calibration plot.
#'  \item \bold{Summary of sample age estimates}: plot credible intervals and Bayes estimate of each sample age on a same graph.
#' }
#'
#' @author Claire Christophe, Anne Philippe, Guillaume Guerin, Sebastian Kreutzer
#'
#' @note Please note that the initial values for all MCMC are currently all the same for all chains since we rely on the automatic
#' initial value generation of JAGS. This is not optimal and will be changed in future. However, it does not affect the quality
#' of the age estimates if the chains have converged.
#'
#' @seealso
#' [rjags], [plot_MCMC], [SCMatrix], [plot_Ages]
#'
#' @references
#' Reimer PJ, Bard E, Bayliss A, Beck JW, Blackwell PC, Bronl Ramsey C, Buck CE, Cheng H, Edwards RL, Friedrich M,
#' Grootes PM, Guilderson TP, Haflidason H, Hajdas I, Hatte C, Heaton TJ, Hoffmann DL, Hogg AG, Hughen KA, Kaiser KF, Kromer B,
#' Manning SW, Niu M, Reimer RW, Richards DA, Scott EM, Southon JR, Staff RA, Turney CSM, van der Plicht J. 2013.
#' IntCal13 ans Marine13 radiocarbon age calibration curves 0-50000 years cal BP. Radiocarbon 55(4)=1869-1887.
#'
#' Hogg AG, Hua Q, Blackwell PG, Niu M, Buck CE, Guilderson TP, Heaton TJ, Palmer JG, Reimer PJ, Reimer RW, Turney CSM, Zimmerman SRH.
#' 2013. SHCal13 Southern Hemisphere calibration, 0-50000 years cal BP. Radiocarbon 55(4):1889-1903
#'
#'
#' @examples
#' ## Load data
#' # OSL data
#' data(DATA1,envir = environment())
#' data(DATA2,envir = environment())
#' Data <- combine_DataFiles(DATA2,DATA1)
#'
#' # 14C data
#  data(DATA_C14,envir = environment())
#' C14Cal <- DATA_C14$C14[1,1]
#' SigmaC14Cal <- DATA_C14$C14[1,2]
#' Names <- DATA_C14$Names[1]
#'
#' # Prior Age
#' prior=rep(c(1,60),3)
#' samplenature=matrix(data=c(1,0,1,0,1,0),ncol=3,nrow=2,byrow=TRUE)
#' SC <- matrix(data=c(1,1,1,0,1,1,0,0,1,0,0,0),ncol=3,nrow=4,byrow=TRUE)
#'
#' ## Age computation of samples
#' Age <- Age_OSLC14(DATA=Data,Data_C14Cal=C14Cal,Data_SigmaC14Cal=SigmaC14Cal,
#'    SampleNames=c("GDB5",Names,"GDB3"),Nb_sample=3,SampleNature=samplenature,
#'    PriorAge=prior,StratiConstraints=SC,Iter=50,n.chains=2)
#' @md
#' @export
Age_OSLC14 <- function(
  DATA,
  Data_C14Cal,
  Data_SigmaC14Cal,
  Nb_sample,
  SampleNames,
  SampleNature,
  PriorAge = rep(c(10, 60), Nb_sample),
  SavePdf = FALSE,
  OutputFileName = c('MCMCplot', 'HPD_Cal14CCurve', "summary"),
  OutputFilePath = c(""),
  SaveEstimates = FALSE,
  OutputTableName = c("DATA"),
  OutputTablePath = c(''),
  StratiConstraints = c(),
  sepSC = c(','),
  BinPerSample = rep(1, sum(SampleNature[1, ])),
  THETA = c(),
  sepTHETA = c(','),
  LIN_fit = TRUE,
  Origin_fit = FALSE,
  distribution = c("cauchy"),
  Model_C14 = c("full"),
  CalibrationCurve = c("AtmosphericNorth"),
  Iter = 50000,
  t = 5,
  n.chains = 3,
  quiet = FALSE) {


  #--- StratiConstraints matrix
  if(length(StratiConstraints)==0){
    StratiConstraints=matrix(data=c(rep(1,Nb_sample),rep(0,Nb_sample*Nb_sample)),ncol=Nb_sample,nrow = (Nb_sample+1),byrow = T)
  }else{
    if(is(StratiConstraints)[1]=="character"){
      SCMatrix=read.csv(StratiConstraints,sep=sepSC)
      StratiConstraints=as.matrix(SCMatrix)
    }
  }

  #--- Calibration curve
  TableauCalib=c()
  if(CalibrationCurve=="AtmosphericNorth"){
    AtmosphericNorth_CalC14<-0
    data(AtmosphericNorth_CalC14,envir = environment())
    TableauCalib=AtmosphericNorth_CalC14
  }else{if(CalibrationCurve=="AtmosphericSouth"){
    AtmosphericSouth_CalC14<-0
    data(AtmosphericSouth_CalC14,envir = environment())
    TableauCalib=AtmosphericSouth_CalC14
  }else{if(CalibrationCurve=="Marine"){
    Marine_CalC14<-0
    data(Marine_CalC14,envir = environment())
    TableauCalib=Marine_CalC14
  }else{
    TableauCalib=read.csv(file=CalibrationCurve,sep=",",dec=".")
  }}}
  AgeBP=rev(TableauCalib[,1])
  CalC14=rev(TableauCalib[,2])
  SigmaCalC14=rev(TableauCalib[,3])

  # #--- C14 prepration: Calibration curve
  # TableauCalib=read.csv(file=paste("inst/extdata/",CalibrationCurve,"_CalC14.csv",sep=""),sep=",",dec=".")
  # AgeBP=rev(TableauCalib[,1])/1000
  # CalC14=rev(TableauCalib[,2])
  # SigmaCalC14=rev(TableauCalib[,3])

  #--- OSL preparation
  #- Index preparation
  CSBinPerSample=cumsum(BinPerSample)
  LengthSample=c()
  for(ns in 1:sum(SampleNature[1,])){
    LengthSample=c(LengthSample,length(DATA$LT[[ns]][,1]))
  }
  CSLengthSample=c()
  CSLengthSample=c(0,cumsum(LengthSample))
  index2=c(0,cumsum(DATA$J))
  #- File preparation
  LT=matrix(data=0,nrow=sum(DATA$J),ncol=(max(DATA$K)+1))
  sLT=matrix(data=0,nrow=sum(DATA$J),ncol=(max(DATA$K)+1))
  IrrT=matrix(data=0,nrow=sum(DATA$J),ncol=(max(DATA$K)))
  for(ns in 1:sum(SampleNature[1,])){
    LT[seq(CSLengthSample[ns]+1,CSLengthSample[ns+1],1),1:length(DATA$LT[[ns]][1,])]<-DATA$LT[[ns]]
    sLT[seq(CSLengthSample[ns]+1,CSLengthSample[ns+1],1),1:length(DATA$sLT[[ns]][1,])]<-DATA$sLT[[ns]]
    IrrT[seq(CSLengthSample[ns]+1,CSLengthSample[ns+1],1),1:length(DATA$ITimes[[ns]][1,])]<-DATA$ITimes[[ns]]
  }
  #- THETA matrix
  if(length(THETA[,1])==0){
    THETA=diag(DATA$ddot_env[2,CSBinPerSample]+(DATA$ddot_env[1,CSBinPerSample])^2*DATA$dLab[2,CSBinPerSample])
  }else{
    if(is(THETA)[1]=="character"){
      errorMatrix=read.csv(THETA,sep=sepTHETA)
      THETA=as.matrix(errorMatrix)
    }
  }

  #---  Index preparation
  ind_OSL=which(SampleNature[1,]==1)
  CS_OSL=cumsum(SampleNature[1,])
  ind_C14=which(SampleNature[2,]==1)
  CS_C14=cumsum(SampleNature[2,])
  ind_change=c(1)
  for(i in 2:(Nb_sample-1)){
    if(SampleNature[1,i]!=SampleNature[1,i+1]){
      ind_change=c(ind_change,i)
    }
  }

  ind_change=c(ind_change,Nb_sample)
  q=length(ind_change)%/%2
  r=length(ind_change)%%2

  ##--- description du model BUG
  BUGModel=c()
  #- Prior
  ModelPrior<-0
  data(ModelPrior,envir = environment())
  BUGPrior=c()

  if(r==1){
    if(SampleNature[1,1]==1){
      BUGPrior=paste(BUGPrior,ModelPrior$Sample1_OSL)
    }else{
      BUGPrior=paste(BUGPrior,ModelPrior$Sample1_C14)
    }
    if(SampleNature[1,2]==1){
      BUGPrior=paste(BUGPrior,ModelPrior$OSL_C14)
    }else{
      BUGPrior=paste(BUGPrior,ModelPrior$C14_OSL)
    }
  }else{
    q=q-1
    if(SampleNature[1,1]==1){
      BUGPrior=paste(BUGPrior,ModelPrior$Sample1_OSL)
    }else{
      BUGPrior=paste(BUGPrior,ModelPrior$Sample1_C14)
    }
    if(SampleNature[1,2]==1){
      BUGPrior=paste(BUGPrior,ModelPrior$OSL_C14)
    }else{
      BUGPrior=paste(BUGPrior,ModelPrior$C14_OSL)
    }
    if(SampleNature[1,Nb_sample]==1){
      BUGPrior=paste(BUGPrior,ModelPrior$OSL)
    }else{BUGPrior=paste(BUGPrior,ModelPrior$C14)}
  }

  #- partie C14
  ModelC14<-0
  data(ModelC14,envir = environment())
  if(Model_C14=="full"){
    BUGModel=paste(ModelC14$full,BUGPrior)
  }else{
    BUGModel=paste(ModelC14$naive,BUGPrior)
  }
  #- partie OSL
  ModelOSL<-0
  data(ModelOSL,envir = environment())
  if(LIN_fit==TRUE){
    cLIN=c('LIN')
  }else{cLIN=c()}
  if(Origin_fit==TRUE){
    cO=c("ZO")
  }else{cO=c()}
  Model_GrowthCurve=c(paste("AgesMultiOSL_EXP",cLIN,cO,sep=""))
  BUGModel=c(paste("model{",ModelOSL[[Model_GrowthCurve]][[distribution]],BUGModel,"}"))

  if(Model_C14=="full"){
    dataList = list('q'=q,"ind_change"=ind_change,"ind_OSL"=ind_OSL,"ind_C14"=ind_C14,"CS_OSL"=CS_OSL,"CS_C14"=CS_C14,
                    'X'=Data_C14Cal,"sigma"=Data_SigmaC14Cal,
                    "xTableauCalib"=AgeBP,"yTableauCalib"=CalC14,"zTableauCalib"=SigmaCalC14,
                    'N'= LT,'sN'=sLT,"IT"=IrrT,
                    "sDlab"=DATA$dLab[1,],
                    'J'=DATA$J,
                    'K'=DATA$K,
                    "ddot"=DATA$ddot_env[1,CSBinPerSample],
                    "Gamma"=THETA,
                    "index"=index2,
                    "BinPerSample"=BinPerSample,
                    "CSBinPerSample"=CSBinPerSample,
                    "xbound"=PriorAge,"StratiConstraints"=StratiConstraints)
  }else{
    dataList = list('q'=q,"ind_change"=ind_change,"ind_OSL"=ind_OSL,"ind_C14"=ind_C14,"CS_OSL"=CS_OSL,"CS_C14"=CS_C14,
                    'X'=Data_C14Cal,"sigma"=Data_SigmaC14Cal,
                    "xTableauCalib"=AgeBP,"yTableauCalib"=CalC14,
                    'N'= LT,'sN'=sLT,"IT"=IrrT,
                    "sDlab"=DATA$dLab[1,],
                    'J'=DATA$J,
                    'K'=DATA$K,
                    "ddot"=DATA$ddot_env[1,CSBinPerSample],
                    "Gamma"=THETA,
                    "index"=index2,
                    "BinPerSample"=BinPerSample,
                    "CSBinPerSample"=CSBinPerSample,
                    "xbound"=PriorAge,"StratiConstraints"=StratiConstraints)
  }

  ##open text connection
  con <-  textConnection(BUGModel)

  jags <-
    rjags::jags.model(
      file = con,
      data = dataList,
      n.chains = n.chains,
      n.adapt = Iter,
      quiet = quiet
    )

  ##close connection
  close(con)

  ##set progress.bar
  if(quiet) progress.bar <- 'none' else progress.bar <- 'text'

  update(jags,Iter)
  echantillon <-
    rjags::coda.samples(jags,
                        c("A", 'Z'),
                        min(Iter, 10000),
                        thin = t,
                        progress.bar = progress.bar)
  U <- summary(echantillon)

  Sample=echantillon[[1]]
  for(i in 2:n.chains){
    Sample=rbind(Sample,echantillon[[i]])
  }

  nom=c()
  for(i in 1:Nb_sample){
    nom=c(nom,paste("A_",SampleNames[i],sep=""))
  }

  ##plot MCMC
  if(SavePdf){
    pdf(file=paste(OutputFilePath,OutputFileName[1],'.pdf',sep=""))
  }

  plot_MCMC(echantillon, sample_names = SampleNames)

  if(SavePdf){
    dev.off()
  }

  Outlier <- SampleNames[ind_C14[which(U$statistics[(Nb_sample+1):(Nb_sample+sum(SampleNature[2,])),1]<1.5)]]

  ##- Gelman and Rubin test of convergency of the MCMC
  CV=gelman.diag(echantillon,multivariate=FALSE)
  cat(paste("\n\n>> Convergence of MCMC for the age parameters <<\n"))
  cat("----------------------------------------------\n")
  cat(paste("Sample name ", " Bayes estimate ", " Uppers credible interval\n"))
  for(i in 1:Nb_sample){
    #cat(paste(" Sample name: ", SampleNames[i],"\n"))
    #cat("---------------------\n")
    cat(paste(paste("A_",SampleNames[i],sep=""),"\t",round(CV$psrf[i,1],2),"\t\t",round(CV$psrf[i,2],2),"\n"))
  }



  cat("\n\n________________________________________________________________________________\n")
  cat(" *** WARNING: following informations are only valid if MCMC chains converged  ***\n")
  cat("________________________________________________________________________________\n")

  # Matrix of results
  rnames=c()
  for(i in 1:Nb_sample){
    rnames=c(rnames,paste("A_",SampleNames[i],sep=""))
  }
  R=matrix(data=NA,ncol=8,nrow=Nb_sample,
           dimnames=list(rnames,c("lower bound at 95%","lower bound at 68%","Bayes estimate",
                                  "upper bound at 68%","upper bound at 95%","",
                                  "Convergencies: Bayes estimate","Convergencies: uppers credible interval")))

  ##- Bayes estimate and credible interval
  cat(paste("\n\n>> Bayes estimates of Age for each sample and credible interval <<\n"))
  AgePlot95=matrix(data=NA,nrow=Nb_sample,ncol=3)
  AgePlot68=matrix(data=NA,nrow=Nb_sample,ncol=3)
  AgePlotMoy=rep(0,Nb_sample)
  for(i in 1:Nb_sample){
    cat("------------------------------------------------------\n")
    #cat(paste(" Sample name: ", SampleNames[i],"\n"))
    #cat("---------------------\n")

    cat(paste("Sample name", "\t","Bayes estimate"," Credible interval: \n"))
    cat(paste(paste("A_",SampleNames[i],sep=""),"\t",round(mean(Sample[,i]),3),'\n'))
    cat("\t\t\t\t\t\t lower bound \t upper bound\n")
    HPD_95=ArchaeoPhases::CredibleInterval(Sample[,i],0.95)
    HPD_68=ArchaeoPhases::CredibleInterval(Sample[,i],0.68)
    cat("\t\t\t\t at level 95% \t",round(c(HPD_95[2]),2),"\t\t",round(c(HPD_95[3]),2),"\n")
    cat("\t\t\t\t at level 68% \t",round(c(HPD_68[2]),2),"\t\t",round(c(HPD_68[3]),2),"\n")
    AgePlot95[i,]=HPD_95
    AgePlot68[i,]=HPD_68
    AgePlotMoy[i]=round(mean(Sample[,i]),3)

    R[i,3]=round(mean(Sample[,i]),3)
    R[i,c(1,5)]=round(HPD_95[2:3],3)
    R[i,c(2,4)]=round(HPD_68[2:3],3)
    R[i,6]=c('')
    R[i,7]=round(CV$psrf[i,1],2)
    R[i,8]=round(CV$psrf[i,2],2)

  }

  cat("\n------------------------------------------------------\n")

  # Representation graphique des resultats
  #        des HPD sur la courbe de calibration
  if(sum(SampleNature[2,])>1){
    couleur=rainbow(Nb_sample)
    par(mfrow=c(1,1),las = 0,mar=c(5,5,2,2))
    xl=c(min(PriorAge[seq(1,(2*Nb_sample-1),2)]),max(PriorAge[seq(2,(2*Nb_sample),2)]))
    plot(xl,xl,col="white",xlab=c("Age (in ka)"),ylab=c("cal C14"),xaxt="n",yaxt="n",cex.lab=1.8)
    axis(2,cex.axis=2)
    axis(1,cex.axis=2)
    polygon(c(AgeBP,rev(AgeBP)),c(CalC14+2*SigmaCalC14,rev(CalC14-2*SigmaCalC14)),col="gray",border="black")
    for(i in ind_C14){
      lines(c(AgePlot95[i,2:3]),rep(Data_C14Cal[CS_C14[i]],2),col=couleur[i],lwd=4)
      lines(AgePlotMoy[i],Data_C14Cal[CS_C14[i]],col="black",lwd=2,type='p')
    }
    legend("topleft",SampleNames[ind_C14],lty=rep(1,Nb_sample),lwd=rep(2,Nb_sample),cex=1,col=couleur[ind_C14])
    if(SavePdf==TRUE){
      dev.print(pdf,file=paste(OutputFilePath,OutputFileName[2],'.pdf',sep=""),width=8,height=10)
    }
  }


  ##CSV output
  if(SaveEstimates){
    df_output <-
      data.frame(
        SampleNames = SampleNames,
        Bayes_estimate = CV$psrf[1:length(SampleNames), 1],
        Upper_Credibility_Interval = CV$psrf[1:length(SampleNames), 2],
        AGE = AgePlotMoy,
        HPD68.MIN = AgePlot68[,2],
        HPD68.MAX = AgePlot68[,3],
        HPD95.MIN = AgePlot95[,2],
        HPD95.MAX = AgePlot95[,3],
        stringsAsFactors = FALSE
      )

    write.table(x = df_output,
                file = paste0(OutputTablePath, OutputTableName, ".csv"),
                sep = ";",
                row.names = FALSE,
                append = FALSE
    )

  }

  # Create return objecty -------------------------------------------------------------------------
  output <- .list_BayLum(
    "Ages" = data.frame(
      SAMPLE = SampleNames,
      AGE = AgePlotMoy,
      HPD68.MIN = AgePlot68[, 2],
      HPD68.MAX = AgePlot68[, 3],
      HPD95.MIN = AgePlot95[, 2],
      HPD95.MAX = AgePlot95[, 3],
      stringsAsFactors = FALSE
    ),
    "Sampling" = echantillon,
    "PriorAge" = PriorAge,
    "StratiConstraints" = StratiConstraints,
    "Model_OSL_GrowthCurve" = Model_GrowthCurve,
    "Model_OSL_Distribution" = distribution,
    "CovarianceMatrix" = THETA,
    "Model_C14" = Model_C14,
    "CalibrationCurve" = CalibrationCurve,
    "Outlier" = Outlier
  )

  # Plot ages -----------------------------------------------------------------------------------
  plot_Ages(object = output, legend.pos = "bottomleft")

  ##TODO: get rid of this
  if(SavePdf){
    dev.print(
      pdf,
      file = paste(OutputFilePath, OutputFileName[3], '.pdf', sep = ""),
      width = 8,
      height = 10
    )
  }


  # Return output -------------------------------------------------------------------------------
  return(output)

}