Monthly moving avg (#363)

* added a monthly moving average function along with unit tests

* cleaned up comments in monthly moving avg cpp code

Author: Paul-Carvalho

NAMESPACE

@@ -12,6 +12,7 @@ export(assignment_column)
 export(bbox)
 export(bin_var)
 export(bycatch)
+export(calc_monthly_moving_avg)
 export(catch_lm)
 export(category_cols)
 export(centroid_to_fsdb)

R/RcppExports.R

@@ -5,6 +5,10 @@ calculate_moving_avg <- function(unique_dates, unique_groups, obs_dates, obs_gro
     .Call(`_FishSET_calculate_moving_avg`, unique_dates, unique_groups, obs_dates, obs_groups, obs_values, window_size, lag, year_lag, temporal, weighted)
 }
 
+calculate_monthly_avg <- function(unique_years, unique_months, unique_groups, obs_years, obs_months, obs_groups, obs_values, window_size, month_lag, year_lag = 0L, weighted = FALSE) {
+    .Call(`_FishSET_calculate_monthly_avg`, unique_years, unique_months, unique_groups, obs_years, obs_months, obs_groups, obs_values, window_size, month_lag, year_lag, weighted)
+}
+
 shift_sort_xcpp <- function(x, ch, y, distance, alts, ab) {
     .Call(`_FishSET_shift_sort_xcpp`, x, ch, y, distance, alts, ab)
 }

R/calc_monthly_moving_avg.R

@@ -0,0 +1,99 @@
+#' Calculate monthly moving window averages
+#'
+#' @param df Primary data frame or data.table containing the observations.
+#' @param name String, the name of the new column to be added to \code{df} 
+#'   containing the calculated moving averages. Defaults to \code{"moving_avg"}.
+#' @param year_col String, the name of the column containing the year (integer).
+#' @param month_col String, the name of the column containing the month (integer).
+#' @param group_cols Character vector, variable(s) from \code{df} that define how to 
+#'   group the data for averaging (e.g. \code{c("VESSEL_ID", "ZONE")}). The moving 
+#'   average is calculated independently for each group defined here.
+#' @param value_col String, the name of the column containing the numeric values 
+#'   to be averaged.
+#' @param window_size Numeric, the size of the moving window in months. Defaults to 3.
+#' @param month_lag Numeric, the number of months to lag the window. Defaults to 1.
+#'   For example, a lag of 1 starts the average from the previous month.
+#' @param year_lag Numeric, the number of years to lag the reference date. Defaults to 0.
+#'   If set to 1, the function looks up the date from exactly one year prior before 
+#'   applying the \code{month_lag} and \code{window_size}.
+#' @param fill_empty_expectation Numeric or \code{NA}. Value used to fill \code{NA} 
+#'   values in the resulting moving average column. Defaults to \code{NA}.
+#' 
+#' @importFrom data.table as.data.table := .SD setnames setorder
+#' 
+#' @returns Returns the original \code{df} as a data.table with the additional 
+#'   column specified by \code{name}.
+#' @export
+
+calc_monthly_moving_avg <- function(df,
+                                    name,
+                                    year_col, 
+                                    month_col, 
+                                    group_cols,
+                                    value_col,
+                                    window_size = 3,
+                                    month_lag = 1,
+                                    year_lag = 0,
+                                    fill_empty_expectation = NA) {
+  
+  dt <- as.data.table(df)
+  
+  # Create unified group ID
+  # Use a temp name to avoid clashes
+  temp_id_col <- "TEMP_UNIFIED_ID_XYZ"
+  dt[, (temp_id_col) := do.call(paste, c(.SD, sep = "_")), .SDcols = group_cols]
+  
+  # Prepare unique vectors
+  unique_times <- unique(dt[, .(get(year_col), get(month_col))])
+  setnames(unique_times, c("yr", "mth"))
+  setorder(unique_times, yr, mth)
+  
+  unique_groups <- unique(dt[[temp_id_col]])
+  
+  # Call C++ function
+  res_mat <- calculate_monthly_avg(
+    unique_years = unique_times$yr,
+    unique_months = unique_times$mth,
+    unique_groups = unique_groups,
+    
+    obs_years = dt[[year_col]],
+    obs_months = dt[[month_col]],
+    obs_groups = dt[[temp_id_col]],
+    obs_values = dt[[value_col]],
+    
+    window_size = window_size,
+    month_lag = month_lag,
+    year_lag = year_lag,
+    weighted = FALSE
+  )
+  
+  # Convert to data table
+  n_times <- nrow(unique_times)
+  n_groups <- length(unique_groups)
+  
+  lookup_table <- data.table(
+    yr = rep(unique_times$yr, times = n_groups),
+    mth = rep(unique_times$mth, times = n_groups),
+    group_id = rep(unique_groups, each = n_times),
+    val_placeholder = as.vector(res_mat) # Flattens matrix column by column
+  )
+  
+  # Rename columns to match the input df and the desired target column name
+  setnames(lookup_table, 
+           old = c("yr", "mth", "group_id", "val_placeholder"), 
+           new = c(year_col, month_col, temp_id_col, name))
+  
+  dt <- merge(dt, lookup_table, by = c(year_col, month_col, temp_id_col), all.x = TRUE)
+  
+  # Handle missing expectations
+  if (!is.na(fill_empty_expectation)) {
+    # Replace NAs in the NEW column with the provided value
+    # We use get() to refer to the column by its string name
+    dt[is.na(get(name)), (name) := fill_empty_expectation]
+  }
+  
+  # Cleanup
+  dt[, (temp_id_col) := NULL]
+  
+  return(dt)
+}

R/globals.R

@@ -26,5 +26,7 @@ utils::globalVariables(c(
   'model_name',
   'policy_name',
   'rename',
-  'V1'
+  'V1',
+  'yr',
+  'mth'
 ))

man/calc_monthly_moving_avg.Rd

@@ -30,6 +30,27 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// calculate_monthly_avg
+NumericMatrix calculate_monthly_avg(IntegerVector unique_years, IntegerVector unique_months, CharacterVector unique_groups, IntegerVector obs_years, IntegerVector obs_months, CharacterVector obs_groups, NumericVector obs_values, int window_size, int month_lag, int year_lag, bool weighted);
+RcppExport SEXP _FishSET_calculate_monthly_avg(SEXP unique_yearsSEXP, SEXP unique_monthsSEXP, SEXP unique_groupsSEXP, SEXP obs_yearsSEXP, SEXP obs_monthsSEXP, SEXP obs_groupsSEXP, SEXP obs_valuesSEXP, SEXP window_sizeSEXP, SEXP month_lagSEXP, SEXP year_lagSEXP, SEXP weightedSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< IntegerVector >::type unique_years(unique_yearsSEXP);
+    Rcpp::traits::input_parameter< IntegerVector >::type unique_months(unique_monthsSEXP);
+    Rcpp::traits::input_parameter< CharacterVector >::type unique_groups(unique_groupsSEXP);
+    Rcpp::traits::input_parameter< IntegerVector >::type obs_years(obs_yearsSEXP);
+    Rcpp::traits::input_parameter< IntegerVector >::type obs_months(obs_monthsSEXP);
+    Rcpp::traits::input_parameter< CharacterVector >::type obs_groups(obs_groupsSEXP);
+    Rcpp::traits::input_parameter< NumericVector >::type obs_values(obs_valuesSEXP);
+    Rcpp::traits::input_parameter< int >::type window_size(window_sizeSEXP);
+    Rcpp::traits::input_parameter< int >::type month_lag(month_lagSEXP);
+    Rcpp::traits::input_parameter< int >::type year_lag(year_lagSEXP);
+    Rcpp::traits::input_parameter< bool >::type weighted(weightedSEXP);
+    rcpp_result_gen = Rcpp::wrap(calculate_monthly_avg(unique_years, unique_months, unique_groups, obs_years, obs_months, obs_groups, obs_values, window_size, month_lag, year_lag, weighted));
+    return rcpp_result_gen;
+END_RCPP
+}
 // shift_sort_xcpp
 NumericMatrix shift_sort_xcpp(NumericMatrix x, NumericMatrix ch, NumericVector y, NumericMatrix distance, int alts, int ab);
 RcppExport SEXP _FishSET_shift_sort_xcpp(SEXP xSEXP, SEXP chSEXP, SEXP ySEXP, SEXP distanceSEXP, SEXP altsSEXP, SEXP abSEXP) {
@@ -49,6 +70,7 @@ END_RCPP
 
 static const R_CallMethodDef CallEntries[] = {
     {"_FishSET_calculate_moving_avg", (DL_FUNC) &_FishSET_calculate_moving_avg, 10},
+    {"_FishSET_calculate_monthly_avg", (DL_FUNC) &_FishSET_calculate_monthly_avg, 11},
     {"_FishSET_shift_sort_xcpp", (DL_FUNC) &_FishSET_shift_sort_xcpp, 6},
     {NULL, NULL, 0}
 };

src/RcppExports.cpp

@@ -0,0 +1,133 @@
+#include <Rcpp.h>
+#include <map>
+#include <vector>
+#include <string>
+#include <algorithm>
+
+using namespace Rcpp;
+
+// Helper to convert Year/Month to an absolute month index (0-based from Year 0)
+// e.g., 2023, 1 -> 24276
+int get_abs_month(int year, int month) {
+  return year * 12 + (month - 1);
+}
+
+// [[Rcpp::export]]
+NumericMatrix calculate_monthly_avg(IntegerVector unique_years,
+                                    IntegerVector unique_months,
+                                    CharacterVector unique_groups,
+                                    IntegerVector obs_years,
+                                    IntegerVector obs_months,
+                                    CharacterVector obs_groups,
+                                    NumericVector obs_values,
+                                    int window_size,
+                                    int month_lag,
+                                    int year_lag = 0,
+                                    bool weighted = false) {
+  
+  int n_times = unique_years.size();
+  int n_groups = unique_groups.size();
+  int n_obs = obs_years.size();
+  
+  // --- Map Setup ---
+  // Map "Absolute Month" -> Row Index
+  std::map<int, int> time_to_row_map;
+  for(int i = 0; i < n_times; ++i) {
+    int abs_m = get_abs_month(unique_years[i], unique_months[i]);
+    time_to_row_map[abs_m] = i;
+  }
+  
+  // Map "Group Name" -> Column Index
+  std::map<String, int> group_to_col_map;
+  for(int i = 0; i < n_groups; ++i) {
+    group_to_col_map[unique_groups[i]] = i;
+  }
+  
+  // --- Data Aggregation (Handle multiple obs per month/group) ---
+  // Initialize matrices
+  NumericMatrix values_matrix(n_times, n_groups);
+  std::fill(values_matrix.begin(), values_matrix.end(), 0.0);
+  
+  IntegerMatrix count_matrix(n_times, n_groups);
+  std::fill(count_matrix.begin(), count_matrix.end(), 0);
+  
+  for(int i = 0; i < n_obs; ++i) {
+    int abs_m = get_abs_month(obs_years[i], obs_months[i]);
+    
+    // Safety check: ensure this time/group exists in our definition vectors
+    if (time_to_row_map.find(abs_m) != time_to_row_map.end() && 
+        group_to_col_map.find(obs_groups[i]) != group_to_col_map.end()) {
+      
+      int row = time_to_row_map[abs_m];
+      int col = group_to_col_map[obs_groups[i]];
+      
+      if (!NumericVector::is_na(obs_values[i])) {
+        values_matrix(row, col) += obs_values[i];
+        count_matrix(row, col)++;
+      }
+    }
+  }
+  
+  // Calculate raw monthly means (normalize the aggregated sums)
+  for (int r = 0; r < n_times; ++r) {
+    for (int c = 0; c < n_groups; ++c) {
+      if (count_matrix(r, c) > 0) {
+        values_matrix(r, c) = values_matrix(r, c) / count_matrix(r, c);
+      } else {
+        values_matrix(r, c) = NA_REAL;
+      }
+    }
+  }
+  
+  // --- Sliding Window Calculation ---
+  
+  NumericMatrix result_matrix(n_times, n_groups);
+  std::fill(result_matrix.begin(), result_matrix.end(), NA_REAL);
+  
+  for (int c = 0; c < n_groups; ++c) {
+    for (int r = 0; r < n_times; ++r) {
+      
+      // Identify the anchor time (Current Time - Year Lag)
+      int current_abs = get_abs_month(unique_years[r], unique_months[r]);
+      int anchor_abs = current_abs - (year_lag * 12);
+      
+      double sum_product = 0.0;
+      double sum_divisor = 0.0;
+      
+      // Define Window (based on Month Lag)
+      // Loop through the offsets for the window size
+      for (int w = 0; w < window_size; ++w) {
+        int lag_amount = month_lag + w; 
+        int target_abs = anchor_abs - lag_amount;
+        
+        // Check if this specific target month exists in our data
+        auto it = time_to_row_map.find(target_abs);
+        
+        if (it != time_to_row_map.end()) {
+          int target_row = it->second;
+          double val = values_matrix(target_row, c);
+          
+          if (!NumericVector::is_na(val)) {
+            if (weighted) {
+              int weight = count_matrix(target_row, c);
+              sum_product += val * weight;
+              sum_divisor += weight;
+            } else {
+              sum_product += val;
+              sum_divisor += 1.0;
+            }
+          }
+        }
+      }
+      
+      if (sum_divisor > 0) {
+        result_matrix(r, c) = sum_product / sum_divisor;
+      }
+    }
+  }
+  
+  // Set Names
+  colnames(result_matrix) = unique_groups;
+  
+  return result_matrix;
+}