Release 1.3.1: polish, safety, and benchmark updates.

- idx_to_letter bounds check [1, 26]
- _is_equal_mindist(std::string) in SAX numerosity reduction
- Remove unused HOT-SAX state; fix cosine_dist roxygen
- Extend bench_discords.R with RRA and SAX cases
- NEWS and version bump summarizing Phases 1-4

Author: seninp

DESCRIPTION

@@ -1,5 +1,5 @@
 Package: jmotif
-Version: 1.3.0
+Version: 1.3.1
 Encoding: UTF-8
 Title: Time Series Analysis Toolkit Based on Symbolic Aggregate Discretization, i.e. SAX
 Description: Implements time series z-normalization, SAX, HOT-SAX, VSM, SAX-VSM, RePair, and RRA

NEWS.md

@@ -1,3 +1,18 @@
+# jmotif 1.3.1
+
+* Memory and performance (Phases 1–4):
+  - Fix memory leaks in RePair grammar (`rule_record` by-value), RRA, and
+    `cosine_sim` (Valgrind-clean on long runs).
+  - HOT-SAX: precompute z-normed sliding windows (~1.4× on ecg0606, ~17× on
+    longer series); identical discord positions and distances.
+  - SAX sliding window: `_paa2` in-place fractional PAA, `_znorm_slice`,
+    unified `_sax_via_window` path (~2.3× on bag construction).
+  - RRA: fused normalized distance, precomputed `w_size` z-norm windows
+    (~1.3× on ecg0606); identical results and distance-call counts.
+* Safety: `idx_to_letter()` validates index range [1, 26].
+* Internal: use `_is_equal_mindist(std::string)` in SAX numerosity reduction;
+  remove unused HOT-SAX state; clarify `cosine_dist()` roxygen.
+
 # Version 1.3.0
 * Cross-implementation alignment with the saxpy (Python) reference and the
   Matrix Profile / MASS convention:

R/jmotif.R

@@ -70,10 +70,10 @@ min_dist <- function(str1, str2, alphabet_size, compression_ratio = 1) {
   }
 }
 
-#' Computes the cosine similarity between numeric vectors
+#' Computes the cosine distance between numeric vectors (1 minus cosine similarity).
 #'
-#' @param m the data matrix
-#' @return Returns the cosine similarity
+#' @param m the data matrix (rows are vectors).
+#' @return A \code{dist} object of pairwise cosine distances.
 #' @export
 #' @examples
 #' a <- c(2, 1, 0, 2, 0, 1, 1, 1)

bench_discords.R

@@ -0,0 +1,145 @@
+#!/usr/bin/env Rscript
+# Benchmark discord discovery — primarily HOT-SAX (Phase 2 z-norm precompute).
+#
+# Usage:
+#   Rscript bench_discords.R [label] [out_csv]
+#
+# Example before/after on poptiplex:
+#   git checkout 37353dc && R CMD INSTALL --library=$HOME/R/library .
+#   Rscript inst/site/bench_discords.R pre-phase2
+#   git checkout master && R CMD INSTALL --library=$HOME/R/library .
+#   Rscript inst/site/bench_discords.R post-phase2
+#   Rscript inst/site/bench_discords.R --compare bench_discords_results.csv
+
+suppressPackageStartupMessages(library(jmotif))
+
+args <- commandArgs(trailingOnly = TRUE)
+
+if (length(args) >= 1 && args[[1]] == "--compare") {
+  csv <- if (length(args) >= 2) args[[2]] else "bench_discords_results.csv"
+  if (!file.exists(csv)) stop("CSV not found: ", csv)
+  d <- read.csv(csv, stringsAsFactors = FALSE)
+  agg <- aggregate(
+    elapsed_sec ~ label + benchmark,
+    data = d,
+    FUN = function(x) c(mean = mean(x), sd = sd(x), min = min(x), max = max(x))
+  )
+  flat <- data.frame(
+    label = agg$label,
+    benchmark = agg$benchmark,
+    mean_sec = agg$elapsed_sec[, "mean"],
+    sd_sec = agg$elapsed_sec[, "sd"],
+    min_sec = agg$elapsed_sec[, "min"],
+    max_sec = agg$elapsed_sec[, "max"],
+    stringsAsFactors = FALSE
+  )
+  labels <- unique(flat$label)
+  if (length(labels) >= 2) {
+    wide <- reshape(flat[, c("label", "benchmark", "mean_sec")],
+                    idvar = "benchmark", timevar = "label", direction = "wide")
+    pre <- grep("pre", labels, value = TRUE)[1]
+    post <- grep("post", labels, value = TRUE)[1]
+    if (is.na(pre)) pre <- labels[1]
+    if (is.na(post)) post <- labels[length(labels)]
+    pre_col <- paste0("mean_sec.", pre)
+    post_col <- paste0("mean_sec.", post)
+    if (pre_col %in% names(wide) && post_col %in% names(wide)) {
+      wide$speedup <- wide[[pre_col]] / wide[[post_col]]
+      cat("\n=== Speedup (", pre, " / ", post, " = times faster) ===\n", sep = "")
+      print(wide[, c("benchmark", pre_col, post_col, "speedup")], row.names = FALSE)
+    }
+  }
+  cat("\n=== Summary by label ===\n")
+  print(flat, row.names = FALSE)
+  quit(save = "no", status = 0)
+}
+
+label <- if (length(args) >= 1) args[[1]] else format(Sys.time(), "%Y%m%d-%H%M")
+out_csv <- if (length(args) >= 2) args[[2]] else "bench_discords_results.csv"
+
+benchmark_one <- function(name, expr, n = 5, warmup = 2) {
+  for (i in seq_len(warmup)) eval(expr, envir = parent.frame())
+  gc(verbose = FALSE)
+
+  times <- numeric(n)
+  for (i in seq_len(n)) {
+    gc(verbose = FALSE)
+    t0 <- proc.time()
+    eval(expr, envir = parent.frame())
+    times[i] <- (proc.time() - t0)[["elapsed"]]
+  }
+
+  data.frame(
+    label = label,
+    benchmark = name,
+    n = n,
+    elapsed_sec = times,
+    timestamp = format(Sys.time(), "%Y-%m-%d %H:%M:%S"),
+    stringsAsFactors = FALSE
+  )
+}
+
+x <- ecg0606
+x_medium <- rep(ecg0606, 3L)
+
+cat("=== jmotif discord benchmarks ===\n")
+cat("label:", label, "\n")
+cat("package:", as.character(packageVersion("jmotif")), "\n")
+cat("ecg0606 length:", length(x), "\n")
+cat("medium series length:", length(x_medium), "\n\n")
+
+cases <- list(
+  list(name = "hotsax_ecg0606_4disc",
+       expr = quote(find_discords_hotsax(x, 100, 4, 4, 0.01, 4)),
+       n = 10, warmup = 3),
+  list(name = "hotsax_ecg0606_5disc",
+       expr = quote(find_discords_hotsax(x, 100, 4, 4, 0.01, 5)),
+       n = 10, warmup = 3),
+  list(name = "hotsax_medium_4disc",
+       expr = quote(find_discords_hotsax(x_medium, 100, 4, 4, 0.01, 4)),
+       n = 5, warmup = 2),
+  list(name = "brute_ecg0606_4disc",
+       expr = quote(find_discords_brute_force(x, 100, 4, 0.01)),
+       n = 5, warmup = 2),
+  list(name = "wordbag_ecg0606",
+       expr = quote(series_to_wordbag(x, 60, 6, 6, "exact", 0.01)),
+       n = 10, warmup = 3),
+  list(name = "rra_ecg0606_4disc",
+       expr = quote(find_discords_rra(x, 100, 4, 4, "none", 0.01, 4)),
+       n = 5, warmup = 2),
+  list(name = "sax_medium_w160",
+       expr = quote(sax_via_window(x_medium, 160, 4, 4, "none", 0.01)),
+       n = 5, warmup = 2)
+)
+
+results <- NULL
+for (case in cases) {
+  cat("Running", case$name, "...\n")
+  flush.console()
+  chunk <- benchmark_one(case$name, case$expr, n = case$n, warmup = case$warmup)
+  print(aggregate(elapsed_sec ~ 1, data = chunk, FUN = mean))
+  results <- if (is.null(results)) chunk else rbind(results, chunk)
+}
+
+agg <- aggregate(elapsed_sec ~ benchmark, data = results, FUN = function(v) {
+  c(mean = mean(v), sd = sd(v), min = min(v), max = max(v))
+})
+summary_df <- data.frame(
+  benchmark = agg$benchmark,
+  mean_sec = round(agg$elapsed_sec[, "mean"], 4),
+  sd_sec = round(agg$elapsed_sec[, "sd"], 4),
+  min_sec = round(agg$elapsed_sec[, "min"], 4),
+  max_sec = round(agg$elapsed_sec[, "max"], 4)
+)
+
+cat("=== Run summary (label =", label, ") ===\n")
+print(summary_df, row.names = FALSE)
+
+if (file.exists(out_csv)) {
+  existing <- read.csv(out_csv, stringsAsFactors = FALSE)
+  combined <- rbind(existing, results)
+} else {
+  combined <- results
+}
+write.csv(combined, out_csv, row.names = FALSE)
+cat("\nAppended to", out_csv, "(", nrow(combined), "rows )\n")

inst/site/bench_discords.R

@@ -103,7 +103,13 @@ cases <- list(
        n = 5, warmup = 2),
   list(name = "wordbag_ecg0606",
        expr = quote(series_to_wordbag(x, 60, 6, 6, "exact", 0.01)),
-       n = 10, warmup = 3)
+       n = 10, warmup = 3),
+  list(name = "rra_ecg0606_4disc",
+       expr = quote(find_discords_rra(x, 100, 4, 4, "none", 0.01, 4)),
+       n = 5, warmup = 2),
+  list(name = "sax_medium_w160",
+       expr = quote(sax_via_window(x_medium, 160, 4, 4, "none", 0.01)),
+       n = 5, warmup = 2)
 )
 
 results <- NULL

jmotif.R

@@ -0,0 +1,84 @@
+#' @useDynLib jmotif
+#' @importFrom Rcpp sourceCpp
+NULL
+
+#' Generates a SAX MinDist distance matrix (i.e. the "lookup table") for a given alphabet size.
+#'
+#' @param a_size the desired alphabet size (a value between 2 and 20, inclusive)
+#' @return Returns a distance matrix (for SAX minDist) for a specified alphabet size
+#' @export
+#' @references Lonardi, S., Lin, J., Keogh, E., Patel, P.,
+#' Finding motifs in time series.
+#' In Proc. of the 2nd Workshop on Temporal Data Mining (pp. 53-68).
+#' @examples
+#' sax_distance_matrix(5)
+sax_distance_matrix <- function(a_size) {
+  if (a_size > 1 && a_size <= 20) {
+    cutlines <- jmotif::alphabet_to_cuts(a_size)[2:a_size]
+    distance_matrix <- matrix(rep(0, a_size * a_size), byrow = T, nrow = a_size, ncol = a_size)
+    i <- 1
+    while (i <= a_size) {
+      # the min_dist for adjacent symbols are 0, so we start with i+2
+      j <- i + 2;
+      while (j <= a_size) {
+        # square the distance now for future use
+        distance_matrix[i,j] <- (cutlines[i] - cutlines[j - 1]) * (cutlines[i] - cutlines[j - 1])
+        # the distance matrix is symmetric
+        distance_matrix[j,i] <- distance_matrix[i,j]
+        j <- j + 1
+      }
+      i <- i + 1
+    }
+    distance_matrix
+  } else {
+    stop(paste("unable to get a distance matrix for the alphabet size",a_size))
+  }
+}
+
+#' Computes the mindist value for two strings
+#'
+#' @param str1 the first string
+#' @param str2 the second string
+#' @param alphabet_size the used alphabet size
+#' @param compression_ratio the distance compression ratio
+#' @return Returns the distance between strings
+#' @export
+#' @references Lonardi, S., Lin, J., Keogh, E., Patel, P.,
+#' Finding motifs in time series.
+#' In Proc. of the 2nd Workshop on Temporal Data Mining (pp. 53-68).
+#' @examples
+#' str1 <- c('a', 'b', 'c')
+#' str2 <- c('c', 'b', 'a')
+#' min_dist(str1, str2, 3)
+min_dist <- function(str1, str2, alphabet_size, compression_ratio = 1) {
+  if (length(str1) != length(str2)) {
+    stop("error: the strings must have equal length")
+  }else{
+    if ( any(letters_to_idx(str1) > alphabet_size) |
+        any(letters_to_idx(str2) > alphabet_size)) {
+           stop('error: some symbol(s) in the string(s) exceed(s)
+                the alphabet size!');
+    }else{
+      dist_table <- sax_distance_matrix(alphabet_size)
+      dist <- 0
+      dist <- sqrt(
+         compression_ratio *
+          sum(diag(dist_table[letters_to_idx(str1), letters_to_idx(str2)]) ^ 2)
+        )
+      dist
+    }
+  }
+}
+
+#' Computes the cosine distance between numeric vectors (1 minus cosine similarity).
+#'
+#' @param m the data matrix (rows are vectors).
+#' @return A \code{dist} object of pairwise cosine distances.
+#' @export
+#' @examples
+#' a <- c(2, 1, 0, 2, 0, 1, 1, 1)
+#' b <- c(2, 1, 1, 1, 1, 0, 1, 1)
+#' sim <- cosine_dist(rbind(a,b))
+cosine_dist <- function(m) {
+  stats::as.dist(1 - m %*% t(m) / (sqrt(rowSums(m ^ 2) %*% t(rowSums(m ^ 2)))))
+}

src/hot-sax.cpp

@@ -11,7 +11,6 @@ discord_record find_best_discord_hotsax(std::vector<double>* ts, int w_size, dou
 
   double best_so_far_distance = 0;
   int best_so_far_index = -1;
-  std::string best_so_far_word = "";
 
   VisitRegistry outerRegistry(ts->size() - w_size, seed);
 
@@ -99,7 +98,6 @@ discord_record find_best_discord_hotsax(std::vector<double>* ts, int w_size, dou
           (nnDistance == best_so_far_distance && candidate_idx < best_so_far_index))){
         best_so_far_distance = nnDistance;
         best_so_far_index = candidate_idx;
-        best_so_far_word = it->second;
         //Rcout << "updated discord record: "<< best_so_far_word << " at " << best_so_far_index <<
         //  " nnDistance " << best_so_far_distance << "\n";
       }

src/sax.cpp

@@ -271,7 +271,7 @@ std::unordered_map<int, std::string> _sax_via_window(
       if (strategy_exact && old_str==curr_str) {
         continue;
       }
-      else if (strategy_mindist && is_equal_mindist(old_str, curr_str) ) {
+      else if (strategy_mindist && _is_equal_mindist(old_str, curr_str) ) {
         continue;
       }
     }

src/string.cpp

@@ -11,7 +11,10 @@
 //' idx_to_letter(2)
 // [[Rcpp::export]]
 char idx_to_letter(int idx) {
-  return LETTERS[idx-1];
+  if (idx < 1 || idx > 26) {
+    stop("'idx' is out of range [1, 26]");
+  }
+  return LETTERS[idx - 1];
 }
 
 //' Get the index for an ASCII letter.
@@ -85,18 +88,17 @@ bool is_equal_mindist(CharacterVector a, CharacterVector b) {
   return true;
 }
 
-/* bool _is_equal_mindist(std::string a, std::string b) {
+bool _is_equal_mindist(std::string a, std::string b) {
   if(a.length() != b.length()){
     return false;
-  }else{
-    for(unsigned i=0; i<a.length(); i++){
-      if( abs(a[i] - b[i]) > 1 ){
-        return false;
-      }
+  }
+  for(unsigned i=0; i<a.length(); i++){
+    if( abs(a[i] - b[i]) > 1 ){
+      return false;
     }
   }
   return true;
-}*/
+}
 
 int _count_spaces(std::string *s) {
   int count = 0;

test_str.R

@@ -0,0 +1,18 @@
+test_that("idx to letter", {
+  expect_true(is.character(idx_to_letter(3)))
+  expect_equal(idx_to_letter(3), "c")
+  expect_error(idx_to_letter(0))
+  expect_error(idx_to_letter(27))
+})
+
+test_that("letter to idx", {
+  expect_true(is.numeric(letter_to_idx('c')))
+  expect_equal(letter_to_idx('c'), 3)
+})
+
+test_that("letters to idx", {
+  expect_true(is.numeric(letters_to_idx(c('a', 'b', 'c', 'a'))))
+  expect_true(is.vector(letters_to_idx(c('a', 'b', 'c', 'a'))))
+  expect_equal(length(letters_to_idx(c('a', 'b', 'c', 'a'))), 4)
+  expect_equal(letters_to_idx(c('a', 'b', 'c', 'a'))[3], 3)
+})