diff --git a/packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py b/packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py index a51c1a3e..4fd07351 100644 --- a/packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py +++ b/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 diff --git a/packages/loop_interpolation/src/loop_interpolation/_finite_difference_interpolator.py b/packages/loop_interpolation/src/loop_interpolation/_finite_difference_interpolator.py index cedad674..73a949a3 100644 --- a/packages/loop_interpolation/src/loop_interpolation/_finite_difference_interpolator.py +++ b/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( diff --git a/packages/loop_interpolation/src/loop_interpolation/_native.py b/packages/loop_interpolation/src/loop_interpolation/_native.py index 02c77ed8..4e2ac2f6 100644 --- a/packages/loop_interpolation/src/loop_interpolation/_native.py +++ b/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( diff --git a/packages/loop_interpolation/tests/test_native_orthogonal_and_limits.py b/packages/loop_interpolation/tests/test_native_orthogonal_and_limits.py new file mode 100644 index 00000000..4871b8a6 --- /dev/null +++ b/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) diff --git a/tests/unit/modelling/test_native_fault_frame.py b/tests/unit/modelling/test_native_fault_frame.py new file mode 100644 index 00000000..42ebfbe3 --- /dev/null +++ b/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