whitsmooth-rust 0.1.3
Batched robust Whittaker smoother (divided differences) implemented in Rust (PyO3).
License: MIT
Requires Python: >=3.9

# whitsmooth_rust

[![PyPI](https://img.shields.io/pypi/v/whitsmooth_rust)](https://pypi.org/project/whitsmooth_rust/)
[![Python](https://img.shields.io/pypi/pyversions/whitsmooth_rust)](https://pypi.org/project/whitsmooth_rust/)
[![License](https://img.shields.io/pypi/l/whitsmooth_rust)](https://pypi.org/project/whitsmooth_rust/)
[![Wheels](https://github.com/MarcYin/whitsmooth_rust_repo/actions/workflows/wheels.yml/badge.svg)](https://github.com/MarcYin/whitsmooth_rust_repo/actions/workflows/wheels.yml)

A **fast, compiled** (no JIT warmup) implementation of a **batched robust Whittaker smoother** with **divided differences**
in **Rust**, exposed to Python via **PyO3**.

It is designed for data shaped as many independent time series:

- `S` series (samples / pixels / bands / locations)
- `T` time points per series

and supports **time- and series-varying weights** (e.g. missing data, quality flags, IRLS robust weights).

This package provides both **float64** and **float32** implementations.

---

## Installation

### Install from PyPI (recommended)

```bash
python -m pip install --upgrade pip
python -m pip install whitsmooth_rust
```

Wheels are built in CI for:

- **Python**: CPython 3.9+
- **Platforms**: Linux (manylinux x86_64), macOS (x86_64 + arm64), Windows (x86_64)

### Build locally (from source)

If you are on an unsupported platform (or want to build the latest code), install Rust and `maturin`:

```bash
python -m pip install maturin
git clone https://github.com/MarcYin/whitsmooth_rust_repo.git
cd whitsmooth_rust_repo
maturin develop --release
```

This builds and installs the extension into the current Python environment.

---

## Mathematical background

We solve a Whittaker-type penalized least squares problem for each series `s`:

$$
\min_{z_s}\; \sum_{t=1}^T w_{s,t}(y_{s,t}-z_{s,t})^2 + \lambda \lVert D z_s \rVert^2
$$

- `y_{s,t}`: observation
- `w_{s,t} >= 0`: weight (0 means missing/ignored)
- `D`: **divided-difference** operator of order `d` defined on a strictly increasing sampling grid `x`
- `λ` (lambda): smoothing strength

The normal equations are:

$$
(\mathrm{diag}(w_s) + \lambda D^\top D + \text{ridge} I) z_s = w_s \odot y_s
$$

`P = D^T D` is **symmetric positive semidefinite** and **banded** with half-bandwidth `k = 2d`.
We store only the **lower band** `Pb[j, i] = P[i, i-j]`.

We solve each system using an **SPD banded Cholesky** factorization:

$$
A = L L^\top
$$

which exploits symmetry and banded structure.

---

## Robust IRLS (iteratively reweighted least squares)

Robust smoothing down-weights outliers by iterating:

1. Solve Whittaker system with current weights
2. Compute residuals `r = y - z`
3. Estimate per-series scale (default: MAD)
4. Compute standardized residuals `u = r / scale`
5. Update robust weights `w_rob(u)` and repeat

Final weights per entry are:

$$
w_{total} = w_{base} \cdot w_{rob}
$$

Supported robust weight functions:

- **Tukey / Bisquare** (hard rejection outside `c`)
- **Huber** (soft clipping)
- **Cauchy**
- **Welsch**
- **Fair**
- **Hampel** (3-parameter redescending)

Supported scale estimators:

- `mad` (default): `scale = 1.4826 * median(|r|)`
- `huber`: iterative Huber M-scale (more stable for heavy tails)

---

## x normalization (recommended)

Because divided differences depend on the units/spacing of `x`, a fixed `λ` can mean different smoothness
when sampling density changes. You can normalize `x` to approximately unit spacing:

- `normalize="mean"`: average step becomes ~1
- `normalize="median"`: median step becomes ~1 (robust)

This makes `λ` more comparable across datasets.

---

## Handling near-duplicate x (recommended for float32)

If your `x` has **very small gaps** (near-duplicate sampling positions), the divided-difference penalty can become
ill-conditioned (especially in float32). By default, the robust solvers use `merge_x_tol=1e-2` (set `merge_x_tol=None`
to disable). Options:

- Use `robust_whittaker_irls_f64` (most stable).
- Tune `merge_x_tol` to merge adjacent points within a tolerance (in input-x units).
- Increase `eps` and/or `ridge` for float32.

Long runs of **missing values** (especially missing tails) reduce the data constraints at the edges and can also make
the float32 solve ill-conditioned. Mitigations: use `robust_whittaker_irls_f64`, increase `ridge`/`eps`, use
`weighting="huber"`, or reduce `iterations`.

---

## API overview

All functions use **series-first layout** for speed:

- `y_st`: shape `(S, T)`
- `w_base_st`: shape `(S, T)` or omitted

### Build penalty band

- `build_pb_f64(x, d=2, normalize=None, eps=None) -> (pb, k)`
- `build_pb_f32(x, d=2, normalize=None, eps=None) -> (pb, k)`

Returns:
- `pb`: array of shape `(k+1, T)` storing lower band of `P = D^T D`
- `k = 2d`

### Robust IRLS smoother

- `robust_whittaker_irls_f64(...) -> z_st` (or `(z_st, w_total_st)` if `return_weights=True`)
- `robust_whittaker_irls_f32(...) -> z_st` (or `(z_st, w_total_st)`)

Parameters (both dtypes):
- `x`: strictly increasing sampling positions (length `T`)
- `y_st`: `(S,T)` observations (NaNs allowed; they will be treated as missing)
- `w_base_st`: `(S,T)` base weights (optional; if provided, NaNs in y are forced to 0 weight)
- `lam`: smoothing parameter λ (default: `10`)
- `d`: difference order (default: `2`)
- `iterations`: IRLS iterations (default: `1`; typical 3..10)
- `weighting`: one of `"tukey"|"huber"|"cauchy"|"welsch"|"fair"|"hampel"` (default: `"tukey"`)
- `scale`: `"mad"` or `"huber"` (default: `"mad"`)
- `ridge`: small diagonal stabilizer (default: `1e-10`; f32 may need larger values on difficult grids)
- `normalize`: None / `"mean"` / `"median"` (default: `"mean"`)
- `eps`: numerical floor for stability (default: `1e-10`; f32 may need larger values on difficult grids)
- `merge_x_tol`: tolerance (in input-x units) to merge near-duplicate sampling points by base-weighted averaging before solving (default: `1e-2`, set to `None` to disable)
- `parallel`: parallelize across series (Rayon). Recommended when S is large. (default: `True`)
- `tuning`: optional tuple
  - for Tukey/Huber/Cauchy/Welsch/Fair: `(0,0,c)`
  - for Hampel: `(a,b,c)`
- `return_weights`: return final `w_total_st` as well (default: `False`)

---

## Usage examples

### Float64

```python
import numpy as np
import whitsmooth_rust

rng = np.random.default_rng(0)
T, S = 2000, 1000
x = np.cumsum(rng.uniform(0.8, 1.2, size=T)).astype(np.float64)

y = rng.normal(size=(S, T)).astype(np.float64)
y[2, T//2] += 25.0
y[4, 10:25] = np.nan

w = rng.uniform(0.2, 1.0, size=(S, T)).astype(np.float64)
w[~np.isfinite(y)] = 0.0

z, wtot = whitsmooth_rust.robust_whittaker_irls_f64(
    x, y, w,
    lam=50.0, d=2, iterations=3,
    weighting="tukey",
    scale="mad",
    ridge=1e-10,
    normalize="mean",
    eps=1e-12,
    parallel=True,
    tuning=None,
    return_weights=True
)
```

### Float32

```python
x32 = x.astype(np.float32)
y32 = y.astype(np.float32)
w32 = w.astype(np.float32)

z32 = whitsmooth_rust.robust_whittaker_irls_f32(
    x32, y32, w32,
    lam=50.0, d=2, iterations=6,
    weighting="tukey",
    scale="mad",
    ridge=1e-6,
    normalize="mean",
    eps=1e-6,
    parallel=True,
    tuning=None,
    return_weights=False
)
```

### If your data is (T,S)
Transpose when calling:

```python
z_st = whitsmooth_rust.robust_whittaker_irls_f64(x, y_ts.T, w_ts.T, ...)
z_ts = z_st.T
```

---

## Notes on performance

- The solver is **SPD banded Cholesky**, per series, with half-bandwidth `k=2d`.
- Parallelization is across series `S`.
- Float32 can be significantly faster and reduce memory footprint for very large `S`.
- For highly irregular / near-duplicate `x`, prefer `merge_x_tol` or `robust_whittaker_irls_f64` (float32 may fallback to float64 internally).

---

## License
MIT

---

## Single-pass (non-robust) Whittaker solver

If you don't need IRLS robustness, you can run a single solve (much cheaper):

- `whittaker_solve_f64(x, y_st, w_base_st=None, lam=10.0, d=2, ridge=1e-10, normalize="mean", eps=1e-10, parallel=True) -> z_st`
- `whittaker_solve_f32(x, y_st, w_base_st=None, lam=10.0, d=2, ridge=1e-10, normalize="mean", eps=1e-10, parallel=True) -> z_st`

These solve:
\[
(\mathrm{diag}(w_s) + \lambda D^\top D + \text{ridge} I) z_s = w_s \odot y_s
\]
with no robust reweighting.

---

## Benchmark script

This repo includes `benchmark.py` which times:

- single solve f64 / f32
- robust IRLS (Tukey) f64 / f32

Across multiple `S` values and plots runtime.

After installing the extension:

```bash
python benchmark.py
```

Optional: if you have SciPy installed, the benchmark can also compare against a sparse reference (looping over series).