R/create.DN.R

In datanugget: Create, Optimize, and Refine Data Nuggets

# This function creates the datanuggets.
# It returns a list of two objects:
# 1. Data Nuggets is a dataframe with cols DN no., DN centers, Scale and Weight.
# 2. Data Nugget Assignments is the vector containing the data nugget assignment of each observation in x.


# Function inputs
# x: original dataset.
# center.method: the method used for choosing data nugget centers.
# R: the number of observations to sample from the data matrix when creating the initial data nugget centers.
# delete.percent: proportion of data points to be deleted at each iteration.
# DN.num1: the number of initial data nugget centers to create.
# DN.num2: the number of final data nuggets to create.
# dist.metric: pairwise distance measure (e.g., "euclidean" or "manhattan").
# seed: random seed for replication.
# no.cores: number of cores used for parallel processing.
# make.pbs: logical; whether to show a progress bar while the function runs.



create.DN = function(x,
                     center.method = "mean",
                     R = 5000,
                     delete.percent = .1,
                     DN.num1 = 10^4,
                     DN.num2 = 2000,
                     dist.metric = "euclidean", 
                     seed = 291102,
                     no.cores = (parallel::detectCores() - 1),
                     make.pbs = FALSE){
  
  
  ## -------------------------------
  ## Argument checks
  ## -------------------------------
  
  
  # check x
  if (!any(class(x) %in% c("matrix", "data.frame", "data.table"))){
    stop('x must be of class "matrix", "data.frame", or "data.table"')
  }
  
  
  
  # check center.method
  if (!(center.method %in% c("mean", "random", "original"))){
    stop('center.method must be "mean" or "random" or "original"')
  }
  
  
  
  # check R
  if (!(class(R) %in% c("numeric", "integer"))){
    stop('R must be of class "numeric" or "integer"')
  }
  
  
  
  # make sure R is between 100 and 10000
  if (R < 100 | R > 10000){
    stop("R must be within [100,10000]")
  }
  
  
  
  # check delete.percent 
  if (!is.numeric(delete.percent)){
    stop('delete.percent must be of class "numeric"')
  }
  
  
  
  # make sure delete.percent is between 0 and 1
  if (delete.percent <= 0 | delete.percent >= 1){
    stop("delete.percent must be within (0,1)")
  }
  
  
  
  # check DN.num1
  if (!(class(DN.num1) %in% c("numeric", "integer"))){
    stop('DN.num1 must be of class "numeric" or "integer"')
  }
  
  
  
  # check DN.num2 
  if (!(class(DN.num2) %in% c("numeric", "integer"))){
    stop('DN.num2 must be of class "numeric" or "integer"')
  }
  
  
  
  # make sure DN.num2 <= DN.num1
  if (DN.num2 > DN.num1){
    stop("DN.num2 must be less than DN.num1")
  }
  
  
  
  # check dist.metric 
  if (dist.metric != "euclidean" & dist.metric != "manhattan"){
    stop('dist.metric must be "euclidean" or "manhattan"')
  }
  
  
  
  # check seed 
  if (!is.numeric(seed)){
    stop('seed must be of class "numeric"')
  }
  
  
  
  # check no.cores 
  if (!(class(no.cores) %in% c("numeric", "integer"))){
    stop('no.cores must be of class "numeric" or "integer"')
  }
  
  
  
  # check make.pbs 
  if (!is.logical(make.pbs)){
    stop("make.pbs must be TRUE OR FALSE")
  }
  
  
  
  ## ----------------------------------
  ## Pre processing and Initialization
  ## ----------------------------------
  

  # convert the input to matrix
  x_mat <- as.matrix(x)
  
  
  # check x elements
  if (!is.numeric(x_mat)) stop("x must contain only numeric columns")
  
  
  # storage mode set to double
  storage.mode(x_mat) <- "double"
  
  
  # number of observations in the data
  obs.num <- nrow(x_mat)
  
  
  # number of columns in the data
  n.col <- ncol(x_mat)
  
  
  # check if x has atleast two observations
  if (obs.num < 2) stop("x must contain at least two observations")
  
  
  # create random sample ordering
  set.seed(seed)
  RS.loc <- sample(obs.num, obs.num, replace = FALSE)

  
  # check if user wants to use parallel processing
  use.parallel <- no.cores > 1
  cl <- NULL
  
  if (use.parallel){
    
    # load the relevant packages
    for(pkg in c("parallel", "foreach")){
      if (!requireNamespace(pkg, quietly = TRUE)) {
        stop("Package '", pkg, "' is required for parallel processing.")
      }
    }
    
    
    have.snow <- requireNamespace("doSNOW", quietly = TRUE)
    have.doparallel <- requireNamespace("doParallel", quietly = TRUE)
    
    
    if (make.pbs && !have.snow){
      stop("Package 'doSNOW' is required for parallel progress bars. Set mak.pbs = FALSE to use doParallel.")
    }
    
    
    if (!have.snow && !have.doparallel){
      stop("Package 'doParallel' or 'doSNOW' is required for parallel processing.")
    }
    
    
    # create the cluster for parallel processing
    cl <- parallel::makeCluster(no.cores)
    parallel::clusterExport(cl, "create.DNcenters", envir = environment(create.DN))
    
    # stop cluster on exit
    on.exit(try(parallel::stopCluster(cl), silent = TRUE), add = TRUE)
    
    
    # engage the cluster for parallel processing
    if (make.pbs || !have.doparallel) doSNOW::registerDoSNOW(cl)
    else doParallel::registerDoParallel(cl)
    
  }
  
  `%dopar%` <- foreach::`%dopar%`
  
  ## ------------------------------------------
  ## Helper functions
  ## -------------------------------------------
  
  # create new environment
  clean.env <- new.env(parent = globalenv())
  
  # Function to call create.DNcenters and get the DN centers for each split
  get_centers_for_split <- function(tmp.RS, DN_num_for_split,
                                    delete.percent, dist.metric){
    
    
    out <- create.DNcenters(RS = tmp.RS, 
                            delete.percent = delete.percent, 
                            DN.num = DN_num_for_split, 
                            dist.metric = dist.metric, 
                            make.pbs = FALSE)
    
    attr(out, "kept") <- NULL
    return(as.data.frame(out))
  }
  
  environment(get_centers_for_split) <- clean.env
  
  
  # Function to assign each observation to the nearest DN center
  assign_chunk <- function(chunk, centers, metric){
    
    as.integer(Rfast::dista(chunk, centers, type = metric, k = 1L, index = TRUE))
    
  }
  
  environment(assign_chunk) <- clean.env
  
  
  ## ------------------------------------------
  ## Create initial DN centers of size DN.num1
  ## -------------------------------------------
  

  
  # split x into chunks to run across CPU cores 
  RS.splits <- ceiling(obs.num / R)
  chunks <- parallel::splitIndices(obs.num, RS.splits)
  chunks <- chunks[lengths(chunks) > 0L]
  RS.splits <- length(chunks)
  
  # divide DN.num1 centers evenly across the chunks 
  m.each <- rep(DN.num1 %/% RS.splits, RS.splits)
  rem <- DN.num1 - sum(m.each)
  if (rem > 0L) m.each[seq_len(rem)] <- m.each[seq_len(rem)] + 1L
  m.each <- pmax(1L, pmin(m.each, lengths(chunks)))
  
  
  # create initial set of data nugget centers
  message("Creating initial set of data nugget centers... ") 
  
  
  # check if user wants to use parallel processing
  if (use.parallel){
    
    opts <- NULL
    
    # check if user wants a progress bar
    if (make.pbs){
      
      # initialize progress bar
      pb <- txtProgressBar(min = 0, max = RS.splits)
      
      # close on exit
      on.exit(close(pb), add = TRUE)
      
      # update the progress bar
      progress <- function(n){utils::setTxtProgressBar(pb, n)}
      opts <- list(progress = progress)
      
    }
    
    # calculate the DN centers
    DN.centers1 <- foreach::foreach(
      i = seq_len(RS.splits),
      .combine = rbind,
      .inorder = TRUE,
      .export = c("create.DNcenters"),
      .packages = c("Rfast"),
      .options.snow = opts
    ) %dopar% {
      
      # reduce the batch to size DN.num1/RS.splits
      get_centers_for_split(tmp.RS = x_mat[RS.loc[chunks[[i]]], , drop = FALSE], 
                            DN_num_for_split = m.each[i], 
                            delete.percent = delete.percent, 
                            dist.metric = dist.metric)
    }
  }else {
    
    if (make.pbs) {
      pb <- utils::txtProgressBar(min = 0, max = RS.splits)
      on.exit(close(pb), add = TRUE)
    }
    
    parts <- vector("list", RS.splits)
    for (i in seq_len(RS.splits)) {
      
      parts[[i]] <- get_centers_for_split(
        tmp.RS = x_mat[RS.loc[chunks[[i]]], , drop = FALSE],
        DN_num_for_split = m.each[i],
        delete.percent = delete.percent,
        dist.metric = dist.metric)
        if (make.pbs) utils::setTxtProgressBar(pb, i)
    }
    DN.centers1 <- do.call(rbind, parts)
  }
  
  DN.centers1 <- as.matrix(DN.centers1)
  message("completed!")
  
  
  ## ---------------------------------------------------------
  ## Create final DN centers (reduce from DN.num1 to DN.num2)
  ## --------------------------------------------------------
  
  
  # create final set of data nugget centers(DN.num1 to DN.num2)
  message("Creating final set of data nugget centers...")
  
  
  DN.data <- create.DNcenters(RS = DN.centers1,
                              delete.percent = delete.percent,
                              DN.num = DN.num2,
                              dist.metric = dist.metric,
                              make.pbs = make.pbs)
  
  
  
  DN.data_mat <- as.matrix(DN.data)
  n.nuggets   <- nrow(DN.data_mat)
  
  message("completed!")
  
  
  
  ## -----------------------------------------
  ## Assign observations to the nearest DN
  ## -----------------------------------------
  
  message("Assigning observations to the nearest data nuggets...")
  
  if (use.parallel){
    
    n.chunks <- min(obs.num, max(1L, as.integer(no.cores) * 2L))
    idx <- parallel::splitIndices(obs.num, n.chunks)
    idx <- idx[lengths(idx) > 0L]
    res <- parallel::parLapply(cl, idx, function(ii)
      assign_chunk(x_mat[ii, , drop = FALSE], DN.data_mat, dist.metric))
    DN.assignments <- integer(obs.num)
    
    for (k in seq_along(idx)) DN.assignments[idx[[k]]] <- res[[k]]
  } else{
    
    DN.assignments <- assign_chunk(x_mat, DN.data_mat, dist.metric)
  }
  
  DN.assignments <- as.integer(DN.assignments)
  
  message("completed!")
  
  
  
  ## ---------------------------------------------------------------------------------------------------
  ## Recalculating the Data nugget centers based on the given method and calculating Data nugget scales
  ## ---------------------------------------------------------------------------------------------------
  
  
  message("Recalculating the data nugget centers and calculating the data nugget scales...")
  
  grp    <- factor(DN.assignments, levels = seq_len(n.nuggets))
  n.i    <- tabulate(DN.assignments, nbins = n.nuggets) # DN weights
  
  
  g.mean <- colMeans(x_mat)
  xc     <- sweep(x_mat, 2L, g.mean, "-")
  
  fill_groups <- function(S, k, p) {
    out <- matrix(0, nrow = k, ncol = p)
    out[as.integer(rownames(S)), ] <- S
    out
  }
  
  S1 <- fill_groups(rowsum(xc, group = grp, reorder = TRUE), n.nuggets, n.col)
  S2 <- fill_groups(rowsum(xc * xc, group = grp, reorder = TRUE), n.nuggets, n.col)
  
  n.safe <- pmax(n.i, 1L)
  mean.c <- S1 / n.safe   # centred means
  ss <- S2 - n.safe * mean.c^2.  # within-nugget sums of squares
  ss[ss < 0] <- 0     # guard tiny negative round-off
  var.mat <- ss / pmax(n.i - 1L, 1L)
  scale.v <- unname(rowMeans(var.mat))
  scale.v[n.i < 2L] <- NA_real_    
  
  centers <- switch(
    center.method,
    mean     = sweep(mean.c, 2L, g.mean, "+"),
    original = DN.data_mat,
    random   = {
      idx.by.nugget <- split(seq_len(obs.num), grp)
      t(vapply(seq_len(n.nuggets), function(i) {
        ii <- idx.by.nugget[[i]]
        if (length(ii) == 0L) DN.data_mat[i, ]
        else x_mat[ii[sample.int(length(ii), 1L)], ]
      }, numeric(n.col)))
    })
  centers <- unname(as.matrix(centers))
  
  
  if (center.method == "mean" && any(n.i == 0L)) {
    centers[n.i == 0L, ] <- DN.data_mat[n.i == 0L, , drop = FALSE]
  }
  
  message("completed!")
  
  
  ## -----------------------------------------
  ## Create Data nugget information 
  ## -----------------------------------------
  
  
  # Create the data nugget information data frame
  DN.information <- data.frame(seq_len(n.nuggets), centers, 
                               n.i, check.names = FALSE) 

  
  # give column names for the centers and scale
  colnames(DN.information) <- c("Data Nugget",
                                paste0("Center", seq_len(n.col)),
                                "Weight")
  
  multi <- scale.v[n.i > 1L & is.finite(scale.v)]
  floor.scale <- if (length(multi) > 0L) min(multi) / 100 else .Machine$double.eps
  if (!is.finite(floor.scale) || floor.scale <= 0) floor.scale <- .Machine$double.eps
  scale.v[!is.finite(scale.v)] <- floor.scale
  
  DN.information$Scale     <- scale.v
  rownames(DN.information) <- seq_len(n.nuggets)
  
  
  # create output dataframe
  output <- list("Data Nuggets" = DN.information,
                "Data Nugget Assignments" = DN.assignments)
  
  
  # assign the data nugget class to the output
  class(output) <- "datanugget"
  
  # return the data nugget
  return(output)
  
}