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) +