Repository navigation
pyspi 3.0.0 - #103
Open
willedibam wants to merge 123 commits into
Open
pyspi 3.0.0#103willedibam wants to merge 123 commits into
willedibam wants to merge 123 commits into
Conversation
…s, parallel compute, new SPIs (LaggedCorrelation, CrossPairwiseDistance)
Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
- basic.py: fix SpearmanR/KendallTau label duplication (self.labels += [x] mutated shared class attribute; changed to self.labels = self.labels + [x]) - utils.py: stub is_jpype_jvm_available() and check_optional_deps() to always return False/empty — JIDT replaced with pure-numpy equivalents, no JVM needed - infotheory.py: implement AIS-based auto-embedding for TransferEntropy (gaussian and kraskov estimators); replaces JIDT MAX_CORR_AIS path that was returning NaN with 'unhandled estimator' warning - calculator.py: related fixes for numpy-only operation Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
CrossPairwiseDistance._dist is RMSE-by-construction (sqrt(mean((a-b)**2))), so the identifier now reflects that explicitly: xpdist_<metric>_tau-<K>_<stat>_rmse. Mirrors the _rmse suffix used by PairwiseDistance and DynamicTimeWarping when normalise=True. No behavior change.
JIDT/Java replacements have been numpy-pure since 7dca294 but the legacy gates remained: - calculator.compute() excluded kernel/symbolic/auto-embed TE from the fork worker pool under a "JVM does not survive fork()" rationale that no longer applies; reunify the parallel path. - config YAMLs declared dependencies: [java] on eight infotheory SPI classes, causing Calculator to print spurious "install java" warnings and skip 24 SPIs in fast_config.yaml. - utils.is_jpype_jvm_available() is a dead stub; fold check_optional_deps down to an empty dict and drop the "Checking optional dependencies..." print. - __init__.py's OMP_NUM_THREADS=1 comment referenced JPype; rewrite to reflect the real reason (fork-worker BLAS oversubscription). tests/test_SPIs.py previously failed collection when any current SPI key (e.g. dtw_constraint-sakoe-chiba_radius-auto) lacked a frozen-baseline entry. Compute the intersection of current vs baseline keys, skip the symmetric differences with a logged summary, and parametrize only the shared set.
Replace setup.py + requirements.txt with a self-contained pyproject.toml. Bump numpy>=2.0, scipy>=1.11, pandas>=2.0, scikit-learn>=1.3, mne>=1.6, pyEDM>=1.15. Drop sktime (replace with lazy aeon import for .ts files, which is the only place it was used). Drop future (unused). mne migration: port wavelet.py and misc.py from the removed mne.connectivity submodule to the mne-connectivity package. The new API returns Connectivity objects rather than tuples; call .get_data(output='dense') and .freqs on them. PSI now calls phase_slope_index per-band because the new API aggregates internally. numpy 2 fixes: - infotheory: use stdlib math.factorial (np.math was removed). - causal: np.max instead of Python max on 2-D slice (int() on 1-D arrays was permissive in numpy 1.x). pyEDM 2.5 fix: - causal: EmbedDimension now requires a target arg even for self- simplex; pass columns as target. pandas 3 fix: - basic: df.corr().values is read-only; copy before fill_diagonal. - calculator: defensively copy SPI returns before filling diagonal so any SPI returning a view doesn't break the fork worker. cdt still imports torch.utils.data at module load without declaring torch as a Requires-Dist, so torch stays an explicit dependency. pyEDM 2.5 no longer needs pkg_resources at runtime, so the setuptools<80 pin is no longer required.
Replaces the single-dataset z-score correctness test with a 3-dataset
ATOL/RTOL tolerance test (CML7 + new VAR1 + new Kuramoto). Baselines
come from upstream pyspi 2.0.1 run in an ephemeral uv venv (10 trials,
seed 42, element-wise mean per SPI). Drift is logged to a session-end
summary table via the spi_warning_logger fixture rather than failing
the suite — version-bump noise stays visible without blocking CI.
- tests/generate_benchmark_datasets.py: fixed-seed VAR(1) and Kuramoto
generators (M=7, T=100) matching cml7.npy dimensions.
- pyspi/data/{var1_7,kuramoto_7}.npy: frozen synthetic datasets.
- tests/{VAR1,Kuramoto}_benchmark_tables.pkl: upstream-computed
per-SPI mean/std across 10 trials.
- tests/test_SPIs.py: parametrized over (dataset, SPI), passes when
abs_diff <= ATOL (1e-6) OR abs_diff <= RTOL * |ref| (1e-2). SPIs
absent on either side are skipped with a one-line per-dataset report.
- tests/conftest.py: logger now records (max_abs, max_rel) instead of
a single z-score; summary table widened to 90 cols.
842 cases pass. Most drift appears on Kuramoto spectral/wavelet SPIs
(coherent phase data saturates them near 1.0 so any BLAS/RNG noise
crosses the 1% threshold) and on ill-conditioned covariance SPIs
(prec*, coint_johansen*) where baseline and fork disagree by large
factors that are numerical not algorithmic.
Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
CrossPairwiseDistance is a bespoke SPI written for the CML case study. The previous per-lag cost was a plain RMSE over the untrimmed overlap, which did not line up with how DynamicTimeWarping(normalise=True) scores the same pair. The new cost realises one specific monotone warping path — pair x[k+s] with y[k] over the overlap, then stutter the dropped boundary samples against the opposite series' corner — and normalises by sqrt(T). Because the result is a valid DTW path cost and DTW minimises over all paths, the identity dtw_rmse <= xpdist now holds by construction, which is the property the case study relies on. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
Two estimator-correctness fixes to infotheory.py: 1. Gaussian entropy: introduce _gaussian_log_det(cov, ridge_rel=1e-8) with epsilon = ridge_rel * mean(diag(cov)) and NaN fallback on persistent singularity. Refactor GaussianEntropyCalculator and the module-level _gaussian_entropy_from_data helper to use it. Fixes xme_gaussian_k10 overflow on Kuramoto (was -inf from slogdet sign<=0 -> 1.8e+308 via np.nan_to_num in regression test) and reduces 30+ unit drift on cce_gaussian and di_gaussian. Rationale matches JIDT's NOISE_LEVEL_TO_ADD=1e-8: cov of x+noise equals Sigma + sigma**2 * I, so ridging Sigma directly is equivalent in expectation and numerically stabler than perturbing samples. 2. KernelTECalculator: DYN_CORR_EXCL was stored via setProperty but never consumed. The compute now queries kd-tree index lists and filters neighbours by |j - i| > dce when dce > 0, matching JIDT's 2*dce+1-point exclusion window. Samples whose counts go to zero after filtering are dropped from the mean; returns NaN if all samples degenerate. Fast path (dce None or <= 0) is byte-identical to the prior implementation. 839/839 regression tests pass. Kuramoto:xme_gaussian_k10 drift drops from 1.8e+308 to 5.9; cce_gaussian 30+ -> 24.9; di_gaussian 30+ -> 17.6. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
- config.yaml: add LaggedCorrelation (18 variants, pearson/spearman/kendall
x tau/max_tau x squared) and CrossPairwiseDistance (2 variants), which
exist as code in pyspi/statistics/{basic,distance}.py but had no YAML
entries in the fork. The fork's config.yaml is now the global SPI
superset that downstream case configs subset from.
- config.yaml + fast_config.yaml: comment out
InformationGeometricConditionalIndependence. cdt's IGCI.predict_proba
takes no kwargs, so the 'kernel' variant has no computational backend
in the fork (and arguably upstream too).
Upstream-inherited commented-out blocks are retained:
TransferEntropy.kernel+MAX_CORR_AIS (JIDT AIS overload bug) and
PartialCoherence (no backing class). Those are deliberate upstream
exclusions, not fork choices.
Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
Python 3.12 venvs no longer bundle setuptools by default (PEP 632), and setuptools>=81 dropped pkg_resources from the default install. pyEDM's LoadData.py still does `import pkg_resources` at module load, so without this pin a fresh venv fails test collection with ModuleNotFoundError. Mirrors the equivalent pin in mts-spi-study-cluster. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
- pyspi/benchmarked_config.yaml: 292-SPI set (~5s-per-SPI threshold at M=8 T=3200 and M=16 T=800). Mirrors configs/pyspi-v2/benchmarked_config.yaml in mts-spi-study-cluster; kept in-tree here as the fork-local reference. - uv.lock: pin resolved deps so downstream repos (incl. Gadi venv rebuilds) reproduce the exact environment.
Calculator.compute() gains n_jobs, checkpoint_dir, resume, mp_context, and progress kwargs. Defaults reproduce prior serial behaviour bit-for-bit; PYSPI_N_JOBS env var still honoured as a fallback. Drops the fork-only _PARALLEL_CALC global path in favour of a portable spawn-default executor. - pyspi/_parallel.py: ProcessPoolExecutor + multiprocessing.shared_memory for the dataset; SPIs bucketed by _cache_namespace so a within-class cache (covariance, spectral_mv/bv, barycenter, coint) is reused across every variant in the bucket; per-SPI failures yield NaN matrices; atomic per-SPI .npy checkpoints enable transparent resume. - pyspi/__main__.py: thin CLI (python -m pyspi compute --data ... --config ... --n-jobs ... --checkpoint-dir ...). - pyspi/statistics/*: six _cache_namespace tags on the relevant SPI bases. - tests/test_parallel.py: parity at n_jobs=2,3; checkpoint resume; failure isolation; CLI smoke. All 47 test_calc.py and 846 test_SPIs.py regression tests against the frozen benchmark tables (CML7, Kuramoto, VAR1) pass unchanged. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
…ness
Calculator and Data both gain verbose= (default True) so library users can
silence the per-SPI init dump, the dataset detrend/normalise banner, and
the compute() completion summary. Dependency warnings still print
unconditionally — they affect correctness, not noise.
YAML parsing extracted into pyspi.calculator.load_spis_from_yaml(), a
plain function reused by Calculator and the parallel worker initializer
(workers no longer instantiate a throwaway Calculator with redirected
stdout — direct loader call).
Data.__init__ now initialises _data/n_processes/etc. explicitly so the
empty-instance state used by add_process() is no longer load-bearing on
hasattr() checks.
Test suite layout:
tests/test_calc.py -> tests/test_calculator.py
tests/test_SPIs.py -> tests/test_regression.py (marked slow)
tests/smoke_test.py -> tests/test_smoke.py (now auto-collected;
8 tests, was 0)
pyproject.toml registers the 'slow' marker and defaults pytest to
-m 'not slow', so the 12-minute regression sweep is opt-in via -m slow
while every-commit pytest stays at ~40s.
bench/bench_compute.py is the non-notebook replacement for the sister-repo
walltime sweep. Outputs tidy CSV + JSON metadata. Two intended workflows
documented in the module docstring: amortized-config re-cutting (n_jobs=1)
and parallel-speedup characterisation (n_jobs swept).
Net: 71/71 fast tests pass (~38s). Slow regression suite intact (run with
'pytest -m slow').
Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
The CML7/VAR1/Kuramoto baseline tables are pandas DataFrames + numpy arrays only — no closures or lambdas — so stdlib pickle suffices and dill becomes one less optional dep. Verified by loading all three .pkl files via pickle.load() and by collecting all 839 regression tests successfully. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
Logging: - New pyspi/_logging.py: a 'pyspi' logger with a colorama-tinted ColorFormatter. __init__.py attaches a NullHandler (library default: silent unless the app configures a handler). - Calculator(verbose=) maps to a process-global pyspi logger level: True -> INFO (config load, per-SPI init, compute summary), False -> WARNING (problems only). Dependency-exclusion notices now log at WARNING so they surface regardless of verbose. - All informational print()/Fore calls in calculator.py and data.py routed through the logger. SPI-failure reporting stays on warnings.warn (unchanged semantics; tests assert pytest.warns). - CLI gains --quiet. Per-SPI progress: - Parallel runs now tick the tqdm bar once per SPI, not once per cache bucket. Workers post (key, failed) to a manager queue after each SPI; the main loop drains it for progress while still harvesting results (and dead-worker detection) from the futures. A stuck SPI is now visible as a stalled bar. Data.__init__ drops the short-lived verbose kwarg added earlier this branch — logging supersedes it. 71/71 fast tests pass; 839 slow regression tests unaffected. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
mp_context default is now platform-aware: fork on Linux, spawn on macOS/Windows. Measured head-to-head (252-SPI benchmarked90 config, M=16 T=800): fork reaches 2.50x at n_jobs=4 / 3.65x at n_jobs=8, vs spawn's 1.60x / 1.95x. spawn re-imports the heavy pyspi stack and re-instantiates every SPI per worker (~2.5s tax); fork inherits that state via copy-on-write. Bug fix: run_parallel created the SharedMemory block and Manager *outside* the try, so a failure constructing either leaked the shared memory. Both are now created inside the try with None-guarded cleanup. bench/bench_compute.py upgraded to the canonical benchmark suite: - (M, T, n_jobs) grid with presets (headline/scaling/parallel/amortized) - incremental JSON output, --resume, --array-index for PBS arrays - environment block (pyspi git sha, dep versions + fingerprint, platform) - per-cell wall time, per-SPI timings, peak RSS, failed-SPI count - prints a speedup-vs-n_jobs summary Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
run_benchmark.pbs runs bench_compute.py on Gadi (parametrised via -v: PRESET, CONFIG, REPEATS, PYSPI_DIR, VENV). README documents the presets, output schema, and the n_jobs / amortized-config invariance. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
cdt (causal discovery toolbox, transitive dep) autosets SETTINGS.NJOBS to cpu_count() at import. Inside an n_jobs worker pool that nests process x thread parallelism and oversubscribes. _worker_init now pins cdt.SETTINGS.NJOBS=1 after the SPI modules import. Note: this fixes cdt-based causal SPIs but NOT ConvergentCrossMapping, whose internal threading is pyEDM-side. CCM-containing configs (the full config.yaml) should be benchmarked with n_jobs=1; the benchmarked_* configs exclude CCM and parallelise cleanly. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
- _parallel.py: pin BLAS/OpenMP/cdt/torch in workers; _pin_blas_env sets thread env vars before spawn so workers start single-threaded. causal.py gates pyEDM parallelism (parallel/numProcess) on PYSPI_PIN_BACKENDS. - Tag ConvergentCrossMapping with _cache_namespace so its variants share one worker and cache instead of recomputing. - bench/cut_config.py: cut benchmarked<N>[_amortized] configs from a bench_compute JSON (amortized or raw cost) — replaces the notebook step. - bench/run_benchmark.pbs + run_benchmark_physics.pbs: M/T grid override and PBS-array (--array-index) support for config-cutting runs at real data sizes. - test_parallel.py: cover CCM/pinning, parametrise fork+spawn. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
USYD Physics has no walltime cap on the user's queue; size the PBS defaults to whole-node + a week so heavy bench cells (full config at large M,T) fit in one sequential job without overrides. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
bench_compute.py: write one JSON per (M, T, n_jobs) cell named
<label>_M<M>_T<T>_n<n>.json. M/T/n_jobs/environment now live at the top
level of each file (self-contained); --output-dir + --label replace
the old single --output path. --resume skips per-cell files.
cut_config.py: reads one per-cell JSON directly; --m/--t cell-picker
removed (file IS the cell).
run_benchmark[_physics].pbs (generic, env-var driven): rewritten to use
--output-dir + --label and to redirect stdout/stderr to bench/logs/.
bench/physics/run_bench_{light,m32,m64}.pbs: dedicated submission scripts
for the M={4,8,16}/{32}/{64} x T={200,400,800,1600,3200} grid. Each is a
single qsub away.
Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
Add MXX size-applicability tags to every SPI in the bundled configs (plus the auxiliary spi_labels.csv reference table). Calculator now splits per-config labels off the constructor kwargs and merges them into the SPI's label list (handling Mxx overrides cleanly), so that downstream consumers can filter SPIs by intended M-range.
- bench_compute: psutil per-cell RSS (was process high-watermark accumulating across cells); failed_spis list; per-SPI category/labels - cut_config: handle per-config labels via _split_config_params; preserve dropped variants/classes as commented blocks by default - analyse_cells: anchor-stability Jaccard, cumulative-cost elbow, per-SPI log-linear scaling fits, extrapolation to a target cell - forecast_cell: predict cell wall-time at (M, T) from scaling fits - backfill_metadata: retrofit category/labels into older cell JSONs - physics PBS scripts: write into bench/results/cells/ to match the new layout (cells/ for per-cell JSONs, analysis/ for derived data)
Cluster sweep on the USYD Physics queue, full config.yaml at n_jobs=1:
M=4,8 across T={200,400,800,1600,3200}; M=16,32 across T={200,400,800,1600};
plus M=64 at T={200,400}. Backfilled with per-SPI category/labels at
current pyspi sha.
Analysis (bench/results/analysis/) holds the cross-cell artefacts that
feed the cut decisions: long_costs.csv, scaling fits, Jaccard matrices
at p{80,90,95,99}, cumulative-cost plots, kept-set drift heatmaps,
extrapolations to (M=64, T=3200), and a human-readable report.md.
Scaling-fit validation against the 2 unseen M=64 cells: median |log
error| ~14% per SPI; total cell-wall predictions within 13% of
measured. The (M=64, T=3200) extrapolations should be good to ~25-30%.
Re-cut from physics_config_M16_T800_n1.json (pyspi badfb4e) via bench.cut_config --mode amortized. Anchored at M=16,T=800 because kept-set Jaccard across 20 cells (M=4..64, T=200..3200) is >=0.95 at p>=90, so the anchor choice is empirically immaterial in the operating range. Counts: p80 keeps 262/328, p90 295/328, p95 312/328, p99 325/328. Dropped variants are commented out (not deleted) so the YAML carries the cut provenance. p98/p97 are deliberately omitted because the 9 CCM-with-EmbedDim variants share a cache namespace and tie in amortized cost; cuts that fall inside the tied block would partially drop the group with zero time saved. Also includes the previously-untracked non-amortized benchmarked configs (80/85/90) and the 85_amortized variant for completeness.
cut_config.amortized_costs was bucketing by _cache_namespace and
spreading the bucket's total cost evenly over every member. But two
classes have caches keyed on a constructor param within the namespace,
so members of different sub-buckets don't actually share work:
- Barycenter: cache is per (mode, pair) — variants only share within
the same mode (euclidean / dtw / sgddtw / softdtw).
- NonparametricSpectralBivariate (5 directed-spectral classes): each
class's measure call dominates the shared per-pair Connectivity
build, so amortization should be per class.
Both classes now declare a _cache_subkey property; amortized_costs and
analyse_cells group by (namespace, *subkey) so the cost of one bucket
isn't smeared across independent caches. Other cache-grouped classes
(Covariance/Precision, CCM, Cointegration, spectral_mv) have the same
bug shape but it doesn't move any cut at the percentiles in use, so
left as-is to keep the change surgical.
Re-cut benchmarked{80,90,95,99}_amortized at M=16,T=800:
- p95 / p99: unchanged (cuts above the affected SPIs).
- p90: 7 swap in / 7 swap out. Now drops hsic, softdtw, te_kraskov_DCE
(fixed-embedding) — actually expensive at large (M, T) — and keeps
the 12 cheap Barycenter variants (euclidean/dtw/sgddtw modes) that
were wrongly dropped together with softdtw. Predicted wall time at
M=64, T=3200 drops from ~11h to ~4h.
- p80: 12 swap each way. Drops dcoh + gpdcoh classes wholly (they
belong above the threshold per-class); keeps the same 12 cheap bary
variants.
Runtime parallel scheduler unaffected: it buckets by _cache_namespace
for worker placement (correct: cache lives per data instance), the
subkey is purely an amortization concept for cut decisions.
… classes Cache-aware corollary to the previous commit: once the expensive shared computation runs for one member of a cache bucket, the rest of the bucket is free at compute time (they're just cheap post-lookup transforms). Dropping a strict subset wastes information for zero time saving. New rule in cut_config.py: if any member of a cache bucket is kept, snap-promote the whole bucket. Kept count may exceed the percentile target by the size of the promoted partial buckets (header records both the target and the snap delta). To make the snap rule visible across all cache classes, declare _cache_subkey on the remaining ones (the previous commit only did Barycenter + spectral_bv since they were the ones biting in shipped cuts): - Covariance/Precision: (estimator,) - ConvergentCrossMapping: (embedding_dimension,) - Cointegration: (method, det_order, k_ar_diff) or (method, autolag, maxlag, trend) depending on method - NonparametricSpectralMultivariate: (class_name,) - GroupDelay / PhaseSlopeIndex override: (class_name, fmin, fmax) matching their narrower cache key Re-cut at M=16,T=800: only p80 changes (262 -> 267; +5 from snapping DirectedCoherence, which had 1/6 variants kept at the boundary). p90/p95/p99 unchanged (no partial buckets existed there).
Delete five configs that no code path references after the recent cache-aware re-cut: - benchmarked_config.yaml (original benchmarked cut, superseded) - benchmarked80_config.yaml (non-amortized; raw-cost variant - benchmarked85_config.yaml unused) - benchmarked90_config.yaml - benchmarked85_amortized_config.yaml (85% tier dropped; p95 covers it) Verified zero load-time references and that the four bundled subsets (fast/sonnet/fabfour, all wired into Calculator subset=) still load. Updated two example invocations in bench/README.md and bench/run_benchmark.pbs that pointed at benchmarked90_config.yaml to benchmarked90_amortized_config.yaml.
Adds the newly-completed (M=64, T=800) per-cell JSON (cell_wall ~19.6h at the full config, n_jobs=1), backfilled with per-SPI category/labels. Re-ran bench.analyse_cells over the now-21-cell corpus; scaling fits tightened (median |log error| per SPI: ~14% -> ~7%; total cell-wall prediction at M=64,T=800 within +10% of observed). All artefacts in bench/results/analysis/ refreshed accordingly.
…over labels Three defects left by the phase-regression repair. *Units.* `C.frequencies` is in Hz, so the fitted slope is radians per Hz and `slope/(2*pi)` is a delay in **seconds**. A true 4-sample lag came back as 4.0, 2.0 and 1.0 at fs = 1, 2 and 4, while the API, the identifier and every other lagged SPI in pyspi count samples. The delay is scaled by fs; no change at the shipped fs=1. *Fit quality.* `rvalue` stored the signed regression r symmetrically, but the fit is of the phase of C_ij and phase(C_ji) = -phase(C_ij), so reversing the process order turned +0.99997 into -0.99997 at the mirrored position -- for a statistic declared symmetric. It is now |r|, which is orientation-free, and the variant is labelled undirected and unsigned to match. *Labels.* A structural trait now replaces `directed`/`undirected` in the merged label set rather than sitting beside a stale one. `gd_*` carried the class's `antisymmetric` and the config's `directed` simultaneously, so `filter_spis` answered both ways for the same SPI; six SPIs were in that state. `gd_*_delay` values are unchanged at fs=1, which is what every bundled config uses. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…eshold `sigonly` is a pointwise amplitude cut -- keep the lags with `|r(l)| > 1.96/sqrt(T)` -- and the name is a historical misnomer. It is not an inferential test: the band is not inflated for the series' own autocorrelation, and it is not corrected for being applied at every one of the ~T/2 lags in the window. Under the null it therefore keeps something almost surely; no seed in 400 produced an empty surviving set at T=600. Both facts are now stated on the class and at the call site, and the Bartlett attribution is dropped -- citing it for the 1/sqrt(T) scale invited exactly the reading the name already encourages. When nothing does clear the cut, the statistic is now 0. Falling back to the unfiltered window reported the largest of ~T/2 sample correlations under the null -- 0.128 for `max` on the test fixture, 0.191 squared -- which is the opposite of what a threshold is for. Tests exercise the returned SPI values rather than the cached lag profile, and cover the no-survivor case for all four enabled combinations, the sign of the association for positive and negative coupling, and the cut itself reconstructed independently from the profile. The CHANGELOG said eight `xcorr_*` SPIs change; `full` contains six. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Reverts the `_rmse` -> `_norm-rootT` rename for non-Euclidean pairwise metrics. The suffix is pyspi's house label for "divided by sqrt(T)", shared with `xpdist` and `DynamicTimeWarping(normalise=True)`, and renaming it for one of the three introduced a second convention for a quantity no bundled config computes. It is documented instead: literal RMSE for `euclidean`, the normalisation and nothing more for cityblock, cosine, chebyshev, canberra and braycurtis. `CrossPairwiseDistance` is documented as what it is -- the root-T-normalised cost of a lag-shifted alignment path with the boundary samples stuttered against the opposite series' endpoint -- and explicitly *not* a conventional path-length RMSE, since the divisor is the record length while the path has T + t steps. The stutter is what makes the alignment a valid DTW path, so `dtw_rmse <= xpdist` holds by construction; that is now a randomised test over five seeds rather than a claim in a comment. `tau` was validated as `int(tau) < 0`, which accepted 1.7 (truncated to 1) and True (silently 1) and raised an opaque conversion error on nan/inf. Validation moves to `utils.require_int`, now shared with the transfer-entropy embedding parameters so the rule is one rule: integral, not boolean, finite, in range. Tests add the definition written out by hand, symmetry and process permutation, bivariate/multivariate equality, and the tau=0 identity with `pdist_euclidean_rmse`. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
pyspi used exactly four functions from `cdt.causality.pairwise`: the ANM independence score (HSIC of a GP residual against the cause), the conditional distribution similarity statistic, the RECI regression-error score and IGCI. They are transcribed into `pyspi/lib/pairwise_causal.py` on NumPy/SciPy/scikit-learn, with the original algorithm attributions and cdt's MIT licence retained at `pyspi/lib/LICENSE-cdt.txt` and shipped in the wheel. **No value changes.** Verified against cdt 0.6 over 60 mixed random pairs (continuous, random-walk, quantised, independent, various lengths): max |cdt - local| is exactly 0 for all four scores. Five of those pairs are frozen into `tests/data/fixtures/cdt_pairwise_reference.npz` so the evidence outlives the dependency, and the slow baseline suite shows no `anm`/`cds`/`reci`/`igci` drift. The ANM score itself was already pyspi's own function, monkey-patched over cdt's (which regressed the wrong way round -- cdt PR #155); only the HSIC came from cdt. Torch drove no pyspi computation. It entered because cdt eagerly imports its Torch-backed models at package load; `InterDependenceScore` is NumPy, and the comment in `_parallel` saying otherwise is corrected along with the now-dead cdt/torch thread pinning. Measured, rather than assumed: - Installed size 1.1 GB -> 546 MB (torch 370 MB plus triton and cdt). - `import cdt.causality.pairwise` cost 1.5-1.7 s, paid once per process -- which under `compute(n_jobs=N)` was once per worker. - Per-call runtime is unchanged: hsic 0.78 -> 0.75 ms, cds 5.38 -> 5.40 ms, reci 0.387 -> 0.386 ms, igci 1.00 -> 0.89 ms at T=800. The win is install size and start-up, not throughput. Tests are three kinds kept separate: independent formulas (an explicit double sum for HSIC, `lstsq` for RECI, a direct order-statistic sum for the IGCI entropy), invariants and known directions, and the frozen cdt fixture as regression evidence. CDS's directional claim is asserted only where it holds -- it goes the wrong way on a cubic pair in 8 of 8 seeds, being one heuristic feature of the Jarfo model rather than a consistent estimator. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
… policy *Benchmark data seeds no longer depend on `n_jobs`.* The n_jobs sweep exists to time the same problem at different worker counts, and seeding on n_jobs handed each column of that sweep a different dataset -- so a scaling curve compared runs on data that were never the same. `n_jobs` stays in the cell *identity*, because it is what the cell measures. *Resume identity gains the algorithm and loses the repeat count.* It now covers `COMPUTATION_VERSION` and the pyspi git sha, so a cell measured before an output-changing change is not reused after it. `repeats` is dropped from it: reuse is already allowed when the stored run has at least as many repeats as requested, so folding the requested count into the identity would discard a good 10-repeat cell the moment someone asked for 5. *`Data.add_process` could create a duplicate name.* It appended `proc-<index>` unconditionally, so a caller who named their processes `["a", "proc-2"]` and appended twice produced a second `proc-2` -- which `to_frame()` cannot stack, and which the constructor would have rejected. It now falls past any name already in use. *The CLI fails on partial failure by default.* `--fail-on-error` is replaced by `--allow-partial`, inverting the default. The previous argument -- that a few SPIs legitimately fail on real data, so a non-zero exit would be the normal case -- no longer holds for the shipped configs, which now run clean on every bundled fixture and on the default smoke data; and it was the weaker consideration against a pipeline silently ingesting a table with failed columns. The partial table is still written, and the failed identifiers still go to stderr and into the NPZ. *`setuptools>=77`* in build-requires, matching the PEP 639 licence metadata already in use; older setuptools would build a wheel with no licence metadata. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
* `tests/README.md` said drift is reported and not enforced. It is enforced -- tolerance violations fail, an all-NaN baseline is rejected, and two gates require that no SPI raises and none produces an empty column. The per-module tolerance split it described is also gone. * The CLI's `--mp-context` help claimed spawn is the default everywhere; `_default_mp_context` returns fork on Linux, spawn elsewhere. * The nats standardisation note claimed the bits/nats split made "correlating, ranking or thresholding" wrong across estimators. Pearson and Spearman between two columns are invariant to a positive rescaling and were never affected; what was affected is any comparison of magnitudes -- a shared threshold, a difference or ratio, a cross-estimator "which found the most" ranking. * `tools/measure_reproducibility.py` compared only entries finite in both runs, so an SPI whose NaN *pattern* moved would have been scored as reproducing exactly. It compares the masks first now; re-measured, still 325/325. * `test_bivariate_rejects_indices_passed_positionally` was `xfail(strict=True)` against a bug fixed earlier in this branch, and expected the old `TypeError`. It is now a passing test asserting the message that names the signature. `test_state_integrity`'s and `test_structural_traits`' docstrings counted open findings that no longer exist. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Every changed statistic was accounted for before regenerating, not after: * **Renamed, 3 SPIs on all three fixtures.** `te_kraskov_*_k-max-*` and `gc_gaussian_k-max-*` gain the auto-embedding method. Their values change too, because `MAX_CORR_AIS` now selects the source embedding instead of hard-coding it to (1, 1). * **Changed, 4 SPIs on kuramoto_M7_T100 only.** The `xcorr_*_sig-True` variants, max |delta| 0.195. Traced to the mechanism, not assumed: 16 of 42 pairs have no lag clearing `1.96/sqrt(100) = 0.196`, and those now return 0 instead of the unfiltered-window statistic. The old values at exactly those pairs ran 0.109-0.195 -- every one of them below the cut the same code had just applied. * **Nothing else moved.** In particular `gd_*` is unchanged, confirming the sample scaling is the no-op at fs=1 that it was predicted to be, and no `anm`/`cds`/`reci`/`igci` value moved when cdt was removed. `COMPUTATION_VERSION` becomes `3.0.0.r2` -- `<release>.r<revision>`, a counter within a release, since only equality is ever tested. One bump for the final output-changing state, covering the KSG per-occurrence dither, full MAX_CORR_AIS, the group-delay scaling and |r|, and the cross-correlation threshold. Checkpoints and run digests from any earlier state of this branch are now refused rather than resumed. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
… changes Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
**Boundary.** The AIS scorer computed `N_eff = T - (dim-1)*delay`, but `_ksg_ais` aligns a one-step-ahead future against the embedding and so spends a sample on the shift as well as on the lookback: the aligned arrays have `T - 1 - (dim-1)*delay` rows. At T=5, kNN=4, dimension 1 the guard therefore saw N=5 (usable neighbours 4, exactly k) and passed, while the real N is 4 (usable 3). `_ksg_mi_general` had no guard of its own, so `tree.query(..., k=5)` padded with infinities and the invalid candidate scored a finite **0.0**. Three changes: the scorer uses the aligned count, `_ksg_mi_general` validates the arrays it was actually handed, and `_select_embedding` raises when no candidate is scorable instead of returning (1, 1) -- an embedding the scorer had just rejected, handed to the estimator as though it had been selected. No shipped value moves: at T=100 with k_max=10, tau_max=4 the smallest aligned N is 63 against k=4, so no candidate was near the boundary. The slow suite passes against the existing baselines, so no regeneration and no COMPUTATION_VERSION bump -- this is stricter rejection of inputs that were already invalid. **Oracle.** The Gaussian TE test built its source design from lags 1..l with a lookback of `l*l_tau`, one sample older throughout than the alignment `_te_build_embeddings` declares (lags 0..l-1, lookback `(l-1)*l_tau`). It passed at rel=2e-3 because on a smoothly autocorrelated source the two lag sets carry nearly the same predictive information. The oracle is rebuilt to the declared alignment and the fixture changed to a *white* source driving the target at lag 1, where the two designs differ by far more than the tolerance -- the test asserts that gap explicitly, so the fixture cannot quietly stop discriminating. Tolerance tightens from 2e-3 to 1e-6, the residual being the documented 1e-8 covariance ridge. The implementation was correct; the evidence was not. Also adds an O(N^2) multivariate KSG reference and uses it to check `_ksg_ais` candidate scores at four (dimension, delay) settings and the selected embedding against an independently computed argmax. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
`int()` and `float()` accept far more than the parameter contracts do:
`int(2.7)` is 2, `int(True)` is 1, `float(True)` is 1.0, and `float("inf")`
sails through a `<= 0` check. Each turns a plainly wrong argument into a
plausible one, silently.
Integral, non-boolean, at or above their minimum: `prop_k` (>=1), numeric
`dyn_corr_excl` (>=0, with `None` and the exact string `"AUTO"` as the only
alternatives -- `"auto"`, `"AUTO "` and `"10"` are now rejected), `n` on
`CausalEntropy`/`DirectedInfo` (>=1), the TE embedding and search parameters
(already strict, unchanged), `LaggedCorrelation.tau` (>=0),
`sakoe_chiba_radius` (>=1) and `CrossPairwiseDistance.tau` (>=0).
Finite real, non-boolean, strictly positive: `kernel_width` and
`sakoe_chiba_ratio`. These are genuinely continuous, so a fractional value
stays legitimate -- a test asserts that too, so the strictness cannot quietly
break the case it exists to allow.
One shared pair of validators in `utils` (`require_int`,
`require_positive_float`), so the rule is one rule rather than a per-class
habit. The validated value is what reaches the identifier, the cache key and
the computation: `prop_k=np.int64(6)` gives `mi_kraskov_NN-6` and a cache key
of plain `int`s, not a numpy type stringified into the name.
Internal calculated indices and the calculators' own `setProperty` paths are
untouched -- they are not public parameters.
No output changes: this is stricter rejection of inputs that were already
invalid, so no baseline regeneration and no COMPUTATION_VERSION bump.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
**The closed form was wrong.** With `x == y` exactly, every joint L-infinity distance equals the marginal one, so the kth joint radius *is* the kth marginal radius and strict counting gives `n_x = n_y = k - 1`, not k. The degenerate estimate is therefore `psi(N) - psi(k)`, not `psi(k) - 2*psi(k+1) + psi(N)` -- which is off by `1/k` and did not match the numbers printed beside it (4.234 vs a measured 4.7341 at N=400). `psi(N) - psi(4)` reproduces all four measurements to four decimals. Corrected in the code comment, the test and the CHANGELOG, and the test now asserts the estimate is *away* from that branch rather than only near the true entropy. **The scope claim was too broad.** Per-occurrence dither changes only paths containing byte-identical normalised coordinate occurrences -- duplicate or exactly collinear processes, a source equal to its target. It does not change ordinary KSG results over distinct coordinates, and this is why no bundled baseline moved even though appending the replica index changed the dither key for every column: only the integer neighbour counts enter the estimate, and a 1e-8 perturbation does not move which points fall inside the radius. Now a test: the same continuous pair returns the identical 15 significant figures under four different dither seeds. **CDT/Torch references finished off** -- README install note and worker-pinning paragraph, the benchmark dependency fingerprint, the PBS comment, and the IDS docstring's `torch.Tensor` parameter types. The two vendored MIT notices are now in the wheel's `dist-info/licenses` as well as in package data. **Remaining claims corrected without renaming anything:** `/sqrt(T)` does not make arbitrary metrics comparable across record lengths (only the Euclidean norm grows as sqrt(T)); the Bartlett attribution is gone from `sigonly`, which is a historical pointwise amplitude threshold; the cross-correlation `max` docstring says "largest retained positive correlation" rather than "strongest association"; `~328` becomes 325 in three places; `tests/README.md` now describes enforced drift, nats, and the single remaining xfail -- recorded explicitly as a fixture/low-data finding on `var1_M3_T100`, not a proven universal defect. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
**IGCI is Information-Geometric Causal Inference**, not "conditional independence" -- it tests no conditional independence. `InformationGeometric CausalInference` is the accurate name; the old one is retained as a subclass that emits a `DeprecationWarning`, so existing configs and scripts keep working with an identical identifier and identical values. Its score is a difference of two entropies and therefore exactly antisymmetric (`A[i,j] == -A[j,i]`, asserted on a fixture). Reporting it as `Unsigned` was not cosmetic: `Calculator._rmmin` shifts every column it believes unsigned by that column's minimum, which on an antisymmetric matrix moves both orientations equally and destroys the sign, and `set_group` folds the directions together through `abs()`. It is now signed and antisymmetric, and stays commented out of the bundled configs -- this is metadata correctness, not a claim that the heuristic is reliable. **ANM was labelled `linear`** while fitting a Gaussian process and testing independence with an RBF-kernel HSIC. Now `nonlinear`. **Structural traits win over config-declared directedness.** The rule was one-sided -- `antisymmetric`/`asymmetric` displaced `directed`/`undirected` -- so a family-level `directed` could still sit on a statistic that declares itself undirected. `gd_*_rvalue` is symmetric by construction (it stores |r|) and came back both `undirected` and `directed` when loaded through YAML, so `filter_spis` answered either way for the same SPI. The SPI's own trait, taken from its instance labels before the family and per-config labels are merged, now displaces the others in both directions. Signedness continues to follow `issigned()`. No shipped label changes: the 325 SPIs in `full` already had consistent traits, and IGCI is disabled. Tests cover direct construction and YAML loading for all three GroupDelay statistics, plus a sweep asserting no SPI in any of the eight bundled configs carries two directedness traits. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
`run_digest` hashed `run_spec` verbatim, and `run_spec` carries `config` (the name) and `configfile` (the *absolute resolved path*). The same config contents over the same data therefore digested differently in a source checkout, in an installed wheel and in a temporary directory -- so a perfectly reusable checkpoint was refused. A false negative rather than unsafe reuse, but it defeats the point of a content hash. Location-only fields (`config`, `configfile`, `dataset_name`) are excluded from the digest and stay in `run_spec` as provenance. What is hashed is what the run computes: preprocessing flags, process count and names, observation count, the SPI identifier set, the config *contents*, `COMPUTATION_VERSION`, and the dataset bytes in their native dtype. Tests pin both directions -- identical contents at two paths under two names digest the same; changed contents, changed data, permuted process order, changed preprocessing and a changed computation version each digest differently. The data and permutation cases needed a non-degenerate fixture: `arange(60).reshape(3, 20)` z-scores to three identical rows, so permuting them is a no-op and a constant shift washes out. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
A compact suite over nine deterministic M=3 fixtures at T=64 and T=256 from two fixed seeds -- independent and correlated Gaussian, a directionally coupled VAR, nonlinear additive coupling, binary and four-level quantised, a noisy delayed oscillation, heavy-tailed with deterministic outliers, and duplicate plus near-collinear processes. Twenty representative SPIs are routed to the fixtures they can say something about; nothing runs the full config, and only the CCM representative is marked slow. It enforces the contract that a statistic either produces a defensible result or fails explicitly, and it found four things. **A real defect.** KSG per-coordinate rescaling invariance failed on quantised data. `x` and `1000*x` normalise to columns differing by one ulp (1.1e-16), which flipped the dither key, re-drew the whole dither, and on tied data a different dither separates the ties differently: MI moved by 0.0375 nats under a rescaling the estimator is supposed to be invariant to. The key is now taken from the column rounded to 12 decimals. This is the only change here that can move a valid output -- on tied or quantised input -- so COMPUTATION_VERSION becomes 3.0.0.r3. Tie-free values are unaffected (only the integer neighbour counts enter the estimate) and no baseline moved; the frozen fixtures are continuous. **A test of mine that was vacuous.** `test_the_dither_realisation_does_not_move_a_continuous_estimate` patched the module attribute `_KNN_NOISE_SEED`, which is a *default argument* bound when `_knn_condition` was defined -- so it asserted that a function returns the same value when called twice with the same arguments. It now passes the seed explicitly and compares against the brute-force reference. **Two assumptions of mine that were wrong, recorded rather than asserted away.** ANM is *not* scale-invariant: the independence test standardises, but the fit is a GaussianProcessRegressor on raw values with a length scale bounded away from zero, so rescaling moves the residual (0.27 on this fixture under `4x + 3`). That matches cdt 0.6 exactly and is a property of the method; pyspi's default z-scoring collapses it to ~4e-5, and a test now pins both halves. RECI prefers the *wrong* direction on a saturating tanh pair at both seeds -- its guarantee needs the cause near-uniform after min-max scaling -- so direction is asserted for ANM only, and RECI's correct case stays on the cubic pair in test_pairwise_causal. **One documented exception.** `gd_delay` is entirely NaN on the four structureless fixtures, which is correct: group delay is defined only where the coherence is significant. Listed in `EXPECTED_EMPTY`, and the suite fails if a listed entry starts producing values. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
# Conflicts: # README.md
# Conflicts: # README.md
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Major release. Brings
mainin line withwilledibam/pyspiv3branch. Breaking changes — see "Migrating from 2.x" inCHANGELOG.md, which is the full record; below is a summary of changes.Headline changes
jpype.cdtandtorchremoved. The four pairwise causal functions pyspi used are transcribed on NumPy/SciPy/scikit-learn, bit-identical to cdt 0.6 on the reference pairs. Installed size roughly halves.pyproject.tomlonly;setup.pyandrequirements.txtare gone.Breaking API changes
Calculator(subset=..., configfile=...)→Calculator(config=...)(bundled name or path to a YAML).normalise=→zscore=onCalculatorandData.pyspi/configs/and renamed (full,fast,sonnet,fabfour,benchmarked_p*).Calculator.save()/pyspi.load_table()as.npz; pickle and parquet output are dropped.Bundled configs
Selected by name via
Calculator(config=...); the default isfull.fullbenchmarked_p99benchmarked_p95benchmarked_p90benchmarked_p80fastsonnetfabfourThe
benchmarked_p<N>sets keep the fastest N% of SPIs by measured cost (M=16, T=800, single job). They existed in 2.x only as files; they are now reachable by name.Numerical changes
fullgoes from 328 to 322 SPIs (7 removed, 1 added), each with a stated reason.kernelandsymbolicvalues change by a factor of ln 2.dcoh_*,coint_aeg_*,psi_wavelet_*,ccm_E-None_*,gd_*,xcorr_*,kraskov, and others) are tabulated with causes under "SPI set changes".Dependencies
Python >= 3.10, numpy >= 2.0, pandas >= 2.1, pyEDM >= 2.5, spectral-connectivity >= 1.1,< 3. No setuptools pin.
Testing
.npz.Relation to open PRs and issues
CalculatorFramefix; the setuptools pin is no longer needed).E=Noneselects argmax of rho).Reviewing
pyspi/configs/, ~6.5k); testing pipelines (tests/, ~7.3k); env dependencies (uv.lock, ~3.3k); and benchmarking suites (bench/, ~2.1k).CHANGELOG.md, thenpyspi/calculator.pyandpyspi/data.py, thenpyspi/statistics/, thentests/.After merge
v3.0.0and publish the release.