diff --git a/COMPAT.md b/COMPAT.md index 11ee6262..99e60e1e 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. 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 cbf4d112..5cd08331 100644 --- a/LoopStructural/modelling/core/_model_serializer.py +++ b/LoopStructural/modelling/core/_model_serializer.py @@ -10,9 +10,11 @@ delegate to the staticmethods here. """ +import copy import json import pathlib +import numpy as np import pandas as pd from ...geometry import BoundingBox @@ -24,19 +26,109 @@ logger = getLogger(__name__) +# 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, 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: + """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). + """ + # 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 + } + 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": parameters, + "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 _SUPPORTED_FEATURE_TYPES.get(feature_type) != kind: + raise LoopValueError( + f"Feature '{feature.name}' ({kind}, created as '{feature_type}') cannot be " + f"stored in a recipe yet: recipes support {sorted(_SUPPORTED_FEATURE_TYPES)}." + ) + 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 +160,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 +225,55 @@ 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: 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 " + "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_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 + ) + if "value" not in parameters: + raise LoopValueError( + f"Unconformity '{feature_name}' has no 'value' in its recipe parameters" + ) + feature = add(base, parameters["value"]) + if parameters.get("onlap"): + # add_onlap_unconformity returns the unconformity applied to the other + # features; the model holds its inverse + feature = model.get_feature_by_name(feature_name) + 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 @@ -145,7 +283,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/LoopStructural/modelling/core/geological_model.py b/LoopStructural/modelling/core/geological_model.py index cabaff0c..03a4d4d7 100644 --- a/LoopStructural/modelling/core/geological_model.py +++ b/LoopStructural/modelling/core/geological_model.py @@ -708,7 +708,18 @@ 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) + record_params = dict(params) + if feature_type == "foliation" and record_params.get("faults") is None: + # faults=None means every fault in the model when the feature is built, which is not + # recoverable from the feature list order once index= has been used: record the + # faults that are actually used (in model order, as add_faults reads them) + record_params["faults"] = [ + f.name for f in self.features if f.type == FeatureType.FAULT + ] + feature = FeatureBuilderRegistry.create(feature_type, self, name, **params) + if feature is not None: + ModelSerializer.record_creation(feature, feature_type, record_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/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_faults_unconformities.py b/tests/unit/modelling/test_recipe_faults_unconformities.py new file mode 100644 index 00000000..e9a20278 --- /dev/null +++ b/tests/unit/modelling/test_recipe_faults_unconformities.py @@ -0,0 +1,244 @@ +"""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) + + +# --------------------------------------------------------------- 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 + + +# ------------------------------------------------------------------ resolved names, not order +def test_a_model_built_with_index_records_the_faults_it_used_and_rebuilds_identically(): + """A foliation built before a fault that is then inserted ahead of it has no faults. + + Rebuilding in list order would give the foliation the fault (faults=None means every fault + in the model at that time), so the recipe records the faults that were actually used. + """ + model = GeologicalModel([0, 0, 0], [10, 10, 10]) + model.data = fault_data() + model.create_and_add_foliation("strati", nelements=1000) + model.create_and_add_fault( + "fault", + 1, + nelements=1000, + buffer=0.5, + major_axis=10, + minor_axis=3, + intermediate_axis=10, + index=0, + ) + model.update() + assert [f.name for f in model.features] == ["fault", "strati"] + entries = json.loads(model.to_recipe_json())["model"]["features"] + assert entries[1]["build"]["parameters"]["faults"] == [] + rebuilt = rebuild(model) + points = np.random.default_rng(2).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, + ) + + +def test_an_unconformity_without_a_value_is_refused_with_a_loop_error(): + model, _ = two_group_model(False) + recipe = model.to_recipe_dict() + del recipe["model"]["features"][1]["build"]["parameters"]["value"] + with pytest.raises(LoopValueError, match="no 'value'"): + GeologicalModel.from_recipe_dict(recipe) + + +def test_an_onlap_unconformity_is_mapped_to_the_feature_in_the_model(): + model, _ = two_group_model(True) + rebuilt = rebuild(model) + assert [f.name for f in rebuilt.features] == [f.name for f in model.features] + for name in (f.name for f in rebuilt.features): + assert rebuilt[name] is rebuilt.features[rebuilt.feature_name_index[name]] diff --git a/tests/unit/modelling/test_recipe_limits.py b/tests/unit/modelling/test_recipe_limits.py new file mode 100644 index 00000000..b488603d --- /dev/null +++ b/tests/unit/modelling/test_recipe_limits.py @@ -0,0 +1,178 @@ +"""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 +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", "index"} & set(parameters) + # faults are stored as the names of the faults that were used, never as features + assert parameters["faults"] == [] + + +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_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() + + +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="only foliations can be rebuilt"): + 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] 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..a6ac2bb7 --- /dev/null +++ b/tests/unit/modelling/test_recipe_reproduces_field.py @@ -0,0 +1,46 @@ +"""A model recipe should rebuild a model that gives the same scalar field. + +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 + +from LoopStructural import GeologicalModel +from LoopStructural.datasets import load_claudius + +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, "regularisation": 10.0, "damp": False}, id="several-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) +