Metadata-Version: 2.4
Name: rank_preserving_calibration
Version: 0.11.0
Summary: Rank-preserving calibration of multiclass probabilities: exact QP, Dykstra's projections, KL and soft-penalty variants.
Keywords: calibration,machine learning,multiclass,isotonic,dykstra,quadratic programming
Author: Gaurav Sood
Author-email: Gaurav Sood <contact@gsood.com>
License-Expression: MIT
License-File: LICENSE
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Programming Language :: Python :: 3.14
Classifier: Topic :: Scientific/Engineering :: Mathematics
Requires-Dist: numpy>=1.20
Requires-Dist: scipy>=1.12
Requires-Dist: clarabel>=0.9
Requires-Python: >=3.12
Project-URL: Homepage, https://github.com/finite-sample/rank-preserving-calibration
Project-URL: Source, https://github.com/finite-sample/rank-preserving-calibration
Project-URL: Issues, https://github.com/finite-sample/rank-preserving-calibration/issues
Project-URL: Documentation, https://finite-sample.github.io/rank-preserving-calibration/
Description-Content-Type: text/markdown

## Rank Preserving Calibration of Multiclass Probabilities

[![Python application](https://github.com/finite-sample/rank-preserving-calibration/actions/workflows/ci.yml/badge.svg)](https://github.com/finite-sample/rank-preserving-calibration/actions/workflows/ci.yml)
[![PyPI version](https://img.shields.io/pypi/v/rank-preserving-calibration.svg)](https://pypi.org/project/rank-preserving-calibration/)
[![Documentation](https://img.shields.io/badge/docs-github.io-blue)](https://finite-sample.github.io/rank-preserving-calibration/)
[![PyPI Downloads](https://static.pepy.tech/badge/rank-preserving-calibration)](https://pepy.tech/projects/rank-preserving-calibration)
[![Python 3.12+](https://img.shields.io/badge/python-3.12+-blue.svg)](https://www.python.org/downloads/)
[![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT)

Survey statisticians and machine learning practitioners often need to adjust the predicted class probabilities from a classifier so that they match known population totals (column marginals). Simple post-hoc methods that apply separate logit shifts or raking to each class can scramble the ranking of individuals within a class when there are three or more classes. This package implements a rank-preserving calibration procedure that projects probabilities onto the intersection of two convex sets:

1. **Row-simplex**: each row sums to one and all entries are non-negative.
2. **Isotonic column marginals**: within each class, values are non-decreasing when instances are sorted by their original scores for that class, and the sum of each column equals a user-supplied target.

The algorithm uses Dykstra's alternating projection method in Euclidean geometry. When the specified column totals are feasible, the procedure returns a matrix that preserves cross-person discrimination within each class, matches the desired totals, and remains a valid probability distribution for each instance. That claim is checked against an independent convex solver in the test suite rather than asserted.

**Feasibility is exact, not approximate.** Every row sums to one, so the grand total is fixed at `N` and the targets must satisfy `sum(M) == N`. When they do, the intersection is *never* empty — the constant matrix `Q[i, j] = M[j] / N` satisfies all three constraint sets — so a failure is a conditioning or iteration-budget problem, never infeasibility. When they do not, no matrix satisfies both constraint sets and there is no "closest point satisfying both" to return: the solver warns, then raises `CalibrationError`. Rescale first:

```python
import numpy as np

N = 3
M = np.array([1.2, 0.9, 1.1])  # any targets
M = M * N / M.sum()  # now sums to N exactly
```

### New: Nearly Isotonic Calibration

This package now supports **nearly isotonic** constraints that allow small violations of strict monotonicity when appropriate:

- **Epsilon-slack constraints**: Allow z[i+1] ≥ z[i] - ε instead of strict z[i+1] ≥ z[i]
- **Penalised rank violations**: `calibrate_soft` (squared error) and `calibrate_kl_soft` (KL divergence) charge λ per unit of downward step instead of forbidding it

These relaxed constraints can provide better balance between rank preservation and probability calibration when strict isotonic constraints are too restrictive.

## Which solver?

`calibrate()` picks for you, and the default is right in every case measured. The two
are mathematically equivalent — they solve the same convex problem, and the test suite
asserts they agree — so the choice is practical:

| method | how | when |
|---|---|---|
| `"qp"` (default) | sparse quadratic program, interior point | always |
| `"dykstra"` | alternating projections in numpy | reference implementation the QP is tested against |

The default is not a close call. Measured at `J=4` on feasible targets:

| N | `dykstra` | `qp` |
|---|---|---|
| 25 | 0.03 s | 0.002 s |
| 50 | 9.6 s | 0.004 s |
| 100 | 10.5 s | 0.01 s |
| 400 | does not converge in 60,000 iterations | 0.12 s |
| 1600 | — | 1.4 s |
| 6400 | — | 50 s |

The reason is the algorithm class rather than the problem: the two constraint sets meet
at a shallow angle, which is exactly where first-order methods crawl and where an
interior-point method is unaffected. A sparse QP alone is not enough — OSQP, also
first-order, is *worse* than Dykstra here, returning 592 rank violations at `N=400`
while reporting that it ran.

**Above roughly `N = 6400` the exact solve becomes expensive again**, since
interior-point factorisation cost grows steeply. Relax the constraint rather than wait:

```python
import numpy as np
from rank_preserving_calibration import calibrate

P = np.array([[0.6, 0.3, 0.1], [0.2, 0.5, 0.3], [0.1, 0.2, 0.7]])
M = np.array([1.0, 1.0, 1.0])

result = calibrate(P, M, nearly={"mode": "epsilon", "eps": 0.05})
```

which permits adjacent within-class decreases of at most `eps` and is handled exactly,
by lowering the isotonic bound rather than by penalising violations.

## Installation

```bash
pip install rank-preserving-calibration
```

Runtime dependencies are `numpy`, `scipy` and `clarabel` (Apache-2.0, ~2.5 MB). The
examples also need `matplotlib` and `scikit-learn`.

## Usage

### Basic Usage

```python
import numpy as np
from rank_preserving_calibration import calibrate

P = np.array(
    [
        [0.6, 0.3, 0.1],
        [0.2, 0.5, 0.3],
        [0.1, 0.2, 0.7],
    ]
)

# Target column sums, e.g. population class frequencies. Must sum to the
# number of rows (3 in this example) for perfect feasibility.
M = np.array([1.0, 1.0, 1.0])

result = calibrate(P, M)

print("Adjusted probabilities:\n", result.Q)
print("Converged:", result.converged)
print("Iterations:", result.iterations)
print("Max row error:", result.max_row_error)
print("Max column error:", result.max_col_error)
print("Rank violations:", result.max_rank_violation)
```

### Nearly Isotonic Usage

```python
from rank_preserving_calibration import calibrate

# Epsilon-slack: allow small rank violations (recommended above N ~ 6400)
result = calibrate(P, M, nearly={"mode": "epsilon", "eps": 0.05})

# Penalise rank violations instead of bounding them. This changes the
# objective, not the constraint set, so it lives in its own function.
from rank_preserving_calibration import calibrate_soft

soft = calibrate_soft(P, M, lam_m=100.0, lam_r=1.0)
```

The returned `CalibrationResult` contains the calibrated matrix `Q` with the same shape as `P`. Each row of `Q` sums to one, the column sums match `M`, and within each column the entries are sorted in non-decreasing order according to the order implied by the original `P`.

## Evaluation and Metrics

After calibration, it's important to validate that the constraints are satisfied and understand the impact on prediction quality. This package provides comprehensive metrics for evaluation:

### Constraint Validation

```python
from rank_preserving_calibration import feasibility_metrics, isotonic_metrics

# Check constraint satisfaction
feasibility = feasibility_metrics(result.Q, M)
print(f"Max row error: {feasibility['row']['max_abs_error']}")
print(f"Max column error: {feasibility['col']['max_abs_error']}")

# Check rank preservation
isotonic = isotonic_metrics(result.Q, P)
print(f"Max rank violation: {isotonic['max_rank_violation']}")
print(f"Violation mass: {isotonic['total_violation_mass']}")
```

### Calibration Quality Assessment

```python
import numpy as np
from sklearn.metrics import log_loss

from rank_preserving_calibration import distance_metrics

# True class labels for the three rows above, if you have them.
y_true = np.array([0, 1, 2])

# Measure calibration changes
distances = distance_metrics(result.Q, P)
print(f"Frobenius distance: {distances['frobenius']}")
print(f"Max change: {distances['max_abs']}")

# Evaluate with labeled data (if available)
if y_true is not None:
    original_nll = log_loss(y_true, P)
    calibrated_nll = log_loss(y_true, result.Q)
    print(f"NLL improvement: {original_nll - calibrated_nll}")
```

### Available Metrics

| Function | Purpose |
| --- | --- |
| `feasibility_metrics(Q, M)` | Validate row (simplex) and column (marginal) constraints |
| `isotonic_metrics(Q, P)` | Check rank preservation and measure violations |
| `distance_metrics(Q, P)` | Quantify changes between original and calibrated probabilities |
| `tie_group_variance(Q, P)` | Assess handling of tied predictions (useful for `ties='group'`) |
| `top_label_ece(y, probs)` | Expected calibration error for top predictions |
| `classwise_ece(y, probs)` | Per-class calibration error analysis |
| `sharpness_metrics(probs)` | Prediction confidence and entropy analysis |
| `auc_deltas(y, P, Q)` | One-vs-rest AUC changes after calibration |

For log loss and the Brier score, use `sklearn.metrics.log_loss` and
`sklearn.metrics.brier_score_loss`; this package only adds the metrics
that are specific to calibration.

### Complete Evaluation Workflow

```python
import numpy as np
from rank_preserving_calibration import (
    calibrate_dykstra,
    feasibility_metrics,
    isotonic_metrics,
    distance_metrics,
    top_label_ece,
)

# Calibrate
result = calibrate_dykstra(P, M)

# 1. Validate constraints
feasibility = feasibility_metrics(result.Q, M)
isotonic = isotonic_metrics(result.Q, P)
print(f"Converged: {result.converged}")
print(f"Row constraint satisfied: {feasibility['row']['max_abs_error'] < 1e-6}")
print(f"Rank preserved: {isotonic['max_rank_violation'] < 1e-6}")

# 2. Assess calibration impact
distances = distance_metrics(result.Q, P)
print(f"Average change per probability: {distances['mean_abs']:.4f}")

# 3. Evaluate predictive quality (if labels available)
if y_true is not None:
    ece_before = top_label_ece(y_true, P)
    ece_after = top_label_ece(y_true, result.Q)
    print(f"Calibration error before: {ece_before['ece']:.3f}")
    print(f"Calibration error after: {ece_after['ece']:.3f}")
```

## Functions

### `calibrate_dykstra(P, M, **kwargs)`

Calibrate using Dykstra's alternating projections (recommended). Supports both strict and nearly isotonic constraints.

### `calibrate_soft(P, M, lam_m, lam_r)`

Minimise `0.5||Q - P||² + 0.5·lam_m·||Qᵀ1 - M||² + lam_r·Σ(downward steps)` with rows on the simplex, as a quadratic program solved exactly by Clarabel.

### `calibrate_kl(P, M, R=None, A=None)` and `calibrate_kl_soft(P, M, lam, R=None, A=None)`

Minimise KL divergence from a reference `R` with ranks taken from an anchor `A`, either strictly or with a λ penalty; both are solved exactly as exponential-cone programs.

### `create_test_case(case_type, N, J, **kwargs)` (in `tests.data_helpers`)

Generate synthetic test data for various scenarios used in testing.

## Arguments

| Parameter | Type | Description |
| --- | --- | --- |
| `P` | `ndarray` of shape `[N, J]` | Base multiclass probabilities or non-negative scores. Rows will be projected to the simplex. |
| `M` | `ndarray` of shape `[J]` | Target column totals (e.g. population class frequencies). The sum of `M` should equal the number of rows `N` for exact feasibility. |
| `max_iters` | `int` | Maximum number of projection iterations (Dykstra, default `3000`). |
| `tol` | `float` | Relative convergence tolerance (Dykstra, default `1e-7`). |
| `verbose` | `bool` | If `True`, prints convergence diagnostics. |
| `nearly` | `dict` | Nearly isotonic slack: `{"mode": "epsilon", "eps": 0.05}`. Any other mode raises. |

## Returns

### CalibrationResult

`calibrate`, `calibrate_qp` and `calibrate_dykstra` return a `CalibrationResult` object with the following attributes:

* `Q`: NumPy array of shape `[N, J]` containing the calibrated probabilities. Each row sums to one, each column approximately sums to the corresponding entry of `M`, and within each column the values are non-decreasing according to the ordering induced by `P`.
* `converged`: boolean indicating whether the solver met the tolerance criteria.
* `iterations`: number of iterations performed.
* `max_row_error`: maximum absolute deviation of row sums from 1.
* `max_col_error`: maximum absolute deviation of column sums from `M`.
* `max_rank_violation`: maximum violation of monotonicity (should be 0 up to numerical tolerance).
* `final_change`: final relative change between iterations.

## Algorithm Notes

* **Dykstra's Method**: Uses alternating projections with memory terms to ensure convergence to the intersection of constraint sets. Rows are projected onto the simplex via the algorithm of Duchi et al., and columns are projected via the pool-adjacent-violators algorithm followed by an additive shift to match column totals. This is the recommended method for most applications.

* **Nearly Isotonic Extensions**:
  - **Epsilon-slack (Dykstra)**: Projects onto the convex set {z : z[i+1] ≥ z[i] - ε} using coordinate transformation. Maintains theoretical convergence guarantees.
  - **Rank penalty (`calibrate_soft`)**: Minimises 0.5||Q - P||² + λ∑max(0, z[i] - z[i+1]) per column, alongside a marginal penalty, as one quadratic program solved by an interior-point method. Checked against an independent convex solver.

## Examples

The examples in [`docs/examples`](docs/examples) are scripts that run on every docs build and are rendered, with their output and figures, in the [documentation](https://finite-sample.github.io/rank-preserving-calibration/). Run one locally with `uv run --group docs python docs/examples/plot_label_shift.py`.

### When to Use Nearly Isotonic Calibration

**Use Nearly Isotonic When:**
- Model predictions have good discrimination but need marginal calibration
- Some predictions are already well-calibrated
- Small rank violations are acceptable in your domain
- You want to preserve model confidence where possible

**Use Strict Isotonic When:**
- Rank order is critical (regulatory, safety applications)
- Model predictions have clear monotonic relationship
- Conservative approach is preferred

## Testing

```bash
python -m pytest tests/ -v
```

## License

This software is released under the terms of the MIT license.

<!-- adjacent:start -->

## 🔗 Adjacent Repositories

- [finite-sample/calibre](https://github.com/finite-sample/calibre) — Advanced Calibration Models
- [finite-sample/optimal-classification-cutoffs](https://github.com/finite-sample/optimal-classification-cutoffs) — Script for calculating the optimal cut-off for max. F1-score, etc.
- [finite-sample/fairlex](https://github.com/finite-sample/fairlex) — Leximin Calibration
- [finite-sample/pyppur](https://github.com/finite-sample/pyppur) — pyppur: Python Projection Pursuit Unsupervised (Dimension) Reduction To Min. Reconstruction Loss or DIstance DIstortion

<!-- adjacent:end -->
