From 7cfd0050e6b49f774aca8a18b0b9658b5badb934 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Wilhelm=20S=C3=B6derkvist=20Vermelin?= Date: Tue, 13 Oct 2026 09:48:46 +0200 Subject: [PATCH 1/7] bench: fair tsfresh comparison in one script tsfresh_compare.py replaces minimal_bench.py, efficient_bench.py, 1brc/benchmark.py and the stale efficient.json. tsfresh now runs exactly the calculators, parameters and attributes tsfx runs for the tier (derived from .tsfx-config.toml, checked by feature count), both 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. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 9 + python_benchmarks/1brc/README.md | 22 ++- python_benchmarks/1brc/benchmark.py | 72 -------- python_benchmarks/efficient.json | 173 ------------------- python_benchmarks/efficient_bench.py | 74 --------- python_benchmarks/minimal_bench.py | 54 ------ python_benchmarks/tsfresh_compare.py | 238 +++++++++++++++++++++++++++ 7 files changed, 264 insertions(+), 378 deletions(-) delete mode 100644 python_benchmarks/1brc/benchmark.py delete mode 100644 python_benchmarks/efficient.json delete mode 100644 python_benchmarks/efficient_bench.py delete mode 100644 python_benchmarks/minimal_bench.py create mode 100644 python_benchmarks/tsfresh_compare.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 64a30cf..ff8de62 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 diff --git a/python_benchmarks/1brc/README.md b/python_benchmarks/1brc/README.md index 0d70445..b12b75c 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, one timed run each): + +| rows | TSFX | TSFresh | +| ---- | ------ | ------- | +| 10M | 2.9 s | 51.5 s | +| 50M | 21.9 s | 194.4 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() From cce5e12ddfc135ef7130acf71edefa12ec4bc17d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Wilhelm=20S=C3=B6derkvist=20Vermelin?= Date: Tue, 13 Oct 2026 09:48:46 +0200 Subject: [PATCH 2/7] perf(wavelet): number_cwt_peaks without quadratic sorting and matching The noise estimate selected the 10th percentile by sorting a len/20 window per ridge line, and every local maximum was compared with every ridge line. Select the two order statistics instead, match by binary search over the lines' sorted end columns (keeping np.argmin's first-line tie-break), and end lines with one partition instead of repeated Vec::remove. Results are bit-identical (snapshot and 11,857 number_cwt_peaks values on 1BRC and edge-case series). On 50M 1BRC rows: number_cwt_peaks 10.9 s -> 2.6 s, Efficient tier 21.6 s -> 13.7 s. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 10 +++++ python_benchmarks/1brc/README.md | 6 +-- src/feature_extractors/wavelet.rs | 69 ++++++++++++++++++++----------- 3 files changed, 57 insertions(+), 28 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ff8de62..5a66021 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -72,6 +72,16 @@ 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. +### 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. + ## 0.5.0 (2026-10-12) ### Changed diff --git a/python_benchmarks/1brc/README.md b/python_benchmarks/1brc/README.md index b12b75c..b82393b 100644 --- a/python_benchmarks/1brc/README.md +++ b/python_benchmarks/1brc/README.md @@ -17,12 +17,12 @@ python python_benchmarks/tsfresh_compare.py --dataset 1brc --rows 10_000_000 `--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, one timed run each): +stations, best of two timed runs): | rows | TSFX | TSFresh | | ---- | ------ | ------- | -| 10M | 2.9 s | 51.5 s | -| 50M | 21.9 s | 194.4 s | +| 10M | 2.3 s | 52.8 s | +| 50M | 13.0 s | 203.6 s | The full billion rows has not been measured. See the docstring of `tsfresh_compare.py` for how the comparison is kept fair. 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). From 9d4062daec2ff9ff64af53bfe26298339a877bc6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Wilhelm=20S=C3=B6derkvist=20Vermelin?= Date: Tue, 13 Oct 2026 09:48:46 +0200 Subject: [PATCH 3/7] fix(ar_coefficient): exact fit with 2k + 1 values, as statsmodels 0.15 statsmodels 0.15 fits an AR(k) model with exactly as many observations as parameters (0.14 raised ZeroDivisionError, which tsfresh turned into NaN, or 0 for coeff = k). CI's fuzz job did not pin statsmodels, so it picked up 0.15 and failed on ar_coefficient([1.0], k=0): tsfresh 1.0, tsfx 0.0. Fit from 2k + 1 values on, and pin statsmodels==0.15.0 in the fuzz job. Co-Authored-By: Claude Opus 5.5 --- .github/workflows/CI.yml | 2 +- CHANGELOG.md | 9 +++++++++ python_tests/test_tsfresh_parity.py | 6 +++++- src/feature_extractors/regression.rs | 28 +++++++++++++++++++++------- 4 files changed, 36 insertions(+), 9 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 4040ca3..ac1dd27 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: diff --git a/CHANGELOG.md b/CHANGELOG.md index 5a66021..35e8027 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -72,6 +72,15 @@ 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 + +- **`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 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/regression.rs b/src/feature_extractors/regression.rs index 2584835..64ffc54 100644 --- a/src/feature_extractors/regression.rs +++ b/src/feature_extractors/regression.rs @@ -254,9 +254,13 @@ 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 @@ -277,10 +281,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; @@ -583,4 +587,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()); + } } From 9ab8370ef54f25e6ee67b212b6f315500bf5522d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Wilhelm=20S=C3=B6derkvist=20Vermelin?= Date: Tue, 13 Oct 2026 10:05:51 +0200 Subject: [PATCH 4/7] perf: faster ADF, permutation entropy and change quantiles Bit-identical (snapshot plus ~12,900 series of 1BRC and edge-case data per feature). Efficient tier on 100M 1BRC rows: 35.3 s -> 25.5 s. - augmented_dickey_fuller: the exact-fit check computed the fit's scale (O(nobs * p) per candidate lag) for every fit. Compare against the cheap bound max|y| + sum_j |b_j| max|X_j| first; the exact scale is only needed when the residual is within twice the tolerance of that bound, so the decision is unchanged. The Jacobi SVD also sums its three dot products in one pass (same order). 16.6 s -> ~10 s. - permutation_entropy: count patterns as base-`dimension` integers (same order as the rank sequences) instead of allocating a Vec per window; dimensions too large for u64 keep the Vec path. 3.7 s -> 2.1 s. - change_quantiles: compute each (ql, qh, isabs) corridor's changes once and aggregate a copy (the median reorders it). 3.0 s -> 2.5 s. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 12 +++++++ src/feature_extractors/change.rs | 40 ++++++++++++++-------- src/feature_extractors/entropy.rs | 46 +++++++++++++++---------- src/feature_extractors/regression.rs | 51 ++++++++++++++++++++-------- 4 files changed, 103 insertions(+), 46 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 35e8027..23ef241 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -90,6 +90,18 @@ All notable changes to TSFX are documented here. 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. ## 0.5.0 (2026-10-12) 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 64ffc54..0478fb9 100644 --- a/src/feature_extractors/regression.rs +++ b/src/feature_extractors/regression.rs @@ -63,6 +63,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 +100,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 { @@ -209,9 +225,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; } From ad216889d9aed5ba6bacca4fd4612ad1a901225a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Wilhelm=20S=C3=B6derkvist=20Vermelin?= Date: Tue, 13 Oct 2026 12:55:12 +0200 Subject: [PATCH 5/7] ci: build free-threaded and PyPy wheels on Linux again With pyo3's abi3 feature, --find-interpreter builds only the abi3 wheel, so 0.5.0 shipped no wheels for free-threaded CPython or PyPy (0.4.0 had cp314t and pp311 on every Linux target) and those installs fell back to the sdist. Name the interpreters instead: python3.11 (abi3), python3.14t and pypy3.11. maturin has their configuration built in, so they need not be installed in the build images (checked locally for both). Co-Authored-By: Claude Opus 5.5 --- .github/workflows/CI.yml | 10 ++++++++-- CHANGELOG.md | 3 +++ Cargo.toml | 4 ++-- 3 files changed, 13 insertions(+), 4 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index ac1dd27..415a399 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -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 23ef241..ab20be2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -74,6 +74,9 @@ All notable changes to TSFX are documented here. ### 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 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" From e43a286f6a8077028d48dddaa6919280e4e2ca91 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Wilhelm=20S=C3=B6derkvist=20Vermelin?= Date: Tue, 13 Oct 2026 13:54:35 +0200 Subject: [PATCH 6/7] test(fuzz): check exact AR designs, skip non-finite ADF designs - _ar_designs left out the exactly determined designs (n = 2k + 1) that statsmodels 0.15 now fits, so their conditioning went unchecked: the thorough profile found [6519, 1, 0, 0, 1] (k = 2, condition number 6e7), where statsmodels is off by 2e-9 of the largest coefficient (tsfx gives the exact 6518 to within 3e-12). - NaN or +-inf inside an ADF regression design can make LAPACK's SVD loop forever in statsmodels' adfuller (autolag=None), which hung the thorough fuzz run. Only a non-finite last value, which enters no design, is still compared; the parity tests check tsfx's NaN for the others. Co-Authored-By: Claude Opus 5.5 --- python_tests/test_fuzz_tsfresh.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) 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 ] From 978f2519a82de0c8556229d2f910e0c8f6a67bc8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Wilhelm=20S=C3=B6derkvist=20Vermelin?= Date: Tue, 13 Oct 2026 13:54:35 +0200 Subject: [PATCH 7/7] perf(regression): solve well-conditioned fits from the Gram matrix The Householder QR in Ols::fit_prefixes sweeps the whole design once per column; for ADF on long series (42 columns x 11k rows, 3.7 MB) that was memory-bandwidth bound and most of ADF's time. Build X^T X and X^T y in one pass over cache-sized row blocks, with column scaling, and use the Cholesky factor, whose leading blocks serve every lag prefix as the QR's did. Not bit-identical: the QR path is still used, unchanged, when the scaled condition number bound reaches 1e3, when the residual (a difference, y'y - |c|^2) has an estimated relative error above 1e-12, or when a fit could be exact. On 1BRC and edge-case series, ADF results agree with the QR path's within 1e-12 relative with the same chosen lags, AR coefficients within 2e-11 of the largest coefficient. ADF alone on 100M 1BRC rows: 12.5 s -> 4.3 s. tsfresh comparison on 1BRC: 39x at 10M rows (1.3 s vs 52.3 s), 27x at 50M (7.1 s vs 192.7 s). Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 12 ++ python_benchmarks/1brc/README.md | 4 +- src/feature_extractors/regression.rs | 177 ++++++++++++++++++++++++++- 3 files changed, 190 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ab20be2..b5d5d39 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -105,6 +105,18 @@ All notable changes to TSFX are documented here. 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) diff --git a/python_benchmarks/1brc/README.md b/python_benchmarks/1brc/README.md index b82393b..e114a79 100644 --- a/python_benchmarks/1brc/README.md +++ b/python_benchmarks/1brc/README.md @@ -21,8 +21,8 @@ stations, best of two timed runs): | rows | TSFX | TSFresh | | ---- | ------ | ------- | -| 10M | 2.3 s | 52.8 s | -| 50M | 13.0 s | 203.6 s | +| 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. diff --git a/src/feature_extractors/regression.rs b/src/feature_extractors/regression.rs index 0478fb9..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 @@ -138,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 { @@ -164,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$. @@ -286,6 +452,9 @@ fn jacobi_svd(mut a: Vec>) -> (Vec>, Vec, Vec>) /// 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 @@ -542,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).