Feat/optimization

.gitignore

@@ -25,6 +25,9 @@ src/MultisetMedian.cpp
 src/SlidingWelford.cpp
 src/SlidingWelfordRing.cpp
 src/SlidingMoments.cpp
+src/SlidingMean.cpp
+src/SlidingCovariance.cpp
+src/SlidingMomentsPrefix.cpp
 
 inst/include/
 

CMakeLists.txt

@@ -13,6 +13,7 @@ add_library(rolling_core
     src_core/MonotonicMin.cpp
     src_core/SlidingWelford.cpp
     src_core/SlidingWelfordRing.cpp
+    src_core/SlidingMean.cpp
     src_core/SlidingMoments.cpp
     src_core/SlidingCovariance.cpp
 )

R/rolling_metrics.R

@@ -1,16 +1,3 @@
-# Count non-NA values inside each rolling window using vectorised cumsum.
-.count_non_na <- function(x, window_size) {
-  cum <- cumsum(!is.na(x))
-  lagged <- c(rep(0L, window_size), cum)[seq_along(x)]
-  cum - lagged
-}
-
-.apply_min_periods <- function(result, x, window_size, min_periods) {
-  if (min_periods == 0L || length(x) == 0L) return(result)
-  result[.count_non_na(x, window_size) < min_periods] <- NA_real_
-  result
-}
-
 .check_window <- function(window_size) {
   if (!is.finite(window_size) || window_size < 1L)
     stop("window_size must be a positive finite integer.", call. = FALSE)
@@ -33,6 +20,9 @@
 #' @param min_periods Minimum number of non-\code{NA} observations required in
 #'   a window to return a result. Defaults to \code{window_size} (pandas
 #'   semantics). Positions with fewer non-\code{NA} values yield \code{NA}.
+#' @param method \code{"stable"} (default) uses the Welford online algorithm.
+#'   \code{"fast"} uses a prefix-sum approach (faster, but susceptible to
+#'   catastrophic cancellation when values are large and variance is small).
 #'
 #' @return
 #' A numeric vector with rolling sample variance values. Entries are
@@ -45,13 +35,17 @@
 #' @examples
 #' x <- as.double(c(1, 2, 3, 4))
 #' rolling_variance(x, 3L)
-rolling_variance <- function(x, window_size, min_periods = window_size) {
+rolling_variance <- function(x, window_size, min_periods = window_size,
+                             method = "stable") {
   x <- as.double(x)
   .check_window(window_size)
   mp <- .check_min_periods(min_periods, window_size)
-  result <- .Call("rolling_variance_c", x, as.integer(window_size),
-                  PACKAGE = "robustrolling")
-  .apply_min_periods(result, x, window_size, mp)
+  if (method == "fast")
+    .Call("rolling_variance_fast_c", x, as.integer(window_size), as.integer(mp),
+          PACKAGE = "robustrolling")
+  else
+    .Call("rolling_variance_c", x, as.integer(window_size), as.integer(mp),
+          PACKAGE = "robustrolling")
 }
 
 #' @title Rolling Maximum
@@ -76,9 +70,8 @@ rolling_max <- function(x, window_size, min_periods = window_size) {
   x <- as.double(x)
   .check_window(window_size)
   mp <- .check_min_periods(min_periods, window_size)
-  result <- .Call("rolling_max_c", x, as.integer(window_size),
-                  PACKAGE = "robustrolling")
-  .apply_min_periods(result, x, window_size, mp)
+  .Call("rolling_max_c", x, as.integer(window_size), as.integer(mp),
+        PACKAGE = "robustrolling")
 }
 
 #' @title Rolling Minimum
@@ -103,9 +96,8 @@ rolling_min <- function(x, window_size, min_periods = window_size) {
   x <- as.double(x)
   .check_window(window_size)
   mp <- .check_min_periods(min_periods, window_size)
-  result <- .Call("rolling_min_c", x, as.integer(window_size),
-                  PACKAGE = "robustrolling")
-  .apply_min_periods(result, x, window_size, mp)
+  .Call("rolling_min_c", x, as.integer(window_size), as.integer(mp),
+        PACKAGE = "robustrolling")
 }
 
 #' @title Rolling Median
@@ -131,9 +123,8 @@ rolling_median <- function(x, window_size, min_periods = window_size) {
   x <- as.double(x)
   .check_window(window_size)
   mp <- .check_min_periods(min_periods, window_size)
-  result <- .Call("rolling_median_c", x, as.integer(window_size),
-                  PACKAGE = "robustrolling")
-  .apply_min_periods(result, x, window_size, mp)
+  .Call("rolling_median_c", x, as.integer(window_size), as.integer(mp),
+        PACKAGE = "robustrolling")
 }
 
 #' @title Rolling Mean
@@ -145,6 +136,10 @@ rolling_median <- function(x, window_size, min_periods = window_size) {
 #' @param window_size Positive integer window length.
 #' @param min_periods Minimum number of non-\code{NA} observations required in
 #'   a window to return a result. Defaults to \code{window_size}.
+#' @param assume_finite If \code{TRUE}, assumes the input contains no
+#'   \code{NA} values and uses a faster SIMD prefix-sum path. Passing
+#'   \code{TRUE} when \code{NA}s are present produces incorrect results.
+#'   Defaults to \code{FALSE}.
 #'
 #' @return A numeric vector with rolling mean values.
 #'
@@ -153,13 +148,13 @@ rolling_median <- function(x, window_size, min_periods = window_size) {
 #' @examples
 #' x <- as.double(c(1, 2, 3, 4))
 #' rolling_mean(x, 3L)
-rolling_mean <- function(x, window_size, min_periods = window_size) {
+rolling_mean <- function(x, window_size, min_periods = window_size,
+                         assume_finite = FALSE) {
   x <- as.double(x)
   .check_window(window_size)
   mp <- .check_min_periods(min_periods, window_size)
-  result <- .Call("rolling_mean_c", x, as.integer(window_size),
-                  PACKAGE = "robustrolling")
-  .apply_min_periods(result, x, window_size, mp)
+  .Call("rolling_mean_c", x, as.integer(window_size), as.integer(mp),
+        as.logical(assume_finite), PACKAGE = "robustrolling")
 }
 
 #' @title Rolling Skewness
@@ -172,6 +167,9 @@ rolling_mean <- function(x, window_size, min_periods = window_size) {
 #' @param window_size Positive integer window length.
 #' @param min_periods Minimum number of non-\code{NA} observations required in
 #'   a window to return a result. Defaults to \code{window_size}.
+#' @param method \code{"stable"} (default) uses Terriberry's online algorithm.
+#'   \code{"fast"} uses a prefix-sum approach (faster, but susceptible to
+#'   catastrophic cancellation when values are large and variance is small).
 #'
 #' @return A numeric vector with rolling skewness values.
 #'
@@ -180,13 +178,17 @@ rolling_mean <- function(x, window_size, min_periods = window_size) {
 #' @examples
 #' x <- as.double(c(1, 2, 3, 4, 5))
 #' rolling_skewness(x, 3L)
-rolling_skewness <- function(x, window_size, min_periods = window_size) {
+rolling_skewness <- function(x, window_size, min_periods = window_size,
+                             method = "stable") {
   x <- as.double(x)
   .check_window(window_size)
   mp <- .check_min_periods(min_periods, window_size)
-  result <- .Call("rolling_skewness_c", x, as.integer(window_size),
-                  PACKAGE = "robustrolling")
-  .apply_min_periods(result, x, window_size, mp)
+  if (method == "fast")
+    .Call("rolling_skewness_fast_c", x, as.integer(window_size), as.integer(mp),
+          PACKAGE = "robustrolling")
+  else
+    .Call("rolling_skewness_c", x, as.integer(window_size), as.integer(mp),
+          PACKAGE = "robustrolling")
 }
 
 #' @title Rolling Kurtosis
@@ -199,6 +201,9 @@ rolling_skewness <- function(x, window_size, min_periods = window_size) {
 #' @param window_size Positive integer window length.
 #' @param min_periods Minimum number of non-\code{NA} observations required in
 #'   a window to return a result. Defaults to \code{window_size}.
+#' @param method \code{"stable"} (default) uses Terriberry's online algorithm.
+#'   \code{"fast"} uses a prefix-sum approach (faster, but susceptible to
+#'   catastrophic cancellation when values are large and variance is small).
 #'
 #' @return A numeric vector with rolling excess kurtosis values.
 #'
@@ -207,13 +212,17 @@ rolling_skewness <- function(x, window_size, min_periods = window_size) {
 #' @examples
 #' x <- as.double(c(1, 2, 3, 4, 5))
 #' rolling_kurtosis(x, 4L)
-rolling_kurtosis <- function(x, window_size, min_periods = window_size) {
+rolling_kurtosis <- function(x, window_size, min_periods = window_size,
+                             method = "stable") {
   x <- as.double(x)
   .check_window(window_size)
   mp <- .check_min_periods(min_periods, window_size)
-  result <- .Call("rolling_kurtosis_c", x, as.integer(window_size),
-                  PACKAGE = "robustrolling")
-  .apply_min_periods(result, x, window_size, mp)
+  if (method == "fast")
+    .Call("rolling_kurtosis_fast_c", x, as.integer(window_size), as.integer(mp),
+          PACKAGE = "robustrolling")
+  else
+    .Call("rolling_kurtosis_c", x, as.integer(window_size), as.integer(mp),
+          PACKAGE = "robustrolling")
 }
 
 #' @title Rolling Covariance
@@ -241,9 +250,8 @@ rolling_cov <- function(x, y, window_size, min_periods = window_size) {
   if (length(x) != length(y)) stop("x and y must have the same length.", call. = FALSE)
   .check_window(window_size)
   mp <- .check_min_periods(min_periods, window_size)
-  result <- .Call("rolling_cov_c", x, y, as.integer(window_size),
-                  PACKAGE = "robustrolling")
-  .apply_min_periods(result, x, window_size, mp)
+  .Call("rolling_cov_c", x, y, as.integer(window_size), as.integer(mp),
+        PACKAGE = "robustrolling")
 }
 
 #' @title Rolling Correlation
@@ -271,7 +279,6 @@ rolling_cor <- function(x, y, window_size, min_periods = window_size) {
   if (length(x) != length(y)) stop("x and y must have the same length.", call. = FALSE)
   .check_window(window_size)
   mp <- .check_min_periods(min_periods, window_size)
-  result <- .Call("rolling_cor_c", x, y, as.integer(window_size),
-                  PACKAGE = "robustrolling")
-  .apply_min_periods(result, x, window_size, mp)
+  .Call("rolling_cor_c", x, y, as.integer(window_size), as.integer(mp),
+        PACKAGE = "robustrolling")
 }

include/MonotonicMax.hpp

@@ -28,4 +28,5 @@ class MonotonicMax : public RollingMetric<MonotonicMax> {
   std::size_t window_size_;
   std::size_t current_tick_;
   std::deque<Element> deque_;
+  std::deque<std::size_t> non_nan_ticks_;
 };

include/MonotonicMin.hpp

@@ -28,4 +28,5 @@ class MonotonicMin : public RollingMetric<MonotonicMin> {
   std::size_t window_size_;
   std::size_t current_tick_;
   std::deque<Element> deque_;
+  std::deque<std::size_t> non_nan_ticks_;
 };

include/RollingMetric.hpp

@@ -21,7 +21,7 @@ template <typename Derived> class RollingMetric {
   void skip() { static_cast<Derived *>(this)->skip_impl(); }
 
   void process_batch(const double *input_data, std::size_t length,
-                     double *output_data) {
+                     double *output_data, std::size_t min_periods = 0) {
     for (std::size_t i = 0; i < length; ++i) {
       if (std::isnan(input_data[i])) {
         this->skip();
@@ -30,6 +30,8 @@ template <typename Derived> class RollingMetric {
         this->update(input_data[i]);
         output_data[i] = this->get_value();
       }
+      if (this->current_size() < min_periods)
+        output_data[i] = std::numeric_limits<double>::quiet_NaN();
     }
   }
 

include/SlidingCovariance.hpp

@@ -8,6 +8,7 @@ class SlidingCovariance {
 public:
   explicit SlidingCovariance(std::size_t window_size);
   void update(double x, double y);
+  std::size_t current_size() const;
   double get_covariance() const;
   double get_correlation() const;
   double get_mean_x() const;

include/SlidingMean.hpp

@@ -0,0 +1,32 @@
+#pragma once
+#include "RollingMetric.hpp"
+#include <cmath>
+#include <cstddef>
+#include <limits>
+#include <stdexcept>
+#include <vector>
+
+class SlidingMean : public RollingMetric<SlidingMean> {
+  friend class RollingMetric<SlidingMean>;
+
+public:
+  explicit SlidingMean(std::size_t window_size);
+  double get_mean() const;
+  void fast_mean_batch(const double *in, std::size_t n, double *out,
+                       std::size_t min_periods) const;
+  void fast_mean_batch_finite(const double *in, std::size_t n, double *out,
+                              std::size_t min_periods) const;
+
+private:
+  void update_impl(double value);
+  void skip_impl();
+  double get_value_impl() const;
+  std::size_t current_size_impl() const;
+
+  std::size_t window_size_;
+  std::vector<double> buffer_;
+  std::size_t head_{0};
+  std::size_t n_written_{0};
+  std::size_t count_{0};
+  double sum_{0.0};
+};

include/SlidingMoments.hpp

@@ -8,6 +8,7 @@ class SlidingMoments {
 public:
   explicit SlidingMoments(std::size_t window_size);
   void update(double x);
+  void update_skewness_only(double value);
   void reset();
   std::size_t current_size() const;
   double get_skewness() const; // adjusted Fisher-Pearson

include/SlidingMomentsPrefix.hpp

@@ -0,0 +1,20 @@
+#pragma once
+#include <cmath>
+#include <cstddef>
+#include <limits>
+#include <vector>
+
+class SlidingMomentsPrefix {
+public:
+  explicit SlidingMomentsPrefix(std::size_t window_size);
+
+  void variance_batch(const double *in, std::size_t n, double *out,
+                      std::size_t min_periods) const;
+  void skewness_batch(const double *in, std::size_t n, double *out,
+                      std::size_t min_periods) const;
+  void kurtosis_batch(const double *in, std::size_t n, double *out,
+                      std::size_t min_periods) const;
+
+private:
+  std::size_t window_size_;
+};

inst/tinytest/test_rolling_moments.R

@@ -148,3 +148,49 @@ expect_equal(robustrolling::rolling_kurtosis(numeric(0), 4L), numeric(0))
 # Input validation
 expect_true(is.double(robustrolling::rolling_kurtosis(1:5, 4L)))
 expect_error(robustrolling::rolling_kurtosis(as.double(1:5), 0L))
+
+# ---- method = "fast" (SlidingMomentsPrefix) ----------------------------------
+
+x_fast <- as.double(c(1, 3, -1, 5, 2, 8, 0, 4))
+
+# fast skewness matches stable on well-conditioned data
+for (k in c(4L, 5L, 6L)) {
+  expect_equal(
+    robustrolling::rolling_skewness(x_fast, k, min_periods = 1L, method = "fast"),
+    robustrolling::rolling_skewness(x_fast, k, min_periods = 1L),
+    tolerance = 1e-9
+  )
+}
+
+# fast kurtosis matches stable on well-conditioned data
+for (k in c(4L, 5L, 6L)) {
+  expect_equal(
+    robustrolling::rolling_kurtosis(x_fast, k, min_periods = 1L, method = "fast"),
+    robustrolling::rolling_kurtosis(x_fast, k, min_periods = 1L),
+    tolerance = 1e-9
+  )
+}
+
+# warmup NAs — skewness needs 3 obs, kurtosis needs 4
+x_lin <- as.double(1:10)
+expect_true(all(is.na(robustrolling::rolling_skewness(x_lin, 5L, method = "fast")[1:4])))
+expect_true(all(is.na(robustrolling::rolling_kurtosis(x_lin, 5L, method = "fast")[1:4])))
+
+# min_periods respected
+x_mp <- as.double(c(1, 2, 3, 4, 5))
+out_mp <- robustrolling::rolling_skewness(x_mp, 5L, min_periods = 3L, method = "fast")
+expect_true(is.na(out_mp[1]))
+expect_true(is.na(out_mp[2]))
+expect_false(is.na(out_mp[3]))
+
+# NA in input handled correctly
+x_na2 <- as.double(c(1, 2, NA_real_, 4, 5, 6, 7, 8))
+expect_equal(
+  robustrolling::rolling_skewness(x_na2, 5L, min_periods = 1L, method = "fast"),
+  robustrolling::rolling_skewness(x_na2, 5L, min_periods = 1L),
+  tolerance = 1e-9
+)
+
+# empty input
+expect_equal(robustrolling::rolling_skewness(numeric(0), 3L, method = "fast"), numeric(0))
+expect_equal(robustrolling::rolling_kurtosis(numeric(0), 4L, method = "fast"), numeric(0))