feat: warn when inequality constraints are supplied to a solver that ignores them

packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py

@@ -46,6 +46,8 @@ def _drop_nan_rows(A, B, idc, w, name):
 class DiscreteInterpolator(GeologicalInterpolator):
     """ """
 
+    _inequality_ignored_warned = False
+
     def __init__(self, support, data=None, c=None, up_to_date=False):
         """
         Base class for a discrete interpolator e.g. piecewise linear or finite difference which is
@@ -1008,6 +1010,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 +1237,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,
@@ -1375,6 +1385,15 @@ def solve_system(
             self._materialise_matrix_free_regularisation_blocks()
 
         Q, bounds = self._assemble_inequality_system(timing)
+        if solver_choice != "admm" and Q.shape[0] > 0 and not self._inequality_ignored_warned:
+            self._inequality_ignored_warned = True  # once per interpolator, not once per solve
+            # only ADMM (and the native solver, handled above) uses the inequality rows
+            solver_name = solver_choice if isinstance(solver_choice, str) else "custom"
+            logger.warning(
+                f"{Q.shape[0]} inequality constraint rows were supplied but solver "
+                f"'{solver_name}' does not use them, so they are ignored. Use "
+                "solver='admm' to enforce them."
+            )
         # Legacy compatibility: ignore deprecated backend switch.
         solver_kwargs.pop("backend", None)
 
@@ -1396,6 +1415,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)

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}

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:
@@ -95,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")
@@ -320,6 +324,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 +424,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 +442,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 +452,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 +471,32 @@ def _solve_linsys(system_matrix, rhs, kwargs):
             if dx_norm <= model_update_tol * model_scale:
                 break
 
+    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}); {consequence}. {advice}"
+        )
+
     if return_history:
         return model, history
     return model