Age_OSLC14.R

#' @title Bayesian analysis for age estimation of OSL measurements and C-14 ages of various samples
#'
#' @description This function computes an age of OSL data consisting of at least two samples and calibrate
#' C-14 ages of samples to get an age (in ka).\cr
#' Ages 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 simultaneously
#' (with output of [Generate_DataFile] or [Generate_DataFile_MG] or both of them using [combine_DataFiles]).
#' Samples, for which data is available in several BIN files, can be analysed.\cr
#' For C-14 data, the user can choose one of the following radiocarbon calibration curve:
#' Northern or Southern Hemisphere or marine atmospheric.
#'
#' @param DATA Two types of inputs are supported.
#' (1): a list of objects: `LT`, `sLT`, `ITimes`, `dLab`, `ddot_env`, `regDose`, `J`,
#' `K`, `Nb_measurement`, provided by the function [Generate_DataFile], [Generate_DataFile_MG] or [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
#' increasing ages. See the details section to for more informations.
#' (2): an object of class `BayLum.list` which is provided by the output of [Age_OSLC14]. When input of class `BayLum.list` is identified, no new JAGS model is created. Instead, the JAGS model specified within the `BayLum.list` is extended. Useful for when convergence was not originally achieved and a complete restart is not desirable.
#'
#' @param Data_C14Cal [numeric] vector: corresponding to C-14 age estimate
#' (in years, conversion in ka is automatically done in the function).
#' If there is stratigraphic relations between samples, `Data_C14Cal` must be ordered
#' by order of increasing ages.
#'
#' @param Data_SigmaC14Cal [numeric]: corresponding to the error of C-14 age estimates.
#'
#' @param Nb_sample [numeric]: number of samples (OSL data and C-14 age),
#' (`Nb_sample>3`, at least to sample of OSL data and one sample of C-14 age).
#'
#' @param SampleNames [character]: sample names for both OSL data and C14 data.
#' The length of this vector is equal to `Nb_sample`.
#' If there is stratigraphic relation, this vector must be ordered by increasing
#' order (to mix OSL samples and C-14 ages if it is needed).
#'
#' @param SampleNature [matrix]: states the nature of the sample.
#' Row number of `SampleNature` matrix is equal to `2` and column number is equal
#' to `Nb_sample`. First line of the matrix: `SampleNature[1,i]` states if sample
#' whose number ID is equal to `i`, is an OSL sample `1` or not `0`.
#' Second line of the matrix: `SampleNature[2,i]` states if sample whose number ID
#' is equal to `i`, is an C-14 sample `1` or not `0.
#'
#' @param PriorAge [numeric] (with default): lower and upper bounds for age parameter of each sample **(in ka)**.
#'  Note that, `length(PriorAge) = 2*Nb_sample` sand `PriorAge[2i-1,2i]` corresponds
#'  to the lower and upper bounds of sample whose number ID is equal to `i`.
#'
#' @param SavePdf [logical] (with default): if `TRUE` save graphs in pdf file
#' named `OutputFileName` in folder `OutputFilePath`.
#'
#' @param OutputFileName [character] (with default): name of the pdf file that
#' will be generated by the function if `SavePdf=TRUE`, `length(OutputFileName)=3`,
#' see **PLOT OUTPUT** in **Value** section for more informations.
#'
#' @param OutputFilePath [character] (with default): path to the pdf file that will
#' be generated by the function if `SavePdf=TRUE`. If it is not equal to "", it must be terminated by "/".
#'
#' @param SaveEstimates [logical] (with default): if TRUE save Bayes' estimates,
#' credible interval at level 68% and 95%  and the result of the Gelman en Rubin
#' test of convergence, in a CSV-table named `OutputFileName` in folder `OutputFilePath`.
#'
#' @param OutputTableName [character] (with default): name of the table that will
#' be generated by the function if `SaveEstimates=TRUE`.
#'
#' @param OutputTablePath [character] (with default): path to the table that will
#' be generated by the function if `SaveEstimates=TRUE`. If it is not equal to "",
#' it must be terminated by "/".
#'
#' @param StratiConstraints [matrix] or [character] (with default): input object for the
#' stratigraphic relation between samples. If there is stratigraphic relation between
#' samples see the details section for instructions regarding how to correctly fill `StratiConstraints`,
#' the user can refer to a [matrix] or to a CSV-file [character]. Otherwise, default value is suitable.
#'
#' @param sepSC [character] (with default): if `StratiConstraints` is character,
#' indicate column separator in `StratiConstraints` CSV-file.
#'
#' @param BinPerSample [numeric] (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]` corresponds to the number of BIN files for the sample whose
#' number ID is equal to `i`.  For more information to fill this vector, we refer
#' to details in [Generate_DataFile] or in [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 `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 `THETA` is character,
#' indicate column separator in `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 p[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 `"cauchy"`, `"gaussian"`, `"lognormal_A"` and `"lognormal_M"`,
#' see details for more informations.
#'
#' @param Model_C14 [character] (with default): if `"full"`, error on estimate calibration
#' curve is taken account, for C-14 samples.  If `"naive"` this error is not taken account in the age estimate.
#'
#' @param CalibrationCurve [character] (with default): calibration curve chosen, for C-14 samples.
#' Allowed inputs are
#' \itemize{
#'   \item \bold{"Intcal13"} or \bold{"Intcal13"} for Northern Hemisphere atmospheric radiocarbon calibration curve,
#'   \item \bold{"Marine13"} or \bold{"Marine13"} for Marine radiocarbon calibration curve,
#'   \item \bold{"SHCal13"} or \bold{"SHCal20"} 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 "XC-14.age", the third to "Error".
#'   The decimal of this file must be a dot, and the separator must be a comma. }
#' }
#'
#' @param adapt [integer] (with default): the number of iterations used in the adaptive phase of the simulation (see [runjags::run.jags]]).
#' @param burnin [integer] (with default): the number of iterations used to "home in" on the stationary posterior distribution. These are not used for assessing convergence (see [runjags::run.jags]]).
#' @param Iter [integer] (with default): the number of iterations to run and who will be used to assess convergence and ages (see [runjags::run.jags]]).
#'
#' @param t [numeric] (with default): 1 every `t` iterations of the MCMC is
#' considered for sampling the posterior distribution (for more information see [[rjags::jags.model].
#'
#' @param n.chains [numeric] (with default): number of independent chains for the model
#' (for more information see [[rjags::jags.model]).
#'
#' @param jags_method [character] (with default): select which method to use in order to call JAGS, supported are  `"rjags"` (the default), `rjparallel`, `simple`, `interruptible`, `parallel`, and `snow` (for more information about each of these possibilities, see [runjags::run.jags]])
#'
#' @param autorun [logical] (with default): choose to automate JAGS processing. JAGS model will be automatically extended until convergence is reached (for more information see [runjags::autorun.jags]).
#'
#' @param quiet [logical] (with default): enables/disables [rjags] messages
#'
#' @param roundingOfValue [integer] (with default):  Integer indicating the number of decimal places to be used, default = 3.
#
#' @param ... further arguments that can be passed to control the Bayesian process, see details
#' for supported arguments
#'
#' @details
#'
#' Note that there are three types of arguments in the previous list.
#' There are arguments for information concerning only OSL samples: \code{DATA}, \code{BinPerSample}, \code{THETA},
#' \code{sepTHETA}, \code{LIN_fit}, \code{Origin_fit}, \code{distribution}.
#'
#' There are arguments for information concerning only C14 samples: \code{Data_C14Cal}, \code{Data_SigmaC14Cal},
#' \code{Model_C14}, \code{CalibrationCurve}.
#'
#' There are arguments for information 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
#'
#' **Supported `...` arguments**
#' \tabular{lllll}{
#' ARGUMENT \tab INPUT \tab METHOD \tab DEFAULT \tab DESCRIPTION\cr
#' `max.time` \tab [character] \tab `rjparallel` \tab `Inf` \tab maximum allowed processing time, e.g.,
#' `10m` for 10 minutes (cf. [runjags::autorun.jags])\cr
#' `interactive` \tab [logical] \tab `rjparallel` \tab `FALSE` \tab enable/disable interactive mode (cf. [runjags::autorun.jags])\cr
#' `startburnin` \tab [integer] \tab `rjparallel` \tab  `4000` \tab number of burn-in iterations (cf. [runjags::autorun.jags]) \cr
#' `startsample` \tab [integer] \tab `rjparallel` \tab `10000` \tab total number of samples to assess convergence
#' (cf. [runjags::autorun.jags]) \cr
#' `inits` \tab named [list] \tab `rjparallel` \tab `NA` \tab fine control over random seeds and random number generator [runjags::autorun.jags]
#' }
#'
#' **How to fill `StratiConstraints?**\cr
#'
#' If there are 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 corresponds 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 construct
#' the \code{StratiConstraints} matrix.
#'
#' The user can also refer to a csv file that contains the relation between samples as defined above.
#' The user must be careful about which separator is used in the csv file using the argument \code{sepSC}.\cr
#'
#' **How to fill `THETA` 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]} contains 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]} contains common error between samples whose number ID are equal to \code{i} and \code{j}.
#' }
#' Note that \code{THETA[i,j]} is a symmetric matrix.
#'
#' The user can also refer to a .csv file that contains the errors as defined above.\cr
#'
#' **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}
#'   }
#' }
#'
#' **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
#' }
#'
#' **More precision on `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 developed by Christen, JA (1994), assume multiplicative effect of errors to address the
#' problem of outliers. In addition, to not penalise 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 convergence for the age estimate for each sample.
#' A result close to one is expected.\cr
#' In addition, the user must visually assess the convergence of the trajectories by looking at the graph
#' generated by the function (see \bold{PLOT OUTPUT} for more informations).\cr
#' If both convergences (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 iterations (\code{Iter})
#' or be more precise on the \code{PriorAge} parameter to reach convergence.
#'   \item \bold{Credible intervals and Bayes estimates}: prints the Bayes' estimates, the credible intervals at 95% and 68% for
#' the age parameters for each sample.
#'  \item \bold{JAGS model output}: returns the JAGS model with class "runjags".
#' }
#'
#' \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, Frederik Harly Baumgarten
#'
#' @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
#' [runjags], [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
#' \dontrun{
#' 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 = 20,
#'  burnin = 20,
#'  adapt = 20,
#'  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("IntCal20"),
    Iter = 10000,
    burnin = 4000,
    adapt = 1000,
    t = 5,
    n.chains = 3,
    jags_method = "rjags",
    autorun = FALSE,
    quiet = FALSE,
    roundingOfValue = 3,
    ...
) {
  if (inherits(DATA, "BayLum.list")) {
    ## reattach mcmc-list to runjags_object
    DATA$runjags_object$mcmc <- DATA$Sampling
      
    ind_OSL <- which(DATA$runjags_object$args$SampleNature[1,] == 1)
    CS_OSL <- cumsum(DATA$runjags_object$args$SampleNature[1,])
    ind_C14 <- which(DATA$runjags_object$args$SampleNature[2,] == 1)
    CS_C14 <- cumsum(DATA$runjags_object$args$SampleNature[2,])


    AgeBP = rev(DATA$runjags_object$args$TableauCalib[, 1])
    CalC14 = rev(DATA$runjags_object$args$TableauCalib[, 2])
    SigmaCalC14 = rev(DATA$runjags_object$args$TableauCalib[, 3])

    results_runjags <-
      runjags::extend.JAGS(
        runjags.object = DATA$runjags_object,
        adapt = adapt,
        burnin = burnin,
        sample = Iter,
        thin = t,
        method = jags_method,
        silent.jags = quiet,
        ...
      )

    # storing the arguments used for the orignal BayLum run (as to not lose them when results are processed)
    results_runjags$args <- list(
      "Model_OSL_GrowthCurve" = DATA$runjags_object$args$Model_OSL_GrowthCurve,
      "Model_OSL_Distribution" = DATA$runjags_object$args$Model_OSL_Distribution,
      "PriorAge" = DATA$runjags_object$args$PriorAge,
      "StratiConstraints" = DATA$runjags_object$args$StratiConstraints,
      "CovarianceMatrix" = DATA$runjags_object$args$CovarianceMatrix,
      "Model_C14" = DATA$runjags_object$args$Model_C14,
      "TableauCalib" = DATA$runjags_object$args$TableauCalib,
      "Outlier" = DATA$runjags_object$args$Outlier,
      "SampleNature" = DATA$runjags_object$args$SampleNature,
      "Data_C14Cal" = DATA$runjags_object$args$Data_C14Cal,
      "Nb_sample" = DATA$runjags_object$args$Nb_sample
    )
  }

  if(!inherits(DATA, "BayLum.list")) {
    #--- 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 %in%  c(
      "IntCal13", "IntCal20", "Marine13", "Marine20", "SHCal13" , "SHCal20")) {
      TableauCalib = get(data(list = CalibrationCurve, envir = environment()))
    } else {
      TableauCalib = read.csv(file = CalibrationCurve,
                              sep = ",",
                              dec = ".")
    }


    AgeBP = rev(TableauCalib[, 1])
    CalC14 = rev(TableauCalib[, 2])
    SigmaCalC14 = rev(TableauCalib[, 3])

    # #--- C14 preparation: 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 (q == 0) {
        stop(
          "[Age_OSLC14()] If you see this message, you are probably trying to run the model with a small number of samples.
           You can still use the function, but the C-14 sample cannot be the first sample.",
          call. = FALSE
        )

      }

      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
      )
    }

    # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
    # JAGS RUN --------------------- START
    # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#

    ##further settings provided eventually
    process_settings <- modifyList(x = list(
      max.time = Inf,
      interactive = FALSE,
      startburnin = 4000,
      startsample = 10000,
      inits = NA

    ), val = list(...))

    if (autorun == FALSE) {
      ##as input a text file is wanted, so we have to cheat a little bit
      temp_file <- tempfile(fileext = ".txt")
      writeLines(BUGModel, con = temp_file)

      ##run JAGS
      results_runjags <-
        runjags::run.JAGS(
          model = temp_file,
          data = dataList,
          n.chains = n.chains,
          monitor = c("A", "Z"),
          adapt = adapt,
          burnin = burnin,
          sample = Iter,
          silent.jags = quiet,
          method = jags_method,
          thin = t
        )
    }

    if (autorun == TRUE) {
      ##as input a text file is wanted, so we have to cheat a little bit
      temp_file <- tempfile(fileext = ".txt")
      writeLines(BUGModel, con = temp_file)

      ##run the auto processing
      results_runjags <-
        runjags::autorun.jags(
          model = temp_file,
          data = dataList,
          n.chains = n.chains,
          monitor = c("A", "Z"),
          adapt = adapt,
          silent.jags = quiet,
          method = jags_method,
          inits = process_settings$inits,
          max.time = process_settings$max.time,
          interactive = process_settings$interactive,
          startburnin = process_settings$startburnin,
          startsample = process_settings$startsample
        )
    }

    results_runjags$args <- list(
      "Model_OSL_GrowthCurve" = Model_GrowthCurve,
      "Model_OSL_Distribution" = distribution,
      "PriorAge" = PriorAge,
      "StratiConstraints" = StratiConstraints,
      "CovarianceMatrix" = THETA,
      "Model_C14" = Model_C14,
      "TableauCalib" = TableauCalib,
      "SampleNature" = SampleNature,
      "Data_C14Cal" = Data_C14Cal,
      "Nb_sample" = Nb_sample
    )
  }

##extract mcmc list
echantillon <- results_runjags$mcmc
U <- summary(echantillon)

##remove mcmc-list from runjags output to reduce output object size
results_runjags$mcmc <- list("MCMC-list is not here. Go to first level -> object named *Sampling*")

##combine chains into one data.frame
Sample <- as.data.frame(runjags::combine.mcmc(echantillon))

## attach "A_" to sample names for the A-parameter
nom = c()
for (i in 1:results_runjags$args$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[(results_runjags$args$Nb_sample + 1):(results_runjags$args$Nb_sample + sum(results_runjags$args$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:results_runjags$args$Nb_sample) {
  #cat(paste(" Sample name: ", SampleNames[i],"\n"))
  #cat("---------------------\n")
  cat(paste(
    paste("A_", SampleNames[i], sep = ""),
    "\t",
    round(CV$psrf[i, 1], roundingOfValue),
    "\t\t",
    round(CV$psrf[i, 2], roundingOfValue),
    "\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:results_runjags$args$Nb_sample) {
  rnames = c(rnames, paste("A_", SampleNames[i], sep = ""))
}
R = matrix(
  data = NA,
  ncol = 8,
  nrow = results_runjags$args$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 = results_runjags$args$Nb_sample, ncol = 3)
AgePlot68 = matrix(data = NA, nrow = results_runjags$args$Nb_sample, ncol = 3)
AgePlotMoy = rep(0, results_runjags$args$Nb_sample)
for (i in 1:results_runjags$args$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]), roundingOfValue),
    '\n'
  ))
  cat("\t\t\t\t\t\t lower bound \t upper bound\n")
  HPD_95 = ArchaeoPhases::CredibleInterval(Sample[, i], 0.95, roundingOfValue =
                                             roundingOfValue)
  HPD_68 = ArchaeoPhases::CredibleInterval(Sample[, i], 0.68, roundingOfValue =
                                             roundingOfValue)
  cat(
    "\t\t\t\t at level 95% \t",
    round(c(HPD_95[2]), roundingOfValue),
    "\t\t",
    round(c(HPD_95[3]), roundingOfValue),
    "\n"
  )
  cat(
    "\t\t\t\t at level 68% \t",
    round(c(HPD_68[2]), roundingOfValue),
    "\t\t",
    round(c(HPD_68[3]), roundingOfValue),
    "\n"
  )
  AgePlot95[i, ] = HPD_95
  AgePlot68[i, ] = HPD_68
  AgePlotMoy[i] = round(mean(Sample[, i]), roundingOfValue)

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

}

cat("\n------------------------------------------------------\n")
R[, c(7, 8)] <- round(CV$psrf[1:results_runjags$args$Nb_sample, ], roundingOfValue)

# Representation graphique des resultats
#        des HPD sur la courbe de calibration
if (sum(results_runjags$args$SampleNature[2, ]) > 1) {
  couleur = rainbow(results_runjags$args$Nb_sample)
  par(mfrow = c(1, 1),
      las = 0,
      mar = c(5, 5, 2, 2))
  xl = c(min(results_runjags$args$PriorAge[seq(1, (2 * results_runjags$args$Nb_sample - 1), 2)]), max(results_runjags$args$PriorAge[seq(2, (2 *
                                                                               results_runjags$args$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(results_runjags$args$Data_C14Cal[CS_C14[i]], 2),
          col = couleur[i],
          lwd = 4)
    lines(
      AgePlotMoy[i],
      results_runjags$args$Data_C14Cal[CS_C14[i]],
      col = "black",
      lwd = 2,
      type = 'p'
    )
  }
  legend(
    "topleft",
    SampleNames[ind_C14],
    lty = rep(1, results_runjags$args$Nb_sample),
    lwd = rep(2, results_runjags$args$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 == TRUE) {
  write.csv(R, file = c(
    paste(OutputTablePath, "Estimates", OutputTableName, ".csv", sep = "")
  ))
}

.list_BayLum <- function(..., originator = sys.call(which = -1)[[1]]){
  ## set list
  l <- list(...)

  ## update originators
  attr(l, "class") <- "BayLum.list"
  attr(l, "originator") <- as.character(originator)

  return(l)

}
# Create return object -------------------------------------------------------------------------
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" = results_runjags$args$PriorAge,
  "StratiConstraints" = results_runjags$args$StratiConstraints,
  "Model_OSL_GrowthCurve" = results_runjags$args$Model_OSL_GrowthCurve,
  "Model_OSL_Distribution" = results_runjags$args$Model_OSL_Distribution,
  "CovarianceMatrix" = results_runjags$args$CovarianceMatrix,
  "Model_C14" = results_runjags$args$Model_C14,
  "CalibrationCurve" = results_runjags$args$CalibrationCurve,
  "Outlier" = Outlier,
  "runjags_object" = results_runjags
)

# 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)

}