find_centers.R

#' Initialize centroids for `cluster_xmap()`
#' 
#' @param xmap an `qm_xmap` class object generated by `read_xmap()``
#' @param qnt an `qm_qnt` class object generated by `read_qnt()``
#' @param fine_phase A character vector to specify fine grained phases which tend to comprise multi-phase pixels e.g., c('Pl', 'Amp')
#' @param saveas File name to save result. FALSE if not saving.
#'
#' @importFrom data.table fwrite
#' @importFrom dplyr bind_cols
#' @importFrom dplyr bind_rows
#' @importFrom dplyr filter
#' @importFrom dplyr group_by
#' @importFrom dplyr mutate
#' @importFrom dplyr summarise
#' @importFrom dplyr ungroup
#' @importFrom pipeR pipeline
#' @importFrom stats lsfit
#' @importFrom stats qnbinom
#' @importFrom stats median
#' @importFrom tidyr spread
#'
#' @export
find_centers <- function(
  xmap,
  qnt,
  fine_phase = NULL,
  saveas = 'centers0.csv'
) {

  quantified <- names(xmap) %in% qnt$elm$elint
  xmap_df <- as.data.frame(lapply(xmap, unlist, use.names = FALSE))
  
  pipeline({
    tidy_epma(qnt = qnt, xmap = xmap, cluster = NULL)
    group_by(phase)
    filter(all(is.na(map)) | !is.na(map))
    ungroup
    mutate( # Let weighting for least squares to be 0 for phases those who listed in fine_phase
      w = if(is.null(fine_phase)) 1 else as.integer(phase %nin% fine_phase)
    )
    group_by(elint) # Peform least squares and estimate 99% prediction interval
    mutate(
      .tmp = is.finite(map),
      map_est =
        pkint * lsfit(pkint[.tmp], map[.tmp], w[.tmp], intercept = FALSE)$coef,
      map_fnite = NULL, 
      pi_L = qnbinom(0.005, map_est, 0.5),
      pi_H = qnbinom(0.995, map_est + 1, 0.5)
    )
    group_by(id)
    mutate(within_pi = all((pi_L <= map) & (map <= pi_H)))  # pi = prediction interval
    group_by(elint, phase)
    mutate(n_within_pi = sum(within_pi))
    ungroup
    select(elint, map, map_est, within_pi, n_within_pi, phase, id, nr)
    bind_rows(
      if(!all(quantified)) {
        pipeline({
          .
          filter(elint == elint[1])
          select(-elint, -map, -map_est)
          bind_cols(
            lapply(xmap_df[!quantified], `[`, .$nr)
          )
          gather(
            elint, map, -within_pi, -n_within_pi, -phase, -id, -nr
          )
          mutate(map_est = map)
        })
      }
    )
    group_by(phase)
    mutate(not_mapped_phase = all(is.na(map)))
    ungroup
    mutate(map = ifelse(not_mapped_phase | within_pi == 0, map_est, map))
    filter(not_mapped_phase | within_pi | n_within_pi == 0)
    group_by(elint, phase)
    summarise(map = median(map))
    ungroup
    spread(elint, map)
    # guess mapping intensity of certain phases in case
    # target elements are not analyzed by reference point analysis
    x ~ { 
      miss <- Reduce(`|`, lapply(x, is.na))
      if(all(!miss)) return(x)
      x[miss, -1] <- pipeline({
        map2(
          list(t(xmap_df[names(x[-1])])),
          pmap(x[miss, -1], c),
          `-`
        )
        map(square)
        map(colSums, na.rm = TRUE)
        map(which.min)
        map(function(i) xmap_df[i, ])
        bind_rows
        `[`(names(x[-1]))
        map2(x[miss, -1], function(y, x) ifelse(is.na(x), y, x))
        bind_cols
      })
      x
    }
    save4qm(nm = saveas, saving = is.character(saveas))
  })

}

#' (Deprecated) Use find_centers
#' @param wd working directory which contains .qnt directory. Use current directory if NULL
#' @param dir_map path to the directory which contains mapping data.
#' @param phase_fine fine grained phases which tend to be appear in multi-phase pixels e.g., c('Pl', 'Amp')
#' @param qnt an object generated by qnt_load or a path to the .RDS file created by qnt_load
#' @param qltmap an object generated by qltmap_load or a path to the .RDS file created by qltmap_load
#' @param saving file name to save. FALSE if not saving.
#' @export
qltmap_cls_centers <- function(
  wd = NULL,
  dir_map,
  phase_fine = NULL,
  qnt = read_qnt(paste0(wd, '/.qnt')),
  qltmap = read_xmap(dir_map),
  saving = 'centers_initial0.csv'
) {
  .Deprecated(new = 'find_centers')
  cd <- getwd()
  on.exit(cd)
  setwd(wd)
  find_centers(fine_phase = phase_fine, qnt = qnt, xmap = qltmap, saveas = saving)
}