c++

Author: ConnorB

NAMESPACE

@@ -33,6 +33,7 @@ export(ks_meso_station_activity)
 export(ks_meso_stations)
 export(ks_meso_vars)
 export(pCO2_to_xCO2)
+export(rcpp_calc_exceedance_prob)
 export(read_ghcnh)
 export(read_ks_meso)
 export(read_shp)

NEWS.md

@@ -20,4 +20,5 @@
 * `has_units()` has been removed. The `units` package is no longer a dependency (#TBD).
 * `ks_meso_fw13()` can now retrieve Kansas Mesonet fire weather data in FW13 format for one station and date range (#TBD).
 * `ks_meso_most_recent()` can now retrieve the most recently ingested Kansas Mesonet timestamp for each station at a requested interval (#TBD).
+* `rcpp_calc_exceedance_prob()` now provides a C++ implementation of flow exceedance probability calculations (#TBD).
 * `read_ghcnh()` is now exported, validates file inputs, can suppress progress messages with `quiet = TRUE`, and handles files that lack `Station_name` or `Station_ID` columns (#TBD).

R/RcppExports.R

@@ -1,6 +1,10 @@
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
+cpp_calc_exceedance_prob <- function(flow, rm_zero = FALSE) {
+    .Call(`_preMetabolizer_cpp_calc_exceedance_prob`, flow, rm_zero)
+}
+
 flag_z <- function(x, width = 5L, threshold = 3.0, return_z = FALSE) {
     .Call(`_preMetabolizer_flag_z`, x, width, threshold, return_z)
 }

R/calc_exceedance_prob.R

@@ -78,3 +78,49 @@ calc_exceedance_prob <- function(flow, rm.zero = FALSE) {
 
   return(full_exceedance_prob)
 }
+
+#' Calculate Flow Exceedance Probabilities with C++
+#'
+#' Calculates exceedance probabilities using the same Weibull plotting position
+#' method and return shape as [calc_exceedance_prob()], but delegates ranking to
+#' a C++ implementation.
+#'
+#' @inheritParams calc_exceedance_prob
+#'
+#' @return A numeric vector of exceedance probabilities. If `rm.zero = TRUE`,
+#'   the returned vector has the same length as the input with `NA` at positions
+#'   of zero or negative values.
+#'
+#' @seealso [calc_exceedance_prob()]
+#'
+#' @examples
+#' rcpp_calc_exceedance_prob(c(10, 5, 0, 15, 8, NA, 0, 20))
+#'
+#' @export
+rcpp_calc_exceedance_prob <- function(flow, rm.zero = FALSE) {
+  if (!is.numeric(flow)) {
+    stop("'flow' must be a numeric vector.")
+  }
+
+  if (any(is.infinite(flow))) {
+    stop("'flow' contains infinite values. These are not allowed.")
+  }
+
+  if (!is.logical(rm.zero)) {
+    stop("'rm.zero' must be a logical value (TRUE or FALSE).")
+  }
+
+  if (length(flow) == 0) {
+    warning("'flow' is an empty vector. Returning an empty numeric vector.")
+    return(numeric(0))
+  }
+
+  if (rm.zero && length(flow[flow > 0 & !is.na(flow)]) == 0) {
+    warning(
+      "All non-NA values in 'flow' were zero. Returning a vector of NAs with the original length."
+    )
+    return(rep(NA, length(flow)))
+  }
+
+  cpp_calc_exceedance_prob(flow, rm.zero)
+}

_pkgdown.yml

@@ -57,6 +57,7 @@ reference:
       - calc_cv
       - calc_mode
       - calc_exceedance_prob
+      - rcpp_calc_exceedance_prob
       - flag_z
 
   - title: Unit conversion

man/rcpp_calc_exceedance_prob.Rd

@@ -10,6 +10,18 @@ Rcpp::Rostream<true>&  Rcpp::Rcout = Rcpp::Rcpp_cout_get();
 Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
 #endif
 
+// cpp_calc_exceedance_prob
+NumericVector cpp_calc_exceedance_prob(NumericVector flow, bool rm_zero);
+RcppExport SEXP _preMetabolizer_cpp_calc_exceedance_prob(SEXP flowSEXP, SEXP rm_zeroSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< NumericVector >::type flow(flowSEXP);
+    Rcpp::traits::input_parameter< bool >::type rm_zero(rm_zeroSEXP);
+    rcpp_result_gen = Rcpp::wrap(cpp_calc_exceedance_prob(flow, rm_zero));
+    return rcpp_result_gen;
+END_RCPP
+}
 // flag_z
 SEXP flag_z(NumericVector x, int width, double threshold, bool return_z);
 RcppExport SEXP _preMetabolizer_flag_z(SEXP xSEXP, SEXP widthSEXP, SEXP thresholdSEXP, SEXP return_zSEXP) {
@@ -26,6 +38,7 @@ END_RCPP
 }
 
 static const R_CallMethodDef CallEntries[] = {
+    {"_preMetabolizer_cpp_calc_exceedance_prob", (DL_FUNC) &_preMetabolizer_cpp_calc_exceedance_prob, 2},
     {"_preMetabolizer_flag_z", (DL_FUNC) &_preMetabolizer_flag_z, 4},
     {NULL, NULL, 0}
 };

src/RcppExports.cpp

@@ -0,0 +1,60 @@
+#include <Rcpp.h>
+#include <algorithm>
+#include <vector>
+using namespace Rcpp;
+
+struct FlowValue {
+  double value;
+  int output_index;
+};
+
+// [[Rcpp::export]]
+NumericVector cpp_calc_exceedance_prob(NumericVector flow, bool rm_zero = false) {
+  int n_flow = flow.size();
+  NumericVector result(n_flow, NA_REAL);
+  std::vector<FlowValue> values;
+  values.reserve(n_flow);
+
+  for (int i = 0; i < n_flow; ++i) {
+    double value = flow[i];
+    if (NumericVector::is_na(value)) {
+      continue;
+    }
+    if (rm_zero && value <= 0) {
+      continue;
+    }
+    values.push_back({value, i});
+  }
+
+  int n = values.size();
+  if (n == 0) {
+    return result;
+  }
+
+  std::sort(
+    values.begin(),
+    values.end(),
+    [](const FlowValue& a, const FlowValue& b) {
+      return a.value > b.value;
+    }
+  );
+
+  int i = 0;
+  while (i < n) {
+    int j = i + 1;
+    while (j < n && values[j].value == values[i].value) {
+      ++j;
+    }
+
+    double rank = (static_cast<double>(i + 1) + static_cast<double>(j)) / 2.0;
+    double exceedance_probability = rank / (static_cast<double>(n) + 1.0);
+
+    for (int k = i; k < j; ++k) {
+      result[values[k].output_index] = exceedance_probability;
+    }
+
+    i = j;
+  }
+
+  return result;
+}

src/calc_exceedance_prob.cpp

@@ -27,3 +27,29 @@ test_that("calc_exceedance_prob errors on non-numeric input", {
 test_that("calc_exceedance_prob errors on Inf values", {
   expect_snapshot(error = TRUE, calc_exceedance_prob(c(1, Inf, 3)))
 })
+
+test_that("rcpp_calc_exceedance_prob matches R implementation", {
+  flow <- c(10, 5, 0, 15, 8, NA, 0, 20, 10)
+
+  expect_equal(
+    rcpp_calc_exceedance_prob(flow),
+    calc_exceedance_prob(flow)
+  )
+  expect_equal(
+    rcpp_calc_exceedance_prob(flow, rm.zero = TRUE),
+    calc_exceedance_prob(flow, rm.zero = TRUE)
+  )
+})
+
+test_that("rcpp_calc_exceedance_prob handles ties with average ranks", {
+  result <- rcpp_calc_exceedance_prob(c(10, 20, 20, 30))
+
+  expect_equal(result, c(0.8, 0.5, 0.5, 0.2))
+})
+
+test_that("rcpp_calc_exceedance_prob preserves NA positions", {
+  result <- rcpp_calc_exceedance_prob(c(10, NA, 30))
+
+  expect_equal(result[2], NA_real_)
+  expect_equal(result[1], 2 / 3)
+})