feat(native): solve structural frames, use the cg iteration limit, and keep no stale field

packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py

@@ -327,6 +327,9 @@ def reset(self):
 
         """
         self.constraints = {}
+        # Gradient orthogonal constraints of a native setup, as raw data:
+        # name -> (points, vectors, w). See add_gradient_orthogonal_constraints.
+        self._native_orthogonal = {}
         # setup_interpolator sets these two values again when it ends.
         self._native_setup = False
         self._setup_constraint_names = None

packages/loop_interpolation/src/loop_interpolation/_finite_difference_interpolator.py

@@ -518,6 +518,19 @@ def add_gradient_orthogonal_constraints(
         -------
 
         """
+        if getattr(self, "_native_setup", False) and np.all(np.asarray(b) == 0):
+            # The native core reads raw data, not sparse rows. With b = 0 this
+            # constraint is a tangent constraint (gradient . vector = 0), so keep
+            # the points and the vectors for the core (see _native._add_data).
+            # A structural frame adds these for the orthogonal coordinates.
+            if points.shape[0] > 0:
+                self._native_orthogonal[name] = (
+                    np.array(points[:, :3], dtype=float),
+                    np.array(vectors[:, :3], dtype=float),
+                    w,
+                )
+                self.up_to_date = False
+            return
         if points.shape[0] > 0:
             # calculate unit vector for orientation data
             node_idx, inside = self.support.position_to_cell_corners(

packages/loop_interpolation/src/loop_interpolation/_native.py

@@ -94,6 +94,9 @@
 _PAIR_SEED = 0
 _PAIR_CHUNK = 100_000
 
+# Iteration limit for each node when solver_kwargs has no maxiter (as cg)
+MAXITER_PER_NODE = 10
+
 _REGULARISATION_KEYS = ("dxx", "dyy", "dzz", "dxy", "dxz", "dyz")
 _BORDER_KEYS = ("dx", "dy", "dz")
 
@@ -207,9 +210,16 @@ def check_supported(interpolator, setup_kwargs, solver_kwargs) -> None:
         logger.info("solver='native': column scaling (apply_scaling_matrix) is skipped.")
 
 
-def _core_settings(tol, solver_kwargs) -> dict:
-    """Make the settings dict of the core from ``tol`` and ``solver_kwargs``."""
+def _core_settings(tol, solver_kwargs, n_nodes=None) -> dict:
+    """Make the settings dict of the core from ``tol`` and ``solver_kwargs``.
+
+    Without ``maxiter``, the limit is ``MAXITER_PER_NODE`` x ``n_nodes``, as the
+    ``cg`` matrix solver (scipy: 10 x the number of unknowns). The default of
+    the core (5000) stops too early on a large grid.
+    """
     settings = {}
+    if n_nodes and solver_kwargs.get("maxiter") is None:
+        settings["max_iter"] = int(MAXITER_PER_NODE * n_nodes)
     if tol is not None:
         settings["tol"] = float(tol)
     for key, value in solver_kwargs.items():
@@ -313,6 +323,21 @@ def _add_data(problem, interpolator) -> None:
             name="tangent",
         )
 
+    # gradient orthogonal constraints of a native setup (for example the
+    # orthogonal coordinates of a structural frame): the same rows as tangents
+    orthogonal = getattr(interpolator, "_native_orthogonal", None)
+    if not isinstance(orthogonal, dict):
+        orthogonal = {}
+    for name, (points, vectors, w) in orthogonal.items():
+        _add_weighted(
+            problem.add_tangent,
+            points,
+            vectors,
+            _weights(w, points.shape[0]),
+            drop_zero_norm=True,
+            name=name,
+        )
+
     # interface: one weight (ipw)
     interface = _rows(interpolator.get_interface_constraints())
     if interface.shape[0] > 1:
@@ -486,16 +511,17 @@ def solve(interpolator, tol, solver_kwargs, x0=None, ridge_factor=None) -> bool:
 
     Notes
     -----
-    ``interpolator.c`` changes only when the solve converged. When the solve
-    does not converge, ``c`` keeps its old value. The partial field is in
-    ``interpolator.latest_native_report["field"]``. A converged report has no
-    ``field`` key. The matrix solvers (for example ``cg``) differ here. They
+    ``interpolator.c`` changes when the solve converged. When the solve does
+    not converge, ``c`` keeps its old value, if it is a field of this grid.
+    The partial field is in ``interpolator.latest_native_report["field"]``.
+    A converged report has no ``field`` key. If the old ``c`` is missing or
+    has another size (the grid changed since the last solve), ``c`` gets the
+    partial field, as with the matrix solvers (for example ``cg``), which
     write the last iterate to ``c`` also when they do not converge.
 
-    The defaults of the core are ``tol=1e-8`` and ``maxiter=5000``. The
-    ``cg`` matrix solver uses ``rtol=1e-8``. Its ``maxiter`` is the scipy
-    default, 10 times the number of nodes. So the native solver can stop
-    earlier on large grids. Set ``maxiter`` to get the same limit.
+    The default of the core is ``tol=1e-8``. The ``cg`` matrix solver uses
+    ``rtol=1e-8``. Without ``maxiter`` in ``solver_kwargs``, the limit is 10
+    times the number of nodes, the scipy default of ``cg``.
     """
     started = perf_counter()
     solver_kwargs = dict(solver_kwargs or {})
@@ -524,7 +550,7 @@ def solve(interpolator, tol, solver_kwargs, x0=None, ridge_factor=None) -> bool:
         problem.set_start(np.ascontiguousarray(x0, dtype=np.float64))
 
     # tol: the solver_kwargs value wins over the tol argument.
-    settings = _core_settings(solver_kwargs.get("tol", tol), solver_kwargs)
+    settings = _core_settings(solver_kwargs.get("tol", tol), solver_kwargs, support.n_nodes)
     field, report = problem.solve(settings)
 
     # The core can return float32. ``c`` is always float64.
@@ -541,6 +567,11 @@ def solve(interpolator, tol, solver_kwargs, x0=None, ridge_factor=None) -> bool:
     else:
         # Keep the old ``c``. Give the partial field in the report.
         report["field"] = field
+        old = getattr(interpolator, "c", None)
+        if old is None or np.shape(old) != (support.n_nodes,) or not np.any(np.isfinite(old)):
+            # There is no old field of this grid (a first solve, or the grid
+            # changed since the last solve), so a caller cannot use the old c.
+            interpolator.c = field
     interpolator.latest_native_report = report
     if not converged:
         logger.warning(

packages/loop_interpolation/tests/test_native_orthogonal_and_limits.py

@@ -0,0 +1,92 @@
+"""Gradient orthogonal constraints, the iteration limit and a solve that does
+not converge, with solver="native".
+
+All tests skip while the native core (loop_interpolation_core) is not installed.
+"""
+
+import numpy as np
+import pytest
+
+# Skip the whole file when the native core is not installed.
+pytest.importorskip("loop_interpolation_core")
+
+from loop_interpolation import FiniteDifferenceInterpolator as FDI
+from loop_interpolation import StructuredGrid, _native
+
+FIELD_TOL = 1e-4  # fraction of the field range
+
+
+def _make_interpolator(solver, b=0.0):
+    """Make a 15^3 grid with value data and a gradient orthogonal constraint."""
+    origin = np.array([-0.1, -0.1, -0.1])
+    maximum = np.array([1.1, 1.1, 1.1])
+    nsteps = np.array([15, 15, 15])
+    grid = StructuredGrid(origin=origin, nsteps=nsteps, step_vector=(maximum - origin) / nsteps)
+    interpolator = FDI(grid)
+    interpolator.solver = solver
+    rng = np.random.default_rng(0)
+    xyz = rng.uniform(0.1, 0.9, size=(60, 3))
+    values = xyz[:, 2] + 0.3 * xyz[:, 0]
+    interpolator.set_value_constraints(np.column_stack([xyz, values, np.ones(60)]))
+    interpolator.setup()
+    # as a structural frame: the gradient is orthogonal to x at the barycentres
+    points = grid.barycentre[::3]
+    vectors = np.tile([1.0, 0.0, 0.0], (points.shape[0], 1))
+    interpolator.add_gradient_orthogonal_constraints(
+        points, vectors, w=1.0, b=b, name="frame_orthogonal"
+    )
+    return interpolator
+
+
+def _range(c):
+    return float(np.max(c) - np.min(c))
+
+
+def test_a_native_setup_keeps_orthogonal_constraints_as_raw_data():
+    interpolator = _make_interpolator("native")
+    assert "frame_orthogonal" in interpolator._native_orthogonal
+    assert "frame_orthogonal" not in interpolator.constraints
+    # the core can solve it: no NotImplementedError
+    _native.check_supported(interpolator, interpolator._setup_kwargs, {})
+
+
+def test_orthogonal_constraints_give_the_same_field_as_cg():
+    native = _make_interpolator("native")
+    assert native.solve_system("native")
+    matrix = _make_interpolator("cg")
+    assert "frame_orthogonal" in matrix.constraints
+    assert matrix.solve_system("cg", solver_kwargs={"rtol": 1e-12, "atol": 0.0})
+    assert np.max(np.abs(native.c - matrix.c)) <= FIELD_TOL * _range(matrix.c)
+    # the constraint has an effect: the gradient in x is small
+    gradient = native.support.evaluate_gradient(native.support.barycentre[::3], native.c)
+    assert np.median(np.abs(gradient[:, 0])) < 0.1 * np.median(np.abs(gradient[:, 2]))
+
+
+def test_an_orthogonal_constraint_with_b_is_not_supported():
+    interpolator = _make_interpolator("native", b=0.5)
+    assert "frame_orthogonal" in interpolator.constraints
+    with pytest.raises(NotImplementedError, match="add_constraints_to_least_squares"):
+        _native.check_supported(interpolator, interpolator._setup_kwargs, {})
+
+
+def test_the_default_iteration_limit_follows_the_grid():
+    assert _native._core_settings(None, {}, 1000)["max_iter"] == 10 * 1000
+    assert _native._core_settings(None, {"maxiter": 7}, 1000)["max_iter"] == 7
+    assert "max_iter" not in _native._core_settings(None, {}, None)
+
+
+def test_a_solve_that_does_not_converge_replaces_a_field_of_another_grid():
+    interpolator = _make_interpolator("native")
+    # the field of an old, smaller grid
+    interpolator.c = np.zeros(10)
+    assert not interpolator.solve_system("native", solver_kwargs={"maxiter": 1})
+    assert interpolator.c.shape == (interpolator.support.n_nodes,)
+    assert "field" in interpolator.latest_native_report
+
+
+def test_a_solve_that_does_not_converge_keeps_a_field_of_this_grid():
+    interpolator = _make_interpolator("native")
+    assert interpolator.solve_system("native")
+    good = interpolator.c.copy()
+    assert not interpolator.solve_system("native", solver_kwargs={"maxiter": 1})
+    assert np.array_equal(interpolator.c, good)

tests/unit/modelling/test_native_fault_frame.py

@@ -0,0 +1,53 @@
+"""A fault frame with solver="native": all three coordinates use the native core.
+
+The coordinates 1 and 2 have gradient orthogonal constraints. The native core
+gets them as raw data, so it solves the full frame.
+The tests skip while the native core (loop_interpolation_core) is not installed.
+"""
+
+import numpy as np
+import pandas as pd
+import pytest
+
+pytest.importorskip("loop_interpolation_core")
+
+from LoopStructural import GeologicalModel
+
+
+def _fault_model(solver):
+    data = pd.DataFrame(
+        {
+            "X": [500.0, 500.0, 500.0, 500.0],
+            "Y": [100.0, 500.0, 900.0, 500.0],
+            "Z": [-500.0, -500.0, -500.0, -500.0],
+            "val": [0.0, 0.0, 0.0, np.nan],
+            "gx": [np.nan, np.nan, np.nan, 1.0],
+            "gy": [np.nan, np.nan, np.nan, 0.0],
+            "gz": [np.nan, np.nan, np.nan, 0.0],
+            "feature_name": "f1",
+        }
+    )
+    model = GeologicalModel([0, 0, -1000], [1000, 1000, 0])
+    model.data = data
+    fault = model.create_and_add_fault("f1", 50, nelements=20_000, solver=solver)
+    fault.builder.up_to_date()
+    return fault
+
+
+def test_all_coordinates_of_a_fault_use_the_native_solver():
+    fault = _fault_model("native")
+    for coordinate in range(3):
+        builder = fault[coordinate].builder
+        assert builder.build_arguments["solver"] == "native"
+        assert builder.interpolator.latest_native_report["converged"]
+
+
+def test_the_native_fault_frame_is_the_cg_fault_frame():
+    native = _fault_model("native")
+    matrix = _fault_model("cg")
+    points = np.random.default_rng(0).random((200, 3)) * [1000, 1000, 1000] + [0, 0, -1000]
+    for coordinate in range(3):
+        a = native[coordinate].evaluate_value(points)
+        b = matrix[coordinate].evaluate_value(points)
+        scale = np.nanmax(np.abs(b))
+        assert np.nanmax(np.abs(a - b)) <= 1e-3 * scale, coordinate