R/proximate_read_cal.R

In proximetricsR: Spectral Preprocessing and Chemometric Calibration of NIR Sensors

#' @title Read model parameters from ProxiMate .cal files
#' @aliases proximate_read_cal
#' @aliases predict.read_cal
#' @description
#'
#' Reads the metadata and model parameters from one or more \code{.cal} files
#' generated by BUCHI ProxiMate sensors. The function extracts the preprocessing
#' recipe, regression method, PLS weights, loadings, scores, intercepts, and
#' bias terms required to project new spectra into the score space and produce
#' predictions. Spectral regression coefficients are not retrieved directly;
#' predictions are computed in the score space via \code{\link{predict.read_cal}}.
#'
#' @usage
#' proximate_read_cal(file, ignore_version = FALSE)
#'
#' \method{predict}{read_cal}(object, newdata, get_comp = c("optimal", "all"),
#'         get_scores = FALSE, bias_index = 1, ...)
#'
#' @param file a character vector of \code{.cal} file paths.
#' @param ignore_version a logical. If \code{FALSE} (default), files with no
#' version information or created with NIRWise PLUS prior to version 1.0 raise
#' an error. If \code{TRUE}, such files are read with a warning instead;
#' predictions from these files may deviate from those produced on the instrument.
#' @param object an object of class \code{read_cal} as returned by
#' \code{proximate_read_cal()}.
#' @param newdata a matrix of new spectral data to predict from. Column names
#' must be coercible to the wavelengths used in the model.
#' @param get_comp a character string. Either \code{"optimal"} (default) to
#' return predictions only for the optimal number of components, or \code{"all"}
#' to return predictions for every available component.
#' @param get_scores a logical indicating whether PLS scores should be returned
#' alongside predictions. Default is \code{FALSE}.
#' @param bias_index the index of the bias to be applied in the list of biases.
#' These are generated in NIRWise PLUS based on the number of files containing
#' the calibration data. Default = 1.
#' @param ... not currently used.
#'
#' @return
#'
#' For \code{proximate_read_cal()}, a list of class \code{"read_cal"} with the following
#' elements:
#'
#' \itemize{
#'     \item \strong{\code{summary}:} a data.frame describing each model:
#'     \itemize{
#'          \item \strong{\code{Property}:} name of the response variable.
#'          \item \strong{\code{Preprocessing}:} sequence of preprocessing steps applied (without parameters).
#'          \item \strong{\code{Method}:} regression method used.
#'          \item \strong{\code{Factors}:} number of PLS components used.
#'          \item \strong{\code{Cross-validation}:} number of cross-validation segments.
#'          A value of 0 indicates no cross-validation was used.
#'          \item \strong{\code{Auto-skip}:} logical indicating whether automatic outlier
#'          removal (auto-delete) was applied during calibration.
#'       }
#'     \item \strong{\code{meta_param}:} a list with one element per model containing
#'     the preprocessing recipe (\code{precipe}), the indices of automatically
#'     removed observations (\code{auto_skip}), and a logical indicating whether
#'     sample aggregation was applied (\code{aggregate}).
#'     \item \strong{\code{file_info}:} a list with one element per model containing
#'     the file paths of the spectral data used for calibration (\code{files})
#'     and the indices of manually skipped observations per file
#'     (\code{skipped_indices}).
#'     \item \strong{\code{models}:} a list with one element per model containing all
#'     parameters required for prediction: wavelengths, preprocessing recipe,
#'     number of factors, mean-centering vector, scores, score scale factors,
#'     PLS weights, loadings, biases, intercept, and target values.
#' }
#'
#' For \code{predict.read_cal()}, a list with the following elements:
#'
#' \itemize{
#'     \item \strong{\code{predictions}:} predicted values for each model in
#'     \code{object}.
#'     \item \strong{\code{distances}:} scaled score distances for each sample and
#'     model, which can be used to assess how well a new sample is represented
#'     by the model.
#'     \item \strong{\code{scores}:} only returned when \code{get_scores = TRUE}. The
#'     projection of new samples into the PLS score space.
#' }
#'
#' @seealso
#' \code{\link{proximate_recalibrate_nax}},
#' \code{\link{proximate_read_nax}}
#'
#' @author Leonardo Ramirez-Lopez and Claudio Orellano
#' @export
proximate_read_cal <- function(file, ignore_version = FALSE) {
  fformat <- tools::file_ext(file)
  fformat <- unique(fformat)

  allowedformat <- c("cal")
  # res <- 2

  if (length(fformat) > 1 || !all(fformat %in% allowedformat)) {
    stop("Ivalid input file format")
  }

  if (fformat == "cal") {
    nms <- c(
      "Property",
      "Preprocessing",
      "Method",
      "Factors",
      "Cross-validation",
      "Auto-skip"
    )

    cal_summary <- data.frame(matrix(NA, length(file), length(nms)))
    colnames(cal_summary) <- nms
    models <- file_info <- meta_param <- NULL
    for (i in file) {
      is_protected <- grepRaw(
        "^BUCHI [A-Z a-z]{1,}Licensed",
        readBin(i, what = "raw", n = 50)
      ) |> length() > 0

      if (is_protected) {
        message("Protected cal file(s). These will not be read.")
        return(list(
          summary = "Protected",
          meta_param = NULL
        ))
        break
      }

      bf <- readLines(i, warn = FALSE, encoding = "UTF-8") |>
        iconv("windows-1252", "UTF-8")

      p_line <- grep("Selector1:Y", bf)
      p_line <- p_line[length(p_line)]
      property <- strsplit(bf[p_line], "\t| ")[[1]][2]
      regression <- strsplit(bf[grep("Model1:Type", bf)], "\t| ")[[1]][2]

      # extract file paths
      tsvs <- bf[grep("Files:File", bf)]
      # tsvs <- tsvs[grep("\\.tsv", tsvs)]
      # tsvs <- tsvs[grep("#TRUE#", tsvs)]
      # tsvs <- gsub("\\\"", "", tsvs)
      tsvs <- strsplit(tsvs, "\t", fixed = TRUE)
      tsvs <- do.call("rbind", tsvs)
      tsvs <- gsub("Files:", "", tsvs)
      tsvs <- gsub("^\\.", "", tsvs)
      tsvs <- gsub("\\", "/", tsvs, fixed = TRUE)
      tsvs <- as.data.frame(tsvs)
      colnames(tsvs) <- c("File", "Path")
      file_info[[which(i == file)]] <- list(files = tsvs)

      cv <- bf[grep("Model1:CrossValidate", bf)]
      if (length(cv) != 0) {
        cv <- as.numeric(strsplit(cv, "\t| ")[[1]][2])
      } else {
        cv <- 0
      }

      factors <- as.numeric(strsplit(bf[grep("Model1:Factors", bf)], "\t| ")[[1]][2])

      zero <- bf[grep("Zero", bf)] |> strsplit(split = "\t| |,")
      zero <- zero[[1]][-1] |> as.numeric()


      scores_idx <- grep("Scores\t", bf)

      if (length(scores_idx) > 1) {
        scores_idx <- scores_idx[grep("^Scores", bf[scores_idx])]
      }

      scores <- bf[scores_idx] |> strsplit(split = "\t| |;")
      scores <- scores[[1]][-1]
      scores <- lapply(seq_along(scores),
        FUN = function(x, i) {
          as.numeric(strsplit(x[i], split = "\t| |,")[[1]])
        },
        x = scores
      )
      scores <- do.call("rbind", scores) |> t()


      pls_weights <- bf[grep("Weights\t", bf)[1]] |> strsplit(split = "\t| |;")
      pls_weights <- pls_weights[[1]][-1]
      pls_weights <- lapply(seq_along(pls_weights),
        FUN = function(x, i) {
          as.numeric(strsplit(x[i], split = "\t| |,")[[1]])
        },
        x = pls_weights
      )
      pls_weights <- do.call("rbind", pls_weights)

      loadings <- bf[grep("Loads\t", bf)[1]] |> strsplit(split = "\t| |;")
      loadings <- loadings[[1]][-1]
      loadings <- lapply(seq_along(loadings),
        FUN = function(x, i) {
          as.numeric(strsplit(x[i], split = "\t| |,")[[1]])
        },
        x = loadings
      )
      loadings <- do.call("rbind", loadings)

      biases <- bf[grep("^Bias\t", bf)[1]] |> strsplit(split = "\t| |;")
      biases <- biases[[1]][-1]

      targets <- bf[grep("^Target\t", bf)[1]] |> strsplit(split = "\t| |,")
      targets <- targets[[1]][-1] |> as.numeric()

      biases <- lapply(seq_along(biases),
        FUN = function(x, i) {
          as.numeric(strsplit(x[i], split = "\t| |,")[[1]])
        },
        x = biases
      )
      biases <- do.call("rbind", biases)


      center <- bf[grep("^Center\t", bf)[1]] |> strsplit(split = "\t| |;")
      center <- center[[1]][-1] |> as.numeric()

      wavelengths <- bf[grep("^Pretreat[0-9].Wavelengths\t", bf)[1]] |> strsplit(split = "\t| |,")
      wavelengths <- wavelengths[[1]][-1] |> as.numeric()

      scale_scores <- bf[grep("^Scale\t", bf)[1]] |> strsplit(split = "\t| |,")
      scale_scores <- scale_scores[[1]][-1] |> as.numeric()


      colnames(loadings) <- colnames(pls_weights) <- names(zero) <- wavelengths
      if (nrow(biases) == 1 & ncol(biases) == ncol(scores)) {
        biases <- t(biases)
      }
      biases <- biases[seq_len(ncol(scores)), , drop = FALSE]
      rownames(biases) <- colnames(scores) <- rownames(pls_weights) <- rownames(loadings) <- paste0("ncomp_", seq_len(nrow(loadings)))
      colnames(biases) <- seq_len(ncol(biases))
      auto_skip <- bf[grep("Model[0-9]+:AutoDelete", bf)]
      manual_skip <- bf[grep("Model[0-9]+:Delete", bf)]
      manual_skip <- manual_skip[-grep("Auto", manual_skip)]
      manual_skip <- gsub("Model1:Delete[0-9]+\t", "", manual_skip)

      skipped <- sapply(
        tsvs$File,
        FUN = function(x) as.numeric(),
        simplify = TRUE,
        USE.NAMES = TRUE
      )

      if (length(manual_skip) == 1) {
        if (manual_skip == "<new deletes>") {
          manual_skip <- NULL
        }
      }


      if (length(manual_skip) > 0) {
        manual_skip <- lapply(
          manual_skip,
          FUN = function(x) {
            strsplit(x, "\t|,")[[1]]
          }
        )
        manual_skip <- do.call("c", manual_skip)
        # file_idx <- gsub("File", "", tsvs$File)
        # paste0(file_idx, collapse = "|")
        manual_skip <- do.call("rbind", strsplit(manual_skip, ".", fixed = TRUE))
        # skipped <- list()
        for (j in unique(manual_skip[, 1])) {
          skipped[[paste0("File", j)]] <- as.numeric(manual_skip[manual_skip[, 1] == j, 2]) + 1
        }
      }

      file_info[[which(i == file)]]$skipped_indices <- skipped


      if (length(auto_skip) != 0) {
        auto_skip <- as.numeric(strsplit(auto_skip, "\t| ")[[1]][-1])
      } else {
        auto_skip <- NULL
      }

      version_line <- grep("^Version\t", bf, value = TRUE, perl = TRUE)
      has_version <- length(version_line) > 0
      supported <- has_version && any(grepl("^Version\t[1-9][0-9]*\\.", bf, perl = TRUE))

      if (!supported) {
        if (!ignore_version) {
          if (has_version) {
            stop(
              "File '", basename(i), "' was created with NIRWise PLUS prior to ",
              "version 1.0 and is not supported. Re-save the model in NIRWise PLUS ",
              "version 1.0 or later, or set ignore_version = TRUE to read it anyway."
            )
          } else {
            stop(
              "File '", basename(i), "' has no version information and cannot be read. ",
              "Set ignore_version = TRUE to attempt reading it anyway."
            )
          }
        } else {
          if (has_version) {
            warning(
              "File '", basename(i), "' was created with NIRWise PLUS prior to version 1.0. ",
              "Predictions may deviate from those produced on the instrument.",
              call. = FALSE
            )
          } else {
            warning(
              "File '", basename(i), "' has no version information and its origin is unknown. ",
              "Predictions may deviate from those produced on the instrument.",
              call. = FALSE
            )
          }
        }
      }

      pret <- grep(
        "Pretreat1(?!.*Means)(?!.*Wavelengths)(?!.*Targets)",
        bf,
        perl = TRUE
      )
      pret <- bf[pret]
      ## FIXME: include other pretreatments of NIRWise PLUS and SX PLUS
      pspline <- grep("Spline", pret)
      if (length(pspline) > 1) {
        pspline <- pspline[1] # FIXME: it can be various splines in a file?? :/
      }
      psnvt <- grep("SNVT", pret)
      pmsc <- grep("MSC", pret)
      pder <- grep("DG", pret)
      psmooth <- grep("SMOOTH", pret, ignore.case = TRUE)
      porder <- order(c(pspline, psnvt, pmsc, pder, psmooth))

      do_aggregation <- any(grepl("AVG", pret, ignore.case = TRUE))

      preprocessing <- NULL
      if (length(pspline) > 0) {
        pspline <- strsplit(pret[pspline], "\t| ")[[1]]
        pspline <- as.numeric(pspline[(length(pspline) - 2):length(pspline)])
        preprocessing[["spline"]] <- prep_resample(
          grid = c(pspline[1], pspline[2], pspline[3])
        )
      }

      if (length(psnvt) > 0) {
        preprocessing[["snvt"]] <- prep_snv()
      }

      if (length(pmsc) > 0) {
        stop(paste0("Sorry, '", basename(i), "' model with MSC. This pretreatment is not yet functional... coming soon"))
      }

      if (length(pder) > 0) {
        pder <- strsplit(pret[pder], "\t| ")[[1]]
        pder <- as.numeric(pder[(length(pder) - 2):length(pder)])
        preprocessing[["der"]] <- prep_derivative(
          m = pder[1],
          w = as.integer(pder[2]) * 2L - 1L,
          p = as.integer(pder[3]) * 2L + 1L,
          algorithm = "nwp"
        )
      }

      if (length(psmooth) > 0) {
        psmooth <- strsplit(pret[psmooth], "\t| ")[[1]]
        psmooth <- as.numeric(psmooth[(length(psmooth) - 2):length(psmooth)])
        preprocessing[["smooth"]] <- prep_smooth(
          w = as.integer(psmooth[1]) * 2L + 1L,
          algorithm = "moving-average"
        )
      }

      preprocessing <- preprocessing[porder]
      preprocessing$device <- "proximate"
      mrecipe <- do.call(preprocess_recipe, preprocessing)
      if (length(mrecipe) == 0) {
        mrecipe$preprocessing_order <- "none"
      }
      cal_summary[which(i == file), "Property"] <- property
      cal_summary[which(i == file), "Preprocessing"] <- mrecipe$preprocessing_order
      cal_summary[which(i == file), "Method"] <- regression
      cal_summary[which(i == file), "Factors"] <- factors
      cal_summary[which(i == file), "Cross-validation"] <- cv
      cal_summary[which(i == file), "Auto-skip"] <- !is.null(auto_skip)

      meta_param[[which(i == file)]] <- list(
        precipe = mrecipe,
        auto_skip = auto_skip,
        aggregate = do_aggregation
      )

      models[[which(i == file)]] <- list(
        Wavelengths = wavelengths,
        Preprocessing = mrecipe,
        Factors = factors,
        Zero = zero,
        Scores = scores,
        Scale_scores = scale_scores,
        Weights = pls_weights,
        Loadings = loadings,
        Biases = biases,
        Intercept = center,
        Targets = targets
      )
    }

    names(models) <- names(meta_param) <- cal_summary$Property

    mcalibration <- list(
      summary = cal_summary,
      meta_param = meta_param,
      file_info = file_info,
      models = models
    )
    class(mcalibration) <- c("read_cal", "list")
    return(mcalibration)
  }
}


#' @aliases proximate_read_cal
#' @export
predict.read_cal <- function(object, newdata, get_comp = c("optimal", "all"),
                             get_scores = FALSE, bias_index = 1, ...) {
  get_comp <- match.arg(get_comp)
  nms <- names(object$models)

  scores <- results_mahal <- results <- list()
  for (i in seq_along(object$models)) {
    spc_pp <- process(newdata, object$models[[i]]$Preprocessing)

    if (get_comp == "optimal") {
      maxf <- object$models[[i]]$Factors
    }
    if (get_comp == "all") {
      maxf <- ncol(object$models[[i]]$Scores)
    }

    spc_pp <- spc_pp[, colnames(spc_pp) %in% object$models[[i]]$Wavelengths, drop = FALSE]

    zrv <- which(colSums(abs(object$models[[i]]$Loadings) == 0) == nrow(object$models[[i]]$Loadings))
    if (length(zrv) > 0) {
      # this fixes issue # 20
      if (all(!names(zrv) %in% colnames(spc_pp))) {
        addm <- matrix(0, nrow(spc_pp), length(zrv), dimnames = list(rownames(spc_pp), names(zrv)))
        spc_pp <- cbind(addm, spc_pp)[, colnames(object$models[[i]]$Loadings), drop = FALSE]
      }
    }

    spc_cent <- sweep(spc_pp, 2, FUN = "-", STATS = object$models[[i]]$Zero)
    # projection_matrix <- t(object$models[[i]]$Weights) %*% solve(object$models[[i]]$Loadings %*% t(object$models[[i]]$Weights))
    # spc_cent %*% projection_matrix

    ij_scores <- matrix(NA, nrow(spc_pp), maxf)
    for (j in 1:maxf) {
      ij_scores[, j] <- spc_cent %*% object$models[[i]]$Weights[j, ]
      # FIXME: check the correct BIAS to take from the columns of object$models[[i]]$Biases
      ij_scores[, j] <- ij_scores[, j] + sum(object$models[[i]]$Biases[j, bias_index])
      spc_cent <- spc_cent - (ij_scores[, j, drop = FALSE] %*% object$models[[i]]$Loadings[j, , drop = FALSE])
    }

    p_val <- object$models[[i]]$Intercept + t(apply(ij_scores, MARGIN = 1, FUN = cumsum))

    colMeans(object$models[[i]]$Scores)
    sc <- sweep(
      ij_scores,
      MARGIN = 2,
      FUN = "/",
      STATS = object$models[[i]]$Scale_scores[seq_len(ncol(ij_scores))]
    )

    mahal <- sc * NA
    for (j in 1:maxf) {
      mahal[, j] <- rowMeans(sc[, seq_len(j), drop = FALSE]^2)
    }

    if (get_comp == "optimal") {
      p_val <- p_val[, maxf, drop = FALSE]
      mahal <- mahal[, maxf, drop = FALSE]
      ith_names <- paste0("ncomp_", maxf)
    } else {
      ith_names <- paste0("ncomp_", seq_len(maxf))
    }

    colnames(mahal) <- colnames(p_val) <- ith_names
    rownames(mahal) <- rownames(p_val) <- seq_len(nrow(p_val))

    results_mahal[[i]] <- mahal
    results[[i]] <- p_val
    if (get_scores) {
      dimnames(ij_scores) <- list(
        seq_len(nrow(ij_scores)),
        paste0("ncomp", seq_len(col(ij_scores)))
      )
      scores[[i]] <- ij_scores
    }
  }
  names(results) <- names(results_mahal) <- nms
  preds <- list(
    predictions = results,
    distances = results_mahal
  )
  if (get_scores) {
    names(scores) <- nms
    preds$scores <- scores
  }
  return(preds)
}

#' @title Calculate the ASCII Value Difference Between Two Strings
#' @description
#' This function calculates the difference between two strings based on their
#' ASCII values.
#' @param s1 A character string.
#' @param s2 A character string.
#' @details
#' The function ensures that both strings are of the same length by padding them
#' with spaces if necessary. It then computes the difference between the ASCII
#'  values of corresponding characters in the strings.
#' @return The absolute difference between the ASCII values of the characters
#' in the two strings.
#' @keywords internal
string_diff <- function(s1, s2) {
  # some examples
  # proximetricsR:::string_diff("abc", "abd")
  # proximetricsR:::string_diff("1001", "1000") == proximetricsR:::string_diff("1000", "1001")
  # proximetricsR:::string_diff("1010", "1000") == 256
  # proximetricsR:::string_diff("1100", "1000") == 256^2
  # proximetricsR:::string_diff("1111", "1000") == 256^2 + 256 + 1
  # proximetricsR:::string_diff("100a", "1000") == 49
  # proximetricsR:::string_diff("100", "1000") == 16
  # proximetricsR:::string_diff("2000", "1000") == 256^3
  # proximetricsR:::string_diff("2000", "1001") == 256^3 - 1

  # Ensure that the strings are of the same length
  max_len <- max(nchar(s1), nchar(s2))
  s1 <- sprintf("%-*s", max_len, s1)
  s2 <- sprintf("%-*s", max_len, s2)

  # Initialize the result
  result <- 0

  # Calculate the difference between the ASCII values of the characters
  for (i in seq_len(max_len)) {
    # We multiply by 256 here to ensure that each character's difference
    # occupies
    # a different byte in the final number.
    result <- result * 256 + (utf8ToInt(substring(s1, i, i)) - utf8ToInt(substring(s2, i, i)))
  }

  abs(result)
}

#' @title Function to locate the serial number index
#' @description
#' This function locates the index of a given serial number in a list of serial
#' numbers.
#' @param serial_numbers A vector of serial numbers.
#' @param serialnumber A single serial number to be located within the list.
#' @details
#' The function first checks if the list of serial numbers is empty and
#' returns 1 if true. If the serial number is found in the list, it returns the
#' last index where the serial number appears. If the serial number is not
#' found, it calculates the ASCII value difference between the given serial
#' number and each element in the list using the `string_diff` function. If
#' all differences are larger than 256^3, it returns 1. Otherwise, it returns
#'  the index of the element with the smallest ASCII value difference.
#' @return The index of the serial number in the list, or the index of the
#' closest match based on ASCII value difference.
#' @keywords internal


locate_serialnumber_index <- function(serial_numbers, serialnumber) {
  # EXAMPLES:
  #   serial_numbers <- c("12345", "67890", "ABCDE")
  #   locate_serialnumber_index(serial_numbers, "12345")
  #   serial_numbers <- c("1001", "1000", "1010", "1100", "1111", "100a", "100", "2000")
  #   locate_serialnumber_index(serial_numbers, "1001") # Should return 1
  #   locate_serialnumber_index(serial_numbers, "2000") # Should return 8
  #   # Should return the index with smallest difference, 6
  #   locate_serialnumber_index(serial_numbers, "100b")

  # seealso the following function: string_diff()

  # Return 1 if the list is empty
  if (length(serial_numbers) == 0) {
    return(1)
  }

  # Return the last index if the serialnumber is found in the list
  if (serialnumber %in% serial_numbers) {
    return(max(which(serial_numbers == serialnumber)))
  }

  # Calculate the string difference between the serialnumber and each element
  # in the list
  differences <- sapply(serial_numbers, string_diff, s2 = serialnumber)

  # If all differences are larger than 256^3, return 1
  if (all(differences >= 256^3)) {
    return(1)
  }

  # Return the index of the element with the smallest difference in terms of
  # ASCII values
  min_diff <- min(differences)
  return(max(which(differences == min_diff)))
}

#' @noRd
#' @export
print.proximate_read_cal <- function(x, ...) {
  cat(.bold_italic("BUCHI ProxiMate calibration model imported with"), "proximate_read_cal()\n")
  cat("---\n")
  print(x$summary)
  invisible(x)
}