AgeS_Computation.R

#' @title Bayesian analysis for OSL age estimation of various samples
#'
#' @description This function computes the age (in ka) of at least two samples
#' according to the model developed in Combès and Philippe (2017),
#' based on outputs of [Generate_DataFile] or [Generate_DataFile_MG]
#' or both of them using [combine_DataFiles].\cr
#' Samples, for which data is available in several BIN files, can be analysed.\cr
#' Single-grain or Multi-grain OSL measurements can be analysed simultaneously.
#'
#' @param DATA (**required**) Two inputs are possible:
#' (1): DATA [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 informations 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 [AgeS_Computation]. When input of class `BayLum.list` is identified, no new JAGS model is created. Instead, the JAGS model specified by the [AgeS_Computation] output is extended. Useful for when convergence was not originally achieved and a complete restart is not desirable.
#'
#' @param SampleNames [character] vector: names of samples. The length of this vector is equal to `Nb_sample`.
#'
#' @param Nb_sample [integer]: number of samples, `Nb_sample>1`.
#'
#' @param PriorAge [numeric] vector (with default): lower and upper bounds for age parameter of each sample (in ka).
#' Note that, \code{length(PriorAge)=2*Nb_sample}
#' and `PriorAge[2i-1,2i]` corresponds to the lower and upper bounds of sample whose number ID is equal to `i`.
#'
#' @param BinPerSample [integer] vector (with default): vector with the number of BIN files per sample.
#' The length of this vector is equal to `Nb_sample`.
#' `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 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`;
#' \code{length(OutputFileName)=2}, 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 [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 \code{OutputFileName} in folder \code{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 THETA [numeric] [matrix] or [character] (with default): input object for systematic and individual error.
#' 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 errors are considered.
#'
#' @param sepTHETA [character] (with default): if `THETA` is character,
#' indicate column separator in `THETA` CSV-file.
#'
#' @param StratiConstraints  numeric 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 (numeric matrix) or to a csv file (character).
#' If there is no stratigraphic relation default value is suitable.
#'
#' @param sepSC [character] (with default): if \code{StratiConstraints} is character,
#' indicate column separator in \code{StratiConstraints} .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.
#' See details section 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 section for more informations on the proposed growth curves.
#'
#' @param distribution [character] (with default): type of distribution that defines
#' how individual equivalent dose values are distributed around the palaeodose.
#' Allowed inputs are `"cauchy"`, `"gaussian"`, `"lognormal_A"` and `"lognormal_M"`,
#' see details section for more informations.
#'
#' @param model [character] (*optional*): allows to provide a custom model to the function as text string. Please note that if this option is chosen the parameter `distribution` is ignored and no safety net is applied. If the function crashes it is up to the user.
#'
#' @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 which will be used to assess convergence and ages (see [runjags::run.jags]).
#'
#' @param t [integer] (with default): 1 every \code{t} iterations of the MCMC is considered for sampling the posterior distribution.
#' (for more information see [runjags::run.jags]).
#'
#' @param n.chains [integer] (with default): number of independent chains for the model (for more information see [runjags::run.jags]).
#'
#' @param jags_method [character] (with default): select which method to use in order to call JAGS. jags_methods `"rjags"` (the default) and `"rjparallel"` have been tested. (for more information about these possibilities and others, 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. 1) When creating a new JAGS model, see details
#' for supported arguments. 2) If extending a JAGS model see arguments of [runjags::extend.JAGS].
#'
#'
#' @details **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 is stratigraphic relations between samples,
#' *informations in DATA must be ordered by order of increasing ages*.
#' To do this the user can either fill right `Names` in [Generate_DataFile]
#' or in [Generate_DataFile_MG] (as it is indicated in Details section of these function),
#' or ordered by order of increasing ages outputs of [Generate_DataFile]
#' or [Generate_DataFile_MG] in [combine_DataFiles]
#'
#' 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+A,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) functions
#' 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 take care about the separator 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 follows.
#' \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 symetric matrix.
#'
#' The user can also refer to a .csv file that contains the errors as defined above.\cr
#'
#' Alternatively you can use the function [create_ThetaMatrix].
#'
#' **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.
#' }
#'
#' @return
#' **NUMERICAL OUTPUT**
#'
#' \enumerate{
#' \item **A list of type `BayLum.list` containing the following objects:**
#'  \itemize{
#'   \item \bold{Sampling}: that corresponds to a sample of the posterior distributions
#'  of the age (in ka), palaeodose (in Gy) and equivalent dose dispersion (in Gy) parameters for each sample;
#'   \item \bold{Model_GrowthCurve}: stating which dose response fitting option was chosen;
#'   \item \bold{Distribution}: stating which distribution was chosen to model the dispersion of
#'  individual equivalent dose values around the palaeodose of the sample;
#'   \item \bold{PriorAge}: stating the priors used for the age parameter (in ka);
#'   \item \bold{StratiConstraints}: stating the stratigraphic relations between samples considered in the model;
#'   \item \bold{CovarianceMatrix}: stating the covariance matrix of error used in the model, highlighting common errors between samples or not.
#'   \item \bold{model}: returns the model that was used for the Bayesian modelling as a [character]
#'   \item \bold{JAGS model output}: returns the JAGS model with class "runjags".
#'  }
#'   \item\bold{The Gelman and Rubin test of convergency}: prints the result of the Gelman and Rubin test of convergence for
#' the age, palaeodose and equivalent dose dispersion parameters 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 being more precise on the \code{PriorAge} parameter (for example specify if it is a young sample \code{c(0.01,10)} an old sample \code{c(10,100)}),
#' or changing the parameter \code{distribution} or the growth curve, to reach convergence.
#'   \item \bold{Credible intervals and Bayes estimates}: prints the Bayes estimates, the credible intervals at 95% and 68% for
#' the age, palaeodose and equivalent dose dispersion parameters for each sample.
#' }
#'
#' **PLOT OUTPUT**
#'
#' \enumerate{
#'  \item\bold{MCMC trajectories}: A graph with the MCMC trajectories and posterior distributions of the age, palaeodose and equivalent dose dispersion parameters
#' is displayed, there is one page per sample. \cr
#' The first line of the figure corresponds to the age parameter, the second to the palaeodose parameter
#' and the third to the equivalent dose dispersion parameter.
#' 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{Summary of sample age estimates}: plot credible intervals and Bayes estimate of each sample age on a same graph.
#' }
#'
#' To give the results in a publication, we recommend to give the Bayes' estimate of the parameters as well as the credible interval at 95% or 68%.
#'
#' @author Claire Christophe, Anne Philippe, Guillaume Guérin, 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 [Generate_DataFile], [Generate_DataFile_MG], [rjags], [plot_MCMC], [SCMatrix], [Age_Computation], [Palaeodose_Computation], [plot_Ages], [create_ThetaMatrix], [runjags]
#'
#' @references
#' Combes, Benoit and Philippe, Anne, 2017.
#' Bayesian analysis of multiplicative Gaussian error for multiple ages estimation in optically stimulated luminescence dating.
#' Quaternary Geochronology (39, 24-34)
#'
#' Combes, B., Philippe, A., Lanos, P., Mercier, N., Tribolo, C., Guerin, G., Guibert, P., Lahaye, C., 2015.
#' A Bayesian central equivalent dose model for optically stimulated luminescence dating.
#' Quaternary Geochronology 28, 62-70. doi:10.1016/j.quageo.2015.04.001
#'
#' @examples
#' ## Age computation of samples GDB5 and GDB3,
#' ## load data
#' data(DATA2) # GD85
#' data(DATA1) # GD83
#'
#' ## produce DataFile
#' Data <- combine_DataFiles(DATA2, DATA1)
#'
#' ## without common error, assuming GDB5 age younger than GDB3 age
#' Nb_sample <- 2
#' SC <- matrix(
#'   data = c(1,1,0,1,0,0),
#'   ncol = 2,
#'   nrow = (Nb_sample+1),
#'   byrow = TRUE)
#'
#'\dontrun{
#' ## run standard
#' Age <- AgeS_Computation(
#'   DATA = Data,
#'   Nb_sample = Nb_sample,
#'   SampleNames = c("GDB5","GDB3"),
#'   PriorAge = rep(c(1,100), 2),
#'   StratiConstraints = SC,
#'   Iter = 100,
#'   quiet = FALSE,
#'   jags_method = "rjags"
#' )
#'
#' ## extend model
#' Age_extended <- AgeS_Computation(
#'   DATA = Age,
#'   Nb_sample = Nb_sample,
#'   SampleNames = c("GDB5","GDB3"),
#'   PriorAge = rep(c(1,100), 2),
#'   StratiConstraints = SC,
#'   adapt = 0,
#'   burnin = 500,
#'   Iter = 1000,
#'   quiet = FALSE,
#'   jags_method = "rjags"
#' )
#'
#' ## autorun mode
#' Age <- AgeS_Computation(
#'   DATA = Data,
#'   Nb_sample = Nb_sample,
#'   SampleNames = c("GDB5","GDB3"),
#'   PriorAge = rep(c(1,100), 2),
#'   StratiConstraints = SC,
#'   Iter = 10000,
#'   quiet = FALSE,
#'   jags_method = "rjags",
#'   autorun = TRUE)
#'
#' ## parallel mode
#' Age <- AgeS_Computation(
#'   DATA = Data,
#'   Nb_sample = Nb_sample,
#'   SampleNames = c("GDB5","GDB3"),
#'   PriorAge = rep(c(1,100), 2),
#'   StratiConstraints = SC,
#'   Iter = 10000,
#'   quiet = FALSE,
#'   jags_method = "rjparallel")
#' }
#'
#' @md
#' @export
AgeS_Computation <- function(
    DATA,
    SampleNames,
    Nb_sample,
    PriorAge = rep(c(0.01, 100), Nb_sample),
    BinPerSample = rep(1, Nb_sample),
    SavePdf = FALSE,
    OutputFileName = c('MCMCplot', "summary"),
    OutputFilePath = c(""),
    SaveEstimates = FALSE,
    OutputTableName = c("DATA"),
    OutputTablePath = c(''),
    THETA = c(),
    sepTHETA = c(','),
    StratiConstraints = c(),
    sepSC = c(','),
    LIN_fit = TRUE,
    Origin_fit = FALSE,
    distribution = c("cauchy"),
    model = NULL,
    Iter = 10000,
    burnin = 4000,
    adapt = 1000,
    t = 5,
    n.chains = 3,
    jags_method = "rjags",
    autorun = FALSE,
    quiet = FALSE,
    roundingOfValue = 3,
    ...
) {
  #---check to see if DATA input is a runjags-object and extend if so ####
  if (inherits(DATA, "BayLum.list")) {
      # reattach mcmc-list which was removed to reduce object size
      DATA$runjags_object$mcmc <- DATA$Sampling
      
     # extend via runjags
    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 original BayLum run (as to not lose them when results are processed)
    results_runjags$args <- list(
      "Model_GrowthCurve" = DATA$runjags_object$args$Model_GrowthCurve,
      "Distribution" = DATA$runjags_object$args$Distribution,
      "PriorAge" = DATA$runjags_object$args$PriorAge,
      "StratiConstraints" = DATA$runjags_object$args$StratiConstraints,
      "CovarianceMatrix" = DATA$runjags_object$args$CovarianceMatrix,
      "model" = DATA$runjags_object$args$model
    )
  }
  #---check to see if DATA input is a DataFile and run JAGS ####
  if (!inherits(DATA, "BayLum.list")) {
    ##---Index preparation ####
    CSBinPerSample <- cumsum(BinPerSample)
    LengthSample <- c()
    for (ns in 1:Nb_sample) {
      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:Nb_sample) {
      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 ####
    ## cover character case
    if (is(THETA)[1] == "character")
      THETA <- as.matrix(read.csv(THETA, sep = sepTHETA, header = FALSE))

    ## cover data.frame case (users should not do it but they do)
    if (inherits(THETA, "data.frame"))
      THETA <- as.matrix(THETA)

    if (length(THETA[,1]) == 0) {
      THETA <- diag(DATA$ddot_env[2, CSBinPerSample] +
                      (DATA$ddot_env[1, CSBinPerSample]) ^ 2 * DATA$dLab[2, CSBinPerSample])
    }

    ##JAGS will crash with a runtime error if the dimension of the theta matrix does not match the
    ##number of samples
    if (sum(dim(THETA)) %% Nb_sample != 0)
      stop(
        "[AgeS_Computation()] The number of samples does not match the dimension of the THETA-matrix!",
        call. = 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)
      }
    }

    ##---BUG file selection ####
    Model_AgeS <- 0
    data(Model_AgeS, 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("AgesMultiCS2_EXP", cLIN, cO, sep = ""))

    # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
    ##---JAGS RUN --------------------- START ####
    # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
    ##set variables
    dataList <- list(
      'N' = LT,
      'sN' = sLT,
      "IT" = IrrT,
      "sDlab" = DATA$dLab[1,],
      'J' = DATA$J,
      'K' = DATA$K,
      "I" = Nb_sample,
      "ddot" = DATA$ddot_env[1, CSBinPerSample],
      "Gamma" = THETA,
      "xbound" = PriorAge,
      "StratiConstraints" = StratiConstraints,
      "index" = index2,
      "BinPerSample" = BinPerSample,
      "CSBinPerSample" = CSBinPerSample
    )

    ##select model
    if (is.null(model))
      model <- Model_AgeS[[Model_GrowthCurve]][[distribution]]

    ##there are two main ways of running JAGS: single-run vs auto-run
    ##(1) if user selects to do single-run:
    if (!autorun) {
      ##a text file is wanted as input, so we have to cheat a little bit
      temp_file <- tempfile(fileext = ".txt")
      writeLines(model, con = temp_file)

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

    ##(2) if user selects to do auto-run:
    if (autorun) {
      ##further settings provided eventually
      process_settings <- modifyList(x = list(
        max.time = Inf,
        interactive = FALSE,
        startburnin = 4000,
        startsample = 10000,
        inits = NA

      ), val = list(...))

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

      ##run the autoprocessing
      results_runjags <-
        runjags::autorun.JAGS(
          model = temp_file,
          data = dataList,
          n.chains = n.chains,
          monitor = c("A", "D", "sD"),
          adapt = adapt,
          startburnin = process_settings$startburnin,
          startsample = process_settings$startsample,
          silent.jags = quiet,
          method = jags_method,
          thin = t,
          inits = process_settings$inits,
          max.time = process_settings$max.time,
          interactive = process_settings$interactive
        )
    }
    # storing the arguments used for the BayLum-run this way,
    # because it allows us an easy way to code the storage of arguments when extending a JAGS-model.
    results_runjags$args <- list(
      "Model_GrowthCurve" = Model_GrowthCurve,
      "Distribution" = distribution,
      "PriorAge" = PriorAge,
      "StratiConstraints" = StratiConstraints,
      "CovarianceMatrix" = THETA,
      "model" = model
    )
    # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
    # JAGS RUN --------------------- END
    # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
  }

  #---processing of JAGS results
  ##extract mcmc list from runjags object
  echantillon <- results_runjags$mcmc
  
  ##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))

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

  ##try makes sure that the function runs
  try(plot_MCMC(echantillon, sample_names = SampleNames))

  if (SavePdf) {
    dev.off()
  }

  #---Gelman and Rubin test of convergence of the MCMC ####
  CV <- gelman.diag(echantillon, multivariate = FALSE)
  cat(paste(
    "\n\n>> Results of the Gelman and Rubin criterion of convergence <<\n"
  ))
  for (i in 1:Nb_sample) {
    cat("----------------------------------------------\n")
    cat(paste(" Sample name: ", SampleNames[i], "\n"))
    cat("---------------------\n")
    cat(paste("\t\t", "Point estimate", "Uppers confidence interval\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(paste(
      paste("D_", SampleNames[i], sep = ""),
      "\t",
      round(CV$psrf[(Nb_sample + i), 1], roundingOfValue),
      "\t\t",
      round(CV$psrf[(Nb_sample + i), 2], roundingOfValue),
      "\n"
    ))
    cat(paste(
      paste("sD_", SampleNames[i], sep = ""),
      "\t",
      round(CV$psrf[(2 * Nb_sample + i), 1], roundingOfValue),
      "\t\t",
      round(CV$psrf[(2 * Nb_sample + i), 2], roundingOfValue),
      "\n"
    ))
  }

  cat(
    "\n\n---------------------------------------------------------------------------------------------------\n"
  )
  cat(
    " *** WARNING: The following information are only valid if the MCMC chains have converged  ***\n"
  )
  cat(
    "---------------------------------------------------------------------------------------------------\n\n"
  )

  #---print results ####
  ##Matrix of results
  rnames <- rep(NA, 3 * Nb_sample)
  for (i in 1:Nb_sample) {
    rnames[i] = c(paste("A_", SampleNames[i], sep = ""))
    rnames[Nb_sample + i] = c(paste("D_", SampleNames[i], sep = ""))
    rnames[2 * Nb_sample + i] = c(paste("sD_", SampleNames[i], sep = ""))
  }
  R = matrix(
    data = NA,
    ncol = 8,
    nrow = 3 * 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: Point estimate",
        "Convergencies: uppers confidence interval"
      )
    )
  )

  ##Bayes estimate and credible interval
  cat(
    paste(
      "\n\n>> Bayes estimates of Age, Palaeodose and its dispersion 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(
      "Parameter",
      "\t",
      "Bayes estimate",
      "\t",
      " 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)

    cat(paste(
      "\nParameter",
      "\t",
      "Bayes estimate",
      "\t",
      " Credible interval \n"
    ))
    cat(paste(
      paste(" D_", SampleNames[i], sep = ""),
      "\t",
      round(mean(sample[, (Nb_sample + i)]), roundingOfValue),
      '\n'
    ))
    cat("\t\t\t\t\t\t lower bound \t upper bound\n")
    HPD_95 = ArchaeoPhases::CredibleInterval(sample[, (Nb_sample + i)], 0.95, roundingOfValue =
                                               roundingOfValue)
    HPD_68 = ArchaeoPhases::CredibleInterval(sample[, (Nb_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"
    )

    R[(Nb_sample + i), 3] = round(mean(sample[, (Nb_sample + i)]), roundingOfValue)
    R[(Nb_sample + i), c(1, 5)] = round(HPD_95[2:3], roundingOfValue)
    R[(Nb_sample + i), c(2, 4)] = round(HPD_68[2:3], roundingOfValue)

    cat(paste(
      "\nParameter",
      "\t",
      "Bayes estimate",
      "\t",
      " Credible interval \n"
    ))
    cat(paste(
      paste("sD_", SampleNames[i], sep = ""),
      "\t",
      round(mean(sample[, (2 * Nb_sample + i)]), roundingOfValue),
      '\n'
    ))
    cat("\t\t\t\t\t\t lower bound \t upper bound\n")
    HPD_95 = ArchaeoPhases::CredibleInterval(echantillon[[1]][, (2 * Nb_sample +
                                                                   i)], 0.95, roundingOfValue = roundingOfValue)
    HPD_68 = ArchaeoPhases::CredibleInterval(echantillon[[1]][, (2 * Nb_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"
    )

    R[(2 * Nb_sample + i), 3] = round(mean(sample[, (2 * Nb_sample + i)]), roundingOfValue)
    R[(2 * Nb_sample + i), c(1, 5)] = round(HPD_95[2:3], roundingOfValue)
    R[(2 * Nb_sample + i), c(2, 4)] = round(HPD_68[2:3], roundingOfValue)

    # R[(3*(i-1)+1):(3*i),6]=c('','','')
    # R[(3*(i-1)+1):(3*i),7]=round(CV$psrf[c(i,i+Nb_sample,i+2*Nb_sample),1],2)
    # R[(3*(i-1)+1):(3*i),8]=round(CV$psrf[c(i,i+Nb_sample,i+2*Nb_sample),2],2)


  }
  cat("\n----------------------------------------------\n")
  R[, c(7, 8)] <- round(CV$psrf, roundingOfValue)

  ##check if there is stratigraphic order between sample.
  ##if(sum(StratiConstraints[2:Nb_sample,1:Nb_sample]>0)){
  ##}

  #---print csv table ####
  if (SaveEstimates == TRUE) {
    write.csv(R, file = c(
      paste(
        OutputTablePath,
        "Estimates",
        OutputTableName,
        ".csv",
        sep = ""
      )
    ))
  }

  #---Create return object ####
  .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)

  }

  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,
    "Model_GrowthCurve" = results_runjags$args$Model_GrowthCurve,
    "Distribution" = results_runjags$args$Distribution,
    "PriorAge" = results_runjags$args$PriorAge,
    "StratiConstraints" = results_runjags$args$StratiConstraints,
    "CovarianceMatrix" = results_runjags$args$CovarianceMatrix,
    "model" = results_runjags$model,
    "runjags_object" = results_runjags
  )

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

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

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