From fd1bb137d6a551a3fb988eecfdc8e40c03fe55df Mon Sep 17 00:00:00 2001 From: Pratik Bhadane Date: Tue, 7 Apr 2026 23:38:36 +0530 Subject: [PATCH] Dtw algo (#9) * feat: implement Dynamic Time Warping (DTW) functionality - Added DTW distance computation and optimal warping path functions in Rust. - Introduced corresponding Python bindings for DTW, DTW_DISTANCE, and BATCH_DTW. - Enhanced WASM support with a new dtw_distance function. - Included comprehensive unit tests for DTW functionality, validating against the dtaidistance library and ensuring mathematical properties. * chore: update ferro-ta version to 1.1.4 - Bumped version number of ferro-ta to 1.1.4 in uv.lock and Cargo.lock files. - Ensured consistency across package dependencies for the updated version. --- Cargo.lock | 52 +++--- Cargo.toml | 4 +- conda/meta.yaml | 2 +- crates/ferro_ta_core/Cargo.toml | 2 +- crates/ferro_ta_core/README.md | 2 +- crates/ferro_ta_core/src/statistic.rs | 200 ++++++++++++++++++++++++ docs/api_manifest.json | 70 ++++++++- docs/changelog.rst | 2 +- docs/support_matrix.rst | 2 +- pyproject.toml | 2 +- python/ferro_ta/indicators/statistic.py | 109 +++++++++++++ src/statistic/dtw.rs | 98 ++++++++++++ src/statistic/mod.rs | 4 + tests/unit/indicators/test_statistic.py | 178 +++++++++++++++++++++ uv.lock | 2 +- wasm/Cargo.lock | 4 +- wasm/Cargo.toml | 2 +- wasm/package.json | 2 +- wasm/src/lib.rs | 12 ++ 19 files changed, 703 insertions(+), 46 deletions(-) create mode 100644 src/statistic/dtw.rs diff --git a/Cargo.lock b/Cargo.lock index 104e6ca..1dd0dfd 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -34,9 +34,9 @@ checksum = "940b3a0ca603d1eade50a4846a2afffd5ef57a9feac2c0e2ec2e14f9ead76000" [[package]] name = "arc-swap" -version = "1.9.0" +version = "1.9.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a07d1f37ff60921c83bdfc7407723bdefe89b44b98a9b772f225c8f9d67141a6" +checksum = "6a3a1fd6f75306b68087b831f025c712524bcb19aad54e557b1129cfa0a2b207" dependencies = [ "rustversion", ] @@ -67,9 +67,9 @@ checksum = "37b2a672a2cb129a2e41c10b1224bb368f9f37a2b16b612598138befd7b37eb5" [[package]] name = "cc" -version = "1.2.57" +version = "1.2.59" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "7a0dd1ca384932ff3641c8718a02769f1698e7563dc6974ffd03346116310423" +checksum = "b7a4d3ec6524d28a329fc53654bbadc9bdd7b0431f5d65f1a56ffb28a1ee5283" dependencies = [ "find-msvc-tools", "shlex", @@ -207,7 +207,7 @@ checksum = "48c757948c5ede0e46177b7add2e67155f70e33c07fea8284df6576da70b3719" [[package]] name = "ferro_ta" -version = "1.1.3" +version = "1.1.4" dependencies = [ "criterion", "ferro_ta_core", @@ -222,7 +222,7 @@ dependencies = [ [[package]] name = "ferro_ta_core" -version = "1.1.3" +version = "1.1.4" dependencies = [ "criterion", "serde", @@ -279,9 +279,9 @@ checksum = "8f42a60cbdf9a97f5d2305f08a87dc4e09308d1276d28c869c684d7777685682" [[package]] name = "js-sys" -version = "0.3.91" +version = "0.3.94" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b49715b7073f385ba4bc528e5747d02e66cb39c6146efb66b781f131f0fb399c" +checksum = "2e04e2ef80ce82e13552136fabeef8a5ed1f985a96805761cbb9a2c34e7664d9" dependencies = [ "once_cell", "wasm-bindgen", @@ -289,9 +289,9 @@ dependencies = [ [[package]] name = "libc" -version = "0.2.183" +version = "0.2.184" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b5b646652bf6661599e1da8901b3b9522896f01e736bad5f723fe7a3a27f899d" +checksum = "48f5d2a454e16a5ea0f4ced81bd44e4cfc7bd3a507b61887c99fd3538b28e4af" [[package]] name = "log" @@ -595,9 +595,9 @@ checksum = "dc897dd8d9e8bd1ed8cdad82b5966c3e0ecae09fb1907d58efaa013543185d0a" [[package]] name = "rustc-hash" -version = "2.1.1" +version = "2.1.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "357703d41365b4b27c590e3ed91eabb1b663f07c4c084095e60cbed4362dff0d" +checksum = "94300abf3f1ae2e2b8ffb7b58043de3d399c73fa6f4b73826402a5c457614dbe" [[package]] name = "rustversion" @@ -729,9 +729,9 @@ dependencies = [ [[package]] name = "wasm-bindgen" -version = "0.2.114" +version = "0.2.117" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "6532f9a5c1ece3798cb1c2cfdba640b9b3ba884f5db45973a6f442510a87d38e" +checksum = "0551fc1bb415591e3372d0bc4780db7e587d84e2a7e79da121051c5c4b89d0b0" dependencies = [ "cfg-if", "once_cell", @@ -742,9 +742,9 @@ dependencies = [ [[package]] name = "wasm-bindgen-macro" -version = "0.2.114" +version = "0.2.117" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "18a2d50fcf105fb33bb15f00e7a77b772945a2ee45dcf454961fd843e74c18e6" +checksum = "7fbdf9a35adf44786aecd5ff89b4563a90325f9da0923236f6104e603c7e86be" dependencies = [ "quote", "wasm-bindgen-macro-support", @@ -752,9 +752,9 @@ dependencies = [ [[package]] name = "wasm-bindgen-macro-support" -version = "0.2.114" +version = "0.2.117" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "03ce4caeaac547cdf713d280eda22a730824dd11e6b8c3ca9e42247b25c631e3" +checksum = "dca9693ef2bab6d4e6707234500350d8dad079eb508dca05530c85dc3a529ff2" dependencies = [ "bumpalo", "proc-macro2", @@ -765,18 +765,18 @@ dependencies = [ [[package]] name = "wasm-bindgen-shared" -version = "0.2.114" +version = "0.2.117" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "75a326b8c223ee17883a4251907455a2431acc2791c98c26279376490c378c16" +checksum = "39129a682a6d2d841b6c429d0c51e5cb0ed1a03829d8b3d1e69a011e62cb3d3b" dependencies = [ "unicode-ident", ] [[package]] name = "web-sys" -version = "0.3.91" +version = "0.3.94" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "854ba17bb104abfb26ba36da9729addc7ce7f06f5c0f90f3c391f8461cca21f9" +checksum = "cd70027e39b12f0849461e08ffc50b9cd7688d942c1c8e3c7b22273236b4dd0a" dependencies = [ "js-sys", "wasm-bindgen", @@ -840,18 +840,18 @@ dependencies = [ [[package]] name = "zerocopy" -version = "0.8.47" +version = "0.8.48" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "efbb2a062be311f2ba113ce66f697a4dc589f85e78a4aea276200804cea0ed87" +checksum = "eed437bf9d6692032087e337407a86f04cd8d6a16a37199ed57949d415bd68e9" dependencies = [ "zerocopy-derive", ] [[package]] name = "zerocopy-derive" -version = "0.8.47" +version = "0.8.48" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "0e8bc7269b54418e7aeeef514aa68f8690b8c0489a06b0136e5f57c4c5ccab89" +checksum = "70e3cd084b1788766f53af483dd21f93881ff30d7320490ec3ef7526d203bad4" dependencies = [ "proc-macro2", "quote", diff --git a/Cargo.toml b/Cargo.toml index 630a9f5..f0c2b45 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -5,7 +5,7 @@ resolver = "2" [package] name = "ferro_ta" -version = "1.1.3" +version = "1.1.4" edition = "2021" description = "Rust-powered Python technical analysis library with a TA-Lib-compatible API" license = "MIT" @@ -30,7 +30,7 @@ ndarray = "0.16" rayon = "1.10" log = "0.4" pyo3-log = "0.12" -ferro_ta_core = { path = "crates/ferro_ta_core", version = "1.1.3", features = ["serde"] } +ferro_ta_core = { path = "crates/ferro_ta_core", version = "1.1.4", features = ["serde"] } [dev-dependencies] criterion = { version = "0.8", features = ["html_reports"] } diff --git a/conda/meta.yaml b/conda/meta.yaml index 1524ede..fdaa6c6 100644 --- a/conda/meta.yaml +++ b/conda/meta.yaml @@ -1,5 +1,5 @@ {% set name = "ferro-ta" %} -{% set version = "1.1.3" %} +{% set version = "1.1.4" %} package: name: {{ name|lower }} diff --git a/crates/ferro_ta_core/Cargo.toml b/crates/ferro_ta_core/Cargo.toml index 6a5b703..afe4370 100644 --- a/crates/ferro_ta_core/Cargo.toml +++ b/crates/ferro_ta_core/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "ferro_ta_core" -version = "1.1.3" +version = "1.1.4" edition = "2021" description = "Pure Rust core indicator library — no PyO3, no numpy dependency" license = "MIT" diff --git a/crates/ferro_ta_core/README.md b/crates/ferro_ta_core/README.md index 80caa1c..ba9a98b 100644 --- a/crates/ferro_ta_core/README.md +++ b/crates/ferro_ta_core/README.md @@ -13,7 +13,7 @@ PyO3, NumPy, or Python runtime dependency, which makes it a good fit for: ```toml [dependencies] -ferro_ta_core = "1.1.3" +ferro_ta_core = "1.1.4" ``` ## Design diff --git a/crates/ferro_ta_core/src/statistic.rs b/crates/ferro_ta_core/src/statistic.rs index 694429c..9fc2219 100644 --- a/crates/ferro_ta_core/src/statistic.rs +++ b/crates/ferro_ta_core/src/statistic.rs @@ -247,6 +247,118 @@ pub fn correl(real0: &[f64], real1: &[f64], timeperiod: usize) -> Vec { result } +// --------------------------------------------------------------------------- +// Dynamic Time Warping (DTW) +// --------------------------------------------------------------------------- + +/// Internal helper: build the full DTW accumulated-cost matrix. +/// +/// Local cost: `|s1[i] - s2[j]|` (Euclidean / L1 for 1-D series). +/// This matches the convention used by `dtaidistance.dtw.distance()`. +/// +/// Out-of-band cells (Sakoe-Chiba constraint) are set to `f64::INFINITY`. +fn dtw_matrix(s1: &[f64], s2: &[f64], window: Option) -> Vec> { + let n = s1.len(); + let m = s2.len(); + let mut dp = vec![vec![f64::INFINITY; m]; n]; + for i in 0..n { + // Window convention matches dtaidistance: window=w means |i-j| < w. + // None = unconstrained (full matrix). + let (j_lo, j_hi) = match window { + None => (0, m), + Some(w) => { + let lo = i.saturating_sub(w.saturating_sub(1)); + let hi = i.saturating_add(w).min(m); + (lo, hi) + } + }; + for j in j_lo..j_hi { + // Squared Euclidean local cost — matches dtaidistance convention. + // The final sqrt is applied only once at the top level (not per-step). + let cost = (s1[i] - s2[j]).powi(2); + let prev = if i == 0 && j == 0 { + 0.0 + } else if i == 0 { + dp[0][j - 1] + } else if j == 0 { + dp[i - 1][0] + } else { + dp[i - 1][j - 1].min(dp[i - 1][j]).min(dp[i][j - 1]) + }; + dp[i][j] = cost + prev; + } + } + dp +} + +/// Compute the Dynamic Time Warping distance between two 1-D series. +/// +/// Returns the accumulated Euclidean cost along the optimal warping path. +/// Uses `|s1[i] - s2[j]|` as the local cost, matching `dtaidistance` convention. +/// +/// # Arguments +/// * `s1` - First time series. +/// * `s2` - Second time series. +/// * `window` - Optional Sakoe-Chiba band width. `None` = unconstrained. +/// +/// Returns `f64::NAN` if either input is empty. +pub fn dtw_distance(s1: &[f64], s2: &[f64], window: Option) -> f64 { + if s1.is_empty() || s2.is_empty() { + return f64::NAN; + } + let dp = dtw_matrix(s1, s2, window); + // sqrt applied once at the end — matches dtaidistance.dtw.distance() convention. + dp[s1.len() - 1][s2.len() - 1].sqrt() +} + +/// Compute the DTW distance and the optimal warping path between two 1-D series. +/// +/// The warping path is a `Vec<(usize, usize)>` of `(i, j)` index pairs, +/// starting at `(0, 0)` and ending at `(n-1, m-1)`, monotonically non-decreasing. +/// +/// # Arguments +/// * `s1` - First time series. +/// * `s2` - Second time series. +/// * `window` - Optional Sakoe-Chiba band width. `None` = unconstrained. +/// +/// Returns `(f64::NAN, vec![])` if either input is empty. +pub fn dtw_path(s1: &[f64], s2: &[f64], window: Option) -> (f64, Vec<(usize, usize)>) { + if s1.is_empty() || s2.is_empty() { + return (f64::NAN, vec![]); + } + let dp = dtw_matrix(s1, s2, window); + let dist = dp[s1.len() - 1][s2.len() - 1].sqrt(); + + // Backtrace from (n-1, m-1) to (0, 0) + let mut path = Vec::new(); + let (mut i, mut j) = (s1.len() - 1, s2.len() - 1); + path.push((i, j)); + while i > 0 || j > 0 { + let (ni, nj) = match (i, j) { + (0, _) => (0, j - 1), + (_, 0) => (i - 1, 0), + _ => { + let diag = dp[i - 1][j - 1]; + let up = dp[i - 1][j]; + let left = dp[i][j - 1]; + let best = diag.min(up).min(left); + if best == diag { + (i - 1, j - 1) + } else if best == up { + (i - 1, j) + } else { + (i, j - 1) + } + } + }; + i = ni; + j = nj; + path.push((i, j)); + } + path.reverse(); + (dist, path) +} + #[cfg(test)] mod tests { use super::*; @@ -259,4 +371,92 @@ mod tests { assert!(v.abs() < 1e-10); } } + + #[test] + fn dtw_identical_series_is_zero() { + let a = vec![1.0, 2.0, 3.0, 4.0, 5.0]; + assert_eq!(dtw_distance(&a, &a, None), 0.0); + } + + #[test] + fn dtw_known_shifted_series() { + // [0,1,2] vs [1,2,3]: DTW uses squared Euclidean local cost + final sqrt. + // Optimal path (0,0)→(1,0)→(2,1)→(2,2), accumulated cost = 1+0+0+1 = 2, sqrt(2). + // Matches dtaidistance.dtw.distance([0,1,2],[1,2,3]) = 1.4142... + let a = vec![0.0, 1.0, 2.0]; + let b = vec![1.0, 2.0, 3.0]; + let expected = 2.0_f64.sqrt(); + let result = dtw_distance(&a, &b, None); + assert!( + (result - expected).abs() < 1e-12, + "got {result}, expected {expected}" + ); + } + + #[test] + fn dtw_known_even_shift() { + // [0,2,4] vs [1,3,5]: diagonal path, squared costs 1+1+1=3, sqrt(3). + // Matches dtaidistance.dtw.distance([0,2,4],[1,3,5]) = 1.7320... + let a = vec![0.0, 2.0, 4.0]; + let b = vec![1.0, 3.0, 5.0]; + let expected = 3.0_f64.sqrt(); + let result = dtw_distance(&a, &b, None); + assert!( + (result - expected).abs() < 1e-12, + "got {result}, expected {expected}" + ); + } + + #[test] + fn dtw_single_element() { + let a = vec![3.0]; + let b = vec![7.0]; + assert_eq!(dtw_distance(&a, &b, None), 4.0); + } + + #[test] + fn dtw_empty_returns_nan() { + assert!(dtw_distance(&[], &[1.0, 2.0], None).is_nan()); + assert!(dtw_distance(&[1.0, 2.0], &[], None).is_nan()); + } + + #[test] + fn dtw_path_endpoints() { + let a = vec![1.0, 2.0, 3.0, 4.0]; + let b = vec![1.5, 2.5, 3.5, 4.5]; + let (_, path) = dtw_path(&a, &b, None); + assert_eq!(path.first(), Some(&(0, 0))); + assert_eq!(path.last(), Some(&(3, 3))); + } + + #[test] + fn dtw_path_is_monotone() { + let a = vec![1.0, 3.0, 2.0, 5.0, 4.0]; + let b = vec![2.0, 1.0, 4.0, 3.0, 6.0]; + let (_, path) = dtw_path(&a, &b, None); + for k in 1..path.len() { + assert!(path[k].0 >= path[k - 1].0); + assert!(path[k].1 >= path[k - 1].1); + } + } + + #[test] + fn dtw_path_distance_matches_distance_only() { + let a = vec![1.0, 4.0, 2.0, 8.0, 3.0]; + let b = vec![2.0, 3.0, 7.0, 4.0, 5.0]; + let d1 = dtw_distance(&a, &b, None); + let (d2, _) = dtw_path(&a, &b, None); + assert!((d1 - d2).abs() < 1e-12); + } + + #[test] + fn dtw_window_constrained_ge_unconstrained() { + // window convention matches dtaidistance: Some(w) means |i-j| < w. + // A narrow window restricts warping, so constrained distance >= unconstrained. + let a: Vec = (0..20).map(|x| x as f64).collect(); + let b: Vec = (0..20).map(|x| x as f64 + 3.0).collect(); + let d_full = dtw_distance(&a, &b, None); + let d_narrow = dtw_distance(&a, &b, Some(3)); + assert!(d_narrow >= d_full - 1e-12); + } } diff --git a/docs/api_manifest.json b/docs/api_manifest.json index 59bf48d..0ff3fea 100644 --- a/docs/api_manifest.json +++ b/docs/api_manifest.json @@ -1,8 +1,8 @@ { "surfaces": { "python": { - "indicator_count": 208, - "method_count": 464, + "indicator_count": 211, + "method_count": 467, "categories": [ "aggregation", "alerts", @@ -131,6 +131,13 @@ "doc": "", "params": [] }, + { + "name": "BATCH_DTW", + "category": "statistic", + "module": "ferro_ta.indicators.statistic", + "doc": "", + "params": [] + }, { "name": "BBANDS", "category": "overlap", @@ -656,6 +663,20 @@ "doc": "", "params": [] }, + { + "name": "DTW", + "category": "statistic", + "module": "ferro_ta.indicators.statistic", + "doc": "", + "params": [] + }, + { + "name": "DTW_DISTANCE", + "category": "statistic", + "module": "ferro_ta.indicators.statistic", + "doc": "", + "params": [] + }, { "name": "DX", "category": "momentum", @@ -3283,6 +3304,13 @@ "doc": "", "params": [] }, + { + "name": "BATCH_DTW", + "category": "statistic", + "module": "ferro_ta.indicators.statistic", + "doc": "", + "params": [] + }, { "name": "BETA", "category": "statistic", @@ -3297,6 +3325,20 @@ "doc": "", "params": [] }, + { + "name": "DTW", + "category": "statistic", + "module": "ferro_ta.indicators.statistic", + "doc": "", + "params": [] + }, + { + "name": "DTW_DISTANCE", + "category": "statistic", + "module": "ferro_ta.indicators.statistic", + "doc": "", + "params": [] + }, { "name": "LINEARREG", "category": "statistic", @@ -4735,7 +4777,7 @@ ] }, "rust_core": { - "public_function_count": 349, + "public_function_count": 351, "functions": [ { "module": "aggregation", @@ -6187,6 +6229,16 @@ "function": "correl", "file": "statistic.rs" }, + { + "module": "statistic", + "function": "dtw_distance", + "file": "statistic.rs" + }, + { + "module": "statistic", + "function": "dtw_path", + "file": "statistic.rs" + }, { "module": "statistic", "function": "linearreg", @@ -6485,7 +6537,7 @@ ] }, "wasm_node": { - "export_count": 221, + "export_count": 222, "exports": [ "ad", "adosc", @@ -6546,6 +6598,7 @@ "digital_price", "donchian", "drawdown_series", + "dtw_distance", "dx", "early_exercise_premium", "ema", @@ -6712,9 +6765,9 @@ } }, "parity_summary": { - "python_indicator_count": 207, - "wasm_export_count": 221, - "common_python_wasm_count": 91, + "python_indicator_count": 210, + "wasm_export_count": 222, + "common_python_wasm_count": 92, "common_python_wasm": [ "ad", "adosc", @@ -6743,6 +6796,7 @@ "dema", "detect_breaks_cusum", "donchian", + "dtw_distance", "dx", "ema", "ht_dcperiod", @@ -6817,6 +6871,7 @@ "asin", "atan", "batch_apply", + "batch_dtw", "beta", "cdl2crows", "cdl3blackcrows", @@ -6886,6 +6941,7 @@ "cosh", "div", "drawdown", + "dtw", "exp", "feature_matrix", "floor", diff --git a/docs/changelog.rst b/docs/changelog.rst index 5ab8048..e5098e8 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -1,7 +1,7 @@ Release Notes ============= -These docs track package version ``1.1.3``. +These docs track package version ``1.1.4``. 1.1.0-audit (2026-03-28) ------------------------ diff --git a/docs/support_matrix.rst b/docs/support_matrix.rst index 6908699..14020ae 100644 --- a/docs/support_matrix.rst +++ b/docs/support_matrix.rst @@ -180,7 +180,7 @@ For source builds, packaging details, and platform notes, see Release status -------------- -These docs track package version ``1.1.3``. +These docs track package version ``1.1.4``. - Release notes by version: :doc:`changelog` - Canonical project changelog: `CHANGELOG.md `_ diff --git a/pyproject.toml b/pyproject.toml index fc896dc..59084e3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "maturin" [project] name = "ferro-ta" -version = "1.1.3" +version = "1.1.4" description = "Rust-powered Python technical analysis library with a TA-Lib-compatible API" readme = "README.md" license = { text = "MIT" } diff --git a/python/ferro_ta/indicators/statistic.py b/python/ferro_ta/indicators/statistic.py index b52ef09..a441d90 100644 --- a/python/ferro_ta/indicators/statistic.py +++ b/python/ferro_ta/indicators/statistic.py @@ -12,19 +12,33 @@ LINEARREG_ANGLE — Linear Regression Angle (degrees) TSF — Time Series Forecast BETA — Beta CORREL — Pearson's Correlation Coefficient (r) +DTW — Dynamic Time Warping (distance + warping path) +DTW_DISTANCE — Dynamic Time Warping distance only (faster) +BATCH_DTW — Batch DTW: N series vs 1 reference, in parallel """ from __future__ import annotations +from typing import Optional + import numpy as np from numpy.typing import ArrayLike +from ferro_ta._ferro_ta import ( + batch_dtw as _batch_dtw, +) from ferro_ta._ferro_ta import ( beta as _beta, ) from ferro_ta._ferro_ta import ( correl as _correl, ) +from ferro_ta._ferro_ta import ( + dtw as _dtw, +) +from ferro_ta._ferro_ta import ( + dtw_distance as _dtw_distance, +) from ferro_ta._ferro_ta import ( linearreg as _linearreg, ) @@ -247,6 +261,98 @@ def CORREL(real0: ArrayLike, real1: ArrayLike, timeperiod: int = 30) -> np.ndarr _normalize_rust_error(e) +def DTW( + series1: ArrayLike, + series2: ArrayLike, + window: Optional[int] = None, +) -> tuple[float, np.ndarray]: + """Dynamic Time Warping — distance and optimal warping path. + + Parameters + ---------- + series1 : array-like + First time series. + series2 : array-like + Second time series (may differ in length from series1). + window : int, optional + Sakoe-Chiba band width. ``None`` (default) = unconstrained. + + Returns + ------- + distance : float + DTW distance (accumulated Euclidean cost along the optimal path). + path : numpy.ndarray, shape (N, 2) + Warping path as ``(i, j)`` index pairs from ``(0, 0)`` to + ``(len(series1)-1, len(series2)-1)``. + """ + try: + return _dtw(_to_f64(series1), _to_f64(series2), window) + except ValueError as e: + _normalize_rust_error(e) + + +def DTW_DISTANCE( + series1: ArrayLike, + series2: ArrayLike, + window: Optional[int] = None, +) -> float: + """Dynamic Time Warping distance only (faster — no path reconstruction). + + Parameters + ---------- + series1 : array-like + First time series. + series2 : array-like + Second time series (may differ in length from series1). + window : int, optional + Sakoe-Chiba band width. ``None`` (default) = unconstrained. + + Returns + ------- + float + DTW distance (accumulated Euclidean cost along the optimal path). + """ + try: + return _dtw_distance(_to_f64(series1), _to_f64(series2), window) + except ValueError as e: + _normalize_rust_error(e) + + +def BATCH_DTW( + matrix: ArrayLike, + reference: ArrayLike, + window: Optional[int] = None, +) -> np.ndarray: + """Batch Dynamic Time Warping — N series vs 1 reference, computed in parallel. + + Parameters + ---------- + matrix : array-like, shape (N, L) + N time series of length L. Each row is compared against ``reference``. + reference : array-like, shape (L,) + The reference series. + window : int, optional + Sakoe-Chiba band width. ``None`` (default) = unconstrained. + + Returns + ------- + numpy.ndarray, shape (N,) + DTW distance from each row of ``matrix`` to ``reference``. + """ + try: + mat = np.ascontiguousarray(matrix, dtype=np.float64) + if mat.ndim != 2: + from ferro_ta.core.exceptions import FerroTAInputError + + raise FerroTAInputError( + f"matrix must be a 2-D array, got {mat.ndim}-D.", + suggestion="Pass a 2-D NumPy array of shape (N, L).", + ) + return _batch_dtw(mat, _to_f64(reference), window) + except ValueError as e: + _normalize_rust_error(e) + + __all__ = [ "STDDEV", "VAR", @@ -257,4 +363,7 @@ __all__ = [ "TSF", "BETA", "CORREL", + "DTW", + "DTW_DISTANCE", + "BATCH_DTW", ] diff --git a/src/statistic/dtw.rs b/src/statistic/dtw.rs new file mode 100644 index 0000000..30b469b --- /dev/null +++ b/src/statistic/dtw.rs @@ -0,0 +1,98 @@ +use ndarray::Array2; +use numpy::{IntoPyArray, PyArray1, PyArray2, PyReadonlyArray1, PyReadonlyArray2}; +use pyo3::exceptions::PyValueError; +use pyo3::prelude::*; +use rayon::prelude::*; + +/// Dynamic Time Warping — distance and optimal warping path between two 1-D series. +/// +/// Returns a tuple `(distance, path)` where `path` is a NumPy array of shape +/// `(N, 2)` containing `(i, j)` index pairs from `(0, 0)` to `(n-1, m-1)`. +/// +/// Local cost: `|series1[i] - series2[j]|` (Euclidean, matches `dtaidistance`). +#[pyfunction] +#[pyo3(signature = (series1, series2, window = None))] +pub fn dtw<'py>( + py: Python<'py>, + series1: PyReadonlyArray1<'py, f64>, + series2: PyReadonlyArray1<'py, f64>, + window: Option, +) -> PyResult<(f64, Bound<'py, PyArray2>)> { + let s1 = series1.as_slice()?; + let s2 = series2.as_slice()?; + if s1.is_empty() || s2.is_empty() { + return Err(PyValueError::new_err( + "series1 and series2 must not be empty", + )); + } + let (dist, path) = ferro_ta_core::statistic::dtw_path(s1, s2, window); + let n = path.len(); + let flat: Vec = path.iter().flat_map(|&(i, j)| [i, j]).collect(); + let arr = + Array2::from_shape_vec((n, 2), flat).map_err(|e| PyValueError::new_err(e.to_string()))?; + Ok((dist, arr.into_pyarray(py))) +} + +/// Dynamic Time Warping — distance only (faster, no path reconstruction). +/// +/// Returns the accumulated Euclidean cost along the optimal warping path. +/// Use this when you only need the distance, not the alignment path. +#[pyfunction] +#[pyo3(signature = (series1, series2, window = None))] +pub fn dtw_distance<'py>( + _py: Python<'py>, + series1: PyReadonlyArray1<'py, f64>, + series2: PyReadonlyArray1<'py, f64>, + window: Option, +) -> PyResult { + let s1 = series1.as_slice()?; + let s2 = series2.as_slice()?; + if s1.is_empty() || s2.is_empty() { + return Err(PyValueError::new_err( + "series1 and series2 must not be empty", + )); + } + Ok(ferro_ta_core::statistic::dtw_distance(s1, s2, window)) +} + +/// Batch Dynamic Time Warping — compute DTW distance from each row of a 2-D matrix +/// to a single reference series, in parallel. +/// +/// Parameters +/// ---------- +/// matrix : np.ndarray, shape (N, L) +/// N time series of length L. Each row is compared against `reference`. +/// reference : np.ndarray, shape (L,) +/// The reference series. +/// window : int, optional +/// Sakoe-Chiba band width. `None` = unconstrained. +/// +/// Returns +/// ------- +/// np.ndarray, shape (N,) +/// DTW distances, one per row. +#[pyfunction] +#[pyo3(signature = (matrix, reference, window = None))] +pub fn batch_dtw<'py>( + py: Python<'py>, + matrix: PyReadonlyArray2<'py, f64>, + reference: PyReadonlyArray1<'py, f64>, + window: Option, +) -> PyResult>> { + let mat = matrix.as_array(); + let ref_slice = reference.as_slice()?; + + if ref_slice.is_empty() { + return Err(PyValueError::new_err("reference must not be empty")); + } + + let (n_rows, _) = mat.dim(); + let rows: Vec> = (0..n_rows).map(|i| mat.row(i).to_vec()).collect(); + + let result: Vec = rows + .par_iter() + .map(|series| ferro_ta_core::statistic::dtw_distance(series, ref_slice, window)) + .collect(); + + Ok(result.into_pyarray(py)) +} diff --git a/src/statistic/mod.rs b/src/statistic/mod.rs index 7fb2fb4..bcef6e7 100644 --- a/src/statistic/mod.rs +++ b/src/statistic/mod.rs @@ -4,6 +4,7 @@ mod beta; pub(crate) mod common; mod correl; +mod dtw; mod linearreg; mod stddev; mod var; @@ -23,5 +24,8 @@ pub fn register(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(pyo3::wrap_pyfunction!(self::linearreg::tsf, m)?)?; m.add_function(pyo3::wrap_pyfunction!(self::beta::beta, m)?)?; m.add_function(pyo3::wrap_pyfunction!(self::correl::correl, m)?)?; + m.add_function(pyo3::wrap_pyfunction!(self::dtw::dtw, m)?)?; + m.add_function(pyo3::wrap_pyfunction!(self::dtw::dtw_distance, m)?)?; + m.add_function(pyo3::wrap_pyfunction!(self::dtw::batch_dtw, m)?)?; Ok(()) } diff --git a/tests/unit/indicators/test_statistic.py b/tests/unit/indicators/test_statistic.py index 96e2953..8487785 100644 --- a/tests/unit/indicators/test_statistic.py +++ b/tests/unit/indicators/test_statistic.py @@ -1,10 +1,14 @@ """Unit tests for ferro_ta.indicators.statistic""" import numpy as np +import pytest from ferro_ta.indicators.statistic import ( + BATCH_DTW, BETA, CORREL, + DTW, + DTW_DISTANCE, LINEARREG, LINEARREG_ANGLE, LINEARREG_INTERCEPT, @@ -308,3 +312,177 @@ class TestTSF: expected = _naive_linearreg(_A, timeperiod=14, x_value=14.0) result = TSF(_A, timeperiod=14) np.testing.assert_allclose(result, expected, equal_nan=True) + + +# --------------------------------------------------------------------------- +# DTW — Dynamic Time Warping +# --------------------------------------------------------------------------- + +dtai = pytest.importorskip("dtaidistance", reason="dtaidistance not installed") + +_DTW_RNG = np.random.default_rng(42) + + +class TestDTW: + # --- Validation against dtaidistance (SOTA reference) --- + + def test_distance_matches_dtaidistance_random(self): + """Core correctness: our distance == dtaidistance on 20 random pairs.""" + for _ in range(20): + n = int(_DTW_RNG.integers(5, 50)) + a = _DTW_RNG.random(n) + b = _DTW_RNG.random(n) + expected = dtai.dtw.distance(a, b) + actual = DTW_DISTANCE(a, b) + np.testing.assert_allclose( + actual, expected, rtol=1e-9, err_msg=f"Mismatch on series length {n}" + ) + + def test_distance_matches_dtaidistance_unequal_length(self): + """Handles unequal-length series correctly.""" + for _ in range(10): + a = _DTW_RNG.random(int(_DTW_RNG.integers(5, 30))) + b = _DTW_RNG.random(int(_DTW_RNG.integers(5, 30))) + expected = dtai.dtw.distance(a, b) + actual = DTW_DISTANCE(a, b) + np.testing.assert_allclose(actual, expected, rtol=1e-9) + + def test_path_distance_matches_dtaidistance(self): + """DTW() path variant: returned distance matches dtaidistance.""" + a = _DTW_RNG.random(20) + b = _DTW_RNG.random(25) + expected = dtai.dtw.distance(a, b) + dist, _ = DTW(a, b) + np.testing.assert_allclose(dist, expected, rtol=1e-9) + + def test_path_matches_dtaidistance_warping_path(self): + """Warping path matches dtaidistance.dtw.warping_path() on same-length series.""" + for _ in range(10): + n = int(_DTW_RNG.integers(5, 20)) + a = _DTW_RNG.random(n) + b = _DTW_RNG.random(n) + expected_path = dtai.dtw.warping_path(a, b) + _, actual_path = DTW(a, b) + actual_pairs = [tuple(int(x) for x in row) for row in actual_path] + assert actual_pairs == expected_path, ( + f"Path mismatch for n={n}:\n ours={actual_pairs}\n dtai={expected_path}" + ) + + def test_window_constrained_matches_dtaidistance(self): + """Sakoe-Chiba window matches dtaidistance window parameter.""" + a = _DTW_RNG.random(30) + b = _DTW_RNG.random(30) + for w in [3, 8, 15]: + expected = dtai.dtw.distance(a, b, window=w) + actual = DTW_DISTANCE(a, b, window=w) + np.testing.assert_allclose( + actual, expected, rtol=1e-9, err_msg=f"Mismatch at window={w}" + ) + + def test_batch_matches_dtaidistance(self): + """BATCH_DTW matches calling dtaidistance per-row.""" + ref = _DTW_RNG.random(20) + matrix = _DTW_RNG.random((8, 20)) + batch_result = BATCH_DTW(matrix, ref) + for i in range(8): + expected = dtai.dtw.distance(matrix[i], ref) + np.testing.assert_allclose( + batch_result[i], + expected, + rtol=1e-9, + err_msg=f"Batch mismatch at row {i}", + ) + + # --- Mathematical properties --- + + def test_identical_distance_is_zero(self): + a = np.array([1.0, 2.0, 3.0, 4.0, 5.0]) + dist, _ = DTW(a, a) + assert dist == pytest.approx(0.0, abs=1e-10) + + def test_symmetry(self): + a, b = _DTW_RNG.random(20), _DTW_RNG.random(20) + assert DTW_DISTANCE(a, b) == pytest.approx(DTW_DISTANCE(b, a), rel=1e-10) + + def test_triangle_inequality(self): + a, b, c = _DTW_RNG.random(15), _DTW_RNG.random(15), _DTW_RNG.random(15) + assert DTW_DISTANCE(a, c) <= DTW_DISTANCE(a, b) + DTW_DISTANCE(b, c) + 1e-9 + + # --- Known hardcoded values --- + + def test_known_shifted_series(self): + # [0,1,2] vs [1,2,3]: optimal path (0,0)→(1,0)→(2,1)→(2,2) + # Squared costs: 1+0+0+1=2, sqrt(2). Verified against dtaidistance. + a = np.array([0.0, 1.0, 2.0]) + b = np.array([1.0, 2.0, 3.0]) + np.testing.assert_allclose(DTW_DISTANCE(a, b), np.sqrt(2.0), rtol=1e-9) + + def test_known_single_element(self): + # sqrt((3-7)^2) = sqrt(16) = 4.0 + np.testing.assert_allclose( + DTW_DISTANCE(np.array([3.0]), np.array([7.0])), 4.0, rtol=1e-9 + ) + + def test_known_constant_series(self): + assert DTW_DISTANCE(np.full(10, 5.0), np.full(10, 5.0)) == pytest.approx( + 0.0, abs=1e-12 + ) + + # --- Path structural guarantees --- + + def test_path_starts_at_origin(self): + _, path = DTW(_DTW_RNG.random(10), _DTW_RNG.random(10)) + assert tuple(int(x) for x in path[0]) == (0, 0) + + def test_path_ends_at_corner(self): + _, path = DTW(_DTW_RNG.random(7), _DTW_RNG.random(9)) + assert tuple(int(x) for x in path[-1]) == (6, 8) + + def test_path_is_monotone(self): + _, path = DTW(_DTW_RNG.random(20), _DTW_RNG.random(20)) + for k in range(1, len(path)): + assert path[k][0] >= path[k - 1][0] + assert path[k][1] >= path[k - 1][1] + + def test_path_steps_unit_size(self): + _, path = DTW(_DTW_RNG.random(15), _DTW_RNG.random(12)) + for k in range(1, len(path)): + di = int(path[k][0]) - int(path[k - 1][0]) + dj = int(path[k][1]) - int(path[k - 1][1]) + assert di in (0, 1) and dj in (0, 1) + assert not (di == 0 and dj == 0) + + # --- DTW_DISTANCE == DTW distance --- + + def test_distance_only_matches_full(self): + a, b = _DTW_RNG.random(25), _DTW_RNG.random(25) + d_full, _ = DTW(a, b) + np.testing.assert_allclose(DTW_DISTANCE(a, b), d_full, rtol=1e-10) + + # --- Batch --- + + def test_batch_single_row(self): + ref = np.array([1.0, 2.0, 3.0]) + result = BATCH_DTW(np.array([[1.0, 2.0, 3.0]]), ref) + assert result[0] == pytest.approx(0.0, abs=1e-10) + + def test_batch_matches_single_calls(self): + ref = _DTW_RNG.random(20) + matrix = _DTW_RNG.random((8, 20)) + batch = BATCH_DTW(matrix, ref) + for i in range(8): + np.testing.assert_allclose( + batch[i], DTW_DISTANCE(matrix[i], ref), rtol=1e-10 + ) + + # --- Edge cases --- + + def test_empty_series_raises(self): + with pytest.raises((ValueError, Exception)): + DTW(np.array([]), np.array([1.0, 2.0])) + + def test_window_constrained_ge_unconstrained(self): + a, b = _DTW_RNG.random(20), _DTW_RNG.random(20) + d_full = DTW_DISTANCE(a, b) + d_narrow = DTW_DISTANCE(a, b, window=2) + assert d_narrow >= d_full - 1e-9 diff --git a/uv.lock b/uv.lock index b9bea22..ffe567a 100644 --- a/uv.lock +++ b/uv.lock @@ -950,7 +950,7 @@ wheels = [ [[package]] name = "ferro-ta" -version = "1.1.3" +version = "1.1.4" source = { editable = "." } dependencies = [ { name = "numpy" }, diff --git a/wasm/Cargo.lock b/wasm/Cargo.lock index 0bda7ce..3f73689 100644 --- a/wasm/Cargo.lock +++ b/wasm/Cargo.lock @@ -49,11 +49,11 @@ checksum = "9330f8b2ff13f34540b44e946ef35111825727b38d33286ef986142615121801" [[package]] name = "ferro_ta_core" -version = "1.1.3" +version = "1.1.4" [[package]] name = "ferro_ta_wasm" -version = "1.1.3" +version = "1.1.4" dependencies = [ "ferro_ta_core", "js-sys", diff --git a/wasm/Cargo.toml b/wasm/Cargo.toml index 1856108..5730162 100644 --- a/wasm/Cargo.toml +++ b/wasm/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "ferro_ta_wasm" -version = "1.1.3" +version = "1.1.4" edition = "2021" description = "WebAssembly bindings for ferro-ta technical analysis indicators" license = "MIT" diff --git a/wasm/package.json b/wasm/package.json index 49225af..8368d94 100644 --- a/wasm/package.json +++ b/wasm/package.json @@ -1,6 +1,6 @@ { "name": "ferro-ta-wasm", - "version": "1.1.3", + "version": "1.1.4", "description": "WebAssembly bindings for ferro-ta technical analysis indicators", "main": "node/ferro_ta_wasm.js", "module": "web/ferro_ta_wasm.js", diff --git a/wasm/src/lib.rs b/wasm/src/lib.rs index d988a42..c2ebfa9 100644 --- a/wasm/src/lib.rs +++ b/wasm/src/lib.rs @@ -1653,6 +1653,18 @@ pub fn correl(real0: &Float64Array, real1: &Float64Array, timeperiod: usize) -> from_vec(ferro_ta_core::statistic::correl(&to_vec(real0), &to_vec(real1), timeperiod)) } +/// Dynamic Time Warping distance between two series. +/// +/// Returns the accumulated Euclidean cost along the optimal warping path. +/// Pass `window` as `0` for unconstrained (no Sakoe-Chiba band). +#[wasm_bindgen] +pub fn dtw_distance(series1: &Float64Array, series2: &Float64Array, window: usize) -> f64 { + let s1 = to_vec(series1); + let s2 = to_vec(series2); + let w = if window == 0 { None } else { Some(window) }; + ferro_ta_core::statistic::dtw_distance(&s1, &s2, w) +} + // =========================================================================== // Streaming / Stateful API // ===========================================================================