Promote dev to main for 0.6.3

CHANGELOG.md

@@ -4,6 +4,44 @@ All notable changes to this project are documented in this file.
 
 The format is based on *Keep a Changelog*, and this project follows *Semantic Versioning*.
 
+## [0.6.3] - 2026-03-19
+
+### Added
+
+- Public power-fit problem export via `pyvoro2.powerfit.build_power_fit_problem(...)`, including stable `PowerFitProblem`, `PowerFitBounds`, and `PowerFitPredictions` objects for external evaluation or solver experiments.
+- Public result packaging via `pyvoro2.powerfit.build_power_fit_result(...)`, so externally computed candidate weights can be turned back into a normal `PowerWeightFitResult` with standard residuals, reports, and algebraic diagnostics.
+- `PowerFitObjectiveBreakdown` and `PowerWeightFitResult.objective_breakdown`, including explicit hard-bound satisfaction and per-penalty summaries for debugging or research workflows.
+- A small advanced problem-API regression suite covering problem export, read-only arrays, result repackaging, and status-detail reporting.
+
+### Changed
+
+- Shared public power-fit dataclasses now live in `pyvoro2.powerfit.types`, while pure public weight/radius conversions live in `pyvoro2.powerfit.transforms`.
+- Mathematical formulas for prediction, algebraic diagnostics, hard-bound conversion, gauge canonicalization, and objective evaluation are now centralized in `pyvoro2.powerfit.problem` instead of being split across solver/report helpers.
+- The standalone native solver, active-set layer, and report layer now consume that shared problem/evaluation machinery instead of maintaining their own private copies of the same formulas.
+- `PairBisectorConstraints` and the new problem-definition objects now expose read-only NumPy arrays to match their frozen-dataclass semantics more honestly.
+- `PowerWeightFitResult.status` is now open-ended and pairs with `status_detail` so externally packaged results can preserve solver-specific termination text without expanding the native status vocabulary.
+
+### Fixed
+
+- Internal formula drift risks between low-level fitting, active-set post-processing, and report serialization were removed by routing them all through one public problem/evaluation layer.
+
+## [0.6.2] - 2026-03-19
+
+### Added
+
+- `PowerWeightFitResult.edge_diagnostics` with theorem-facing edge-space arrays and summaries: `alpha`, `beta`, `z_obs`, `z_fit`, algebraic residuals, and weighted/unweighted algebraic inconsistency metrics.
+- Fit reports and per-constraint record exporters now serialize the new algebraic edge diagnostics alongside the existing measurement-space predictions and residuals.
+- A medium-size robust-fit regression test that checks the native Huber/ADMM path on a deterministic sparse-outlier benchmark rather than only on tiny two-point cases.
+
+### Changed
+
+- The iterative `solver='admm'` backend now splits on the predicted measurement variable itself (`y = beta + alpha (w_i - w_j)`) instead of on raw weight differences, so the linear solve uses the scientifically natural `alpha^2` scaling of the graph Laplacian.
+- The ADMM backend now warm-starts from the quadratic analytic fit when that warm start is identifiable, improving convergence on robust Huber fits without expanding the public optimizer API.
+
+### Fixed
+
+- Native Huber fitting via `FitModel(mismatch=HuberLoss(...))` and `solver='admm'` no longer stalls or blows up on the medium-size sparse-outlier benchmark family used during the JCAM paper work; the corrected native path now matches the expected measurement-space residual quality and converges in a small number of iterations on the representative failure case.
+
 ## [0.6.1] - 2026-03-16
 
 ### Added

docs/guide/powerfit.md

@@ -123,7 +123,9 @@ The result contains:
 - fitted `weights` and shifted `radii`,
 - predicted separator locations in both fraction and position form,
 - residuals in the chosen measurement space,
-- solver/termination metadata,
+- `edge_diagnostics` with algebraic edge-space quantities such as `z_obs`, `z_fit`, and weighted inconsistency summaries,
+- `objective_breakdown` with mismatch / penalty / regularization totals for the packaged candidate weights,
+- solver/termination metadata including optional `status_detail`,
 - and explicit infeasibility reporting for contradictory hard constraints.
 
 For example, if hard interval or equality restrictions cannot all hold
@@ -141,7 +143,47 @@ Both low-level fits and active-set results also provide `to_records(...)` helper
 that turn per-constraint diagnostics into plain Python rows for downstream
 packages, table exporters, or custom reporting.
 
-For radii output, 0.6.1 makes the gauge choice explicit:
+
+### Measurement-space vs algebraic edge-space diagnostics
+
+`PowerWeightFitResult` now exposes two complementary diagnostic views.
+
+Measurement-space quantities live in the same space as the chosen separator
+targets:
+
+- `target`, `predicted`, `residuals`
+
+Algebraic edge-space quantities live in the implied difference model
+
+\[
+ y = \beta + \alpha (w_i - w_j),
+\]
+
+with
+
+\[
+ z_{\mathrm{obs}} = \frac{y_{\mathrm{target}} - \beta}{\alpha},
+ \qquad
+ z_{\mathrm{fit}} = w_i - w_j.
+\]
+
+The edge diagnostics expose `alpha`, `beta`, `z_obs`, `z_fit`, the algebraic
+residual `z_obs - z_fit`, and edge weights
+
+\[
+ \omega = \mathrm{confidence} \cdot \alpha^2.
+\]
+
+The exported `weighted_rmse` is defined explicitly as
+
+\[
+ \sqrt{\mathrm{mean}(\omega r^2)},
+\]
+
+not as `sqrt(sum(w r^2) / sum(w))`. That distinction matters when you compare
+results to other code that uses a normalized weighted RMSE convention.
+
+For radii output, 0.6.3 makes the gauge choice explicit:
 
 - by default, `weights_to_radii(...)` uses the minimal additive shift that makes
   all returned radii non-negative;
@@ -327,6 +369,45 @@ text = pv.dumps_report_json(solve_report, sort_keys=True)
 pv.write_report_json(solve_report, 'solve_report.json', sort_keys=True)
 ```
 
+## Native Huber fit on sparse outliers
+
+A short robust-fitting example looks like this:
+
+```python
+model = pv.FitModel(mismatch=pv.HuberLoss(delta=0.03))
+fit = pv.fit_power_weights(points, constraints, model=model, solver='admm')
+
+print(fit.status, fit.rms_residual)
+print(fit.edge_diagnostics.weighted_rmse)
+```
+
+This is still the native `pyvoro2` solver path. The robust part comes from the
+measurement-space Huber objective, while `edge_diagnostics` lets you inspect
+the theorem-facing algebraic residuals directly.
+
+## Advanced problem export / repackaging
+
+For research workflows or external solvers, `pyvoro2` now exposes the resolved
+inverse problem itself:
+
+```python
+resolved = pv.powerfit.resolve_pair_bisector_constraints(points, raw_constraints)
+problem = pv.powerfit.build_power_fit_problem(resolved, model=model)
+
+weights = some_external_solver(problem)
+result = pv.powerfit.build_power_fit_result(
+    problem,
+    weights,
+    solver='external',
+    status='external_failure',
+    status_detail='candidate iterate only',
+)
+```
+
+This keeps `fit_power_weights(...)` solver-owned while still giving downstream
+code a stable public export of the mathematics, prediction formulas, objective
+evaluation, and result packaging.
+
 ## Current scope
 
 The current implementation supports both **3D** domains through `pyvoro2` and
@@ -341,9 +422,9 @@ The main current restriction is geometric, not algebraic:
 - 2D currently supports `Box` and rectangular `RectangularCell`;
 - there is **no** planar oblique-periodic `PeriodicCell` yet.
 
-### Objective-model scope for 0.6.1
+### Objective-model scope for 0.6.3
 
-The 0.6.1 line still keeps the built-in objective family compact:
+The 0.6.3 line still keeps the built-in objective family compact:
 
 - mismatch terms: `SquaredLoss`, `HuberLoss`
 - hard feasibility: `Interval`, `FixedValue`
@@ -356,7 +437,7 @@ hard-feasibility checks, residual diagnostics, and solver behavior easy to
 reason about.
 
 Additional mismatch or penalty families should wait until downstream packages
-validate a concrete need for them. In particular, 0.6.1 does **not** try to
+validate a concrete need for them. In particular, 0.6.3 still does **not** try to
 freeze an open-ended callback API for arbitrary user-defined objectives.
 
 ## Worked example notebooks

src/pyvoro2/__about__.py

@@ -5,4 +5,4 @@
 """
 
 # Keep this as a simple assignment so scikit-build-core can extract it via regex.
-__version__ = '0.6.1'
+__version__ = '0.6.3'

src/pyvoro2/__init__.py

@@ -52,6 +52,7 @@
     ReciprocalBoundaryPenalty,
     L2Regularization,
     FitModel,
+    AlgebraicEdgeDiagnostics,
     ConstraintGraphDiagnostics,
     ConnectivityDiagnostics,
     ConnectivityDiagnosticsError,
@@ -114,6 +115,7 @@
     'ReciprocalBoundaryPenalty',
     'L2Regularization',
     'FitModel',
+    'AlgebraicEdgeDiagnostics',
     'ConstraintGraphDiagnostics',
     'ConnectivityDiagnostics',
     'ConnectivityDiagnosticsError',

src/pyvoro2/powerfit/__init__.py

@@ -22,6 +22,11 @@
     SelfConsistentPowerFitResult,
     solve_self_consistent_power_weights,
 )
+from .problem import (
+    PowerFitProblem,
+    build_power_fit_problem,
+    build_power_fit_result,
+)
 from .realize import (
     RealizedPairDiagnostics,
     UnaccountedRealizedPair,
@@ -35,16 +40,18 @@
     dumps_report_json,
     write_report_json,
 )
-from .solver import (
+from .solver import ConnectivityDiagnosticsError, fit_power_weights
+from .transforms import radii_to_weights, weights_to_radii
+from .types import (
+    AlgebraicEdgeDiagnostics,
     ConnectivityDiagnostics,
-    ConnectivityDiagnosticsError,
     ConstraintGraphDiagnostics,
     HardConstraintConflict,
     HardConstraintConflictTerm,
+    PowerFitBounds,
+    PowerFitObjectiveBreakdown,
+    PowerFitPredictions,
     PowerWeightFitResult,
-    fit_power_weights,
-    radii_to_weights,
-    weights_to_radii,
 )
 
 __all__ = [
@@ -59,12 +66,19 @@
     'ReciprocalBoundaryPenalty',
     'L2Regularization',
     'FitModel',
+    'AlgebraicEdgeDiagnostics',
     'ConstraintGraphDiagnostics',
     'ConnectivityDiagnostics',
     'ConnectivityDiagnosticsError',
     'HardConstraintConflictTerm',
     'HardConstraintConflict',
+    'PowerFitBounds',
+    'PowerFitPredictions',
+    'PowerFitObjectiveBreakdown',
+    'PowerFitProblem',
     'PowerWeightFitResult',
+    'build_power_fit_problem',
+    'build_power_fit_result',
     'RealizedPairDiagnostics',
     'UnaccountedRealizedPair',
     'UnaccountedRealizedPairError',

src/pyvoro2/powerfit/active.py

@@ -2,25 +2,26 @@
 
 from __future__ import annotations
 
-from dataclasses import dataclass, replace
+from dataclasses import dataclass
 from typing import Literal
 
 import numpy as np
 
 from .constraints import PairBisectorConstraints, resolve_pair_bisector_constraints
 from .model import FitModel
 from .realize import RealizedPairDiagnostics, match_realized_pairs
+from .problem import (
+    _build_active_set_connectivity_diagnostics,
+    build_power_fit_problem,
+    build_power_fit_result,
+)
 from .solver import (
     ConnectivityDiagnostics,
     PowerWeightFitResult,
     _apply_connectivity_policy,
-    _build_active_set_connectivity_diagnostics,
-    _connected_components,
-    _difference_identifying_mask,
-    _predict_measurements,
     fit_power_weights,
-    weights_to_radii,
 )
+from .transforms import weights_to_radii
 from ..diagnostics import TessellationDiagnostics as TessellationDiagnostics3D
 from ..domains import Box as Box3D, OrthorhombicCell, PeriodicCell
 from ..planar.diagnostics import TessellationDiagnostics as TessellationDiagnostics2D
@@ -323,6 +324,7 @@ def solve_self_consistent_power_weights(
             raise ValueError('active0 must have shape (m,)')
 
     warnings_list = list(resolved.warnings)
+    full_problem = build_power_fit_problem(resolved, model=model)
     add_streak = np.zeros(m, dtype=np.int64)
     drop_streak = np.zeros(m, dtype=np.int64)
     toggle_count = np.zeros(m, dtype=np.int64)
@@ -509,10 +511,10 @@ def solve_self_consistent_power_weights(
         n_added = int(np.count_nonzero((~active) & new_active))
         n_removed = int(np.count_nonzero(active & (~new_active)))
 
-        pred_fraction, pred_position, pred = _predict_measurements(
-            weights_eval,
-            resolved,
-        )
+        pred_all = full_problem.predict(weights_eval)
+        pred_fraction = np.asarray(pred_all.fraction, dtype=np.float64)
+        pred_position = np.asarray(pred_all.position, dtype=np.float64)
+        pred = np.asarray(pred_all.measurement, dtype=np.float64)
         target = (
             resolved.target_fraction
             if resolved.measurement == 'fraction'
@@ -616,6 +618,7 @@ def solve_self_consistent_power_weights(
             final_fit,
             active_constraints,
             final_weights,
+            model=model,
             r_min=r_min,
             weight_shift=weight_shift,
         )
@@ -631,10 +634,10 @@ def solve_self_consistent_power_weights(
             unaccounted_pair_check=unaccounted_pair_check,
         )
         warnings_list.extend(final_realized.warnings)
-        pred_fraction, pred_position, pred = _predict_measurements(
-            final_fit.weights,
-            resolved,
-        )
+        pred_all = full_problem.predict(final_fit.weights)
+        pred_fraction = np.asarray(pred_all.fraction, dtype=np.float64)
+        pred_position = np.asarray(pred_all.position, dtype=np.float64)
+        pred = np.asarray(pred_all.measurement, dtype=np.float64)
     else:
         final_realized = last_diag
         pred_fraction = np.full(m, np.nan, dtype=np.float64)
@@ -754,12 +757,11 @@ def _active_alignment_components(
     constraints: PairBisectorConstraints,
     model: FitModel,
 ) -> list[list[int]]:
-    effective_mask = _difference_identifying_mask(constraints, model)
-    return _connected_components(
-        constraints.n_points,
-        constraints.i[effective_mask],
-        constraints.j[effective_mask],
-    )
+    problem = build_power_fit_problem(constraints, model=model)
+    return [
+        list(component)
+        for component in problem.connectivity.effective_graph.connected_components
+    ]
 
 
 def _self_consistent_gauge_policy_description() -> str:
@@ -863,35 +865,24 @@ def _rebuild_fit_with_weights(
     constraints: PairBisectorConstraints,
     weights: np.ndarray,
     *,
+    model: FitModel,
     r_min: float,
     weight_shift: float | None,
 ) -> PowerWeightFitResult:
-    radii, shift = weights_to_radii(
+    problem = build_power_fit_problem(constraints, model=model)
+    return build_power_fit_result(
+        problem,
         weights,
+        solver=fit.solver,
+        status=fit.status,
+        status_detail=fit.status_detail,
+        converged=fit.converged,
+        n_iter=fit.n_iter,
+        warnings=fit.warnings,
+        canonicalize_gauge=False,
         r_min=r_min,
         weight_shift=weight_shift,
     )
-    pred_fraction, pred_position, pred = _predict_measurements(weights, constraints)
-    target = (
-        constraints.target_fraction
-        if constraints.measurement == 'fraction'
-        else constraints.target_position
-    )
-    residuals = pred - target
-    rms = float(np.sqrt(np.mean(residuals * residuals))) if residuals.size else 0.0
-    mx = float(np.max(np.abs(residuals))) if residuals.size else 0.0
-    return replace(
-        fit,
-        weights=np.asarray(weights, dtype=np.float64).copy(),
-        radii=radii,
-        weight_shift=shift,
-        predicted=pred,
-        predicted_fraction=pred_fraction,
-        predicted_position=pred_position,
-        residuals=residuals,
-        rms_residual=rms,
-        max_residual=mx,
-    )
 
 
 def _empty_realized_pair_diagnostics(