R/read_XSYG2R.R

In Luminescence: Comprehensive Luminescence Dating Data Analysis

#' @title Import XSYG files into R
#'
#' @description Imports XSYG-files produced by a Freiberg Instruments lexsyg reader into R.
#'
#' @details
#' **How does the import function work?**
#'
#' The function uses the `'XML'` package to parse the file structure. Each
#' sequence is subsequently translated into an [Luminescence::RLum.Analysis-class] object.
#'
#' **General structure XSYG format**
#'
#' ```
#' <?xml?>
#' <Sample>
#'   <Sequence>
#'     <Record>
#'       <Curve name="first curve" />
#'       <Curve name="curve with data">x0 , y0 ; x1 , y1 ; x2 , y2 ; x3 , y3</Curve>
#'     </Record>
#'   </Sequence>
#' </Sample>
#' ```
#'
#' So far, each
#' XSYG file can only contain one `<Sample></Sample>`, but multiple
#' sequences.
#'
#' Each record may comprise several curves.
#'
#' **TL curve recalculation**
#'
#' On the FI lexsyg device TL curves are recorded as time against count values.
#' Temperature values are monitored on the heating plate and stored in a
#' separate curve (time vs. temperature). If the option
#' `recalculate.TL.curves = TRUE` is chosen, the time values for each TL
#' curve are replaced by temperature values.
#'
#' Practically, this means combining two matrices (Time vs. Counts and Time vs.
#' Temperature) with different row numbers by their time values. Three cases
#' are considered:
#'
#' 1. HE: Heating element
#' 2. PMT: Photomultiplier tube
#' 3. Interpolation is done using the function [approx]
#'
#' CASE (1): `nrow(matrix(PMT))` > `nrow(matrix(HE))`
#'
#' Missing temperature values from the heating element are calculated using
#' time values from the PMT measurement.
#'
#' CASE (2): `nrow(matrix(PMT))` < `nrow(matrix(HE))`
#'
#' Missing count values from the PMT are calculated using time values from the
#' heating element measurement.
#'
#' CASE (3): `nrow(matrix(PMT))` == `nrow(matrix(HE))`
#'
#' A new matrix is produced using temperature values from the heating element
#' and count values from the PMT.
#'
#' **Note:**
#' Please note that due to the recalculation of the temperature
#' values based on values delivered by the heating element, it may happen that
#' multiple count values exists for each temperature value and temperature
#' values may also decrease during heating, not only increase.
#'
#' **Count linearity correction**
#'
#' If the argument `auto_linearity_correction = TRUE`, the function offers a,
#' theoretically, automated linearity correction of the PMT signal using
#' the internal function [Luminescence::correct_PMTLinearity]. The critical parameter
#' is the count-pair resolution, which is provided to the function
#' automatically depending on the detector. Currently, the following
#' settings are used:
#'
#' \tabular{ll}{
#' DETECTOR \tab COUNT-PAIR-RESOLUTION\cr
#' UVVIS \tab 18 ns\cr
#' NIR50 \tab 70 ns\cr
#' NIR40 \tab 70 ns \cr
#' ETPMT \tab 25 ns
#' }
#'
#' Unfortunately, not all XSYG files provide correct information on
#' the used detector, depending on software version. Hence, in case
#' of doubt, please verify the settings and conduct a manual correction
#' if required.
#'
#' **Advanced file import**
#'
#' To allow for a more efficient usage of the function, instead of single path
#' to a file just a directory can be passed as input. In this particular case
#' the function tries to extract all XSYG-files found in the directory and import
#' them all. Using this option internally the function constructs as list of
#' the XSYG-files found in the directory. Please note no recursive detection
#' is supported as this may lead to endless loops.
#'
#' @param file [character] or [list] (**required**):
#' name of one or multiple XSYG files (URLs are supported); it can be the path
#' to a directory, in which case the function tries to detect and import all
#' XSYG files found recursively from the given directory.
#'
#' @param recalculate.TL.curves [logical] (*with default*):
#' if set to `TRUE`, TL curves are returned as temperature against count values
#' (see details for more information) Note: The option overwrites the time vs.
#' count TL curve. Select `FALSE` to import the raw data delivered by the
#' lexsyg. Works for TL curves and spectra.
#'
#' @param n_records [numeric] (*with default*):
#' number of records to be imported; by default the function attempts to import
#' all records.
#'
#' @param fastForward [logical] (*with default*):
#' if `TRUE` for a more efficient data processing only a list of
#' [Luminescence::RLum.Analysis-class] objects is returned.
#'
#' @param import [logical] (*with default*):
#' if set to `FALSE`, only the XSYG file structure is shown.
#'
#' @param pattern [character] (*with default*):
#' regular expression pattern passed to [list.files] to construct a list of
#' files to read (used only when a path is provided).
#'
#' @param auto_linearity_correction [logical] (*with default*): enable/disable
#' automatic count linearity correction. Because the information on the detectors
#' are not consistent and are not consistently stored, this option should be used
#' with caution.
#'
#' @param verbose [logical] (*with default*):
#' enable/disable output to the terminal.
#'
#' @param txtProgressBar [logical] (*with default*):
#' enable/disable the progress bar during import. Ignored if `verbose = FALSE`.
#'
#' @return
#' **Using the option `import = FALSE`**
#'
#' A list consisting of two elements is shown:
#' - [data.frame] with information on file.
#' - [data.frame] with information on the sequences stored in the XSYG file.
#'
#' **Using the option `import = TRUE` (default)**
#'
#' A list is provided, the list elements
#' contain: \item{Sequence.Header}{[data.frame] with information on the
#' sequence.} \item{Sequence.Object}{[Luminescence::RLum.Analysis-class]
#' containing the curves.}
#'
#' @note
#' This function is a beta version as the XSYG file format is not yet
#' fully specified. Thus, further file operations (merge, export, write) should
#' be done using the functions provided with the package `'XML'`.
#'
#' **So far, no image data import is provided!** \cr
#' Corresponding values in the XSXG file are skipped.
#'
#' @section Function version: 0.8.3
#'
#' @author
#' Sebastian Kreutzer, F2.1 Geophysical Parametrisation/Regionalisation, LIAG - Institute for Applied Geophysics (Germany)\cr
#' Marco Colombo, Institute of Geography, Heidelberg University (Germany)
#'
#' @seealso `'XML'`, [Luminescence::RLum.Analysis-class], [Luminescence::RLum.Data.Curve-class],
#' [approx], [Luminescence::correct_PMTLinearity]
#'
#' @references
#' Grehl, S., Kreutzer, S., Hoehne, M., 2013. Documentation of the
#' XSYG file format. Unpublished Technical Note. Freiberg, Germany
#'
#' **Further reading**
#'
#' XML: [https://en.wikipedia.org/wiki/XML]()
#'
#' @keywords IO
#'
#' @examples
#'
#' ##(1) import XSYG file to R (uncomment for usage)
#'
#' #FILE <- file.choose()
#' #temp <- read_XSYG2R(FILE)
#'
#' ##(2) additional examples for pure XML import using the package XML
#' ##    (uncomment for usage)
#'
#'   ##import entire XML file
#'   #FILE <- file.choose()
#'   #temp <- XML::xmlRoot(XML::xmlTreeParse(FILE))
#'
#'   ##search for specific subnodes with curves containing 'OSL'
#'   #getNodeSet(temp, "//Sample/Sequence/Record[@@recordType = 'OSL']/Curve")
#'
#' ##(2) How to extract single curves ... after import
#' data(ExampleData.XSYG, envir = environment())
#'
#' ##grep one OSL curves and plot the first curve
#' OSLcurve <- get_RLum(OSL.SARMeasurement$Sequence.Object, recordType="OSL")[[1]]
#'
#' ##(3) How to see the structure of an object?
#' structure_RLum(OSL.SARMeasurement$Sequence.Object)
#'
#' @export
read_XSYG2R <- function(
  file,
  recalculate.TL.curves = TRUE,
  n_records = NULL,
  fastForward = FALSE,
  import = TRUE,
  pattern = "\\.xsyg",
  auto_linearity_correction = FALSE,
  verbose = TRUE,
  txtProgressBar = TRUE
) {
  .set_function_name("read_XSYG2R")
  on.exit(.unset_function_name(), add = TRUE)

  ##TODO: this function should be reshaped:
  ##  - metadata from the sequence should go into the info slot of the RLum.Analysis object
  ##   >> however, the question is whether this works with subsequent functions
  ##  - currently not all metadata are supported, it should be extended
  ##  - the should be a mode importing ALL metadata
  ##  - xlum should be general, xsyg should take care about subsequent details

  .validate_logical_scalar(verbose)
  .validate_class(pattern, "character")
  file <- .validate_file(file, pattern = pattern, recursive = TRUE,
                         throw.error = FALSE, verbose = verbose)
  if (length(file) == 0)
    return(NULL)
  .validate_positive_scalar(n_records, int = TRUE, null.ok = TRUE)
  .validate_logical_scalar(fastForward)
  .validate_logical_scalar(import)

  # Self Call -----------------------------------------------------------------------------------
  # Option (a): Input is a list, every element in the list will be treated as file connection
  # with that many file can be read in at the same time
  # Option (b): The input is just a path, the function tries to grep ALL xsyg/XSYG files in the
  # directory and import them, if this is detected, we proceed as list
  if (inherits(file, "list")) {
    temp.return <- lapply(seq_along(file), function(x) {
      read_XSYG2R(
        file = file[[x]],
        recalculate.TL.curves = recalculate.TL.curves,
        fastForward = fastForward,
        import = import,
        verbose = verbose,
        txtProgressBar = txtProgressBar
      )
    })

    if (!fastForward)
      return(temp.return)
    if (import)
      return(unlist(temp.return, recursive = FALSE))

    return(as.data.frame(data.table::rbindlist(temp.return)))
  }

  ## Integrity checks -------------------------------------------------------
  .validate_logical_scalar(auto_linearity_correction)
  .validate_logical_scalar(txtProgressBar)

  ## don't show the progress bar if not verbose
  if (!verbose)
    txtProgressBar <- FALSE

  # (0) config --------------------------------------------------------------
  #version.supported <- c("1.0")

  #additional functions
  #get spectrum values
  # TODO: This function could be written also in C++, however, not necessary due to a low demand
  get_XSYG.spectrum.values <- function(curve.node){
    ##1st grep wavelength table
    wavelength <- XML::xmlAttrs(curve.node)["wavelengthTable"]

    ##string split (scan is considerably faster by a factor of 2)
    wavelength <- na.exclude(scan(text = wavelength, what = numeric(), sep = ";", quiet = TRUE))

    ##2nd grep time values
    raw_curve_node <- XML::xmlValue(curve.node)
    curve.node <- unlist(strsplit(raw_curve_node, split = "];", fixed = TRUE))
    curve.node <- strsplit(curve.node, split = ",[", fixed = TRUE)
    curve.node.time <- as.numeric(vapply(curve.node, function(x) x[1], character(1)))

    ##3rd grep count values
    curve.node.count <- vapply(curve.node, function(x) {
      x[length(x)]
    }, character(1))

    ## remove last bracket
    curve.node.count[length(curve.node.count)] <- gsub(
      pattern = "]", replacement = "", x = curve.node.count[length(curve.node.count)], fixed = TRUE)

    ##4th combine to spectrum matrix
    spectrum.matrix <- matrix(NA, nrow = length(wavelength), ncol = length(curve.node.time))
    for(i in 1:length(curve.node.time)) {
      tmp <- scan(text = curve.node.count[i], what = numeric(), sep = "|", quiet = TRUE)
      if(length(tmp) == length(wavelength))
        spectrum.matrix[,i] <- tmp
    }

    ## remove NA values from matrix
    id_NA <- colSums(is.na(spectrum.matrix)) != nrow(spectrum.matrix)
    spectrum.matrix <- spectrum.matrix[, id_NA]

    ##change row names (rows are wavelength)
    rownames(spectrum.matrix) <- round(wavelength, digits = 3)

    ##change column names (columns are time/temp values)
    colnames(spectrum.matrix) <- round(curve.node.time[id_NA], digits=3)

    return(spectrum.matrix)
  }

  # (1) Integrity checks ----------------------------------------------------
  ##set HUGE for larger nodes
  HUGE <- 524288

  ##parse XML tree using the package XML
  temp <- try(
    XML::xmlRoot(XML::xmlTreeParse(file, useInternalNodes = TRUE, options = HUGE, error = NULL)),
    silent = TRUE)

  ##show error
  if(inherits(temp, "try-error")){
    if(verbose)
      .throw_message("XML file not readable, nothing imported, NULL returned")
    return(NULL)
  }

  # (2) Further file processing ---------------------------------------------
  ##==========================================================================##
  ##SHOW STRUCTURE
  if(!import){
    ##sample information
    temp.sample <- as.data.frame(XML::xmlAttrs(temp), stringsAsFactors = FALSE)

    ##grep sequences files

    ##set data.frame
    names.sequence.header <- names(XML::xmlAttrs(temp[[1]]))
    temp.sequence.header <- data.frame(t(1:length(names.sequence.header)),
                                       stringsAsFactors = FALSE)
    colnames(temp.sequence.header) <- names.sequence.header

    ##fill information in data.frame
    for(i in 1:XML::xmlSize(temp)){
      temp.sequence.header[i,] <- t(XML::xmlAttrs(temp[[i]]))
    }

      ##additional option for fastForward == TRUE
      if(fastForward){
        ##change column header
        temp.sample <- t(temp.sample)
        colnames(temp.sample) <- paste0("sample::", colnames(temp.sample))
        output <- cbind(temp.sequence.header, temp.sample)

      }else{
        output <-  list(Sample = temp.sample, Sequences = temp.sequence.header)
      }

    return(output)
  }

  ## ========================================================================
  ## IMPORT XSYG FILE

  if (verbose) {
      cat("\n[read_XSYG2R()] Importing ...")
      cat("\n path: ", dirname(file))
      cat("\n file: ", .shorten_filename(basename(file)))
      cat("\n")
    }

  ## set n_records
  n_records <- min(XML::xmlSize(temp), n_records)

  ## initialize the progress bar
  if (txtProgressBar) {
    pb <- txtProgressBar(min = 0, max = n_records, char = "=", style = 3)
  }

  ## loop over the entire sequence up to the number of records to be read
  output <- lapply(1:n_records, function(x) {
    sequence <- temp[[x]]
    sequence.header <- as.data.frame(XML::xmlAttrs(sequence),
                                     stringsAsFactors = FALSE)

    ## account for non set value
    if (length(sequence.header) > 0)
        colnames(sequence.header) <- ""

    ###----------------------------------------------------------------------
    ##LOOP
    ##read records >> records are combined to one RLum.Analysis object
    sequence.object <- unlist(lapply(seq_len(XML::xmlSize(sequence)), function(i) {
      ## sequence becomes the record
      record <- sequence[[i]]

      ## the XSYG file might be broken due to a machine error during the measurement
      recordType <- try(XML::xmlAttrs(record)["recordType"], silent = TRUE)
      
      if (any(inherits(recordType, "try-error")))
          return(NULL) # nocov

      ## create a fallback, the function should not fail
      if (is.null(recordType) || is.na(recordType))
        recordType <- "not_set"

      ## correct record type in depending on the stimulator
      xml.size <- XML::xmlSize(record)
      if (recordType == "OSL" &&
          XML::xmlAttrs(record[[xml.size]])["stimulator"] %in%
              c("ir_LED_850", "ir_LD_850")) {
          recordType  <- "IRSL"
      }

      ## get all record attributes
      attrs_record <- XML::xmlAttrs(record)
      names(attrs_record) <- gsub("^name$", "recordName", names(attrs_record))

      header.position <- as.integer(as.character(sequence.header["position", ]))
      header.name <- as.character(sequence.header["name", ])

      ## loop 3rd level
      lapply(1:xml.size, function(j) {
        curve <- record[[j]]
        attrs <- XML::xmlAttrs(curve)

        ## all curves after the first in a record are marked with a leading
        ## underscore: this should make it easier to identify the curve to
        ## analyze (the first) from the others
        recordType <- paste0(ifelse(j == 1, "", "_"), recordType)

          ##get curveType
          temp.sequence.object.curveType <- as.character(attrs["curveType"])

          ##get detector
          temp.sequence.object.detector <- as.character(attrs["detector"])

          ## combine attributes
          attrs_comb <- as.list(c(attrs, attrs_record))
          attrs_comb <- attrs_comb[!duplicated(names(attrs_comb))]

        ## get additional information
        temp.sequence.object.info <- modifyList(attrs_comb,
                                                list(position = header.position,
                                                     sequenceName = header.name))

        ## TL curve recalculation ===========================================
        if (recalculate.TL.curves) {
            ##TL curve heating values is stored in the 3rd curve of every set
            if (recordType == "TL" && j == 1) {
              #grep values from PMT measurement or spectrometer
              if(!"Spectrometer" %in% temp.sequence.object.detector){
                temp.sequence.object.curveValue.PMT <- src_get_XSYG_curve_values(XML::xmlValue(
                  record[[j]]))

                time.values <- temp.sequence.object.curveValue.PMT[, 1]
              }else{
                temp.sequence.object.curveValue.spectrum <-
                  get_XSYG.spectrum.values(curve)

                ##get time values which are stored in the row labels
                time.values <- as.numeric(
                  colnames(temp.sequence.object.curveValue.spectrum))
              }

              ## round values (1 digit is technical resolution of the heating element)
              time.values <- round(time.values, digits = 1)

              #grep values from heating element
              temp.sequence.object.curveValue.heating.element <-
                src_get_XSYG_curve_values(XML::xmlValue(record[[3]]))

              temp.element <- temp.sequence.object.curveValue.heating.element[, 1]

              ## reduce matrix values to values of the detection
              temp.sequence.object.curveValue.heating.element <-
                  temp.sequence.object.curveValue.heating.element[
                      temp.element >= min(time.values) &
                      temp.element <= max(time.values), ,
                      drop = FALSE]

              ## calculate corresponding heating rate, this makes only sense
              ## for linear heating, therefore it has to be the maximum value

              ##remove 0 values (not measured) and limit to peak
              heating.rate.values <- temp.sequence.object.curveValue.heating.element[
                temp.sequence.object.curveValue.heating.element[,2] > 0 &
                  temp.sequence.object.curveValue.heating.element[,2] <=
                  max(temp.sequence.object.curveValue.heating.element[,2]),,drop = FALSE]

              heating.rate <- (heating.rate.values[length(heating.rate.values[,2]), 2] -
                                 heating.rate.values[1,2])/
                (heating.rate.values[length(heating.rate.values[,1]), 1] -
                   heating.rate.values[1,1])

              ##round values
              heating.rate <- round(heating.rate, digits=1)

              ##add to info element
              temp.sequence.object.info <- c(temp.sequence.object.info,
                                             RATE = heating.rate)

              ##PERFORM RECALCULATION
              ##check which object contains more data
              if(!"Spectrometer" %in% temp.sequence.object.detector){
                ##CASE (1)
                if(nrow(temp.sequence.object.curveValue.PMT) >
                   nrow(temp.sequence.object.curveValue.heating.element)){
                  temp.sequence.object.curveValue.heating.element.i <- approx(
                    x = temp.sequence.object.curveValue.heating.element[,1],
                    y = temp.sequence.object.curveValue.heating.element[,2],
                    xout = time.values,
                    rule = 2)

                  temperature.values <-
                    temp.sequence.object.curveValue.heating.element.i$y

                  count.values <-
                    temp.sequence.object.curveValue.PMT[,2]

                  ##CASE (2)
                }else if((nrow(temp.sequence.object.curveValue.PMT) <
                          nrow(temp.sequence.object.curveValue.heating.element))){
                  temp.sequence.object.curveValue.PMT.i <- approx(
                    x = time.values,
                    y = temp.sequence.object.curveValue.PMT[,2],
                    xout = temp.sequence.object.curveValue.heating.element[,1],
                    rule = 2)

                  temperature.values <-
                    temp.sequence.object.curveValue.heating.element[,2]

                  count.values <- temp.sequence.object.curveValue.PMT.i$y

                  ##CASE (3)
                }else{
                  temperature.values <-
                    temp.sequence.object.curveValue.heating.element[,2]

                  count.values <- temp.sequence.object.curveValue.PMT[,2]
                }

                ##combine as matrix
                curve <- cbind(temperature.values, count.values)

                ##set curve identifier
                temp.sequence.object.info$curveDescripter <- "Temperature [\u00B0C]; Counts [a.u.]"

              }else{
                ##CASE (1) here different approach. in contrast to the PMT measurements, as
                ##         usually the resolution should be much, much lower for such measurements
                ##         Otherwise we would introduce some pseudo signals, as we have to
                ##         take care of noise later one

                if (length(time.values) !=
                    nrow(temp.sequence.object.curveValue.heating.element)) {
                  temp.sequence.object.curveValue.heating.element.i <- approx(
                    x = temp.sequence.object.curveValue.heating.element[,1],
                    y = temp.sequence.object.curveValue.heating.element[,2],
                    xout = time.values,
                    rule = 2,
                    ties = mean,
                    na.rm = FALSE)

                  temperature.values <-
                    temp.sequence.object.curveValue.heating.element.i$y

                  ##check for duplicated values and if so, increase this
                  idx.dup <- which(duplicated(temperature.values))
                  if (length(idx.dup) > 0) {
                    temperature.values[idx.dup] <- temperature.values[idx.dup] + 1
                    .throw_warning("Temperature values are found to be ",
                                   "duplicated and increased by 1 K.")
                  }
                  ##CASE (2)  (equal)
                }else{
                  ##CASE (2)  (equal)
                  temperature.values <-
                    temp.sequence.object.curveValue.heating.element[,2]
                }

                ##reset values of the matrix
                colnames(temp.sequence.object.curveValue.spectrum) <- temperature.values
                curve <- temp.sequence.object.curveValue.spectrum

                ##change curve descriptor
                temp.sequence.object.info$curveDescripter <- "Temperature [\u00B0C]; Wavelength [nm]; Counts [1/ch]"
              }
            }##endif
          }##endif recalculate.TL.curves == TRUE


          ## Cleanup info objects -------------------------------------------
          if ("curveType" %in% names(temp.sequence.object.info))
            temp.sequence.object.info[["curveType"]] <- NULL

          ## Set RLum.Data-objects ------------------------------------------
          if (!"Spectrometer" %in% temp.sequence.object.detector) {
            out.class <- "RLum.Data.Curve"
            if (!inherits(curve, "matrix")) {
              curve <- src_get_XSYG_curve_values(XML::xmlValue(curve))
            }
          } else {
            out.class <- "RLum.Data.Spectrum"
            if (!inherits(curve, "matrix")) {
              curve <- get_XSYG.spectrum.values(curve)
            }
          }

          set_RLum(
              class = out.class,
              originator = "read_XSYG2R",
              recordType = paste0(recordType,
                                  " (", temp.sequence.object.detector,")"),
              curveType = temp.sequence.object.curveType,
              data = curve,
              info = temp.sequence.object.info)
        })
    }), use.names = FALSE)

      ##if the XSYG file is broken we get NULL as list element
      if (!is.null(sequence.object)) {
        ##set RLum.Analysis object
        sequence.object <- set_RLum(
          originator = "read_XSYG2R",
          class = "RLum.Analysis",
          records = sequence.object,
          protocol = as.character(sequence.header["protocol", 1]),
          info = list(file = file)
        )

        ##set parent uid of RLum.Anlaysis as parent ID of the records
        sequence.object <- .set_pid(sequence.object)

        ##update progress bar
        if (txtProgressBar) {
          setTxtProgressBar(pb, x)
        }

        ##merge output and return values
        if(fastForward){
          return(sequence.object)
        }
        return(list(Sequence.Header = sequence.header,
                    Sequence.Object = sequence.object))
      }
  }) ##end loop for sequence list

  ## close ProgressBar
  if (txtProgressBar)
    close(pb)

  ## show output information
  num.removed <- length(output[sapply(output, is.null)])
  if (num.removed > 0)
    .throw_warning(num.removed, " incomplete sequence(s) removed")

  if (verbose) {
    cat(sprintf("\t >> %d of %d sequence(s) loaded successfully\n",
                n_records - num.removed, n_records))
  }

  ##get rid of the NULL elements (as stated before ... invalid files)
  output <- .rm_NULL_elements(output)

  ## account for linearity
  if (auto_linearity_correction) {
    ## set look-up table for common FI PMTs
    count_pair_res <- c(
      UVVIS = 18,
      NIR50 = 70,
      NIR40 = 70,
      ETPMT = 25)

    ## correct only curves that can be corrected
    output <- lapply(output, \(x) {
      x@records <- lapply(x@records, \(y) {
        detector <- y@info$detector
        if (!is.null(detector) && any(detector %in% names(count_pair_res))) {
          y <- correct_PMTLinearity(y, PMT_pulse_pair_resolution = count_pair_res[detector])
        }
        return(y)
      })
      return(x)
    })
  }

  ## return object
  return(output)
}