fix: store build parameters in model recipes, and refuse what they cannot rebuild

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 |

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,99 @@
 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,
+            # 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,
+        }
+
+    @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 +150,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 +215,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

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(

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]:

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)