From ef50446491ff35fd4729bfa50fb8847d11841d54 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Sat, 31 Oct 2026 13:20:58 +1030 Subject: [PATCH 1/4] feat: use the native solver of LoopStructural when it is installed - New setting Solver: Native (the default) or Default. With Native, and when loop-interpolation-core is installed, the interpolators get solver='native'. Without the core, the plugin uses the default solver. - The native solver also solves the inequality constraints, so ADMM is not set for them when native is used. - A feature that the native solver cannot solve (for example a fold or a region mask) is built again with the default solver (cg, or ADMM for inequalities). The plugin uses the default solver for that feature from then on, and writes one log message. --- loopstructural/gui/dlg_settings.py | 26 ++++++ loopstructural/main/model_manager.py | 73 +++++++++++++++-- loopstructural/main/solver_choice.py | 63 +++++++++++++++ loopstructural/toolbelt/preferences.py | 2 + tests/qgis/test_native_solver.py | 108 +++++++++++++++++++++++++ tests/unit/test_solver_choice.py | 51 ++++++++++++ 6 files changed, 316 insertions(+), 7 deletions(-) create mode 100644 loopstructural/main/solver_choice.py create mode 100644 tests/qgis/test_native_solver.py create mode 100644 tests/unit/test_solver_choice.py diff --git a/loopstructural/gui/dlg_settings.py b/loopstructural/gui/dlg_settings.py index d719475..a7461a6 100644 --- a/loopstructural/gui/dlg_settings.py +++ b/loopstructural/gui/dlg_settings.py @@ -14,6 +14,7 @@ from qgis.PyQt import uic from qgis.PyQt.QtCore import QUrl from qgis.PyQt.QtGui import QDesktopServices, QIcon +from qgis.PyQt.QtWidgets import QComboBox, QLabel # project from loopstructural.__about__ import ( @@ -23,6 +24,7 @@ __uri_tracker__, __version__, ) +from loopstructural.main import solver_choice from loopstructural.toolbelt import PlgLogger, PlgOptionsManager from loopstructural.toolbelt.preferences import PlgSettingsStructure @@ -84,6 +86,7 @@ def __init__(self, parent): self.n_elements_auto_check_box.toggled.connect( lambda checked: self.n_elements_spin_box.setEnabled(not checked) ) + self._add_solver_combo() # load previously saved settings self.load_settings() @@ -100,6 +103,7 @@ def apply(self): settings.interpolator_npw = self.npw_spin_box.value() settings.interpolator_cpw = self.cpw_spin_box.value() settings.interpolator_regularisation = self.regularisation_spin_box.value() + settings.interpolator_solver = self.solver_combo.currentData() settings.version = __version__ debug_dir_text = ( self.le_debug_directory.text() if hasattr(self, "le_debug_directory") else "" @@ -131,9 +135,31 @@ def load_settings(self): self.regularisation_spin_box.setValue(settings.interpolator_regularisation) self.cpw_spin_box.setValue(settings.interpolator_cpw) self.npw_spin_box.setValue(settings.interpolator_npw) + index = self.solver_combo.findData(settings.interpolator_solver) + self.solver_combo.setCurrentIndex(max(index, 0)) if hasattr(self, "le_debug_directory"): self.le_debug_directory.setText(settings.debug_directory or "") + def _add_solver_combo(self): + """Add the choice of the solver below the number of elements.""" + self.solver_combo = QComboBox(self) + for choice in solver_choice.CHOICES: + self.solver_combo.addItem(solver_choice.LABELS[choice], choice) + if solver_choice.native_available(): + tooltip = ( + "The native solver is much faster. A feature that it cannot solve " + "(for example a fold) is solved with the default solver." + ) + else: + tooltip = ( + "The native solver is not installed, so the default solver is used. " + "Install it with: pip install loopstructural[native]" + ) + self.solver_combo.setToolTip(tooltip) + row = self.gridLayout_solver.rowCount() + self.gridLayout_solver.addWidget(QLabel("Solver", self), row, 0) + self.gridLayout_solver.addWidget(self.solver_combo, row, 1) + def reset_settings(self): """Reset settings in the UI and persisted settings to plugin defaults.""" default_settings = PlgSettingsStructure() diff --git a/loopstructural/main/model_manager.py b/loopstructural/main/model_manager.py index 346b7e8..cfecfdf 100644 --- a/loopstructural/main/model_manager.py +++ b/loopstructural/main/model_manager.py @@ -29,7 +29,7 @@ from LoopStructural import GeologicalModel from loopstructural.toolbelt.preferences import PlgOptionsManager -from ..main import constraints, parametric_fault +from ..main import constraints, parametric_fault, solver_choice from ..main.data_types import FaultEntry, StratigraphyEntry from ..main.helpers import qgisAttributeIsNone from ..main.interpolation_size import DataSummary, suggest_nelements, summarise_data @@ -248,6 +248,11 @@ def __init__(self, debug_manager=None): self.nelements_used: Dict[str, tuple] = {} # feature name -> DataSummary of the data of the last build self._data_summaries: Dict[str, DataSummary] = {} + # features that the native solver could not solve; they use the + # default solver (see `_solve_feature`) + self._native_failed: Set[str] = set() + # features with constraints that need an inequality solver + self._needs_inequalities: Set[str] = set() @property def uses_map_data(self) -> bool: @@ -325,6 +330,8 @@ def reset(self): self.nelements_overrides = {} self.nelements_used = {} self._data_summaries = {} + self._native_failed = set() + self._needs_inequalities = set() self.dem_function = lambda x, y: 0 self._topology_dirty = False self._data_dirty = False @@ -1015,12 +1022,64 @@ def interpolator_arguments(self, name, data) -> dict: else: nelements, automatic = self.default_nelements(name) self.nelements_used[name] = (nelements, automatic) - return { + arguments = { 'nelements': nelements, 'npw': settings.interpolator_npw, 'cpw': settings.interpolator_cpw, 'regularisation': settings.interpolator_regularisation, } + if self._use_native(name): + arguments['solver'] = solver_choice.NATIVE + return arguments + + def _use_native(self, name) -> bool: + """Return True if the feature `name` uses the native solver. + + Not for a feature that the native solver could not solve before, see + `_solve_feature`. + """ + if name in self._native_failed: + return False + settings = PlgOptionsManager.get_plg_settings() + return solver_choice.use_native(settings.interpolator_solver) + + def _inequality_solver(self, name, solver): + """Return the solver to pass for a constraint type that needs `solver` (or None). + + The native solver also solves the inequality constraints, so with it + `interpolator_arguments` gives the solver and this returns None. + """ + if not solver: + return None + self._needs_inequalities.add(name) + return None if self._use_native(name) else solver + + def _solve_feature(self, feature, rebuild=True): + """Build `feature`. If the native solver cannot solve it, build it with + the default solver, and use the default solver for it from now on. + + `rebuild` builds the feature again; otherwise only if it is out of date. + """ + try: + if rebuild: + feature.builder.update() + else: + feature.builder.up_to_date() + except solver_choice.NATIVE_ERRORS as error: + if not solver_choice.uses_native(feature): + raise + name = feature.name + self._native_failed.add(name) + fallback = solver_choice.fallback_solver(name in self._needs_inequalities) + if self._debug_manager is not None: + self._debug_manager.log( + f"The native solver cannot solve '{name}' ({error}). " + f"It is solved with '{fallback}'.", + log_level=1, + ) + for builder in solver_choice.builders(feature): + builder.update_build_arguments({'solver': fallback}) + feature.builder.update() def default_nelements(self, name) -> tuple: """Return (number of elements, True if automatic) for the feature @@ -1238,7 +1297,7 @@ def _extra_kwargs(self, name) -> dict: return {} kwargs = {} for layer_data in spec['data'].values(): - solver = constraints.solver_for(layer_data.get('type')) + solver = self._inequality_solver(name, constraints.solver_for(layer_data.get('type'))) if solver: kwargs['solver'] = solver return kwargs @@ -1941,7 +2000,7 @@ def update_feature(self, feature_name: str): raise ValueError(f"Feature '{feature_name}' not found in the model.") # Allow UI to react to a feature update self._emit('model_update_started') - feature.builder.update() + self._solve_feature(feature) # Notify observers and include feature name for interested listeners self._emit('feature_updated', feature_name) self._emit('model_update_finished') @@ -2003,7 +2062,7 @@ def update_all_features( feature = self.model.get_feature_by_name(feature_name) if feature is not None: self._report_progress(f"Solving '{feature_name}'") - feature.builder.update() + self._solve_feature(feature) finally: self._progress_callback = None else: @@ -2020,7 +2079,7 @@ def update_all_features( for f in self.model.features: if not f.name.startswith('__'): self._report_progress(f"Solving '{f.name}'") - f.builder.up_to_date() + self._solve_feature(f, rebuild=False) finally: self._progress_callback = None if dbg is not None: @@ -2300,7 +2359,7 @@ def _foliation_data( dfs.append(constraints.add_weight(rows, df, layer_data)) else: dfs.append(constraints.constraint_rows(df, layer_data, name)) - solver = constraints.solver_for(layer_type) + solver = self._inequality_solver(name, constraints.solver_for(layer_type)) if solver: kwargs['solver'] = solver except Exception as e: diff --git a/loopstructural/main/solver_choice.py b/loopstructural/main/solver_choice.py new file mode 100644 index 0000000..7c22634 --- /dev/null +++ b/loopstructural/main/solver_choice.py @@ -0,0 +1,63 @@ +"""Choose the solver of the interpolators: the native (Rust) solver or the default. + +The native solver is ``solver='native'`` of loop-interpolation 0.2 and later. +It needs the package loop-interpolation-core (``pip install +loopstructural[native]``). It does not support all settings, for example a +fold interpolator or a region mask. For these, the model manager builds the +feature again with the default solver (see `fallback_solver`). + +The module does not import QGIS, so the tests run in the fast tests/unit/ job. +""" + +from importlib.util import find_spec + +DEFAULT = 'default' +NATIVE = 'native' +CHOICES = (DEFAULT, NATIVE) +LABELS = { + DEFAULT: "Default (Python)", + NATIVE: "Native (fast, if it is installed)", +} + +# The errors that tell that the native solver cannot solve a feature +NATIVE_ERRORS = (NotImplementedError, ImportError) + + +def native_available() -> bool: + """Return True if loop-interpolation has the native solver and its core is installed.""" + try: + return ( + find_spec('loop_interpolation_core') is not None + and find_spec('loop_interpolation._native') is not None + ) + except (ImportError, ValueError): + return False + + +def use_native(setting: str) -> bool: + """Return True if the setting asks for the native solver and it is installed.""" + return setting == NATIVE and native_available() + + +def fallback_solver(needs_inequalities: bool): + """Return the solver for a feature that the native solver cannot solve. + + Inequality constraints need ADMM. Otherwise the default of LoopStructural (cg). + """ + return 'admm' if needs_inequalities else 'cg' + + +def builders(feature): + """Return the builders of a feature: one, or the three of a structural frame.""" + builder = getattr(feature, 'builder', None) + if builder is None: + return [] + return list(getattr(builder, 'builders', None) or [builder]) + + +def uses_native(feature) -> bool: + """Return True if a builder of the feature uses the native solver.""" + return any( + (getattr(b, 'build_arguments', None) or {}).get('solver') == NATIVE + for b in builders(feature) + ) diff --git a/loopstructural/toolbelt/preferences.py b/loopstructural/toolbelt/preferences.py index 196ed4e..84c6682 100644 --- a/loopstructural/toolbelt/preferences.py +++ b/loopstructural/toolbelt/preferences.py @@ -36,6 +36,8 @@ class PlgSettingsStructure: interpolator_regularisation: float = 1.0 interpolator_cpw: float = 1.0 interpolator_npw: float = 1.0 + # 'native' or 'default'; see `main.solver_choice` + interpolator_solver: str = 'native' separate_dock_widgets: bool = False # map units; see `fault_topology_calc.fault_pair_direction` fault_end_tolerance: float = 10.0 diff --git a/tests/qgis/test_native_solver.py b/tests/qgis/test_native_solver.py new file mode 100644 index 0000000..4fe479a --- /dev/null +++ b/tests/qgis/test_native_solver.py @@ -0,0 +1,108 @@ +"""Pytest tests for the native solver in the model manager. + +The features are small stand-ins: the tests check the choice of the solver and +the fallback, not LoopStructural. +""" + +import pytest + +from loopstructural.main import model_manager as model_manager_module +from loopstructural.main import solver_choice +from loopstructural.main.model_manager import GeologicalModelManager +from loopstructural.toolbelt.preferences import PlgSettingsStructure + + +class _DebugManager: + def __init__(self): + self.messages = [] + + def log(self, message, *args, **kwargs): + self.messages.append(message) + + +class _Builder: + def __init__(self, solver): + self.build_arguments = {'solver': solver} + self.solvers = [] + + def update_build_arguments(self, arguments): + self.build_arguments.update(arguments) + + def update(self): + solver = self.build_arguments.get('solver') + self.solvers.append(solver) + if solver == solver_choice.NATIVE: + raise NotImplementedError("solver='native' does not support the setting 'region'") + + up_to_date = update + + +class _Feature: + def __init__(self, name, solver=solver_choice.NATIVE): + self.name = name + self.builder = _Builder(solver) + + +def _use(monkeypatch, solver, available=True): + settings = PlgSettingsStructure(interpolator_solver=solver) + monkeypatch.setattr( + model_manager_module.PlgOptionsManager, 'get_plg_settings', staticmethod(lambda: settings) + ) + monkeypatch.setattr(solver_choice, 'native_available', lambda: available) + + +@pytest.fixture +def manager(): + return GeologicalModelManager(debug_manager=_DebugManager()) + + +def test_native_is_passed_when_it_is_chosen_and_installed(manager, monkeypatch): + _use(monkeypatch, solver_choice.NATIVE) + assert manager.interpolator_arguments('s1', None)['solver'] == 'native' + + +def test_no_solver_is_passed_without_native(manager, monkeypatch): + _use(monkeypatch, solver_choice.DEFAULT) + assert 'solver' not in manager.interpolator_arguments('s1', None) + _use(monkeypatch, solver_choice.NATIVE, available=False) + assert 'solver' not in manager.interpolator_arguments('s1', None) + + +def test_native_also_solves_the_inequalities(manager, monkeypatch): + _use(monkeypatch, solver_choice.NATIVE) + assert manager._inequality_solver('s1', 'admm') is None + _use(monkeypatch, solver_choice.DEFAULT) + assert manager._inequality_solver('s1', 'admm') == 'admm' + assert manager._inequality_solver('s1', None) is None + + +def test_a_feature_that_native_cannot_solve_uses_the_default(manager, monkeypatch): + _use(monkeypatch, solver_choice.NATIVE) + feature = _Feature('s1') + manager._solve_feature(feature) + assert feature.builder.solvers == ['native', 'cg'] + assert 's1' in manager._native_failed + assert any('native' in message for message in manager._debug_manager.messages) + # the next build of the feature does not try native again + assert 'solver' not in manager.interpolator_arguments('s1', None) + + +def test_the_fallback_of_inequalities_is_admm(manager, monkeypatch): + _use(monkeypatch, solver_choice.NATIVE) + manager._inequality_solver('s1', 'admm') + feature = _Feature('s1') + manager._solve_feature(feature, rebuild=False) + assert feature.builder.solvers == ['native', 'admm'] + + +def test_an_error_without_native_is_raised(manager, monkeypatch): + _use(monkeypatch, solver_choice.DEFAULT) + + class _Failing(_Builder): + def update(self): + raise NotImplementedError("something else") + + feature = _Feature('s1', solver='cg') + feature.builder = _Failing('cg') + with pytest.raises(NotImplementedError): + manager._solve_feature(feature) diff --git a/tests/unit/test_solver_choice.py b/tests/unit/test_solver_choice.py new file mode 100644 index 0000000..f43956d --- /dev/null +++ b/tests/unit/test_solver_choice.py @@ -0,0 +1,51 @@ +"""Pytest tests for the choice of the solver (main/solver_choice.py). + +The module does not import QGIS, so the tests run in the fast tests/unit/ job. +""" + +from types import SimpleNamespace + +from loopstructural.main import solver_choice + + +def _builder(solver=None): + return SimpleNamespace(build_arguments={'solver': solver} if solver else {}) + + +def test_native_is_used_only_if_it_is_installed(monkeypatch): + monkeypatch.setattr(solver_choice, 'native_available', lambda: True) + assert solver_choice.use_native(solver_choice.NATIVE) + assert not solver_choice.use_native(solver_choice.DEFAULT) + monkeypatch.setattr(solver_choice, 'native_available', lambda: False) + assert not solver_choice.use_native(solver_choice.NATIVE) + + +def test_native_available_is_a_bool(): + assert isinstance(solver_choice.native_available(), bool) + + +def test_fallback_solver(): + assert solver_choice.fallback_solver(needs_inequalities=True) == 'admm' + assert solver_choice.fallback_solver(needs_inequalities=False) == 'cg' + + +def test_builders_of_a_feature_and_of_a_frame(): + one = _builder() + assert solver_choice.builders(SimpleNamespace(builder=one)) == [one] + parts = [_builder(), _builder(), _builder()] + frame = SimpleNamespace(builder=SimpleNamespace(builders=parts)) + assert solver_choice.builders(frame) == parts + assert solver_choice.builders(SimpleNamespace()) == [] + + +def test_uses_native(): + assert solver_choice.uses_native(SimpleNamespace(builder=_builder('native'))) + assert not solver_choice.uses_native(SimpleNamespace(builder=_builder('cg'))) + frame = SimpleNamespace( + builder=SimpleNamespace(builders=[_builder('cg'), _builder('native'), _builder('cg')]) + ) + assert solver_choice.uses_native(frame) + + +def test_every_choice_has_a_label(): + assert set(solver_choice.LABELS) == set(solver_choice.CHOICES) From 459027658e9f71530b7736786f91798c55386723 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Sun, 1 Nov 2026 15:10:18 +1030 Subject: [PATCH 2/4] build: require LoopStructural[native] 1.11 and fix the pin bump workflow - LoopStructural 1.11.0 allows loop-interpolation 0.2, which has the native solver. Require LoopStructural[native]~=1.11.0 in the three requirement files. - The pin bump workflow read the key "LoopStructural" from the release manifest, but the key of the root package is ".", so it always stopped. It also did not update loopstructural/requirements.txt (the file that qpip reads) and did not match a pin with an extra. --- .github/workflows/bump-loopstructural-pin.yml | 17 ++++++++++------- loopstructural/requirements.txt | 2 +- requirements/embedded.txt | 2 +- requirements/testing.txt | 2 +- 4 files changed, 13 insertions(+), 10 deletions(-) diff --git a/.github/workflows/bump-loopstructural-pin.yml b/.github/workflows/bump-loopstructural-pin.yml index 56d135d..2ac31f5 100644 --- a/.github/workflows/bump-loopstructural-pin.yml +++ b/.github/workflows/bump-loopstructural-pin.yml @@ -22,8 +22,9 @@ jobs: # plus loop_common/loop_interpolation sub-packages), so `releases/latest` # can return whichever component's release was published last, not # necessarily the LoopStructural package. Read the manifest instead — - # it maps package name to version unambiguously. - latest=$(curl -fsS https://raw.githubusercontent.com/Loop3D/LoopStructural/master/.release-please-manifest.json | jq -r '.["LoopStructural"] // empty') + # it maps package path to version unambiguously. LoopStructural is the + # root package, so its key is ".". + latest=$(curl -fsS https://raw.githubusercontent.com/Loop3D/LoopStructural/master/.release-please-manifest.json | jq -r '.["."] // empty') if ! [[ "$latest" =~ ^[0-9]+\.[0-9]+\.[0-9]+$ ]]; then echo "::error::Could not read a valid LoopStructural version from the release manifest (got '${latest}')." exit 1 @@ -36,11 +37,13 @@ jobs: latest="${{ steps.latest.outputs.version }}" latest_minor="${latest%.*}" changed=0 - for f in requirements/embedded.txt requirements/testing.txt; do - current=$(grep -oP 'LoopStructural~=\K[0-9]+\.[0-9]+\.[0-9]+' "$f") + # loopstructural/requirements.txt is the file that QGIS (qpip) reads. + # The pin can have an extra, for example LoopStructural[native]~=1.11.0. + for f in loopstructural/requirements.txt requirements/embedded.txt requirements/testing.txt; do + current=$(grep -oP 'LoopStructural(\[[^]]*\])?~=\K[0-9]+\.[0-9]+\.[0-9]+' "$f") current_minor="${current%.*}" if [ "$latest_minor" != "$current_minor" ]; then - sed -i "s/LoopStructural~=[0-9.]\+/LoopStructural~=${latest}/" "$f" + sed -i -E "s/(LoopStructural(\[[^]]*\])?)~=[0-9.]+/\1~=${latest}/" "$f" changed=1 fi done @@ -54,8 +57,8 @@ jobs: title: "Bump LoopStructural pin to ${{ steps.latest.outputs.version }}" body: | LoopStructural released a new minor version (${{ steps.latest.outputs.version }}). - This bumps the `~=` pin in `requirements/embedded.txt` and - `requirements/testing.txt` accordingly. + This bumps the `~=` pin in `loopstructural/requirements.txt`, + `requirements/embedded.txt` and `requirements/testing.txt` accordingly. Verify the `qgis-compat` CI job passes before merging. branch: bump-loopstructural-pin diff --git a/loopstructural/requirements.txt b/loopstructural/requirements.txt index 498d635..cf97248 100644 --- a/loopstructural/requirements.txt +++ b/loopstructural/requirements.txt @@ -1,6 +1,6 @@ pyvistaqt pyvista -LoopStructural~=1.8.0 +LoopStructural[native]~=1.11.0 pyqtgraph loopsolver geopandas diff --git a/requirements/embedded.txt b/requirements/embedded.txt index 1567159..231939e 100644 --- a/requirements/embedded.txt +++ b/requirements/embedded.txt @@ -1,4 +1,4 @@ pyvista>=0.4 #3D visualisation/vtk meshing -LoopStructural~=1.8.0 #3D geological modelling +LoopStructural[native]~=1.11.0 #3D geological modelling geoh5py>=0.9 #saving models to hdf5 dill>=0.4.1 #saving models to python diff --git a/requirements/testing.txt b/requirements/testing.txt index efabef7..5b0ef25 100644 --- a/requirements/testing.txt +++ b/requirements/testing.txt @@ -4,7 +4,7 @@ pytest-cov>=7.1.0 packaging>=23 pyvista>=0.4 #3D visualisation/vtk meshing -LoopStructural~=1.8.0 #3D geological modelling +LoopStructural[native]~=1.11.0 #3D geological modelling geoh5py>=0.9 #saving models to hdf5 dill>=0.4.1 #saving models to python shapely From c814c4686335610193d38ed3763d165f0a1e539e Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Sun, 1 Nov 2026 15:24:51 +1030 Subject: [PATCH 3/4] feat: fall back to the default solver for each coordinate, not each feature The coordinates 1 and 2 of a fault frame have orthogonal rows that the native core cannot read, so all the faults used cg. Now the plugin finds the coordinates that the core refuses (with the support check of loop-interpolation), solves only those with the default solver, and keeps native for the others. It records the coordinates, so the next build does not try native for them again. If the check finds no coordinate, the full feature uses the default solver, as before. The tests for the ADMM solver of inequality rows now turn the native core off, because with it the native solver solves the inequalities. --- loopstructural/main/model_manager.py | 70 +++++++++++++++++++--------- loopstructural/main/solver_choice.py | 29 ++++++++++++ tests/qgis/test_detach_feature.py | 17 ++++++- tests/qgis/test_manual_foliations.py | 14 +++++- tests/qgis/test_native_solver.py | 69 +++++++++++++++++++++++++++ 5 files changed, 175 insertions(+), 24 deletions(-) diff --git a/loopstructural/main/model_manager.py b/loopstructural/main/model_manager.py index cfecfdf..54d043e 100644 --- a/loopstructural/main/model_manager.py +++ b/loopstructural/main/model_manager.py @@ -251,6 +251,9 @@ def __init__(self, debug_manager=None): # features that the native solver could not solve; they use the # default solver (see `_solve_feature`) self._native_failed: Set[str] = set() + # feature name -> indexes of the builders (coordinates) that use the + # default solver because the native solver could not solve them + self._native_fallbacks: Dict[str, Set[int]] = {} # features with constraints that need an inequality solver self._needs_inequalities: Set[str] = set() @@ -331,6 +334,7 @@ def reset(self): self.nelements_used = {} self._data_summaries = {} self._native_failed = set() + self._native_fallbacks = {} self._needs_inequalities = set() self.dem_function = lambda x, y: 0 self._topology_dirty = False @@ -1055,31 +1059,53 @@ def _inequality_solver(self, name, solver): return None if self._use_native(name) else solver def _solve_feature(self, feature, rebuild=True): - """Build `feature`. If the native solver cannot solve it, build it with - the default solver, and use the default solver for it from now on. + """Build `feature`. A coordinate that the native solver cannot solve is + built with the default solver, and uses it from now on. + + Each coordinate of a structural frame has its own builder. For a fault, + the coordinates 1 and 2 have orthogonal rows that the native core cannot + read, but the coordinate 0 can use it. If the check of loop-interpolation + finds no coordinate, the full feature uses the default solver. `rebuild` builds the feature again; otherwise only if it is out of date. """ - try: - if rebuild: - feature.builder.update() - else: - feature.builder.up_to_date() - except solver_choice.NATIVE_ERRORS as error: - if not solver_choice.uses_native(feature): - raise - name = feature.name - self._native_failed.add(name) - fallback = solver_choice.fallback_solver(name in self._needs_inequalities) - if self._debug_manager is not None: - self._debug_manager.log( - f"The native solver cannot solve '{name}' ({error}). " - f"It is solved with '{fallback}'.", - log_level=1, - ) - for builder in solver_choice.builders(feature): - builder.update_build_arguments({'solver': fallback}) - feature.builder.update() + name = feature.name + builders = solver_choice.builders(feature) + fallback = solver_choice.fallback_solver(name in self._needs_inequalities) + # the coordinates that failed in an earlier build + for index in self._native_fallbacks.get(name, ()): + if index < len(builders): + builders[index].update_build_arguments({'solver': fallback}) + build = feature.builder.update if rebuild else feature.builder.up_to_date + # one coordinate (or more) changes to the fallback at each attempt + for _attempt in range(len(builders) + 1): + try: + build() + return + except solver_choice.NATIVE_ERRORS as error: + if not solver_choice.uses_native(feature): + raise + failed = solver_choice.unsupported_builders(builders) + if not failed: + self._native_failed.add(name) + failed = list(range(len(builders))) + self._native_fallbacks.setdefault(name, set()).update(failed) + if self._debug_manager is not None: + which = ( + f"the coordinates {sorted(failed)} of '{name}'" + if len(builders) > 1 + else f"'{name}'" + ) + self._debug_manager.log( + f"The native solver cannot solve {which} ({error}). " + f"They are solved with '{fallback}'.", + log_level=1, + ) + for index in failed: + builders[index].update_build_arguments({'solver': fallback}) + # the coordinates that are built already stay as they are + build = feature.builder.up_to_date + build() def default_nelements(self, name) -> tuple: """Return (number of elements, True if automatic) for the feature diff --git a/loopstructural/main/solver_choice.py b/loopstructural/main/solver_choice.py index 7c22634..25408a9 100644 --- a/loopstructural/main/solver_choice.py +++ b/loopstructural/main/solver_choice.py @@ -55,6 +55,35 @@ def builders(feature): return list(getattr(builder, 'builders', None) or [builder]) +def unsupported_builders(builders): + """Return the indexes of the builders that use the native solver but that + its core cannot solve, for example the coordinates 1 and 2 of a fault + frame, which have orthogonal rows that the core cannot read. + + It uses the support check of loop-interpolation. Returns an empty list if + the check is not available or finds no problem. + """ + try: + from loop_interpolation._native import check_supported + except ImportError: + return [] + failed = [] + for index, builder in enumerate(builders): + if (getattr(builder, 'build_arguments', None) or {}).get('solver') != NATIVE: + continue + interpolator = getattr(builder, 'interpolator', None) + if interpolator is None: + continue + try: + check_supported(interpolator, getattr(interpolator, '_setup_kwargs', {}), {}) + except NotImplementedError: + failed.append(index) + except Exception: + # the check is for one solve; any other problem is not a reason to change + continue + return failed + + def uses_native(feature) -> bool: """Return True if a builder of the feature uses the native solver.""" return any( diff --git a/tests/qgis/test_detach_feature.py b/tests/qgis/test_detach_feature.py index d174609..3248d8f 100644 --- a/tests/qgis/test_detach_feature.py +++ b/tests/qgis/test_detach_feature.py @@ -15,6 +15,7 @@ from loopstructural.main.data_manager import ModellingDataManager from loopstructural.main.model_manager import GeologicalModelManager +from loopstructural.main import solver_choice def _contact(x): @@ -138,7 +139,9 @@ def test_added_rows_stay_when_a_feature_is_detached(self, manager): assert sorted(manager.detached[name]['X']) == [1.0, 2.0] assert _built_x(manager) == [1.0, 2.0, 7.0] - def test_an_inequality_row_gives_the_solver(self, manager): + def test_an_inequality_row_gives_the_solver(self, manager, monkeypatch): + # the default solver; with the native solver, see the next test + monkeypatch.setattr(solver_choice, 'native_available', lambda: False) name = _group_name(manager) spec = self._extra(7.0) spec['data']['values'].update( @@ -148,6 +151,18 @@ def test_an_inequality_row_gives_the_solver(self, manager): manager.update_foliation_features() assert manager._captured[-1][2]['solver'] == 'admm' + def test_the_native_solver_also_solves_an_inequality_row(self, manager, monkeypatch): + monkeypatch.setattr(solver_choice, 'native_available', lambda: True) + monkeypatch.setattr(solver_choice, 'use_native', lambda setting: True) + name = _group_name(manager) + spec = self._extra(7.0) + spec['data']['values'].update( + {'type': 'Inequality', 'lower_field': 'v', 'upper_field': 'v'} + ) + manager.extra_constraints[name] = spec + manager.update_foliation_features() + assert manager._captured[-1][2]['solver'] == 'native' + def test_a_bad_row_is_skipped_and_the_other_rows_are_used(self, manager): name = _group_name(manager) spec = self._extra(7.0) diff --git a/tests/qgis/test_manual_foliations.py b/tests/qgis/test_manual_foliations.py index 8e27583..9343cdd 100644 --- a/tests/qgis/test_manual_foliations.py +++ b/tests/qgis/test_manual_foliations.py @@ -13,6 +13,7 @@ from shapely.geometry import Point from loopstructural.main.model_manager import GeologicalModelManager +from loopstructural.main import solver_choice from loopstructural.toolbelt.preferences import PlgSettingsStructure @@ -169,13 +170,24 @@ def test_a_constant_z_source_gives_one_z(self, manager): data, _ = self._rows(manager, layer) assert (data['Z'] == 12.5).all() - def test_a_pairwise_inequality_needs_the_admm_solver(self, manager): + def test_a_pairwise_inequality_needs_the_admm_solver(self, manager, monkeypatch): + # the default solver; the native solver also solves inequalities + monkeypatch.setattr(solver_choice, 'native_available', lambda: False) df = _points(order=lambda p: 1 if p.x < 50 else 2) layer = {'type': 'Pairwise Inequality', 'pair_field': 'order', 'df': df} data, kwargs = self._rows(manager, layer) assert kwargs == {'solver': 'admm'} assert 'pair_id' in data + def test_with_the_native_solver_a_pairwise_inequality_needs_no_admm( + self, manager, monkeypatch + ): + monkeypatch.setattr(solver_choice, 'use_native', lambda setting: True) + df = _points(order=lambda p: 1 if p.x < 50 else 2) + layer = {'type': 'Pairwise Inequality', 'pair_field': 'order', 'df': df} + _data, kwargs = self._rows(manager, layer) + assert 'solver' not in kwargs + def test_a_feature_from_value_and_normal_layers_solves_without_a_column(self, manager): values = _value_layer() normals = { diff --git a/tests/qgis/test_native_solver.py b/tests/qgis/test_native_solver.py index 4fe479a..2a4ad8a 100644 --- a/tests/qgis/test_native_solver.py +++ b/tests/qgis/test_native_solver.py @@ -106,3 +106,72 @@ def update(self): feature.builder = _Failing('cg') with pytest.raises(NotImplementedError): manager._solve_feature(feature) + + +class _FrameBuilder: + """Three coordinates. A coordinate that the native core cannot solve raises.""" + + def __init__(self, refused): + self.builders = [_Builder(solver_choice.NATIVE) for _ in range(3)] + self.refused = set(refused) + self.built = set() + + def _build(self, everything): + for index, builder in enumerate(self.builders): + if index in self.built and not everything: + continue + solver = builder.build_arguments.get('solver') + builder.solvers.append(solver) + if solver == solver_choice.NATIVE and index in self.refused: + raise NotImplementedError( + "solver='native' does not support the setting " + "'add_constraints_to_least_squares'" + ) + self.built.add(index) + + def update(self): + self.built = set() + self._build(everything=True) + + def up_to_date(self): + self._build(everything=False) + + +class _Fault: + def __init__(self, name, refused=(1, 2)): + self.name = name + self.builder = _FrameBuilder(refused) + + +def test_only_the_coordinates_that_native_cannot_solve_use_the_default(manager, monkeypatch): + _use(monkeypatch, solver_choice.NATIVE) + fault = _Fault('f1') + monkeypatch.setattr(solver_choice, 'unsupported_builders', lambda builders: [1, 2]) + manager._solve_feature(fault) + solvers = [b.build_arguments['solver'] for b in fault.builder.builders] + assert solvers == ['native', 'cg', 'cg'] + # the coordinate 0 is built one time, with native + assert fault.builder.builders[0].solvers == ['native'] + assert manager._native_fallbacks['f1'] == {1, 2} + # the feature still gets native, so the coordinate 0 uses it next time + assert manager.interpolator_arguments('f1', None)['solver'] == 'native' + + +def test_the_next_build_uses_the_recorded_coordinates(manager, monkeypatch): + _use(monkeypatch, solver_choice.NATIVE) + manager._native_fallbacks['f1'] = {1, 2} + fault = _Fault('f1') + manager._solve_feature(fault) + # no native attempt for the coordinates 1 and 2 + assert [b.solvers for b in fault.builder.builders] == [['native'], ['cg'], ['cg']] + + +def test_without_a_coordinate_from_the_check_the_full_feature_uses_the_default( + manager, monkeypatch +): + _use(monkeypatch, solver_choice.NATIVE) + fault = _Fault('f1') + monkeypatch.setattr(solver_choice, 'unsupported_builders', lambda builders: []) + manager._solve_feature(fault) + assert [b.build_arguments['solver'] for b in fault.builder.builders] == ['cg'] * 3 + assert 'f1' in manager._native_failed From 4a4b185b39337fea04e1466e0f54afdbd2791e66 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Sun, 1 Nov 2026 15:30:30 +1030 Subject: [PATCH 4/4] fix: give the native solver the iteration limit of cg, and fall back when it does not converge The native core stops after 5 000 iterations by default; cg allows 10 x the number of nodes. On a large grid the native solve stopped before it converged. loop-interpolation then keeps the old field, which had the size of the old grid, so the surfaces failed ("cannot assign 51000 values to array of shape 87232"). - Before each solve, give each native builder maxiter = 10 x its number of nodes. - After the build, a coordinate whose native solve did not converge is solved again with the default solver, and uses it from then on. The native limit is not given to the default solver. --- loopstructural/main/model_manager.py | 65 ++++++++++++++++++++-------- loopstructural/main/solver_choice.py | 31 +++++++++++++ tests/qgis/test_native_solver.py | 42 ++++++++++++++++++ tests/unit/test_solver_choice.py | 19 ++++++++ 4 files changed, 140 insertions(+), 17 deletions(-) diff --git a/loopstructural/main/model_manager.py b/loopstructural/main/model_manager.py index 54d043e..57a7123 100644 --- a/loopstructural/main/model_manager.py +++ b/loopstructural/main/model_manager.py @@ -1067,6 +1067,10 @@ def _solve_feature(self, feature, rebuild=True): read, but the coordinate 0 can use it. If the check of loop-interpolation finds no coordinate, the full feature uses the default solver. + A native solve that does not converge is also solved again with the + default solver: loop-interpolation then keeps the old field, which can + have the size of an old grid. + `rebuild` builds the feature again; otherwise only if it is out of date. """ name = feature.name @@ -1075,13 +1079,14 @@ def _solve_feature(self, feature, rebuild=True): # the coordinates that failed in an earlier build for index in self._native_fallbacks.get(name, ()): if index < len(builders): - builders[index].update_build_arguments({'solver': fallback}) + self._set_fallback(builders[index], fallback) + self._set_native_limits(builders) build = feature.builder.update if rebuild else feature.builder.up_to_date # one coordinate (or more) changes to the fallback at each attempt for _attempt in range(len(builders) + 1): try: build() - return + break except solver_choice.NATIVE_ERRORS as error: if not solver_choice.uses_native(feature): raise @@ -1089,23 +1094,49 @@ def _solve_feature(self, feature, rebuild=True): if not failed: self._native_failed.add(name) failed = list(range(len(builders))) - self._native_fallbacks.setdefault(name, set()).update(failed) - if self._debug_manager is not None: - which = ( - f"the coordinates {sorted(failed)} of '{name}'" - if len(builders) > 1 - else f"'{name}'" - ) - self._debug_manager.log( - f"The native solver cannot solve {which} ({error}). " - f"They are solved with '{fallback}'.", - log_level=1, - ) - for index in failed: - builders[index].update_build_arguments({'solver': fallback}) + self._use_fallback(name, builders, failed, fallback, f"not supported ({error})") # the coordinates that are built already stay as they are build = feature.builder.up_to_date - build() + else: + build() + unconverged = solver_choice.unconverged_builders(builders) + if unconverged: + self._use_fallback(name, builders, unconverged, fallback, "did not converge") + feature.builder.up_to_date() + + @staticmethod + def _set_fallback(builder, fallback): + # the native iteration limit is not for the other solvers + builder.update_build_arguments({'solver': fallback, 'solver_kwargs': {}}) + + @staticmethod + def _set_native_limits(builders): + """Give each native builder the iteration limit of its grid. + + Only a change of the limit marks the builder out of date. + """ + for builder in builders: + if (builder.build_arguments or {}).get('solver') != solver_choice.NATIVE: + continue + limits = solver_choice.native_solver_kwargs(builder) + if limits and builder.build_arguments.get('solver_kwargs') != limits: + builder.update_build_arguments({'solver_kwargs': limits}) + + def _use_fallback(self, name, builders, indexes, fallback, reason): + """Solve the builders `indexes` of the feature `name` with `fallback` from now on.""" + self._native_fallbacks.setdefault(name, set()).update(indexes) + if self._debug_manager is not None: + which = ( + f"the coordinates {sorted(indexes)} of '{name}'" + if len(builders) > 1 + else f"'{name}'" + ) + self._debug_manager.log( + f"Native solver, {which}: {reason}. Solved with '{fallback}'.", + log_level=1, + ) + for index in indexes: + self._set_fallback(builders[index], fallback) def default_nelements(self, name) -> tuple: """Return (number of elements, True if automatic) for the feature diff --git a/loopstructural/main/solver_choice.py b/loopstructural/main/solver_choice.py index 25408a9..ae0165e 100644 --- a/loopstructural/main/solver_choice.py +++ b/loopstructural/main/solver_choice.py @@ -21,6 +21,8 @@ # The errors that tell that the native solver cannot solve a feature NATIVE_ERRORS = (NotImplementedError, ImportError) +# The iteration limit of a native solve, for each node of the grid (as cg) +NATIVE_ITERATIONS_PER_NODE = 10 def native_available() -> bool: @@ -84,6 +86,35 @@ def unsupported_builders(builders): return failed +def native_solver_kwargs(builder): + """Return the solver_kwargs for a native solve of `builder`, or None. + + The core stops after 5 000 iterations by default. The cg solver of scipy + allows 10 x the number of nodes. With the default, a large grid can stop + before it converges, so give the native solver the same limit as cg. + """ + support = getattr(getattr(builder, 'interpolator', None), 'support', None) + n_nodes = getattr(support, 'n_nodes', None) + if not n_nodes: + return None + return {'maxiter': int(NATIVE_ITERATIONS_PER_NODE * n_nodes)} + + +def unconverged_builders(builders): + """Return the indexes of the builders whose native solve did not converge. + + Then loop-interpolation keeps the old field of the interpolator, which can + have the size of an old grid. + """ + return [ + index + for index, builder in enumerate(builders) + if (getattr(builder, 'build_arguments', None) or {}).get('solver') == NATIVE + and getattr(builder, 'interpolator', None) is not None + and not getattr(builder.interpolator, 'up_to_date', True) + ] + + def uses_native(feature) -> bool: """Return True if a builder of the feature uses the native solver.""" return any( diff --git a/tests/qgis/test_native_solver.py b/tests/qgis/test_native_solver.py index 2a4ad8a..29bf37d 100644 --- a/tests/qgis/test_native_solver.py +++ b/tests/qgis/test_native_solver.py @@ -4,6 +4,8 @@ the fallback, not LoopStructural. """ +from types import SimpleNamespace + import pytest from loopstructural.main import model_manager as model_manager_module @@ -175,3 +177,43 @@ def test_without_a_coordinate_from_the_check_the_full_feature_uses_the_default( manager._solve_feature(fault) assert [b.build_arguments['solver'] for b in fault.builder.builders] == ['cg'] * 3 assert 'f1' in manager._native_failed + + +class _GridBuilder(_Builder): + """A builder with an interpolator. The native solve converges if `converges`.""" + + def __init__(self, converges=True, n_nodes=1000): + super().__init__(solver_choice.NATIVE) + self.converges = converges + self.interpolator = SimpleNamespace( + support=SimpleNamespace(n_nodes=n_nodes), up_to_date=False + ) + self.kwargs_at_solve = [] + + def update(self): + solver = self.build_arguments.get('solver') + self.solvers.append(solver) + self.kwargs_at_solve.append(dict(self.build_arguments.get('solver_kwargs') or {})) + self.interpolator.up_to_date = solver != solver_choice.NATIVE or self.converges + + up_to_date = update + + +def test_the_native_solve_gets_the_iteration_limit_of_its_grid(manager, monkeypatch): + _use(monkeypatch, solver_choice.NATIVE) + feature = _Feature('s1') + feature.builder = _GridBuilder(n_nodes=87232) + manager._solve_feature(feature) + assert feature.builder.kwargs_at_solve == [{'maxiter': 872320}] + + +def test_a_native_solve_that_does_not_converge_is_solved_with_the_default(manager, monkeypatch): + _use(monkeypatch, solver_choice.NATIVE) + feature = _Feature('Group_0') + feature.builder = _GridBuilder(converges=False) + manager._solve_feature(feature) + assert feature.builder.solvers == ['native', 'cg'] + # the native limit is not given to cg + assert feature.builder.kwargs_at_solve[-1] == {} + assert manager._native_fallbacks['Group_0'] == {0} + assert any('did not converge' in m for m in manager._debug_manager.messages) diff --git a/tests/unit/test_solver_choice.py b/tests/unit/test_solver_choice.py index f43956d..040fee1 100644 --- a/tests/unit/test_solver_choice.py +++ b/tests/unit/test_solver_choice.py @@ -49,3 +49,22 @@ def test_uses_native(): def test_every_choice_has_a_label(): assert set(solver_choice.LABELS) == set(solver_choice.CHOICES) + + +def _native_builder(n_nodes=1000, up_to_date=True, solver='native'): + interpolator = SimpleNamespace(support=SimpleNamespace(n_nodes=n_nodes), up_to_date=up_to_date) + return SimpleNamespace(build_arguments={'solver': solver}, interpolator=interpolator) + + +def test_the_native_limit_is_ten_iterations_for_each_node(): + assert solver_choice.native_solver_kwargs(_native_builder(87232)) == {'maxiter': 872320} + assert solver_choice.native_solver_kwargs(SimpleNamespace()) is None + + +def test_unconverged_builders(): + builders = [ + _native_builder(), + _native_builder(up_to_date=False), + _native_builder(up_to_date=False, solver='cg'), + ] + assert solver_choice.unconverged_builders(builders) == [1]