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 new file mode 100644 index 00000000..34928370 --- /dev/null +++ b/packages/loop_interpolation/tests/test_inequality_enforcement.py @@ -0,0 +1,142 @@ +"""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. + +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 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 + +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="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_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 " + "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) +