Added the site selction functionalities to pecan

docker/depends/pecan_package_dependencies.csv

@@ -281,6 +281,7 @@
 "parallel","*","modules/assim.batch","Imports",FALSE
 "parallel","*","modules/data.atmosphere","Suggests",FALSE
 "parallel","*","modules/data.remote","Imports",FALSE
+"pathviewr","*","modules/assim.sequential","Suggests",FALSE
 "PEcAn.allometry","*","base/all","Suggests",TRUE
 "PEcAn.assim.batch","*","base/all","Depends",TRUE
 "PEcAn.assim.batch","*","modules/rtm","Imports",TRUE

modules/assim.sequential/DESCRIPTION

@@ -48,6 +48,7 @@ Suggests:
     magic (>= 1.5.0),
     methods,
     Metrics,
+    pathviewr,
     PEcAn.benchmark,
     PEcAn.data.land,
     PEcAn.data.remote,

modules/assim.sequential/NAMESPACE

@@ -34,6 +34,7 @@ export(downscale_metrics)
 export(downscale_qsub_main)
 export(dwtmnorm)
 export(ensemble_downscale)
+export(faiss_kmeans)
 export(get_ensemble_weights)
 export(hop_test)
 export(interactive.plotting.sda)
@@ -64,6 +65,7 @@ export(sda.qsub.job.submission)
 export(sda_assemble)
 export(sda_weights_site)
 export(simple.local)
+export(site_selction)
 export(stack_covariates_2_geotiff)
 export(subset_ensemble)
 export(tobit.model)

modules/assim.sequential/R/faiss_kmeans.R

@@ -0,0 +1,61 @@
+#' Perform K-means clustering using the FAISS Python library
+#'
+#' @param data_matrix Numeric matrix: Input dataset where each row is a sample and each column is a feature.
+#' @param k Integer: The desired number of clusters (centroids).
+#' @param niter Integer: Number of K-means training iterations (default: 20).
+#' @param verbose Logical: Whether to display clustering progress and loss details during training (default: TRUE).
+#' @param n_cores Integer: Number of CPU threads to use for execution (default: number of available cores minus 1).
+#'
+#' @return A FAISS K-means object containing cluster centroids and assignment details.
+#' @export
+#'
+#' @examples
+#' \dontrun{
+#'   mat <- matrix(runif(1000 * 10), nrow = 1000, ncol = 10)
+#'   res <- faiss_kmeans(data_matrix = mat, k = 5, niter = 25)
+#' }
+faiss_kmeans <- function(data_matrix, k, niter = 20, verbose = TRUE, n_cores = parallel::detectCores() - 1) {
+  # Make sure we have the faiss package before using it.
+  if (!reticulate::py_module_available("faiss")) {
+    PEcAn.logger::logger.info("Python library faiss is not installed!")
+    return(0)
+  }
+
+  # Import faiss object.
+  faiss <- reticulate::import("faiss")
+  # Set OpenMP threads for FAISS CPU execution
+  faiss$omp_set_num_threads(as.integer(n_cores))
+
+  # Reformat data to meet FAISS requirements
+  data_mat <- as.matrix(data_matrix)
+  mode(data_mat) <- "single"
+  d <- ncol(data_mat)
+
+  # Initialize FAISS K-means object
+  kmeans <- faiss$Kmeans(
+    d       = as.integer(d),
+    k       = as.integer(k),
+    niter   = as.integer(niter),
+    verbose = verbose
+  )
+
+  # Train K-means on matrix
+  kmeans$train(data_mat)
+
+  # Compute cluster assignments.
+  index <- faiss$IndexFlatL2(as.integer(d))
+  index$add(kmeans$centroids)
+
+  # Search nearest centroid for each point (explicit integer for k)
+  res <- index$search(data_mat, 1L)
+
+  # Total within-cluster variance (Total WCSS)
+  total_within_ss <- sum(as.vector(res[[1]]))
+
+  list(
+    centroids = kmeans$centroids,
+    cluster_ids = as.vector(res[[2]]) + 1, # Convert 0-indexed Python to 1-indexed R
+    distances = as.vector(res[[1]]),
+    total_within_ss = total_within_ss
+  )
+}

modules/assim.sequential/R/site_selection.R

@@ -0,0 +1,126 @@
+#' Perform Stratified Site Selection
+#'
+#' Selects sampling sites stratified by land cover class `lc`, scaling sample 
+#' sizes proportionally while guaranteeing a minimum sample size per class. 
+#' Optional normalization is applied prior to clustering.
+#'
+#' @param df Data frame or matrix containing feature columns, a land cover column `lc`, 
+#'   and a weight column `weights`, where each row is a sample and each column is a feature.
+#' @param tot.sample.size Integer. Total target number of sample points to select across all land cover classes.
+#' @param min.sample.size Integer. Minimum number of sample points allocated to each land cover class. (Default = 50).
+#' @param normalize Logical. Whether to apply min-max scaling to numeric columns in `df`. (Default = TRUE).
+#' @param max.cluster.num Integer. Maximum number of clusters to consider if calculating optimal cluster count. (Default = 6).
+#' @param elbow.operation Logical. Whether to automatically calculate the proper cluster number based on the elbow location. (Default = TRUE).
+#' If FALSE, the `max.cluster.num` will be used directly as the cluster size.
+#'
+#' @return A list of cluster allocations and sample indices grouped by land cover class.
+#' @export
+site_selction <- function (df, tot.sample.size, min.sample.size = 50, normalize = TRUE, max.cluster.num = 6, elbow.operation = TRUE) {
+  # grab column names from df.
+  col.names <- colnames(df)
+  # make sure we have both land cover (lc) and weights (weight) layer in the df object.
+  if (!grepl("lc", col.names) & !grepl("weights", col.names)) {
+    PEcAn.logger::logger.info("Please add land cover (lc) and weights (weight) columns to the df!")
+    return(0)
+  }
+  # remove NAs if there is any.
+  if (any(!complete.cases(df))) {
+    non.na.inds <- which(complete.cases(df))
+    df <- df[non.na.inds]
+  } else {
+    non.na.inds <- NULL
+  }
+  # normalize df.
+  if (normalize) {
+    # Define min-max scaling function
+    min_max_scale <- function(x) {
+      (x - min(x, na.rm = TRUE)) / (max(x, na.rm = TRUE) - min(x, na.rm = TRUE))
+    }
+
+    # Apply to all numeric columns in df
+    df <- as.data.frame(lapply(df, function(col) {
+      if (is.numeric(col)) min_max_scale(col) else col
+    }))
+  }
+  # calculate sample size per each land cover class.
+  # land cover operation.
+  LC.inds <- sort(unique(df$lc)) %>% 
+    purrr::map(function(i) {
+      which(df$lc == (i))
+    })
+  proportion <- LC.inds %>% purrr::map(length) %>% unlist
+  proportion <- proportion/sum(proportion)
+  tot.points.num <- smooth.div(tot.sample.size, (tot.sample.size - min.sample.size*length(LC.inds)) * proportion + min.sample.size) - pre.selected.LC.num
+  # loop over land cover classes.
+  clusters <- vector("list", length = length(LC.inds))
+  for (i in seq_along(clusters)) {
+    # prepare datasets for the clustering of the current land cover class.
+    temp.df <- df[LC.inds[[i]], ]
+    temp.weights <- temp.df$weights
+    # if we want to calculate the proper cluster size based on the elbow location of total_within_ss.
+    if (elbow.operation) {
+      total_within_ss <- c()
+      for (j in seq_len(max.cluster.num)) {
+        temp.kmeans <- faiss_kmeans(data_matrix = temp.df, k = j)
+        total_within_ss <- c(total_within_ss, temp.kmeans$total_within_ss)
+      }
+      size <- pathviewr::find_curve_elbow(df)
+    } else {
+      size <- max.cluster.num # default cluster size.
+    }
+    # calculate sample size for each cluster.
+    clust.sample.size <- smooth.div(tot.points.num[i],  rep(tot.points.num[i]/size, size))
+    # conduct kmeans clustering with the selected cluster size.
+    temp.kmeans <- faiss_kmeans(data_matrix = temp.df, k = size)
+    # loop over clusters.
+    temp.cluster.inds <- vector("list", length = size)
+    for (j in seq_len(size)) {
+      temp.inds <- which(temp.kmeans$cluster_ids == j)
+      # weighted sampling.
+      samples <- sample(x = temp.inds, size = clust.sample.size[j], 
+                        replace = F, prob = temp.weights[temp.inds])
+      # below we need to convert to the actual index.
+      # if we observe NAs in the original data frame.
+      if (!is.null(non.na.inds)) {
+        temp.cluster.inds[[j]] <- non.na.inds[LC.inds[[i]][samples]]
+      } else {
+        temp.cluster.inds[[j]] <- LC.inds[[i]][samples]
+      }
+    }
+    # store results.
+    clusters[[i]] <- list(cluster.ids = temp.cluster.inds)
+  }
+  return(clusters)
+}
+
+#' Smoothly Distribute a Total Amount Across Vector Elements
+#'
+#' Distributes an integer total (`tot`) across the elements of a numeric vector (`v`). 
+#' by alternating ceiling and floor operations to prevent rounding bias, while ensuring 
+#' that the returned vector sums exactly to `tot`.
+#'
+#' @param tot Numeric/Integer. The target sum total to distribute.
+#' @param v Numeric vector containing fractional values or proportion targets.
+#'
+#' @return A numeric vector of the same length as `v` containing integer values that sum to `tot`.
+smooth.div <- function(tot, v) {
+  vec <- c()
+  ceil <- TRUE
+  remain <- tot
+
+  for (i in seq_along(v)) {
+    if (i == length(v)) {
+      vec <- c(vec, remain)
+      break
+    }
+    if (ceil) {
+      temp <- ceiling(v[i])
+    } else {
+      temp <- floor(v[i])
+    }
+    ceil <- !ceil
+    remain <- remain - temp
+    vec <- c(vec, temp)
+  }
+  vec
+}