From e15f0fd082d61aeb90427fe7c297f9008235abbd Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 29 Oct 2026 14:50:02 +1030 Subject: [PATCH 1/4] test: pin the recipe and inequality-constraint problems (strict xfail) Recipe: to_recipe_json / from_recipe_json do not store a feature's build parameters (nelements, interpolatortype, regularisation), so a model built with anything but the defaults rebuilds with the defaults and gives a different field (differences of 12 to 60 on a field of range about 560). The default-parameter case passes; the three others are strict xfails. Inequality: ADMM enforces the bounds with the arguments the existing builder test uses, but with its default arguments returns essentially the unconstrained solution (4 of 24 bounds violated, worst 0.056), and with solver cg the inequality data are dropped without a warning. A control (the unconstrained solution does violate bounds) and the working ADMM configuration pass; the default-ADMM case and the cg-warning expectation are strict xfails. Co-Authored-By: Claude Sonnet 5.5 --- .../tests/test_inequality_enforcement.py | 122 ++++++++++++++++++ .../modelling/test_recipe_reproduces_field.py | 64 +++++++++ 2 files changed, 186 insertions(+) create mode 100644 packages/loop_interpolation/tests/test_inequality_enforcement.py create mode 100644 tests/unit/modelling/test_recipe_reproduces_field.py diff --git a/packages/loop_interpolation/tests/test_inequality_enforcement.py b/packages/loop_interpolation/tests/test_inequality_enforcement.py new file mode 100644 index 00000000..b0c4a018 --- /dev/null +++ b/packages/loop_interpolation/tests/test_inequality_enforcement.py @@ -0,0 +1,122 @@ +"""Inequality constraints must be enforced, or at least not silently ignored. + +``InterpolatorBuilder.add_inequality_constraints`` takes rows ``x, y, z, lower, upper`` and the +documented meaning is ``lower <= f(x, y, z) <= upper``. Inequality constraints are kept apart +from the equality system and are only used by the ADMM solver. + +The scenario is a sparse model on a gaussian dome whose true field is known: one borehole (three +contact points on the 0.3, 0.5 and 0.7 horizons) and two orientations. Thirty random points whose +true value is clearly above or below the 0.5 horizon become one-sided bounds on it: +``0.5 <= f`` above, ``f <= 0.5`` below. Without the bounds this data leaves 4 of the 24 bounded +points on the wrong side of 0.5. + +What was observed (not diagnosed): + +* ADMM with ``nmajor=10, admm_weight=0.01, maxiter=50`` (the arguments the existing builder test + uses) enforces every bound, and the error against the truth falls. +* ADMM with its default arguments, or with only some of those arguments, returns essentially the + unconstrained solution: 4 of 24 bounds still violated, worst by 0.056. +* With the ``cg`` solver the inequality data have no effect, and nothing is logged. + +The last two are strict xfails. When they start to pass they XPASS, which fails the run: remove +the markers then. +""" +import logging + +import numpy as np +import pytest +from loop_common.geometry import BoundingBox +from loop_interpolation import InterpolatorBuilder +from scipy.optimize import brentq + +TOLERANCE = 1e-3 +LEVEL = 0.5 + + +def truth(xyz): + r2 = (xyz[:, 0] - 0.5) ** 2 + (xyz[:, 1] - 0.5) ** 2 + return xyz[:, 2] - 0.2 * np.exp(-r2 / (2 * 0.2**2)) + + +def truth_gradient(xyz): + r2 = (xyz[:, 0] - 0.5) ** 2 + (xyz[:, 1] - 0.5) ** 2 + e = 0.2 * np.exp(-r2 / (2 * 0.2**2)) / 0.2**2 + return np.c_[e * (xyz[:, 0] - 0.5), e * (xyz[:, 1] - 0.5), np.ones(len(xyz))] + + +def sparse_data(): + contacts = [] + for level in (0.3, 0.5, 0.7): + z = brentq(lambda z, level=level: truth(np.array([[0.5, 0.5, z]]))[0] - level, 0.0, 1.0) + contacts.append([0.5, 0.5, z, level]) + orientation_xyz = np.array([[0.4, 0.4, 0.5], [0.6, 0.6, 0.5]]) + normals = np.c_[orientation_xyz, truth_gradient(orientation_xyz)] + + rng = np.random.default_rng(9) + candidates = rng.uniform(0.1, 0.9, (30, 3)) + rows = [] + for point, value in zip(candidates, truth(candidates)): + if value > LEVEL + 0.05: + rows.append([*point, LEVEL, 10.0]) # f >= 0.5 + elif value < LEVEL - 0.05: + rows.append([*point, -10.0, LEVEL]) # f <= 0.5 + return np.array(contacts), normals, np.array(rows) + + +def builder(with_inequality=True): + contacts, normals, inequality = sparse_data() + bounding_box = BoundingBox(np.array([0.0, 0.0, 0.0]), np.array([1.0, 1.0, 1.0])) + built = InterpolatorBuilder( + "FDI", bounding_box=bounding_box, nelements=4000, buffer=0.1 + ) + built.add_value_constraints(contacts).add_normal_constraints(normals) + if with_inequality: + built.add_inequality_constraints(inequality) + return built + + +def violations(interpolator): + """How far each bounded point is outside its bounds (0 where satisfied).""" + _, _, inequality = sparse_data() + value = interpolator.evaluate_value(inequality[:, :3]) + below = np.maximum(inequality[:, 3] - value, 0.0) + above = np.maximum(value - inequality[:, 4], 0.0) + return below + above + + +def test_the_unconstrained_solution_violates_some_bounds(): + """The scenario is only meaningful if the data alone get some bounded points wrong.""" + interpolator = builder(with_inequality=False).setup_interpolator().solve().build() + assert np.count_nonzero(violations(interpolator) > TOLERANCE) >= 3 + + +def test_admm_with_the_arguments_of_the_existing_builder_test_enforces_the_bounds(): + interpolator = ( + builder() + .setup_interpolator() + .use_solver("admm", nmajor=10, admm_weight=0.01, maxiter=50) + .solve() + .build() + ) + assert violations(interpolator).max() <= TOLERANCE + + +@pytest.mark.xfail( + strict=True, + reason="ADMM with its default arguments returns essentially the unconstrained solution " + "here (4 of 24 bounds violated, worst 0.056); cause not diagnosed", +) +def test_admm_with_default_arguments_enforces_the_bounds(): + interpolator = builder().setup_interpolator().use_solver("admm").solve().build() + assert violations(interpolator).max() <= TOLERANCE + + +@pytest.mark.xfail( + strict=True, + reason="expectation, not a documented contract: with solver 'cg' the inequality data are " + "dropped without a warning", +) +def test_cg_with_inequality_constraints_warns_that_they_are_not_used(caplog): + with caplog.at_level(logging.WARNING): + builder().setup_interpolator().use_solver("cg").solve().build() + assert any("inequality" in record.getMessage().lower() for record in caplog.records) diff --git a/tests/unit/modelling/test_recipe_reproduces_field.py b/tests/unit/modelling/test_recipe_reproduces_field.py new file mode 100644 index 00000000..27ecac30 --- /dev/null +++ b/tests/unit/modelling/test_recipe_reproduces_field.py @@ -0,0 +1,64 @@ +"""A model recipe should rebuild a model that gives the same scalar field. + +``GeologicalModel.to_recipe_json`` / ``from_recipe_json`` are documented as a way to describe and +rebuild a model. The existing recipe tests check the names of features, the data and the +stratigraphic column, but not the field, and the field depends on how each feature was built: +the number of elements, the interpolator type and the regularisation weight. + +Those build parameters are not stored: a feature entry in the recipe has only ``name``, ``kind``, +``faults`` and ``regions``. A model built with the defaults rebuilds exactly, because the rebuild +uses the defaults. A model built with anything else rebuilds with the defaults and gives a +different field (on the Claudius data, differences of 12 to 62 on a field of range about 560). + +The non-default cases are strict xfails. When the recipe stores the build parameters they will +XPASS, which fails the run: remove the markers then. +""" +import numpy as np +import pytest + +from LoopStructural import GeologicalModel +from LoopStructural.datasets import load_claudius + +RECIPE_DROPS_BUILD_PARAMETERS = pytest.mark.xfail( + strict=True, + reason="to_recipe_json stores only name, kind, faults and regions for a feature, not the " + "build parameters (nelements, interpolatortype, regularisation, ...), so the rebuilt " + "model uses the defaults", +) + +CASES = [ + pytest.param({}, id="defaults"), + pytest.param( + {"nelements": 2000}, id="nelements", marks=RECIPE_DROPS_BUILD_PARAMETERS + ), + pytest.param( + {"interpolatortype": "PLI"}, id="interpolatortype", marks=RECIPE_DROPS_BUILD_PARAMETERS + ), + pytest.param( + {"regularisation": 10.0}, id="regularisation", marks=RECIPE_DROPS_BUILD_PARAMETERS + ), +] + + +def build_model(**build_parameters): + data, bb = load_claudius() + model = GeologicalModel(bb[0, :], bb[1, :]) + model.data = data + model.stratigraphic_column.clear(basement=False) + model.create_and_add_foliation("strati", **build_parameters) + model.update() + return model, bb + + +@pytest.mark.parametrize("build_parameters", CASES) +def test_rebuilt_model_gives_the_same_field(build_parameters): + model, bb = build_model(**build_parameters) + + rebuilt = GeologicalModel.from_recipe_json(model.to_recipe_json()) + rebuilt.update() + + points = np.random.default_rng(0).uniform(bb[0], bb[1], (300, 3)) + original = model.evaluate_feature_value("strati", points) + again = rebuilt.evaluate_feature_value("strati", points) + np.testing.assert_allclose(again, original, atol=1e-6) + From 90e03fb89efc6f88d363d5bcc0ee84cf9f87bcd8 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 29 Oct 2026 15:34:43 +1030 Subject: [PATCH 2/4] fix: scale the inequality rows with the equality system in the ADMM solve solve_system column-scales the equality system (apply_scaling_matrix, on by default) and solves for y where c = S @ y, but passed the inequality matrix Q to ADMM unscaled. The bounds were therefore applied to y instead of to the field Q @ c. With scaling on, even a large penalty left bounds violated (admm_weight=10: 10 of 24). Pass Q @ S so Q @ y' = Q @ c. With this fixed, ADMM behaves like ordinary ADMM: admm_weight 1 and 10 now satisfy all 24 bounds (admm_weight=1 stops at iteration 90). The default penalty (0.01) is still too weak for the test scenario (primal residual 0.056 after 200 iterations against a tolerance of 0.0025) and the solver returns without a warning; that strict xfail stays, with its reason rewritten. The cg warning xfail is unchanged. Co-Authored-By: Claude Sonnet 5.5 --- .../_discrete_interpolator.py | 5 +++ .../tests/test_inequality_enforcement.py | 42 ++++++++++++++----- 2 files changed, 36 insertions(+), 11 deletions(-) diff --git a/packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py b/packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py index a51c1a3e..4818d0d7 100644 --- a/packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py +++ b/packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py @@ -1396,6 +1396,11 @@ def solve_system( if x0 is not None and scaling_matrix is not None: # the solver works on y where c = S @ y x0 = x0 / scaling_matrix.diagonal() + if scaling_matrix is not None and Q.shape[0] > 0: + # the solver works on y where c = S @ y, so the inequality rows must act on y + # too: Q @ c = (Q @ S) @ y. Without this the bounds are applied to y, not to the + # field. + Q = (Q @ scaling_matrix).tocsr() if callable(solver_choice): self.up_to_date = self._solve_with_callable(solver_choice, A, b, timing) diff --git a/packages/loop_interpolation/tests/test_inequality_enforcement.py b/packages/loop_interpolation/tests/test_inequality_enforcement.py index b0c4a018..34928370 100644 --- a/packages/loop_interpolation/tests/test_inequality_enforcement.py +++ b/packages/loop_interpolation/tests/test_inequality_enforcement.py @@ -10,16 +10,25 @@ ``0.5 <= f`` above, ``f <= 0.5`` below. Without the bounds this data leaves 4 of the 24 bounded points on the wrong side of 0.5. -What was observed (not diagnosed): - -* ADMM with ``nmajor=10, admm_weight=0.01, maxiter=50`` (the arguments the existing builder test - uses) enforces every bound, and the error against the truth falls. -* ADMM with its default arguments, or with only some of those arguments, returns essentially the - unconstrained solution: 4 of 24 bounds still violated, worst by 0.056. +Diagnosis: + +* ``solve_system`` column-scales the equality system (``apply_scaling_matrix``, on by default) and + solves for ``y`` where ``c = S @ y``. The inequality matrix ``Q`` was passed to ADMM unscaled, + so the bounds were applied to ``y`` instead of to the field ``Q @ c``. With scaling on, even a + large penalty left bounds violated (``admm_weight=10``: 10 of 24). Fixed by passing ``Q @ S``. +* With that fixed ADMM behaves like ordinary ADMM. Its default penalty ``admm_weight=0.01`` is too + weak for this problem: after the default 200 outer iterations the primal residual is 0.056 + against a tolerance of 0.0025 (still falling after 5000 iterations: 0.0145), so 4 of 24 bounds + are still violated. ``admm_weight=0.1`` stops at iteration 150, ``admm_weight=1`` at iteration + 90 with every bound satisfied. The solver reaches ``nmajor`` without meeting its tolerance and + returns without a warning. +* ``nmajor=10, maxiter=50`` (the arguments of the existing builder test) enforces the bounds + only by accident: it is an unconverged, truncated solve, and ``maxiter=100`` with the same + ``nmajor`` violates them again. * With the ``cg`` solver the inequality data have no effect, and nothing is logged. -The last two are strict xfails. When they start to pass they XPASS, which fails the run: remove -the markers then. +The default-penalty and ``cg`` cases are strict xfails. When they start to pass they XPASS, which +fails the run: remove the markers then. """ import logging @@ -103,14 +112,25 @@ def test_admm_with_the_arguments_of_the_existing_builder_test_enforces_the_bound @pytest.mark.xfail( strict=True, - reason="ADMM with its default arguments returns essentially the unconstrained solution " - "here (4 of 24 bounds violated, worst 0.056); cause not diagnosed", + reason="the default penalty admm_weight=0.01 is too weak here: after 200 iterations the primal " + "residual (0.056) is far above the tolerance (0.0025), 4 of 24 bounds are violated, and the " + "solver returns without a warning", ) -def test_admm_with_default_arguments_enforces_the_bounds(): +def test_admm_with_the_default_penalty_enforces_the_bounds(): interpolator = builder().setup_interpolator().use_solver("admm").solve().build() assert violations(interpolator).max() <= TOLERANCE +@pytest.mark.parametrize("apply_scaling_matrix", [True, False]) +@pytest.mark.parametrize("admm_weight", [1.0, 10.0]) +def test_admm_enforces_the_bounds_with_and_without_column_scaling(apply_scaling_matrix, admm_weight): + """The inequality rows must be scaled with the equality system (they were not).""" + built = builder() + built.interpolator.apply_scaling_matrix = apply_scaling_matrix + interpolator = built.setup_interpolator().use_solver("admm", admm_weight=admm_weight).solve().build() + assert violations(interpolator).max() <= TOLERANCE + + @pytest.mark.xfail( strict=True, reason="expectation, not a documented contract: with solver 'cg' the inequality data are " From 5f2bff99b12b3c2a4af6fc3d146e514d213b8faa Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 29 Oct 2026 16:01:11 +1030 Subject: [PATCH 3/4] feat: warn when ADMM uses all its iterations without converging admm_solve returned after nmajor iterations with no indication of whether it had converged, so an under-penalised or under-iterated solve left inequality constraints violated silently (with the default admm_weight=0.01 on the sparse dome, 4 of 24 bounds, worst 0.056). The loop now records how it ended. If it used every iteration without meeting the primal and dual tolerances it logs a warning with both residuals and their tolerances and suggests a larger admm_weight or nmajor. Stopping because the model stopped changing (model_update_tol, which the caller asks for) is not reported. Co-Authored-By: Claude Sonnet 5.5 --- .../loopsolver/admm_solver.py | 21 +++++ .../tests/test_inequality_enforcement.py | 80 ++++++++++++++++++- 2 files changed, 99 insertions(+), 2 deletions(-) diff --git a/packages/loop_interpolation/src/loop_interpolation/loopsolver/admm_solver.py b/packages/loop_interpolation/src/loop_interpolation/loopsolver/admm_solver.py index e033ea85..f1617a84 100644 --- a/packages/loop_interpolation/src/loop_interpolation/loopsolver/admm_solver.py +++ b/packages/loop_interpolation/src/loop_interpolation/loopsolver/admm_solver.py @@ -5,11 +5,14 @@ from dataclasses import dataclass import numpy as np +from loop_common.logging import get_logger from scipy.sparse import csr_matrix, diags, vstack from scipy.sparse.linalg import LinearOperator, cg, lsmr, lsqr from .admm_method import ADMM +logger = get_logger(__name__) + @dataclass class Config: @@ -320,6 +323,11 @@ def _solve_linsys(system_matrix, rhs, kwargs): ) history = [] + # how the major loop ended: it converged, it stopped because the model stopped changing (the + # caller asked for that with model_update_tol), or it used all of nmajor iterations + converged = False + stopped_on_model_update = False + last_residuals = None inner_x0 = None active_rows_prev = None use_cached_normal_eq = solver_name == "cg" and not active_set @@ -415,6 +423,8 @@ def _solve_linsys(system_matrix, rhs, kwargs): eps_pri = float(n_ie_sqrt * admm_abs_tol + admm_rel_tol * max(qx_norm, z_norm)) eps_dual = float(n_ie_sqrt * admm_abs_tol + admm_rel_tol * u_norm) + last_residuals = (primal_norm, eps_pri, dual_norm, eps_dual) + if return_history: history.append( { @@ -431,6 +441,7 @@ def _solve_linsys(system_matrix, rhs, kwargs): ) if (_i + 1) >= min_iterations and primal_norm <= eps_pri and dual_norm <= eps_dual: + converged = True break if model_update_tol > 0.0 and (_i + 1) >= min_iterations: @@ -440,6 +451,7 @@ def _solve_linsys(system_matrix, rhs, kwargs): history[-1]["stopped_on_model_update"] = True history[-1]["dx_norm"] = dx_norm history[-1]["model_norm"] = model_norm_before + stopped_on_model_update = True break if adaptive_rho_enabled and (_i + 1) >= min_iterations: @@ -458,6 +470,15 @@ def _solve_linsys(system_matrix, rhs, kwargs): if dx_norm <= model_update_tol * model_scale: break + if Q.shape[0] > 0 and not converged and not stopped_on_model_update and last_residuals: + primal_norm, eps_pri, dual_norm, eps_dual = last_residuals + logger.warning( + f"ADMM used all {nmajor} iterations without converging: primal residual " + f"{primal_norm:.3g} (tolerance {eps_pri:.3g}), dual residual {dual_norm:.3g} " + f"(tolerance {eps_dual:.3g}). The inequality constraints may still be violated. " + f"Try a larger admm_weight (the final value used was {rho:g}) or a larger nmajor." + ) + if return_history: return model, history return model diff --git a/packages/loop_interpolation/tests/test_inequality_enforcement.py b/packages/loop_interpolation/tests/test_inequality_enforcement.py index 34928370..5b94ba8e 100644 --- a/packages/loop_interpolation/tests/test_inequality_enforcement.py +++ b/packages/loop_interpolation/tests/test_inequality_enforcement.py @@ -113,8 +113,8 @@ def test_admm_with_the_arguments_of_the_existing_builder_test_enforces_the_bound @pytest.mark.xfail( strict=True, reason="the default penalty admm_weight=0.01 is too weak here: after 200 iterations the primal " - "residual (0.056) is far above the tolerance (0.0025), 4 of 24 bounds are violated, and the " - "solver returns without a warning", + "residual (0.056) is far above the tolerance (0.0025) and 4 of 24 bounds are violated (the " + "solver now warns about this, see the non-convergence tests below)", ) def test_admm_with_the_default_penalty_enforces_the_bounds(): interpolator = builder().setup_interpolator().use_solver("admm").solve().build() @@ -140,3 +140,79 @@ def test_cg_with_inequality_constraints_warns_that_they_are_not_used(caplog): with caplog.at_level(logging.WARNING): builder().setup_interpolator().use_solver("cg").solve().build() assert any("inequality" in record.getMessage().lower() for record in caplog.records) + + +# ------------------------------------------------------------------ non-convergence is reported +def not_converged_warnings(caplog): + return [r for r in caplog.records if "without converging" in r.getMessage()] + + +def test_the_default_penalty_that_does_not_converge_says_so(caplog): + with caplog.at_level(logging.WARNING): + builder().setup_interpolator().use_solver("admm").solve().build() + (warning,) = not_converged_warnings(caplog) + message = warning.getMessage() + assert "ADMM used all 200 iterations" in message + assert "primal residual" in message and "tolerance" in message + assert "admm_weight" in message + + +def test_a_penalty_that_converges_does_not_warn(caplog): + with caplog.at_level(logging.WARNING): + interpolator = ( + builder().setup_interpolator().use_solver("admm", admm_weight=1.0).solve().build() + ) + assert not not_converged_warnings(caplog) + assert violations(interpolator).max() <= TOLERANCE + + +def small_problem(): + """Data say the unknowns sum to 3; the first two are bounded to [0, 0.2].""" + from scipy import sparse + + n = 6 + a = sparse.vstack( + [sparse.identity(n, format="csr") * 0.2, sparse.csr_matrix(np.ones((1, n)))] + ).tocsr() + b = np.r_[np.zeros(n), 3.0] + q = sparse.csr_matrix(np.eye(n)[:2]) + bounds = np.array([[0.0, 0.2, 1.0], [0.0, 0.2, 1.0]]) + return a, b, q, bounds, np.zeros(n) + + +def solve_small(**kwargs): + from loop_interpolation.loopsolver import admm_solve + + a, b, q, bounds, x0 = small_problem() + return admm_solve(a, b, q, bounds, x0, linsys_solver="lsmr", **kwargs) + + +def test_solver_converged_in_time_does_not_warn(caplog): + with caplog.at_level(logging.WARNING): + x = solve_small(admm_weight=1.0, nmajor=200) + assert not not_converged_warnings(caplog) + np.testing.assert_allclose(x[:2], 0.2, atol=1e-3) + + +def test_solver_out_of_iterations_warns_and_reports_the_residuals(caplog): + with caplog.at_level(logging.WARNING): + x = solve_small(admm_weight=1e-4, nmajor=20) + (warning,) = not_converged_warnings(caplog) + assert "all 20 iterations" in warning.getMessage() + assert "the final value used was 0.0001" in warning.getMessage() + # and the bounds really are violated, so the warning is true + assert x[:2].max() > 0.2 + 0.1 + + +def test_solver_stopped_by_model_update_tolerance_does_not_warn(caplog): + """Stopping because the model stopped changing was asked for; it is not a failure.""" + with caplog.at_level(logging.WARNING): + solve_small(admm_weight=1e-4, nmajor=20, model_update_tol=10.0) + assert not not_converged_warnings(caplog) + + +def test_a_solver_with_fewer_iterations_than_the_minimum_warns(caplog): + """It cannot be declared converged before min_iterations, so it used all of its iterations.""" + with caplog.at_level(logging.WARNING): + solve_small(admm_weight=1.0, nmajor=3) + assert len(not_converged_warnings(caplog)) == 1 From 54e9d43f5be1ac2b7b5196471ecd158837f6b447 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 29 Oct 2026 18:06:20 +1030 Subject: [PATCH 4/4] fix: warn about ADMM non-convergence once per solve and describe the residuals accurately The constant-norm polish re-solves with ADMM every iteration, so the warning repeated N+1 times. The first solve stays quiet when a polish follows and only the last polish iteration reports. When only the dual residual missed its tolerance the message no longer says the inequality constraints may be violated. warn_on_nonconvergence=False silences it. Co-Authored-By: Claude Sonnet 5.5 --- .../_discrete_interpolator.py | 8 ++++ .../loop_interpolation/_solver_strategy.py | 1 + .../loopsolver/admm_solver.py | 24 ++++++++++-- .../tests/test_inequality_enforcement.py | 38 +++++++++++++++++++ 4 files changed, 68 insertions(+), 3 deletions(-) diff --git a/packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py b/packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py index 4818d0d7..699fcd17 100644 --- a/packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py +++ b/packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py @@ -1008,6 +1008,9 @@ def _run_constant_norm_polish( polish_solver_kwargs = dict(stable_solver_kwargs) polish_solver_kwargs["x0"] = lambda _support, x0=x0_seed: np.array(x0, copy=True) polish_solver_kwargs["constant_norm_iterations"] = 0 + # warn about non-convergence once per solve: only the final polish iteration reports + # (the earlier ones are intermediate, warm-started steps) + polish_solver_kwargs["warn_on_nonconvergence"] = i == int(iterations) - 1 solved = self.solve_system("admm", solver_kwargs=polish_solver_kwargs) if not solved: @@ -1232,6 +1235,11 @@ def _solve_with_admm( logger.warning("No inequality constraints, using lsmr") return self.solve_system("lsmr", solver_kwargs=solver_kwargs) + if constant_norm_iterations > 0 and constant_norm_weight > 0.0: + # the polish re-solves follow and the last one reports non-convergence for the + # whole solve, so this first solve stays quiet + solver_kwargs = {**solver_kwargs, "warn_on_nonconvergence": False} + try: c, history, ok = _solver_strategy.solve_with_admm( A, diff --git a/packages/loop_interpolation/src/loop_interpolation/_solver_strategy.py b/packages/loop_interpolation/src/loop_interpolation/_solver_strategy.py index a9f447f1..e34dfa78 100644 --- a/packages/loop_interpolation/src/loop_interpolation/_solver_strategy.py +++ b/packages/loop_interpolation/src/loop_interpolation/_solver_strategy.py @@ -183,6 +183,7 @@ def extract_admm_kwargs(solver_kwargs: dict, admm_solve) -> tuple[str, dict]: "inner_maxiter_schedule": solver_kwargs.pop("inner_maxiter_schedule", "linear"), "model_update_tol": solver_kwargs.pop("model_update_tol", 0.0), "return_history": solver_kwargs.pop("return_history", False), + "warn_on_nonconvergence": solver_kwargs.pop("warn_on_nonconvergence", True), } supported_optional = set(inspect.signature(admm_solve).parameters.keys()) admm_kwargs = {k: v for k, v in admm_kwargs.items() if k in supported_optional} diff --git a/packages/loop_interpolation/src/loop_interpolation/loopsolver/admm_solver.py b/packages/loop_interpolation/src/loop_interpolation/loopsolver/admm_solver.py index f1617a84..0ab48459 100644 --- a/packages/loop_interpolation/src/loop_interpolation/loopsolver/admm_solver.py +++ b/packages/loop_interpolation/src/loop_interpolation/loopsolver/admm_solver.py @@ -98,6 +98,7 @@ def admm_solve( model_update_tol: float = 0.0, min_iterations: int = 5, return_history: bool = False, + warn_on_nonconvergence: bool = True, ): if A.shape[1] != x0.shape[0]: raise ValueError("Number of columns in interpolation matrix does not match x0") @@ -470,13 +471,30 @@ def _solve_linsys(system_matrix, rhs, kwargs): if dx_norm <= model_update_tol * model_scale: break - if Q.shape[0] > 0 and not converged and not stopped_on_model_update and last_residuals: + if ( + warn_on_nonconvergence + and Q.shape[0] > 0 + and not converged + and not stopped_on_model_update + and last_residuals + ): primal_norm, eps_pri, dual_norm, eps_dual = last_residuals + if primal_norm > eps_pri: + consequence = "the inequality constraints may not be fully satisfied" + advice = ( + f"Try a larger admm_weight (the final value used was {rho:g}) or a larger nmajor." + ) + else: + # the bounds hold to tolerance; only the dual residual (optimality) missed + consequence = ( + "the inequality constraints are satisfied to tolerance but the solution may " + "not be fully converged" + ) + advice = f"Try a larger nmajor (the final admm_weight was {rho:g})." logger.warning( f"ADMM used all {nmajor} iterations without converging: primal residual " f"{primal_norm:.3g} (tolerance {eps_pri:.3g}), dual residual {dual_norm:.3g} " - f"(tolerance {eps_dual:.3g}). The inequality constraints may still be violated. " - f"Try a larger admm_weight (the final value used was {rho:g}) or a larger nmajor." + f"(tolerance {eps_dual:.3g}); {consequence}. {advice}" ) if return_history: diff --git a/packages/loop_interpolation/tests/test_inequality_enforcement.py b/packages/loop_interpolation/tests/test_inequality_enforcement.py index 5b94ba8e..1b6c1f80 100644 --- a/packages/loop_interpolation/tests/test_inequality_enforcement.py +++ b/packages/loop_interpolation/tests/test_inequality_enforcement.py @@ -155,6 +155,21 @@ def test_the_default_penalty_that_does_not_converge_says_so(caplog): assert "ADMM used all 200 iterations" in message assert "primal residual" in message and "tolerance" in message assert "admm_weight" in message + # the primal residual missed its tolerance here, so the bounds really may be violated + assert "may not be fully satisfied" in message + + +def test_the_constant_norm_polish_warns_once_not_once_per_iteration(caplog): + """The polish re-solves with ADMM each iteration; the warning is for the whole solve.""" + with caplog.at_level(logging.WARNING): + ( + builder() + .setup_interpolator() + .use_solver("admm", constant_norm_iterations=3, constant_norm_weight=0.1) + .solve() + .build() + ) + assert len(not_converged_warnings(caplog)) == 1 def test_a_penalty_that_converges_does_not_warn(caplog): @@ -204,6 +219,29 @@ def test_solver_out_of_iterations_warns_and_reports_the_residuals(caplog): assert x[:2].max() > 0.2 + 0.1 +def test_dual_only_miss_does_not_claim_the_bounds_are_violated(caplog): + """admm_weight=1 with 20 iterations meets the primal tolerance but not the dual one.""" + with caplog.at_level(logging.WARNING): + interpolator = ( + builder() + .setup_interpolator() + .use_solver("admm", admm_weight=1.0, nmajor=20) + .solve() + .build() + ) + (warning,) = not_converged_warnings(caplog) + message = warning.getMessage() + assert "satisfied to tolerance" in message + assert "may not be fully satisfied" not in message + assert violations(interpolator).max() <= TOLERANCE # and the message is true + + +def test_warning_can_be_silenced_per_solve(caplog): + with caplog.at_level(logging.WARNING): + solve_small(admm_weight=1e-4, nmajor=20, warn_on_nonconvergence=False) + assert not not_converged_warnings(caplog) + + def test_solver_stopped_by_model_update_tolerance_does_not_warn(caplog): """Stopping because the model stopped changing was asked for; it is not a failure.""" with caplog.at_level(logging.WARNING):