perf: reduce SAX sliding-window allocations (Phase 3).

Use _paa2 with in-place segment sums, _znorm_slice for window z-norm
without copying, and a single _sax_via_window path shared with sax_via_window.
Removes deprecated _paa.

Author: seninp

inst/include/jmotif.h

@@ -326,8 +326,9 @@ std::vector<double> _alphabet_to_cuts(int a_size);
 int _count_spaces(std::string *s);
 double _mean(std::vector<double> *ts, int *start, int *end);
 std::vector<double> _znorm(const std::vector<double>& ts, double threshold);
-std::vector<double> _paa(std::vector<double> ts, int paa_num); // deprecated
-std::vector<double> _paa2(std::vector<double> ts, int paa_num);
+void _znorm_slice(const std::vector<double>& ts, int start, int end,
+                  double threshold, std::vector<double>& out);
+std::vector<double> _paa2(const std::vector<double>& ts, int paa_num);
 double _euclidean_dist(std::vector<double>* seq1, std::vector<double>* seq2);
 double _early_abandoned_dist(std::vector<double>* seq1, std::vector<double>* seq2,
                              double upper_limit);
@@ -337,7 +338,7 @@ std::string _series_to_string(std::vector<double> ts, int a_size);
 bool _is_equal_mindist(std::string a, std::string b);
 //
 std::unordered_map<int, std::string> _sax_via_window(
-    std::vector<double> ts, int w_size, int paa_size, int a_size,
+    const std::vector<double>& ts, int w_size, int paa_size, int a_size,
     std::string nr_strategy, double n_threshold);
 
 #endif

src/paa.cpp

@@ -29,146 +29,54 @@ NumericVector paa(NumericVector ts, int paa_num) {
   return wrap(_paa2(Rcpp::as< std::vector<double> >(ts), paa_num));
 }
 
-std::vector<double> _paa(std::vector<double> ts, int paa_num) {
+std::vector<double> _paa2(const std::vector<double>& ts, int paa_num) {
 
-  // fix the length
   int len = ts.size();
 
-  // check for the trivial case
-  if (len == paa_num) {
-    std::vector<double> res(ts);
-    return res;
-  }
-  else {
-    // if the number of points in a segment is even
-    if (len % paa_num == 0) {
-      std::vector<double> res(paa_num, 0);
-      int inc = len / paa_num;
-      for (int i = 0; i < len; i++) {
-        int idx = i / inc; // the spot
-        res[idx] = res[idx] + ts[i];
-      }
-      double dl = (double) inc;
-      for (int i = 0; i < paa_num; i++) {
-        res[i] = res[i] / dl;
-      }
-      return res;
-    }else{
-      // if the number of points in a segment is odd
-      std::vector<double> res(paa_num);
-      for (int i = 0; i < len * paa_num; i++) {
-        int idx = i / len; // the spot
-        int pos = i / paa_num; // the col spot
-        res[idx] = res[idx] + ts[pos];
-      }
-      double dl = (double) len;
-      for (int i = 0; i < paa_num; i++) {
-        res[i] = res[i] / dl;
-      }
-      return res;
-    }
-  }
-}
-
-std::vector<double> _paa2(std::vector<double> ts, int paa_num) {
-
-  std::vector<double> tempVector;
-
-  // fix the length
-  int len = ts.size();
-  // Rcout << "len " << len << std::endl;
-
   if(len < paa_num){
     stop("'paa_num' size is invalid");
   }
 
-  // check for the trivial case
   if (len == paa_num) {
-    std::vector<double> res(ts);
-    return res;
+    return ts;
   }
-  else {
 
-    // Rcpp::Rcout.precision(5);
-    // Rcout << std::fixed;
+  std::vector<double> res(paa_num, 0.0);
+  double points_per_segment = (double) len / (double) paa_num;
 
-    std::vector<double> res(paa_num, 0.0);
-    double points_per_segment = (double) len / (double) paa_num;
-    // Rcout << "points per seg " << points_per_segment << std::endl;
+  std::vector<double> breaks(paa_num + 1, 0);
+  for(int i = 0; i < paa_num + 1; i++){
+    breaks[i] = i * points_per_segment;
+  }
+  breaks[paa_num] = (double) len;
 
-    std::vector<double> breaks(paa_num + 1, 0);
-    for(int i = 0; i < paa_num + 1; i++){
-      breaks[i] = i * points_per_segment;
-    }
-    // The final break is paa_num * (len/paa_num), which should be exactly len
-    // but in IEEE-754 often rounds to len + epsilon (e.g. 7*(29/7) =
-    // 29.000000000000004). That makes floor(seg_end) == len for the last
-    // segment, so frac_end collapses to ~1e-15 and the last sample is dropped
-    // from the segment sum while the divisor stays the same -- corrupting the
-    // final PAA value. Snap it to len so the boundary is handled exactly,
-    // matching the saxpy / Java fractional PAA.
-    breaks[paa_num] = (double) len;
-    // Rcout << "ts length: " << len << ", breaks: ";
-    // for(auto it=breaks.begin(); it<breaks.end(); ++it){
-    //   Rcout << *it << ", ";
-    //}
-    //Rcout << std::endl;
+  for(int i = 0; i < paa_num; i++){
 
-    for(int i = 0; i < paa_num; i++){
+    double seg_start = breaks[i];
+    double seg_end = breaks[i+1];
 
-      double seg_start = breaks[i];
-      double seg_end = breaks[i+1];
-      // Rcout << " * seg_start " << seg_start << ", end " << seg_end << std::endl;
+    double frac_begin = ceil(seg_start) - seg_start;
+    double frac_end = seg_end - floor(seg_end);
 
-      double frac_begin = ceil(seg_start) - seg_start;
-      double frac_end = seg_end - floor(seg_end);
-      // Rcout << " ** frac_begin " << frac_begin << ", frac_end " << frac_end << std::endl;
+    int full_begin = (int)floor(seg_start);
+    int full_end = (int)ceil(seg_end);
+    if(full_end > len) {
+      full_end = len;
+    }
 
-      int full_begin = floor(seg_start);
-      int full_end = ceil(seg_end);
-      if(full_end > len) {
-        full_end = len;
+    double sum_of_elems = 0.0;
+    for (int j = full_begin; j < full_end; j++) {
+      double v = ts[j];
+      if (j == full_begin && frac_begin > 0) {
+        v *= frac_begin;
       }
-      // Rcout << " *** full_begin " << full_begin << ", full_end " << full_end << std::endl;
-
-      // the following taken from
-      // https://stackoverflow.com/questions/421573/best-way-to-extract-a-subvector-from-a-vector
-      // but it may leak memory?
-      std::vector<double>::const_iterator first = ts.begin() + full_begin;
-      std::vector<double>::const_iterator last = ts.begin() +  full_end;
-      std::vector<double> segment(first, last);
-      //std::vector<double> segment(ts.begin() + full_begin, ts.begin() +  full_end);
-
-      // Rcout << "segment ";
-      // for(auto it=segment.begin(); it<segment.end(); ++it){
-        // Rcout << *it << ", ";
-      // }
-      // Rcout << std::endl;
-
-      if(frac_begin > 0){
-        segment[0] = segment[0] * frac_begin;
+      if (j == full_end - 1 && frac_end > 0) {
+        v *= frac_end;
       }
-
-      if(frac_end > 0){
-        segment[segment.size()-1] = segment[segment.size()-1] * frac_end;
-      }
-      // Rcout << "adj_segment ";
-      // for(auto it=segment.begin(); it<segment.end(); ++it){
-        // Rcout << *it << ", ";
-      // }
-      // Rcout << std::endl;
-
-      double sum_of_elems = 0.0;
-      for (double e : segment)
-        sum_of_elems += e;
-      // Rcout << " **** sum " << sum_of_elems << std::endl;
-
-      segment.swap(tempVector);
-
-      res[i] = sum_of_elems / points_per_segment;
-      // Rcout << " **** res[" << i << "] " << res[i] << std::endl;
-
+      sum_of_elems += v;
     }
-    return res;
+
+    res[i] = sum_of_elems / points_per_segment;
   }
+  return res;
 }

src/sax.cpp

@@ -236,20 +236,21 @@ std::string _series_to_string(std::vector<double> ts, int a_size) {
 }
 
 std::unordered_map<int, std::string> _sax_via_window(
-    std::vector<double> ts, int w_size, int paa_size, int a_size,
+    const std::vector<double>& ts, int w_size, int paa_size, int a_size,
     std::string nr_strategy, double n_threshold) {
 
   std::unordered_map<int, std::string> idx2word;
   bool strategy_exact = is_equal_str("exact", nr_strategy);
   bool strategy_mindist = is_equal_str("mindist", nr_strategy);
 
   std::string old_str;
-
-  //Rcout << "series length " << ts.size() << ", window " << w_size << std::endl;
+  std::vector<double> znorm_buf;
+  znorm_buf.reserve(w_size);
+  std::vector<double> paa_buf;
+  paa_buf.reserve(paa_size);
 
   for (unsigned i = 0; i <= ts.size() - w_size; i++) {
 
-    // check if NA is encountered
     int idx = i + w_size - 1;
     if(R_IsNA(ts[idx])) {
       size_t size = std::snprintf(nullptr, 0, "encountered an Na and stopped processing at %i", i + w_size - 1) + 1;
@@ -261,17 +262,10 @@ std::unordered_map<int, std::string> _sax_via_window(
       break;
     }
 
-    // subseries extraction
-    // std::vector<double>::const_iterator first = ts.begin() + i;
-    // std::vector<double>::const_iterator last = ts.begin() + i + w_size;
-    // std::vector<double> subSection(first, last);
-    std::vector<double> subSection = _subseries(&ts, i, i + w_size);
-
-    subSection = _znorm(subSection, n_threshold);
-
-    subSection = _paa(subSection, paa_size);
+    _znorm_slice(ts, (int)i, (int)i + w_size, n_threshold, znorm_buf);
+    paa_buf = _paa2(znorm_buf, paa_size);
 
-    std::string curr_str = _series_to_string(subSection, a_size);
+    std::string curr_str = _series_to_string(paa_buf, a_size);
 
     if (!(old_str.empty())) {
       if (strategy_exact && old_str==curr_str) {
@@ -287,7 +281,6 @@ std::unordered_map<int, std::string> _sax_via_window(
     old_str = curr_str;
   }
 
-
   return idx2word;
 }
 
@@ -319,74 +312,11 @@ std::map<int, std::string> sax_via_window(
     forward_exception_to_r(ex);
   }
 
-  std::vector<double> series = Rcpp::as< std::vector<double> > (ts);
-  bool strategy_exact = is_equal_str("exact", nr_strategy);
-  bool strategy_mindist = is_equal_str("mindist", nr_strategy);
-
-  std::map<int, std::string> idx2word;
-
-  // std::ofstream bw("test_outR.txt");
-
-  std::string old_str("");
-
-  for (unsigned i = 0; i <= series.size() - w_size; i++) {
-
-    // check if NA is encountered
-    int idx = i + w_size - 1;
-    if(R_IsNA(ts[idx])) {
-      size_t size = std::snprintf(nullptr, 0, "encountered an Na and stopped processing at %i", i + w_size - 1) + 1;
-      std::unique_ptr<char[]> buf( new char[ size ] );
-      std::snprintf( buf.get(), size, "encountered an Na and stopped processing at %i", i + w_size - 1 );
-      Rcpp::warning(
-        std::string( buf.get(), buf.get() + size - 1 )
-      );
-      break;
-    }
-
-    NumericVector subSection = subseries(ts, i, i + w_size);
-    std::vector<double>::const_iterator first = series.begin() + i;
-    std::vector<double>::const_iterator last = series.begin() + i + w_size;
-    std::vector<double> sub_section(first, last);
-
-    // bw << i << "\t[" ;
-    // for (auto i = subSection.begin(); i != subSection.end(); ++i)
-    //   bw << *i << ',';
-    // bw << "]\t" ;
-
-    sub_section = _znorm(sub_section, n_threshold);
-
-    sub_section = _paa(sub_section, paa_size);
-    // bw << "[" ;
-    // for (auto i = subSection.begin(); i != subSection.end(); ++i)
-    //   bw << *i << ',';
-    // bw << "]\t" ;
-
-    std::string curr_str = _series_to_string(sub_section, a_size);
-    // bw << curr_str << "\t" ;
-
-    if (!(0 == old_str.length())) {
-      if ( strategy_exact
-             && old_str == curr_str ) {
-        // std::cout << std::endl;
-        // bw << "skipped\n" ;
-        continue;
-      }
-      else if (strategy_mindist
-                 && is_equal_mindist(old_str, curr_str) ) {
-        // std::cout << std::endl;
-        continue;
-      }
-    }
-
-    idx2word.insert(std::make_pair(i, curr_str));
-
-    old_str = curr_str;
-
-    // std::cout << " " << curr_str << std::endl;
-    // bw << "kept\n" ;
-  }
-  // bw.close();
-  return idx2word;
+  std::vector<double> series = Rcpp::as<std::vector<double>>(ts);
+  std::string nr = Rcpp::as<std::string>(nr_strategy);
+  std::unordered_map<int, std::string> um = _sax_via_window(
+      series, w_size, paa_size, a_size, nr, n_threshold);
+  return std::map<int, std::string>(um.begin(), um.end());
 }
 
 //' Discretize a time series with SAX using chunking (no sliding window).

src/znorm.cpp

@@ -53,3 +53,33 @@ std::vector<double> _znorm(const std::vector<double>& ts, double threshold) {
   return diff;
 
 }
+
+void _znorm_slice(const std::vector<double>& ts, int start, int end,
+                  double threshold, std::vector<double>& out) {
+  int len = end - start;
+  out.resize(len);
+
+  double sum = 0.0;
+  for (int i = start; i < end; i++) {
+    sum += ts[i];
+  }
+  double mean = sum / len;
+
+  double sq_sum = 0.0;
+  for (int i = start; i < end; i++) {
+    double d = ts[i] - mean;
+    sq_sum += d * d;
+  }
+  double stdev = std::sqrt(sq_sum / len);
+
+  if (stdev < threshold) {
+    for (int i = start; i < end; i++) {
+      out[i - start] = ts[i];
+    }
+    return;
+  }
+
+  for (int i = start; i < end; i++) {
+    out[i - start] = (ts[i] - mean) / stdev;
+  }
+}