Fit power weights from separator observations

pyvoro2 can solve a fixed-site inverse problem for power/Laguerre tessellations: fit power weights so that selected pairwise separators land at desired locations along the connector lines between sites. A normalized connector fraction may be any finite real value; restricting it to the segment between the two sites is a separate model constraint.

New code should begin with the concise fixed-observation surface in pyvoro2.inverse. Advanced objective models, realization checks, reports, and the experimental active-set outer loop live in pyvoro2.inverse.separator.

The API is geometry-first and domain-agnostic. The same high-level functions work with supported 3D domains and planar pyvoro2.planar domains. Downstream code decides:

• which site pairs are observed or proposed;
• which periodic image shift belongs to each observation;
• the target separator location;
• and the confidence of each observation.

pyvoro2 then provides the mathematical and geometric layers:

• resolve and validate separator observations;
• fit power weights under a configurable convex model;
• expose graph, connectivity, and hard-feasibility diagnostics;
• compute the resulting power tessellation;
• detect which requested pairs and periodic images are realized;
• and optionally run a realization-aware active-set outer loop.

The fixed-observation inverse fit and the geometric realization check answer different questions. For the API-independent derivation, see Inverse fitting from separator observations.

For namespace selection and lifecycle status, begin with Choosing an API. The glossary defines gauge, component offsets, representation shift, realized face, and active set. Users migrating from v0.6.3 should also read the v0.6.3–v0.8 migration guide.

The high-level resolver, observation container, fit result, fit entry point, and neutral transforms are stable. Advanced models, problem and operator views, report/realization helpers, and layered convenience views are provisional. Active-set refinement is experimental. The optional explicit sparse quadratic backend is provisional and supports large static sparse observation graphs only.

Canonical downstream integration

External IDs follow the same contract throughout separator resolution, fitting, active-set refinement, realization records, and reports: provide one unique non-negative integer per input site. Python integers and NumPy integer scalars are accepted. With index_mode='id', raw observation endpoints must be those exact integer IDs; floats, numeric strings, and booleans are rejected rather than converted.

The same exact policy applies to direct observation indices, periodic shift components, search and iteration counts, and active-set hysteresis counts. Public flags and masks accept only Python or NumPy Booleans. Confidence, targets, model parameters, regularization references, solver tolerances, r_min, and optional weight_shift must use finite real numeric values in their documented positive, non-negative, or signed ranges. Validation happens before casting or solver work.

Resolved observations own read-only C-contiguous copies of their numerical and mask data. Model scalars are canonical built-in floats, L2Regularization.reference is an owned read-only float64 copy even when its strength is zero, and FitModel.penalties is an owned tuple. Mutating a caller array or list after construction therefore does not change these values.

The repository-owned examples/chemvoro_workflow.py script is the canonical chemistry-neutral downstream example. It uses only current canonical imports and keeps application metadata outside pyvoro2 in an external-ID-keyed sidecar:

import numpy as np
import pyvoro2 as pv
import pyvoro2.inverse as inverse
import pyvoro2.inverse.separator as separator

points = np.array([[0.1, 0.5, 0.5], [0.9, 0.5, 0.5]])
site_ids = np.array([205, 101], dtype=int)
metadata_by_id = {
    205: {'label': 'left-site', 'source_row': 0},
    101: {'label': 'right-site', 'source_row': 1},
}
cell = pv.PeriodicCell(
    vectors=((1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0))
)

observations = inverse.resolve_separator_observations(
    points,
    [(205, 101, 0.5, (-1, 0, 0))],
    ids=site_ids,
    index_mode='id',
    domain=cell,
    image='given_only',
)
fit = inverse.fit_weights_from_separators(
    points,
    observations,
    connectivity_check='diagnose',
)
weights = fit.state.mathematical_weights

result = pv.compute(
    points,
    domain=cell,
    ids=site_ids,
    mode='power',
    weights=weights,
    include_empty=True,
    return_faces=True,
    return_face_shifts=True,
)
boundaries_by_input = result.require_boundaries()

realized = separator.match_realized_pairs(
    points,
    domain=cell,
    weights=weights,
    constraints=observations,
)
fit_rows = fit.to_records(observations, use_ids=True)
fit_report = fit.to_report(observations, use_ids=True)

result.ids, result.cell_measures, result.empty_mask, and the collections returned by require_boundaries() share input-site order. A downstream package can therefore build its own ID-labelled rows with an explicit None if empty else measure policy without reading raw backend order. Boundary records preserve external neighbor IDs and adjacent_shift; fit and realization records use external IDs when use_ids=True.

Inspect fit.identification before comparing fitted state across observation components. fit.state.global_representation_shift records only the common backend representation shift; it is not a fitted physical quantity. The complete executable example also demonstrates same-image and wrong-image realization reporting and JSON-friendly report export. The public deterministic paper-style ladder is in examples/paper_regressions.py; both scripts and their run instructions are described in examples/README.md.

Geometry of one pair

For a pair of sites i and j, choose one specific image q_j of site j. In a nonperiodic domain, q_j = p_j. Let

• d = ||q_j - p_i||,
• z = w_i - w_j,

where w are the fitted power weights.

Then the separator position along the connector is affine in z:

\[ t(z) = \frac{1}{2} + \frac{z}{2 d^2} \]

for normalized fraction, and

\[ s(z) = \frac{d}{2} + \frac{z}{2 d} \]

for absolute position measured from site i.

The connector line extends beyond both sites, so neither the geometry nor the measurement type inherently restricts t to [0, 1]. Add an explicit between-sites restriction when that is part of the observation model.

This is why pyvoro2 exposes the measurement type explicitly: a loss in fraction-space and a loss in position-space are different optimization problems.

Step 1: resolve separator observations once

import numpy as np
import pyvoro2 as pv
import pyvoro2.inverse as inverse
import pyvoro2.inverse.separator as separator

points = np.array([[0.0, 0.0, 0.0], [2.0, 0.0, 0.0]], dtype=float)
box = pv.Box(((-5, 5), (-5, 5), (-5, 5)))

observations = inverse.resolve_separator_observations(
    points,
    [(0, 1, 0.25)],
    measurement='fraction',
    domain=box,
)

Each raw tuple is (i, j, value[, shift]), where shift=(na, nb, nc) is the integer lattice image applied to site j.

An explicit shift names that image and is used unchanged, even if another image is nearer. If a periodic observation omits the shift and image='nearest', pyvoro2 certifies the minimum-image displacement for the exact binary64 values supplied. Exact ties use stable endpoint orientation, so reversing an ordered pair negates its shift and displacement. For separator inference, the stable keys are the resolved internal site indices in the fixed point ordering. External IDs remain metadata and do not participate in exact- tie geometry; relabeling them cannot change the selected image. This does not define a point-array permutation invariant or a public tie mode. The public image_search=1 default remains available as a performance hint: it seeds a bounded incumbent search but cannot change a successful shift, displacement, or distance. Certification that exceeds the private bounded resource contract fails explicitly; increasing image_search is not a correctness remedy.

The resolved SeparatorObservations object stores the validated pair indices, shifts, connector geometry, and targets in both fraction and position form. Resolver-created observations are also bound to the exact caller-order points, domain representation, dimension/count, and ID provenance used here. The public direct constructor remains useful for already-resolved row data; a valid directly constructed object is source-unbound rather than being assigned fabricated points or a domain.

When such an already resolved object is passed to a fitting function with points, those points establish or verify exact source points. Leaving domain=None makes no additional domain assertion and does not erase a domain already bound by the resolver. On an unbound object, the same call establishes an exact no-domain source; supplying a non-None domain establishes or verifies that exact representation. The observation object's owned IDs remain authoritative. Realization and active-set operations always verify the full source they use.

Every valid object has source-independent row IDs and an ordered observation-set fingerprint. Row identity includes the complete canonical row model but not warnings. Subsetting retains row IDs and input indices, duplicate rows remain distinct, and reordering changes the set fingerprint. A later verified source binding never changes those identities.

Step 2: define the fitting model

model = separator.FitModel(
    mismatch=separator.SquaredLoss(),
    feasible=separator.Interval(0.0, 1.0),
    penalties=(
        separator.ExponentialBoundaryPenalty(
            lower=0.0,
            upper=1.0,
            margin=0.05,
            strength=1.0,
            tau=0.01,
        ),
    ),
)

The model separates three ideas:

• mismatch=: how target-vs-predicted separator locations are scored,
• feasible=: hard admissible sets such as an interval or fixed value,
• penalties=: soft penalties such as outside-interval or near-boundary repulsion.

Built-in pieces currently include:

• SquaredLoss()
• HuberLoss(delta=...)
• Interval(lower, upper)
• FixedValue(value)
• SoftIntervalPenalty(lower, upper, strength=...)
• ExponentialBoundaryPenalty(...)
• ReciprocalBoundaryPenalty(...)
• L2Regularization(...)

For residual

\[ e_r(w)=\beta_r+\alpha_r(w_{i_r}-w_{j_r})-\mathrm{target}_r, \]

the squared loss is \(\frac12e_r^2\). Huber loss is

\[ \ell_\delta(e)= \begin{cases} \frac12e^2, & |e|\le\delta,\\ \delta\left(|e|-\frac12\delta\right), & |e|>\delta. \end{cases} \]

Row confidence multiplies only the mismatch loss. L2Regularization contributes

\[ \frac{\lambda}{2}\lVert w-w^{\mathrm{ref}}\rVert_2^2, \]

so its gradient is \(\lambda(w-w^{\mathrm{ref}})\) and its Hessian contribution is \(\lambda I\). A reference with strength=0 contributes nothing to the objective; the existing output-alignment policy may still use it to choose otherwise unidentified component offsets.

SoftIntervalPenalty(a, b, s) contributes

\[ s\left(\max(a-y,0)^2+\max(y-b,0)^2\right), \]

while ExponentialBoundaryPenalty(a, b, m, s, tau) contributes

\[ s\left[ \exp\left(\frac{a+m-y}{\tau}\right) +\exp\left(\frac{y-(b-m)}{\tau}\right) \right]. \]

Those existing strengths are not rescaled by a one-half factor. For one inward distance \(d\), ReciprocalBoundaryPenalty uses

\[ q(d)= \begin{cases} 0, & d\ge m,\\ s(1/d-1/m), & \epsilon<d<m,\\ s\left[(1/\epsilon-1/m)-(d-\epsilon)/\epsilon^2\right], & d\le\epsilon, \end{cases} \]

and contributes \(q(y-a)+q(b-y)\). The linear branch is used at \(d=\epsilon\) and the inactive branch at \(d=m\). Value and first derivative are continuous at \(\epsilon\); the derivative jumps at \(m\), so no unique derivative is claimed there. The parameters require upper > lower, margin > 0, 0 < epsilon < margin, 2 * margin <= upper - lower, and strength >= 0.

A scalar penalty with strength=0 is mathematically absent. Its value and derivatives are exactly zero without evaluating dangerous branch expressions; it does not force ADMM, hide the quadratic operator, change graph coupling, or change the fitted solution. An exactly zero named report component may remain.

For hard lower bound \(a\), prediction \(y\), and upper bound \(b\), classification uses

\[ v=\max(a-y,y-b,0), \qquad t=10^{-12}+64\epsilon_{64}\max(|a|,|y|,|b|), \]

where \(\epsilon_{64}\) is float64 machine epsilon. The row is satisfied exactly when \(v\le t\). This roundoff policy is also used by the hard-feasibility graph check and is separate from the ADMM stopping tolerances.

Step 3: fit power weights

fit = inverse.fit_weights_from_separators(
    points,
    observations,
    model=model,
)

The default is a certified direct solve with dense NumPy linear algebra. Select sparse SciPy linear algebra explicitly when the fixed observation graph is large and local:

python -m pip install "pyvoro2[sparse]"

sparse_fit = inverse.fit_weights_from_separators(
    points,
    observations,
    solver='direct',
    linear_backend='sparse',
)
print(sparse_fit.solver_termination.solver)          # direct
print(sparse_fit.solver_termination.linear_backend)  # sparse

solver='direct' accepts a purely quadratic model. Huber mismatch, hard constraints, and positive-strength scalar penalties require solver='admm'. Whenever a component solve is required, explicit ADMM runs the iterative method, including for a purely quadratic model. linear_backend='dense' never imports SciPy; linear_backend='sparse' explicitly requires it. Neither route changes backend at a site-count threshold.

All four quadratic method/backend combinations use one continuous-objective success contract. A floating objective that rounds to zero is not sufficient: the returned binary64 weights must either make every active source term exactly zero or pass the universal forward objective-gap certificate. Exact small-case helpers may propose or evaluate candidates, but coordinatewise rounding of an exact optimum is not treated as a separate discrete-optimality proof. A continuous optimum that cannot be represented accurately by the returned binary64 vector therefore produces status='numerical_failure' rather than a weaker meaning of optimal.

The linear backend changes only matrix storage and weight-system linear algebra. Both backends use the same observation rows, periodic image labels, effective graph, component anchors, final component-alignment policy, and quadratic success certificate. The method and backend are inspectable through fit.solver, fit.linear_backend, and the corresponding fields on fit.solver_termination. A degenerate fit that needs no component solve, such as an empty observation set or an all-zero-confidence model with only singleton components, reports solver='none', linear_backend=None, and n_iter=0. When ADMM has completed iterations but final quadratic certification fails, n_iter retains the completed iteration count in the structured numerical failure result.

ADMM scalar updates are themselves certified. Mismatch-only squared and Huber rows keep their vectorized proximal path. A row that needs a positive-strength scalar penalty uses a private bracketed solver with exact branch breakpoints, rigorous one-sided derivative enclosures, scaled binary64 exponential evaluation, and safeguarded Newton/ordered-float bisection. It returns a coordinate only with proved exact point signs or an adjacent numeric-float sign bracket. Adjacent endpoints are selected by a direct termwise objective difference. A bracket, evaluation, expansion, or iteration failure produces status='numerical_failure' with no stale weights or objective breakdown; status_detail identifies the original and component-local observation rows and retains scalar iteration, expansion, candidate, bracket, derivative, localization, and fallback evidence. A failure during ADMM attempt k normally records only the k - 1 iterations that completed.

The v0.8 migration is:

Removed call	Current call
solver='auto'	solver='direct', linear_backend='dense'
solver='analytic'	solver='direct', linear_backend='dense'
solver='sparse'	solver='direct', linear_backend='sparse'
solver='admm'	solver='admm', linear_backend='dense'

The removed ADMM keywords max_iter, rho, tol_abs, and tol_rel become admm_max_iter, admm_rho, admm_abs_tol, and admm_rel_tol. The active-set wrapper uses the same names with a fit_ prefix.

The result contains:

• fitted weights and shifted radii,
• predicted separator locations in both fraction and position form,
• residuals in the chosen measurement space,
• edge_diagnostics with quantities such as z_obs, z_fit, and weighted difference-space inconsistency summaries,
• objective_breakdown with mismatch, penalty, and regularization totals for the packaged candidate weights, plus hard_max_violation and the scale-aware hard_max_tolerance used for hard-bound classification,
• solver/termination metadata including optional status_detail,
• and explicit infeasibility reporting for contradictory hard constraints.

Read a fit through its scientific layers

The layers are deliberately one-directional rather than one monolithic result:

resolved observations + sites
            |
            v
fixed-observation fit
    |-- fitted state and identification
    |-- observation predictions and objective
    |-- graph and operator diagnostics
    `-- fixed-solver termination
            |
            v  (explicit forward realization request)
realized geometry and requested-image matching
            |
            v  (experimental outer loop only)
active-set path and outer termination

A small algebraic residual does not imply that the requested pair is a realized face. Realization and the active-set path are therefore separate objects and separate lifecycle layers.

SeparatorFitResult keeps all existing flat fields and adds lightweight layered views. The views reference existing arrays; they do not copy fitted or observation data.

state = fit.state
print(state.mathematical_weights)
print(state.backend_radii, state.global_representation_shift)

observation_fit = fit.observation_view(observations)
print(observation_fit.targets, observation_fit.confidence)
print(observation_fit.predictions, observation_fit.residuals)

identification = fit.identification
print(identification.effective_observation_components)
print(identification.relative_component_offsets_identified_by_data)
print(identification.component_offsets_selected_by_objective)
print(identification.component_alignment_policy)

termination = fit.solver_termination
print(
    termination.status,
    termination.solver,
    termination.linear_backend,
    termination.converged,
)

The identification view always reports global_geometric_gauge_identified_by_data == False: separator differences do not identify one common additive constant. The informative observation graph contains only positive-confidence separator rows. Its connected components are reported by effective_observation_components, and relative_component_offsets_identified_by_data is true exactly when this graph is connected. A positive L2 regularization term guarantees selection of otherwise free component offsets and is reported by component_offsets_selected_by_objective. Other supported scalar penalties are not classified as selecting offsets because they may have flat regions or zero strength.

identification.unconstrained_sites reports sites isolated in that informative graph. This can differ from the compatibility diagnostic fit.connectivity.unconstrained_points, which retains its candidate-graph meaning. A site mentioned only by zero-confidence rows is candidate-connected but observationally unconstrained. Hard restrictions and penalties apply independently of mismatch confidence and may constrain or bound offsets, but they are not separator-observation data and never add informative graph edges. An exact hard equality may fix an offset in a particular model; the current identification view deliberately does not attempt to summarize that separate constraint-identifiability question.

With connectivity_check='none', connectivity-derived identification values are None; accessing the view does not rebuild diagnostics that the caller disabled.

The state view uses global_representation_shift for the common shift used to form non-negative backend radii. This backend representation choice selects a representative within the global geometric gauge; it is distinct from independent component offsets and is not information recovered from separator observations. The compatibility flat field fit.weight_shift has exactly this meaning. Likewise, the compatibility field fit.connectivity.gauge_policy contains the same string as the canonical component_alignment_policy, despite its historical name.

The complete mapping is:

Scientific layer	Canonical access	Existing flat fields or objects
Fitted state and backend representation	fit.state	weights, radii, weight_shift
Identification and component alignment	fit.identification	connectivity and its effective components, offset flags, effective-graph isolated sites, and gauge_policy
Observation-space fit	fit.observation_view(observations)	measurement, target, predicted*, residuals, residual summaries, used_shifts; confidence comes from observations
Objective contributions	fit.objective	objective_breakdown
Algebraic diagnostics	fit.algebraic	edge_diagnostics, connectivity
Fixed-solver termination	fit.solver_termination	status, status_detail, solver, linear_backend, n_iter, converged, hard_feasible, conflict, warnings
Requested-image matching and realized geometry	realized.requested_image_matching, realized.geometry	all RealizedPairDiagnostics fields
Experimental outer-loop termination and path	result.outer_termination, result.path	active-set termination fields, active_mask, marginal_constraints, history, path_summary

The observation accessor applies one exact association rule before presenting observation-owned arrays beside fit predictions:

• two source-unbound objects match only when they have the exact same observation model;
• two source-bound objects match only when they have the exact same canonical source;
• a bound/unbound pair and two differently bound sources are rejected.

Length alone is never enough. Private observation origins and source bindings survive subsets, shallow and deep copies, same-version pickle round trips, dataclasses.replace(...), and copy.replace(...) where available. A source-inconsistent replacement raises. Reports always use the authoritative origin held by the result or diagnostic rather than borrowing provenance from an arbitrary supplied object.

For example, if hard interval or equality restrictions cannot all hold simultaneously, the fit returns:

• status == 'infeasible_hard_constraints'
• hard_feasible == False
• weights is None
• conflict with a compact contradiction witness
• conflicting_constraint_indices for the participating rows

instead of pretending the issue is merely slow convergence.

Both low-level fits and active-set results also provide to_records(...) helpers that turn per-constraint diagnostics into plain Python rows for downstream packages, table exporters, or custom reporting. Every observation-aligned row contains its stable row_id in addition to the existing fields.

Measurement-space and difference-space diagnostics

SeparatorFitResult exposes two complementary diagnostic views.

Measurement-space quantities live in the same space as the chosen separator targets:

• target, predicted, residuals

Difference-space quantities live in the implied weight-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 difference-space 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, the API makes the global representation shift explicit:

• by default, weights_to_radii(...) uses the minimal common additive shift that makes all returned radii non-negative;
• r_min= remains available as a compatibility-oriented convenience when a specific minimum radius is required;
• weight_shift= lets downstream code request one explicit common shift.

Both conversion directions require finite inputs and finite representable results. weights_to_radii(...) also requires a finite, non-negative r_min. If squaring a radius or applying the common shift would overflow, the transform raises ValueError; it does not clip or substitute a fallback value.

One common shift of every weight is the geometric gauge of the complete power diagram. Disconnected observation graphs introduce an additional and different ambiguity: each informative component can be shifted independently without changing its observed separator equations, but relative shifts between components can change the complete realized tessellation.

The current solver therefore chooses and reports a component-alignment policy:

• the numerical decomposition follows a model-coupling mask that includes positive-confidence rows and rows touched by hard restrictions or positive-strength penalties;
• without positive regularization or a supplied reference, disconnected model-coupling components are centered to mean zero;
• if a zero-strength regularization reference is supplied, disconnected model-coupling components are aligned to the reference mean;
• positive L2 regularization selects weights relative to its zero or supplied reference;
• connectivity_check='none'|'diagnose'|'warn'|'raise' controls whether the observationally unidentified component offsets are reported, warned about, or raised.

The model-coupling mask is exposed under the historical problem-field name offset_identifying_constraint_mask; despite that name, it is a solver decomposition detail and is not an identifiability claim. Hard restrictions may bound offsets, and penalties or numerical conventions may choose a returned representative without guaranteeing a unique optimum. None of those values is information identified by disconnected separator observations. Inspect realization results when cross-component competition matters.

The experimental active-set wrapper carries offsets from one iterate to the next only for true zero-L2 gauge components. Its final post-refit alignment is kept only when exact binary64-input checks prove that all within-component weight differences are unchanged. Positive L2 removes this gauge, so the final certified L2 solution is returned without alignment to the previous outer iterate.

Connectivity is computed on the graph of site unknowns, not on a graph of periodic images. A periodic shift changes the geometry of one observation row and of realized-boundary matching, but it does not create an additional fitted unknown. Interpenetrating periodic nets therefore remain disconnected unless an observation actually couples their site indices.

Inspect the observation graph and quadratic normal operator

Advanced workflows can inspect the fixed-observation mathematics directly from the public problem. The graph and operator views are provisional and are exported only from pyvoro2.inverse.separator.

import numpy as np
import pyvoro2.inverse as inverse
import pyvoro2.inverse.separator as separator

points = np.array(
    [[0.0, 0.0], [2.0, 0.0], [5.0, 0.0]],
    dtype=float,
)
observations = inverse.resolve_separator_observations(
    points,
    [(0, 1, 0.20), (1, 2, 0.70), (0, 2, 0.40)],
    confidence=[1.0, 0.5, 2.0],
)
model = separator.FitModel(
    regularization=separator.L2Regularization(
        strength=0.25,
        reference=np.array([1.0, -1.0, 2.0]),
    )
)
problem = separator.build_power_fit_problem(observations, model=model)

graph = problem.observation_graph
operator = problem.quadratic_operator
B = graph.incidence_dense()
L_obs = operator.observation_laplacian_dense()
A = operator.regularized_normal_matrix_dense()
b_obs = operator.observation_rhs
b = operator.regularized_normal_rhs

fit = inverse.fit_weights_from_separators(
    points,
    observations,
    model=model,
    connectivity_check='diagnose',
)
assert np.allclose(B.T @ fit.weights, problem.predict_difference(fit.weights))
assert np.allclose(
    graph.beta + graph.alpha * (B.T @ fit.weights),
    fit.predicted,
)
assert np.allclose(A @ fit.weights, b)

# Optional conversion; SciPy is also used by linear_backend='sparse'.
try:
    B_sparse = graph.incidence_sparse(format='csc')
    L_obs_sparse = operator.observation_laplacian_sparse(format='csr')
except ImportError:
    B_sparse = L_obs_sparse = None

For n sites and m observations, B.shape == (n, m). Column r has +1 at graph.site_i[r] and -1 at graph.site_j[r], hence B.T @ weights == weights[site_i] - weights[site_j]. Every observation remains a column: repeated measurements and observations for different periodic images of the same site pair are not deduplicated. graph.observation_indices maps those columns back to the resolved input indices, while graph.requested_shifts retains image identity.

The view names the two systems separately:

\[ \rho_r=c_r\alpha_r^2, \qquad q_r=c_r\alpha_r(y_r^{\mathrm{obs}}-\beta_r), \]

\[ L_{\mathrm{obs}}=B\operatorname{diag}(\rho)B^\mathsf{T}, \qquad b_{\mathrm{obs}}=Bq, \]

and, for L2 strength \(\lambda\) and reference \(w^{\mathrm{ref}}\),

\[ A=L_{\mathrm{obs}}+\lambda I, \qquad b=b_{\mathrm{obs}}+\lambda w^{\mathrm{ref}}. \]

There is no extra factor of two. Zero-confidence rows stay in B and in all row-facing arrays, but their informative mask is false and rho is zero, so they add nothing to \(L_{\mathrm{obs}}\) or \(b_{\mathrm{obs}}\) and do not connect informative components. The implementation constructs \(q\) directly with scale-safe products; z_obs remains diagnostic and is not required to be representable for a finite normal system.

Without positive L2 regularization, the observation-Laplacian nullity is the number of informative components, including isolated sites. One common null direction is global geometric gauge; additional component constants are unidentified offsets that can affect the complete realized diagram. Positive L2 regularization removes those null directions from A, but it does not make the separator observations themselves connected or observationally identify the offsets.

The repository benchmark harness benchmarks/benchmark_sparse_separator.py covers small dense-favorable and medium/large molecular-shaped locality graphs, including disconnected static components. It records dense/sparse assembly, direct-solve and complete-fit times, matrix storage, and numerical agreement. This scalability support is for large static geometries. It does not provide trajectory processing, MD frame reuse, prepared solvers across changing frames, parallel tessellation, GPU/distributed execution, or scalable all-pairs observation construction.

The fixed normal system is available only for SquaredLoss with no positive-strength scalar penalties. Zero-strength penalties are absent, so they do not hide this view. Huber mismatch and positive-strength scalar-penalty models still expose problem.observation_graph, but problem.quadratic_operator raises ValueError rather than claiming to represent their full objective. Hard interval or equality restrictions may coexist with the quadratic view; they remain separately visible through problem.bounds, and operator.normal_equations_characterize_fit is false because a constrained optimum need not solve the unconstrained normal equation.

Step 4: check geometric realization

A requested pairwise separator is not automatically a realized face in the full power tessellation. After fitting, you can ask which requested pairs became real neighbors.

realized = separator.match_realized_pairs(
    points,
    domain=box,
    weights=fit.state.mathematical_weights,
    constraints=observations,
    return_boundary_measure=True,
    return_tessellation_diagnostics=True,
    unaccounted_pair_check='warn',
)

weights= is the preferred realization input. radii= remains accepted as a backend-compatible representation route; supply exactly one. The common shift used to form radii is a geometric representation choice, not another fitted scientific variable.

This returns purely geometric diagnostics:

• whether each pair is realized at all,
• whether it is realized with the same requested periodic shift,
• whether only some other image is realized,
• whether one of the endpoint cells is empty,
• an optional boundary measure of the matched boundary (face area in 3D, edge length in 2D),
• any realized unordered site pairs that were absent from the candidate set, exposed through unaccounted_pairs,
• and optional tessellation-wide diagnostics.

Realization remains an explicit, separate computation. Read periodic-image matching and optional geometry independently:

matching = realized.requested_image_matching
print(matching.any_realization, matching.same_requested_shift)
print(matching.another_periodic_shift, matching.realized_shifts)

geometry = realized.geometry
print(geometry.endpoint_i_empty, geometry.endpoint_j_empty)
print(geometry.boundary_measure, geometry.tessellation_diagnostics)

The fit result never computes or owns a tessellation automatically.

Optional: refine the active set

For sparse or noisy candidate sets, the useful high-level workflow is often:

1. fit on a current active set;
2. run the actual power tessellation;
3. keep or re-add observations according to realized support;
4. repeat until active and realized sets agree.

This is a practical realization-aware outer algorithm. It is not part of the exact graph/Laplacian theory of the fixed-observation inner fit, and its termination status and path diagnostics should be inspected explicitly.

The explicitly experimental separator API provides this as:

result = separator.solve_self_consistent_power_weights(
    points,
    observations,
    domain=box,
    model=model,
    options=separator.ActiveSetOptions(
        add_after=1,
        drop_after=2,
        relax=0.5,
        max_iter=25,
        cycle_window=8,
    ),
    fit_solver='admm',
    fit_linear_backend='dense',
    return_history=True,
    return_boundary_measure=True,
    return_tessellation_diagnostics=True,
)

The solver is generic:

• it never invents candidate pairs,
• it never silently changes the user-supplied periodic image,
• it uses realized faces rather than any domain-specific contact logic,
• it supports hysteresis, under-relaxation, cycle detection, and marginal-pair reporting.

Reading the final diagnostics

solve_self_consistent_power_weights(...) returns both a final low-level fit and rich per-constraint diagnostics.

Useful fields include:

• result.constraints: the resolved pair set used throughout the solve,
• result.active_mask: final active-set membership,
• result.realized: optional realized-face matching diagnostics, including unaccounted_pairs when the final tessellation realizes candidate-absent pairs; it is None when the final fit has no usable weights,
• result.connectivity: final candidate-graph and active-graph connectivity diagnostics plus the component-alignment policy used for disconnected components,
• result.path_summary: compact optimization-path diagnostics that answer questions such as whether the fit-active graph was ever disconnected, whether active-component offsets were ever not identified by the pairwise data, and whether the tessellation ever realized pairs that were absent from the candidate set,
• result.history: optional per-iteration rows; each row distinguishes the fit-active mask (n_active_fit) from the post-toggle mask used for the next iteration (n_active), and also records fit-active component counts and the number of realized pairs absent from the candidate set on that iteration,
• result.diagnostics: optional per-constraint targets, predictions, residuals, endpoint-empty flags, boundary measure, toggle counts, and generic status labels,
• result.rms_residual_all / result.max_residual_all: summaries over all candidate constraints, or None when final weights are unavailable,
• result.tessellation_diagnostics: final tessellation-wide checks,
• result.marginal_constraints: indices of toggling / cycle / wrong-shift pairs.

The layered aliases make the inner/outer boundary explicit:

final_inner_fit = result.inner_fit
final_realization = result.final_realization
candidate_diagnostics = result.candidate_diagnostics
outer_termination = result.outer_termination
path = result.path

path.active_mask, path.marginal_constraint_indices, path.history, and path.summary share the existing active-set result data. The outer termination is experimental and does not change the exact fixed-observation meaning of final_inner_fit. In particular, outer and inner convergence are distinct:

assert result.converged == (result.termination == 'self_consistent')
final_state_available = result.final_state_available
unavailable_reason = result.final_state_unavailable_reason
final_inner_converged = result.final_refit_converged

A final optimal or max_iter fit with complete finite weights has an available state; all final geometry, predictions, residuals, and candidate records are recomputed from those exact weights. A weighted max_iter result remains non-converged at the inner layer. If a final fit has no usable weights, final_realization, candidate_diagnostics, the candidate residual summaries, and optional tessellation diagnostics are None. Historical path data, connectivity, warnings, cycle metadata, and path-derived marginal indices remain inspectable. No preceding realization is presented as the final one.

Transient path diagnostics are intentionally inspectable rather than noisy: final-state connectivity_check= / unaccounted_pair_check= policies still control warnings or exceptions, while result.path_summary and result.history expose optimization-path events without turning every transient component split into a default warning.

Status labels are intentionally generic, for example:

• stable_active
• stable_inactive
• toggled_active
• toggled_inactive
• realized_other_shift
• active_unrealized
• cycle_member

Exporting diagnostics as plain records

Downstream packages often want rows rather than structured NumPy-heavy result objects. The power-fitting package now exposes lightweight record exporters:

rows = result.to_records(use_ids=True)
fit_rows = result.fit.to_records(result.constraints, use_ids=True)
realized_rows = (
    None
    if result.realized is None
    else result.realized.to_records(result.constraints, use_ids=True)
)
if result.fit.conflict is not None:
    conflict_rows = result.fit.conflict.to_records(ids=result.constraints.ids)

These helpers keep the core API numerical while making it straightforward to feed results into custom logs, JSON encoders, or dataframe construction in a downstream package.

Full report bundles

When downstream code wants a single nested object rather than several row sets, use the report helpers or the corresponding result methods:

fit_report = fit.to_report(observations, use_ids=True)
realized_report = realized.to_report(observations, use_ids=True)
solve_report = result.to_report(use_ids=True)

The standalone helpers are also exported:

fit_report = separator.build_fit_report(fit, observations, use_ids=True)
solve_report = separator.build_active_set_report(result, use_ids=True)

These report bundles stay plain-Python and JSON-native. They are useful when a downstream package wants a complete diagnostic payload for logging, caching, or UI work without manually unpacking NumPy-heavy result objects. Every family retains its existing kind and adds the same versioned envelope:

assert fit_report['schema'] == {
    'name': 'pyvoro2.inverse.separator.report',
    'version': 1,
}
assert fit_report['producer'] == {
    'name': 'pyvoro2',
    'version': pv.__version__,
}

source = fit_report['source']
observation_set = fit_report['observation_set']

observation_set has exactly fingerprint, measurement, n_rows, and row_ids. source has exactly binding, fingerprint, dimension, n_points, points, domain, and ids. Resolver-backed reports use binding='bound' and retain exact caller-order points, the exact domain record, and ID provenance. A report produced by the valid public row-only chain uses:

{
    'binding': 'unbound',
    'fingerprint': None,
    'dimension': 2,
    'n_points': 3,
    'points': None,
    'domain': None,
    'ids': None,
}

An unbound null domain is different from a bound source whose domain is {'kind': 'none'}. Bound domains retain exactly one of none, planar_box, planar_rectangular_cell, spatial_box, spatial_orthorhombic_cell, and spatial_periodic_cell.

Report sections map to the same scientific layers:

Report section	Layer
fit weights, radii, weight_shift	state
fit connectivity	identification and graph context for algebraic diagnostics
fit constraints, fit_records, and residual summaries in summary	observations
fit objective_breakdown	objective
fit edge_diagnostics	algebraic diagnostics
fit summary, conflict, and warnings	fixed solver termination
realized records, unrealized	requested-image matching
realized unaccounted_pairs, tessellation_diagnostics, and optional values in records	realized geometry
active-set availability, fit, realized, diagnostics, summary, history, and path_summary	final-layer availability, final inner fit, optional final realization/candidate diagnostics, outer termination, and active-set path

The exact fit objective_breakdown fields are total, mismatch, penalties_total, penalty_terms, regularization, hard_constraints_satisfied, hard_max_violation, and hard_max_tolerance. The last two are respectively the maximum raw violation and the maximum rowwise tolerance actually used, or zero when no hard bounds exist.

To serialize them directly:

text = separator.dumps_report_json(solve_report, sort_keys=True)
separator.write_report_json(solve_report, 'solve_report.json', sort_keys=True)

dumps_report_json(...) rejects NaN and infinity. Fit, realized, and active reports round-trip exactly through JSON, including active no-weights failures. The active report's availability block records weights, realization, records, and reason. When unavailable, all three flags are false, reason is the final fit status, and weights-dependent report sections and summary values are JSON null. The nested fit still reports its own status and convergence while the active summary retains the outer stop reason.

Native Huber fit on sparse outliers

A short robust-fitting example looks like this:

model = separator.FitModel(mismatch=separator.HuberLoss(delta=0.03))
fit = inverse.fit_weights_from_separators(
    points,
    observations,
    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 mathematical difference-space residuals directly.

Advanced problem export and result packaging

For research workflows or external solvers, pyvoro2 now exposes the resolved inverse problem itself:

observations = inverse.resolve_separator_observations(points, raw_observations)
problem = separator.build_power_fit_problem(observations, model=model)

weights = some_external_solver(problem)
result = separator.build_power_fit_result(
    problem,
    weights,
    solver='external',
    status='external_failure',
    status_detail='candidate iterate only',
)

This keeps fit_weights_from_separators(...) solver-owned while giving downstream code a public export of the mathematics, prediction formulas, objective evaluation, and result packaging.

The result builder rejects status='optimal' or converged=True when any reported soft-objective component or total is NaN or infinite. Native solver paths turn such outcomes into status='numerical_failure'. Direct problem.evaluate_objective(...) may still return positive infinity for hard infeasibility or a genuine extended-real objective. Failure of the optional direct ADMM warm start falls back to the reference or zero initialization and does not by itself end the solve.

Current scope

The current implementation supports both 3D domains through pyvoro2 and 2D planar domains through pyvoro2.planar. The shared solver vocabulary is intentionally dimension-safe: fitting is phrased in terms of separator observations and generic boundary measure rather than chemistry-specific or 3D-only semantics.

The v0.7 compatibility package, broad top-level separator exports, and five historical core aliases are absent in v0.8. See the migration guide, architecture, and API lifecycle.

The main current restriction is geometric, not algebraic:

• 3D supports Box, OrthorhombicCell, and triclinic PeriodicCell;
• 2D currently supports Box and rectangular RectangularCell;
• there is no planar oblique-periodic PeriodicCell yet.

Current objective-model scope

The built-in objective family remains compact:

• mismatch terms: SquaredLoss, HuberLoss
• hard feasibility: Interval, FixedValue
• soft penalties: SoftIntervalPenalty, ExponentialBoundaryPenalty, ReciprocalBoundaryPenalty
• regularization: L2Regularization

That set is broad enough for the current generic inverse workflow while keeping 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, pyvoro2 does not try to freeze an open-ended callback API for arbitrary user-defined objectives.

Worked example notebooks

Four focused notebooks complement the guide:

• 04_powerfit presents the canonical external-ID, weight-first periodic workflow and then introduces advanced objective and active-set features.

• 06_powerfit_reports shows how to export low-level fits, realized-pair diagnostics, and self-consistent active-set results as rows or JSON-friendly reports.

• 07_powerfit_infeasibility shows how contradictory hard restrictions are reported through status, is_infeasible, conflict, and report bundles.
• 08_powerfit_active_path shows how to inspect transient active-set path diagnostics separately from the final-state report objects.

These examples are aimed at downstream packages that want to keep the solver API numerical while still producing human-readable logs, cached payloads, or UI views.