API Reference¶

This section contains the complete API documentation for the rank_preserving_calibration package.

Main Calibration Functions¶

Robust rank-preserving multiclass probability calibration.

This module provides numerically stable implementations of rank-preserving calibration by Dykstra’s alternating projections.

rank_preserving_calibration.calibration.calibrate_dykstra(P: ndarray, M: ndarray, max_iters: int = 3000, tol: float = 1e-07, feasibility_tol: float = 0.1, verbose: bool = False, callback: CallbackFunction = None, detect_cycles: bool = False, cycle_window: int = 10, nearly: dict | None = None, ties: str = 'stable', row_atol: float = 1e-08, col_atol: float = 1e-08) → CalibrationResult[source]¶

Calibrate using Dykstra’s alternating projections.

Projects multiclass probabilities onto the intersection of:

row simplex: {rows ≥ 0, rows sum to 1} and
column-wise isotone-by-score + fixed column sums: {nondecreasing in score order; column sum = M_j}.

This is the recommended default method for rank-preserving calibration. The algorithm uses exact Euclidean projections via Pool Adjacent Violators (PAV) followed by uniform shifts to satisfy sum constraints.

Parameters:

P – Input probability matrix of shape (N, J). Each row represents predicted class probabilities for one instance. Rows need not sum to 1 initially.
M – Target column sums of shape (J,). Should sum to approximately N for feasibility.
max_iters – Maximum number of iterations. Default 3000 is usually sufficient.
tol – Convergence tolerance for relative change in solution. Default 1e-7.
feasibility_tol – Tolerance for feasibility warnings when sum(M) differs from N.
verbose – If True, enables debug logging.
callback – Optional function called each iteration as callback(iter, change, Q). Should return False to terminate early.
detect_cycles – If True, detects and breaks cycles in the solution sequence.
cycle_window – Number of iterations to look back for cycle detection.
nearly – Optional dict for nearly-isotonic constraints. Use {“mode”: “epsilon”, “eps”: 0.01} to allow small isotonicity violations.
ties – How to handle tied scores. “stable” preserves input order, “group” pools equal-score instances.
row_atol – Absolute tolerance on |row sum - 1| for convergence. Absolute, not relative – see is_feasible. The 1e-8 default is what alternating projections can actually reach: each sweep ends on a column projection, which perturbs the row sums, so no finite iterate satisfies both constraint sets exactly. It is still three orders tighter than the 1e-5 the previous np.allclose check allowed.
col_atol – Absolute tolerance on |column sum - M_j| for declaring convergence. Being absolute, it does not loosen as M_j grows with N, which the previous check did.

Returns:

Q: Calibrated probability matrix of shape (N, J)
converged: Always True (failures raise CalibrationError instead)
iterations: Number of iterations performed
max_row_error: Maximum absolute row sum error
max_col_error: Maximum absolute column sum error
max_rank_violation: Maximum rank order violation
final_change: Final relative change in solution

Return type:

CalibrationResult object containing

Raises:

CalibrationError – If inputs are invalid, the algorithm fails to converge, or other errors occur.
ValueError – If ties parameter is not “stable” or “group”

Examples

Basic calibration:

>>> import numpy as np
>>> from rank_preserving_calibration import calibrate_dykstra
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> M = np.array([0.8, 0.8, 0.4])  # target column sums; P's are 1.0, 0.7, 0.3
>>> result = calibrate_dykstra(P, M)
>>> print(f"Converged: {result.converged}")
Converged: True
>>> result.Q.sum(axis=1).round(6).tolist()
[1.0, 1.0]
>>> result.Q.sum(axis=0).round(6).tolist()
[0.8, 0.8, 0.4]

With nearly-isotonic constraints:

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

Notes

Converges to the exact intersection of the constraint sets
Preserves ranking within each class (column) by original model scores
Memory complexity is O(N*J) for the probability matrices
Time complexity per iteration is O(N*J*log(N)) due to sorting
Raises CalibrationError on convergence failure rather than returning an unreliable result

Result Classes¶

class rank_preserving_calibration.calibration.CalibrationResult(Q: ndarray, converged: bool, iterations: int, max_row_error: float, max_col_error: float, max_rank_violation: float, final_change: float)[source]¶

Bases: object

Result container for rank-preserving calibration algorithms.

Returned by calibrate_dykstra() containing the calibrated probability matrix and diagnostic information about convergence and constraint satisfaction.

Variables:

Q (numpy.ndarray) – Calibrated probability matrix of shape (N, J) where rows sum to 1 and columns preserve rank ordering from original scores.
converged (bool) – True if algorithm converged within specified tolerance.
iterations (int) – Number of iterations performed before termination.
max_row_error (float) – Maximum absolute error in row sum constraint (should be ≈0).
max_col_error (float) – Maximum absolute error in column sum constraint.
max_rank_violation (float) – Maximum rank-order violation across all columns.
final_change (float) – Final relative change in solution between iterations.

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import calibrate_dykstra
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> M = np.array([0.8, 0.8, 0.4])
>>> result = calibrate_dykstra(P, M)
>>> if result.converged:
...     print(f"Calibration successful in {result.iterations} iterations")
Calibration successful in 2 iterations
>>> print(f"Max rank violation: {result.max_rank_violation:.6f}")
Max rank violation: 0.000000

Q: ndarray¶

converged: bool¶

iterations: int¶

max_row_error: float¶

max_col_error: float¶

max_rank_violation: float¶

final_change: float¶

Nearly Isotonic Functions¶

Nearly isotonic regression utilities.

This module provides “relaxed” isotonic constraints that allow small violations of monotonicity—useful when strict isotonicity is too restrictive.

Exports:

project_near_isotonic_euclidean: ε-slack projection (exact L2), with an optional sum target.
prox_near_isotonic: exact proximal operator for λ * sum (z_i - z_{i+1})_+.
prox_near_isotonic_with_sum: same prox with an exact sum constraint, imposed by translation.

rank_preserving_calibration.nearly.project_near_isotonic_euclidean(v: ndarray, eps: float, sum_target: float | None = None, weights: ndarray | None = None) → ndarray[source]¶

Project v onto the set { z : z_{i+1} >= z_i - eps } in L2.

Reduction: define w_i = v_i + i * eps. Then the constraint becomes standard isotonic on w. Let w* = isotonic(w) (weighted if weights are provided). The exact projection is z* = w* - i * eps. If sum_target is given, apply a uniform shift so 1^T z* = sum_target. That is exact for uniform weights only: the sum constraint’s multiplier shifts entry i by mu / w_i.

Parameters:

v – Input vector of shape (n,) to project.
eps – Slack parameter (>= 0). Allows z[i+1] >= z[i] - eps rather than the strict z[i+1] >= z[i].
sum_target – If provided, apply an exact uniform shift so the result sums to this value. Requires uniform weights.
weights – Positive weights of shape (n,) for a weighted L2 projection; cannot be combined with sum_target unless uniform.

Returns:

Projected vector of shape (n,) satisfying the near-isotonic constraint, and the sum constraint if one was requested.

Raises:

ValueError – If eps is negative, the weights are invalid, or sum_target is combined with non-uniform weights.

Notes

Idempotent: applying the projection twice gives the same result.

rank_preserving_calibration.nearly.prox_near_isotonic(y: ndarray, lam: float) → ndarray[source]¶

Compute the prox of the λ-penalty nearly-isotonic term exactly.

Solves:

minimize_z  0.5 * ||z - y||^2 + lam * sum_i (z_i - z_{i+1})_+

which charges every downward step and leaves upward ones free; lam = 0 returns y, and lam large enough returns the isotonic regression of y.

By dynamic programming along the chain, as in Johnson’s (2013) algorithm for the fused lasso. With f_i the optimal cost of the first i terms as a function of z_i, the message passed on is m_i(b) = min_a f_i(a) + lam * (a - b)_+, whose derivative is f_i’ clipped to [-lam, 0] (the fused lasso’s symmetric penalty clips to [-lam, lam]). f_i’ is piecewise linear, so its knots are kept in a deque and each step costs amortised O(1). Going back, z_i is z_{i+1} clipped to the points where f_i’ equals -lam and 0.

No maintained Python package solves this problem; the R package neariso (Hoefling, Tibshirani & Tibshirani) implements their path algorithm.

Parameters:

y – Input vector of shape (n,).
lam – Penalty weight on downward steps, non-negative.

Returns:

The minimiser z, of shape (n,).

Raises:

ValueError – If lam is negative.

rank_preserving_calibration.nearly.prox_near_isotonic_with_sum(y: ndarray, lam: float, sum_target: float) → ndarray[source]¶

Compute the prox under a sum constraint.

Solves:

minimize_z  0.5 * ||z - y||^2 + lam * sum_i (z_i - z_{i+1})_+
subject to  sum(z) = sum_target

The penalty depends only on differences, so it is translation-invariant, and the constrained prox is the unconstrained one shifted by a constant.

Parameters:

y – Input vector of shape (n,).
lam – Penalty weight on downward steps.
sum_target – Value the result must sum to.

Returns:

The constrained minimiser.

Evaluation and Metrics¶

Calibration quality metrics: feasibility, isotonicity, distances, scoring.

rank_preserving_calibration.metrics.kl_divergence(Q: ndarray, R: ndarray, eps: float = 1e-15) → float[source]¶

Compute KL divergence KL(Q||R) = sum Q * log(Q/R).

The KL divergence measures how much Q differs from R, where R is the reference distribution. Uses convention 0 * log(0/x) = 0.

Parameters:

Q – Probability matrix of shape (N, J).
R – Reference probability matrix of shape (N, J).
eps – Small constant for numerical stability.

Returns:

KL divergence value (non-negative).

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import kl_divergence
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> Q = np.array([[0.6, 0.3, 0.1], [0.4, 0.4, 0.2]])
>>> kl = kl_divergence(Q, P)
>>> print(f"KL divergence: {kl:.4f}")
KL divergence: 0.0550

rank_preserving_calibration.metrics.reverse_kl_divergence(Q: ndarray, R: ndarray, eps: float = 1e-15) → float[source]¶

Compute reverse KL divergence KL(R||Q) = sum R * log(R/Q).

The reverse KL divergence has different properties than forward KL: - Forward KL(Q||R) is mean-seeking (Q spreads over modes of R) - Reverse KL(R||Q) is mode-seeking (Q focuses on main mode of R)

Parameters:

Q – Probability matrix of shape (N, J).
R – Reference probability matrix of shape (N, J).
eps – Small constant for numerical stability.

Returns:

Reverse KL divergence value (non-negative).

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import reverse_kl_divergence
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> Q = np.array([[0.6, 0.3, 0.1], [0.4, 0.4, 0.2]])
>>> rkl = reverse_kl_divergence(Q, P)

rank_preserving_calibration.metrics.symmetrized_kl(Q: ndarray, R: ndarray, eps: float = 1e-15) → float[source]¶

Compute symmetrized KL divergence (Jensen-Shannon related).

Symmetrized KL = 0.5 * (KL(Q||R) + KL(R||Q))

This is symmetric in Q and R, unlike the standard KL divergence.

Parameters:

Q – Probability matrix of shape (N, J).
R – Reference probability matrix of shape (N, J).
eps – Small constant for numerical stability.

Returns:

Symmetrized KL divergence value (non-negative).

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import symmetrized_kl
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> Q = np.array([[0.6, 0.3, 0.1], [0.4, 0.4, 0.2]])
>>> skl = symmetrized_kl(Q, P)

rank_preserving_calibration.metrics.feasibility_metrics(Q: ndarray, M: ndarray | None = None) → dict[str, Any][source]¶

Compute feasibility metrics for calibrated probability matrix.

Analyzes how well the probability matrix Q satisfies the row and column sum constraints required by rank-preserving calibration.

Parameters:

Q – Calibrated probability matrix of shape (N, J). Should have rows summing to 1 and non-negative entries.
M – Target column sums of shape (J,). If provided, analyzes column sum constraint satisfaction.

Returns:

“row”: Row constraint metrics – max/mean absolute errors and the minimum value
”col”: Column constraint metrics (if M provided) – various error norms

Return type:

Dictionary containing feasibility metrics

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import feasibility_metrics
>>> Q = np.array([[0.6, 0.3, 0.1], [0.4, 0.4, 0.2]])
>>> M = np.array([1.0, 0.7, 0.3])
>>> metrics = feasibility_metrics(Q, M)
>>> print(f"Max row error: {metrics['row']['max_abs_error']:.6f}")
Max row error: 0.000000
>>> print(f"Max col error: {metrics['col']['max_abs_error']:.6f}")
Max col error: 0.000000

rank_preserving_calibration.metrics.isotonic_metrics(Q: ndarray, P: ndarray) → dict[str, Any][source]¶

Analyze isotonic (rank-preserving) properties of calibrated probabilities.

Verifies that calibrated probabilities are nondecreasing when ordered by original model scores within each class. Computes violation statistics and flatness measures.

Parameters:

Q – Calibrated probability matrix of shape (N, J).
P – Original probability matrix of shape (N, J) used for score ordering.

Returns:

“max_rank_violation”: Maximum rank-order violation across all columns
”mass_weighted_violation”: Sum of violation magnitudes weighted by mass
”flat_fraction”: Fraction of adjacent pairs that are equal (flatness)
”per_class_max_violation”: List of max violations per class
”per_class_mass_violation”: List of mass-weighted violations per class

Return type:

Dictionary containing isotonic metrics

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import isotonic_metrics
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> Q = np.array([[0.65, 0.25, 0.1], [0.35, 0.45, 0.2]])
>>> metrics = isotonic_metrics(Q, P)
>>> print(f"Max violation: {metrics['max_rank_violation']:.6f}")
Max violation: 0.000000

rank_preserving_calibration.metrics.tie_group_variance(Q: ndarray, P: ndarray) → dict[str, Any][source]¶: Within-equal-score group variance of Q; for ties=’group’ this should be ~0.

rank_preserving_calibration.metrics.distance_metrics(Q: ndarray, P: ndarray) → dict[str, float][source]¶

Compute distance metrics between original and calibrated probabilities.

Measures how much the calibration procedure has changed the probability estimates using various matrix norms and distance measures.

Parameters:

Q – Calibrated probability matrix of shape (N, J).
P – Original probability matrix of shape (N, J).

Returns:

“frobenius”: Frobenius norm ||Q - P||_F
”frobenius_sq”: Squared Frobenius norm ||Q - P||²_F
”mean_abs”: Mean absolute difference
”max_abs”: Maximum absolute difference

Return type:

Dictionary containing distance metrics

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import distance_metrics
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> Q = np.array([[0.65, 0.25, 0.1], [0.35, 0.45, 0.2]])
>>> distances = distance_metrics(Q, P)
>>> print(f"Frobenius distance: {distances['frobenius']:.4f}")
Frobenius distance: 0.1000

rank_preserving_calibration.metrics.top_label_ece(y: ndarray, probs: ndarray, n_bins: int = 15) → dict[source]¶

Compute Expected Calibration Error for top-label predictions.

Bins predictions by confidence (max probability) and measures the gap between confidence and accuracy in each bin. ECE is the weighted average of these gaps, and MCE is the maximum gap.

Parameters:

y – True class labels as integers of shape (N,).
probs – Probability matrix of shape (N, J).
n_bins – Number of confidence bins (default 15).

Returns:

“ece”: Expected Calibration Error (weighted average gap)
”mce”: Maximum Calibration Error (max gap across bins)
”bins”: List of (count, mean_confidence, accuracy) per bin

Return type:

Dictionary containing

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import top_label_ece
>>> y = np.array([0, 1, 2])
>>> probs = np.array([[0.8, 0.1, 0.1], [0.2, 0.7, 0.1], [0.1, 0.2, 0.7]])
>>> ece_result = top_label_ece(y, probs)
>>> print(f"ECE: {ece_result['ece']:.4f}")
ECE: 0.2667

rank_preserving_calibration.metrics.classwise_ece(y: ndarray, probs: ndarray, n_bins: int = 15) → dict[source]¶

Compute class-wise Expected Calibration Error.

For each class j, bins rows by their predicted probability of j and compares the mean prediction in each bin with the share of rows in it that are class j, weighting each bin by its share of all rows. That is the class-j ECE; the overall figure is its mean over classes (Kull et al., 2019, eq. 4).

Parameters:

y – True class labels as integers of shape (N,).
probs – Probability matrix of shape (N, J).
n_bins – Number of equal-width probability bins per class (default 15).

Returns:

“ece”: Mean of the class-j ECEs
”mce”: Largest bin gap in any class
”per_class”: Per-class ECE, MCE and (count, mean prediction, frequency) bin tables

Return type:

Dictionary containing

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import classwise_ece
>>> y = np.array([0, 1, 2, 2])
>>> probs = np.array([[0.7, 0.2, 0.1], [0.2, 0.7, 0.1],
...                   [0.2, 0.2, 0.6], [0.6, 0.2, 0.2]])
>>> round(classwise_ece(y, probs, n_bins=2)["ece"], 4)
0.2167

rank_preserving_calibration.metrics.sharpness_metrics(probs: ndarray) → dict[str, float][source]¶

Compute sharpness metrics for probability predictions.

Sharpness measures how confident (peaked) the predictions are, independent of whether they are correct. Sharper predictions are more informative.

Parameters:

probs – Probability matrix of shape (N, J).

Returns:

“mean_entropy”: Average entropy across predictions (lower = sharper)
”mean_max_prob”: Average of maximum probabilities (higher = sharper)
”mean_margin”: Average difference between top two probabilities

Return type:

Dictionary containing

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import sharpness_metrics
>>> probs = np.array([[0.9, 0.05, 0.05], [0.5, 0.3, 0.2]])
>>> sharp = sharpness_metrics(probs)
>>> print(f"Mean max prob: {sharp['mean_max_prob']:.3f}")
Mean max prob: 0.700

rank_preserving_calibration.metrics.column_variance(Q: ndarray) → dict[str, Any][source]¶

Compute per-column variance of a probability matrix.

Useful for measuring how “spread out” the calibrated probabilities are within each class. Lower variance may indicate a “flatter” solution.

Parameters:

Q – Probability matrix of shape (N, J).

Returns:

“per_column”: List of variance values for each column
”mean”: Mean variance across columns
”min”: Minimum column variance
”max”: Maximum column variance

Return type:

Dictionary containing

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import column_variance
>>> Q = np.array([[0.6, 0.3, 0.1], [0.4, 0.4, 0.2]])
>>> var = column_variance(Q)
>>> print(f"Mean column variance: {var['mean']:.4f}")
Mean column variance: 0.0050

rank_preserving_calibration.metrics.informativeness_ratio(Q: ndarray, P: ndarray) → dict[str, Any][source]¶

Compare variance of Q vs P to measure information preservation.

The informativeness ratio measures how much of the original probability spread is preserved after calibration. A ratio close to 1 means the calibration preserved the discriminative structure; a ratio near 0 indicates a “flat” solution with little discrimination between samples.

Parameters:

Q – Calibrated probability matrix of shape (N, J).
P – Original probability matrix of shape (N, J).

Returns:

“total_ratio”: Var(Q) / Var(P)
”per_column_ratio”: Variance ratio for each column
”mean_column_ratio”: Mean of per-column ratios
”interpretation”: String describing the result

Return type:

Dictionary containing

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import informativeness_ratio
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> Q = np.array([[0.5, 0.35, 0.15], [0.5, 0.35, 0.15]])  # Flat
>>> ratio = informativeness_ratio(Q, P)
>>> print(f"Total ratio: {ratio['total_ratio']:.3f}")
Total ratio: 0.487

rank_preserving_calibration.metrics.auc_deltas(y: ndarray, P: ndarray, Q: ndarray) → dict[str, Any][source]¶: One-vs-rest AUC before vs after; useful to check rank effects of plateaus.

Exception Classes¶

class rank_preserving_calibration.calibration.CalibrationError[source]¶

Bases: Exception

Raised when calibration fails due to invalid inputs or numerical issues.

Complete Module Documentation¶

Rank-preserving calibration of multiclass probabilities.

This package provides robust implementations of rank-preserving calibration algorithms: an exact sparse QP, Dykstra’s alternating projections, KL-divergence variants, and soft-penalty trade-offs.

Quick start¶

>>> import numpy as np
>>> from rank_preserving_calibration import calibrate_dykstra
>>> from rank_preserving_calibration import feasibility_metrics, isotonic_metrics
>>> # Toy data
>>> rng = np.random.default_rng(42)
>>> P = rng.dirichlet(np.ones(4), size=100)  # N x J predicted probs, rows sum to 1
>>> # Target class counts: shifted, but still summing to N, since each row sums to 1
>>> M = P.sum(axis=0) * np.array([1.2, 1.0, 0.9, 0.9])
>>> M *= len(P) / M.sum()
>>> # Calibrate
>>> res = calibrate_dykstra(P, M)
>>> res.converged
True
>>> # Check invariants: rows still sum to 1, columns hit M, order is kept
>>> fm = feasibility_metrics(res.Q, M)
>>> bool(fm["row"]["max_abs_error"] < 1e-6 and fm["col"]["max_abs_error"] < 1e-6)
True
>>> isotonic_metrics(res.Q, P)["max_rank_violation"]
0.0

exception rank_preserving_calibration.CalibrationError[source]

Bases: Exception

Raised when calibration fails due to invalid inputs or numerical issues.

class rank_preserving_calibration.CalibrationResult(Q: ndarray, converged: bool, iterations: int, max_row_error: float, max_col_error: float, max_rank_violation: float, final_change: float)[source]

Bases: object

Result container for rank-preserving calibration algorithms.

Returned by calibrate_dykstra() containing the calibrated probability matrix and diagnostic information about convergence and constraint satisfaction.

Variables:

Q (numpy.ndarray) – Calibrated probability matrix of shape (N, J) where rows sum to 1 and columns preserve rank ordering from original scores.
converged (bool) – True if algorithm converged within specified tolerance.
iterations (int) – Number of iterations performed before termination.
max_row_error (float) – Maximum absolute error in row sum constraint (should be ≈0).
max_col_error (float) – Maximum absolute error in column sum constraint.
max_rank_violation (float) – Maximum rank-order violation across all columns.
final_change (float) – Final relative change in solution between iterations.

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import calibrate_dykstra
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> M = np.array([0.8, 0.8, 0.4])
>>> result = calibrate_dykstra(P, M)
>>> if result.converged:
...     print(f"Calibration successful in {result.iterations} iterations")
Calibration successful in 2 iterations
>>> print(f"Max rank violation: {result.max_rank_violation:.6f}")
Max rank violation: 0.000000

Q: ndarray

converged: bool

iterations: int

max_row_error: float

max_col_error: float

max_rank_violation: float

final_change: float

class rank_preserving_calibration.IPFResult(Q: ndarray, converged: bool, iterations: int, max_row_error: float, max_col_error: float, final_change: float)[source]

Bases: object

Result from Iterative Proportional Fitting.

Variables:

Q (numpy.ndarray) – Probability matrix after IPF.
converged (bool) – Whether both the relative change and the returned matrix’s largest column error fell below tol.
iterations (int) – Number of iterations.
max_row_error (float) – Maximum row sum error.
max_col_error (float) – Maximum column sum error.
final_change (float) – Final relative change between iterations.

Q: ndarray

converged: bool

iterations: int

max_row_error: float

max_col_error: float

final_change: float

class rank_preserving_calibration.KLCalibrationResult(Q: ndarray, converged: bool, iterations: int, kl_divergence: float, max_row_error: float, max_col_error: float, max_rank_violation: float, final_change: float)[source]

Bases: object

Result container for KL-divergence rank-preserving calibration.

Variables:

Q (numpy.ndarray) – Calibrated probability matrix of shape (N, J) where rows sum to 1 and columns preserve rank ordering from anchor scores.
converged (bool) – True if the row and column errors are both at most 1e-7.
iterations (int) – Number of iterations performed before termination.
kl_divergence (float) – Final KL(Q||R) divergence value.
max_row_error (float) – Maximum absolute error in row sum constraint.
max_col_error (float) – Maximum absolute error in column sum constraint.
max_rank_violation (float) – Maximum rank-order violation across all columns.
final_change (float) – Final relative change in solution between iterations.

Q: ndarray

converged: bool

iterations: int

kl_divergence: float

max_row_error: float

max_col_error: float

max_rank_violation: float

final_change: float

class rank_preserving_calibration.KLParetoResult(solutions: list[tuple[float, ~numpy.ndarray]]=<factory>, kl_values: list[float] = <factory>, rank_violations: list[float] = <factory>, lambda_path: list[float] = <factory>)[source]

Bases: object

Result container for KL Pareto frontier computation.

Variables:

solutions (list[tuple[float, numpy.ndarray]]) – List of (lambda, Q) pairs along the Pareto frontier.
kl_values (list[float]) – KL divergence values for each solution.
rank_violations (list[float]) – Total rank violation for each solution.
lambda_path (list[float]) – Lambda values used in the sweep.

solutions: list[tuple[float, ndarray]]

kl_values: list[float]

rank_violations: list[float]

lambda_path: list[float]

class rank_preserving_calibration.SoftCalibrationResult(Q: ndarray, converged: bool, iterations: int, objective: float, fit_term: float, marginal_term: float, rank_term: float, max_row_error: float, max_col_error: float, max_rank_violation: float)[source]

Bases: object

Result from soft-constraint calibration.

Variables:

Q (numpy.ndarray) – Calibrated probability matrix.
converged (bool) – Whether the solver solved the program and every row of Q is a probability distribution to within 1e-7.
iterations (int) – Interior-point iterations.
objective (float) – The minimised objective at Q.
fit_term (float) – ||Q - P||^2 at Q.
marginal_term (float) – ||col_sums(Q) - M||^2 at Q.
rank_term (float) – Summed downward steps at Q.
max_row_error (float) – Maximum row sum error.
max_col_error (float) – Maximum column sum error.
max_rank_violation (float) – Maximum rank violation.

Q: ndarray

converged: bool

iterations: int

objective: float

fit_term: float

marginal_term: float

rank_term: float

max_row_error: float

max_col_error: float

max_rank_violation: float

class rank_preserving_calibration.TwoStageResult(Q: ndarray, converged: bool, ipf_iterations: int, projection_iterations: int, max_row_error: float, max_col_error: float, max_rank_violation: float, ipf_result: ndarray)[source]

Bases: object

Result from two-stage IPF-based calibration.

Variables:

Q (numpy.ndarray) – Calibrated probability matrix.
converged (bool) – Whether the algorithm converged.
ipf_iterations (int) – Number of IPF iterations performed.
projection_iterations (int) – Number of isotonic projection iterations.
max_row_error (float) – Maximum row sum error.
max_col_error (float) – Maximum column sum error.
max_rank_violation (float) – Maximum rank violation.
ipf_result (numpy.ndarray) – The intermediate IPF result before isotonic projection.

Q: ndarray

converged: bool

ipf_iterations: int

projection_iterations: int

max_row_error: float

max_col_error: float

max_rank_violation: float

ipf_result: ndarray

rank_preserving_calibration.auc_deltas(y: ndarray, P: ndarray, Q: ndarray) → dict[str, Any][source]: One-vs-rest AUC before vs after; useful to check rank effects of plateaus.

rank_preserving_calibration.calibrate(P: np.ndarray, M: np.ndarray, *, method: str = 'auto', nearly: dict | None = None, **kwargs) → CalibrationResult[source]

Calibrate a probability matrix, preserving within-class rankings.

Finds the matrix closest to P in Frobenius norm whose rows are probability distributions, whose column sums equal M, and whose columns are non-decreasing in the order of the original scores – so the ranking of individuals within each class survives calibration.

Parameters:

P – Probability matrix of shape (N, J). The within-column orderings of P are what get preserved.
M – Target column sums of shape (J,). Must sum to N; rows summing to 1 fix the grand total, so any other value is infeasible. Rescale with M = M * N / M.sum().
method – "auto" (default, currently "qp"), "qp" or "dykstra". See the module docstring for the trade-offs.
nearly – None for strict rank preservation, or {"mode": "epsilon", "eps": e} to permit adjacent decreases of at most e. Worth reaching for above the size where the exact solve becomes expensive.
**kwargs – Passed through to the chosen solver.

Returns:

CalibrationResult holding the calibrated matrix and its diagnostics.

Raises:

CalibrationError – If method is unknown, the inputs are invalid, or the chosen solver fails.

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import calibrate
>>> P = np.array([[0.7, 0.3], [0.4, 0.6], [0.5, 0.5]])
>>> M = np.array([1.8, 1.2])          # must sum to N = 3
>>> result = calibrate(P, M)
>>> result.converged
True

Ranking within each class is preserved exactly:

>>> float(result.max_rank_violation) < 1e-9
True

The solvers agree, because they solve the same problem:

>>> other = calibrate(P, M, method="dykstra")
>>> bool(np.allclose(result.Q, other.Q, atol=1e-6))
True

rank_preserving_calibration.calibrate_dykstra(P: ndarray, M: ndarray, max_iters: int = 3000, tol: float = 1e-07, feasibility_tol: float = 0.1, verbose: bool = False, callback: CallbackFunction = None, detect_cycles: bool = False, cycle_window: int = 10, nearly: dict | None = None, ties: str = 'stable', row_atol: float = 1e-08, col_atol: float = 1e-08) → CalibrationResult[source]

Calibrate using Dykstra’s alternating projections.

Projects multiclass probabilities onto the intersection of:

row simplex: {rows ≥ 0, rows sum to 1} and
column-wise isotone-by-score + fixed column sums: {nondecreasing in score order; column sum = M_j}.

Parameters:

P – Input probability matrix of shape (N, J). Each row represents predicted class probabilities for one instance. Rows need not sum to 1 initially.
M – Target column sums of shape (J,). Should sum to approximately N for feasibility.
max_iters – Maximum number of iterations. Default 3000 is usually sufficient.
tol – Convergence tolerance for relative change in solution. Default 1e-7.
feasibility_tol – Tolerance for feasibility warnings when sum(M) differs from N.
verbose – If True, enables debug logging.
callback – Optional function called each iteration as callback(iter, change, Q). Should return False to terminate early.
detect_cycles – If True, detects and breaks cycles in the solution sequence.
cycle_window – Number of iterations to look back for cycle detection.
nearly – Optional dict for nearly-isotonic constraints. Use {“mode”: “epsilon”, “eps”: 0.01} to allow small isotonicity violations.
ties – How to handle tied scores. “stable” preserves input order, “group” pools equal-score instances.
row_atol – Absolute tolerance on |row sum - 1| for convergence. Absolute, not relative – see is_feasible. The 1e-8 default is what alternating projections can actually reach: each sweep ends on a column projection, which perturbs the row sums, so no finite iterate satisfies both constraint sets exactly. It is still three orders tighter than the 1e-5 the previous np.allclose check allowed.
col_atol – Absolute tolerance on |column sum - M_j| for declaring convergence. Being absolute, it does not loosen as M_j grows with N, which the previous check did.

Returns:

Q: Calibrated probability matrix of shape (N, J)
converged: Always True (failures raise CalibrationError instead)
iterations: Number of iterations performed
max_row_error: Maximum absolute row sum error
max_col_error: Maximum absolute column sum error
max_rank_violation: Maximum rank order violation
final_change: Final relative change in solution

Return type:

CalibrationResult object containing

Raises:

CalibrationError – If inputs are invalid, the algorithm fails to converge, or other errors occur.
ValueError – If ties parameter is not “stable” or “group”

Examples

Basic calibration:

>>> import numpy as np
>>> from rank_preserving_calibration import calibrate_dykstra
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> M = np.array([0.8, 0.8, 0.4])  # target column sums; P's are 1.0, 0.7, 0.3
>>> result = calibrate_dykstra(P, M)
>>> print(f"Converged: {result.converged}")
Converged: True
>>> result.Q.sum(axis=1).round(6).tolist()
[1.0, 1.0]
>>> result.Q.sum(axis=0).round(6).tolist()
[0.8, 0.8, 0.4]

With nearly-isotonic constraints:

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

Notes

Converges to the exact intersection of the constraint sets
Preserves ranking within each class (column) by original model scores
Memory complexity is O(N*J) for the probability matrices
Time complexity per iteration is O(N*J*log(N)) due to sorting
Raises CalibrationError on convergence failure rather than returning an unreliable result

rank_preserving_calibration.calibrate_ipf(P: ndarray, M: ndarray, max_iters: int = 100, tol: float = 1e-08, verbose: bool = False) → IPFResult[source]

Iterative Proportional Fitting (raking) to match target marginals.

IPF alternately scales rows and columns to match target constraints: - Row scaling: Ensure rows sum to 1 (probability simplex) - Column scaling: Scale columns to match target marginals M

Unlike Dykstra, IPF preserves the relative structure of probabilities within rows (ratios between classes). This can produce less “flat” solutions but may not perfectly satisfy all constraints simultaneously.

Parameters:

P – Input probability matrix of shape (N, J).
M – Target column sums of shape (J,). Should sum to N.
max_iters – Maximum number of row/column scaling iterations.
tol – Convergence tolerance for relative change.
verbose – If True, enables debug logging.

Returns:

IPFResult with the scaled probability matrix.

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import calibrate_ipf
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> M = np.array([0.8, 0.8, 0.4])
>>> result = calibrate_ipf(P, M)
>>> result.Q.sum(axis=0).round(6).tolist()
[0.8, 0.8, 0.4]

rank_preserving_calibration.calibrate_kl(P: ndarray, M: ndarray, R: ndarray | None = None, A: ndarray | None = None, tol: float = 1e-10, verbose: bool = False) → KLCalibrationResult[source]

Calibrate by minimising KL divergence under rank preservation.

Solves exactly, as an exponential-cone program:

min KL(Q||R)  s.t.  rows of Q on the simplex, column sums = M,
                    each column non-decreasing in the order of A

The anchor-reference decoupling separates what divergence is measured from (R) and what fixes the rank order (A).

Parameters:

P – Input probability matrix of shape (N, J).
M – Target column sums of shape (J,). Must sum to N.
R – Reference distribution for KL divergence. If None, uses P.
A – Anchor for rank ordering. If None, uses P.
tol – Solver gap and feasibility tolerance.
verbose – Print the solver log.

Returns:

KLCalibrationResult with calibrated matrix and diagnostics.

Raises:

CalibrationError – If the inputs are invalid or the solver fails.

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import calibrate_kl
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> M = np.array([0.8, 0.8, 0.4])
>>> result = calibrate_kl(P, M)
>>> round(result.max_col_error, 6), result.max_rank_violation < 1e-8
(0.0, True)

With anchor-reference decoupling:

>>> A = np.array([[0.6, 0.3, 0.1], [0.4, 0.4, 0.2]])
>>> result = calibrate_kl(P, M, R=P, A=A)

rank_preserving_calibration.calibrate_kl_pareto(P: ndarray, M: ndarray, lambda_grid: list[float] | ndarray | None = None, R: ndarray | None = None, A: ndarray | None = None, tol: float = 1e-10, verbose: bool = False) → KLParetoResult[source]

Compute the frontier of KL divergence against rank violation.

Solves calibrate_kl_soft() at each λ, reporting the whole path rather than a single tuned point. Because each solve is exact, KL is non-decreasing and the rank penalty non-increasing along the path.

Parameters:

P – Input probability matrix of shape (N, J).
M – Target column sums of shape (J,).
lambda_grid – λ values to solve at. Default: 20 points from 0.01 to 1000.
R – Reference for KL divergence. Default: P.
A – Anchor for rank ordering. Default: P.
tol – Solver gap and feasibility tolerance.
verbose – Print the solver logs.

Returns:

KLParetoResult with solutions along the frontier, in increasing λ.

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import calibrate_kl_pareto
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> M = np.array([0.8, 0.8, 0.4])
>>> result = calibrate_kl_pareto(P, M)
>>> len(result.solutions) == len(result.kl_values) == len(result.lambda_path)
True

rank_preserving_calibration.calibrate_kl_soft(P: ndarray, M: ndarray, lam: float = 1.0, R: ndarray | None = None, A: ndarray | None = None, tol: float = 1e-10, verbose: bool = False) → KLCalibrationResult[source]

Soft KL calibration with a λ-weighted rank penalty.

Solves exactly, as an exponential-cone program:

min KL(Q||R) + λ·V_A(Q)  s.t.  Q ∈ C(M)

Where:

KL(Q||R) = sum q log(q / r) is the divergence from the reference R
V_A(Q) sums every adjacent downward step, (q_lower - q_upper)_+, within each column in the order of the anchor A
C(M) is rows on the simplex with column sums equal to M

At λ = 0 this is the KL projection onto C(M), which IPF also computes; for λ large enough the penalty is exact and the result equals calibrate_kl().

Parameters:

P – Input probability matrix of shape (N, J).
M – Target column sums of shape (J,). Must sum to N.
lam – Rank penalty weight. Larger = more isotonic.
R – Reference for KL divergence. Default: P.
A – Anchor for rank ordering. Default: P.
tol – Solver gap and feasibility tolerance.
verbose – Print the solver log.

Returns:

KLCalibrationResult with the calibrated matrix.

Raises:

CalibrationError – If the inputs are invalid, λ is negative, R or A has the wrong shape or R a negative entry, or the solver fails.

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import calibrate_kl_soft
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> M = np.array([0.8, 0.8, 0.4])
>>> strong = calibrate_kl_soft(P, M, lam=10.0)
>>> weak = calibrate_kl_soft(P, M, lam=0.0)
>>> bool(strong.max_rank_violation <= weak.max_rank_violation)
True

rank_preserving_calibration.calibrate_ovr_isotonic(y: ndarray, probs: ndarray) → dict[str, Any][source]

Calibrate multiclass probabilities by one-vs-rest isotonic regression.

Fits a separate isotonic regression per class on the binary problem of that class against the rest, then normalises each row to sum to 1. This is the conventional approach, used by scikit-learn among others, and is included here as the baseline the rank-preserving solvers are measured against.

Warning

This does not preserve rank, despite living in a package named for it. The final row normalisation divides every entry by a row-specific total, and two people with the same class-j score but different row totals come out in a different order than they went in. Measured on columns that are perfectly isotonic before normalisation, 57% of adjacent within-class pairs invert.

If you rank individuals by their probability of a given class, use calibrate_dykstra() instead, which projects onto the constraint set rather than rescaling rows.

The row normalisation is also not free statistically: it partially undoes the per-class calibration it just performed, which is why one-vs-rest isotonic ranks last of seven methods on NLL in Dimitriadis-style benchmark comparisons of multiclass calibrators.

Parameters:

y – True class labels as integers of shape (N,).
probs – Original probability matrix of shape (N, J).

Returns:

A dictionary containing the calibrated probabilities under key 'Q', with rows summing to 1.

Raises:

CalibrationError – If y is not N integer labels in [0, J).

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import calibrate_ovr_isotonic
>>> rng = np.random.default_rng(0)
>>> P = rng.dirichlet(np.ones(3), size=50)
>>> y = rng.integers(0, 3, 50)
>>> Q = calibrate_ovr_isotonic(y, P)["Q"]
>>> bool(np.allclose(Q.sum(axis=1), 1.0))
True

rank_preserving_calibration.calibrate_qp(P: ndarray, M: ndarray, *, nearly: dict | None = None, row_atol: float = 1e-08, col_atol: float = 1e-08, solver_tol: float = 1e-10, feasibility_tol: float = 0.1, verbose: bool = False) → CalibrationResult[source]

Calibrate by solving the projection exactly as a sparse QP.

This is the recommended solver. It minimises ||Q - P||_F^2 over matrices with unit row sums, column sums M, non-negative entries, and columns that are non-decreasing in the order of the original scores – the same problem calibrate_dykstra() solves by alternating projections, but by interior point, which does not stall on this geometry.

Parameters:

P – Probability matrix of shape (N, J). Rows need not already sum to 1; the within-column orderings of P are what get preserved.
M – Target column sums of shape (J,). Must sum to N: rows summing to 1 fix the grand total, so any other total makes the problem infeasible.
nearly – None for strict isotonicity, or {"mode": "epsilon", "eps": e} to permit adjacent decreases of at most e. Handled exactly, by lowering the isotonic bound rather than by penalising violations.
row_atol – Absolute tolerance on |row sum - 1| when reporting convergence.
col_atol – Absolute tolerance on |column sum - M_j|.
solver_tol – Interior-point convergence tolerance. The 1e-10 default is tighter than the solver’s own default because rank preservation is the guarantee here: at the looser default the result carries isotonic violations of order 1e-9 to 1e-7. The cost is negligible.
feasibility_tol – Retained for signature compatibility with the other solvers; feasibility itself is exact and is checked in validation.
verbose – Print solver progress.

Returns:

CalibrationResult with the calibrated matrix and its diagnostics.

Raises:

CalibrationError – If the inputs are invalid, the optional solver dependencies are missing, or the solver does not reach an optimum.

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import calibrate_qp
>>> P = np.array([[0.7, 0.3], [0.4, 0.6], [0.5, 0.5]])
>>> M = np.array([1.8, 1.2])          # sums to N = 3
>>> result = calibrate_qp(P, M)
>>> result.converged
True
>>> bool(np.allclose(result.Q.sum(axis=1), 1.0))
True
>>> bool(np.allclose(result.Q.sum(axis=0), M))
True

Ranking within each class is preserved exactly:

>>> float(result.max_rank_violation) < 1e-9
True

rank_preserving_calibration.calibrate_soft(P: ndarray, M: ndarray, lam_m: float = 1.0, lam_r: float = 10.0, *, tol: float = 1e-10, verbose: bool = False) → SoftCalibrationResult[source]

Calibrate with soft penalties on the marginals and on rank violations.

Solves the convex problem:

minimise  0.5 ||Q - P||^2 + 0.5 lam_m ||Q^T 1 - M||^2
          + lam_r * sum_j sum_i (q_(i),j - q_(i+1),j)_+
subject to each row of Q on the probability simplex,

where q_(i),j is column j in the order of P’s column j, so the last term charges every adjacent downward step. Rows stay exact probability distributions; the marginals and the ranking are traded off against staying close to P.

It is a quadratic program, solved by Clarabel’s interior-point method as calibrate_qp is: variables Q, the marginal gaps s = Q^T 1 - M, and a hinge v >= q_lower - q_upper, v >= 0 per adjacent pair. An ADMM with an exact rank prox reached the same optimum at N = 30 but stalled at its iteration cap from N = 300, 25-90x slower and up to 2e-2 above the optimum.

Parameters:

P – Input probability matrix of shape (N, J).
M – Target column sums of shape (J,).
lam_m – Marginal penalty weight. 0 ignores the marginals.
lam_r – Rank penalty weight. 0 ignores the ranking; large values approach strict isotonicity.
tol – Solver gap and feasibility tolerance.
verbose – Print the solver log.

Returns:

SoftCalibrationResult with the calibrated matrix.

Raises:

CalibrationError – If P or M are malformed, a weight is negative, or the solver fails.

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import calibrate_soft
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> M = np.array([0.8, 0.8, 0.4])
>>> result = calibrate_soft(P, M, lam_m=100.0, lam_r=1.0)
>>> result.converged
True
>>> np.round(result.Q.sum(axis=0), 2).tolist()
[0.8, 0.8, 0.4]

rank_preserving_calibration.calibrate_two_stage(P: ndarray, M: ndarray, ipf_max_iters: int = 100, ipf_tol: float = 1e-08, proj_max_iters: int = 100, proj_tol: float = 1e-08, preserve_marginals: bool = False, verbose: bool = False) → TwoStageResult[source]

Two-stage calibration: IPF followed by isotonic projection.

This approach first uses IPF (raking) to approximately match target marginals while preserving the relative probability structure, then applies isotonic projection to restore rank ordering.

The two-stage approach can produce less “flat” solutions compared to direct Dykstra projection when distribution shifts are large, because IPF preserves probability ratios rather than doing Euclidean projection.

Parameters:

P – Input probability matrix of shape (N, J).
M – Target column sums of shape (J,). Should sum to N.
ipf_max_iters – Maximum IPF iterations in stage 1.
ipf_tol – Convergence tolerance for IPF.
proj_max_iters – Maximum iterations for isotonic projection in stage 2.
proj_tol – Convergence tolerance for projection.
preserve_marginals – If True, re-apply column scaling after isotonic projection to maintain marginals (may re-introduce rank violations).
verbose – If True, enables debug logging.

Returns:

TwoStageResult with calibrated matrix and diagnostics.

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import calibrate_two_stage
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2], [0.2, 0.3, 0.5]])
>>> M = np.array([1.0, 1.2, 0.8])
>>> result = calibrate_two_stage(P, M)
>>> print(f"Converged: {result.converged}")
Converged: True

rank_preserving_calibration.classwise_ece(y: ndarray, probs: ndarray, n_bins: int = 15) → dict[source]

Compute class-wise Expected Calibration Error.

Parameters:

y – True class labels as integers of shape (N,).
probs – Probability matrix of shape (N, J).
n_bins – Number of equal-width probability bins per class (default 15).

Returns:

“ece”: Mean of the class-j ECEs
”mce”: Largest bin gap in any class
”per_class”: Per-class ECE, MCE and (count, mean prediction, frequency) bin tables

Return type:

Dictionary containing

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import classwise_ece
>>> y = np.array([0, 1, 2, 2])
>>> probs = np.array([[0.7, 0.2, 0.1], [0.2, 0.7, 0.1],
...                   [0.2, 0.2, 0.6], [0.6, 0.2, 0.2]])
>>> round(classwise_ece(y, probs, n_bins=2)["ece"], 4)
0.2167

rank_preserving_calibration.column_variance(Q: ndarray) → dict[str, Any][source]

Compute per-column variance of a probability matrix.

Useful for measuring how “spread out” the calibrated probabilities are within each class. Lower variance may indicate a “flatter” solution.

Parameters:

Q – Probability matrix of shape (N, J).

Returns:

“per_column”: List of variance values for each column
”mean”: Mean variance across columns
”min”: Minimum column variance
”max”: Maximum column variance

Return type:

Dictionary containing

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import column_variance
>>> Q = np.array([[0.6, 0.3, 0.1], [0.4, 0.4, 0.2]])
>>> var = column_variance(Q)
>>> print(f"Mean column variance: {var['mean']:.4f}")
Mean column variance: 0.0050

rank_preserving_calibration.compare_calibration_methods(P: ndarray, M: ndarray, results: dict[str, ndarray]) → dict[str, Any][source]

Compare multiple calibration results on flatness and constraint satisfaction.

Useful for deciding which calibration approach is best for a given problem by comparing solutions from different methods (QP, Dykstra, soft, two-stage).

Parameters:

P – Original probability matrix.
M – Target column sums.
results – Dictionary mapping method names to calibrated Q matrices.

Returns:

Dictionary with per-method metrics and rankings.

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import (
...     calibrate_dykstra, calibrate_soft, compare_calibration_methods
... )
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> M = np.array([0.9, 0.8, 0.3])
>>> Q_dykstra = calibrate_dykstra(P, M).Q
>>> Q_soft = calibrate_soft(P, M).Q
>>> comparison = compare_calibration_methods(P, M, {
...     "dykstra": Q_dykstra, "soft": Q_soft
... })

rank_preserving_calibration.distance_metrics(Q: ndarray, P: ndarray) → dict[str, float][source]

Compute distance metrics between original and calibrated probabilities.

Measures how much the calibration procedure has changed the probability estimates using various matrix norms and distance measures.

Parameters:

Q – Calibrated probability matrix of shape (N, J).
P – Original probability matrix of shape (N, J).

Returns:

“frobenius”: Frobenius norm ||Q - P||_F
”frobenius_sq”: Squared Frobenius norm ||Q - P||²_F
”mean_abs”: Mean absolute difference
”max_abs”: Maximum absolute difference

Return type:

Dictionary containing distance metrics

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import distance_metrics
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> Q = np.array([[0.65, 0.25, 0.1], [0.35, 0.45, 0.2]])
>>> distances = distance_metrics(Q, P)
>>> print(f"Frobenius distance: {distances['frobenius']:.4f}")
Frobenius distance: 0.1000

rank_preserving_calibration.feasibility_metrics(Q: ndarray, M: ndarray | None = None) → dict[str, Any][source]

Compute feasibility metrics for calibrated probability matrix.

Analyzes how well the probability matrix Q satisfies the row and column sum constraints required by rank-preserving calibration.

Parameters:

Q – Calibrated probability matrix of shape (N, J). Should have rows summing to 1 and non-negative entries.
M – Target column sums of shape (J,). If provided, analyzes column sum constraint satisfaction.

Returns:

“row”: Row constraint metrics – max/mean absolute errors and the minimum value
”col”: Column constraint metrics (if M provided) – various error norms

Return type:

Dictionary containing feasibility metrics

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import feasibility_metrics
>>> Q = np.array([[0.6, 0.3, 0.1], [0.4, 0.4, 0.2]])
>>> M = np.array([1.0, 0.7, 0.3])
>>> metrics = feasibility_metrics(Q, M)
>>> print(f"Max row error: {metrics['row']['max_abs_error']:.6f}")
Max row error: 0.000000
>>> print(f"Max col error: {metrics['col']['max_abs_error']:.6f}")
Max col error: 0.000000

rank_preserving_calibration.flatness_bound(P: ndarray, M: ndarray) → dict[str, Any][source]

Compute theoretical bound on expected solution flatness.

This function estimates how flat the calibrated solution is likely to be based on the distribution shift magnitude. When the shift is large relative to the variance in P, we expect Q to be flat.

The bound is based on the observation that calibration seeks Q closest to P subject to constraints. When constraints force Q far from P, the solution tends toward the uniform point Q_ij = M_j/N which satisfies all column constraints by construction.

Parameters:

P – Original probability matrix of shape (N, J).
M – Target column sums of shape (J,).

Returns:

“shift_magnitude”: Normalized L2 shift
”variance_p”: Total variance of P
”expected_variance_ratio”: Estimated Var(Q)/Var(P)
”flatness_risk”: Score from 0 to 1 indicating risk of flat solution
”recommendation”: String with practical advice

Return type:

Dictionary containing

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import flatness_bound
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> M = np.array([0.5, 0.8, 0.7])  # Large shift
>>> bound = flatness_bound(P, M)
>>> print(f"Flatness risk: {bound['flatness_risk']:.2f}")
Flatness risk: 0.32

rank_preserving_calibration.flatness_metrics(Q: ndarray, P: ndarray | None = None, M: ndarray | None = None) → dict[str, Any][source]

Measure how “flat” (uninformative) the calibrated solution is.

A perfectly flat solution has equal values within each column when sorted by original scores, meaning all samples receive the same probability for each class. This indicates loss of discriminative information.

Parameters:

Q – Calibrated probability matrix of shape (N, J).
P – Original probability matrix of shape (N, J). If provided, enables variance ratio computation.
M – Target column sums of shape (J,). If provided, enables comparison to the uniform solution Q_ij = M_j/N.

Returns:

“column_variance”: Per-column variance of Q values
”mean_column_variance”: Average variance across columns
”total_variance”: Total variance of all Q entries
”entropy_per_row”: Mean entropy of each row
”mean_max_prob”: Mean of maximum probability in each row
”variance_ratio”: Var(Q) / Var(P), lower means flatter. Only when P is provided.
”distance_to_uniform”: Frobenius distance from the uniform solution. Only when M is provided.

Return type:

Dictionary containing flatness metrics

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import flatness_metrics
>>> Q = np.array([[0.5, 0.3, 0.2], [0.4, 0.4, 0.2]])
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> metrics = flatness_metrics(Q, P)
>>> print(f"Variance ratio: {metrics['variance_ratio']:.3f}")
Variance ratio: 0.289

rank_preserving_calibration.informativeness_ratio(Q: ndarray, P: ndarray) → dict[str, Any][source]

Compare variance of Q vs P to measure information preservation.

Parameters:

Q – Calibrated probability matrix of shape (N, J).
P – Original probability matrix of shape (N, J).

Returns:

“total_ratio”: Var(Q) / Var(P)
”per_column_ratio”: Variance ratio for each column
”mean_column_ratio”: Mean of per-column ratios
”interpretation”: String describing the result

Return type:

Dictionary containing

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import informativeness_ratio
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> Q = np.array([[0.5, 0.35, 0.15], [0.5, 0.35, 0.15]])  # Flat
>>> ratio = informativeness_ratio(Q, P)
>>> print(f"Total ratio: {ratio['total_ratio']:.3f}")
Total ratio: 0.487

rank_preserving_calibration.is_feasible(Q: ndarray, M: ndarray, *, row_atol: float = 1e-10, col_atol: float = 1e-10) → bool[source]

Check that Q has unit row sums and column sums equal to M.

The tolerances are absolute. This function exists because that was not always true: every feasibility check in this module used to call np.allclose(x, target, atol=...) without passing rtol, and np.allclose keeps its default rtol=1e-5. The effective bound was atol + 1e-5 * |target|, so a stated atol=1e-12 on row sums behaved as 1e-5, and the column check – measured against M_j, which grows with N – became looser as the problem got larger. Convergence was therefore reported least reliably exactly where it mattered most.

Parameters:

Q – Candidate matrix of shape (N, J).
M – Target column sums of shape (J,).
row_atol – Absolute tolerance on |row_sum - 1|.
col_atol – Absolute tolerance on |col_sum - M_j|.

Returns:

True if every row sum and every column sum is within tolerance.

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import is_feasible
>>> Q = np.full((4, 2), 0.5)
>>> is_feasible(Q, np.array([2.0, 2.0]))
True

A row error of 1e-6 is a real violation, and is reported as one:

>>> Q_bad = Q.copy()
>>> Q_bad[0, 0] += 1e-6
>>> is_feasible(Q_bad, np.array([2.0, 2.0]))
False

rank_preserving_calibration.isotonic_metrics(Q: ndarray, P: ndarray) → dict[str, Any][source]

Analyze isotonic (rank-preserving) properties of calibrated probabilities.

Verifies that calibrated probabilities are nondecreasing when ordered by original model scores within each class. Computes violation statistics and flatness measures.

Parameters:

Q – Calibrated probability matrix of shape (N, J).
P – Original probability matrix of shape (N, J) used for score ordering.

Returns:

“max_rank_violation”: Maximum rank-order violation across all columns
”mass_weighted_violation”: Sum of violation magnitudes weighted by mass
”flat_fraction”: Fraction of adjacent pairs that are equal (flatness)
”per_class_max_violation”: List of max violations per class
”per_class_mass_violation”: List of mass-weighted violations per class

Return type:

Dictionary containing isotonic metrics

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import isotonic_metrics
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> Q = np.array([[0.65, 0.25, 0.1], [0.35, 0.45, 0.2]])
>>> metrics = isotonic_metrics(Q, P)
>>> print(f"Max violation: {metrics['max_rank_violation']:.6f}")
Max violation: 0.000000

rank_preserving_calibration.kl_divergence(Q: ndarray, R: ndarray, eps: float = 1e-15) → float[source]

Compute KL divergence KL(Q||R) = sum Q * log(Q/R).

The KL divergence measures how much Q differs from R, where R is the reference distribution. Uses convention 0 * log(0/x) = 0.

Parameters:

Q – Probability matrix of shape (N, J).
R – Reference probability matrix of shape (N, J).
eps – Small constant for numerical stability.

Returns:

KL divergence value (non-negative).

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import kl_divergence
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> Q = np.array([[0.6, 0.3, 0.1], [0.4, 0.4, 0.2]])
>>> kl = kl_divergence(Q, P)
>>> print(f"KL divergence: {kl:.4f}")
KL divergence: 0.0550

rank_preserving_calibration.marginal_shift_metrics(P: ndarray, M: ndarray) → dict[str, Any][source]

Quantify the distribution shift between P’s empirical marginals and M.

Larger shifts generally require more aggressive calibration adjustments, which can lead to flatter solutions. This function helps diagnose whether flatness is expected given the shift magnitude.

Parameters:

P – Original probability matrix of shape (N, J).
M – Target column sums of shape (J,).

Returns:

“empirical_marginals”: Column sums of P
”target_marginals”: M values
”marginal_diff”: Difference (M - empirical)
”l1_shift”: L1 norm of the difference
”l2_shift”: L2 norm of the difference
”max_shift”: Maximum absolute shift in any column
”relative_l2_shift”: L2 shift divided by sqrt(N)
”per_column_relative_shift”: Relative shift per column

Return type:

Dictionary containing shift metrics

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import marginal_shift_metrics
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> M = np.array([0.8, 0.8, 0.4])
>>> metrics = marginal_shift_metrics(P, M)
>>> print(f"L2 shift: {metrics['l2_shift']:.3f}")
L2 shift: 0.245

rank_preserving_calibration.project_near_isotonic_euclidean(v: ndarray, eps: float, sum_target: float | None = None, weights: ndarray | None = None) → ndarray[source]

Project v onto the set { z : z_{i+1} >= z_i - eps } in L2.

Parameters:

v – Input vector of shape (n,) to project.
eps – Slack parameter (>= 0). Allows z[i+1] >= z[i] - eps rather than the strict z[i+1] >= z[i].
sum_target – If provided, apply an exact uniform shift so the result sums to this value. Requires uniform weights.
weights – Positive weights of shape (n,) for a weighted L2 projection; cannot be combined with sum_target unless uniform.

Returns:

Projected vector of shape (n,) satisfying the near-isotonic constraint, and the sum constraint if one was requested.

Raises:

ValueError – If eps is negative, the weights are invalid, or sum_target is combined with non-uniform weights.

Notes

Idempotent: applying the projection twice gives the same result.

rank_preserving_calibration.project_near_kl_isotonic(v: ndarray, eps_slack: float, sum_target: float | None = None, weights: ndarray | None = None, eps: float = 1e-300) → ndarray[source]

Project v onto near-isotonic set with multiplicative slack under KL.

The constraint is: z[i+1] >= z[i] * (1 - eps_slack)

This is the KL analog of the Euclidean nearly-isotonic projection. With u[i] = z[i] / c**i and c = 1 - eps_slack the constraint becomes plain isotonicity in u, and KL(z||v) becomes a KL isotonic regression in u with weights c**i, solved by geometric-mean pooling.

Parameters:

v – Input vector (positive values).
eps_slack – Multiplicative slack parameter (0 to 1). eps_slack=0 is strict isotonic, eps_slack=0.1 allows 10% violation.
sum_target – If provided, the projection also sums to this. Exact by multiplicative rescaling, which holds only for uniform weights.
weights – Optional weights for weighted projection; cannot be combined with sum_target unless uniform.
eps – Numerical stability constant.

Returns:

Near-isotonic projection in KL sense.

Raises:

ValueError – If eps_slack is outside [0, 1), or sum_target is combined with non-uniform weights.

rank_preserving_calibration.prox_kl_near_isotonic(y: ndarray, lam: float, rho: float = 1.0, max_iters: int = 1000, tol: float = 1e-08, eps: float = 1e-300, return_info: bool = False) → ndarray | tuple[ndarray, dict][source]

Proximal operator for λ-penalty nearly-isotonic term under KL.

Solves:: min_z KL(z||y) + λ * sum_i (z_i - z_{i+1})_+

where (x)_+ = max(0, x) penalizes rank violations.

Uses ADMM in the dual space. This is more experimental than the epsilon-slack approach.

Parameters:

y – Input vector (positive values).
lam – Penalty weight for rank violations.
rho – ADMM penalty parameter.
max_iters – Maximum ADMM iterations.
tol – Convergence tolerance.
eps – Numerical stability constant.
return_info – If True, also return convergence info dict.

Returns:

z when return_info is False, otherwise (z, info_dict).

Raises:

ValueError – If lam is negative.

rank_preserving_calibration.prox_near_isotonic(y: ndarray, lam: float) → ndarray[source]

Compute the prox of the λ-penalty nearly-isotonic term exactly.

Solves:

minimize_z  0.5 * ||z - y||^2 + lam * sum_i (z_i - z_{i+1})_+

which charges every downward step and leaves upward ones free; lam = 0 returns y, and lam large enough returns the isotonic regression of y.

No maintained Python package solves this problem; the R package neariso (Hoefling, Tibshirani & Tibshirani) implements their path algorithm.

Parameters:

y – Input vector of shape (n,).
lam – Penalty weight on downward steps, non-negative.

Returns:

The minimiser z, of shape (n,).

Raises:

ValueError – If lam is negative.

rank_preserving_calibration.prox_near_isotonic_with_sum(y: ndarray, lam: float, sum_target: float) → ndarray[source]

Compute the prox under a sum constraint.

Solves:

minimize_z  0.5 * ||z - y||^2 + lam * sum_i (z_i - z_{i+1})_+
subject to  sum(z) = sum_target

The penalty depends only on differences, so it is translation-invariant, and the constrained prox is the unconstrained one shifted by a constant.

Parameters:

y – Input vector of shape (n,).
lam – Penalty weight on downward steps.
sum_target – Value the result must sum to.

Returns:

The constrained minimiser.

rank_preserving_calibration.reverse_kl_divergence(Q: ndarray, R: ndarray, eps: float = 1e-15) → float[source]

Compute reverse KL divergence KL(R||Q) = sum R * log(R/Q).

The reverse KL divergence has different properties than forward KL: - Forward KL(Q||R) is mean-seeking (Q spreads over modes of R) - Reverse KL(R||Q) is mode-seeking (Q focuses on main mode of R)

Parameters:

Q – Probability matrix of shape (N, J).
R – Reference probability matrix of shape (N, J).
eps – Small constant for numerical stability.

Returns:

Reverse KL divergence value (non-negative).

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import reverse_kl_divergence
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> Q = np.array([[0.6, 0.3, 0.1], [0.4, 0.4, 0.2]])
>>> rkl = reverse_kl_divergence(Q, P)

rank_preserving_calibration.sharpness_metrics(probs: ndarray) → dict[str, float][source]

Compute sharpness metrics for probability predictions.

Sharpness measures how confident (peaked) the predictions are, independent of whether they are correct. Sharper predictions are more informative.

Parameters:

probs – Probability matrix of shape (N, J).

Returns:

“mean_entropy”: Average entropy across predictions (lower = sharper)
”mean_max_prob”: Average of maximum probabilities (higher = sharper)
”mean_margin”: Average difference between top two probabilities

Return type:

Dictionary containing

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import sharpness_metrics
>>> probs = np.array([[0.9, 0.05, 0.05], [0.5, 0.3, 0.2]])
>>> sharp = sharpness_metrics(probs)
>>> print(f"Mean max prob: {sharp['mean_max_prob']:.3f}")
Mean max prob: 0.700

rank_preserving_calibration.symmetrized_kl(Q: ndarray, R: ndarray, eps: float = 1e-15) → float[source]

Compute symmetrized KL divergence (Jensen-Shannon related).

Symmetrized KL = 0.5 * (KL(Q||R) + KL(R||Q))

This is symmetric in Q and R, unlike the standard KL divergence.

Parameters:

Q – Probability matrix of shape (N, J).
R – Reference probability matrix of shape (N, J).
eps – Small constant for numerical stability.

Returns:

Symmetrized KL divergence value (non-negative).

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import symmetrized_kl
>>> P = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]])
>>> Q = np.array([[0.6, 0.3, 0.1], [0.4, 0.4, 0.2]])
>>> skl = symmetrized_kl(Q, P)

rank_preserving_calibration.tie_group_variance(Q: ndarray, P: ndarray) → dict[str, Any][source]: Within-equal-score group variance of Q; for ties=’group’ this should be ~0.

rank_preserving_calibration.top_label_ece(y: ndarray, probs: ndarray, n_bins: int = 15) → dict[source]

Compute Expected Calibration Error for top-label predictions.

Bins predictions by confidence (max probability) and measures the gap between confidence and accuracy in each bin. ECE is the weighted average of these gaps, and MCE is the maximum gap.

Parameters:

y – True class labels as integers of shape (N,).
probs – Probability matrix of shape (N, J).
n_bins – Number of confidence bins (default 15).

Returns:

“ece”: Expected Calibration Error (weighted average gap)
”mce”: Maximum Calibration Error (max gap across bins)
”bins”: List of (count, mean_confidence, accuracy) per bin

Return type:

Dictionary containing

Examples

>>> import numpy as np
>>> from rank_preserving_calibration import top_label_ece
>>> y = np.array([0, 1, 2])
>>> probs = np.array([[0.8, 0.1, 0.1], [0.2, 0.7, 0.1], [0.1, 0.2, 0.7]])
>>> ece_result = top_label_ece(y, probs)
>>> print(f"ECE: {ece_result['ece']:.4f}")
ECE: 0.2667