Metadata-Version: 2.4
Name: polars-online
Version: 0.4.1
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: Intended Audience :: Financial and Insurance Industry
Classifier: Operating System :: MacOS
Classifier: Operating System :: Microsoft :: Windows
Classifier: Operating System :: POSIX :: Linux
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Programming Language :: Rust
Classifier: Topic :: Scientific/Engineering :: Mathematics
Classifier: Typing :: Typed
Requires-Dist: polars>=1.34.0,<2
Requires-Dist: numpy>=1.24 ; extra == 'numpy'
Provides-Extra: numpy
License-File: LICENSE
Summary: Streaming / online regression models for Polars
Keywords: polars,online-learning,streaming,regression,time-series,kalman,ewma
Author: Hans Gilde
License-Expression: Apache-2.0
Requires-Python: >=3.12
Description-Content-Type: text/markdown; charset=UTF-8; variant=GFM
Project-URL: Changelog, https://github.com/hgilde/polars-online/blob/main/CHANGELOG.md
Project-URL: Documentation, https://hgilde.github.io/polars-online/
Project-URL: Issues, https://github.com/hgilde/polars-online/issues
Project-URL: Repository, https://github.com/hgilde/polars-online

# polars-online

Online model fitting for [Polars](https://pola.rs) — linear models, streaming
moments, clustering and regime detection. A bank of models that learns one
chunk at a time and predicts every row *before* it learns from it. Rust core,
Python API, and a standalone CLI.

It is built for data that does not fit in memory. Feed it a stream and memory
stays at *state + one chunk* however long the stream runs, and the numbers are
the same whether the stream arrives as one chunk or a thousand.

> **A note on Polars versions.** Two of the three interfaces this rides on
> carry no stability promise from Polars, so `polars>=1.34.0,<2` is measured
> rather than guaranteed. A weekly canary runs the whole suite on the newest
> py-polars, and the response to a red one is decided in advance. Details in
> [Versioning and the Polars pin](#versioning-and-the-polars-pin).

**Contents.** [What you get](#what-you-get) ·
[Install](#install) ·
[Quick start](#quick-start) ·
[How a bank sees a stream](#how-a-bank-sees-a-stream) ·
[Running a bank](#running-a-bank) ·
[Memory](#memory-which-calls-stream) ·
[Saving, loading and serving](#saving-loading-and-serving) ·
[Preparing a stream](#preparing-a-stream) ·
[Reading the fit](#reading-the-fit) ·
[Diagnostics, selection and evaluation](#diagnostics-selection-and-evaluation) ·
[Models](#models) ·
[Parallelism](#parallelism) ·
[Performance](#performance) ·
[Against scikit-learn](#against-scikit-learn) ·
[What this is not](#what-this-is-not) ·
[Versioning and the Polars pin](#versioning-and-the-polars-pin) ·
[Testing](#testing) ·
[Development](#development)

## What you get

**Twenty model families, one set of stream semantics.** A spec's clock, decay,
grouping and warm-up mean the same thing whichever model it names.

| model | what it is |
|---|---|
| [`ewridge`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.ewridge) · [math](#ewridge--ew-ridge-on-sufficient-statistics) | exponentially weighted ridge on sufficient statistics — the workhorse; grids over ridge values, feature sets and halflives come almost free |
| [`rls`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.rls) · [math](#rls--recursive-least-squares) | recursive least squares, in the numerically safe square-root form |
| [`lasso`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.lasso) · [math](#lasso--lasso-path-with-free-λ-selection) | lasso / elastic-net path with online λ selection |
| [`kalman`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.kalman) · [math](#kalman--random-walk-β-dynamic-linear-model) | Kalman filter with random-walk coefficients |
| [`huber`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.huber) · [`quantile`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.quantile) · [math](#huber--quantile--robust-regression) | robust and quantile regression |
| [`sgd`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.sgd) · [math](#sgd--stochastic-gradient-descent) | stochastic gradient descent with squared, Huber, quantile, ε-insensitive, Poisson and logistic losses |
| [`pa`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.pa) · [math](#pa--passive-aggressive-regression) | passive-aggressive regression — no learning rate |
| [`ftrl`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.ftrl) · [math](#ftrl--online-logistic-regression) | FTRL-proximal logistic regression, L1-sparse |
| [`ew_cov`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.ew_cov) · [math](#ew_cov--exponentially-weighted-moments) | running mean, variance, covariance, correlation, partial correlation, Mahalanobis distance and PCA |
| [`holt`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.holt) · [math](#holt--holts-linear-trend) | Holt's linear trend — the no-feature baseline |
| [`kmeans`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.kmeans) · [math](#kmeans--exponentially-weighted-k-means) | exponentially weighted k-means — out-of-sample cluster labels, with a split–merge move that finds a cluster born after seeding |
| [`micro`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.micro) · [math](#micro--density-based-clustering-any-shape) | density-based clustering — DenStream micro-clusters linked into clusters of any shape and number; flags the rows that belong to none |
| [`ew_class`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.ew_class) · [math](#ew_class--gaussian-classification-on-ew_cov-moments) | Gaussian classification — QDA, LDA or naive Bayes on one `ew_cov` state per class; a label column in, out-of-sample posteriors out |
| [`seqtest`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.seqtest) · [math](#seqtest--a-sequential-test-of-a-sign-by-betting) | a sequential test of a sign by betting — an e-process you can read at any row; on its own a column's sign, with `a`/`b` whether one spec of the bank predicts closer than another |
| [`marginal`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.marginal) · [math](#marginal--every-pairs-moments-kept-in-the-state) | every (feature, target) pair's running mean, variance, covariance, correlation, slope and t — O(p·T) per row for a wide set of columns, kept in the state and read back as a frame. Optionally the pair at a set of lags, for a `t` that knows about serial dependence, and binned target moments, for the relations a correlation cannot see |
| [`deco`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.deco) · [math](#deco--one-correlation-for-the-whole-matrix) | one correlation for the whole matrix — Engle & Kelly's equicorrelation, or one per block and per pair of blocks, in O(m) a row |
| [`rcov`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.rcov) · [math](#rcov--a-blocks-realised-covariance-robust-to-noise) | a block's realised covariance, robust to microstructure noise — the Barndorff-Nielsen–Hansen–Lunde–Shephard kernel or Christensen–Kinnebrock–Podolskij pre-averaging, emitted when a group closes |
| [`hmm`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.hmm) · [math](#hmm--which-regime-are-we-in) | a Gaussian hidden Markov model, filtered online — `ew_class` without the labels, with a transition matrix that can be learned |
| [`corrchange`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.corrchange) · [math](#corrchange--has-the-correlation-structure-changed) | has the correlation structure changed — the Wied–Krämer–Dehling constancy test span by span, or the size of a change between two windows against a permutation null |
| [`bocpd`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.bocpd) · [math](#bocpd--how-long-has-this-regime-lasted) | how long has this regime lasted — Adams & MacKay's run-length posterior, so the answer is the age of the regime and not a flag |

**Three ways to run a bank, same numbers from each.** A Python loop over
chunks ([`ModelBank`](https://hgilde.github.io/polars-online/polars_online.html#polars_online.ModelBank)); a Polars query ([`lf.online.fit_predict(specs)`](https://hgilde.github.io/polars-online/namespaces.html#polars_online._frame.LazyFrameOnlineNamespace.fit_predict) is a
`LazyFrame` you `collect`, `sink` or batch like any other); or a file-to-file
job (`po.run(...)` from Python, or the `online` CLI from a TOML with no Python
at all). There is also an expression form for a frame in memory — it cannot
stream, and it says so.

**Time, built in.** A clock column in any unit, decay by halflife in that
unit, a ceiling on gaps, session boundaries, a policy for a clock that runs
backwards. One state per group, row weights, warm-up thresholds. Or no clock
at all: then row order is the clock — and with decay off, the bank is plain
least squares over everything it has seen, in any row order.

**The clock does not have to be a time.** It is any monotone numeric column,
so sort a frame by one of its own features and name that feature as the
clock. The halflife is then a bandwidth in that feature's units, and each row
is fit on the rows before it under weight `0.5 ** (Δx / halflife)`. That is
an exponential kernel, so the fit is a *local* regression in that column —
LOESS in one streaming pass, out of a state that does not grow. Same for the
moments: clock `ew_cov` on a column and its correlations are local in it.
[A clock that is not time](#a-clock-that-is-not-time) has the example.

**Two guarantees.** Predictions are out-of-sample by construction. Output is
chunk-invariant. Both are tests, not intentions.

**Tested the way those guarantees demand.** About 650 Rust tests and 2,200
pytest cases: numpy oracles to ~1e-13, cross-checks against river,
hypothesis-generated adversarial streams, chunk and thread invariance, golden
numbers on every OS, the README's own code blocks executed. CI runs all of it
on macOS, Windows and Linux on every push, and a weekly canary runs it on the
newest Polars. Details under [Testing](#testing).

**State is a file.** Save a bank; load it to keep learning, or to serve
predictions without learning. The file is written atomically, is the same
bytes from every entry point, and carries nothing host-specific.

**Introspection and diagnostics.** Coefficients as a table or as columns.
Residual sigma and z-scores, drift detection, online model selection and
averaging across a grid, streaming R²/IC/hit rate, residual quantiles and
autocorrelation — all out-of-sample, all O(state).

**Parallel by group, deterministic by construction.** Every (spec, group)
pair is one task on a thread pool, the halflives of a grid run side by side,
and the runner overlaps reading, fitting and writing. Rows within a stream go
one at a time — that is the recursion — so the thread count changes the speed
and nothing else, and a test holds it to that. Details under
[Parallelism](#parallelism).

## Install

```sh
pip install polars-online      # or: uv add polars-online
```

Wheels for macOS (arm64, x86_64), Windows x64 and Linux (x64 glibc and musl,
aarch64 glibc) are on PyPI; those and the CLI binaries are attached to each
GitHub release. Python 3.12+.
The wheel is ~19 MB to download and ~59 MB installed: it statically links the
Rust half of Polars, so nothing beyond `polars` itself has to be present at
run time. `numpy` is an optional extra: only `ModelBank.gram()` and the
[`po.gram`](https://hgilde.github.io/polars-online/gram.html#module-polars_online.gram), [`po.corr`](https://hgilde.github.io/polars-online/corr.html#module-polars_online.corr) and `po.sim` toolkits need it.

From source, with [uv](https://docs.astral.sh/uv/) and a stable Rust toolchain:

```sh
uv sync
uv run maturin develop --release -m crates/online-py/Cargo.toml
```

## Quick start

```python
import polars as pl
import polars_online as po

spec = po.spec.ewridge(
    "ridge",
    targets=["ret"], features=["signal_a", "signal_b"],
    clock="ts", halflife=600.0, max_dclock=300.0,   # decay in seconds; gaps capped at 5 min
    group="bond_id",                                 # one model per bond
)

# Fit and predict over a stream; save the fitted state when the run reaches the last row
(
    pl.scan_parquet("ticks/*.parquet")
    .online.fit_predict([spec], save_state="bank.state")
    .filter(pl.col("ridge").struct.field("n_eff") > 100)
    .sink_parquet("fitted.parquet")                  # memory: state + one chunk, end to end
)

# Later: score new rows against the saved state, learning nothing
scored = pl.scan_parquet("today.parquet").online.predict("bank.state").collect()

# The result is one struct column per spec; unnest it into plain columns
flat = scored.online.unnest([spec])   # pred_ret, resid_ret, n_eff, coef_ret_intercept, coef_ret_signal_a, ...
```

Each spec adds one struct column, named after the spec, holding
`pred_<target>`, `resid_<target>`, `n_eff`, `coef` and whatever diagnostics
you switch on. `df.online.fit_predict(specs)` does the same for a frame in
memory, eagerly.

## How a bank sees a stream

These parameters are shared by every model.

### Time and decay

| parameter | meaning |
|---|---|
| `clock` | a monotone **numeric** column — seconds, cumulative volume, anything. `None` means row count. Temporal columns are rejected; cast first (`pl.col("ts").dt.epoch("s")`), so that the units of `halflife`, `max_dclock` and `session_gap` are yours, not the dtype's |
| `halflife` / `lam` | decay in clock units, one or the other. A list of halflives gives one accumulator per value. `halflife=inf` (or `lam=1.0`) turns decay off |
| `max_dclock` | ceiling on the clock step; required with a `clock`. `0` disables decay, `inf` removes the ceiling |
| `on_clock_reset` | what a backwards clock means: `"max"` (default), `"zero"`, `"reset_state"`, or `"error"` to refuse the chunk and leave the bank as it was |
| `session`, `session_gap` | on a session change, apply this clock delta; `"reset"` resets the state, `inf` never applies it |
| `session_shrink`, `long_halflife` | `ewridge` only: at a session change, mix partway back toward a slow-moving twin |

Per-row decay is `λ = 0.5 ** (Δclock / halflife)`.

### A clock that is not time

The clock is any monotone numeric column, and nothing says the column has to
be a time. Sort the frame by one of its features and name that feature as the
clock. `halflife` is then a bandwidth in that feature's units, and every row
is fit on the rows before it under weight `0.5 ** (Δx / halflife)` — an
exponential kernel. The fit is a local linear regression, computed in one
pass with `O(k²)` of state rather than a window of rows.

```python
curve = df.sort("x0")
local = po.spec.ewridge(
    "local",
    targets=["y"],
    features=["x0"],
    clock="x0",                 # the clock is a feature: the decay is a kernel in it
    halflife=0.5,               # the bandwidth, in x0's own units
    max_dclock=1.0,             # a wider gap decays as if it were this wide
    max_rows_between_solves=1,  # refit at every row
    min_periods=10.0,
)
fitted = po.ModelBank([local]).fit_predict(curve).unnest("local")
```

`pred_y` is the fitted curve read at each row's own `x0`, and `coef` is the
line through that row's neighbourhood. These are the batch numbers: against a
kernel-weighted least squares recomputed from scratch at every row they agree
to 1e-12 (`tests/test_ewridge.py`, `test_a_feature_as_the_clock_is_a_local_linear_regression`).

The kernel is one-sided, because a row is scored before it is learned from.
On `sin(x)` with a bandwidth of 0.25 the fit sits 0.08 from the truth where
the best straight line sits 0.39, and it follows the curve with the lag that
looking only backwards implies. The bandwidth trades that lag against noise:
at 1.0 the same stream gives 0.29, most of the way back to the line.

Three things carry over unchanged. `features` need not include the clock
column — with other features the same fit is a varying-coefficient model,
whose coefficients move along the clock. `group` gives one local fit per key
from the same pass. And the clock belongs to every model, not to `ewridge`:
clocked on a feature, `ew_cov` reports moments and correlations local in it,
and `marginal` does the same pair by pair.

### Groups, weights and warm-up

| parameter | meaning |
|---|---|
| `targets`, `features` | column names, ≥1 target; targets share one `X'X`. Columns must be numeric (any width, `Decimal` and `Boolean` included; cast to `f64`) — a String column is refused rather than cast to null. Columns the spec does not name pass through untouched |
| `add_intercept` | default `True` |
| `group` | one state per key |
| `group_close` | `"monotone"` or `"session"`: emit the group's accumulators when it is finished and drop the stream, which is what bounds a bank over an unbounded key space. Read them with `bank.closed_groups()` |
| `weight` | row weight column |
| `min_periods` | in `n_eff` units; outputs are null until it is reached. A list gives one threshold per target. Warm-up gates output, not learning |
| `coef_every` | snapshot the coefficients every N rows (`0` = only on each chunk's last row) |
| `label_delay` | hold each row back from *learning* until the clock has moved this much further on. For a target that is only known later — see [Labels that arrive late](#labels-that-arrive-late) |

`n_eff` is the exponentially weighted observation count: the weight behind
the state that produced *this row's* prediction, measured before the row's
own update and decay. So it is `0` on a stream's first row, lags the row
count by one, saturates at `1 / (1 − λ)`, and means the same thing in every
model — which is what makes `min_periods` portable across a bank.

### Nulls

A null in any feature, or in the weight, skips the row: outputs are null, no
update happens, the clock still advances. A null in one target still emits
that target's `pred`, leaves its `resid` null, and skips only that target's
update. NaN, ±inf and any magnitude above `1e100` count as null, so sentinels
never reach a model.

### Three ways to hold a row back

They differ, and it matters which you use.

- **Weight `0`** — the row is scored, the clock advances, nothing is learned:
  the coefficients are frozen bit for bit. Use it to keep a row's place in
  the stream. Since the clock advances, `n_eff` keeps decaying and can fall
  below `min_periods` if you score for a long stretch this way.
- **A null target** — the feature moments still update while the target's
  cross-moment does not, so the coefficients wander with feature noise. Use
  it only for a label that has not arrived yet.
- **`predict`** — scores every row against the bank exactly as it stands and
  touches nothing: no clock advance, no decay, `n_eff` frozen. Use it to
  serve. It is also the fast path: `ewridge` scores at 1.8–2.9× its learning
  throughput.

### Any row order

A bank is a set of sufficient statistics, so row order reaches the fit only
through decay. With decay off, an `ewridge` with `ridge=0` is ordinary least
squares over every row it has seen, in whatever order they came: forwards,
backwards or shuffled, its coefficients match `numpy.linalg.lstsq` to 2e-13.
Memory use is proportional to the model's state, not to the amount of data
that has passed through it — 6M rows × 20 features from a parquet stream peak
at 1.4 GB, against 3.97 GB for `lstsq` on the same rows, and the frame never
has to fit. Set `solve_every=1000` to solve less often than every
row (1.4 s instead of 11 s there, coefficients at most 1000 rows stale).

A finite halflife with no clock discounts each row by how far back it sits,
so the fit is a weighted least squares of *that* order. One trap: a huge
finite halflife is not `inf`. Its solve cadence defaults to `halflife/50`, so
`halflife=1e12` solves once, at `min_periods`, and never again. Say `inf`, or
set `solve_every`.

### Two guarantees

- **Predictions are out-of-sample.** Every row is predicted from the state
  before its own target is folded in. Nothing here can leak.
- **Chunk invariance.** One chunk or a thousand, with or without a save and
  resume in the middle, the output is bit-identical. The one exception is
  `coef`, which is a *reporting* cadence: it is snapshotted every
  `coef_every` rows and on each chunk's last row, so smaller chunks report it
  more often.

### Mistakes are named

A builder checks each keyword against its type hints — `halflife="10"` says
`spec "m": halflife must be a number or a list of numbers, got str '10'`. A
missing column says which spec wanted it, in what role, and what the frame
has. A spec named like an input column is refused rather than silently
replacing it.

## Running a bank

### In a loop: `ModelBank`

```python
spec = po.spec.ewridge(
    "ridge",
    targets=["y"], features=["x0", "x1", "x2"],
    clock="t", halflife=600.0, max_dclock=300.0,
    group="bond_id", ridge=[1e-6, 0.1], standardize=True,
)
bank = po.ModelBank([spec])

for chunk in lf.collect_batches():        # never materializes the whole stream
    out = bank.fit_predict(chunk)
    ...

bank.save("bank.state")                    # atomic: temp file, then rename
```

A bank says what it holds: `repr(bank)` is
`ModelBank(['ridge'], groups=412, rows_seen=3000000)`, `bank.specs` gives
back the spec dicts, and `bank.groups()` is a frame of every `(spec, group)`
with its row count and last clock value. Groups live until dropped, so a
long-running bank forgets the quiet ones with:

```python
stale = bank.groups().filter(pl.col("last_clock") < now - 30 * 86400)
bank.drop_groups(stale["group"])           # they start cold if they reappear
```

### As a query: `lf.online.fit_predict`

The loop above as a `LazyFrame`. Executing it — `collect()`,
`collect_batches()`, any `sink_*()` — streams the plan's rows through a
fresh bank in `chunk_rows` chunks, so the query stays at *state + one chunk*
however long the stream, and everything after the bank is ordinary polars:

```python
(
    pl.scan_parquet("ticks/*.parquet")
    .online.fit_predict([spec], chunk_rows=100_000)
    .filter(pl.col("ridge").struct.field("n_eff") > 100)  # after the bank: what comes out
    .select("ts", "bond_id", "ridge")                     # pushed into the scan
    .sink_parquet("fitted.parquet")
)
lf.online.predict(bank).collect()                                         # serve
lf.online.fit_predict(load_state="bank.state", save_state="bank.state")   # resume, and save
```

Things worth knowing about the plan:

- **It is pure.** Every execution starts from the same state (the specs',
  or `load_state`), so collecting twice gives the same frame, and `head(n)`
  learns from the first `n` rows and no more.
- **`save_state` writes when the run reaches the last row**, atomically, and
  the same bytes a `ModelBank` or [`po.run`](https://hgilde.github.io/polars-online/polars_online.html#polars_online.run) would write. A run abandoned early
  or ended by a bank error leaves the file untouched. A failure *after* the
  bank does not stop the bank, so the state is written although the query
  failed; if the two must be tied together, `po.run` saves only after its
  output is committed. [docs/STATE-WORKFLOW.md](docs/STATE-WORKFLOW.md) has
  the measurements behind each of these.
- **Filters, selections and `head` after the bank are honoured at the
  source**, and a selection reaches the input scan, so only the columns the
  specs and the query need are read.
- **Filter after the bank, not before, unless the model must skip those
  rows.** A filter after never changes what the bank learns from, and it
  streams; a filter *before* holds several row groups per thread in polars'
  parquet reader — 2.5 GB at 12M rows against 0.78 GB for the same filter
  after ([docs/PERFORMANCE.md](docs/PERFORMANCE.md) §11).

If the model must not learn from some rows, give them weight `0` instead of
filtering them out: they still stream, still come out scored, and no gap
opens in the clock.

```python
(
    lf.with_columns(pl.when(pl.col("venue") == "X").then(1.0).otherwise(0.0).alias("w"))
    .online.fit_predict([po.spec.ewridge("ridge", targets=["y"], features=["x0", "x1"],
                                         clock="t", halflife=600.0, max_dclock=300.0,
                                         weight="w")])
    .sink_parquet("fitted.parquet")
)
```

`df.online.fit_predict(specs)` is the eager twin. `po.fit_predict(frame, ...)`,
[`po.predict(frame, bank)`](https://hgilde.github.io/polars-online/polars_online.html#polars_online.predict) and [`po.unnest(frame, specs)`](https://hgilde.github.io/polars-online/polars_online.html#polars_online.unnest) are the same calls as
plain functions, for a type checker, which cannot see a registered namespace.

### As a job: `po.run` and the `online` CLI

The same bank as a three-stage pipeline — read, fit, write, one chunk in
flight per stage — with the output written through a temporary file and
renamed into place, so a failed run leaves the previous output where it was.

```python
po.run(input="ticks.parquet", output="fitted.parquet",
       specs=[spec], chunk_rows=100_000, save_state="bank.state")   # -> {"rows": ..., "chunks": ...}

po.run("bank.toml", input="today.csv")                              # keywords override the TOML
po.run(input="today.parquet", output="scored.parquet", specs=[spec],
       load_state="bank.state", predict=True)                       # serve: learn nothing
```

`input` is anything py-polars can stream:

- a **path** in parquet, ipc, csv or ndjson — told from the extension, or
  named with `input_format=`; globs and cloud URLs as `pl.scan_*` takes them;
- a **`LazyFrame`**, with whatever scan options its query needs;
- a **`DataFrame`**;
- **any iterable of frames** in stream order — a database cursor, a socket,
  a generator.

`output` is a path in any of the four formats. `keep_columns=[...]` selects
input columns before the bank sees them. `progress(rows, chunks)` is called
after each chunk; raising in it stops the run. CSV cannot hold struct columns,
so there each spec's struct is flattened to `<spec>.<field>` columns and the
`coef` list becomes a JSON string that
`pl.col("ridge.coef").str.json_decode(pl.List(pl.Float64))` reads back
bit-exact.

**When the product of a run is its state, leave `output` out.** An
accumulator-only spec emits `n_eff` a row and nothing else; over a billion
rows that is 8 GB of file written so it can be deleted. Without an `output`
the run writes nothing and `save_state` is required — a run that writes
nothing and saves nothing has done nothing:

```python
wide = po.spec.ew_cov("gram", features=[f"x{i}" for i in range(3)], stats=[], halflife=1000.0)
po.run(input="ticks.parquet", specs=[wide], save_state="gram.state")   # no output at all
```

`no_output=True` says the same thing over a config that names an output, and
is what the CLI's `--no-output` sets.

The CLI is the same pipeline as one binary and one TOML
([examples/bank.toml](examples/bank.toml)), for deployments with no Python:

```sh
online --config bank.toml
online --config bank.toml --resume bank.state --save-state bank.state
online --config bank.toml --resume bank.state --predict --input today.parquet
online --config bank.toml --input ticks.csv --output scored.ndjson
online --config bank.toml --input feed.dat --input-format ipc
online --config bank.toml --no-output --save-state gram.state   # the state is the product
online --config bank.toml --dry-run          # validate and print the output schema
```

`--predict` scores against the resumed state and learns nothing; it drops
the config's `save_state`, so one TOML serves both runs. `--no-output`
suppresses the per-row output; a run needs one or the other.

A TOML spec for a model that learns from no target — `ew_cov`, `kmeans`,
`micro` — may leave `targets` out, and it is filled with `features[0]` the
way the Python builders fill it. The two surfaces then write byte-identical
specs, so a state saved from one resumes under the other. The CLI reads with
polars' own scanners, which on a stable toolchain lack the SIMD CSV parser
py-polars' wheels have, so for a large CSV `po.run` is the faster of the two.
In TOML, a Windows path needs single quotes or forward slashes
(`input = 'C:\data\in.parquet'`), since a backslash in a double-quoted string
starts an escape sequence.

From Rust, the same pipeline is `online_polars::run_config` for a `RunConfig`,
`run_config_on` for a `LazyFrame` or batches the caller already has, and
`run` with a callback instead of an output file.

### The expression form (in memory only)

For a frame that is already in memory, the shortest way to write a model is
as an expression. Features may be expressions, evaluated per group under
`.over`, so a lag never crosses a group boundary:

```python
out = df.with_columns(
    pl.col("y").online.ewridge(
        features=["x0", "x1", pl.col("y").shift(1).alias("y_lag")],
        clock="t", halflife=600.0, max_dclock=300.0,
    ).over("group").alias("fit")
)
```

The numbers are the bank's — the expression *is* the bank, run over the
column polars hands it. And that is the catch: polars gives a stateful user
expression its whole column at once, in either engine, so wrapping the
expression in a lazy query does not make it stream. On the 12M-row file
below, it peaks at **7.3 GB against 1.35 GB** for the plan. Every call
therefore warns with [`polars_online.InMemoryExpressionWarning`](https://hgilde.github.io/polars-online/polars_online.html#polars_online.InMemoryExpressionWarning); using the
expression on a frame in memory on purpose is fine, and one line says so:

```python
import warnings

warnings.filterwarnings("ignore", category=po.InMemoryExpressionWarning)
```

[`po.online(pl.col("y"))`](https://hgilde.github.io/polars-online/polars_online.html#polars_online.online) is the same namespace as a plain function.
[docs/PLAN.md](docs/PLAN.md) §6 has the design and the condition under which
the warning would go away.

## Memory: which calls stream

All of them but the expression form. Peak footprint on one file, `ewridge`
with 20 features, parquet in and out:

| what you write | 3M rows | 12M rows | |
|---|---:|---:|---|
| `lf.online.fit_predict([spec])` | 0.90 GB | 1.35 GB | a query over a stream |
| `for chunk in lf.collect_batches(): bank.fit_predict(chunk)` | 0.80 GB | 1.24 GB | your own loop |
| `po.run(...)`, `online --config` | 0.95 / 0.73 GB | 1.41 / 0.75 GB | file in, file out |
| `pl.col("y").online.ewridge(...)` in `with_columns` | | 7.3 GB | the expression: whole column at once |

The first three are flat; what growth they show is the allocator holding
freed pages, and nearly all of the rest is polars' parquet read-ahead
(`POLARS_ROW_GROUP_PREFETCH_SIZE=1` takes it to 0.31–0.46 GB). Everything
after the bank — filters, joins, group-bys, sinks — streams as polars streams
it. Note that this is polars' rule for the operations *around* the bank too:
a rolling window under `.over("group")` or `group_by=` collects (6.5 GB and
1.7 GB on the same rows, against 0.25–0.28 GB ungrouped), whereas a bank's
`group=` is one accumulator per group and stays O(state).
[docs/PERFORMANCE.md](docs/PERFORMANCE.md) §11 has every measurement.

## Saving, loading and serving

A fitted model is the bank's *state* — one accumulator per `(spec, group)` —
and it travels as one file, written whole or not at all. The same words work
from a bank and from a query:

```python
# From a bank
bank.fit_predict(df)
bank.save("bank.state")                               # atomic: temp file, then rename
bank = po.ModelBank.load("bank.state", specs=[spec])  # specs= checks the file is this model's
bank.fit_predict(today)                               # learn on: the state moves
scored = bank.predict(today)                          # serve: score, learn nothing

# From a query
lf.online.fit_predict([spec], save_state="bank.state").sink_parquet("fitted.parquet")
lf.online.fit_predict(load_state="bank.state", save_state="bank.state").sink_parquet("more.parquet")
served = lf.online.predict("bank.state").collect()
```

`po.run(..., save_state=)`, `po.run(..., load_state=, predict=True)` and the
CLI's `--save-state` / `--resume` / `--predict` read and write the same file,
and the bytes are the same whichever wrote them. `save_bytes()` and
`load_bytes()` do the same in memory, for a checkpoint that lives somewhere
else. Loading names the problem it hits: `FileNotFoundError` for no file yet,
`ValueError` for a file that is not a bank, a newer build's, or another
model's.

What `predict` reports: row `i` carries what `fit_predict` would have
reported had it been the next row of the stream — `pred`, `n_eff`, `sigma`,
`resid_z`, selection, metrics, field for field — with every row scored from
the same state. The target column may be absent (then `resid` is null),
`weight` is not read, a group the bank has never seen scores null, and the
stream's session and clock policies still hold.
[docs/STATE-WORKFLOW.md](docs/STATE-WORKFLOW.md) walks the whole workflow:
fit, save, serve, learn on, with what each step guarantees.

### A state file describes itself

A saved bank can be read by something that knows nothing about it. No specs,
no config, no data — the file carries what it needs:

```python
bank = po.ModelBank.load("bank.state")   # no specs=: the file is enough

bank.specs                # every spec back as the dict its builder made
bank.groups()             # spec, group, rows_processed, last_clock
bank.output_fields()      # {'ridge': ['pred_y', 'resid_y', 'n_eff', 'coef'], ...}
bank.rows_seen()          # rows fed, over every chunk and group
bank.solve_failures()     # per spec, per group
```

That is enough to walk one: list the specs, read each one's model type and
parameters, ask what fields it emits, list its groups. Then the four
diagnostic tables say how it is doing and what it was trained on —
[`last_row()`](#the-last-row) and [`coef()`](#coefficients) for the fit,
[`summary()` and `describe()`](#what-it-was-fed) for the stream behind it.
Each returns every spec by default, with `spec` as the first column, and the
schemas are fixed across specs, so banks from different runs stack with a
plain `concat`.

`bank.specs` is a **copy, and read-only**. The bank's behaviour comes from
the state built at construction, so a list on the Python side could only ever
disagree with it — and used to: editing `bank.specs[0]["features"]` in place
left `coef()` labelling coefficients from a spec the bank was not running.

### Reading a state without this library

`bank.to_json()` is everything `save` writes, as JSON, and `save_json(path)`
puts it in a file. Use it to look at a state, diff two of them, or hand one
to something that is not Python.

It is an export, not a second format — `load` reads msgpack and only
msgpack — and it is faithful, including the values JSON has no literal for.
`NaN` and `±inf` are written as `"nan"`, `"inf"` and `"-inf"`, the same
spelling a spec's `halflife` already uses. That matters more than it sounds:
`halflife=inf` means no decay, so an ordinary state carries an infinity, and
a plain JSON encoder writes it as `null` without saying so. Every export is
read back and checked against the state before you get it, so a state that
could not be carried is an error rather than a file that is quietly wrong.

## Preparing a stream

Two things a stream may need before a bank sees it: a target that is not
known at the row it sits on, and series that do not tick together.

### Labels that arrive late

A target that is a forward quantity — the next five minutes' return, the
next day's fill rate — is not known at the row it sits on. A stream that
learns it there hands the model that much of the future before it predicts
the rows in between. Every "out-of-sample" number after that is
contaminated, and with an autocorrelated feature even a pure noise column
starts to look predictive.

`label_delay` is the fix, and it is one parameter:

```python
spec = po.spec.ewridge("fwd", targets=["ret_5m"], features=["x0", "x1"],
                       clock="ts", max_dclock=3600.0, halflife=1800.0,
                       label_delay=300.0)     # the return takes 5 minutes to be known
```

Each row is **scored** where it sits and **learned from** 300 clock units
later. Everything downstream of the label moves with it: the prediction,
`sigma`, `resid_z`, the metrics, drift, the conformal interval, `n_eff` and
`min_periods` all see only labels that had really arrived.

The clock is the model's own — the raw column capped by `max_dclock`, with
skipped rows' time folded in — which is what makes release depend on the
clock alone and so survive any chunking. With no `clock` column that is one
unit per accepted row, so `label_delay=20` is twenty rows. Two events empty
the buffer: a **reset** drops it (the state those rows would teach is being
thrown away), and a **session change** releases it in order (one session's
clock does not measure time in the next). Rows still waiting when the stream
ends are simply never learned from — their labels never matured.

The buffer lives in the state and is saved with it, so a run that stops
mid-stream resumes with the same rows still waiting. The memory it needs is one
row's values for every row inside the delay, per group.

[`po.prep.embargo`](https://hgilde.github.io/polars-online/prep.html#polars_online.prep.embargo) writes the same thing out as data: every row twice, a
zero-weight prediction at `t` and a lesson at `t + delay`, merged back into
clock order. That is the recipe to reach for when the delay has to be
visible in the frame, or for an engine that is not this one. The native path
is tested against it field by field and agrees to the bit, with three
exceptions: `resid_quantiles`, `emit_autocorr` and `emit_drift`. Those three
take no row weight, so a zero-weight row feeds them as much as its learn copy
does — and in a doubled stream every residual therefore lands twice.
`label_delay` feeds them once.

### Series that tick at their own times

Two series observed at different instants cannot be correlated directly. A
fine common grid attenuates the correlation towards zero (the Epps effect)
and filling forward invents observations that were never made.

[`po.prep.refresh_time`](https://hgilde.github.io/polars-online/prep.html#polars_online.prep.refresh_time) puts them on the grid Barndorff-Nielsen, Hansen, Lunde
and Shephard defined: a point wherever **every** series has ticked at least
once since the last one, each carrying its last observed value.

```python
from polars_online import prep

grid = prep.refresh_time(ticks, series="symbol", names=["AAA", "BBB", "CCC"],
                         time="t", value="px").collect()
```

The input is long — one row per tick, with the series named in a column. The
output has one row per grid point: `time_refresh`, one `<s>_value` per
series, `n_obs_<s>` since the previous point, and `retained_fraction`,
which is how much of the data survived. Look at that last one before
trusting a correlation computed on the result: the grid runs at the pace of
the slowest series, so a fast one loses most of its ticks.

`pairs=True` runs an independent two-series grid per pair instead, which
keeps far more when one series is slow. Rows must be in time order; a
backwards time is an error naming the row, and nothing is interpolated.

The output looks synchronous and is not: each value is up to one of its own
inter-tick intervals old. `n_obs_<s>` is that staleness made visible — the
series with the largest count is the one holding the grid up.

## Reading the fit

What a bank can tell you about its fit — and about what it was fed — with
no data at hand, and the grammar of the field names it writes.

### Coefficients

Two ways, and they agree row for row:

```python
ols = po.spec.ewridge("ols", targets=["y"], features=["x0", "x1"], clock="t",
                      halflife=600.0, max_dclock=300.0, group="bond_id", coef_every=1)

# 1. From a bank -- live, or loaded from a state file with no data at hand.
#    One row per coefficient, with the term it belongs to.
bank = po.ModelBank([ols])
bank.fit_predict(df)
betas = bank.coef()                  # every spec: spec, group, instance, n_eff, ..., term, coef
wide = betas.pivot("term", index=["group", "instance"], values="coef")

# 2. From the output, as columns: the fit as it moved, one row per row.
path = (
    lf.online.fit_predict([ols])
    .online.unnest([ols])            # pred_y, resid_y, n_eff, coef_y_intercept, coef_y_x0, coef_y_x1
    .select("t", "bond_id", "^coef_.*$")
    .collect()
)
```

`bank.coef()` is the fit as of the last row each group learned from, with
`n_eff` for how much weight is behind it. The output's `coef` is the same fit,
snapshotted *after* each row's update — the row's own `pred` comes from the fit
*before* it. It is written every `coef_every` rows, and on the last row of
every chunk. The default, `coef_every=0`, is the chunk end only; the per-row
path above asks for `coef_every=1`, and writes a list of `k` floats on every
row of the output.

Under a grid — several `ridge` values, `feature_sets`, a `lasso_path`,
several targets — the list holds one block per (target × grid point).
`unnest` names each block's columns the way the `pred` fields are named
(`coef_y_x0__r0.5@h500` beside `pred_y__r0.5@h500`), `bank.coef()` carries
the same columns to tell blocks apart (add them to the pivot's `index`), and
`unnest` reads a saved output the same way:
`pl.scan_parquet("fitted.parquet").online.unnest([ols])`. It takes the specs,
a bank, or the path of a saved state. `bank.gram("ols")` gives the EW
accumulators behind the fit, for anything other than our solve.

### The accumulators behind a fit

`bank.gram(spec)` hands back the matrices the model itself solves against,
per group and decay instance. They are a *complete* sufficient statistic, so
a saved state answers questions the run never asked:

| key | what it is |
|---|---|
| `means`, `comoments` | the feature means and the centred `k × k` co-moment matrix |
| `cross_moments`, `target_weights` | per target, the uncentred `E[z·y]` the solve consumes, and the weight behind it |
| `target_means`, `target_vars` | per target, the target's own mean and centred variance |
| `n_eff`, `n_kish`, `target_n_kish` | the accumulated weight, and Kish's effective sample size for the features and per target |

`n_eff` counts weight, not rows. `n_kish = n_eff² / Σw²` is the number of
equally weighted rows the moments are worth, which is what a standard error
divides by; an exponentially weighted window of unit rows settles at `(1 +
λ)/(1 − λ)` whatever the halflife's units. `n_kish` is scale-free: multiply
every weight by the same factor and it does not move, so it does not fall when
a stream goes quiet. `n_eff` is the number that falls, and the one to read for
that.

The target moments are the half that makes the rest usable. Without `Var[y]`
there is no residual variance, no R², no information criterion and no
standard error to be had from a saved Gram:

```python
ols = po.spec.ewridge("ols", targets=["y"], features=["x0", "x1", "x2"],
                      halflife=500.0, ridge=1e-9, standardize=False)
fitted = po.ModelBank([ols])
fitted.fit_predict(df)

g = fitted.gram("ols")[0]
raw = g["comoments"] + np.outer(g["means"], g["means"])   # the solve's pairing
beta = np.linalg.solve(raw, g["cross_moments"][0])
slopes = beta[1:]                                          # column 0 is the intercept
resid_var = g["target_vars"][0] - slopes @ g["comoments"][1:, 1:] @ slopes
r2 = 1 - resid_var / g["target_vars"][0]
```

A state saved by 0.1.x has no `Σw²` and no target moments, and
they cannot be recovered from what it does have — those four keys are `None`
there, for that state's whole remaining life. Reading any of this back needs
`numpy`, which is an optional extra.

`po.gram` is the toolkit for these, so the algebra above is one call:

| call | what it does |
|---|---|
| `merge(grams)` | pools the Grams of disjoint row sets into the Gram of their union, exactly |
| `subset(g, cols)` | the Gram of some of the columns — a sub-block, not a recomputation |
| `solve(g, ridge=, target=, features=)` | the model's own ridge, in original units; a list of ridges is one eigendecomposition |
| `lasso_path(g, lambdas, ...)` | the `lasso` model's coordinate descent, offline, plus per-feature `penalty_weights` |
| `coef_stats(g, coef)` | residual variance, R², standard errors and t, at the Kish sample size |
| `correlation(g)`, `vif(g)`, `condition(g)` | the correlation matrix, variance inflation factors, and Belsley's condition indexes and variance-decomposition proportions |

```python
ols = po.spec.ewridge("ols", targets=["y"], features=["x0", "x1", "x2"],
                      halflife=500.0, ridge=1e-9, standardize=False)
fitted = po.ModelBank([ols])
fitted.fit_predict(df)
g = fitted.gram("ols")[0]

r2 = po.gram.coef_stats(g, po.gram.solve(g, ridge=1e-9))["r2"]          # the block above
ridges = po.gram.solve(g, ridge=[0.0, 0.01, 0.1, 1.0], standardize=True)  # one decomposition
worst = po.gram.condition(g)["kappa"]
```

It is the same arithmetic the models run, so `solve` on a spec's Gram is
that spec's fit and `lasso_path` is that spec's path — the tests hold both
against the models rather than against a second copy of the formula. What it
is not is the same arithmetic to the last bit: the models factorize with
`faer`'s Cholesky and numpy with LAPACK's LU, which round differently in the
last place or two.

Two things to know before reading a disagreement as a bug. `bank.coef()` is
as of the model's last *solve*, which its `solve_every` schedule decides,
while `gram()` is as of the last row. And `merge` pools parts that share a
weighting — shards of a pass, groups being combined — not two halves of a
decayed stream in time order, where each part's weights are relative to its
own last row; the docstring gives the rescaling for that case.

### One row per finished group

A bank keeps one state per group key, for the life of the bank. On a stream
whose key space keeps growing — a day id, a session id, a block number — that
is unbounded memory for state nobody will read again.

`group_close` says when a group is finished. The bank then emits its
accumulators as one row and drops the stream.

```python
blocks = po.spec.ew_cov("cov", features=["x0", "x1"], lam=1.0,
                        group="block", group_close="monotone")
by_block = df.with_columns(block=pl.int_range(pl.len()) // 100)

bank = po.ModelBank([blocks])
bank.fit_predict(by_block)

closed = bank.closed_groups()          # one row per finished block
first = po.gram.from_row(closed.head(1))
corr = po.gram.correlation(first)      # everything in po.gram works on it
```

`"monotone"` means the key column never goes backwards, so a key smaller
than the largest one fed so far is finished; a chunk whose keys are out of
order is refused, naming the row. An integer key column is ordered as
numbers; anything else is ordered as text, so `"9"` comes after `"10"`. Sort
by the same rule the bank reads, or the chunk is refused: a `Categorical`
column sorts by its physical order by default, which is the order the
categories were first seen, so sort it with `pl.col("k").cast(pl.String)` or
cast the column itself. `"session"` closes a group where its `session` value
changes. Either way the last group never closes — nothing proves it is
finished — and stays readable through `gram()`.

The row is the `gram()` a driver would have read at that moment, bit for
bit: one builder makes both. It carries the span's own `rows_fed`,
`rows_learned` and clock range beside the moments, `coef` for a model that
has one, the eigendecomposition for an `ew_cov` with `pca`, and a
`marginal`'s pairs.

A run writes the same rows to a sidecar file, which with no output at all is
the whole shape of an accumulate-only pass — read a stream that does not fit
in memory, write one row per block:

```python
po.run(input=by_block.lazy(), specs=[blocks], closed_groups="blocks.parquet")
```

(`by_block` and `blocks` are from the block above.)
`lf.online.fit_predict(closed_groups=path)` and `online --closed-groups path`
write it too. What has closed and not been read is saved with the state, so
a driver that saves between chunks does not lose rows silently.

### Reading a correlation matrix

`po.gram` solves and diagnoses a design matrix. `po.corr` is its
complement: the arithmetic that comes *after* a correlation matrix, in the
same style — numpy only, pure functions, each held against the paper it
comes from.

```python
r = po.corr.matrix(closed.head(1))       # an array, a gram() dict, or a closed row
fixed, dist, iters = po.corr.nearest(r)  # Higham's nearest correlation matrix
shrunk, alpha = po.corr.shrink(r, alpha=0.2)
lo, hi = po.corr.mp_edge(n=2000, m=50)   # where pure noise puts its eigenvalues
```

| function | what it does |
|---|---|
| `to_z`, `from_z` | Fisher's transform, clipped so a degenerate ±1 is finite |
| `nearest(A, W)` | Higham (2002): the nearest correlation matrix, by alternating projections with Dykstra's correction. Returns the matrix, the distance and the iteration count |
| `shrink(R, target, alpha, x)` | Ledoit–Wolf shrinkage towards a constant-correlation or identity target; the optimal intensity needs the rows, and the docstring says why |
| `equicorr`, `equicorr_row`, `equicorr_loglik` | `deco`'s three quantities, offline |
| `absorption`, `shift` | the absorption ratio and its standardised shift |
| `spectral`, `from_spectral` | the top eigenpairs and the completion back to a correlation matrix |
| `block_means`, `from_blocks` | mean correlation within and between labelled blocks, and back |
| `mp_edge`, `mp_density` | the Marchenko–Pastur edges and density: which eigenvalues are noise |
| `signal_share` | how much of a correlation's movement across blocks is not the sampling floor |
| `loss` | `qlike`, `z_mse` or Engle–Colacito `minvar`, each zero or minimal at the truth |
| `epps_invert` | the correlation at a coarser scale from `ew_cov`'s lagged co-moments |
| `fisher_se` | the standard error of a correlation, with the AR(1) inflation and its caveats |

### The last row

The state file also carries the output row of the last row each group
learned from, so a saved model says how it was doing without its output
frame:

```python
bank = po.ModelBank.load("bank.state", specs=[spec])
last = bank.last_row("ridge")    # one row per group: spec, group, pred_y__r0.000001, ..., n_eff, coef
```

It is the `fit_predict` row field for field — `pred`, `sigma`, the metrics
and the interval when the spec asks for them, `n_eff`, and `coef` when that
row carried it (a chunk's last row does; `bank.coef()` has the coefficients
either way). Fit many models, save each, and comparing them is one `concat`
over the files:

```python
from pathlib import Path
table = pl.concat(
    [po.ModelBank.load(f).last_row() for f in sorted(Path(".").glob("*.state"))],
    how="diagonal_relaxed",          # specs with different fields stack with nulls
)
```

A group that has not learned from a row yet, or a file written by 0.1.x,
gives a row of nulls, and `predict` does not move the row.

### What it was fed

The state file also carries what each group has seen: how many rows, what
became of them, the clock's range, and a count and moments for every input
column. A saved model can say what it was trained on without the data at
hand:

```python
bank = po.ModelBank.load("bank.state", specs=[spec])
fed = bank.summary("ridge")    # one row per group: rows_fed, rows_processed, rows_skipped, rows_learned, ..., clock_min, clock_max
cols = bank.describe("ridge")  # one row per input column per group: column, role, count, null_count, mean, std, min, max
```

`summary` counts rows. `rows_fed` were routed to the group; `rows_processed`
the model accepted, `rows_skipped` it did not (a feature or the weight was
missing); `rows_learned` moved the fit (a weight above zero and a target
present) and `rows_zero_weight` advanced the clock and nothing else. With
them come `weight_sum`, the clock's first and last values, `last_clock`,
and what the clock schedule met: `session_changes`, `clock_backwards`,
`resets`. `describe` is `DataFrame.describe` for each feature, target and
weight column over the rows fed, counting as the models count: a null, a
NaN, an infinity or a magnitude beyond `1e100` is a `null_count`, not a
value. A label column has counts only, and an unsupervised model lists its
features.

Neither decays: they are plain counts over the whole stream, computed in
row order, so they are the same whatever the chunking, to the bit. `predict`
does not move them. A file written before 0.2.0 reports nulls for both — a
count that began partway would read as the whole history — while
`rows_processed` and `last_clock`, which the stream always kept, are
filled in.

### Output field names

You index the result struct by strings, so the names are a contract. The
grammar:

```
pred_{target}{combo}{instance}     combo    = ""            single ridge, no feature sets
resid_{target}{combo}{instance}             | __r{ridge}     ridge grid
sigma_{target}{combo}{instance}             | __{set}        feature sets, single ridge
absresid_q{level}_{target}...               | __{set}_r{ridge}
n_eff{instance}                    instance = ""            single halflife
coef{instance}                              | @h{halflife}   halflife grid
```

Numbers render as plain decimals in `[1e-6, 1e7)` and as compact scientific
outside it. Every name, default and signature is pinned by
`tests/test_api_surface.py` against a checked-in snapshot, so a change is a
reviewable diff and a version bump, never a silent rename of your columns.

You never have to build these strings. [`po.spec.output_index(spec)`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.output_index) lists
every field with the values its name encodes, and [`po.spec.coef_fields(spec)`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.coef_fields)
does the same for the `coef` lists — one row per coefficient with its list,
position, and the column `unnest` gives it — so selecting is a filter, not
string formatting:

```python
grid = po.spec.ewridge("m", targets=["y"], features=["x0", "x1"], clock="t",
                       max_dclock=300.0, halflife=[100.0, 500.0], ridge=[1e-6, 0.5])

idx = po.spec.output_index(grid)
name = idx.filter((pl.col("kind") == "pred") & (pl.col("target") == "y")
                  & (pl.col("ridge") == 0.5) & (pl.col("halflife") == 500.0))["field"].item()
out["m"].struct.field(name)                                    # "pred_y__r0.5@h500"

row = po.spec.coef_fields(grid).filter(
    (pl.col("term") == "x1") & (pl.col("ridge") == 0.5) & (pl.col("halflife") == 500.0)
).row(0, named=True)
out["m"].struct.field(row["field"]).list.get(row["position"])  # field "coef@h500", position 5
```

Both tables come from the same Rust code that renders the names, so they
cannot drift from the strings, and the index carries each field's `dtype` —
the schema the bank declares to polars before the first row is read. One
sharp edge: avoid `__` and `@` in target names and feature-set labels if you
parse field names downstream, since a target named `y__r0.5` renders like a
ridge grid on `y`.

## Diagnostics, selection and evaluation

Opt-in outputs, all derived from state the models already keep and all read
*before* each row, so they are as out-of-sample as the predictions they
describe:

| flag | adds | meaning |
|---|---|---|
| `emit_sigma` | `sigma_<slot>` | EW standard deviation of that slot's out-of-sample residuals |
| `emit_resid_z` | `resid_z_<slot>` | `resid / sigma` — how surprising the row was, in units of the model's own recent error |
| `emit_selected` | `selected_<t>`, `pred_<t>__selected` | online model selection across ridge values, feature sets and halflives, by lowest EW out-of-sample error |
| `emit_averaged` | `pred_<t>__averaged` | `softmax(−eta · EW error)` blend over the same slots — hedges where `emit_selected` commits |
| `emit_drift` | `drift_<slot>` | Page-Hinkley break detection on the residual stream; `drift_action="reset"` also restarts the stream |
| `emit_metrics` | `ic_<slot>`, `r2_<slot>`, `hit_rate_<slot>` | the numbers [`po.eval`](https://hgilde.github.io/polars-online/eval.html#module-polars_online.eval) computes, kept in O(state) beside the model |
| `resid_quantiles` | `absresid_q<p>_<slot>` | P² quantiles of \|resid\| — a distribution-free interval where `sigma` gives a Gaussian one |
| `emit_autocorr` | `autocorr_<slot>` | EW residual autocorrelation; non-zero means the model is mis-specified |
| `conformal=0.9` | `lo_<slot>`, `hi_<slot>`, `coverage_<slot>` | an adaptive conformal interval at that coverage, distribution-free, and the coverage it has actually delivered |

`emit_metrics` on a `sgd` or `ftrl` fit with `loss="logistic"` reads
differently, because `pred` is a probability and `y` a 0/1 label rather than
a signed target: `hit_rate` is accuracy at a 0.5 threshold, not sign
agreement (the sign test always agrees on two positive numbers — it read
1.0 for every such fit before this was fixed), `r2` is the Brier skill score
against the running base rate, and `ic` the point-biserial correlation
between the probability and the label, both under their usual names. There
is no streaming log loss (a `ln` result cannot go into the state without
losing cross-platform reproducibility — see `docs/PLAN.md` §11a);
`po.eval.metrics(..., binary=True)` adds it over a collected frame instead.

Drift detection complements the halflife rather than replacing it: decay
forgets smoothly and always; a detector notices a break and says so, within
a couple of rows of a sign flip.

`conformal` is the interval to use when the residuals are not Gaussian. It
tracks the `coverage` quantile of `|resid|` directly: the radius grows by
`conformal_rate · sigma · coverage` on a miss, and shrinks by
`conformal_rate · sigma · (1 − coverage)` on a hit. So its long-run coverage
is the number you asked for, whatever the residuals do, with an error that
shrinks like `1/T`.

`sigma` gives a Gaussian interval instead. On Gaussian residuals the two
agree; on fat-tailed or heteroskedastic ones the Gaussian interval
over-covers by several points where this one lands on target. The conformal
radius starts at `sigma · Φ⁻¹(1 − α/2)` and is null until then. It is read
before the row, like everything else, and adds three numbers to the state per
slot.

```python
ci = po.spec.ewridge("ci", targets=["y"], features=["x0", "x1"], clock="t",
                     max_dclock=300.0, halflife=500.0, conformal=0.9)
band = po.ModelBank([ci]).fit_predict(df).unnest("ci")
held = band.select(((pl.col("lo_y") <= df["y"]) & (df["y"] <= pl.col("hi_y"))).mean())
```

After the fact, `po.eval` reads the output frame:

```python
po.eval.metrics(out, "ridge", by=["bond_id"])                       # R², IC, hit rate, MSE
po.eval.rolling_metrics(out, "ridge", clock="t", window=3600.0)     # per clock window
po.eval.compare_specs(out, ["ridge", "kalman"])                     # one table, many specs
po.eval.seqtest(out, a="kalman", b="ridge", by=["bond_id"])        # is kalman closer? evidence per row
```

`compare_specs` says which spec had the lower error over the frame.
`seqtest` says how sure you can be, at every row, and is the same test the
[`seqtest`](#seqtest--a-sequential-test-of-a-sign-by-betting) model runs
inside a bank, streaming.

Those four all need the whole frame. When the output is never materialised —
fifty slots over a billion rows — reduce each chunk instead and keep ten
doubles per key:

```python
ridge = po.spec.ewridge("ridge", targets=["y"], features=["x0", "x1"],
                        clock="t", max_dclock=300.0, halflife=500.0, group="bond_id")
scoring = po.ModelBank([ridge])

running = None
for chunk in df.iter_slices(100):
    part = po.eval.sums(scoring.fit_predict(chunk), "ridge", by=["bond_id"])
    running = part if running is None else po.eval.merge_sums(running, part)

po.eval.from_sums(running, min_obs=10)   # R², IC, hit rate, MSE and RMSE
```

`from_sums(sums(df))` is `metrics(df)`, and `merge_sums` over any split of
the rows is `sums` over all of them — both held by tests. The sums are
**centred** (weighted means and centred second moments, merged with a
parallel-axis term) rather than raw `Σy` and `Σy²`: a target sitting around
1e8 with unit spread destroys the raw form's variance entirely, and this one
does not notice. `weight=` names a column to weight the rows by.

### Data whose truth is known

A regime detector is a claim about a stream, and a claim needs a stream
whose answer is written down. [`po.sim.regimes`](https://hgilde.github.io/polars-online/sim.html#module-polars_online.sim) produces one from a seed,
with the awkward parts included — series that tick at their own times,
prices observed with noise, autocorrelated returns, a volatility that moves
with the regime, an intraday pattern and a volume clock — and hands back the
truth beside the data.

```python
out = po.sim.regimes(4, states=[0.2, 0.7],
                     transition=[[0.98, 0.02], [0.02, 0.98]],
                     n_blocks=8, rows_per_block=500,
                     phi=0.3, noise=0.01, async_rates=[1.0, 1.0, 0.4, 0.4],
                     seed=0)
rows, truth = out["rows"], out["truth_blocks"]
```

`rows` is what a consumer sees: **levels** `x_1 … x_m` (so
`po.prep.refresh_time` and then `.diff()` apply), a clock, a session and an
optional volume. `truth_rows` gives the block, state, volatility multiplier
and interpolation fraction per bar; `truth_blocks` gives each block's true
correlation matrix. Two calls with the same seed are byte-identical.

Under `async_rates` a bar where a series drew no tick carries `null`, so
both the sparse and the previous-tick forms are one line away. `durations`
makes each state last exactly as long as it says; `design="smooth"`
interpolates the matrix across a boundary instead of stepping.

## Models

All accumulators are exponentially weighted **means**, not sums, so they stay
bounded over arbitrarily long runs; second moments are kept **centered** (a
weighted Welford update), so the variance is right even when features sit on
a large offset. `z` denotes `[1, x]` when an intercept is configured, `w` the
row weight, `λ` the row's decay.

Each model links to its builder in the API reference, whose docstring lists
every keyword with its default, and to its Rust source under
`crates/online-core/src/`, where the module comment states the recursion.

### `ewridge` — EW ridge on sufficient statistics

*API:* [`po.spec.ewridge`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.ewridge) — *Rust:* [`ewridge.rs`](crates/online-core/src/ewridge.rs) — *Outputs:* [fields](docs/OUTPUTS.md#ewridge)

```
W'   = λW + w                       S' = (λW·S + w·z zᵀ) / W'
W_j' = λW_j + w                     r_j' = (λW_j·r_j + w·z·y_j) / W_j'
solve:  (S + ridge·D) β_j = r_j     D = I minus the intercept slot
```

O(k²) per row; Cholesky solves on a schedule (`solve_every` in clock units,
default `halflife/50`, every row for `halflife=inf` and for `lam`;
`max_rows_between_solves` caps it in rows). Ridge values and named
`feature_sets` are expanded at solve time from the same accumulator, so grids
are nearly free. `coef_prior` shrinks toward a stated belief rather than toward
zero; because `S` is a mean, a plain `ridge` is a permanent per-observation
penalty, and the fading warm start ("start at yesterday's fit") is
`ridge_decay`. With `standardize`, the solve is done in correlation form and
unscaled afterwards, dropping near-zero-variance features rather than
blowing up.

The O(k²) per row is the update of `S`, and at a thousand features it is
most of the cost. `gram_block_rows=256` holds 256 rows back and folds them
into `S` with one matrix product instead of 256 rank-one updates. Measured
single-threaded, that is 5.1× the rows per second at 256 features, 6.6× at
1,000 and 5.9× at 2,000; a solve every 512 rows brings each down to about
4×, because the solve costs the same either way. `n_eff`, the timing of
every prediction and chunk invariance do not change. The coefficients
agree with the row-by-row fit to rounding, not to the bit: the merged sum
is the same sum in a different order. The held rows travel in the state
file, so a save mid-block resumes on the same block. It is refused with
`window`, and where a solve happens every row (`solve_every=0` or
`max_rows_between_solves=1`), because the block is merged before each
solve and would never hold more than one row. `docs/PERFORMANCE.md` §18
has the table.

### `rls` — recursive least squares

*API:* [`po.spec.rls`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.rls) — *Rust:* [`rls.rs`](crates/online-core/src/rls.rs) — *Outputs:* [fields](docs/OUTPUTS.md#rls)

```
A ← λA + w zzᵀ       b_j ← λb_j + w y_j z        β_j = A⁻¹ b_j
A₀ = ridge·I         b₀ = ridge·coef_prior
```

Coefficients move every row with no solve staleness. The state is the
Cholesky factor of `A`, updated by Givens rotations. That is the square-root
form: `O(k²)` per row, the same as the textbook `P` recursion, without
either of its two failure modes. `P` accumulates rounding asymmetry by a
factor of `1/λ` per row, and one extreme row can cancel it and freeze a
coefficient for good. `ridge` sets
`A₀ = ridge·I` and, unlike `ewridge`, penalizes the intercept too.
Algebraically identical to `ewridge(ridge_decay=True)` solved every row; a
test holds them to <1e-9. A row with any null target is predict-only for all targets, since the
factor is shared.

### `lasso` — lasso path with free λ selection

*API:* [`po.spec.lasso`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.lasso) — *Rust:* [`lasso.rs`](crates/online-core/src/lasso.rs) — *Outputs:* [fields](docs/OUTPUTS.md#lasso)

Coordinate descent on the standardized statistics, warm-started along the
path and across solves:

```
ρ_i = c_i − Σ_{j≠i} C_ij β_j
β_i = soft(ρ_i, λ·l1_ratio) / (C_ii + λ(1 − l1_ratio))
```

`l1_ratio < 1` gives elastic net. Predictions for every path point are
computed anyway, so `lam_selected_<target>` — the argmin of an EW
out-of-sample squared error over the path — adds no work of its own. It is
reported as it stood *before* the row, like every other output.

### `kalman` — random-walk-β dynamic linear model

*API:* [`po.spec.kalman`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.kalman) — *Rust:* [`kalman.rs`](crates/online-core/src/kalman.rs) — *Outputs:* [fields](docs/OUTPUTS.md#kalman)

```
β_j ← Φβ_j    P_j ← ΦP_jΦ + Q·Δclock    Φ = diag(2^(−Δclock/r_i))
s   = zᵀP_j z + R_j/w                   k   = P_j z / s
β_j ← β_j + k(y_j − zᵀβ_j)              P_j ← P_j − k zᵀP_j
```

Process noise comes from a per-factor **coefficient halflife** on
standardized features, `q_i = σ²(ln2 / h_i)²`, matching the steady-state
gain of EW-RLS; `coef_halflife` may be a scalar or one value per slot, and
`inf` pins a coefficient. Observation noise defaults to the EW residual
variance. Standardization is internal and on by default. With
`standardize=False`, `q=0` and a fixed `obs_var`, this is exactly a Bayesian
linear regression (it reproduces river's `BayesianLinearRegression` to
3.6e-15).

**Reverting coefficients.** By default `Φ = I` and a coefficient is a random
walk: once a slope has been learned it stays until new rows move it. With
`revert_halflife`, a slope decays toward zero with the clock instead —
halved every `r_i` clock units while nothing is observed, and pulled back
toward zero by the same factor before each update. That is a mean-reverting
(AR(1)) prior: a regressor that is only occasionally active is forgotten
between its bursts rather than kept at its last value, and a stale effect
cannot persist through a run of null targets. The reversion acts in the
standardized coordinates, so "zero" means "no effect" for a slope and "the
target averages zero" for the intercept; a scalar applies to every slot
including the intercept, and a list with `inf` in the first slot exempts it.
The steady-state prior variance of a reverting slot is `q_i·Δclock/(1−φ_i²)`
instead of growing without bound. `predict` propagates the coefficients by
the same `Φ` over the distance from the last learned row (capped by
`max_dclock`), so a prediction far past the data is the intercept alone.

```python
revert = po.spec.kalman(
    "k",
    targets=["y"],
    features=["signal_a", "signal_b"],
    coef_halflife=100.0,
    revert_halflife=[float("inf"), 50.0, 50.0],  # intercept stays; slopes revert
    halflife=200.0,
    clock="t",
    max_dclock=10.0,
)
out = po.ModelBank([revert]).fit_predict(df)
```

### `huber` / `quantile` — robust regression

*API:* [`po.spec.huber`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.huber) and [`po.spec.quantile`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.quantile) — *Rust:* [`robust.rs`](crates/online-core/src/robust.rs) — *Outputs:* [huber](docs/OUTPUTS.md#huber), [quantile](docs/OUTPUTS.md#quantile)

IRLS reweighting on the ridge update, using each row's *prior* residual so
the reweighting stays out-of-sample. Huber: `w = min(1, δσ/|r|)`. Quantile at
level τ: the check loss's own IRLS weight, `2τσ/|r|` above the fit and
`2(1−τ)σ/|r|` below it, with `|r|` floored at `quantile_eps · σ` so a residual
near zero cannot blow the weight up. Weights are per target, so `S` is per
target here.

### `sgd` — stochastic gradient descent

*API:* [`po.spec.sgd`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.sgd) — *Rust:* [`sgd.rs`](crates/online-core/src/sgd.rs) — *Outputs:* [fields](docs/OUTPUTS.md#sgd)

```
eta = zᵀβ        p = link(eta)        gᵢ = (dL/d eta)·zᵢ·w + l2·βᵢ        βᵢ -= lrᵢ·gᵢ
```

| loss | link | `dL/d eta` |
|---|---|---|
| `squared` | identity | `p − y` |
| `huber` | identity | `clamp(p − y, ±delta)` |
| `quantile` | identity | `1{y < p} − τ` |
| `epsilon_insensitive` | identity | 0 inside the tube, else `sign(p − y)` |
| `poisson` | log | `p − y` |
| `logistic` | sigmoid | `p − y` |

O(k) per row and no solves — the cheap baseline, and the only model here that
takes count targets (`loss="poisson"`). The learning rate is `constant`,
`inv_scaling` (`lr/(1+n_eff)^power`) or `adagrad`, whose accumulator decays
on the clock so an adapted rate re-opens after a long gap. `clip_gradient`
defaults to `1e3`, because with a log link one large count makes the next
gradient exponentially bigger; it does not bind for identity-link losses.

**Constrained coefficients.** `coef_min` and `coef_max` bound each slope,
and `coef_sum` fixes their total. After every update the slopes are moved to
the nearest point that satisfies all three (the Euclidean projection); the
intercept is never constrained. A bound is a number for every slope or a
list with one entry per feature, and `inf` means no bound on that side.
Portfolio weights that must be long-only and fully invested are
`coef_min=0.0, coef_sum=1.0`; a sign the model must respect is
`coef_min=0.0` alone; a slope pinned at a known value is `coef_min` equal
to `coef_max`. The fit starts from the projected zero (the uniform weights,
on a simplex), and `coef` reports what the projection returned, in the
caller's units even under `scale_features=True`.

```python
weights = po.spec.sgd(
    "w",
    targets=["y"],
    features=["signal_a", "signal_b", "x0"],
    halflife=200.0,
    learning_rate=0.01,
    coef_min=0.0,
    coef_sum=1.0,
    coef_every=1,
)
fit = po.ModelBank([weights]).fit_predict(df)
last = fit["w"].struct.field("coef").drop_nulls()[-1]
assert min(last[1:]) >= 0.0 and abs(sum(last[1:]) - 1.0) < 1e-12
```

### `pa` — passive-aggressive regression

*API:* [`po.spec.pa`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.pa) — *Rust:* [`pa.rs`](crates/online-core/src/pa.rs) — *Outputs:* [fields](docs/OUTPUTS.md#pa)

```
loss = max(0, |y − p| − eps)      s = ‖z‖²
pa    τ = loss / s          pa1  τ = min(c, loss/s)      pa2  τ = loss / (s + 1/(2c))
β    += τ · sign(y − p) · z
```

Each row poses a constraint and the update is the smallest change that
satisfies it — no learning rate to tune. Plain `pa` moves the fit as far as
one bad row demands, so `pa1` is the default. PA keeps no accumulators, so
its coefficients have no halflife; the clock only drives `n_eff`. A row
weight below 1 scales `τ`; above 1 it counts as 1.

`pa` takes the same `coef_min`, `coef_max` and `coef_sum` as `sgd`, with
the projection applied after each update. The step then no longer meets
the row's margin exactly, and a truth outside the set is never reached, so
keep `c` small: each row moves the fit only as far as `c` allows and the
projection takes the rest back.

### `ew_cov` — exponentially weighted moments

*API:* [`po.spec.ew_cov`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.ew_cov) — *Rust:* [`ewcov.rs`](crates/online-core/src/ewcov.rs) — *Outputs:* [fields](docs/OUTPUTS.md#ew_cov)

```
W'   = λW + w        m'ᵢ = (λW·mᵢ + w·xᵢ) / W'      S'ᵢⱼ = (λW·Sᵢⱼ + w·xᵢxⱼ) / W'
varᵢ = Sᵢᵢ − mᵢ²     covᵢⱼ = Sᵢⱼ − mᵢmⱼ             corrᵢⱼ = covᵢⱼ / √(varᵢ·varⱼ)
```

Running mean, variance, std, covariance and correlation of the columns you
name, on the same clock as every model here — one O(k²) update per row,
where a pure-Polars pairwise EW correlation needs O(k²) *passes*. With
`precision_prior` set it also gives `partial_corr`, the correlation between
two columns controlling for all the others, read off `(C + s·prior·I)⁻¹`
(O(k³), paid only when asked for). Values are read from the state before each
row, so an `ew_cov` output can be a feature for that same row without leaking
it.

**`window` — a hard cutoff, not a softer decay.** A halflife of `h` never
forgets: three halflives back still carries 12.5% of the weight. `window=w`
makes that exactly zero — a row older than `w` clock units contributes
nothing:

```
weight(age) = 0.5 ** (age / halflife)   if age <= window
            = 0                          otherwise
```

Note the first line: *inside* the window the weights are still exponential,
so this is not a rolling flat mean and the newest row still dominates. It is
exact, not approximate, because an EW sum contains its own past — everything
at or before a time `u` is `λ^(t−u)` times the accumulator as it stood then,
so subtracting that leaves precisely the rest. The model keeps a ring of
snapshots to do it, which is the one place here where memory grows with a
*window* rather than with the state: about 3 MB per group for a 1,000-row
window over 20 columns, divided by `window_every` if you snapshot less often.

Four things worth knowing before you read the numbers. The guarantee is
**one-sided**: the boundary is the oldest snapshot still inside the window,
so a coarse `window_every` discards a little more than asked, never less.
The clock is the **decayed** one, after `max_dclock` and any `session_gap`.
The edge is a **discontinuity** — a row ageing out drops its whole weight at
once, so the series has small steps an EWMA does not. And it is a
**subtraction**, so precision falls with the fraction discarded: negligible
at `window = 3h`, worse as the window shortens toward the halflife.
`n_eff` becomes the weight inside the window, so `min_periods` now gates on
something that stops growing, and a clock gap longer than `window` empties
it and reports nulls rather than stale numbers.

If your data fits in memory and you only want moments, polars already does
this: `df.rolling("t", period="3h").agg(...)` with an exponential weight is the
same number to 1e-14. The reason to reach for the spec is a stream, a saved
state, or the time: the rolling window recomputes each window at `O(n·W)` where
this is `O(n)`. At 200k rows and a 4,680-row window that is 15.8 s against 14
ms.

Three more outputs come off the same state. `"mahal"` in `stats` adds
`mahal`, the Mahalanobis distance of the row from the running mean,
`√(δᵀ (C + s·prior·I)⁻¹ δ)` with `δ = x − m` — how far the row is from what
the columns have been doing *together*, in standard deviations. A row with
every column in range but in a combination never seen before scores high here
and nowhere else; on Gaussian columns `mahal²` is χ² with k degrees of
freedom, and with one column it is `|z|`. It needs `precision_prior`, and the
prior fades like `partial_corr`'s. `mahal_quantiles=[0.99]` adds
`mahal_q0.99`, a P² quantile of the scores so far, so `mahal > mahal_q0.99`
is a distribution-free "one row in a hundred". `pca=r` adds the top `r`
eigenpairs of the covariance: `pc<j>_var`, `pc<j>_share` of the trace, the
loading on each feature `pc<j>_<feature>`, and the row's score on that
component `pc<j>_score`. The eigendecomposition costs O(k³), so
`pca_every=n` refreshes it every `n` rows and scores the rows in between on
the last loadings; each refresh keeps the sign of the previous one, so a
loading never flips between rows.

`stats=[]` is legal and means accumulate only. The spec learns the same
moments and emits nothing but `n_eff`; its value is its state, read back
with `bank.gram("mv")` or `bank.describe()`. That is the form for a wide
set of columns, where emitting even the means is k values per row nobody
reads. `pca` and `mahal_quantiles` stand without a statistic in the list.

```python
mv = po.spec.ew_cov("mv", features=["x0", "x1", "x2"], clock="t", max_dclock=300.0,
                    halflife=500.0, stats=["mahal"], precision_prior=1e-6,
                    mahal_quantiles=[0.99], pca=1, pca_every=20)
scores = po.ModelBank([mv]).fit_predict(df).unnest("mv")
odd = scores.filter(pl.col("mahal") > pl.col("mahal_q0.99"))   # the joint outliers
first = scores.select("pc0_share", "pc0_x0", "pc0_x1", "pc0_x2", "pc0_score")
```

`lags` accumulates the same co-moments one step further out: how column `a`
now moves with column `b` `ℓ` rows ago. With `W` and `m` the weight and mean
before the row, and both deviations taken against that mean,

```
C_ℓ' = a·C_ℓ + a·b·(x_t − m)(x_{t−ℓ} − m)'
```

with the same `a` and `b` the co-moments use — so lag 0 would be `comoments`
exactly. Lags count **learned rows within the group**, not clock units, and
must be strictly increasing and ≥ 1. Read them from `bank.gram("mv")` as
`lags` and `lag_comoments` (an `(L, k, k)` array), or add `"lagcorr"` to
`stats` to emit `lagcorr_<a>_<b>_l<ℓ>` per lag and **ordered** pair — both
orientations, because a lagged matrix is not symmetric: `a` leading `b` is
not `b` leading `a`.

```python
lagged = po.spec.ew_cov("lagged", features=["x0", "x1"], halflife=500.0,
                        stats=["corr", "lagcorr"], lags=[1, 5])
lead = po.ModelBank([lagged]).fit_predict(df).unnest("lagged")
```

The ring of past rows is emptied on a session change and on a clock gap
beyond `max_dclock` — the two events after which "the row `ℓ` back" no longer
means a row `ℓ` ago — and a zero-weight row ages the matrices without
entering it. Nothing else moves: clearing the ring is not a reset.

### `ftrl` — online logistic regression

*API:* [`po.spec.ftrl`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.ftrl) — *Rust:* [`ftrl.rs`](crates/online-core/src/ftrl.rs) — *Outputs:* [fields](docs/OUTPUTS.md#ftrl)

FTRL-proximal (McMahan et al. 2013) for binary targets, with the accumulators
decayed on the same clock as everything else:

```
β_i = 0 if |z_i| ≤ l1 else −(z_i − sgn(z_i)l1) / ((β + √n_i)/α + l2)
p   = sigmoid(zᵀβ)     g_i = (p − y)·z_i·w
z_i += g_i − ((√(n_i + g_i²) − √n_i)/α)·β_i      n_i += g_i²
```

With `loss="logistic"` (default) `pred` is a probability and `resid = y − p`;
with `loss="squared"` it is the linear prediction — sparse linear regression
with no solves, and L1 support, which `ewridge` does not have.

### `holt` — Holt's linear trend

*API:* [`po.spec.holt`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.holt) — *Rust:* [`holt.rs`](crates/online-core/src/holt.rs) — *Outputs:* [fields](docs/OUTPUTS.md#holt)

The one model that takes no features: it extrapolates the target's own level
and trend.

```
pred     = l + b·Δt
l' = α·y + (1−α)·pred        b' = β·(l' − l)/Δt + (1−β)·b
```

`α` and `β` come from `level_halflife` and `trend_halflife` in clock units;
the trend is per clock unit, so an irregular clock extrapolates the right
distance. `coef` is `[level, trend]` per target; `trend_halflife=inf` pins
the trend at zero, leaving a plain EW level. There is no seasonal term,
because a seasonal index is a `group` on the phase, which the bank already
does. Run it in the same bank as the real model to answer "how much is the
regression actually adding?" — compare `sigma`, or let `emit_selected`
choose.

```python
po.spec.holt("baseline", targets=["y"], clock="t", max_dclock=600.0,
             level_halflife=200.0, trend_halflife=2000.0)
```

### `kmeans` — exponentially weighted k-means

*API:* [`po.spec.kmeans`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.kmeans) — *Rust:* [`cluster/kmeans.rs`](crates/online-core/src/cluster/kmeans.rs) — *Outputs:* [fields](docs/OUTPUTS.md#kmeans)

The one model with no target: it labels each row with the nearest of `k`
centres, read before the row is learned, so the label is out-of-sample like
every prediction here.

```
j*   = argmin_j ‖x − c_j‖²          distances in units of each feature's EW sd
n'_j = λn_j + w                      c'_j = c_j + (w/n'_j)(x − c_j)     for j = j*
```

Each centre is the EW mean of the rows assigned to it — `ew_cov`'s mean
recursion, per cluster. The struct holds `cluster`, `dist` (to the centre),
`dist2` (to the runner-up), `n_eff`, and `coef` = the centres, `k` rows of
`len(features)`.

```python
km = po.spec.kmeans("km", features=["x0", "x1", "x2"], k=3, clock="t",
                    halflife=2000.0, max_dclock=300.0, warm_rows=100)
out = po.ModelBank([km]).fit_predict(df).unnest("km")
po.spec.coef_index(km)        # target = "cluster0".., term = the feature
```

Seeding waits for `warm_rows` rows (default 500), then places the centres
with `seed_rule="lloyd"`: the best of ten k-means++ starts by inertia. One
start lands in the wrong partition a third of the time on five blobs in four
dimensions; the restarts tell them apart. The rows are replayed and the
buffer freed, so the model is O(k·p) from then on.

**What split–merge repairs.** A row far outside its cluster (about four
standard deviations of `dist²` above the typical radius) is scored but not
learned: it is summarised. Every `split_merge_every` rows the two closest centres are
compared, and if they are closer than `split_merge` times the sum of their
radii — two centres in one blob — one is freed and placed on the far rows,
provided enough have gathered to be a cluster's worth. A centre whose blob
vanished decays; once under `dead_frac` of an equal share it is re-placed
the same way. That takes `log2(1/dead_frac)` halflives: 4.3 at the default
0.05, 2 at 0.25. Raise `dead_frac` when regimes change faster than that;
the price is that a cluster lighter than `dead_frac/k` of the stream loses
its centre whenever any row is far. What the move cannot see is one centre
owning two blobs, whose rows are all within its own radius — seeding with
`lloyd` is what prevents it. Set `split_merge=0` for plain sequential
k-means.

### `micro` — density-based clustering, any shape

*API:* [`po.spec.micro`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.micro) — *Rust:* [`cluster/micro.rs`](crates/online-core/src/cluster/micro.rs) — *Outputs:* [fields](docs/OUTPUTS.md#micro)

`kmeans` needs `k` and finds round clusters. `micro` finds clusters of any
shape, does not need their number, flags the rows that belong to none, and
follows clusters that appear and vanish. It is DenStream's micro-clusters
with a linkage step over them.

A summary is a small cluster: a decayed weight `n`, a centre `c` and a
radius `r`, the EW root-mean-square distance of its rows from the centre.
Each row goes to the nearest summary that can take it without its radius
passing `eps`, in units of each feature's EW sd. If none can, the row opens
one. A summary with `n ≥ beta_mu` is established. Every `prune_every` rows
the light summaries are dropped and the established ones are linked:
centres within `L` of each other share a label, and `L` is read from the
spacing the summaries already show unless `macro_link` sets it.

```
n_j  ← λ n_j                                               every summary
j*   = nearest summary that keeps  a r²_j + a b ‖x − c_j‖² ≤ eps² p,
       a = n_j/(n_j + 1),  b = 1/(n_j + 1);  else a new one at x
n_j* ← n_j* + w     c_j* ← c_j* + (w/n_j*)(x − c_j*)     r²_j* ← min(·, eps² p)
```

The struct holds:

| field | meaning |
|---|---|
| `cluster` | the label of the nearest established summary; null while there is none |
| `dist` | the distance to that summary's centre |
| `micro` | the id of the summary this row goes to |
| `outlier` | no established summary takes the row |
| `n_clusters`, `n_micro` | how many of each the state holds |
| `n_eff` | the weight behind the state before this row |
| `coef` | the established summaries, one `[id, label, n, radius, c_1 … c_p]` row each |

All are read before the row is learned. Ids are monotone and never reused; a
label is the smallest id in its chain, so it outlives everything but that
summary.

```python
mc = po.spec.micro("mc", features=["x0", "x1"], eps=0.1, clock="t",
                   halflife=2000.0, max_dclock=300.0, min_periods=50.0)
out = po.ModelBank([mc]).fit_predict(df).unnest("mc")
out.select("cluster", "outlier", "n_clusters", "n_micro").tail(3)
```

**Choosing `eps`.** It is the spread the model should read as *one*
cluster, per standardized coordinate: about 0.07 for two-dimensional
shapes, 0.3 for well-separated Gaussians in twenty dimensions. Both ways to
get it wrong show in the outputs. If nearly every row is an `outlier` and
`cluster` stays null, `eps` is too small: no summary reaches `beta_mu`
before it is pruned. If `n_micro` is about the number of clusters, `eps` is
too coarse: each cluster is one summary, so the derived `L` reads the
spacing *between* clusters and bridges them into one. Lower `eps`, or set
`macro_link=2` to link only summaries that touch.

Measured at 20k rows with the `eps` above: moons, rings and five Gaussians
in twenty dimensions all score ARI 1.000 against the truth, where `kmeans`
cannot follow the first two. Noise drawn uniformly over the box is flagged
`outlier` 94% of the time, real rows 0.3%. A cluster born mid-stream has a
label within 200 rows; one whose rows stop lingers `halflife · log2(n /
beta_mu)`, with `n` the weight it had.

### `ew_class` — Gaussian classification on `ew_cov` moments

*API:* [`po.spec.ew_class`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.ew_class) — *Rust:* [`ewclass.rs`](crates/online-core/src/ewclass.rs) — *Outputs:* [fields](docs/OUTPUTS.md#ew_class)

A label column in place of a numeric target. The model keeps one `ew_cov`
state per class — a weight `n_c`, a mean `μ_c` and a centered covariance
`C_c` — and scores a row by Bayes' rule over Gaussian classes. `covariance`
picks the shape. `"full"` gives each class its own covariance: QDA.
`"shared"` pools them, weighted by the class weights: LDA. `"diagonal"`
keeps only the variances: Gaussian naive Bayes. `precision_prior` is the
ridge that makes a class scoreable from its first row, and it fades the way
`ew_cov`'s does.

```
π_c = n_c / Σ n         r_c = precision_prior · s_c        (s_c: the prior's fade)
M_c = C_c + r_c I  (full)      M = Σ π_c M_c  (shared)      diag(C_c) + r_c  (diagonal)
ℓ_c = ln π_c − ½ ln det M_c − ½ (x − μ_c)ᵀ M_c⁻¹ (x − μ_c)
p_c = exp(ℓ_c − max ℓ) / Σ exp(ℓ − max ℓ)                  class = argmax ℓ
n_c ← λ n_c + w·[y = c]        μ_c, C_c ← weighted Welford on the row's own class
```

The struct holds `class`, the most probable class as a string; one
`p_<class>` per declared class; `n_eff`; and `coef`, the class means in the
order of `classes` (`coef_up_x0` after `unnest`). All are read before the row
is learned, so a row's posterior never saw its own label. A class no row has
carried yet has `p = 0` exactly and null means. A null label scores the row
and learns nothing from it — so a stream whose labels arrive late is scored
by nulling the label and keeping the features. A label the spec does not list
is an error naming the row, the value and the classes. Integer and boolean
columns work as labels through their text: `classes=["0", "1"]`,
`classes=["true", "false"]`.

```python
labelled = df.with_columns(
    pl.when(pl.col("y") > 0).then(pl.lit("up")).otherwise(pl.lit("down")).alias("dir")
)
cl = po.spec.ew_class("cl", features=["x0", "x1", "x2"], label="dir", classes=["down", "up"],
                      covariance="shared", precision_prior=0.1, clock="t",
                      halflife=200.0, max_dclock=300.0, min_periods=20.0)
out = po.ModelBank([cl]).fit_predict(labelled).unnest("cl")
out.select("dir", "class", "p_up", "n_eff").tail(3)
```

**Choosing the shape.** `"full"` is the general case and costs one `k×k`
Cholesky per class per row. `"shared"` factorizes once per row, and is the
right model when the classes differ in location but not in spread — it then
matches `"full"` to a fraction of a percent on the test data, with fewer
parameters to learn. `"diagonal"` is the cheapest and cannot see a
correlation: two classes with the same marginals and opposite correlations
are one class to it. Measured at 400k rows, six features and three classes:
0.9M rows/s full, 1.8M shared, 5M diagonal. On three Gaussian classes with
their own covariances the accuracy sits within 0.001 of the Bayes rate the
generating parameters allow, and the posteriors are calibrated to about 0.01.

### `seqtest` — a sequential test of a sign, by betting

*API:* [`po.spec.seqtest`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.seqtest) — *Rust:* [`seqtest.rs`](crates/online-core/src/seqtest.rs) — *Outputs:* [fields](docs/OUTPUTS.md#seqtest)

Not a regression. A `seqtest` asks whether a column tends to be positive; with
`a` and `b`, it asks instead whether one spec of the bank predicts closer than
another. Either way the answer is evidence you can read at any row, as often as
you like, and act on the first time it is enough. A p-value cannot be used that
way, because peeking at one inflates its error rate. An e-process can, and that
is the whole reason to reach for it. Per target it keeps the wealth of two
gamblers, one betting that the next sign is positive and one that it is
negative. Each stakes the Krichevsky–Trofimov fraction set by the counts so
far, and never bets against its own lead:

```
s = sign(y)                    n⁺, n⁻: the signs counted before this row,  n = n⁺ + n⁻
λ⁺ = max(0, (n⁺ − n⁻) / (n + 1))          λ⁻ = max(0, (n⁻ − n⁺) / (n + 1))
ln E⁺ ← ln E⁺ + ln(1 + λ⁺ s)              ln E⁻ ← ln E⁻ + ln(1 − λ⁻ s)
```

Under the null — given everything so far, the next sign is no more likely
positive than negative — `E⁺` is a nonnegative supermartingale, and Ville's
inequality gives `P(E⁺ ever reaches 1/α) ≤ α`. So `log_e_pos ≥ ln 20`
rejects at the 5% level however many times you looked, and however the
rows depend on each other. No distribution is assumed and the size of the
values is invisible: 60% small gains and 40% huge losses is "positive". Where
the clip never binds the wealth has the closed form `2ⁿ B(n⁺+½, n⁻+½) / π`,
the Beta(½, ½) mixture, and the tests hold the bank to it; the two sides'
average is an e-value for the two-sided question. The struct holds
`log_e_pos_<t>`, `log_e_neg_<t>`, `n_pos_<t>`, `n_neg_<t>` and `n_eff`, all
as they stood **before** the row. A zero, a null or a NaN is a tie: it bets
nothing and counts nothing. A trial is a row, so there is no `weight` and no
`halflife` — a spec that gives them is refused — and `session` or
`on_clock_reset="reset_state"` restarts the test.

With `a` and `b` the test compares two specs of the bank: the sign tested is
`|resid_b| − |resid_a|`, positive when `a` came closer, and the fields are
`log_e_a_<t>`, `log_e_b_<t>`, `wins_a_<t>`, `wins_b_<t>`. The bank runs a
comparison after the specs it names, on the out-of-sample residuals their
structs report; a row where either side is null — warm-up, a skipped row —
is no trial. `a_suffix` and `b_suffix` pick a grid instance (`"@h500"`,
`"__r0.5@h500"`). A comparison inside a bank is chunk-invariant, saved with
the state and streams like everything else; [`po.eval.seqtest`](https://hgilde.github.io/polars-online/eval.html#polars_online.eval.seqtest) is the same
computation over a frame you already have.

```python
common = dict(targets=["y"], features=["x0", "x1"], clock="t", max_dclock=300.0, group="bond_id")
ridge = po.spec.ewridge("ridge", halflife=500.0, **common)
kalman = po.spec.kalman("kalman", halflife=500.0, coef_halflife=100.0, **common)
closer = po.spec.seqtest("closer", targets=["y"], a="kalman", b="ridge", group="bond_id")
out = po.ModelBank([ridge, kalman, closer]).fit_predict(df)
verdict = out.group_by("bond_id").agg(pl.col("closer").struct.field("log_e_a_y").max())
# log_e_a_y >= ln(20): on that bond, kalman beat ridge at the 5% level, read at any row
```

### `marginal` — every pair's moments, kept in the state

*API:* [`po.spec.marginal`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.marginal) — *Rust:* [`marginal.rs`](crates/online-core/src/marginal.rs) — *Outputs:* [fields](docs/OUTPUTS.md#marginal)

A `marginal` is not a regression and not a joint fit. It keeps the
exponentially weighted moments of each (feature, target) pair on its own,
as if every pair were a two-column `ew_cov`. For `p` features and `T`
targets that is O(p·T) per row; one `ew_cov` over all the columns would be
O((p + T)²).

Per target `t`, on a row where `y_t` is present, with `W_t` the weight
behind that target before the row:

```
W'_t = λW_t + w        a = λW_t / W'_t        b = w / W'_t        Q'_t = λ²Q_t + w²
S'_yy = a·S_yy + a·b·(y_t − m_y)²             S'_xx = a·S_xx + a·b·(x_j − m_x)²
S'_xy = a·S_xy + a·b·(x_j − m_x)(y_t − m_y)    m' = m + b·(value − m)
```

That is `ew_cov`'s arithmetic. A pair's correlation is the one an `ew_cov`
over the two columns would report, to the bit. A null target ages its own
pairs (`W_t ← λW_t`) and learns nothing for them. A null feature skips the
row, as everywhere. Weights, the clock, sessions and groups apply as they
do to every model.

Nothing is emitted per row but `n_eff`. The pairs are the state, and
`bank.marginal("pairs")` reads them as a long frame with one row per
(group, instance, feature, target):

| column | meaning |
|---|---|
| `n_eff` | the target's `W_t`, the weight behind its pairs |
| `n_kish` | `W_t² / Q_t`, Kish's effective sample size: the count of equally weighted rows that carry the same information (`(1 + λ)/(1 − λ)` in the limit for unit weights) |
| `mean_x`, `var_x`, `mean_y`, `var_y`, `cov` | the pair's moments, population form |
| `corr` | `cov / √(var_x · var_y)` |
| `beta` | `cov / var_x`, the slope of the target on that feature alone |
| `t` | `corr·√((n_kish − 2)/(1 − corr²))`, the t-statistic of the correlation at the Kish sample size |

`t` is a scale for comparing pairs, not a p-value: the rows are neither
independent nor Gaussian. `corr`, `beta` and `t` are null until the
target's `W_t` reaches `min_periods` (default 3; two rows give ±1 whatever
the data), and where they are undefined — a constant feature, or
`n_kish ≤ 2` for `t`. A bank loaded from a file reports the pairs the bank
that saved it would. One chunk or a thousand gives the same frame to the
bit.

Two views sit on top of that, both off unless asked for.

**`lags=[1, 2, 3, 5, 8]` — is `t` telling the truth?** `t` is built on
`n_kish`, which is the right count for unequal weights and says nothing about
serial dependence. On a smooth stream consecutive rows are nearly the same
observation, so `t` claims evidence that is not there. Lags fix that. The pair's
moments are accumulated at each lag too — the same statistic `ew_cov(lags=)`
computes, to the bit — and `serial_rule` turns them into Bartlett's correction.

| column | meaning |
|---|---|
| `lagcorr_xx`, `lagcorr_yy` | each series' own autocorrelation, one entry per lag |
| `lagcorr_xy`, `lagcorr_yx` | the feature now against the target `ℓ` rows back, and the reverse |
| `n_serial` | `n_kish` divided by `1 + 2·Σ ρ_x(ℓ)·ρ_y(ℓ)` (Bartlett 1935) |
| `t_serial` | the same statistic as `t`, against that count |
| `phi_x`, `phi_y` | the fitted per-row decays, under `serial_rule="geometric"` |

The four lists are `ew_cov`'s `lagcorr` numbers exactly — the lagged
covariance over the two standard deviations, not clamped to `[−1, 1]`, since a
lagged correlation is not bounded by one in a finite sample. A row where the
target is missing ages its weight and holds the lag moments, as it holds the
pair's.

Two *independent* AR(1) series with `φ = 0.9` and `0.8` come out at `t = 2.39`
and `t_serial = 1.03`. The first is a finding; the second is the truth.

`lagcorr_xy` against `lagcorr_yx` is worth having on its own. A feature whose
`lagcorr_yx[0]` beats its `corr` *leads* its target. One whose `lagcorr_xy[0]`
does *follows* it, which is usually a column sampled late.

**`bins=16` — what a correlation cannot see.** Everything above is linear. A
feature can be strongly related to a target with `corr` at zero: a threshold, a
V, a saturation. Bin the feature and keep the target's moments inside each bin,
and all three become visible.

| column | meaning |
|---|---|
| `bin_edges` | the feature's edges, fixed once and never moved |
| `bin_n`, `bin_mean_y`, `bin_var_y` | the target's weight, mean and variance in each bin — the response curve |
| `split_gain` | the fraction of the target's variance removed by the best single cut |
| `split_at` | where that cut falls, in the feature's units |
| `split_gain_t` | the `t` a `corr` would need to match that gain |

`split_gain` is a regression stump's `R²`, so it compares directly with
`corr²` and the difference is the nonlinear surplus. It costs `O(bins)` of
state per pair and one binary search per pair per row, which is why it can run
across ten thousand columns in the pass that gives them `corr`.

Read `split_gain_t` as a ranking, not a p-value. The cut was chosen by
maximising over the candidates, and the statistic does not know that.

Give the edges outright with `bin_edges` — a list per feature, or a dict keyed
by name — and they are exact and comparable across runs; `bins`, `bin_rule` and
`bin_warm_rows` describe learning them and are refused beside it. Otherwise
they are learned from the first `bin_warm_rows` rows (default 1,000), by
weighted quantile or by equal width. A value that carries more than a bin's
share — an indicator's zero — fills a bin of its own and the rest share what is
left, so the 5% of rows that carry the signal are not lost among the zeros.
Those rows are held and replayed, not spent: the histogram is what it would
have been had the edges been known before the first row. A feature keeps only
the bins it can support, so a binary feature has two and a constant one has a
single bin and no split. Each bin's moments are kept the way every accumulator
here is kept, so a target at `1e7` keeps its variance.

Both views ride into `bank.closed_groups()` as `pair_*` columns beside the
others: `pair_split_gain` as a list over the pairs, `pair_lagcorr_xx` and
`pair_bin_n` as lists of lists.

```python
pairs = po.spec.marginal("pairs", targets=["y", "ret"],
                         features=["x0", "x1", "x2", "signal_a", "signal_b"],
                         clock="t", max_dclock=300.0, halflife=500.0, group="bond_id")
bank = po.ModelBank([pairs])
bank.fit_predict(df)                       # the struct holds n_eff alone
table = bank.marginal("pairs")             # group, instance, feature, target, n_eff, n_kish, ..., corr, beta, t
one_bond = bank.marginal("pairs", group="b0")   # 10 rows: five features by two targets

honest = po.spec.marginal("pairs", targets=["y"], features=["x0", "x1"],
                          halflife=500.0,
                          lags=[1, 2, 3, 5, 8], serial_rule="geometric",
                          bins=16)                 # + lagcorr_*, n_serial, split_gain, ...
```

### `corrchange` — has the correlation structure changed?

*API:* [`po.spec.corrchange`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.corrchange) — *Rust:* [`corrchange.rs`](crates/online-core/src/corrchange.rs) — *Outputs:* [fields](docs/OUTPUTS.md#corrchange)

Two tests, because there are two questions.

`kind="monitor"` is the **closed-sample** constancy test of Wied, Krämer and
Dehling (2012), run over consecutive spans of `span_rows` rows. At the last row
of a span, per pair:

```
Q = max_{2≤j≤T} (j/√T)·|ρ̂_j − ρ̂_T| / D̂
```

with `ρ̂_j` the correlation of the span's first `j` rows and `D̂` the
delta-method long-run standard deviation of `ρ̂`. Under the null `Q` converges
to `sup|B|`, a Brownian bridge, so the critical value is the Kolmogorov
quantile — computed from the series, not pinned, and it reproduces the
published 1.3581 at 5%. That published null is the point: the test's size and
power are held to the paper's own tables (`.035` at ρ = 0 and `T = 500`,
`.587` power on a `0.5 → 0.7` break), not to numbers this implementation
happened to produce.

```python
c = po.spec.corrchange("break", features=["x0", "x1"], span_rows=500)
out = df.online.fit_predict([c]).unnest("break")   # stat, crit, flag, since_flag
```

The price is a delay of at most `span_rows` rows: nothing is reported until a
span closes. The paper's own *sequential* form, with a boundary function, is
Wied and Galeano (2013), which has not been read here.

`scalar=True` runs the same CUSUM on the **equicorrelation** of the
standardised row (`deco`'s `u`) — one statistic however many columns, and a
test of its level rather than of a pair.

`kind="window"` asks how *big* the change is instead: `‖vech(R̂_pre − R̂_post)‖`
over two adjacent windows, against a fixed `crit` or a **permutation**
quantile — `n_perm` shuffles of the pooled rows between the windows, in
blocks of `perm_block` so serial dependence does not make the null too
liberal. Not a sign-flip null, which a first reading of the literature
suggests: negating a whole row leaves every correlation exactly where it was.

The flag rate per row is not `alpha` for the window kind: two windows that
slide by one row are almost the same windows, so a statistic above the
quantile stays above it for a run of rows.

### `hmm` — which regime are we in

*API:* [`po.spec.hmm`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.hmm) — *Rust:* [`hmm.rs`](crates/online-core/src/hmm.rs) — *Outputs:* [fields](docs/OUTPUTS.md#hmm)

`ew_class` classifies a row against *labelled* Gaussians. An `hmm` does the
same arithmetic with no labels: the state is hidden, and a transition matrix
carries information from one row to the next. That is the difference between
"which regime does this row look like" and "which regime are we in", and the
second is usually the question.

Hamilton's filter, one row at a time, from the filtered `p` the previous row
left:

```
p1_l   = Σ_k p_k·Π_kl                      the predicted state
f_l    = N(x | μ_l, Σ_l + r_l·I)           the state's density
loglik = ln Σ_l p1_l·f_l                   the row's surprise
p_l   ← p1_l·f_l / Σ                       the filtered state
```

Everything reported is read before the row is learned from. Each state's
accumulator then takes the row at weight `w·p_l`. The responsibilities sum
to `w`, so `n_eff` is the shared recursion untouched — a row splits across
the states rather than counting more than once. The transition matrix is
learned from the **filtered joint of consecutive states**,
`ξ_kl = p_k(t−1)·Π_kl·f_l / Σ`, with a Dirichlet pseudo-count keeping a
never-visited row a distribution.

```python
h = po.spec.hmm("regime", features=["x0", "x1"], k=2, precision_prior=1e-2,
                halflife=500.0, warm_rows=400)
out = df.online.fit_predict([h]).unnest("regime")   # p_0, p_1, p1_0, p1_1, state, loglik
```

What the transition chain adds, measured: on two-dimensional blobs 1.5
apart, a memoryless nearest-centre rule *given the true centres* is 85% right
and the filter is 99%.

`precision_prior` is required — a state's centred co-moments start at zero,
and a zero matrix has no density. Give `means` and `covs` to filter with
known states, and `learn=False` to freeze them. Otherwise it seeds from the
first `warm_rows` learned rows with `kmeans`' rule, and every output is null
until then. **`warm_rows` should span more than one regime**, or the seeds
are two halves of one. `exog_tvtp` drives the matrix from a column instead,
through fixed `tvtp_coef`.

One limitation worth knowing: a single extreme row can be captured by one
state, and in mean form a state with zero responsibility keeps its moments —
so a state that stops winning never forgets, and the mixture is left short
one state. A larger `precision_prior`, given states, or cleaning upstream
are the mitigations.

### `rcov` — a block's realised covariance, robust to noise

*API:* [`po.spec.rcov`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.rcov) — *Rust:* [`rcov.rs`](crates/online-core/src/rcov.rs) — *Outputs:* [fields](docs/OUTPUTS.md#rcov)

A realised covariance over ticks is the sum of outer products of returns.
Over real tick data it is wrong twice: each price is the efficient one plus
a measurement error, and the error's variance accumulates with every tick;
and if the series are not observed together, the correlation is attenuated
towards zero. Both are estimated away by published estimators that are sums
over lags — which is exactly what a stream can accumulate.

`rcov` has no decay and no per-row output but `n_eff`. Its value is the
block, emitted when the group closes, so it needs `group` and `group_close`
and the estimate rides in that row.

```python
r = po.spec.rcov("rk", features=["x0", "x1"], kind="kernel",
                 group="block", group_close="monotone", block_rows=2000)
bank = po.ModelBank([r])
bank.fit_predict(by_block.select("x0", "x1", "block"))
blocks = bank.closed_groups()      # rcov, rcorr, rcov_n, bandwidth_used, ...
```

Rows are **returns**: difference upstream. Three kinds:

| `kind` | what it is |
|---|---|
| `plain` | `Σ x x'`. Equals `n` × an `ew_cov(lam=1)`'s uncentred second moment at close, to the bit — the cross-check, and the reference the other two are measured against |
| `kernel` | the multivariate realised kernel (Barndorff-Nielsen, Hansen, Lunde & Shephard 2011), `Σ_h k(h/(H+1))·Γ̂_h` with Parzen weights and jittered end points |
| `preavg` | the modulated realised covariance (Christensen, Kinnebrock & Podolskij 2010): returns pre-averaged over `k_n = ⌊θ√n⌋` with `g(x) = min(x, 1−x)`, less the residual bias |

Parzen is the only kernel: the Bartlett kernel is not consistent for this
estimator, and Parzen's 0.97 efficiency beats the quadratic spectral's 0.93.
`bandwidth` is a fixed `H`; left out it is `H = ⌈c*·ξ̂^{4/5}·n^{3/5}⌉` with
`c* = 3.5134`, which needs `block_rows` — the ring has to be sized before the
first row and `n` is known only at the close. `block_rows` is a sizing hint, not
a limit: a longer block runs, clipped, and reports `bandwidth_used`.

The closed row carries `rcov` and `rcorr` (`vech` of the upper triangle),
`rcov_n`, `rcov_kind`, `bandwidth_used`, `omega2` and `iv_sparse` (the noise
variance and sparse integrated variance behind the bandwidth), `iq` (a
realised-quarticity proxy, labelled one) and `psd_repaired`. A block too
short to estimate from gives nulls, not an error.

Nothing reads a future row: the jittered *end* point is formed at close from
observations already in state, and a product enters `Γ̂_h` only once both
legs are final. `weight` is taken as 0 or 1 only — a sum over returns has no
fractional row.

### `deco` — one correlation for the whole matrix

*API:* [`po.spec.deco`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.deco) — *Rust:* [`deco.rs`](crates/online-core/src/deco.rs) — *Outputs:* [fields](docs/OUTPUTS.md#deco)

A correlation matrix of `m` series has `m(m−1)/2` free entries. A stream
cannot keep them all moving without O(m²) work a row, and most of them are
estimated from too little data to be worth moving. `deco` (Engle & Kelly
2012) replaces them with their average and estimates that, in O(m) a row.

The row is standardised against the pre-row means and variances of an
`EwDiag`, `r_i = (x_i − m_i)/√v_i`. With `S₁ = Σ r_i` and `S₂ = Σ r_i²` over
`n` features, the row's estimate is their Lemma 2.3:

```
u = (S₁² − S₂) / ((n − 1)·S₂)        = mean of r_i·r_j over i ≠ j, over mean r_i²
```

and the level follows one of two dynamics, on the model's own clock:

```
"ew":      W' = λW + w,  b = w/W'      ρ' = ρ + b·(u − ρ)
"linear":  ρ' = (1 − α − β)·ρ̄' + α·u + β·ρ      (ρ̄ the "ew" level)
```

`"ew"` is `EwCov`'s mean form, so `rho` is exactly what an
`ew_cov(stats=["mean"])` over the `u` sequence would report. `"linear"` is
the paper's eq. 21 with correlation targeting, which needs `alpha` and
`beta` with `alpha + beta < 1`.

| output | meaning |
|---|---|
| `u` | this row's own estimate, read before the row is learned from |
| `rho` | the level as it stood before the row |
| `loglik` | the row's Gaussian log-density in standardised coordinates under that level |
| `n_eff` | as everywhere: the weight behind the state before this row |

Two things to know. `u` is a **downward biased** estimate of the
equicorrelation — the paper says so, and it is a ratio of two averages, so
`E[u]` is about 0.20 for a true 0.30 at six columns. Use it as a signal that
moves with the market's correlation, not as the correlation. And `rho` is
not the same thing as an `ew_cov`'s `corr` over the columns; the mean of a
ratio is not the ratio of means, and the gap is large.

`blocks` maps a name to a subset of the features. The model then estimates
one number per block and one per pair of blocks — the useful middle between
one correlation and all of them — and the outputs become `u_<A>`, `u_<A>_<B>`
and their `rho_*` twins, with one `loglik` over all of them. Every feature
must be in exactly one block, and a block needs at least two.

```python
eq = po.spec.deco("eq", features=["x0", "x1", "x2"],
                  clock="t", max_dclock=300.0, halflife=500.0)
blocked = po.spec.deco("blocks", features=["x0", "x1", "x2", "signal_a"],
                       blocks={"fast": ["x0", "x1"], "slow": ["x2", "signal_a"]},
                       halflife=500.0)
out = df.online.fit_predict([eq, blocked])
```

### `bocpd` — how long has this regime lasted?

*API:* [`po.spec.bocpd`](https://hgilde.github.io/polars-online/spec.html#polars_online.spec.bocpd) — *Rust:* [`bocpd.rs`](crates/online-core/src/bocpd.rs) — *Outputs:* [fields](docs/OUTPUTS.md#bocpd)

Every other detector here answers "has something changed?" with a statistic.
`bocpd` (Adams & MacKay 2007) keeps a posterior over the **run length** — how
many rows since the last break — so the answer carries the age of the regime
with it. "We are forty rows into a regime" is different information from
"something broke".

Their Algorithm 1, with `H = 1/hazard` and `π_r` run `r`'s posterior
predictive for this row:

```
growth:      P(r_t = r+1, x_1:t) = P(r_t-1 = r, x_1:t-1)·π_r·(1 − H)
changepoint: P(r_t = 0,   x_1:t) = Σ_r P(r_t-1 = r, x_1:t-1)·π_r·H
```

Each run keeps its own conjugate sufficient statistics, so slot `r` holds
exactly the `r` rows that hypothesis says came before this one in the run —
and slot 0 holds none, so its predictive is the prior's. That is what makes
"a new run starts here" a hypothesis the data can vote on.

The vector would grow by a slot every row. Two knobs bound it.
`prune_below` drops the runs holding less than that share of the mass.
`max_run` folds every longer run into the last kept one, which takes their
mass and keeps its own statistics — so it caps how much history any run
holds, and `run_mode` saturates one below it.

| output | meaning |
|---|---|
| `p_change` | `P(r_t ≤ 1)` given this row: the alarm |
| `run_mode` | the most likely run length, before the row — so **`t − run_mode` is the row the run began on** |
| `run_mean` | the posterior mean run length, before the row |
| `pred_<f>` | the pre-row predictive mean of each feature, mixed over runs |
| `logscore` | the row's log predictive density under that mixture |
| `n_eff` | as everywhere: the weight behind the state before this row |

**`run_mode` is the answer; `p_change` is the alarm.** The two are not the
same quality of signal. `p_change` is a per-row likelihood ratio, so it is
spiky, and its height depends on the size of the break against the prior
scale. A ten-fold variance step takes it to 0.83 on the row itself. A
four-sigma mean shift with a diffuse prior barely lifts it. A change in
correlation alone never moves it at all. The run length finds all three, one
to three rows later, and dates them to the right row.

It is `P(r ≤ 1)` and not `P(r = 0)` because the changepoint branch and the
growth branch share the same predictive, which makes the normalised mass at
`r = 0` *exactly* `H` on every row whatever the data. Row one of a group
reports nothing at all: `P(r ≤ 1)` is 1 there however the row looks.

`prior_scale` is the prior guess at the variance and is the one parameter
you must set from your data. Too large and the model goes quiet — no row is
ever surprising under a predictive that wide, and a real break is never
found. `prior_nu` and `prior_scale` are `2a` and `2b` in the gamma
parametrisation, which is how Adams and MacKay give their own finance
example (`a = 1`, `b = 1e-4`, `hazard = 250`).

`emission="diag"` (the default) is a normal-inverse-gamma per feature;
`"gaussian"` is a normal-inverse-Wishart over all of them, which costs
`O(runs·d²)` a row and is the one that can see a break in the *correlation*
with the marginals unchanged.

`emission="robust"` weights each row's contribution by
`(π(x)/π(mode))**robust_beta` — in what the run learns *and* in the message
it passes on. A 20-σ row is then atypical under every run, every tempered
likelihood is about 1, and nothing moves. Without it that one row is a
changepoint (`p_change` 0.91) and the run it starts carries the outlier in
its mean. The knob is a trade: a whole new regime is a run of individually
forgiven rows, so at `robust_beta` above about 0.2 nothing is ever detected
again. The default of 0.1 ignores the outlier and still dates a four-sigma
shift to the right row.

```python
b = po.spec.bocpd("regime", features=["ret"], hazard=250.0,
                  prior_nu=2.0, prior_scale=[2e-4], group="bond_id")
out = df.online.fit_predict([b]).unnest("regime")
run_started_at = pl.int_range(pl.len()) - pl.col("run_mode")
```

`hazard_col` reads the hazard per row from a column, declared in the target
slot the way a weight is — a wider prior on a quiet session, a narrower one
across a data release.

## Parallelism

The unit of work is a *stream*: one spec on one group (with no `group`, one
stream per spec). On every chunk, each stream in the bank becomes one task on
the bank's own thread pool (a [rayon](https://github.com/rayon-rs/rayon)
pool, separate from polars') — one flat pool across all specs and all
groups, longest stream first so a few big groups do not leave cores idle at
the tail. Within a stream the rows go one at a time, because each row's
update depends on the last. That is what makes the numbers independent of how
the work is split. It also means a bank with one spec and one group is one
thread's work per chunk — polars' own reading and writing still run in
parallel around it.

So a bank fills the pool with groups, with specs, or with both.

A search over factor sets is a list of specs, one per set. Each spec is its
own accumulator, with its own standardization and its own grid inside, and
each is one task. A null in a factor a spec does not use costs that spec
nothing. (Subsets of one list that should share an accumulator are
`feature_sets` on one spec: one solve each, not one task.) The list runs as
one plan in one pass, with the thread counts set before anything is built:

```python
import os
os.environ["POLARS_ONLINE_MAX_THREADS"] = "8"   # the bank's pool: read at the first bank call
os.environ["POLARS_MAX_THREADS"] = "8"          # polars' readers and writers: read at import

import polars as pl
import polars_online as po
from itertools import product

factors = {"mkt": ["x0"], "mkt-sz": ["x0", "x1"], "mkt-sz-val": ["x0", "x1", "x2"]}

def spec(name, features, standardize):
    return po.spec.ewridge(f"{name}-std{standardize:d}",
                           targets=["y"], features=features, clock="t", max_dclock=300.0,
                           group="bond_id", session="session", session_gap=60.0,
                           halflife=[100.0, 1000.0], ridge=[1e-3, 0.1],   # gridded inside the spec
                           standardize=standardize)

specs = [spec(n, f, s) for (n, f), s in product(factors.items(), [False, True])]

(pl.scan_parquet("ticks.parquet")
   .online.fit_predict(specs, chunk_rows=200_000, save_state="grid.state")
   .sink_parquet("grid.parquet"))

scores = po.eval.compare_specs(pl.read_parquet("grid.parquet"),
                               [s["name"] for s in specs]).sort("r2", descending=True)
```

Every chunk puts 6 × 64 stream tasks on the pool. On 2.56M rows over 64
groups that plan takes 12.3 s at one thread and 2.2 s at fourteen; the
three-factor spec alone goes from 2.5 s to 0.62 s, because with one task
per group the fixed cost of reading and assembling each chunk shows
through. The output is one struct column per spec, which is what
`compare_specs` reads, and one state file holds them all. The same list
runs the same way through `ModelBank`, `po.run` and the CLI.

Where the parallelism comes from, then:

- **Groups.** k=20 over 64 groups: 1.02M, 1.91M, 3.52M, 6.44M and 8.20M
  rows/s at 1, 2, 4, 8 and 14 threads — **8.0×** on a 14-core machine.
- **Specs.** Eight single-group specs in one bank run in 130 ms against
  515 ms one at a time.
- **Halflives.** Each halflife in a grid is its own accumulator, and the
  instances of a stream run alongside each other (except with
  `drift_action="reset"`, which couples them). Ridge and feature-set grids
  are not parallel because they need not be: they share one accumulator and
  are expanded at solve time.
- **The runner.** `po.run` and the CLI are a three-stage pipeline — a reader
  thread, the bank on the calling thread, a writer thread — with one chunk
  in flight per stage; `ONLINE_TIMING=1` prints how long the bank waited on
  each side. Reading and writing are polars' work on polars' pool: parquet
  pages are encoded a column at a time there, NDJSON a slice per thread.
- **Python.** The GIL is released while a chunk is in the bank, so a Python
  reader thread can run ahead of `ModelBank.fit_predict`, and independent
  `po.run` calls in threads of one process share the one pool.
- **The expression form.** Under `.over("group")`, polars runs the groups
  through its own pool — which is why the plugin packs its inputs into one
  struct: the single-input path is parallel, the multi-input one is not
  (12.2M rows/s at 1000 groups).

Thread count is `POLARS_ONLINE_MAX_THREADS` for the bank's pool and
`POLARS_MAX_THREADS` for polars' readers and writers; unset, each is one
thread per core. The bank builds its pool at the first bank call, and polars
builds its own at import, so each must be set before that point — as above,
or in the shell
(`POLARS_ONLINE_MAX_THREADS=8 python fit.py`), which is the form that always
works. Set later, the variable is ignored, and
[`po.thread_pool_size()`](https://hgilde.github.io/polars-online/polars_online.html#polars_online.thread_pool_size) says what took (`pl.thread_pool_size()` for
polars'). A value that is not a count is refused by name at the first bank
call. It changes the speed and nothing else: `tests/test_portability.py`
runs the same stream at 1 and 8 threads in separate processes and requires
identical output. Everything is one process — there is no distributed
execution, by design (see [What this is not](#what-this-is-not)).

Two knobs because the two counts do different things. Polars' also sizes
what its parquet reader holds in flight — it prefetches row groups ahead of
the consumer, so more threads is a bigger pile of decoded rows
([Memory](#memory-which-calls-stream), above) — while the bank's count buys speed and
nothing else. So a run that has to fit in a smaller box keeps polars small
and gives the bank every core:

```python
import os
os.environ["POLARS_MAX_THREADS"] = "4"           # the reader's prefetch is sized from this
os.environ["POLARS_ONLINE_MAX_THREADS"] = "14"   # the bank still has every core

import polars as pl
import polars_online as po

(pl.scan_parquet("ticks.parquet")
   .online.fit_predict([spec], chunk_rows=200_000)
   .sink_parquet("fit.parquet"))
```

On 12M rows over 64 groups, one spec: 14 and 14 takes 2.6 s at a peak of 1.1
GB. 4 and 14 takes the same 2.6 s at 0.8 GB — a third less memory at the same
speed. One shared count of 4 takes 3.9 s at 0.6 GB, and polars alone at one
thread 7.4 s, because reading and writing are then one thread's work. Six specs
split the same way: 10.3 s at 1.5 GB, 11.8 s at 1.2 GB, 16.6 s at 1.0 GB.
(Memory here and below is the peak footprint `/usr/bin/time -l` reports. RSS
reads about 0.7 GB higher, because the memory-mapped input file counts there.)
The pools never wait on each other, because a bank task never calls back into
polars' pool. So oversubscribing both costs no time either: 28 and 28 on 14
cores ran the grid above in 2.18 s against 2.21.

### Chunk size

`chunk_rows` is how many rows the bank takes at a time. It is a keyword on
`lf.online.fit_predict` and `lf.online.predict`, on `po.run`, and on the
CLI (`--chunk-rows`, or `chunk_rows` in the TOML); the default is 100,000.
With `ModelBank.fit_predict(df)` the chunk is whatever frame you pass.

It never changes the numbers. One chunk or a thousand gives the same
output; the one thing that moves is where `coef` lands, because each stream
reports its coefficients on its last row of every chunk. `coef_every` gives
it a cadence that does not move.

What it changes is speed and memory, and two things pull against each
other:

- **A chunk can only run the groups it holds.** If the file is sorted by
  group, a 100k chunk holds one or two groups, so the bank has one or two
  tasks per chunk and most cores sit idle. Bigger chunks fix that.
- **Bigger chunks lose the overlap.** Reading, fitting and writing run side
  by side, a chunk apart. With huge chunks the stages spend more time
  waiting on each other, and the chunks in flight (about three) cost
  memory.

The same 12M rows and 64 groups, one spec, 14 threads:

| `chunk_rows` | groups interleaved | sorted by group | peak memory |
|---|---|---|---|
| 20,000 | 2.7 s | 9.2 s | 1.0 GB |
| 50,000 | 2.5 s | 8.8 s | 0.9 GB |
| 100,000 (default) | 2.4 s | 8.1 s | 1.0 GB |
| 200,000 | 2.6 s | 7.1 s | 1.1 GB |
| 500,000 | 2.8 s | 4.6 s | 1.5 GB |
| 1,000,000 | 3.2 s | 4.2 s | 1.8 GB |
| 2,000,000 | 4.5 s | 6.0 s | 2.4 GB |

(The memory column is the interleaved file's; the sorted file is within
0.1 GB of it, except 2.6 GB at 2M.) So:

- **Groups mixed through the file** — tick data in time order — leave the
  default. Everything from 50k to 500k is within 0.4 s of it.
- **Sorted or clustered by group** — raise it until a chunk spans several
  groups: a few times the rows per group. This file has about 190k rows
  per group, and 1M, five groups a chunk, is twice as fast as the default.
  Past that the overlap goes and memory climbs; the interleaved file shows
  the cost, with 2M nearly twice the default's time.
- **A smaller box** — lower it, but expect little below the default: most
  of the first gigabyte is polars' reader prefetch, not the chunks, and
  `POLARS_MAX_THREADS` or `POLARS_ROW_GROUP_PREFETCH_SIZE` is what shrinks
  that ([Memory](#memory-which-calls-stream), above).

## Performance

Apple M-series, single process, best of 3, 200k rows per run
(`uv run python scripts/benchmark.py --markdown`):

| configuration | notes | rows/sec |
|---|---|---|
| `ewridge` k=5 | 1 target, 1 halflife | 10,306,024 |
| `ewridge` k=20 | 1 target, 1 halflife | 4,122,355 |
| `ewridge` k=50 | 1 target, 1 halflife | 1,040,742 |
| `ewridge` k=20 | 10 targets | 2,340,621 |
| `ewridge` k=20 | 5 halflives | 2,355,548 |
| `rls` | k=20, 1 target | 1,843,129 |
| `kalman` | k=20, 1 target | 2,122,805 |
| `lasso` | k=20, 1 target (3-point path) | 2,060,329 |
| `huber` | k=20, 1 target | 4,175,881 |
| `ftrl` | k=20, 1 target | 6,006,156 |

Targets share one `S` accumulator, so 10 targets cost far less than 10× one.
Each halflife in a grid is its own accumulator, but they run in parallel, so
a 5-halflife grid costs about 2× one rather than 5×. `rls` pays 1.3–2.1× for
the square-root form that keeps it from dying of cancellation on one extreme
row; that is worth it.

The other families, and the options that add a pass, on the same machine
and rows:

| configuration | notes | rows/sec |
|---|---|---|
| `ewridge` + `conformal` | k=20, 90% interval | 4,097,559 |
| `sgd` | k=20, squared loss | 8,961,276 |
| `sgd` | k=20, `coef_min=0`, `coef_sum=1` | 2,481,671 |
| `pa` | k=20 | 11,150,059 |
| `kalman` | k=20, `revert_halflife` | 1,844,755 |
| `ew_cov` | k=20: mean, std, corr (230 statistics) | 2,117,751 |
| `ew_cov` | k=20: mean, mahal, `mahal_q0.99` | 750,903 |
| `ew_class` | k=20, 3 classes, full covariance | 495,968 |
| `ew_class` | k=20, 3 classes, shared covariance | 577,082 |
| `ew_class` | k=20, 3 classes, diagonal | 2,824,922 |
| `kmeans` | 4 features, K=8 | 6,118,711 |
| `kmeans` | k=20, K=8 | 3,103,963 |
| `micro` | 4 features, `eps=1` | 15,060,666 |
| `seqtest` | sign of one column | 22,436,196 |

A conformal interval is free: it reads the residual the model already has.
A simplex constraint sorts `2k` breakpoints per row, so it costs `sgd`
about 4×. `ew_cov` writes 230 numbers a row and still runs at half the
speed of one `ewridge`. The Mahalanobis distance and the full-covariance
`ew_class` each pay for a Cholesky factor of a `k × k` matrix, one per row
for `mahal` and one per *learned* row for `ew_class` — the classes a row
does not touch keep theirs. `kmeans` and `micro` cost a distance to each
centre; `seqtest` a handful of operations.

The correlation families, on the same machine and rows:

| configuration | notes | rows/sec |
|---|---|---|
| `deco` | k=20, one equicorrelation | 2,666,809 |
| `deco` | k=20 in 4 blocks | 1,892,836 |
| `rcov` | 4 features, kernel, blocks of 1000 | 3,864,849 |
| `ew_cov` | k=20: mean, cov, lags 1–5 | 1,148,104 |
| `hmm` | 4 features, K=2 | 1,343,086 |
| `hmm` | k=20, K=2 | 353,258 |
| `bocpd` | 4 features, diagonal | 1,009,075 |
| `bocpd` | 4 features, full covariance | 585,180 |
| `corrchange` | 4 features, monitor, `span_rows=500` | 388,557 |
| `corrchange` | 4 features, window 100, permute every 500 | 187,407 |

`deco` is one number for the whole matrix and costs `O(m)` a row, which is
why it runs at `ew_cov`'s speed and not at a covariance matrix's. `rcov`
accumulates per row and pays for its kernel only when the block closes.
`hmm` factorizes a `k × k` covariance per state per row, which is
`ew_class`'s cost with the classes hidden. `bocpd` costs
`O(runs · d²)`, and the length of the run vector is the whole story —
see below.

**`prune_below` is not a tuning knob on `bocpd`, it is what makes it finite.**
The run vector grows by one entry every row, so with `prune_below = 0` the
model is `O(rows²)`: measured at 1,897 / 947 / 472 rows/s on 5k / 10k / 20k
rows, halving each time the stream doubles. At the default `1e-6` it is
flat in the length of the stream, and the knob is a direct dial on
throughput — 204k, 324k and 687k rows/s at `1e-8`, `1e-6` and `1e-4` on
i.i.d. Gaussian rows. `max_run` is the belt to that pair of braces and
usually never binds. One consequence worth knowing: **`bocpd` is faster on
data that actually breaks**, because a changepoint collapses the posterior
onto a short run — the 1.0M rows/s in the table is on data with regimes,
against 324k on a stationary stream.

`corrchange`'s window kind is the slowest model here, and deliberately: the
permutation null re-draws `n_perm` statistics every `permute_every` rows.
At the default cadence that is a `O(n_perm · window · k²)` job amortized
over 500 rows; `crit` given as a number skips it entirely.

Grouped data goes wider, as [Parallelism](#parallelism) shows: 8.2M rows/s
at k=20 over 64 groups.

**Memory** is three things: the state, the chunks in flight, and whatever
polars' reader prefetches. Three chunks are in flight at once, so
`chunk_rows` is the knob for the middle one ([Chunk size](#chunk-size),
above). The prefetch is usually the largest of the three: on a 14-thread
machine the parquet reader front-loads ~0.7 GB of decoded row groups
whatever the file's length, and `POLARS_ROW_GROUP_PREFETCH_SIZE=1` takes the
CLI to 0.15 GB at the same speed. It is sized from the thread count, so
`POLARS_MAX_THREADS` shrinks it too. Where the time goes, and what to reach
for, is in [docs/PERFORMANCE.md](docs/PERFORMANCE.md).

## Against scikit-learn

The model most people compare this to is
[`SGDRegressor.partial_fit`](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.SGDRegressor.html),
and the honest comparison starts by naming what each side is. `SGDRegressor`
is a first-order stochastic optimiser: its answer depends on the learning
rate, the schedule, the feature scaling and the row order. The primary
regression here is a different algorithm class — `ewridge`, `rls`, `lasso`,
`huber` and `quantile` accumulate sufficient statistics and solve, so there
is no learning rate, and with decay off `ewridge` is ordinary least squares
to 2e-13 of `numpy.linalg.lstsq` in any row order. The counterpart to
`SGDRegressor` here is [`sgd`](#sgd--stochastic-gradient-descent), the cheap
`O(k)` baseline.

Measured on one generated stream, 100,000 rows, `k = 20`, each contender at
its best over a sweep of its own settings (`scripts/sklearn_comparison.py`,
scikit-learn 1.9.0, `docs/PERFORMANCE.md` §19 for the full tables and the
sweeps). The noise ceiling is the R² of the generating signal itself:

| contender | R² stationary | R² drifting | rows/sec | what a prediction saw |
|---|---:|---:|---:|---|
| noise ceiling | 0.9831 | 0.9923 | | |
| `SGDRegressor`, row by row | 0.9829 | 0.9899 | 3,400 | every row before it |
| `SGDRegressor`, batches of 1,000 | 0.9830 | 0.9820 | 2,200,000 | every row before its *batch* |
| `po.spec.sgd` | 0.9826 | 0.9906 | 6,000,000 | every row before it |
| `po.spec.ewridge` | **0.9831** | **0.9907** | 575,000 | every row before it |

Three things to take from it. **Accuracy is not the difference.** Everything
reaches the ceiling on the stationary stream, and the row-by-row contenders
are within 0.001 of each other on the drifting one — `SGDRegressor` and
`po.spec.sgd` at the same constant step give the same number to four
places, because they are the same recursion. The one gap in the table, the
batched 0.9820, is staleness: a prediction made up to 999 rows before its
update. **The speed difference is a difference in semantics**: sklearn's fast
form updates once per 1,000-row batch, while every prediction here is made
from the state as it stands. Asked for that same guarantee — `partial_fit`
per row — `SGDRegressor` runs at 3,400 rows/second, and the cost is Python's
per-row overhead rather than the algorithm. And **a halflife is not the
advantage here**: a constant learning rate forgets too, at about `1/eta`
rows, and on evenly spaced rows that is a halflife. The clock matters when
the rows are not evenly spaced, which this stream is not.

**Where the Gram wins: the first hundred rows of every group.** An exact
solve is right as soon as its `X'X` is full rank, about `k` rows in; a
first-order method needs about `1/eta` rows per direction. 500 groups of 200
rows, each group its own coefficients, `k = 20`, R² by position in the group
— sklearn as one estimator and one scaler per group in a dict, row by row,
the bank as `group="g"`:

| contender | rows 25–50 | rows 50–100 | rows 100–200 | rows/sec |
|---|---:|---:|---:|---:|
| noise ceiling | 0.9896 | 0.9903 | 0.9900 | |
| `SGDRegressor` per group, at its best | 0.7277 | 0.9242 | 0.9789 | 3,342 |
| `po.spec.sgd`, `scale_features=True`, the same step | 0.7182 | 0.9213 | 0.9788 | 18,519,660 |
| `po.spec.ewridge` | **0.9693** | **0.9860** | **0.9882** | 5,130,803 |

The `sgd` row is sklearn's step, and through 0.3.1 it was not this close:
the table found `po.spec.sgd` diverging at the start of every group (R² of
−6.9 at rows 25–50), because `scale_features` standardised a row against
moments from *before* it, and a two-row variance estimate can be tiny by
chance. It now standardises against moments that include the row —
sklearn's scaler order, which is not a leak: the features are known at
prediction time, and the rule is about the target (`docs/PLAN.md` task 74).
The 0.01 of R² left between the two rows is where the *prediction* is
standardised: sklearn's loop standardises the row it predicts against the
moments before it and the row it learns from against the moments including
it, while here one standardised row serves both, so a prediction is on the
same footing as every row the coefficients were learned from; the gap
closes as the fit converges. The condition was few rows per feature, and a
wide fit is that on every row: at `k = 10,000` the scaled fit's predictions
correlated 0.52 with sklearn's before the change and 0.999997 after.

What else is different: a grid of six penalties is 3.8× the work for
sklearn (six estimators) and 1.6× here (one accumulator, six solves);
multiple targets share one `X'X`; `group=` is one state per key rather than a
dict of estimators; standardisation is streaming and cannot leak the way a
`StandardScaler` fitted on the whole frame does; chunk invariance is a test;
and the state is a versioned cross-OS file rather than a pickle.

**Where sklearn wins: a wide row against `ewridge`, by batching.**
`ewridge` keeps a `(k+1)²` matrix, so at `k = 10,000` it carries 860 MB of
state and runs at 53 rows/second, against 0.31 MB and 20,007 for
`SGDRegressor` in batches of 1,000. That cost is the matrix itself — a
rank-1 update moves all 800 MB of it every row, at 85 GB/s, near this
machine's memory bandwidth; emitting coefficients is not it (`coef_every=1`
costs 2.5%), and neither is the solve (0.2%). `gram_block_rows=1024`
touches the matrix once per 1,024 rows instead and buys 7.2× of the
throughput (379 rows/second) and none of the memory. `sgd` is the `O(k)`
answer here — 0.89 MB and 33,271 rows/second at the same width, 45,000 fed
in chunks of 20,000 rows — and it is faster than sklearn's batch with
every prediction made from the state as it stands; `SGDRegressor` asked
for the same, row by row, runs at 2,206 rows/second at that width, and the
two agree (correlation 0.9954 between their predictions unscaled, and
0.999997 with `scale_features=True`, which is sklearn's own recipe of a
scaler in front of the step). Both are run at `learning_rate = 0.2 / k`: an LMS step is
stable only while `eta · |z|² < 2`, and a standardised row has `|z|² ≈ k`.
And the ecosystem is sklearn's: pipelines, `GridSearchCV`, calibration,
and far more use. What this has is the stream.

## What this is not

A model layer, not a stream-processing framework. It expects a frame that is
already aligned — and, when a spec names a `clock`, each group's rows in
clock order — and it keeps O(state) per stream. It deliberately does **not**
provide:

- **connectors or ingestion** — feed it whatever Polars can read;
- **event-time windowing, asof or interval joins** — build features with
  Polars expressions upstream, or with a streaming framework such as
  [Pathway](https://pathway.com);
- **watermarks or late-arrival policy** — `clock`, `max_dclock`,
  `on_clock_reset` and `session` describe time *within* a stream, not
  pipeline lateness. Under a `clock`, a row that arrives out of order is a
  data error, and `on_clock_reset="error"` will say so;
- **distributed execution** — one process, a thread pool across (spec × group).

Those boundaries make the two compose:
[examples/pathway_integration.py](examples/pathway_integration.py) runs a
`ModelBank` as a stateful operator inside a Pathway pipeline — Pathway does
ingestion, event-time alignment and windowing; we do the model. Chunk
invariance means the engine's batching cannot change the numbers, and
`save_bytes`/`load_bytes` let a pipeline checkpoint carry the model state.
Pathway is not a dependency; the example imports it lazily.

## Versioning and the Polars pin

### What is pinned

| py-polars | rust polars | pyo3-polars | pyo3 | Python |
|---|---|---|---|---|
| **>= 1.34.0, < 2** (built and tested against 1.44.1) | 0.55.2 | 0.28 | 0.29 | ≥ 3.12 (`abi3-py312`) |

The Rust `polars` is pinned exactly and linked into the wheel; the runtime
requirement is a range, because the two copies never meet. The floor is
`LazyFrame.collect_batches`, which `po.run` and `lf.online.fit_predict` read
with and py-polars added in 1.34.0; the whole suite passes on 1.34.0, 1.38.1
and 1.44.1 with identical numbers. `ModelBank` and the expression form alone
work from 1.28.1 (tested across 17 releases). `tests/test_scaffold.py`
asserts the pins; the matrix is in
[docs/RELEASE-READINESS.md](docs/RELEASE-READINESS.md).

### Why a mismatch is an error, not a crash

`ModelBank` and the expression plugin move data across the boundary through
the Arrow C Data Interface, the same cross-language ABI pyarrow and DuckDB
use; nothing here uses the version-sensitive types that cross as serialized
query plans. The only thing `ModelBank` asks of the Python side is
`PySeries._export` / `_import`, and a Polars without them fails with a clean
`AttributeError` before any data moves. The plugin loader goes further and
negotiates its ABI, refusing a major it does not know. The pin exists so you
never see those messages, not because something worse waits behind them.

### Which interfaces carry a promise

Polars supports three, and only one carries a guarantee:

- the **expression plugin** — the supported path, with a negotiated handshake;
- **pyo3-polars' extension types** (`ModelBank`) — provided "for
  convenience", with no guarantee beyond the latest definitions working for
  the latest Polars;
- the **IO plugin** (`lf.online.fit_predict`) — documented, but `@unstable`
  in py-polars.

The two that stream are the two without a promise, so a break on a new
Polars is expected maintenance, not a surprise.

### How the pin moves

A weekly job ([`polars-canary.yml`](.github/workflows/polars-canary.yml))
drops the range from `pyproject.toml`, installs the newest py-polars — a 2.0
included, the week it appears — builds the wheel as CI does and runs the
whole suite. Only polars moves in that run, so a red canary means Polars
broke us and nothing else. The response is decided in advance: **cap** the
range at the last release that passed, in a patch release, so no resolver
hands anyone the broken pair; then **fix**, and widen again. Where to look
first: `ModelBank`, then the IO-plugin tests in `tests/test_frame.py`, then
the plugin. The Rust copy of polars moves by hand, together with
pyo3-polars, polars-arrow, polars-parquet and polars-utils, through CI.

### This package's own versioning

Semantic versioning. While pre-1.0 the **minor** version carries breaking
changes, and any change to the numbers a model returns, so pin `~=0.4.0` if
you need stability. Widening the Polars range is a minor release; narrowing
it is breaking. See [CHANGELOG.md](CHANGELOG.md). Output field names are
part of the API ([above](#output-field-names)).

## Testing

The guarantees above are only worth what checks them, so the suite is built
around oracles and invariants rather than expected values typed in by hand.
About 650 Rust tests and 2,200 pytest cases (from some 1,250 test functions),
all green on three OSes; [docs/TESTING.md](docs/TESTING.md) is the ledger of what
each part proves and what it has found.

**Against references.** `ewridge` and `rls` match numpy references in
`tests/reference.py` to 1e-9, `kalman` to ~1e-15, `huber` and `quantile` to
~1e-13, `ftrl` to ~1e-16; `rls` equals `ewridge(ridge_decay=True)` solved
every row to <1e-9. The lasso is checked against the KKT conditions of its
objective rather than a ported solver, which cannot share a bug with it.
[river](https://riverml.xyz) is an independent implementation of several of
the same algorithms. Its FTRL recursion agrees with ours to 1e-12 row for
row; its EW moments agree in closed form and in the limit; its quantile and
Huber models agree statistically. Two convention differences are pinned as
tests
rather than left as surprises.

**Invariants, for every model.** Each of these is checked at the bank, and
where it applies at the expression and CLI levels too:

| | |
|---|---|
| chunk invariance | one chunk, seven, four hundred, one row at a time, and with a save and load in the middle |
| thread invariance | 1 thread against 8 |
| group independence | a group's numbers do not depend on what else is in the bank |
| the paths agree | expression ≡ bank; runner ≡ bank for every input source and format |
| `predict` ≡ `fit_predict` | of the next row, field for field, with every diagnostic on |
| stream semantics | the null policy, warm-up, and the clock |
| `n_eff` | the same recursion in every model (`crates/online-core/tests/model_contract.rs`) |

Hypothesis generates adversarial streams — mixed nulls, duplicate and
long-gap clocks, values at ±1e8, zero weights, tiny groups — and asserts the
strongest one: **changing a row's own target never changes that row's own
prediction.** IC ≈ 0 on pure-noise targets says the same thing from the
other side.

**Fixed numbers.** One golden stream per model in the Rust core, and the
whole pipeline — extraction, fan-out, diagnostics, struct assembly — pinned
to fixed output and compared on every OS, so a divergence in polars'
vectorized paths on another CPU would show.

**Hardening.** What the suite does to a bank on purpose:

| | |
|---|---|
| everything at once | a 30k-row stream with every output switched on, compared by digest across chunkings, a mid-stream save and load, and thread counts |
| weight scale | all weights ×1e±6 changes nothing but `n_eff` |
| parameter edges | `halflife` from `1e-3` to `inf` |
| a corrupt state file | any byte flipped fails cleanly, and never panics |
| concurrent misuse | two threads calling `fit_predict` at once get a clean error |
| copying a bank | `pickle` and `copy.deepcopy` resume bit-exactly |
| across the FFI | memory safety where two copies of Polars share one process |
| sustained load | a 10M-row soak, opt-in with `pytest -m soak` |

**Contracts that are files.** The public API — every name, default and
signature, every output field name — is a checked-in snapshot
(`tests/api_surface.txt`), so a change is a reviewable diff. Every python
block in this README runs. Everything under `examples/` runs unmodified —
the TOML through the real CLI, the Pathway operator end to end.
`docs/VALIDATION.md`, where the defaults were chosen,
is regenerated and compared, so the numbers behind them cannot silently stop
being true. A data file, a large file or generated output that gets tracked
fails a test. Bank files from the previous schema version still load.

**Beyond the suite.** `cargo mutants` runs over the core: the last pass left
8.3% of 2,616 mutants surviving, clustered where only the Python suite
reaches (`cargo test` cannot see it), which is what the golden and contract
tests in Rust were added for. Coverage is 96% of the Python package and 75%
of Rust regions, understated for the same reason.

**Where it runs.** `./scripts/gate.sh` before every commit — `cargo fmt`,
`clippy -D warnings`, `cargo test`, `ruff`, `mypy`, the build, `pytest`,
`sphinx -W`. CI runs the same on ubuntu, windows and macos for every push
and pull request. The release build writes a state file on macOS and
continues the stream from it on Windows and Linux. The weekly canary runs
the suite against the newest py-polars.

Tests generate or download their own data; there are no data files in the
repo. Downloads are cached under `.cache/` and skipped when offline.

## Development

```sh
uv sync                                                # Python env (CPython 3.12)
./scripts/gate.sh                                      # everything CI checks
uv run cargo test --workspace                          # Rust tests
uv run maturin develop --release -m crates/online-py/Cargo.toml
uv run pytest                                          # Python tests
uv run --group docs sphinx-build -W docs/reference docs/_build/html   # API reference
uv run python scripts/validate.py > docs/VALIDATION.md # re-run the [validate] experiments
uv run python scripts/regime_experiments.py all        # the docs/REGIMES.md experiments
uv run python scripts/benchmark.py                     # throughput
```

Prerequisites: [uv](https://docs.astral.sh/uv/) and a stable Rust toolchain
([rustup](https://rustup.rs)). `source scripts/env.sh` (`. .\scripts\env.ps1`
in PowerShell) puts both on the `PATH` for a shell; `.vscode/settings.json`
does it for VS Code's terminal. `cargo` runs via `uv run` because `online-py`
builds against pyo3's `abi3-py312` and needs a 3.12+ interpreter at build
time.

- Every document, and which to read for what: [docs/README.md](docs/README.md)
- A map for coding agents, at the repo root and on the docs site:
  [llms.txt](llms.txt) ([llmstxt.org](https://llmstxt.org))
- API reference: <https://hgilde.github.io/polars-online/> — built from the
  docstrings and published from every green push to `main`
- Design and task list: [docs/PLAN.md](docs/PLAN.md)
- Saving, serving, resuming: [docs/STATE-WORKFLOW.md](docs/STATE-WORKFLOW.md)
- Measured defaults: [docs/VALIDATION.md](docs/VALIDATION.md)
- What the regime detectors find: [docs/REGIMES.md](docs/REGIMES.md)
- Where the time and memory go: [docs/PERFORMANCE.md](docs/PERFORMANCE.md)
- Adding a model: [docs/EXTENDING.md](docs/EXTENDING.md)

## License

Apache-2.0. See [CONTRIBUTING.md](CONTRIBUTING.md) to make changes,
[SECURITY.md](SECURITY.md) to report a vulnerability.

