diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 4040ca3..415a399 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -48,7 +48,7 @@ jobs: HYPOTHESIS_PROFILE: ci run: | # bit-exact parity is checked against these versions - uv pip install --system hypothesis tsfresh==0.21.2 numpy==2.4.4 pandas==3.0.2 scipy==1.17.1 + uv pip install --system hypothesis tsfresh==0.21.2 numpy==2.4.4 pandas==3.0.2 scipy==1.17.1 statsmodels==0.15.0 pytest -m fuzz linux: @@ -75,7 +75,12 @@ jobs: CFLAGS_aarch64_unknown_linux_gnu: "-D__ARM_ARCH=8" with: target: ${{ matrix.platform.target }} - args: --release --out dist --find-interpreter + # One abi3 wheel for CPython >= 3.11, plus version-specific wheels + # for the interpreters abi3 does not cover (free-threaded CPython, + # PyPy). --find-interpreter builds only the abi3 wheel. maturin + # knows these interpreters' configuration, so they need not be + # installed in the build image + args: --release --out dist -i python3.11 -i python3.14t -i pypy3.11 sccache: false manylinux: auto - name: Upload wheels @@ -108,7 +113,8 @@ jobs: uses: PyO3/maturin-action@v1 with: target: ${{ matrix.platform.target }} - args: --release --out dist --find-interpreter + # As for manylinux: abi3 plus free-threaded and PyPy wheels + args: --release --out dist -i python3.11 -i python3.14t -i pypy3.11 sccache: false manylinux: musllinux_1_2 - name: Upload wheels diff --git a/CHANGELOG.md b/CHANGELOG.md index 64a30cf..b5d5d39 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -33,6 +33,15 @@ All notable changes to TSFX are documented here. ### Changed +- **Fair TSFresh benchmark.** `python_benchmarks/tsfresh_compare.py` + replaces `minimal_bench.py`, `efficient_bench.py`, `1brc/benchmark.py` and + `efficient.json`. TSFresh now runs exactly the calculators and parameters + TSFX runs for the tier (derived from `.tsfx-config.toml` and checked by + feature count), both libraries get the same pre-sorted in-memory data with + CSV parsing outside the timer, TSFresh's `n_jobs` matches the Polars + thread pool, and both get a warm-up and repeated runs. Previously TSFX's + time included reading the CSV, the stocks series were not sorted by date + for TSFX, and the feature sets differed. - **Many more output columns by default.** The `Efficient` tier now gives about 750 columns per value column (was about 155), 400 of them from `fft_coefficient`. Existing columns keep their names, but the new ones @@ -63,6 +72,52 @@ All notable changes to TSFX are documented here. scaled to the series or fit) rather than bit for bit, since tsfx does not reproduce numpy's FFT, LAPACK and SIMD `log` arithmetic. +### Fixed + +- **Linux wheels for free-threaded CPython 3.14 and PyPy 3.11 again.** 0.5.0 + shipped only the abi3 wheels, which those interpreters cannot use, so + installing there built tsfx from source (with a Rust toolchain). +- **`ar_coefficient` matches statsmodels 0.15.** A series of exactly + 2k + 1 values (as many observations as AR(k) parameters) now gives the + exact fit, as tsfresh does with statsmodels 0.15, instead of NaN (or 0 for + `coeff = k`), which matched statsmodels 0.14. For example, a single value + with `k = 0` gives that value. The tsfresh fuzz tests in CI now pin + `statsmodels==0.15.0`. + +### Performance + +- **`number_cwt_peaks` is much faster on long series,** with the same + results bit for bit: on 50M rows of 1BRC data (series of about 5,600 + values) it takes 2.6 s instead of 10.9 s. Its noise estimate + now selects the 10th percentile of each window instead of sorting it, and + ridge lines are matched by binary search instead of comparing every pair. + The `Efficient` tier on 50M rows of 1BRC data takes 13.0 s instead of + 21.9 s. +- **`augmented_dickey_fuller`, `permutation_entropy` and + `change_quantiles` are faster,** with the same results bit for bit. + Together they make the `Efficient` tier on 100M rows of 1BRC data take + 25.5 s instead of 35.3 s (−28%), the features alone measured as follows: + - `augmented_dickey_fuller` 16.6 → about 10 s: the check for an exact fit + first compares against a cheap upper bound of the fit's scale, and only + computes the exact scale (a pass over every row and regressor) when the + residual is small enough to matter. + - `permutation_entropy` 3.7 → 2.1 s: each window's pattern is counted as + one integer instead of an allocated list. + - `change_quantiles` 3.0 → 2.5 s: each corridor's changes are computed + once and shared by its aggregations. +- **`augmented_dickey_fuller` and `ar_coefficient` solve well-conditioned + regressions from the Gram matrix.** A Householder QR sweeps the whole + design once per regressor and was limited by memory bandwidth on long + series; the Gram matrix needs one pass. On 100M rows of 1BRC data, + `augmented_dickey_fuller` alone takes 4.4 s instead of 12.6 s. **Not bit + for bit:** results agree with the previous ones to within about + 1e-10 (ADF test statistics and p-values within 1e-12 relative on + 1BRC and edge-case series, the same chosen lags, AR coefficients within + 2e-11 of the fit's largest coefficient). Ill-conditioned, + near-exact and exactly fitting regressions still use the QR decomposition, + so their results, including the deviations from tsfresh documented for + exact fits, are unchanged. + ## 0.5.0 (2026-10-12) ### Changed diff --git a/Cargo.toml b/Cargo.toml index 2d3780b..b0e56a4 100755 --- a/Cargo.toml +++ b/Cargo.toml @@ -15,8 +15,8 @@ crate-type = ["lib", "cdylib"] serde = { version = "1", features = ["derive"] } toml = "1.1" # abi3: one wheel per platform for every CPython >= 3.11, including future -# versions (free-threaded builds, which have no stable ABI, still get -# version-specific wheels) +# versions. Free-threaded CPython and PyPy have no stable ABI: on Linux, CI +# builds version-specific wheels for them (see the -i flags in CI.yml) pyo3 = { version = "0.29", features = ["abi3-py311"] } chrono = "0.4.42" csv = "1.3" diff --git a/python_benchmarks/1brc/README.md b/python_benchmarks/1brc/README.md index 0d70445..e114a79 100644 --- a/python_benchmarks/1brc/README.md +++ b/python_benchmarks/1brc/README.md @@ -1,6 +1,6 @@ # 1 Billion Row Benchmark -To do perform the benchmark you must first generate data (if you are in the root +To perform the benchmark you must first generate data (if you are in the root directory): ```bash @@ -9,14 +9,26 @@ python python_benchmarks/1brc/generate_data.py N where `N` is a number (1_000_000_000 for the proper 1 billion row challenge). -Then run benchmark like so: +Then compare TSFX against TSFresh like so: ```bash -python python_benchmarks/1brc/benchmark.py +python python_benchmarks/tsfresh_compare.py --dataset 1brc --rows 10_000_000 ``` -If `tsfx` runs out of memory you can specify how many threads polars should use, -by setting the environment variable: +`--rows` limits how many measurements are read (default 10 million, `0` for +all). For scale, on an 80-core machine with the `efficient` tier (8,859 +stations, best of two timed runs): + +| rows | TSFX | TSFresh | +| ---- | ------ | ------- | +| 10M | 1.3 s | 52.3 s | +| 50M | 7.1 s | 192.7 s | + +The full billion rows has not been measured. See +the docstring of `tsfresh_compare.py` for how the comparison is kept fair. + +Both libraries use as many threads (TSFresh: processes) as Polars does. To +limit them, for example if you run out of memory, set: ```bash export POLARS_MAX_THREADS=16 diff --git a/python_benchmarks/1brc/benchmark.py b/python_benchmarks/1brc/benchmark.py deleted file mode 100644 index b361af6..0000000 --- a/python_benchmarks/1brc/benchmark.py +++ /dev/null @@ -1,72 +0,0 @@ -import json -import time - -import pandas as pd -import polars as pl -import tsfresh -import tsfx -from tsfx import ExtractionSettings, FeatureSetting - - -def tsfresh_efficient() -> float: - """Benchmark TSFRESH efficient features extraction.""" - print("Benchmarking TSFRESH efficient features extraction...") - print("Loading data...") - df = pd.read_csv( - "./test_data/measurements.txt", - sep=";", - usecols=[0, 1], - names=["location", "temperature"], - ) - df = df.dropna(axis=0, how="any") - print("Data loaded") - with open("./python_benchmarks/efficient.json") as f: - settings = json.load(f) - print("Starting TSFRESH efficient features extraction...") - start = time.time() - fdf = tsfresh.extract_features( - df[["location", "temperature"]], - column_id="location", - column_sort="location", - default_fc_parameters=settings, - ) - end = time.time() - print(f"TSFRESH Efficient extraction took {end - start} s") - print(fdf.head()) - return end - start - - -def tsfx_efficient() -> float: - """Benchmark TSFX efficient features extraction.""" - print("Benchmarking TSFX efficient features extraction...") - print("Loading data...") - lf = pl.scan_csv( - "./test_data/measurements.txt", - separator=";", - has_header=False, - ) - lf = lf.with_columns( - pl.col("column_1").alias("location"), - pl.col("column_2").alias("temperature"), - ) - lf = lf.select(pl.col("location"), pl.col("temperature")) - lf = lf.drop_nulls() - print("Data loaded") - opts = ExtractionSettings( - grouping_cols=["location"], - feature_setting=FeatureSetting.Efficient, - value_cols=["temperature"], - ) - print("Starting TSFX efficient features extraction...") - start = time.time() - fdf = tsfx.extract_features(lf, opts) - end = time.time() - print(f"TSFX Efficient extraction took {end - start} s") - print(fdf.sort(pl.col("location")).head()) - return end - start - - -if __name__ == "__main__": - time1 = tsfx_efficient() - time2 = tsfresh_efficient() - print(f"Speed-up ~ {time2 / time1:.1f}x") diff --git a/python_benchmarks/efficient.json b/python_benchmarks/efficient.json deleted file mode 100644 index f65c9cb..0000000 --- a/python_benchmarks/efficient.json +++ /dev/null @@ -1,173 +0,0 @@ -{ - "variance_larger_than_standard_deviation": null, - "has_duplicate_max": null, - "has_duplicate_min": null, - "has_duplicate": null, - "sum_values": null, - "abs_energy": null, - "mean_abs_change": null, - "mean_change": null, - "median": null, - "mean": null, - "length": null, - "standard_deviation": null, - "variation_coefficient": null, - "variance": null, - "skewness": null, - "kurtosis": null, - "root_mean_square": null, - "absolute_sum_of_changes": null, - "longest_strike_below_mean": null, - "longest_strike_above_mean": null, - "count_above_mean": null, - "count_below_mean": null, - "last_location_of_maximum": null, - "first_location_of_maximum": null, - "last_location_of_minimum": null, - "first_location_of_minimum": null, - "percentage_of_reoccurring_values_to_all_values": null, - "percentage_of_reoccurring_datapoints_to_all_datapoints": null, - "sum_of_reoccurring_values": null, - "sum_of_reoccurring_data_points": null, - "ratio_value_number_to_time_series_length": null, - "maximum": null, - "absolute_maximum": null, - "minimum": null, - "time_reversal_asymmetry_statistic": [ - { "lag": 1 }, - { "lag": 2 }, - { "lag": 3 } - ], - "c3": [{ "lag": 1 }, { "lag": 2 }, { "lag": 3 }], - "cid_ce": [{ "normalize": true }, { "normalize": false }], - "symmetry_looking": [ - { "r": 0.0 }, - { "r": 0.05 }, - { "r": 0.1 }, - { "r": 0.15000000000000002 }, - { "r": 0.2 }, - { "r": 0.25 }, - { "r": 0.30000000000000004 }, - { "r": 0.35000000000000003 }, - { "r": 0.4 }, - { "r": 0.45 }, - { "r": 0.5 }, - { "r": 0.55 }, - { "r": 0.6000000000000001 }, - { "r": 0.65 }, - { "r": 0.7000000000000001 }, - { "r": 0.75 }, - { "r": 0.8 }, - { "r": 0.8500000000000001 }, - { "r": 0.9 }, - { "r": 0.9500000000000001 } - ], - "large_standard_deviation": [ - { "r": 0.05 }, - { "r": 0.1 }, - { "r": 0.15000000000000002 }, - { "r": 0.2 }, - { "r": 0.25 }, - { "r": 0.30000000000000004 }, - { "r": 0.35000000000000003 }, - { "r": 0.4 }, - { "r": 0.45 }, - { "r": 0.5 }, - { "r": 0.55 }, - { "r": 0.6000000000000001 }, - { "r": 0.65 }, - { "r": 0.7000000000000001 }, - { "r": 0.75 }, - { "r": 0.8 }, - { "r": 0.8500000000000001 }, - { "r": 0.9 }, - { "r": 0.9500000000000001 } - ], - "quantile": [ - { "q": 0.1 }, - { "q": 0.2 }, - { "q": 0.3 }, - { "q": 0.4 }, - { "q": 0.6 }, - { "q": 0.7 }, - { "q": 0.8 }, - { "q": 0.9 } - ], - "autocorrelation": [ - { "lag": 0 }, - { "lag": 1 }, - { "lag": 2 }, - { "lag": 3 }, - { "lag": 4 }, - { "lag": 5 }, - { "lag": 6 }, - { "lag": 7 }, - { "lag": 8 }, - { "lag": 9 } - ], - "number_peaks": [ - { "n": 1 }, - { "n": 3 }, - { "n": 5 }, - { "n": 10 }, - { "n": 50 } - ], - "index_mass_quantile": [ - { "q": 0.1 }, - { "q": 0.2 }, - { "q": 0.3 }, - { "q": 0.4 }, - { "q": 0.6 }, - { "q": 0.7 }, - { "q": 0.8 }, - { "q": 0.9 } - ], - "range_count": [ - { "min": -1, "max": 1 }, - { "min": -1000000000000.0, "max": 0 }, - { "min": 0, "max": 1000000000000.0 } - ], - "linear_trend": [{ "attr": "intercept" }, { "attr": "slope" }], - "agg_linear_trend": [ - { "attr": "intercept", "chunk_len": 5, "f_agg": "max" }, - { "attr": "intercept", "chunk_len": 5, "f_agg": "min" }, - { "attr": "intercept", "chunk_len": 5, "f_agg": "mean" }, - { "attr": "intercept", "chunk_len": 5, "f_agg": "var" }, - { "attr": "intercept", "chunk_len": 10, "f_agg": "max" }, - { "attr": "intercept", "chunk_len": 10, "f_agg": "min" }, - { "attr": "intercept", "chunk_len": 10, "f_agg": "mean" }, - { "attr": "intercept", "chunk_len": 10, "f_agg": "var" }, - { "attr": "intercept", "chunk_len": 50, "f_agg": "max" }, - { "attr": "intercept", "chunk_len": 50, "f_agg": "min" }, - { "attr": "intercept", "chunk_len": 50, "f_agg": "mean" }, - { "attr": "intercept", "chunk_len": 50, "f_agg": "var" }, - { "attr": "slope", "chunk_len": 5, "f_agg": "max" }, - { "attr": "slope", "chunk_len": 5, "f_agg": "min" }, - { "attr": "slope", "chunk_len": 5, "f_agg": "mean" }, - { "attr": "slope", "chunk_len": 5, "f_agg": "var" }, - { "attr": "slope", "chunk_len": 10, "f_agg": "max" }, - { "attr": "slope", "chunk_len": 10, "f_agg": "min" }, - { "attr": "slope", "chunk_len": 10, "f_agg": "mean" }, - { "attr": "slope", "chunk_len": 10, "f_agg": "var" }, - { "attr": "slope", "chunk_len": 50, "f_agg": "max" }, - { "attr": "slope", "chunk_len": 50, "f_agg": "min" }, - { "attr": "slope", "chunk_len": 50, "f_agg": "mean" }, - { "attr": "slope", "chunk_len": 50, "f_agg": "var" } - ], - "number_crossing_m": [{ "m": 0 }, { "m": -1 }, { "m": 1 }], - "ratio_beyond_r_sigma": [ - { "r": 0.5 }, - { "r": 1 }, - { "r": 1.5 }, - { "r": 2 }, - { "r": 2.5 }, - { "r": 3 }, - { "r": 5 }, - { "r": 6 }, - { "r": 7 }, - { "r": 10 } - ], - "count_above": [{ "t": 0 }], - "count_below": [{ "t": 0 }], - "mean_n_absolute_max": [{ "number_of_maxima": 7 }] -} diff --git a/python_benchmarks/efficient_bench.py b/python_benchmarks/efficient_bench.py deleted file mode 100644 index cf46064..0000000 --- a/python_benchmarks/efficient_bench.py +++ /dev/null @@ -1,74 +0,0 @@ -import json -import time - -import pandas as pd -import polars as pl -import tsfresh -import tsfx -from tsfx import ExtractionSettings, FeatureSetting - - -def tsfresh_efficient() -> float: - """Benchmark TSFRESH efficient features extraction.""" - print("Benchmarking TSFRESH efficient features extraction...") - print("Loading data...") - df = pd.read_csv("./test_data/all_stocks_5yr.csv") - df = df.dropna(axis=0, how="any") - print("Data loaded") - with open("./python_benchmarks/efficient.json") as f: - settings = json.load(f) - print("Starting TSFRESH efficient features extraction...") - start = time.time() - fdf = tsfresh.extract_features( - df[ - [ - "date", - "open", - "high", - "low", - "close", - "volume", - "Name", - ] - ], - column_id="Name", - column_sort="date", - default_fc_parameters=settings, - ) - end = time.time() - print(f"TSFRESH Efficient extraction took {end - start} s") - print(fdf.head()) - return end - start - - -def tsfx_efficient() -> float: - """Benchmark TSFX efficient features extraction.""" - print("Benchmarking TSFX efficient features extraction...") - print("Loading data...") - lf = pl.scan_csv("test_data/all_stocks_5yr.csv") - lf = lf.drop_nulls() - print("Data loaded") - opts = ExtractionSettings( - grouping_cols=["Name"], - feature_setting=FeatureSetting.Efficient, - value_cols=[ - "open", - "high", - "low", - "close", - "volume", - ], - ) - print("Starting TSFX efficient features extraction...") - start = time.time() - fdf = tsfx.extract_features(lf, opts) - end = time.time() - print(f"TSFX Efficient extraction took {end - start} s") - print(fdf.sort(pl.col("Name")).head()) - return end - start - - -if __name__ == "__main__": - time1 = tsfx_efficient() - time2 = tsfresh_efficient() - print(f"Speed-up ~ {time2 / time1:.1f}x") diff --git a/python_benchmarks/minimal_bench.py b/python_benchmarks/minimal_bench.py deleted file mode 100644 index 1b76619..0000000 --- a/python_benchmarks/minimal_bench.py +++ /dev/null @@ -1,54 +0,0 @@ -import time - -import pandas as pd -import polars as pl -import tsfresh -import tsfx -from tsfresh.feature_extraction.settings import MinimalFCParameters -from tsfx import ExtractionSettings, FeatureSetting - - -def tsfresh_minimal() -> None: - """Benchmark TSFRESH minimal features extraction.""" - print("Benchmarking TSFRESH minimal features extraction...") - print("Loading data...") - df = pd.read_csv("./test_data/all_stocks_5yr.csv") - df = df.dropna(axis=0, how="any") - print("Data loaded") - settings = MinimalFCParameters() - print("Starting TSFRESH minimal features extraction...") - start = time.time() - fdf = tsfresh.extract_features( - df, - column_id="Name", - column_sort="date", - default_fc_parameters=settings, - ) - end = time.time() - print(f"TSFRESH Minimal extraction took {end - start} s") - print(fdf.head()) - - -def tsfx_minimal() -> None: - """Benchmark TSFX minimal features extraction.""" - print("Benchmarking TSFX minimal features extraction...") - print("Loading data...") - lf = pl.scan_csv("test_data/all_stocks_5yr.csv") - lf = lf.drop_nulls() - print("Data loaded") - opts = ExtractionSettings( - grouping_cols=["Name"], - feature_setting=FeatureSetting.Minimal, - value_cols=["open", "high", "low", "close", "volume"], - ) - print("Starting TSFX minimal features extraction...") - start = time.time() - fdf = tsfx.extract_features(lf, opts) - end = time.time() - print(f"TSFX Minimal extraction took {end - start} s") - print(fdf.sort(pl.col("Name")).head()) - - -if __name__ == "__main__": - tsfx_minimal() - tsfresh_minimal() diff --git a/python_benchmarks/tsfresh_compare.py b/python_benchmarks/tsfresh_compare.py new file mode 100644 index 0000000..50aad33 --- /dev/null +++ b/python_benchmarks/tsfresh_compare.py @@ -0,0 +1,238 @@ +"""Head-to-head benchmark of TSFX against TSFresh. + +Usage (from anywhere): + + python python_benchmarks/tsfresh_compare.py --tier efficient + python python_benchmarks/tsfresh_compare.py --dataset 1brc --rows 10_000_000 + +The comparison is kept fair by: + +- Same features: TSFresh runs exactly the feature calculators, parameters and + output attributes that TSFX runs for the chosen tier, derived from the same + `.tsfx-config.toml`. The script checks that both produce the same number of + features per value column before timing anything. +- Same input: the data is read, cleaned, cast to Float64 and sorted by + (id, time) once, outside the timed region. Both libraries then get the same + in-memory rows in the same order, so neither pays for CSV parsing or sorting. +- Same hardware budget: TSFresh's `n_jobs` is set to the size of the Polars + thread pool TSFX runs on. Set `POLARS_MAX_THREADS` to limit both. BLAS is + pinned to one thread per process, as TSFresh recommends for multiprocessing. +- Repeated runs: both get an untimed warm-up on a small subset, then + `--repeats` timed runs; the best and median are reported. +""" + +# ruff: noqa: E402 +import os + +# TSFresh parallelises over processes; multithreaded BLAS inside each worker +# oversubscribes the CPU and slows it down. TSFX does not use BLAS. +for _var in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"): + os.environ.setdefault(_var, "1") + +import argparse +import statistics +import time +import tomllib +import warnings +from collections.abc import Callable +from pathlib import Path + +import numpy as np +import pandas as pd +import polars as pl +import tsfresh +from tsfx import ExtractionSettings, FeatureSetting, extract_features + +ROOT = Path(__file__).resolve().parent.parent +CONFIG_PATH = ROOT / ".tsfx-config.toml" +TIERS = { + "minimal": FeatureSetting.Minimal, + "efficient": FeatureSetting.Efficient, + "comprehensive": FeatureSetting.Comprehensive, +} +# TSFX config names that TSFresh spells differently +RENAMED = { + "absolute_energy": "abs_energy", + "mean_absolute_change": "mean_abs_change", + "percentage_of_reoccurring_values_to_all_datapoints": ( + "percentage_of_reoccurring_datapoints_to_all_datapoints" + ), +} +# TSFX only emits these attributes of the linear-trend fits +TREND_ATTRS = ("intercept", "slope") + + +def active_calculators(tier: FeatureSetting) -> set[str]: + """Feature calculators TSFX actually computes for `tier`. + + Read off TSFX's output columns on a tiny frame, so tier membership never + has to be duplicated here. + """ + df = pl.DataFrame({"id": [0] * 50, "x": np.random.default_rng(0).normal(size=50)}) + settings = ExtractionSettings( + grouping_cols=["id"], + value_cols=["x"], + feature_setting=tier, + config_path=str(CONFIG_PATH), + ) + cols = set(extract_features(df.lazy(), settings).columns) - {"id"} + active = set() + for name in [*tomllib.loads(CONFIG_PATH.read_text()), "length"]: + prefix = f"x__{name}" + # Trend features are named `x__linear_trend_slope` etc. + starts = (f"{prefix}__", *(f"{prefix}_{a}" for a in TREND_ATTRS)) + matched = {c for c in cols if c in (name, prefix) or c.startswith(starts)} + if matched: + active.add(name) + cols -= matched + # Parameterless features without a config section (e.g. `sample_entropy`) + active |= {c.removeprefix("x__") for c in cols} + return active + + +def tsfresh_parameters(tier: FeatureSetting) -> dict: + """TSFresh `default_fc_parameters` matching what TSFX computes for `tier`.""" + config = tomllib.loads(CONFIG_PATH.read_text()) + fc_parameters = {} + for name in sorted(active_calculators(tier)): + params = config.get(name, {}).get("parameters") + if name == "linear_trend": + params = [{"attr": a} for a in TREND_ATTRS] + elif name == "agg_linear_trend": + params = [ + {"attr": a, "chunk_len": p["chunk_size"], "f_agg": p["aggregator"]} + for p in params + for a in TREND_ATTRS + ] + elif name == "mean_n_absolute_max": + params = [{"number_of_maxima": p["n"]} for p in params] + elif name == "cwt_coefficients": + params = [{**p, "widths": tuple(p["widths"])} for p in params] + fc_parameters[RENAMED.get(name, name)] = params + return fc_parameters + + +def load_stocks() -> tuple[pl.DataFrame, str, list[str]]: + """Kaggle S&P 500 daily prices: one series per ticker.""" + value_cols = ["open", "high", "low", "close", "volume"] + df = ( + pl.read_csv(ROOT / "test_data/all_stocks_5yr.csv") + .drop_nulls() + .with_columns(pl.col(value_cols).cast(pl.Float64)) + .sort("Name", "date", maintain_order=True) + .select("Name", *value_cols) + ) + return df, "Name", value_cols + + +def load_1brc(rows: int | None) -> tuple[pl.DataFrame, str, list[str]]: + """First `rows` 1BRC measurements: one series per weather station.""" + # No timestamps: each station's series is in file order, which a stable + # sort by station keeps + df = ( + pl.read_csv( + ROOT / "test_data/measurements.txt", + separator=";", + has_header=False, + new_columns=["location", "temperature"], + schema_overrides={"temperature": pl.Float64}, + n_rows=rows, + ) + .drop_nulls() + .sort("location", maintain_order=True) + ) + return df, "location", ["temperature"] + + +def time_runs(name: str, run: Callable[[], object], repeats: int) -> list[float]: + """Wall-clock seconds of `repeats` calls to `run`.""" + times = [] + for i in range(repeats): + start = time.perf_counter() + run() + times.append(time.perf_counter() - start) + print(f" {name} run {i + 1}/{repeats}: {times[-1]:.2f} s") + return times + + +def main() -> None: + """Time both libraries on one dataset and tier.""" + parser = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + parser.add_argument("--dataset", choices=["stocks", "1brc"], default="stocks") + parser.add_argument("--tier", choices=list(TIERS), default="efficient") + parser.add_argument("--repeats", type=int, default=3) + parser.add_argument( + "--rows", + type=int, + default=10_000_000, + help="1brc only: rows to read (0 for all, which is slow)", + ) + args = parser.parse_args() + + tier = TIERS[args.tier] + n_threads = pl.thread_pool_size() + if args.dataset == "stocks": + df, id_col, value_cols = load_stocks() + else: + df, id_col, value_cols = load_1brc(args.rows or None) + pdf = df.to_pandas() + fc_parameters = tsfresh_parameters(tier) + settings = ExtractionSettings( + grouping_cols=[id_col], + value_cols=value_cols, + feature_setting=tier, + config_path=str(CONFIG_PATH), + ) + + def run_tsfx(data: pl.DataFrame = df) -> pl.DataFrame: + return extract_features(data.lazy(), settings) + + def run_tsfresh(data: pd.DataFrame = pdf) -> pd.DataFrame: + with warnings.catch_warnings(): + # Constant or short series trigger numpy/statsmodels warnings + warnings.simplefilter("ignore") + return tsfresh.extract_features( + data, + column_id=id_col, + default_fc_parameters=fc_parameters, + n_jobs=n_threads, + disable_progressbar=True, + ) + + n_series = df[id_col].n_unique() + print( + f"{args.dataset}: {df.height:,} rows, {n_series:,} series x " + f"{len(value_cols)} value columns, tier={args.tier}, " + f"{n_threads} threads/processes" + ) + + # Warm-up (thread/process pools, lazy imports) and feature-parity check + first_ids = df[id_col].unique(maintain_order=True)[:5] + head = df.filter(pl.col(id_col).is_in(first_ids.implode())) + tsfx_head, tsfresh_head = run_tsfx(head), run_tsfresh(head.to_pandas()) + for col in value_cols: + # TSFX emits one shared `length` column, TSFresh one per value column + n_tsfx = sum(c.startswith(f"{col}__") for c in tsfx_head.columns) + n_tsfx += "length" in tsfx_head.columns + n_tsfresh = sum(c.startswith(f"{col}__") for c in tsfresh_head.columns) + if n_tsfx != n_tsfresh: + raise SystemExit( + f"Feature mismatch for {col!r}: TSFX {n_tsfx} vs TSFresh {n_tsfresh}" + ) + print(f"Both compute {n_tsfresh} features per value column\n") + + print("TSFX") + tsfx_times = time_runs("TSFX", run_tsfx, args.repeats) + print("TSFresh") + tsfresh_times = time_runs("TSFresh", run_tsfresh, args.repeats) + + print(f"\n{'':8}{'best':>10}{'median':>10}") + for name, times in (("TSFX", tsfx_times), ("TSFresh", tsfresh_times)): + print(f"{name:8}{min(times):>9.2f}s{statistics.median(times):>9.2f}s") + best = min(tsfresh_times) / min(tsfx_times) + median = statistics.median(tsfresh_times) / statistics.median(tsfx_times) + print(f"Speed-up: {best:.1f}x (best), {median:.1f}x (median)") + + +if __name__ == "__main__": + main() diff --git a/python_tests/test_fuzz_tsfresh.py b/python_tests/test_fuzz_tsfresh.py index 9e9e562..9958e78 100644 --- a/python_tests/test_fuzz_tsfresh.py +++ b/python_tests/test_fuzz_tsfresh.py @@ -369,6 +369,12 @@ def _deviation(feature: str, x: list[float]) -> bool: return True return not _well_posed(x, _ar_designs(x)) if feature == "augmented_dickey_fuller": + # NaN or ±inf inside a regression design can make LAPACK's SVD loop + # forever in statsmodels (autolag=None); tsfx gives NaN, which the + # parity tests check. A non-finite last value enters no design and is + # compared + if not all(math.isfinite(v) for v in x[:-1]): + return True return not _well_posed(x, _adf_design(x)) or _adf_exact_fit(x) return False @@ -392,13 +398,14 @@ def _well_posed(values: list[float], designs: list) -> bool: def _ar_designs(values: list[float]) -> list: """AutoReg's designs [1, x_{t-1}, ..., x_{t-k}] for the parity config's - k that can be fitted (n >= 2k + 2).""" + k that can be fitted (n >= 2k + 1: statsmodels 0.15 also fits the + exactly determined case).""" x = np.array(values, dtype=float) n = len(x) return [ np.column_stack([np.ones(n - k)] + [x[k - c : n - c] for c in range(1, k + 1)]) for k in (1, 2, 4) - if n >= 2 * k + 2 + if n >= 2 * k + 1 ] diff --git a/python_tests/test_tsfresh_parity.py b/python_tests/test_tsfresh_parity.py index 4c2a00c..7b6ddff 100644 --- a/python_tests/test_tsfresh_parity.py +++ b/python_tests/test_tsfresh_parity.py @@ -1635,8 +1635,12 @@ def _ar_2() -> list[float]: # too short to fit: tsfresh's fallback [nan] * k ([1, 2, 3, 7], 0, 2, math.nan), ([1, 2, 3, 7], 2, 2, 0), + # exactly 2k + 1 values: an exact fit (statsmodels 0.15; 0.14 raised + # where the design had full rank) ([1, 2, 3, 4, 5], 2, 2, 0), - ([1], 0, 0, 0), + ([1], 0, 0, 1), + ([1, 3, 2], 0, 1, 3.5), + ([1, 3, 2], 1, 1, -0.5), ], ) def test_ar_coefficient(values, coeff, k, expected): diff --git a/src/feature_extractors/change.rs b/src/feature_extractors/change.rs index 8ee5f8f..3f3520f 100644 --- a/src/feature_extractors/change.rs +++ b/src/feature_extractors/change.rs @@ -234,28 +234,40 @@ impl Kernel for ChangeQuantiles { .copied() .filter(|v| !v.is_nan()) .collect::>(); - let mut changes = Vec::new(); + // The corridor's changes, shared by the aggregations of the same + // (ql, qh, isabs); the default grid has 3 per corridor and side + let mut corridors: Vec<((f64, f64, bool), Vec)> = Vec::new(); + let mut scratch = Vec::new(); for &(ql, qh, isabs, f_agg) in &self.params { if ql >= qh { out.push(0.0); continue; } - let bins = [ql, qh].map(|q| pandas_quantile(&finite_sorted, q)); - let inside = qcut_first_bin(x, bins); - changes.clear(); - // Only changes that start and end inside the corridor - changes.extend( - (0..x.len().saturating_sub(1)) - .filter(|&i| inside[i] && inside[i + 1]) - .map(|i| { - let d = x[i + 1] - x[i]; - if isabs { d.abs() } else { d } - }), - ); + let key = (ql, qh, isabs); + let changes = match corridors.iter().position(|(k, _)| *k == key) { + Some(i) => &corridors[i].1, + None => { + let bins = [ql, qh].map(|q| pandas_quantile(&finite_sorted, q)); + let inside = qcut_first_bin(x, bins); + // Only changes that start and end inside the corridor + let changes = (0..x.len().saturating_sub(1)) + .filter(|&i| inside[i] && inside[i + 1]) + .map(|i| { + let d = x[i + 1] - x[i]; + if isabs { d.abs() } else { d } + }) + .collect(); + corridors.push((key, changes)); + &corridors.last().unwrap().1 + } + }; out.push(if changes.is_empty() { 0.0 } else { - f_agg.aggregate(&mut changes) + // A fresh copy: the median reorders it + scratch.clear(); + scratch.extend_from_slice(changes); + f_agg.aggregate(&mut scratch) }); } } diff --git a/src/feature_extractors/entropy.rs b/src/feature_extractors/entropy.rs index 900bedd..6be0757 100644 --- a/src/feature_extractors/entropy.rs +++ b/src/feature_extractors/entropy.rs @@ -312,31 +312,43 @@ fn permutation_entropy(x: &[f64], tau: usize, dimension: usize) -> f64 { let less = |a: f64, b: f64| a < b || (b.is_nan() && !a.is_nan()); let equal = |a: f64, b: f64| a == b || (a.is_nan() && b.is_nan()); let windows = (x.len() - dimension) / tau + 1; - // np.unique(axis=0) counts the patterns in lexicographic order - let mut counts: BTreeMap, usize> = BTreeMap::new(); - for k in 0..windows { + // The rank of each value of window k + let ranks = |k: usize| { let window = &x[k * tau..k * tau + dimension]; - let pattern = window - .iter() - .enumerate() - .map(|(i, &v)| { - window - .iter() - .enumerate() - .filter(|&(j, &w)| less(w, v) || (j < i && equal(w, v))) - .count() - }) - .collect(); + window.iter().enumerate().map(move |(i, &v)| { + window + .iter() + .enumerate() + .filter(|&(j, &w)| less(w, v) || (j < i && equal(w, v))) + .count() + }) + }; + let terms = match (dimension as u64).checked_pow(dimension as u32) { + // A pattern as one base-`dimension` number, which orders the + // patterns as their rank sequences: no allocation per window + Some(_) => pattern_terms( + windows, + (0..windows).map(|k| ranks(k).fold(0u64, |code, r| code * dimension as u64 + r as u64)), + ), + None => pattern_terms(windows, (0..windows).map(|k| ranks(k).collect::>())), + }; + -numpy_sum(terms.len(), |i| terms[i]) +} + +/// $p \ln p$ of each pattern's frequency, in increasing pattern order (as +/// `np.unique(axis=0)` counts them, lexicographically). +fn pattern_terms(windows: usize, patterns: impl Iterator) -> Vec { + let mut counts: BTreeMap = BTreeMap::new(); + for pattern in patterns { *counts.entry(pattern).or_insert(0) += 1; } - let terms: Vec = counts + counts .values() .map(|&count| { let p = count as f64 / windows as f64; p * p.ln() }) - .collect(); - -numpy_sum(terms.len(), |i| terms[i]) + .collect() } /// Approximate entropy feature. diff --git a/src/feature_extractors/regression.rs b/src/feature_extractors/regression.rs index 2584835..ae7c832 100644 --- a/src/feature_extractors/regression.rs +++ b/src/feature_extractors/regression.rs @@ -19,6 +19,22 @@ struct Ols { /// an exact fit (see [`Ols::fit`]). const PERFECT_FIT_TOL: f64 = 1e-10; +/// Fits whose column-scaled design has a condition number bound below this +/// are solved from the Gram matrix (see [`Ols::fit_prefixes_gram`]). +const MAX_GRAM_CONDITION: f64 = 1e3; + +/// The largest estimated relative error of a residual sum of squares the +/// Gram matrix path accepts. That residual is computed as a difference, +/// $y^T y - \lVert c \rVert^2$, whose rounding error is about +/// $\varepsilon\,\kappa^2\,y^T y$ for a (scaled) condition number +/// $\kappa$; t-values and information criteria magnify it when few degrees +/// of freedom remain. +const MAX_GRAM_RESIDUAL_ERROR: f64 = 1e-12; + +/// Rows per block when accumulating the Gram matrix, so that a block of every +/// column stays in cache. +const GRAM_BLOCK_ROWS: usize = 512; + /// Information criteria closer than this (relative to their magnitude, at /// least 1) count as equal when choosing the ADF lag. const IC_TIE_TOL: f64 = 1e-9; @@ -48,12 +64,19 @@ impl Ols { /// residual sum of squares is set to 0: the rounding noise would /// otherwise decide t-values and information criteria. /// + /// Well-conditioned fits that are far from exact take a faster path, see + /// [`Ols::fit_prefixes_gram`]; its results agree with this one's to + /// within rounding, not bit for bit. + /// /// `None` if the design is not finite (numpy's SVD raises /// `LinAlgError`); a non-finite `y` propagates into the results. fn fit_prefixes(x: &[Vec], y: &[f64], min_cols: usize) -> Option> { if x.iter().flatten().any(|v| !v.is_finite()) { return None; } + if let Some(fits) = Self::fit_prefixes_gram(x, y, min_cols) { + return Some(fits); + } let (nobs, k) = (y.len(), x.len()); let (r, c) = householder_qr(x, y); // tails[p] = |c[p..]|^2, the part of y no combination of the first @@ -63,6 +86,10 @@ impl Ols { for p in (0..k).rev() { tails[p] = tails[p + 1] + if p < nobs { c[p] * c[p] } else { 0.0 }; } + // For a cheap upper bound of each fit's scale (see below) + let abs_max = |v: &[f64]| v.iter().fold(0.0, |m: f64, e| m.max(e.abs())); + let y_max = abs_max(y); + let col_max: Vec = x.iter().map(|col| abs_max(col)).collect(); let fits = (min_cols..=k) .map(|p| { // The p x p triangular factor of the first p columns @@ -96,18 +123,30 @@ impl Ols { let fitted: f64 = (i..p).map(|j| r_p[j][i] * params[j]).sum(); ssr += (c[i] - fitted) * (c[i] - fitted); } - let scale = (0..nobs) - .map(|row| { - y[row].abs() - + x[..p] - .iter() - .zip(¶ms) - .map(|(col, b)| (col[row] * b).abs()) - .sum::() - }) - .fold(0.0, f64::max); - if (ssr / nobs as f64).sqrt() <= PERFECT_FIT_TOL * scale { - ssr = 0.0; + let rms = (ssr / nobs as f64).sqrt(); + // The scale is at most max |y| + sum_j |b_j| max |X_j| (twice + // that covers rounding), so the O(nobs p) scale is only + // needed when the residual is below the tolerance of that + let bound = y_max + + params + .iter() + .zip(&col_max) + .map(|(b, m)| b.abs() * m) + .sum::(); + if !bound.is_finite() || rms <= 2.0 * PERFECT_FIT_TOL * bound { + let scale = (0..nobs) + .map(|row| { + y[row].abs() + + x[..p] + .iter() + .zip(¶ms) + .map(|(col, b)| (col[row] * b).abs()) + .sum::() + }) + .fold(0.0, f64::max); + if rms <= PERFECT_FIT_TOL * scale { + ssr = 0.0; + } } let tol = s_max * p as f64 * f64::EPSILON; Self { @@ -122,6 +161,134 @@ impl Ols { Some(fits) } + /// [`Ols::fit_prefixes`] from the Gram matrix $X^T X$ and $X^T y$, built + /// in one pass over the rows, and its Cholesky factor $R$ ($R^T R = X^T + /// X$, the same $R$ as a QR decomposition up to signs). A Householder QR + /// sweeps the whole design once per column, which on long series is + /// limited by memory bandwidth; this reads it once. + /// + /// The columns are scaled to unit norm first. Squaring the design squares + /// its condition number, so this path is only taken when every fit is + /// well conditioned (scaled condition number bound below + /// [`MAX_GRAM_CONDITION`]: results within about $10^{6}\,\varepsilon$ + /// relative of the QR path's), leaves a residual large enough that its + /// estimated relative error stays below [`MAX_GRAM_RESIDUAL_ERROR`] and is + /// clearly not an exact fit. `None` otherwise, for the QR path. + fn fit_prefixes_gram(x: &[Vec], y: &[f64], min_cols: usize) -> Option> { + let (nobs, k) = (y.len(), x.len()); + if k == 0 { + return None; + } + let mut gram = vec![vec![0.0; k]; k]; + let mut xty = vec![0.0; k]; + let mut yty = 0.0; + for start in (0..nobs).step_by(GRAM_BLOCK_ROWS) { + let rows = start..(start + GRAM_BLOCK_ROWS).min(nobs); + let y_block = &y[rows.clone()]; + for a in 0..k { + let xa = &x[a][rows.clone()]; + for b in a..k { + gram[a][b] += dot(xa, &x[b][rows.clone()]); + } + xty[a] += dot(xa, y_block); + } + yty += dot(y_block, y_block); + } + if !yty.is_finite() || xty.iter().any(|v| !v.is_finite()) { + return None; + } + // Unit columns: gram[a][b] / (norm_a norm_b), so the diagonal is 1 + let norms: Vec = (0..k).map(|a| gram[a][a].sqrt()).collect(); + if norms.contains(&0.0) { + return None; + } + // Upper triangular Cholesky factor, by rows: r[i][j] for j >= i + let mut r = vec![vec![0.0; k]; k]; + for j in 0..k { + for i in 0..=j { + let g = gram[i][j] / (norms[i] * norms[j]); + let s: f64 = (0..i).map(|l| r[l][i] * r[l][j]).sum(); + if i == j { + let pivot = g - s; + if pivot.is_nan() || pivot <= 0.0 { + return None; + } + r[j][j] = pivot.sqrt(); + } else { + r[i][j] = (g - s) / r[i][i]; + } + } + } + // R^{-1} (upper triangular) by back substitution, column by column; + // |R|_F^2 = k for unit columns, so the condition bound is + // sqrt(k) |R^{-1}|_F + let mut inv = vec![vec![0.0; k]; k]; + let mut inv_sq = 0.0; + for j in 0..k { + inv[j][j] = 1.0 / r[j][j]; + for i in (0..j).rev() { + let s: f64 = (i + 1..=j).map(|l| r[i][l] * inv[l][j]).sum(); + inv[i][j] = -s / r[i][i]; + } + inv_sq += (0..=j).map(|i| inv[i][j] * inv[i][j]).sum::(); + } + let bound_sq = k as f64 * inv_sq; + if bound_sq.is_nan() || bound_sq >= MAX_GRAM_CONDITION * MAX_GRAM_CONDITION { + return None; + } + // c = R^{-T} (scaled X^T y): the first p entries are Q^T y for the + // first p columns + let mut c = vec![0.0; k]; + for j in 0..k { + let s: f64 = (0..j).map(|i| r[i][j] * c[i]).sum(); + c[j] = (xty[j] / norms[j] - s) / r[j][j]; + } + let y_max = y.iter().fold(0.0, |m: f64, e| m.max(e.abs())); + let col_max: Vec = x + .iter() + .map(|col| col.iter().fold(0.0, |m: f64, e| m.max(e.abs()))) + .collect(); + let mut reached = (0..min_cols).map(|i| c[i] * c[i]).sum::(); + let mut fits = Vec::with_capacity(k + 1 - min_cols); + for p in min_cols..=k { + if p > min_cols { + reached += c[p - 1] * c[p - 1]; + } + let ssr = yty - reached; + // bound_sq bounds the squared condition number + if ssr.is_nan() || ssr * MAX_GRAM_RESIDUAL_ERROR < f64::EPSILON * bound_sq * yty { + return None; + } + // beta = R_p^{-1} c, diag((X^T X)^{-1}) = squared rows of R_p^{-1}, + // both scaled back by the column norms + let params: Vec = (0..p) + .map(|i| (i..p).map(|j| inv[i][j] * c[j]).sum::() / norms[i]) + .collect(); + let normalized_cov_diag = (0..p) + .map(|i| (i..p).map(|j| inv[i][j] * inv[i][j]).sum::() / (norms[i] * norms[i])) + .collect(); + // Could this be an exact fit (see fit_prefixes)? Then the QR path + // decides + let bound = y_max + + params + .iter() + .zip(&col_max) + .map(|(b, m)| b.abs() * m) + .sum::(); + if (ssr / nobs as f64).sqrt() <= 2.0 * PERFECT_FIT_TOL * bound || !bound.is_finite() { + return None; + } + fits.push(Self { + params, + ssr, + rank: p, + nobs, + normalized_cov_diag, + }); + } + Some(fits) + } + /// The t-value of parameter `i`: the estimate over its standard error, /// with the error variance estimated as `ssr / (nobs - rank)`. fn tvalue(&self, i: usize) -> f64 { @@ -148,6 +315,21 @@ impl Ols { } } +/// Dot product with eight independent accumulators, so it vectorises (the +/// sum is rounded differently from a sequential one). +fn dot(a: &[f64], b: &[f64]) -> f64 { + let (a_lanes, a_rest) = a.as_chunks::<8>(); + let (b_lanes, b_rest) = b.as_chunks::<8>(); + let mut acc = [0.0; 8]; + for (u, v) in a_lanes.iter().zip(b_lanes) { + for l in 0..8 { + acc[l] += u[l] * v[l]; + } + } + let rest: f64 = a_rest.iter().zip(b_rest).map(|(u, v)| u * v).sum(); + acc.iter().sum::() + rest +} + /// Householder QR decomposition $X = QR$ of the matrix with columns `x` /// (at least as many rows as columns). Returns the columns of the upper /// triangular $R$ (each of length $k$) and $Q^T y$. @@ -209,9 +391,14 @@ fn jacobi_svd(mut a: Vec>) -> (Vec>, Vec, Vec>) let mut rotated = false; for p in 0..k { for q in p + 1..k { - let alpha = dot(&a[p], &a[p]); - let beta = dot(&a[q], &a[q]); - let gamma = dot(&a[p], &a[q]); + // The three dot products in one pass, each summed in the + // same order as `dot` + let (mut alpha, mut beta, mut gamma) = (-0.0f64, -0.0f64, -0.0f64); + for (x, y) in a[p].iter().zip(&a[q]) { + alpha += x * x; + beta += y * y; + gamma += x * y; + } if gamma == 0.0 || gamma.abs() <= f64::EPSILON * (alpha * beta).sqrt() { continue; } @@ -254,13 +441,20 @@ fn jacobi_svd(mut a: Vec>) -> (Vec>, Vec, Vec>) /// # Edge cases /// - Nulls are dropped first; a group with no non-null values gives NaN. /// - `coeff` $> k$ gives NaN. -/// - Fewer than $2k + 2$ values (statsmodels needs more observations than -/// regressors and raises otherwise): NaN for `coeff` $< k$ and -/// 0 for `coeff` $= k$, as tsfresh (whose fallback has only $k$ entries). +/// - Fewer than $2k + 1$ values (statsmodels needs at least as many +/// observations as regressors and raises otherwise): NaN for `coeff` +/// $< k$ and 0 for `coeff` $= k$, as tsfresh (whose fallback has only $k$ +/// entries). +/// - Exactly $2k + 1$ values: the exact fit (no residual degrees of +/// freedom), as statsmodels 0.15; for example the value itself for a +/// single value and $k = 0$. /// - A constant or otherwise collinear series gives the minimum-norm least /// squares solution, as statsmodels' pseudoinverse. For near-collinear /// (ill-conditioned) lag matrices the result is rounding-dependent, in /// statsmodels too. +/// - Well-conditioned fits far from exact are solved from the Gram matrix +/// rather than a QR decomposition (faster on long series), agreeing to +/// within about $10^{-10}$ of the largest coefficient, not bit for bit. /// - **Deliberate deviation from tsfresh:** a series containing NaN or /// $\pm\infty$ gives NaN for every coefficient. tsfresh raises /// (statsmodels' `MissingDataError` is not caught), failing the whole @@ -277,10 +471,10 @@ pub struct ArCoefficient { /// model cannot be fitted. fn ar_fit(x: &[f64], k: usize) -> Option> { let n = x.len(); - // statsmodels needs more observations (n - k) than regressors (k + 1): - // ValueError if fewer, ZeroDivisionError (a degree-of-freedom correction) - // if as many - if n < 2 * k + 2 { + // statsmodels needs at least as many observations (n - k) as regressors + // (k + 1) and raises ValueError otherwise; with exactly as many, 0.15 + // returns the exact fit (0.14 raised ZeroDivisionError) + if n < 2 * k + 1 { return None; } let nobs = n - k; @@ -517,10 +711,16 @@ fn adfuller(x: &[f64], autolag: AdfAutolag) -> [f64; 3] { /// - Nonzero values so small that their squares underflow (below about /// $10^{-150}$) leave the regression numerically rank-deficient; the /// result then depends on rounding, in statsmodels as here. +/// - Well-conditioned regressions that are far from an exact fit (most +/// series) are solved from the Gram matrix instead of a QR decomposition, +/// which is much faster on long series. The two agree to within about +/// $10^{-10}$ relative, not bit for bit; ill-conditioned, near-exact and +/// exactly fitting regressions always use the QR decomposition. /// /// # Cost /// $O(n\,p^2)$ per distinct `autolag`, with $p$ the largest lag (68 for -/// $n = 10^5$): the $p + 1$ candidate fits share one QR decomposition. +/// $n = 10^5$): the $p + 1$ candidate fits share one Gram matrix (or QR +/// decomposition), which needs a single pass over the series. /// /// # tsfresh /// `feature_calculators.augmented_dickey_fuller` (v0.21.2). @@ -583,4 +783,14 @@ mod tests { } assert!(ar_fit(&[1.0, 2.0, 3.0, 7.0], 3).is_none()); } + + #[test] + fn test_ar_fit_exactly_determined() { + // statsmodels 0.15: AutoReg([2.5], lags=0).fit().params + assert_eq!(ar_fit(&[2.5], 0), Some(vec![2.5])); + // 2k + 1 values fit exactly; one fewer cannot be fitted + let params = ar_fit(&[1.0, 3.0, 2.0], 1).unwrap(); + assert!((params[0] - 3.5).abs() < 1e-12 && (params[1] + 0.5).abs() < 1e-12); + assert!(ar_fit(&[1.0, 3.0], 1).is_none()); + } } diff --git a/src/feature_extractors/wavelet.rs b/src/feature_extractors/wavelet.rs index 2d8403e..80e9790 100644 --- a/src/feature_extractors/wavelet.rs +++ b/src/feature_extractors/wavelet.rs @@ -196,22 +196,27 @@ fn ricker(points: usize, a: usize) -> Vec { } /// `scipy.stats.scoreatpercentile(values, per)` (`interpolation_method = -/// "fraction"`); NaN sorts last, as in `np.sort`. -fn score_at_percentile(values: &[f64], per: f64) -> f64 { - let mut sorted = values.to_vec(); - sorted.sort_by(|a, b| match (a.is_nan(), b.is_nan()) { +/// "fraction"`); NaN sorts last, as in `np.sort`. Selects the two order +/// statistics it needs instead of sorting; `buf` is scratch space. +fn score_at_percentile(values: &[f64], per: f64, buf: &mut Vec) -> f64 { + let cmp = |a: &f64, b: &f64| match (a.is_nan(), b.is_nan()) { (true, true) => std::cmp::Ordering::Equal, (true, false) => std::cmp::Ordering::Greater, (false, true) => std::cmp::Ordering::Less, _ => a.partial_cmp(b).unwrap(), - }); - let idx = per / 100.0 * (sorted.len() - 1) as f64; + }; + buf.clear(); + buf.extend_from_slice(values); + let idx = per / 100.0 * (buf.len() - 1) as f64; let i = idx as usize; + let (_, &mut lo, above) = buf.select_nth_unstable_by(i, cmp); if i as f64 == idx { - sorted[i] * 1.0 / 1.0 + lo * 1.0 / 1.0 } else { + // The next order statistic is the smallest value above position i + let hi = *above.iter().min_by(|a, b| cmp(a, b)).unwrap(); let (w0, w1) = ((i + 1) as f64 - idx, idx - i as f64); - (sorted[i] * w0 + sorted[i + 1] * w1) / (w0 + w1) + (lo * w0 + hi * w1) / (w0 + w1) } } @@ -268,18 +273,27 @@ fn count_cwt_peaks(x: &[f64], n: usize) -> usize { for line in ridge_lines.iter_mut() { line.gap += 1; } - let prev_cols: Vec = ridge_lines + // The lines' last columns at the start of the row, sorted, keeping + // the first line at each column (np.argmin picks the first closest) + let mut ends: Vec<(usize, usize)> = ridge_lines .iter() - .map(|l| *l.cols.last().unwrap()) + .enumerate() + .map(|(i, l)| (*l.cols.last().unwrap(), i)) .collect(); + ends.sort_unstable(); + ends.dedup_by_key(|e| e.0); for &col in &max_cols[row] { - // np.argmin: the first closest line - let closest = prev_cols - .iter() - .enumerate() - .min_by_key(|&(i, &c)| (c.abs_diff(col), i)) - .filter(|&(_, &c)| c.abs_diff(col) as f64 <= max_distance) - .map(|(i, _)| i); + // The closest line is at the nearest end column on either side; + // on a tie in distance, the first line wins + let p = ends.partition_point(|&(c, _)| c < col); + let left = p.checked_sub(1).map(|j| ends[j]); + let right = ends.get(p).copied(); + let closest = [left, right] + .into_iter() + .flatten() + .min_by_key(|&(c, i)| (c.abs_diff(col), i)) + .filter(|&(c, _)| c.abs_diff(col) as f64 <= max_distance) + .map(|(_, i)| i); match closest { Some(i) => { let line = &mut ridge_lines[i]; @@ -294,11 +308,12 @@ fn count_cwt_peaks(x: &[f64], n: usize) -> usize { }), } } - for i in (0..ridge_lines.len()).rev() { - if ridge_lines[i].gap > gap_thresh { - final_lines.push(ridge_lines.remove(i)); - } - } + // End lines that skipped too many rows, keeping the others in order + let (ended, alive): (Vec<_>, Vec<_>) = std::mem::take(&mut ridge_lines) + .into_iter() + .partition(|line| line.gap > gap_thresh); + final_lines.extend(ended); + ridge_lines = alive; } // _filter_ridge_lines(cwt, lines, min_length = ceil(n / 4), @@ -307,10 +322,11 @@ fn count_cwt_peaks(x: &[f64], n: usize) -> usize { let window_size = (len as f64 / 20.0).ceil() as usize; let (hf_window, odd) = (window_size / 2, window_size % 2); let row_one = &cwt[0]; - let noise = |ind: usize| { + let mut buf = Vec::new(); + let mut noise = |ind: usize| { let start = ind.saturating_sub(hf_window); let end = (ind + hf_window + odd).min(len); - score_at_percentile(&row_one[start..end], 10.0) + score_at_percentile(&row_one[start..end], 10.0, &mut buf) }; final_lines .iter() @@ -371,7 +387,10 @@ fn count_cwt_peaks(x: &[f64], n: usize) -> usize { /// Irregular series are not affected. /// /// # Cost -/// $O(\text{len} \cdot n^2)$. +/// $O(\text{len} \cdot n^2)$ for the transform and $O(n \cdot \text{len} +/// \log \text{len})$ to join ridge lines, plus a selection over a window of +/// $\lceil \text{len} / 20 \rceil$ for each ridge line's noise: up to +/// $O(\text{len}^2)$ in all, which dominates on long series. /// /// # tsfresh /// `feature_calculators.number_cwt_peaks` (v0.21.2).