fix: address review findings on the derived refresh, native solver and export

loopstructural/gui/modelling/steps/checks.py

@@ -116,12 +116,13 @@ def check_model(data_manager, model_manager=None) -> StepCheck:
         problems.append("The model CRS must be a projected CRS (in metres).")
     changed = data_manager.get_changed_layers()
     if changed:
-        problems.append("Input layers changed after the model was built: " + ", ".join(changed))
+        # "Update data and solve" puts the new data into the features
+        todo.append("Input layers changed after the model was built: " + ", ".join(changed))
     state = model_manager.model_state if model_manager is not None else 'empty'
     if state == 'empty':
-        todo.append("Initialise the model.")
+        todo.append("Build the model.")
     elif state == 'stale':
-        problems.append("Inputs changed after the model was built. Initialise the model again.")
+        problems.append("Inputs changed after the model was built. Rebuild the model.")
     elif state == 'initialized':
         todo.append("Solve the model.")
     return StepCheck(tuple(problems), tuple(todo))

loopstructural/main/data_manager.py

@@ -1404,7 +1404,15 @@ def sync_extra_constraints(self):
                 if entry.get('processed') or entry.get('layer') is None:
                     continue
                 row = dict(entry)
-                row['df'] = qgsLayerToGeoDataFrame(entry['layer'], target_crs=model_crs)
+                try:
+                    row['df'] = qgsLayerToGeoDataFrame(entry['layer'], target_crs=model_crs)
+                except RuntimeError as error:
+                    # The layer was removed from the project: skip this row only
+                    self.logger(
+                        message=f"Could not read the added constraint '{key}' of '{name}': {error}",
+                        log_level=2,
+                    )
+                    continue
                 rows[key] = row
             if rows:
                 extras[name] = {'data': rows, 'use_z_coordinate': True}

loopstructural/main/derived_refresh.py

@@ -89,11 +89,9 @@ def _read_contacts_job(self):
         dm = self.data_manager
         settings = dm.get_widget_settings('basal_contacts_widget', {}) or {}
         geology = dm.get_layer_role(layer_roles.GEOLOGY)
-        faults = (
-            dm.find_layer_by_name(settings['faults_layer'])
-            if settings.get('faults_layer')
-            else dm.get_layer_role(layer_roles.FAULT_TRACES)
-        )
+        faults = dm.get_layer_role(layer_roles.FAULT_TRACES)
+        if faults is None and settings.get('faults_layer'):
+            faults = dm.find_layer_by_name(settings['faults_layer'])
         unit_field = dm.get_layer_role(layer_roles.GEOLOGY_UNIT_FIELD)
         inputs = dm.basal_contacts_inputs()
         _check_contacts_inputs(inputs)
@@ -126,10 +124,12 @@ def _read_thickness_job(self):
             )
 
         def layer(key, role=None):
-            name = settings.get(key)
-            if name:
-                return dm.find_layer_by_name(name)
-            return dm.get_layer_role(role) if role else None
+            # The role first, as in `thickness_inputs`, so that the layer of a
+            # run is the layer that its inputs record.
+            found = dm.get_layer_role(role) if role else None
+            if found is None and settings.get(key):
+                found = dm.find_layer_by_name(settings[key])
+            return found
 
         geology = layer('geology_layer', layer_roles.GEOLOGY)
         structure = layer('structure_layer', layer_roles.STRUCTURE)

loopstructural/main/model_export.py

@@ -5,6 +5,7 @@
 and `get_stratigraphic_unit_names`) and write files with pyvista.
 """
 
+import csv
 import re
 from pathlib import Path
 
@@ -73,7 +74,8 @@ def _check_extension(extension, formats):
 def unit_names_for_ids(model_manager, ids) -> np.ndarray:
     """Return the unit name for each stratigraphic id. An id of -1 has no unit."""
     names = list(model_manager.get_stratigraphic_unit_names())
-    ids = np.asarray(ids, dtype=int)
+    # NaN (a point with no unit) has the id -1
+    ids = np.nan_to_num(np.asarray(ids, dtype=float), nan=-1).astype(int)
     return np.array([names[i] if 0 <= i < len(names) else '' for i in ids], dtype=object)
 
 
@@ -84,6 +86,9 @@ def add_stratigraphy(model_manager, mesh, *, cells=False):
     """
     points = mesh.cell_centers().points if cells else mesh.points
     ids = np.asarray(model_manager.evaluate_stratigraphy_on_points(points))
+    if ids.dtype.kind == 'f':
+        # a point with no unit has the id -1
+        ids = np.nan_to_num(ids, nan=-1).astype(int)
     data = mesh.cell_data if cells else mesh.point_data
     data[STRATIGRAPHY_ID_FIELD] = ids
     # a VTK string array: '' where there is no unit
@@ -108,9 +113,11 @@ def existing_surface_files(folder, extension):
 def _write_csv(path, points, ids, names):
     points = np.asarray(points)
     with open(path, 'w', encoding='utf-8', newline='') as handle:
-        handle.write(f"x,y,z,{STRATIGRAPHY_ID_FIELD},{UNIT_FIELD}\n")
+        writer = csv.writer(handle, lineterminator='\n')
+        writer.writerow(['x', 'y', 'z', STRATIGRAPHY_ID_FIELD, UNIT_FIELD])
         for (x, y, z), unit_id, name in zip(points, ids, names):
-            handle.write(f"{x:.6f},{y:.6f},{z:.6f},{int(unit_id)},\"{name}\"\n")
+            unit_id = -1 if np.isnan(unit_id) else int(unit_id)
+            writer.writerow([f"{x:.6f}", f"{y:.6f}", f"{z:.6f}", unit_id, name])
 
 
 def export_surfaces(

loopstructural/main/model_manager.py

@@ -1104,6 +1104,15 @@ def _solve_feature(self, feature, rebuild=True):
             self._use_fallback(name, builders, unconverged, fallback, "did not converge")
             feature.builder.up_to_date()
 
+    def _forget_native_state(self, name=None):
+        """Let the feature `name` (or every feature) try the native solver again."""
+        if name is None:
+            self._native_failed.clear()
+            self._native_fallbacks.clear()
+            return
+        self._native_failed.discard(name)
+        self._native_fallbacks.pop(name, None)
+
     @staticmethod
     def _set_fallback(builder, fallback):
         # the native iteration limit is not for the other solvers
@@ -1155,9 +1164,10 @@ def _model_area(self) -> float:
         """Return the map (XY) area of the model box, or 0 if the box is not set."""
         bounding_box = getattr(self.model, 'bounding_box', None)
         try:
-            lengths = np.asarray(bounding_box.maximum, dtype=float)[:2] - np.asarray(
-                bounding_box.origin, dtype=float
-            )[:2]
+            lengths = (
+                np.asarray(bounding_box.maximum, dtype=float)[:2]
+                - np.asarray(bounding_box.origin, dtype=float)[:2]
+            )
         except (AttributeError, TypeError, ValueError, IndexError):
             return 0.0
         area = float(np.prod(lengths))
@@ -1422,9 +1432,7 @@ def _group_by_name(self, name):
     def detached_to_dict(self) -> dict:
         """Return `detached` in a form that `json.dump` can write."""
         return {
-            name: {
-                'columns': {str(c): [_json_value(v) for v in data[c]] for c in data.columns}
-            }
+            name: {'columns': {str(c): [_json_value(v) for v in data[c]] for c in data.columns}}
             for name, data in self.detached.items()
         }
 
@@ -1659,7 +1667,13 @@ def mark_data_dirty(self):
     def remove_parametric_fault(self, name: str):
         """Stop building the parametric fault `name`. The model is stale until the next build."""
         if self.parametric_faults.pop(name, None) is not None:
+            self.nelements_overrides.pop(name, None)
+            self.nelements_used.pop(name, None)
+            self._data_summaries.pop(name, None)
+            self._needs_inequalities.discard(name)
+            self._forget_native_state(name)
             self._data_dirty = True
+            self._emit('model_updated')
 
     def _build_parametric_faults(self):
         for name, spec in self.parametric_faults.items():
@@ -1677,13 +1691,22 @@ def _build_parametric_faults(self):
 
     def parametric_faults_to_dict(self) -> dict:
         """Return `parametric_faults` in a form that `json.dump` can write."""
-        return {name: parametric_fault.clean_spec(spec) for name, spec in self.parametric_faults.items()}
+        return {
+            name: parametric_fault.clean_spec(spec) for name, spec in self.parametric_faults.items()
+        }
 
     def parametric_faults_from_dict(self, faults: dict):
         """Replace `parametric_faults` with the specs from `parametric_faults_to_dict`."""
-        self.parametric_faults = {
-            name: parametric_fault.spec_from_dict(spec) for name, spec in (faults or {}).items()
-        }
+        self.parametric_faults = {}
+        for name, spec in (faults or {}).items():
+            # one damaged fault must not stop the restore of the project
+            try:
+                self.parametric_faults[name] = parametric_fault.spec_from_dict(spec)
+            except ValueError as error:
+                if self._debug_manager is not None:
+                    self._debug_manager.log(
+                        f"Skipped the saved parametric fault '{name}': {error}", log_level=2
+                    )
 
     def _get_feature_by_name_or_none(self, name):
         """Non-raising counterpart to `GeologicalModel.get_feature_by_name`.
@@ -1883,6 +1906,8 @@ def update_model(
         self.model.features = []
         self.model.feature_name_index = {}
         self.skipped_constraints = {}
+        # The data is read again, so a feature can use the native solver again
+        self._forget_native_state()
 
         # Notify start (if requested) so UI can react
         if notify_observers:
@@ -1980,6 +2005,8 @@ def refresh_feature_data(self) -> dict:
         """
         updated = []
         needs_initialize = []
+        # The data paths below add the features that need inequalities again
+        self._needs_inequalities.clear()
 
         def refresh(name, data):
             feature = self._get_feature_by_name_or_none(name)
@@ -2000,6 +2027,8 @@ def refresh(name, data):
             if isinstance(current, pd.DataFrame) and current.equals(data):
                 return
             target.add_data_from_data_frame(data)
+            # new data can be solved by the native solver
+            self._forget_native_state(name)
             # add_data_from_data_frame does not reset this flag, and while
             # it is set the builder keeps the constraints of the old data
             target.data_added = False

loopstructural/main/parametric_fault.py

@@ -174,16 +174,34 @@ def clean_spec(spec: dict) -> dict:
     """Return the spec as plain numbers and text, safe for JSON."""
     result = {key: spec[key] for key in FIELDS if key in spec}
     result['name'] = str(result['name']).strip()
-    for key in ('strike', 'dip', 'pitch', 'displacement', 'major_axis', 'intermediate_axis',
-                'minor_axis'):
+    for key in (
+        'strike',
+        'dip',
+        'pitch',
+        'displacement',
+        'major_axis',
+        'intermediate_axis',
+        'minor_axis',
+    ):
         if key in result:
             result[key] = float(result[key])
     result['centre'] = [float(v) for v in result['centre']]
     return result
 
 
 def spec_from_dict(data: Optional[dict]) -> dict:
-    """Read a spec that was saved. The centre is a tuple again."""
-    spec = clean_spec(data)
+    """Read a spec that was saved. The centre is a tuple again.
+
+    Raises
+    ------
+    ValueError
+        If the saved spec is empty or has no name or centre.
+    """
+    if not isinstance(data, dict) or 'name' not in data or 'centre' not in data:
+        raise ValueError("The saved parametric fault has no name or centre.")
+    try:
+        spec = clean_spec(data)
+    except (TypeError, ValueError) as error:
+        raise ValueError(f"The saved parametric fault is not valid: {error}") from error
     spec['centre'] = tuple(spec['centre'])
     return spec

tests/unit/test_model_export.py

@@ -129,8 +129,8 @@ def test_export_block_model_csv(tmp_path):
     )
     lines = path.read_text(encoding='utf-8').splitlines()
     assert lines[0] == 'x,y,z,stratigraphy_id,unit'
-    assert lines[1].endswith('0,"lower"')
-    assert lines[2].endswith('1,"upper"')
+    assert lines[1].endswith('0,lower')
+    assert lines[2].endswith('1,upper')
 
 
 def test_export_cross_section_plane_and_line(tmp_path):
@@ -153,3 +153,21 @@ def test_export_cross_section_needs_input(tmp_path):
         model_export.export_cross_section(
             manager, tmp_path / 'x.vtk', origin=(0, 0, 0), normal=(0, 0, 0), size=1
         )
+
+
+def test_csv_escapes_a_quote_in_a_unit_name(tmp_path):
+    import csv
+
+    path = tmp_path / 'quote.csv'
+    model_export._write_csv(path, [(0, 0, 0)], [3], ['Unit "A"'])
+    with open(path, newline='', encoding='utf-8') as handle:
+        rows = list(csv.reader(handle))
+    assert rows[1][3:] == ['3', 'Unit "A"']
+
+
+def test_nan_id_has_no_unit():
+    import numpy as np
+
+    manager = FakeManager(FakeModel())
+    names = model_export.unit_names_for_ids(manager, np.array([0.0, np.nan]))
+    assert names[1] == ''

tests/unit/test_parametric_fault.py

@@ -146,3 +146,9 @@ def test_the_default_spec_is_in_the_middle_of_the_area(self):
         assert spec['major_axis'] == 1000.0
         assert spec['minor_axis'] == 500.0
         assert parametric_fault.problems({**spec, 'name': 'F'}) == []
+
+
+@pytest.mark.parametrize('saved', [None, {}, {'name': 'F'}, {'centre': [0, 0, 0]}])
+def test_spec_from_dict_rejects_a_damaged_spec(saved):
+    with pytest.raises(ValueError):
+        parametric_fault.spec_from_dict(saved)

tests/unit/test_step_checks.py

@@ -110,7 +110,10 @@ def test_a_geographic_crs_is_a_problem(self, dm):
 class TestStratigraphyStep:
     @pytest.fixture(autouse=True)
     def map_layers(self, dm, request):
-        if request.node.name not in ('test_an_empty_column_is_not_started', 'test_missing_map_layers_are_a_todo'):
+        if request.node.name not in (
+            'test_an_empty_column_is_not_started',
+            'test_missing_map_layers_are_a_todo',
+        ):
             dm.layer_roles.values = {
                 layer_roles.GEOLOGY: object(),
                 layer_roles.STRUCTURE: object(),
@@ -241,12 +244,13 @@ def test_out_of_date_data_is_shown_on_step_2_only(self, dm):
         assert checks.check_stratigraphy(dm).status == Status.PROBLEM
         assert checks.check_model(dm, FakeModelManager('solved')).status == Status.DONE
 
-    def test_changed_layers_are_a_problem(self, dm):
+    def test_changed_layers_are_a_todo(self, dm):
+        # "Update data and solve" fixes this with one click, so it is not a problem
         dm.bounding_box_set = True
         dm.changed_layers = ['geology']
         check = checks.check_model(dm, FakeModelManager('solved'))
-        assert check.status == Status.PROBLEM
-        assert 'geology' in check.problems[0]
+        assert not check.problems
+        assert 'geology' in check.todo[0]
 
 
 class TestViewStep: