R/refine.DN.R

In datanugget: Create, Optimize, and Refine Data Nuggets

# This function refines the datanuggets created by create.DN() function.
# 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.
# DN: a datanugget object.
# EV.tol: percentile at which the eigen value threshold is taken. Default is 0.9.
# max.splits: maximum number of iterations for splitting. Default is 5.
# min.nugget.size: minimum number of observations a data nugget may contain as
#                  a result of a split. Default is 2.
# shape.split: logical; whether to run the shape splitting stage.
# min.shape.size: minimum weight for a data nugget to be shape split. Default is 10.
# delta: a data nugget is elongated when its first two eigenvalues satisfy
#        EV1/EV2 > delta. Default is 2.
# nstart: number of random starts used by kmeans when splitting. Default is 25.
# 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.


refine.DN = function(x,
                     DN,
                     EV.tol = 0.9,
                     max.splits = 5,
                     min.nugget.size = 2,
                     shape.split = FALSE,
                     min.shape.size = 10, 
                     delta = 2,
                     nstart = 25,
                     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 DN
  if (!inherits(DN,"datanugget")){
    stop('DN must be of class "datanugget"')
  }
  
  
  
  # check EV.tol
  if (!is.numeric(EV.tol) || EV.tol >= 1 || EV.tol <= 0){
    stop('EV.tol must be within (0, 1)')
  }
  
  
  
  # check max.splits
  if (!(class(max.splits) %in% c("numeric","integer")) || max.splits < 0){
    stop('max.splits must be non-negative')
  }
  
  
  
  # check min.nugget.size
  if (!is.numeric(min.nugget.size) || min.nugget.size < 2){
    stop('min.nugget.size must be at least 2')
  }
  
  
  
  # check min.shape.size
  if (!is.numeric(min.shape.size) || min.shape.size < 10){
    stop("min.shape.size must be at least 10")
  }
  
  
  
  # check delta
  if (!is.numeric(delta) || delta <= 1) {
    stop("delta must be greater than 1")
  }
  
  
  # check nstart 
  if (!is.numeric(nstart)){
    stop('nstart must be of class "numeric"')
  }
  
  
  
  # 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")
  
  
  # covert max.splits and min.nugget.size to integer
  max.splits <- floor(max.splits)
  min.nugget.size <- floor(min.nugget.size)
  
  
  # 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)
  
  
  # Data nuggets information
  DN.information <- DN[["Data Nuggets"]]
  DN.assignments <- as.integer(DN[["Data Nugget Assignments"]])
  if (length(DN.assignments) != obs.num)
    stop("length of the data nugget assignments must equal nrow(x)")
  DN.ids <- as.integer(DN.information[, "Data Nugget"])
  DN.centers <- as.matrix(DN.information[, 1 + seq_len(n.col), drop = FALSE])
  DN.weight <- as.integer(DN.information[, "Weight"])
  DN.scale <- as.numeric(DN.information[, "Scale"])
  
  
  # next available data nugget id
  next.id <- max(DN.ids) + 1L
  
  
  # for storing two largest eigenvalues of each data nugget
  DN.EV1 <- rep(NA_real_, length(DN.ids))
  DN.EV2 <- rep(NA_real_, length(DN.ids))
  
  
  # 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 make.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)
    
    
    # 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 summarize one data nugget from the observation it contains
  get_nugget_features <- function(assign.obs, n.col){
    
    # no of assigned observation in the nugget
    m <- nrow(assign.obs)
    
    # if empty return a row of all zeroes
    if (m == 0L){
      return(list(center = rep(NA_real_, n.col), weight = 0L,
                  scale = 0, EV1 = 0, EV2 = 0))}
    
    # center 
    ctr <- colMeans(assign.obs)
    
  
    # a single observation has no within-data-nugget variability
    if (m == 1L){
      return(list(center = ctr, weight = 1L, scale = 0, EV1 = 0, EV2 = 0))}
    
    
    # scale
    S <- stats::cov(assign.obs)
    sc <- sum(diag(S)) / n.col
    
    # eigenvalues
    ev <- if (n.col == 1L) as.numeric(S) else 
      eigen(S, symmetric = TRUE, only.values = TRUE)$values
    
    
    return(list(center = ctr, weight = m, scale = sc,
                EV1 = ev[1],
                EV2 = if (length(ev) > 1L) ev[2] else ev[1]))
  }
  
  environment(get_nugget_features) <- clean.env
  
  
  # Function to split one data nugget into two using K-Means algorithm.
  split_nugget <- function(assign.obs, min.nugget.size, nstart, project, seed){
    
    
    # check if the current data nugget can be split
    if (nrow(assign.obs) < 2 * min.nugget.size) return (NULL)
    
    # split on the projection onto the first eigen vector when shape splitting
    if (project){
      eg <- eigen(stats::cov(assign.obs), symmetric = TRUE)
      to.split <- as.numeric(assign.obs %*% eg$vectors[, 1, drop = FALSE])
    } else to.split <- assign.obs
    
    set.seed(seed)
    new.splitted.nuggets <- tryCatch(stats::kmeans(x = to.split, 
                                                   centers = 2, 
                                                   iter.max = 1000, 
                                                   nstart = nstart)$cluster, 
                                     error = function(e) NULL)
    
    if (is.null(new.splitted.nuggets)) return(NULL)
    
    if (min(tabulate(new.splitted.nuggets, nbins = 2)) < min.nugget.size) return(NULL)
    
    return(new.splitted.nuggets)
    
  }
  
  environment(split_nugget) <- clean.env
  
  
  # Function to run the splits over a list of candidate data nuggets
  run_splitting <- function(sub.list, FUN){
    
    if(length(sub.list) == 0) return(list())
    
    # 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 <- utils::txtProgressBar(min = 0, max = length(sub.list))
        
        # update the progress bar
        progress <- function(n){utils::setTxtProgressBar(pb, n)}
        opts <- list(progress = progress)
        
        # close the progress bar upon exit
        on.exit(close(pb), add = TRUE)
      }
      
      
      res <- foreach::foreach(i = seq_along(sub.list), 
                              .inorder = TRUE, 
                              .options.snow = opts
      )%dopar% {
        
        FUN(sub.list[[i]])
      
      }
    }else{
      
      if (make.pbs){
        
        pb <- utils::txtProgressBar(min = 0, max = length(sub.list))
        on.exit(close(pb), add = TRUE)
      }
      
      res <- vector("list", length(sub.list))
      
      for(i in seq_along(sub.list)){
        
        res[[i]] <- FUN(sub.list[[i]])
        
        if (make.pbs) utils::setTxtProgressBar(pb, i)
      }
    }
    
  return(res)  
  }
  
  
  # Function to list the observations belonging to each data nugget
  get_members <- function(){
    split(seq_len(obs.num), factor(DN.assignments, levels = DN.ids))
  }
  
  
  # Function to recalculate the information of data nuggets once a split occurs
  refresh_nuggets_info <- function(which.pos, members){
    
    for (k in which.pos){
      
      tmp.stats <- get_nugget_features(x_mat[members[[k]], , drop = FALSE], n.col)
      
      DN.centers[k, ] <<- tmp.stats$center
      DN.weight[k] <<- tmp.stats$weight
      DN.scale[k] <<- tmp.stats$scale
      DN.EV1[k] <<- tmp.stats$EV1
      DN.EV2[k] <<- tmp.stats$EV2
      
    }
  }
  
  
  # Function to assign observations to new nuggets
  apply_splits <- function(splits){
    
    touched <- integer(0)
    
    for (s in splits){
      
      # reassign observations to the new data nugget
      DN.assignments[s$rows[s$cluster == 2]] <<- next.id
      
      DN.ids <<- c(DN.ids, next.id)
      DN.centers <<- rbind(DN.centers, rep(NA_real_, n.col))
      DN.weight <<- c(DN.weight, NA_real_)
      DN.scale <<- c(DN.scale, NA_real_)
      DN.EV1 <<- c(DN.EV1, NA_real_)
      DN.EV2 <<- c(DN.EV2, NA_real_)
      
      touched <- c(touched, s$pos, length(DN.ids))
      next.id <<- next.id + 1L
    }
    return(touched)
  }
  
  
  # calculate the statistics of the existing nuggets before the split
  tmp.members <- get_members()
  refresh_nuggets_info(seq_along(DN.ids), tmp.members)
  
  
  ## -------------------------------------------------------
  ## Splitting data nuggets according to their max eigenvalue
  ## -------------------------------------------------------
  
  if (max.splits > 0L){
    
    # positive first eigen values
    tmp.EVs <- DN.EV1[is.finite(DN.EV1) & DN.EV1 > 0]
    
    if (length(tmp.EVs) == 0){
      message("No data nugget has a positive first eigen value")
    } else{
      
      # threshold on the maximum eigen values
      EV.threshold <- stats::quantile(tmp.EVs, probs = EV.tol)
      
      message("Splitting data nuggets according to max EV...")
      
      for (split.attempts in seq_len(max.splits)){
        
        # candidate data nuggets
        large.max.EVs <- which(is.finite(DN.EV1) & 
                                 DN.EV1 > EV.threshold & 
                                 DN.weight >= 2 * min.nugget.size)
        
        if (length(large.max.EVs) == 0L) break
        
        message(paste("split #", split.attempts, sep = ""))
        
        tmp.members <- get_members()
        
        # split the candidate data nuggets
        tmp.splits <- run_splitting(
          lapply(large.max.EVs, function(k)
            x_mat[tmp.members[[k]], , drop = FALSE]), 
          function(assign.obs)
            split_nugget(assign.obs, min.nugget.size, nstart, project = FALSE, seed = seed))
        
        
        # keep the nuggets that were split successfully
        tmp.keep <- which(!vapply(tmp.splits, is.null, logical(1)))
        
        if(length(tmp.keep) == 0L) break
        
        tmp.splits <- lapply(tmp.keep, function(j)
          list(pos = large.max.EVs[j], 
               rows = tmp.members[[large.max.EVs[j]]], 
               cluster = tmp.splits[[j]]))
        
        tmp.touched <- apply_splits(tmp.splits)
        
        tmp.members <- get_members()
        refresh_nuggets_info(tmp.touched, tmp.members)
        
      }
      message("complete!")
    }
  }
  
  
  
  
  ## -------------------------------------------------------
  ## Splitting data nuggets according to their shape
  ## -------------------------------------------------------
  
  
  if (shape.split & n.col >= 2L & max.splits > 0L){
    
    message("Splitting data nuggets according to shape...")
    
    for (split.attempts in seq_len(max.splits)){
      
      tmp.finite <- is.finite(DN.scale) & DN.weight > 1L
      
      s.threshold <- if (any(tmp.finite)){
        # threshold for elongated nuggets
        stats::median(DN.scale[tmp.finite]) + 2 * stats::mad(DN.scale[tmp.finite])
      } else Inf
      
      
      # ratio of the first two eigen values
      tmp.ratio <- ifelse(is.finite(DN.EV2) & DN.EV2 > 0 , DN.EV1 / DN.EV2, Inf)
      
      # candidate nuggets
      large.shapes <- which(DN.weight >= max(min.shape.size, 2 * min.nugget.size) & 
                              is.finite(DN.EV1) & DN.EV1 > 0 & 
                              (tmp.ratio > delta | DN.EV1 > s.threshold)) 
      
      if (length(large.shapes) == 0L) break
      
      message(paste("split #", split.attempts, sep = ""))
        
      tmp.members <- get_members()
      
      # split the candidate data nuggets projected on the first eigen vector
      tmp.splits <- run_splitting(
        lapply(large.shapes, function(k)
          x_mat[tmp.members[[k]], , drop = FALSE]), 
        function(assign.obs){
          
          new.split <- split_nugget(assign.obs, min.nugget.size, nstart, project = TRUE, seed = seed)
          
          if (is.null(new.split)) return(NULL)
        
        
          # a shape split is kept if only it makes the two nuggets sufficiently rounder
          tmp.ratios <- vapply(1:2, function(g){
            
            tmp.obs <- assign.obs[new.split == g, , drop = FALSE]
            
            if (nrow(tmp.obs) < 2L) return(NA_real_)
            
            tmp.ev <- eigen(stats::cov(tmp.obs), symmetric = TRUE, 
                            only.values = TRUE)$values
            
            if(length(tmp.ev) < 2L || tmp.ev[2] <= 0) return(Inf)
            
            return(tmp.ev[1]/tmp.ev[2])
          }, numeric(1))
          
          
          if(anyNA(tmp.ratios) || mean(tmp.ratios) >= delta) return(NULL)
          
          return(new.split)
        })
      
      # keep the nuggets that were split successfully
      tmp.keep <- which(!vapply(tmp.splits, is.null, logical(1)))
      
      if (length(tmp.keep) == 0L) break
      
      tmp.splits <- lapply(tmp.keep, function(j)
        list(pos = large.shapes[j], 
             rows = tmp.members[[large.shapes[j]]], 
             cluster = tmp.splits[[j]]))
      
      tmp.touched <- apply_splits(tmp.splits)
      
      tmp.members <- get_members()
      refresh_nuggets_info(tmp.touched, tmp.members)
    }
    
    message("complete")
  }
  
  
  
  ## -----------------------------------------
  ## Final data nugget information
  ## -----------------------------------------
  
  
  
  message("Renumbering the data nuggets...")
  
  # drop any nugget with no observation
  keep <- which(DN.weight > 0L)
  DN.ids <- DN.ids[keep]
  DN.centers <- DN.centers[keep, , drop = FALSE]
  DN.weight <- DN.weight[keep]
  DN.scale <- DN.scale[keep]
  
  
  # renumber the data nuggets 
  renum <- integer(max(DN.ids))
  renum[DN.ids] <- seq_along(DN.ids)
  DN.assignments <- renum[DN.assignments]
  
  # a data nugget of one observation has no within-data-nugget variability
  DN.scale[DN.weight < 2L] <- 0
  
  # Create the data nugget information data frame
  DN.information <- data.frame(seq_along(DN.ids), DN.centers,
                               DN.weight, DN.scale, check.names = FALSE)
  
  colnames(DN.information) <- c("Data Nugget",
                                paste0("Center", seq_len(n.col)),
                                "Weight", "Scale")
  
  rownames(DN.information) <- seq_along(DN.ids)
  
  
  stopifnot(sum(DN.information[, "Weight"]) == obs.num,
            !anyNA(DN.assignments),
            all(DN.information[, "Weight"] > 0),
            all(sort(unique(DN.assignments)) == seq_len(nrow(DN.information))))
  
  message("completed!")
  
  
  # 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)
}