From e15f0fd082d61aeb90427fe7c297f9008235abbd Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 29 Oct 2026 14:50:02 +1030 Subject: [PATCH] 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) +