doc update, error check, etc

Author: rvalavi

DESCRIPTION

@@ -2,7 +2,7 @@ Package: ClassicHCAS
 Type: Package
 Title: The Classic Habitat Condition Assessment System (HCAS).
 Version: 1.0.0
-Date: 2025-10-08
+Date: 2025-10-15
 Author: Roozbeh Valavi, Chris Ware, Eric Lehmann, Kristen Williams, Mike Birchall, Simon Collings, Simon Ferrier, and Tom Harwood
 Maintainer: <roozbeh.valavi@csiro.au>
 Description: The HCAS evaluates habitat condition by comparing observed and 

R/benchmark.R

@@ -1,4 +1,4 @@
-#' Condition benchmarking of target points
+#' Condition benchmarking of target sites
 #'
 #' The HCAS (Habitat Condition Assessment System) benchmarking function evaluates habitat
 #' condition by comparing observed and predicted remote sensing (RS) variables. It integrates
@@ -64,9 +64,9 @@
 #' @param confidence Numeric. The confidence value for LDC methods. See details below..
 #' @param lambda Numeric. The lambda param for LDC Cauchy weighting...
 #' @param exclude_slef Logical. To exclude a benchmark point from assessing itself.
-#' @param drop_features Integer vector. Completely remove the RS variable from the benchmarking process. For
-#' consistency, it is recommended to exclude the same variables used in the histogram step; unless
-#' you have a specific reason not to.
+#' @param drop_features Integer vector of remote sensing variable IDs. Completely remove the RS variable
+#' from the benchmarking process. For consistency, it is recommended to exclude the same variables used
+#' in the histogram step; unless you have a specific reason not to.
 #' @param make_su Logical. To make the uncertainty map or not.
 #' @param ... Additional arguments for writing raster outputs e.g. \code{filename},
 #' \code{overwrite}, and \code{wopt} from terra \code{\link[terra]{predict}}.
@@ -146,6 +146,9 @@ benchmark <- function(
     # drop features from calculation if requested
     if (length(drop_features)) {
         n_vars <- (ncol(samples) - 2) / 2
+        if (any(drop_features > n_vars)) {
+            stop("The feature IDs in drop_features must not be larger than the number of remote sensing variables.")
+        }
         exclude_var <- c(drop_features + 2, drop_features + 2 + n_vars)
     }
 

R/histogram.R

@@ -47,9 +47,9 @@
 #' @param bin_num Integer. Specifies the number of bins for the histogram. It is generally recommended
 #' to use the default value of 650. Adjusting \code{bin_width} is often more effective than changing
 #' \code{bin_num}.
-#' @param drop_features Integer vector. Completely remove the RS variable from the histogram generation. For
-#' consistency, it is recommended to exclude the same variables later in the benchmarking step; unless
-#' you have a specific reason not to.
+#' @param drop_features Integer vector of remote sensing variable IDs. Completely remove the RS variable
+#' from the histogram generation. For consistency, it is recommended to exclude the same variables later
+#' in the benchmarking step; unless you have a specific reason not to.
 #' @param num_threads Integer. Specifies the number of CPU threads to be used for processing. A value
 #' below 1 indicates that all available threads will be utilized (default). Refer to the details section for
 #' more information.
@@ -105,10 +105,10 @@ histogram <- function(
 
 
     # number of observed RS vars
-    num_layers <- (ncol(data_vals) - 2) / 2
+    n_vars <- (ncol(data_vals) - 2) / 2
     # id of obs and mod for saving
-    mod_layers <- seq_len(num_layers) + 2
-    obs_layers <- mod_layers + num_layers
+    mod_layers <- seq_len(n_vars) + 2
+    obs_layers <- mod_layers + n_vars
 
     # get the correct columns for the C++ code
     samples_xy <- data_vals[, 1:2]
@@ -133,6 +133,9 @@ histogram <- function(
 
     # drop features from calculation if requested
     if (length(drop_features)) {
+        if (any(drop_features > n_vars)) {
+            stop("The feature IDs in drop_features must not be larger than the number of remote sensing variables.")
+        }
         observed[, drop_features] <- 0
         modelled[, drop_features] <- 0
     }

man/benchmark.Rd

@@ -84,7 +84,11 @@ Rcpp::NumericMatrix bench_cpp(
     // Cast the values once; cleaner code
     const float32_t xypenalty = static_cast<float32_t>(xy_penalty);
     std::vector<float32_t> xystats(xy_stats.begin(), xy_stats.end());
-
+    // Check the size of stats so indexing doesn't panic!
+    if (xystats.size() != 4) {
+        Rcpp::stop("xy_stats must contain exactly 4 values!");
+    }
+    
     // Normalise XY and apply weight to it
     samples.col(0) = ((samples.col(0).array() - xystats[0]) / xystats[2]) * xypenalty;
     samples.col(1) = ((samples.col(1).array() - xystats[1]) / xystats[3]) * xypenalty;
@@ -189,17 +193,17 @@ Rcpp::NumericMatrix bench_cpp(
     if (make_su)
     {
         // output matrix - both condition and SU values
-        for (const auto& cval : condition_vect) {
-            out_mat(i, 0) = cval.hc;
-            out_mat(i, 1) = std::log(cval.su);
+        for (const auto& v : condition_vect) {
+            out_mat(i, 0) = v.hc;
+            out_mat(i, 1) = std::log(v.su);
             i++;
         }
     }
     else
     {
         // output matrix - only condition values
-        for (const auto& cval : condition_vect) {
-            out_mat(i, 0) = cval.hc;
+        for (const auto& v : condition_vect) {
+            out_mat(i, 0) = v.hc;
             i++;
         }
     }

man/calibrate.Rd

@@ -30,12 +30,12 @@ std::vector<int> combined_Search(
     int64_t cos_scale,
     bool is_geo
 ) {
-    const int size = sample_x.size();
     // Use pair<ENV_Distance, Original_Index>
     std::vector<std::pair<float, int>> dist_idx;
     // Reserve 5000; it's a heuristic, but better than nothing
     dist_idx.reserve(5000);
 
+    const int size = sample_x.size();
     // 1. Perform radius check, and calculate ENV distance for KNN later
     for (int i = 0; i < size; ++i) {
         const int64_t dx = sample_x[i] - query_x;
@@ -108,7 +108,7 @@ inline Condition get_Condition(
         double p_dist = pred_dists[i];
 
         double weight = 1.0;
-        if (p_dist > 0) {
+        if (p_dist > 0.0) {
             weight = 1.0 / (PI_SQ * p_dist * lambda * (1.0 + (p_dist * p_dist) / lambda_sq));
         }
 

man/histogram.Rd

@@ -55,12 +55,10 @@ Rcpp::IntegerMatrix histo_cpp(
     RowMajorMatrix<double> xy = get_XY(xy_vals);
     // define the output matrix;
     RowMajorMatrix<uint64_t> matrix(bin_num, bin_num);
-    // Matrix MUST be initialised with zero
+    // Matrix MUST be initialised with zeros
     matrix.setZero();
 
     const int nr = rs.rows();
-
-
     double scale;
     int64_t squared_dist;
     const double radius_m = radius_km * 1000.0;