Skip to content

Commit 62bb2f8

Browse files
mstorathclaude
andcommitted
feat(detection): faithful L1Potts solver, hardening, tests, docs
Rewrite the L1Potts implementation to match Algorithm 1 of Storath, Weinmann & Unser (2017) explicitly. The forward DP is identical to the paper; backtracking uses explicit parent pointers (uint8 matrix + prev-argmin vector) instead of the paper's in-place float64 subtraction trick. Mathematically equivalent, ~8x less memory. Hardening: reject non-numeric dtypes, NaN/Inf in signal or weights, empty signal, and predict-before-fit; clear ValueError/RuntimeError messages instead of cryptic IndexError/AttributeError. Defensive copies of signal and weights shield the solver from caller mutation. Tests (tests/test_l1potts.py, 71 cases): - functional-value parity with Pelt+CostL1 over 5 random seeds - replication-equivalence for integer weights, including weight ratios up to 10000:1 - positive homogeneity over 7 orders of magnitude in alpha - reversal symmetry of the functional - multi-pen independence on a single fit - edge cases: n=1, K=1, n=2, shape (N,1), pen=+inf, pen=NaN, empty signal, complex/object/bool dtypes, list/tuple/read-only inputs - SyntaxWarning regression net via subprocess + compile() Docs: user-guide page (docs/user-guide/detection/l1potts.md) with usage example and reference, code-reference autodoc stub, mkdocs nav entries. Export L1Potts at top level (rpt.L1Potts). Performance: on N=5000 noisy 1D signal, ~20x faster than Pelt(model="l1", min_size=1, jump=1) with identical functional value. Implementation prepared with assistance from Claude (Anthropic) acting as a coding agent. The algorithm is the cited authors'; the agent's contributions were the parent-pointer backtracking variant, input hardening, and the test suite. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
1 parent 22869a6 commit 62bb2f8

7 files changed

Lines changed: 801 additions & 126 deletions

File tree

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
# L1Potts
2+
3+
::: ruptures.detection.l1potts.L1Potts
4+
rendering:
5+
show_root_heading: true
Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,61 @@
1+
# Penalized L1 Potts segmentation (`L1Potts`)
2+
3+
## Description
4+
5+
The method is implemented in [`L1Potts`][ruptures.detection.l1potts.L1Potts].
6+
It computes the global minimizer of the **L1 Potts functional** for piecewise constant 1D signals:
7+
8+
$$
9+
\min_{u \in \mathbb{R}^N} \;\; \gamma \sum_{i=1}^{N-1} \mathbb{1}(u_i \neq u_{i+1}) \;+\; \sum_{i=1}^{N} w_i \, |f_i - u_i|
10+
$$
11+
12+
where $f$ is the observed signal, $w$ are non-negative per-sample weights, and $\gamma > 0$ is the jump penalty.
13+
14+
The L1 fit makes the estimator robust to heavy-tailed noise and outliers, in contrast to the L2 Potts model used by other ruptures detectors (`Pelt(model="l2")`, `Dynp(model="l2")`).
15+
16+
The implementation is **Algorithm 1 of [[Storath2017]](#Storath2017)**, which solves the problem exactly in $\mathcal{O}(KN)$ time, where $N$ is the number of samples and $K$ the number of distinct values in the signal. The algorithm is much faster than `Pelt(model="l1")` (a 20–30× speedup is typical on a few-thousand-sample noisy signal). It uses a Viterbi-type dynamic program over (level, sample) pairs, where the candidate levels are the unique observed values — the optimal segment level is always one of them, since the weighted L1 median lies in the data.
17+
18+
`L1Potts` accepts only 1D signals. Penalty-only mode (`predict(pen=...)`) is the only supported prediction mode; `n_bkps` and `epsilon` are not.
19+
20+
21+
## Usage
22+
23+
```python
24+
import numpy as np
25+
import matplotlib.pylab as plt
26+
import ruptures as rpt
27+
28+
# creation of data with heavy-tailed (Laplace) noise
29+
n, sigma = 500, 1.0
30+
n_bkps = 3
31+
signal, bkps = rpt.pw_constant(n, 1, n_bkps, noise_std=sigma)
32+
signal = signal.ravel() + np.random.default_rng(0).laplace(scale=sigma, size=n)
33+
34+
# change point detection
35+
algo = rpt.L1Potts().fit(signal)
36+
my_bkps = algo.predict(pen=3.0)
37+
38+
# show results
39+
rpt.show.display(signal, bkps, my_bkps, figsize=(10, 6))
40+
plt.show()
41+
```
42+
43+
To downweight known outlier samples, pass per-sample weights to `fit`:
44+
45+
```python
46+
weights = np.ones(n)
47+
weights[outlier_indices] = 1e-3 # near-zero weight effectively ignores those samples
48+
my_bkps = rpt.L1Potts().fit(signal, weights=weights).predict(pen=3.0)
49+
```
50+
51+
`fit_predict` is also available:
52+
53+
```python
54+
my_bkps = rpt.L1Potts().fit_predict(signal, pen=3.0, weights=weights)
55+
```
56+
57+
58+
## References
59+
60+
<a id="Storath2017">[Storath2017]</a>
61+
Storath, M., Weinmann, A., & Unser, M. (2017). Jump-penalized least absolute values estimation of scalar or circle-valued signals. *Information and Inference: A Journal of the IMA*. Preprint: <https://bigwww.epfl.ch/preprints/storath1602p.pdf>

‎mkdocs.yml‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -49,6 +49,7 @@ nav:
4949
- 'Binary segmentation': user-guide/detection/binseg.md
5050
- 'Bottom-up segmentation': user-guide/detection/bottomup.md
5151
- 'Window sliding segmentation': user-guide/detection/window.md
52+
- 'L1 Potts': user-guide/detection/l1potts.md
5253
- Cost functions:
5354
- 'CostL1': user-guide/costs/costl1.md
5455
- 'CostL2': user-guide/costs/costl2.md
@@ -92,6 +93,7 @@ nav:
9293
- Binseg: code-reference/detection/binseg-reference.md
9394
- BottomUp: code-reference/detection/bottomup-reference.md
9495
- Window: code-reference/detection/window-reference.md
96+
- L1Potts: code-reference/detection/l1potts-reference.md
9597
- Cost functions:
9698
- 'CostL1': code-reference/costs/costl1-reference.md
9799
- 'CostL2': code-reference/costs/costl2-reference.md

‎src/ruptures/__init__.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
"""Offline change point detection for Python."""
22

33
from .datasets import pw_constant, pw_linear, pw_normal, pw_wavy
4-
from .detection import Binseg, BottomUp, Dynp, KernelCPD, Pelt, Window
4+
from .detection import Binseg, BottomUp, Dynp, KernelCPD, L1Potts, Pelt, Window
55
from .exceptions import NotEnoughPoints
66
from .show import display
77

‎src/ruptures/detection/__init__.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,5 +4,6 @@
44
from .bottomup import BottomUp
55
from .dynp import Dynp
66
from .kernelcpd import KernelCPD
7+
from .l1potts import L1Potts
78
from .pelt import Pelt
89
from .window import Window

0 commit comments

Comments
 (0)