From e15f0fd082d61aeb90427fe7c297f9008235abbd Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 29 Oct 2026 14:50:02 +1030 Subject: [PATCH 1/5] 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 150c145b283cd3409759e0977794b9c89a52752f Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 29 Oct 2026 15:23:38 +1030 Subject: [PATCH 2/5] fix: store build parameters in model recipes, and refuse what they cannot rebuild A recipe rebuilt every feature with the default build parameters and as a foliation, so a model built with non-default nelements, interpolatortype or regularisation gave a different field (12 to 60 on a field of range ~560), a fault came back as a foliation, and a model with an unconformity could not be rebuilt at all. - create_and_add_feature records the parameters it was *requested* with on the feature (the element count realised by the interpolator would not reproduce the grid) - the recipe stores them per foliation (build.feature_type, build.parameters) and the reader replays them through create_and_add_feature - the stratigraphic column's group_mapping is stored - feature_recipe_kind tests subclasses first, so faults and unconformities are no longer labelled structural_frame / foliation - to_recipe_* raise LoopValueError for a feature it cannot rebuild yet (non-foliation, custom data=, a parameter JSON cannot hold, no creation record) instead of writing a wrong recipe - recipes written before this (no build entry) still load, with a warning; an older non-foliation entry is refused The three strict xfails from the recipe test now pass and their markers are removed. Co-Authored-By: Claude Sonnet 5.5 --- COMPAT.md | 1 + .../modelling/core/_model_serializer.py | 112 ++++++++++++- .../modelling/core/geological_model.py | 5 +- .../modelling/core/stratigraphic_column.py | 3 + tests/unit/modelling/test_recipe_limits.py | 148 ++++++++++++++++++ .../modelling/test_recipe_reproduces_field.py | 34 +--- 6 files changed, 268 insertions(+), 35 deletions(-) create mode 100644 tests/unit/modelling/test_recipe_limits.py diff --git a/COMPAT.md b/COMPAT.md index 11ee6262..5075c172 100644 --- a/COMPAT.md +++ b/COMPAT.md @@ -37,6 +37,7 @@ unchanged, but the new path is preferred going forward. |---|---|---| | `LoopStructural.modelling.features._feature_converters.add_fold_to_feature` (private module, imported directly by the QGIS plugin) | `GeologicalModel.add_fold_to_feature` | 2026-07-24 | | `LoopStructural.modelling.features._feature_converters.convert_feature_to_structural_frame` | `GeologicalModel.convert_feature_to_structural_frame` | 2026-07-24 | +| `GeologicalModel.to_recipe_json()` / `to_recipe_dict()` / `save_recipe()` silently wrote a lossy recipe for any model: only `name`, `kind`, `faults` for each feature (build parameters dropped, kinds mislabelled), and a model with a fault, unconformity or custom `data=` rebuilt wrongly or not at all | Each foliation entry now stores the parameters it was requested with (`build`), `group_mapping` is stored, and `from_recipe_*` replays them, so a model built with non-default `nelements`, `interpolatortype` or `regularisation` rebuilds to the same field. **A model containing a fault, unconformity or any other non-foliation feature, a feature built with its own `data=`, or a parameter JSON cannot hold, now raises `LoopValueError` from `to_recipe_*` instead of writing a recipe that rebuilds wrongly.** Older recipes (no `build`) still load, with a warning, but an older entry whose kind is not `foliation` is refused. Faults and unconformities are not yet supported | 2026-10-29 | | `GeologicalModel.create_and_add_fault(..., faults=[])` | `GeologicalModel.create_and_add_fault(..., faults=None)` (list built internally, same effective default) | 2026-07-27 | | `GeologicalModel.create_and_add_intrusion(..., intrusion_frame_parameters={}, geometric_scaling_parameters={})` | `GeologicalModel.create_and_add_intrusion(..., intrusion_frame_parameters=None, geometric_scaling_parameters=None)` (dicts built internally, same effective default) | 2026-07-27 | | `GeologicalModel.get_fault_surfaces(faults=[])` | `GeologicalModel.get_fault_surfaces(faults=None)` (list built internally, same effective default) | 2026-07-27 | diff --git a/LoopStructural/modelling/core/_model_serializer.py b/LoopStructural/modelling/core/_model_serializer.py index cbf4d112..70dbe545 100644 --- a/LoopStructural/modelling/core/_model_serializer.py +++ b/LoopStructural/modelling/core/_model_serializer.py @@ -13,6 +13,7 @@ import json import pathlib +import numpy as np import pandas as pd from ...geometry import BoundingBox @@ -24,19 +25,96 @@ logger = getLogger(__name__) +# Parameters of create_and_add_feature that are not stored with a feature's creation record: +# ``index`` is the position (the recipe keeps the order), ``faults`` are other features (the recipe +# keeps their names), and ``data`` is a data frame (see ``custom_data``). +_NOT_RECORDED = ("index", "faults", "data") + +# Feature types a recipe can rebuild from the parameters stored with it. +_SUPPORTED_FEATURE_TYPES = ("foliation",) + + +def _json_safe(value: object, where: str) -> object: + """Return ``value`` as something ``json`` can write, or raise ``LoopValueError``.""" + if value is None or isinstance(value, (bool, int, float, str)): + return value + if isinstance(value, np.generic): + return value.item() + if isinstance(value, np.ndarray): + return value.tolist() + if isinstance(value, (list, tuple)): + return [_json_safe(v, where) for v in value] + if isinstance(value, dict): + return {str(k): _json_safe(v, f"{where}[{k!r}]") for k, v in value.items()} + raise LoopValueError( + f"{where} is a {type(value).__name__}, which cannot be written to a recipe. " + "Only numbers, strings, booleans, None, lists and dictionaries of these can be stored." + ) + + class ModelSerializer: @staticmethod def feature_recipe_kind(feature): + # subclasses first: FaultSegment is a StructuralFrame and UnconformityFeature a + # GeologicalFeature + if isinstance(feature, FaultSegment): + return "fault" + if isinstance(feature, UnconformityFeature): + return "unconformity" if isinstance(feature, GeologicalFeature): return "foliation" if isinstance(feature, StructuralFrame): return "structural_frame" - if isinstance(feature, UnconformityFeature): - return "unconformity" - if isinstance(feature, FaultSegment): - return "fault" return feature.__class__.__name__.lower() + @staticmethod + def record_creation(feature: object, feature_type: str, params: dict) -> None: + """Remember how ``feature`` was created, for the recipe. + + Called by ``GeologicalModel.create_and_add_feature`` with the parameters that were + *requested* (for example ``nelements``, not the element count that was realised). + """ + feature._recipe_creation = { + "feature_type": feature_type, + "parameters": {k: v for k, v in params.items() if k not in _NOT_RECORDED}, + "custom_data": params.get("data") is not None, + } + + @staticmethod + def feature_build_entry(feature: object) -> dict: + """Return the ``build`` part of a feature's recipe entry. + + Raises + ------ + LoopValueError + if the feature cannot be rebuilt from a recipe yet + """ + kind = ModelSerializer.feature_recipe_kind(feature) + creation = getattr(feature, "_recipe_creation", None) + if creation is None: + raise LoopValueError( + f"Feature '{feature.name}' ({kind}) cannot be stored in a recipe: there is no " + "record of how it was created. Recipes store features created with " + "GeologicalModel.create_and_add_feature (or a create_and_add_* method)." + ) + feature_type = creation["feature_type"] + if kind != "foliation" or feature_type not in _SUPPORTED_FEATURE_TYPES: + raise LoopValueError( + f"Feature '{feature.name}' ({kind}, created as '{feature_type}') cannot be " + "stored in a recipe yet: recipes currently support foliations only." + ) + if creation["custom_data"]: + raise LoopValueError( + f"Feature '{feature.name}' was created with its own data=, which a recipe cannot " + "store: a recipe rebuilds a feature from the rows of the model data whose " + "feature_name matches. Put the data in model.data instead." + ) + parameters = { + key: _json_safe(value, f"parameter '{key}' of feature '{feature.name}'") + for key, value in creation["parameters"].items() + } + return {"feature_type": feature_type, "parameters": parameters} + @staticmethod def to_dict(model): result = {} @@ -68,6 +146,7 @@ def to_recipe_dict(model, data_reference=None): if getattr(fault, "name", None) ], "regions": [], + "build": ModelSerializer.feature_build_entry(feature), } recipe["model"]["features"].append(feature_entry) if data_reference is not None: @@ -132,10 +211,27 @@ def from_recipe_dict(cls, recipe): feature_name = feature_entry.get("name") if not isinstance(feature_name, str): raise TypeError("each feature entry must include a string name") - feature_data = model.data.loc[model.data["feature_name"] == feature_name].copy() - if feature_data.empty: - feature_data = None - feature = model.create_and_add_foliation(feature_name, data=feature_data) + build = feature_entry.get("build") + if build is None: + # a recipe written before build parameters were stored + kind = feature_entry.get("kind", "foliation") + if kind != "foliation": + raise LoopValueError( + f"Cannot rebuild feature '{feature_name}' of kind '{kind}' from this " + "recipe: recipes currently support foliations only." + ) + logger.warning( + f"The recipe has no build parameters for feature '{feature_name}' (it was " + "written by an older version); building it with the default parameters" + ) + build = {"feature_type": "foliation", "parameters": {}} + if not isinstance(build, dict) or not isinstance(build.get("parameters", {}), dict): + raise TypeError("a feature's build entry must be a dictionary with parameters") + feature = model.create_and_add_feature( + build.get("feature_type", "foliation"), + feature_name, + **build.get("parameters", {}), + ) if feature is None: raise ValueError(f"Could not recreate feature '{feature_name}' from recipe") feature_map[feature_name] = feature diff --git a/LoopStructural/modelling/core/geological_model.py b/LoopStructural/modelling/core/geological_model.py index cabaff0c..bf844e12 100644 --- a/LoopStructural/modelling/core/geological_model.py +++ b/LoopStructural/modelling/core/geological_model.py @@ -708,7 +708,10 @@ def create_and_add_feature(self, feature_type: str, name: str, **params): feature : BaseFeature the created geological feature, or None if it could not be built """ - return FeatureBuilderRegistry.create(feature_type, self, name, **params) + feature = FeatureBuilderRegistry.create(feature_type, self, name, **params) + if feature is not None: + ModelSerializer.record_creation(feature, feature_type, params) + return feature @public_api(tier="stable") def create_and_add_foliation( diff --git a/LoopStructural/modelling/core/stratigraphic_column.py b/LoopStructural/modelling/core/stratigraphic_column.py index a2d1971f..b1f83c04 100644 --- a/LoopStructural/modelling/core/stratigraphic_column.py +++ b/LoopStructural/modelling/core/stratigraphic_column.py @@ -641,6 +641,7 @@ def to_dict(self): """ return { "elements": [element.to_dict() for element in self.order], + "group_mapping": dict(self.group_mapping), } def update_from_dict(self, data): @@ -658,6 +659,7 @@ def update_from_dict(self, data): else: element = StratigraphicUnit.from_dict(element_data) self.add_element(element) + self.group_mapping = dict(data.get("group_mapping", {})) @classmethod def from_dict(cls, data): @@ -675,6 +677,7 @@ def from_dict(cls, data): else: element = StratigraphicUnit.from_dict(element_data) column.add_element(element) + column.group_mapping = dict(data.get("group_mapping", {})) return column def get_isovalues(self, where: str = 'bottom') -> dict[str, float]: diff --git a/tests/unit/modelling/test_recipe_limits.py b/tests/unit/modelling/test_recipe_limits.py new file mode 100644 index 00000000..741693cf --- /dev/null +++ b/tests/unit/modelling/test_recipe_limits.py @@ -0,0 +1,148 @@ +"""What a recipe can and cannot store, and that it says so instead of rebuilding something wrong.""" +import json + +import numpy as np +import pandas as pd +import pytest + +from LoopStructural import GeologicalModel +from LoopStructural.datasets import load_claudius +from LoopStructural.modelling.core._model_serializer import ModelSerializer +from LoopStructural.utils import LoopValueError + + +def claudius_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", **{"nelements": 1000, **build_parameters}) + return model, data, bb + + +def fault_model(): + data = pd.DataFrame( + [ + [5, 5, 5, 0, 0.70710678, 0.0, 0.70710678, 0, "fault"], + [5, 5, 5, 0, -0.70710678, 0.0, 0.70710678, 1, "fault"], + [8, 5, 5, 0, 0, 0, 1, np.nan, "strati"], + ], + columns=["X", "Y", "Z", "val", "nx", "ny", "nz", "coord", "feature_name"], + ) + model = GeologicalModel([0, 0, 0], [10, 10, 10]) + model.data = data + model.create_and_add_fault( + "fault", 1, nelements=1000, buffer=0.5, major_axis=10, minor_axis=3, intermediate_axis=10 + ) + model.create_and_add_foliation("strati", nelements=1000, faults=[model["fault"]]) + return model + + +# ------------------------------------------------------------------ what is stored +def test_the_requested_parameters_are_stored_not_the_realised_ones(): + model, _, _ = claudius_model(nelements=2000, regularisation=10.0) + entry = json.loads(model.to_recipe_json())["model"]["features"][0] + parameters = entry["build"]["parameters"] + assert entry["build"]["feature_type"] == "foliation" + assert parameters["nelements"] == 2000 + assert parameters["regularisation"] == 10.0 + assert parameters["interpolatortype"] == "FDI" + # the realised element count is different, and storing it would not reproduce the grid + model.update() + assert model["strati"].builder.build_arguments["nelements"] != 2000 + + +def test_data_and_other_features_are_not_stored_as_parameters(): + model, _, _ = claudius_model() + parameters = json.loads(model.to_recipe_json())["model"]["features"][0]["build"]["parameters"] + assert not {"data", "faults", "index"} & set(parameters) + + +def test_group_mapping_survives_a_round_trip(): + model, _, _ = claudius_model() + model.stratigraphic_column.add_unit("unit_0", thickness=60, id=0) + model.stratigraphic_column.group_mapping["Group_0"] = "strati" + rebuilt = GeologicalModel.from_recipe_json(model.to_recipe_json()) + assert rebuilt.stratigraphic_column.group_mapping == {"Group_0": "strati"} + + +def test_a_column_without_a_group_mapping_still_loads(): + """Recipes written before group_mapping was stored have no such key.""" + model, _, _ = claudius_model() + model.stratigraphic_column.add_unit("unit_0", thickness=60, id=0) + recipe = model.to_recipe_dict() + del recipe["model"]["stratigraphic_column"]["group_mapping"] + rebuilt = GeologicalModel.from_recipe_dict(recipe) + assert rebuilt.stratigraphic_column.group_mapping == {} + + +# ------------------------------------------------------------------ kind detection +def test_a_fault_is_recognised_as_a_fault(): + model = fault_model() + assert ModelSerializer.feature_recipe_kind(model["fault"]) == "fault" + assert ModelSerializer.feature_recipe_kind(model["strati"]) == "foliation" + + +def test_an_unconformity_is_recognised_as_an_unconformity(): + model, _, _ = claudius_model() + model.add_unconformity(model["strati"], 0) + kinds = {f.name: ModelSerializer.feature_recipe_kind(f) for f in model.features} + assert kinds["__strati_unconformity"] == "unconformity" + + +# ------------------------------------------------------------------ what is refused +def test_a_model_with_a_fault_is_refused_not_rebuilt_without_the_fault(): + with pytest.raises(LoopValueError, match=r"fault.*foliations only"): + fault_model().to_recipe_json() + + +def test_a_model_with_an_unconformity_is_refused(): + model, _, _ = claudius_model() + model.add_unconformity(model["strati"], 0) + with pytest.raises(LoopValueError, match=r"__strati_unconformity.*unconformity"): + model.to_recipe_json() + + +def test_a_feature_built_with_its_own_data_is_refused(): + model, data, _ = claudius_model() + model.create_and_add_foliation( + "strati_2", data=data[data["feature_name"] == "strati"].iloc[::2], nelements=1000 + ) + with pytest.raises(LoopValueError, match=r"strati_2.*its own data"): + model.to_recipe_json() + + +def test_a_parameter_json_cannot_hold_is_refused(): + model, _, _ = claudius_model(solver=lambda a, b: None) + with pytest.raises(LoopValueError, match=r"parameter 'solver' of feature 'strati'"): + model.to_recipe_json() + + +def test_a_feature_with_no_creation_record_is_refused(): + model, _, _ = claudius_model() + del model["strati"]._recipe_creation + with pytest.raises(LoopValueError, match="no record of how it was created"): + model.to_recipe_json() + + +# ------------------------------------------------------------------ older recipes +def test_an_older_recipe_without_build_parameters_loads_with_a_warning(caplog): + import logging + + model, _, _ = claudius_model() + recipe = model.to_recipe_dict() + del recipe["model"]["features"][0]["build"] + with caplog.at_level(logging.WARNING): + rebuilt = GeologicalModel.from_recipe_dict(recipe) + assert [f.name for f in rebuilt.features] == ["strati"] + assert any("no build parameters" in record.getMessage() for record in caplog.records) + + +def test_an_older_recipe_entry_that_is_not_a_foliation_is_refused(): + model, _, _ = claudius_model() + recipe = model.to_recipe_dict() + entry = recipe["model"]["features"][0] + del entry["build"] + entry["kind"] = "structural_frame" + with pytest.raises(LoopValueError, match="foliations only"): + GeologicalModel.from_recipe_dict(recipe) diff --git a/tests/unit/modelling/test_recipe_reproduces_field.py b/tests/unit/modelling/test_recipe_reproduces_field.py index 27ecac30..a6ac2bb7 100644 --- a/tests/unit/modelling/test_recipe_reproduces_field.py +++ b/tests/unit/modelling/test_recipe_reproduces_field.py @@ -1,17 +1,9 @@ """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. +The field depends on how each feature was built: the number of elements, the interpolator type +and the regularisation weight. The recipe stores the parameters each feature was *requested* +with (not the element count that was realised, which would not reproduce the grid), and the +rebuilt model repeats the request. """ import numpy as np import pytest @@ -19,23 +11,13 @@ 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"), + pytest.param({"interpolatortype": "PLI"}, id="interpolatortype"), + pytest.param({"regularisation": 10.0}, id="regularisation"), 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 + {"nelements": 2000, "regularisation": 10.0, "damp": False}, id="several-parameters" ), ] From 9784a6ab78196e96debf9f0b82d5a3a7a0ee7d4e Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 29 Oct 2026 15:31:08 +1030 Subject: [PATCH 3/5] feat: recipes rebuild faults and unconformities Stage 2 of the recipe fix. Foliations were already rebuilt from their stored parameters; now faults and unconformities are too. - a fault is stored with its creation parameters (displacement, axes, buffer, ...) and an explicit faults list by name: [] (no faults) is kept distinct from the default (all faults) - add_unconformity / add_onlap_unconformity record the feature, value and whether it is onlap on the object that lands in model.features (the inverse, for onlap); the reader replays them in order, which reproduces the regions they add to earlier features - foliation faults restored after the fact are only applied to foliations; faults and unconformities are built with theirs - an entry that refers to a feature that is not an earlier entry is refused - anything else (fold frames, folded foliations, domain faults, intrusions) is still refused when the recipe is written A fault model rebuilds to the identical strata and fault displacement (previously 0.21 and 0.49 different); erode and onlap unconformity models rebuild to identical fields and unit ids. Co-Authored-By: Claude Sonnet 5.5 --- COMPAT.md | 2 +- .../modelling/core/_model_relationships.py | 10 +- .../modelling/core/_model_serializer.py | 62 +++++-- .../test_recipe_faults_unconformities.py | 154 ++++++++++++++++++ tests/unit/modelling/test_recipe_limits.py | 23 +-- 5 files changed, 223 insertions(+), 28 deletions(-) create mode 100644 tests/unit/modelling/test_recipe_faults_unconformities.py diff --git a/COMPAT.md b/COMPAT.md index 5075c172..99e60e1e 100644 --- a/COMPAT.md +++ b/COMPAT.md @@ -37,7 +37,7 @@ unchanged, but the new path is preferred going forward. |---|---|---| | `LoopStructural.modelling.features._feature_converters.add_fold_to_feature` (private module, imported directly by the QGIS plugin) | `GeologicalModel.add_fold_to_feature` | 2026-07-24 | | `LoopStructural.modelling.features._feature_converters.convert_feature_to_structural_frame` | `GeologicalModel.convert_feature_to_structural_frame` | 2026-07-24 | -| `GeologicalModel.to_recipe_json()` / `to_recipe_dict()` / `save_recipe()` silently wrote a lossy recipe for any model: only `name`, `kind`, `faults` for each feature (build parameters dropped, kinds mislabelled), and a model with a fault, unconformity or custom `data=` rebuilt wrongly or not at all | Each foliation entry now stores the parameters it was requested with (`build`), `group_mapping` is stored, and `from_recipe_*` replays them, so a model built with non-default `nelements`, `interpolatortype` or `regularisation` rebuilds to the same field. **A model containing a fault, unconformity or any other non-foliation feature, a feature built with its own `data=`, or a parameter JSON cannot hold, now raises `LoopValueError` from `to_recipe_*` instead of writing a recipe that rebuilds wrongly.** Older recipes (no `build`) still load, with a warning, but an older entry whose kind is not `foliation` is refused. Faults and unconformities are not yet supported | 2026-10-29 | +| `GeologicalModel.to_recipe_json()` / `to_recipe_dict()` / `save_recipe()` silently wrote a lossy recipe for any model: only `name`, `kind`, `faults` for each feature (build parameters dropped, kinds mislabelled), and a model with a fault, unconformity or custom `data=` rebuilt wrongly or not at all | Each foliation entry now stores the parameters it was requested with (`build`), `group_mapping` is stored, and `from_recipe_*` replays them, so a model built with non-default `nelements`, `interpolatortype` or `regularisation` rebuilds to the same field. Foliations, faults and unconformities (erode and onlap) are rebuilt. **A model containing any other feature (fold frame, folded foliation, domain fault, intrusion, ...), a feature built with its own `data=`, or a parameter JSON cannot hold, now raises `LoopValueError` from `to_recipe_*` instead of writing a recipe that rebuilds wrongly.** Older recipes (no `build`) still load, with a warning, but an older entry whose kind is not `foliation` is refused | 2026-10-29 | | `GeologicalModel.create_and_add_fault(..., faults=[])` | `GeologicalModel.create_and_add_fault(..., faults=None)` (list built internally, same effective default) | 2026-07-27 | | `GeologicalModel.create_and_add_intrusion(..., intrusion_frame_parameters={}, geometric_scaling_parameters={})` | `GeologicalModel.create_and_add_intrusion(..., intrusion_frame_parameters=None, geometric_scaling_parameters=None)` (dicts built internally, same effective default) | 2026-07-27 | | `GeologicalModel.get_fault_surfaces(faults=[])` | `GeologicalModel.get_fault_surfaces(faults=None)` (list built internally, same effective default) | 2026-07-27 | diff --git a/LoopStructural/modelling/core/_model_relationships.py b/LoopStructural/modelling/core/_model_relationships.py index 5076864d..fc03db27 100644 --- a/LoopStructural/modelling/core/_model_relationships.py +++ b/LoopStructural/modelling/core/_model_relationships.py @@ -12,6 +12,7 @@ from ...modelling.features import FeatureType, UnconformityFeature from ...utils import getLogger +from ._model_serializer import ModelSerializer logger = getLogger(__name__) @@ -159,6 +160,9 @@ def add_unconformity(model, feature, value, index=None): else: f.add_region(uc_feature) # now add the unconformity to the feature list + ModelSerializer.record_creation( + uc_feature, "unconformity", {"feature": feature.name, "value": value, "onlap": False} + ) model._add_feature(uc_feature, index=index) return uc_feature @@ -195,6 +199,10 @@ def add_onlap_unconformity(model, feature, value, index=None): continue if f != feature: f.add_region(uc_feature) - model._add_feature(uc_feature.inverse(), index=index) + in_model = uc_feature.inverse() + ModelSerializer.record_creation( + in_model, "unconformity", {"feature": feature.name, "value": value, "onlap": True} + ) + model._add_feature(in_model, index=index) return uc_feature diff --git a/LoopStructural/modelling/core/_model_serializer.py b/LoopStructural/modelling/core/_model_serializer.py index 70dbe545..f35b39a6 100644 --- a/LoopStructural/modelling/core/_model_serializer.py +++ b/LoopStructural/modelling/core/_model_serializer.py @@ -25,13 +25,18 @@ logger = getLogger(__name__) -# Parameters of create_and_add_feature that are not stored with a feature's creation record: -# ``index`` is the position (the recipe keeps the order), ``faults`` are other features (the recipe -# keeps their names), and ``data`` is a data frame (see ``custom_data``). +# Parameters of create_and_add_feature that are not stored as they are: ``index`` is the position +# (the recipe keeps the order), ``data`` is a data frame (see ``custom_data``), and ``faults`` are +# other features, stored by name. _NOT_RECORDED = ("index", "faults", "data") -# Feature types a recipe can rebuild from the parameters stored with it. -_SUPPORTED_FEATURE_TYPES = ("foliation",) +# Feature types a recipe can rebuild from the parameters stored with it, and the kind that +# ``feature_recipe_kind`` must report for a feature created as that type. +_SUPPORTED_FEATURE_TYPES = { + "foliation": "foliation", + "fault": "fault", + "unconformity": "unconformity", +} def _json_safe(value: object, where: str) -> object: @@ -74,9 +79,14 @@ def record_creation(feature: object, feature_type: str, params: dict) -> None: Called by ``GeologicalModel.create_and_add_feature`` with the parameters that were *requested* (for example ``nelements``, not the element count that was realised). """ + parameters = {k: v for k, v in params.items() if k not in _NOT_RECORDED} + faults = params.get("faults") + if faults is not None: + # [] means "no faults", which is different from None ("all faults in the model") + parameters["faults"] = [f if isinstance(f, str) else f.name for f in faults] feature._recipe_creation = { "feature_type": feature_type, - "parameters": {k: v for k, v in params.items() if k not in _NOT_RECORDED}, + "parameters": parameters, "custom_data": params.get("data") is not None, } @@ -98,10 +108,10 @@ def feature_build_entry(feature: object) -> dict: "GeologicalModel.create_and_add_feature (or a create_and_add_* method)." ) feature_type = creation["feature_type"] - if kind != "foliation" or feature_type not in _SUPPORTED_FEATURE_TYPES: + if _SUPPORTED_FEATURE_TYPES.get(feature_type) != kind: raise LoopValueError( f"Feature '{feature.name}' ({kind}, created as '{feature_type}') cannot be " - "stored in a recipe yet: recipes currently support foliations only." + f"stored in a recipe yet: recipes support {sorted(_SUPPORTED_FEATURE_TYPES)}." ) if creation["custom_data"]: raise LoopValueError( @@ -218,7 +228,8 @@ def from_recipe_dict(cls, recipe): if kind != "foliation": raise LoopValueError( f"Cannot rebuild feature '{feature_name}' of kind '{kind}' from this " - "recipe: recipes currently support foliations only." + "recipe: it has no build parameters, and only foliations can be rebuilt " + "without them." ) logger.warning( f"The recipe has no build parameters for feature '{feature_name}' (it was " @@ -227,11 +238,30 @@ def from_recipe_dict(cls, recipe): build = {"feature_type": "foliation", "parameters": {}} if not isinstance(build, dict) or not isinstance(build.get("parameters", {}), dict): raise TypeError("a feature's build entry must be a dictionary with parameters") - feature = model.create_and_add_feature( - build.get("feature_type", "foliation"), - feature_name, - **build.get("parameters", {}), - ) + feature_type = build.get("feature_type", "foliation") + parameters = dict(build.get("parameters", {})) + if "faults" in parameters: + fault_names = parameters["faults"] + missing = [n for n in fault_names if n not in feature_map] + if missing: + raise LoopValueError( + f"Feature '{feature_name}' refers to {missing}, which are not earlier " + "entries of the recipe" + ) + parameters["faults"] = [feature_map[n] for n in fault_names] + if feature_type == "unconformity": + base = feature_map.get(parameters.get("feature")) + if base is None: + raise LoopValueError( + f"Unconformity '{feature_name}' refers to '{parameters.get('feature')}', " + "which is not an earlier entry of the recipe" + ) + add = ( + model.add_onlap_unconformity if parameters.get("onlap") else model.add_unconformity + ) + feature = add(base, parameters["value"]) + else: + feature = model.create_and_add_feature(feature_type, feature_name, **parameters) if feature is None: raise ValueError(f"Could not recreate feature '{feature_name}' from recipe") feature_map[feature_name] = feature @@ -241,7 +271,9 @@ def from_recipe_dict(cls, recipe): fault_names = feature_entry.get("faults", []) if not isinstance(fault_names, list): raise TypeError("faults must be a list") - if fault_names: + # restores faults assigned to a foliation after it was built; faults, unconformities + # and features given an explicit faults list were built with them + if fault_names and feature_entry.get("kind", "foliation") == "foliation": feature = feature_map[feature_name] feature.faults = [feature_map[name] for name in fault_names if name in feature_map] diff --git a/tests/unit/modelling/test_recipe_faults_unconformities.py b/tests/unit/modelling/test_recipe_faults_unconformities.py new file mode 100644 index 00000000..2d33fb1d --- /dev/null +++ b/tests/unit/modelling/test_recipe_faults_unconformities.py @@ -0,0 +1,154 @@ +"""Recipes rebuild faults and unconformities, not just foliations.""" +import json + +import numpy as np +import pandas as pd +import pytest + +from LoopStructural import GeologicalModel +from LoopStructural.datasets import load_claudius +from LoopStructural.modelling.core._model_serializer import ModelSerializer +from LoopStructural.utils import LoopValueError + + +def fault_data(): + return pd.DataFrame( + [ + [5, 5, 5, 0, 0.70710678, 0.0, 0.70710678, 0, "fault"], + [5, 5, 5, 0, -0.70710678, 0.0, 0.70710678, 1, "fault"], + [8, 5, 5, 0, 0, 0, 1, np.nan, "strati"], + ], + columns=["X", "Y", "Z", "val", "nx", "ny", "nz", "coord", "feature_name"], + ) + + +def fault_model(**foliation_parameters): + model = GeologicalModel([0, 0, 0], [10, 10, 10]) + model.data = fault_data() + model.create_and_add_fault( + "fault", 1, nelements=1000, buffer=0.5, major_axis=10, minor_axis=3, intermediate_axis=10 + ) + model.create_and_add_foliation("strati", nelements=1000, **foliation_parameters) + model.update() + return model + + +def two_group_model(onlap): + data, bb = load_claudius() + data = data.reset_index() + data.loc[:, "val"] *= -1 + data.loc[:, ["nx", "ny", "nz"]] *= -1 + data.loc[792, "feature_name"] = "strati2" + data.loc[792, ["nx", "ny", "nz"]] = [0, 0, 1] + data.loc[792, "val"] = 0 + model = GeologicalModel(bb[0, :], bb[1, :]) + model.data = data + upper = model.create_and_add_foliation("strati2", nelements=1000) + if onlap: + model.add_onlap_unconformity(upper, 1) + else: + model.add_unconformity(upper, 1) + model.create_and_add_foliation("strati", nelements=1000) + model.update() + return model, bb + + +def rebuild(model): + rebuilt = GeologicalModel.from_recipe_json(model.to_recipe_json()) + rebuilt.update() + return rebuilt + + +def same_field(a, b, name, points): + first, second = a.evaluate_feature_value(name, points), b.evaluate_feature_value(name, points) + assert np.array_equal(np.isnan(first), np.isnan(second)) + np.testing.assert_allclose(second, first, atol=1e-6, equal_nan=True) + + +# ------------------------------------------------------------------------------- faults +def test_a_fault_is_rebuilt_as_a_fault(): + model = fault_model() + rebuilt = rebuild(model) + assert [(f.name, type(f)) for f in rebuilt.features] == [(f.name, type(f)) for f in model.features] + assert type(rebuilt["fault"]).__name__ == "FaultSegment" + + +def test_the_faulted_field_and_the_displacement_are_reproduced(): + model = fault_model() + rebuilt = rebuild(model) + points = np.random.default_rng(0).uniform(0, 10, (300, 3)) + same_field(model, rebuilt, "strati", points) + np.testing.assert_allclose( + rebuilt.evaluate_fault_displacements(points), + model.evaluate_fault_displacements(points), + atol=1e-8, + ) + # the fault does something here, so the test is not vacuous + assert np.abs(model.evaluate_fault_displacements(points)).max() > 0.1 + + +def test_the_fault_parameters_are_stored(): + model = fault_model() + entry = json.loads(model.to_recipe_json())["model"]["features"][0] + assert entry["kind"] == "fault" + parameters = entry["build"]["parameters"] + assert entry["build"]["feature_type"] == "fault" + assert parameters["displacement"] == 1 + assert parameters["major_axis"] == 10 + assert parameters["faults"] == [] # a fault with no faults is not "all faults" + + +def test_a_foliation_given_no_faults_is_not_given_the_faults_when_rebuilt(): + """faults=[] means none; faults=None (the default) means every fault in the model.""" + without = fault_model(faults=[]) + with_faults = fault_model() + points = np.random.default_rng(1).uniform(0, 10, (300, 3)) + # the two models differ, so the test can tell the cases apart + assert not np.allclose( + without.evaluate_feature_value("strati", points), + with_faults.evaluate_feature_value("strati", points), + ) + same_field(without, rebuild(without), "strati", points) + same_field(with_faults, rebuild(with_faults), "strati", points) + + +def test_a_fault_that_refers_to_a_later_feature_is_refused(): + model = fault_model() + recipe = model.to_recipe_dict() + recipe["model"]["features"][1]["build"]["parameters"]["faults"] = ["not_there"] + with pytest.raises(LoopValueError, match="not earlier entries"): + GeologicalModel.from_recipe_dict(recipe) + + +# ------------------------------------------------------------------------ unconformities +@pytest.mark.parametrize("onlap", [False, True], ids=["erode", "onlap"]) +def test_an_unconformity_is_rebuilt_and_gives_the_same_model(onlap): + model, bb = two_group_model(onlap) + rebuilt = rebuild(model) + assert [(f.name, int(f.type)) for f in rebuilt.features] == [ + (f.name, int(f.type)) for f in model.features + ] + points = np.random.default_rng(0).uniform(bb[0], bb[1], (400, 3)) + same_field(model, rebuilt, "strati2", points) + same_field(model, rebuilt, "strati", points) + np.testing.assert_array_equal(rebuilt.evaluate_model(points), model.evaluate_model(points)) + + +@pytest.mark.parametrize("onlap", [False, True], ids=["erode", "onlap"]) +def test_the_unconformity_entry(onlap): + model, _ = two_group_model(onlap) + entry = json.loads(model.to_recipe_json())["model"]["features"][1] + assert entry["kind"] == "unconformity" + assert entry["build"] == { + "feature_type": "unconformity", + "parameters": {"feature": "strati2", "value": 1, "onlap": onlap}, + } + assert ModelSerializer.feature_recipe_kind(model.features[1]) == "unconformity" + + +def test_an_unconformity_that_refers_to_a_missing_feature_is_refused(): + model, _ = two_group_model(False) + recipe = model.to_recipe_dict() + recipe["model"]["features"][1]["build"]["parameters"]["feature"] = "not_there" + with pytest.raises(LoopValueError, match="not an earlier entry"): + GeologicalModel.from_recipe_dict(recipe) diff --git a/tests/unit/modelling/test_recipe_limits.py b/tests/unit/modelling/test_recipe_limits.py index 741693cf..6e9594d5 100644 --- a/tests/unit/modelling/test_recipe_limits.py +++ b/tests/unit/modelling/test_recipe_limits.py @@ -1,4 +1,8 @@ -"""What a recipe can and cannot store, and that it says so instead of rebuilding something wrong.""" +"""What a recipe can and cannot store, and that it says so instead of rebuilding something wrong. + +Recipes rebuild foliations, faults and unconformities. Anything else is refused when the recipe is +written, naming the feature. +""" import json import numpy as np @@ -91,15 +95,12 @@ def test_an_unconformity_is_recognised_as_an_unconformity(): # ------------------------------------------------------------------ what is refused -def test_a_model_with_a_fault_is_refused_not_rebuilt_without_the_fault(): - with pytest.raises(LoopValueError, match=r"fault.*foliations only"): - fault_model().to_recipe_json() - - -def test_a_model_with_an_unconformity_is_refused(): - model, _, _ = claudius_model() - model.add_unconformity(model["strati"], 0) - with pytest.raises(LoopValueError, match=r"__strati_unconformity.*unconformity"): +def test_a_domain_fault_is_refused(): + model = fault_model() + model.create_and_add_domain_fault("domain_fault_data", nelements=500) + # refused, and the message names the feature (it has no creation record: the domain fault in + # model.features is not the object create_and_add_feature returned) + with pytest.raises(LoopValueError, match=r"'domain_fault_data'.*cannot be stored in a recipe"): model.to_recipe_json() @@ -144,5 +145,5 @@ def test_an_older_recipe_entry_that_is_not_a_foliation_is_refused(): entry = recipe["model"]["features"][0] del entry["build"] entry["kind"] = "structural_frame" - with pytest.raises(LoopValueError, match="foliations only"): + with pytest.raises(LoopValueError, match="only foliations can be rebuilt"): GeologicalModel.from_recipe_dict(recipe) From 18436d52ad496b1068d2ad1faceb79db10b499e4 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 29 Oct 2026 15:32:03 +1030 Subject: [PATCH 4/5] test: a model with a fault and an unconformity rebuilds exactly Co-Authored-By: Claude Sonnet 5.5 --- .../test_recipe_faults_unconformities.py | 40 +++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/tests/unit/modelling/test_recipe_faults_unconformities.py b/tests/unit/modelling/test_recipe_faults_unconformities.py index 2d33fb1d..7121e1a1 100644 --- a/tests/unit/modelling/test_recipe_faults_unconformities.py +++ b/tests/unit/modelling/test_recipe_faults_unconformities.py @@ -152,3 +152,43 @@ def test_an_unconformity_that_refers_to_a_missing_feature_is_refused(): recipe["model"]["features"][1]["build"]["parameters"]["feature"] = "not_there" with pytest.raises(LoopValueError, match="not an earlier entry"): GeologicalModel.from_recipe_dict(recipe) + + +# --------------------------------------------------------------- both in one model +def test_a_model_with_a_fault_and_an_unconformity_rebuilds_exactly(): + rows = [ + [5, 5, 5, 0, 0.70710678, 0.0, 0.70710678, 0, "fault"], + [5, 5, 5, 0, -0.70710678, 0.0, 0.70710678, 1, "fault"], + [2, 2, 8, 0, 0.05, 0, 1, np.nan, "upper"], + [8, 2, 8.2, 0, 0.05, 0, 1, np.nan, "upper"], + [5, 8, 8.1, 0, 0, 0.05, 1, np.nan, "upper"], + [2, 5, 4, 0, 0, 0, 1, np.nan, "lower"], + [8, 5, 3.5, 0, 0, 0, 1, np.nan, "lower"], + [5, 2, 3.8, 0, 0, 0, 1, np.nan, "lower"], + [5, 8, 3.9, 0, 0, 0, 1, np.nan, "lower"], + [2, 5, 2, -2, 0, 0, 1, np.nan, "lower"], + ] + model = GeologicalModel([0, 0, 0], [10, 10, 10]) + model.data = pd.DataFrame( + rows, columns=["X", "Y", "Z", "val", "nx", "ny", "nz", "coord", "feature_name"] + ) + model.create_and_add_fault( + "fault", 1, nelements=1000, buffer=0.5, major_axis=10, minor_axis=3, intermediate_axis=10 + ) + upper = model.create_and_add_foliation("upper", nelements=1000) + model.add_unconformity(upper, 0) + model.create_and_add_foliation("lower", nelements=1000) + model.update() + + rebuilt = rebuild(model) + assert [(f.name, type(f)) for f in rebuilt.features] == [(f.name, type(f)) for f in model.features] + points = np.random.default_rng(0).uniform(0, 10, (500, 3)) + same_field(model, rebuilt, "upper", points) + same_field(model, rebuilt, "lower", points) + np.testing.assert_allclose( + rebuilt.evaluate_fault_displacements(points), + model.evaluate_fault_displacements(points), + atol=1e-8, + ) + # the unconformity hides part of each group, so the test is not vacuous + assert 0.1 < np.isnan(model.evaluate_feature_value("upper", points)).mean() < 0.95 From ee1e5c34a1a032ebe4b08c487f8bb09d956df089 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 29 Oct 2026 18:02:55 +1030 Subject: [PATCH 5/5] fix: deep-copy recipe parameters at creation; test mutation and array parameters Co-Authored-By: Claude Sonnet 5.5 --- .../modelling/core/_model_serializer.py | 6 ++++- tests/unit/modelling/test_recipe_limits.py | 27 +++++++++++++++++++ 2 files changed, 32 insertions(+), 1 deletion(-) diff --git a/LoopStructural/modelling/core/_model_serializer.py b/LoopStructural/modelling/core/_model_serializer.py index 70dbe545..5696c5ab 100644 --- a/LoopStructural/modelling/core/_model_serializer.py +++ b/LoopStructural/modelling/core/_model_serializer.py @@ -10,6 +10,7 @@ delegate to the staticmethods here. """ +import copy import json import pathlib @@ -76,7 +77,10 @@ def record_creation(feature: object, feature_type: str, params: dict) -> None: """ feature._recipe_creation = { "feature_type": feature_type, - "parameters": {k: v for k, v in params.items() if k not in _NOT_RECORDED}, + # deep copy: a list/dict/array the caller changes later must not change the recipe + "parameters": { + k: copy.deepcopy(v) for k, v in params.items() if k not in _NOT_RECORDED + }, "custom_data": params.get("data") is not None, } diff --git a/tests/unit/modelling/test_recipe_limits.py b/tests/unit/modelling/test_recipe_limits.py index 741693cf..fe9c673b 100644 --- a/tests/unit/modelling/test_recipe_limits.py +++ b/tests/unit/modelling/test_recipe_limits.py @@ -146,3 +146,30 @@ def test_an_older_recipe_entry_that_is_not_a_foliation_is_refused(): entry["kind"] = "structural_frame" with pytest.raises(LoopValueError, match="foliations only"): GeologicalModel.from_recipe_dict(recipe) + + +# ------------------------------------------------------- parameters are copied, not aliased +def test_changing_a_parameter_after_creation_does_not_change_the_recipe(): + model, _, _ = claudius_model() + values = [1, 2, 3] + weights = {"a": [1.0, 2.0]} + ModelSerializer.record_creation( + model["strati"], "foliation", {"nelements": 1000, "v": values, "w": weights} + ) + before = model.to_recipe_json() + values.append(4) + weights["a"].append(3.0) + weights["b"] = 1 + assert model.to_recipe_json() == before + parameters = json.loads(before)["model"]["features"][0]["build"]["parameters"] + assert parameters["v"] == [1, 2, 3] + assert parameters["w"] == {"a": [1.0, 2.0]} + + +def test_a_numpy_array_parameter_is_stored_as_a_list_and_copied(): + model, _, _ = claudius_model() + vector = np.array([0.0, 1.0, 0.5]) + ModelSerializer.record_creation(model["strati"], "foliation", {"nelements": 1000, "vec": vector}) + vector[0] = 99.0 + parameters = json.loads(model.to_recipe_json())["model"]["features"][0]["build"]["parameters"] + assert parameters["vec"] == [0.0, 1.0, 0.5]