perf(prediction): optional SciPy-free evaluation of the self-consumption tables (#1378)

* perf(prediction): evaluate self-consumption tables without SciPy

The GENETIC inverter (and the GENETIC0 one) evaluates a pickled
scipy.interpolate.RegularGridInterpolator for the self-consumption
probabilities. SciPy is used nowhere else in EOS, yet loading it costs
about 35 MB RSS (21 MB anonymous) on Linux arm64, which every
optimization run pays.

Add utils/gridinterpolator.py:

- LinearGridInterpolator: 2-D linear interpolation with the same steps
  as SciPy's linear 2-D path (interval search, normalised distances,
  corner weights summed in the same order, fill value outside the grid,
  NaN for NaN input).
- load_grid_interpolator(): reads the shipped pickles with a restricted
  unpickler that maps the SciPy class to a stand-in and accepts only
  the NumPy globals needed for the arrays. SciPy is not imported, and
  the data files no longer depend on SciPy's private module paths.

Results are bit-identical to SciPy on Linux (aarch64, checked on 3.6M
points per table including the query pattern of the inverter). The
macOS arm64 SciPy wheel fuses multiply-add, so there SciPy itself
differs by <= 1.1e-16; GENETIC plans for fixed seeds are unchanged on
both platforms.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

* feat(prediction): make the SciPy-free table evaluation opt-in

Add optimization.self_consumption_interpolator (scipy | numpy, default scipy).
The default keeps evaluating the self-consumption tables with SciPy's
RegularGridInterpolator, built from the same grid, values and fill value the
pickle holds, so results are unchanged. With numpy the tables are evaluated by
LinearGridInterpolator and SciPy is never imported. The backend is read on
every evaluation, so a config change applies to the next run.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

---------

Co-authored-by: Christin <christin@BigBrother.bikini.bottom.zone>
Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
Co-authored-by: Andreas <drbacke@gmx.de>

CHANGELOG.md

@@ -65,6 +65,10 @@ compatibility changes and a short test procedure. These changes are not yet a ta
   values. Runtime forecast energy is Wh per slot and prices are currency per Wh.
 - Runtime configuration changes now take priority over environment and file values for
   explicitly updated keys; command-line settings retain higher priority.
+- New setting `optimization.self_consumption_interpolator` (`scipy` by default, `numpy`).
+  With `numpy` the self-consumption probability tables are evaluated without SciPy, so a
+  GENETIC run no longer loads SciPy (about 35 MB less memory). Results are bit-identical to
+  SciPy's Linux builds. The tables are read without SciPy in both cases.
 - Lower idle memory: matplotlib (PDF reports), statsmodels (price fallback forecast), pvlib,
   bs4/lxml, scipy (self-consumption table), the GENETIC0 optimizer and the timezone lookup
   data are imported on first use instead of at server start. Results are unchanged.

docs/_generated/configexample.md

@@ -276,7 +276,8 @@
                "penalties": {
                    "ev_soc_miss": 10
                }
-           }
+           },
+           "self_consumption_interpolator": "scipy"
        },
        "prediction": {
            "hours": 48,

docs/_generated/configoptimization.md

@@ -12,6 +12,7 @@
 | genetic | `EOS_OPTIMIZATION__GENETIC` | `GeneticCommonSettings` | `rw` | `required` | GENETIC optimization algorithm configuration. |
 | genetic0 | `EOS_OPTIMIZATION__GENETIC0` | `Genetic0CommonSettings` | `rw` | `required` | GENETIC0 optimization algorithm configuration. |
 | keys | | `list[str]` | `ro` | `N/A` | The keys of the solution. |
+| self_consumption_interpolator | `EOS_OPTIMIZATION__SELF_CONSUMPTION_INTERPOLATOR` | `<enum 'GridInterpolatorBackend'>` | `rw` | `scipy` | How the inverter's self-consumption probability table is evaluated [scipy | numpy]. 'scipy' (default) uses scipy.interpolate. 'numpy' gives the same values without importing SciPy, which saves about 35 MB RSS after the first optimization run (useful on devices with little RAM); a run takes about 2.5 % longer. |
 :::
 <!-- pyml enable line-length -->
 
@@ -48,7 +49,8 @@
                "penalties": {
                    "ev_soc_miss": 10
                }
-           }
+           },
+           "self_consumption_interpolator": "scipy"
        }
    }
 ```
@@ -91,6 +93,7 @@
                "interval_sec": 3600,
                "horizon": 24
            },
+           "self_consumption_interpolator": "scipy",
            "algorithms": [
                "GENETIC",
                "GENETIC0"

docs/_generated/openapi.md

@@ -1,6 +1,6 @@
 # Akkudoktor-EOS
 
-**Version**: `v0.3.0.dev2609171762020752`
+**Version**: `v0.3.0.dev2610021392659061`
 
 <!-- pyml disable line-length -->
 **Description**: This project provides a comprehensive solution for simulating and optimizing an energy system based on renewable energy sources. With a focus on photovoltaic (PV) systems, battery storage (batteries), load management (consumer requirements), heat pumps, electric vehicles, and consideration of electricity price data, this system enables forecasting and optimization of energy flow and costs over a specified period.

openapi.json

@@ -8,7 +8,7 @@
       "name": "Apache 2.0",
       "url": "https://www.apache.org/licenses/LICENSE-2.0.html"
     },
-    "version": "v0.3.0.dev2609171762020752"
+    "version": "v0.3.0.dev2610021392659061"
   },
   "paths": {
     "/v1/measurement/battery-capacity/{battery_id}": {
@@ -9970,6 +9970,15 @@
         ],
         "title": "GesamtlastRequest"
       },
+      "GridInterpolatorBackend": {
+        "type": "string",
+        "enum": [
+          "scipy",
+          "numpy"
+        ],
+        "title": "GridInterpolatorBackend",
+        "description": "How the self-consumption tables are evaluated."
+      },
       "HTTPValidationError": {
         "properties": {
           "detail": {
@@ -13091,6 +13100,15 @@
                 }
               }
             ]
+          },
+          "self_consumption_interpolator": {
+            "$ref": "#/components/schemas/GridInterpolatorBackend",
+            "description": "How the inverter's self-consumption probability table is evaluated [scipy | numpy]. 'scipy' (default) uses scipy.interpolate. 'numpy' gives the same values without importing SciPy, which saves about 35 MB RSS after the first optimization run (useful on devices with little RAM); a run takes about 2.5 % longer.",
+            "default": "scipy",
+            "examples": [
+              "scipy",
+              "numpy"
+            ]
           }
         },
         "type": "object",
@@ -13131,6 +13149,15 @@
               }
             ]
           },
+          "self_consumption_interpolator": {
+            "$ref": "#/components/schemas/GridInterpolatorBackend",
+            "description": "How the inverter's self-consumption probability table is evaluated [scipy | numpy]. 'scipy' (default) uses scipy.interpolate. 'numpy' gives the same values without importing SciPy, which saves about 35 MB RSS after the first optimization run (useful on devices with little RAM); a run takes about 2.5 % longer.",
+            "default": "scipy",
+            "examples": [
+              "scipy",
+              "numpy"
+            ]
+          },
           "algorithms": {
             "items": {
               "type": "string"

src/akkudoktoreos/optimization/genetic0/genetic0loadinterpolator.py

@@ -1,25 +1,28 @@
 #!/usr/bin/env python
-import pickle
 import threading
 from pathlib import Path
-from typing import TYPE_CHECKING, Optional
+from typing import Optional
 
 import numpy as np
 
-if TYPE_CHECKING:
-    # Unpickling the table imports scipy when it is first loaded.
-    from scipy.interpolate import RegularGridInterpolator
-
 from akkudoktoreos.core.cache import cache_energy_management
 from akkudoktoreos.core.coreabc import SingletonMixin
+from akkudoktoreos.utils.gridinterpolator import (
+    SwitchableGridInterpolator,
+    configured_grid_interpolator_backend,
+    load_grid_interpolator,
+)
 
 
 class SelfConsumptionProbabilityInterpolator:
     def __init__(self, filepath: str | Path):
         self.filepath = filepath
-        # Load the RegularGridInterpolator
-        with open(self.filepath, "rb") as file:
-            self.interpolator: RegularGridInterpolator = pickle.load(file)  # noqa: S301
+        # The table is a pickled SciPy RegularGridInterpolator. It is read
+        # without SciPy; optimization.self_consumption_interpolator decides
+        # whether SciPy (default) or NumPy evaluates it.
+        self.interpolator = SwitchableGridInterpolator(
+            load_grid_interpolator(self.filepath), configured_grid_interpolator_backend
+        )
 
     def _generate_points(
         self, load_1h_power: float, pv_power: float
@@ -97,7 +100,7 @@ class Genetic0LoadInterpolator(SelfConsumptionProbabilityInterpolator, Singleton
         super().__init__(filename)
 
 
-# Created on first use: loading the table imports scipy.
+# Created on first use.
 _genetic0_load_interpolator: Optional[Genetic0LoadInterpolator] = None
 _genetic0_load_interpolator_lock = threading.Lock()
 

src/akkudoktoreos/optimization/optimization.py

@@ -15,6 +15,7 @@ from akkudoktoreos.optimization.genetic.geneticsettings import (
     normalize_genetic_settings,
 )
 from akkudoktoreos.utils.datetimeutil import DateTime
+from akkudoktoreos.utils.gridinterpolator import GridInterpolatorBackend
 
 
 class OptimizationAlgorithm(StrEnum):
@@ -66,6 +67,20 @@ class OptimizationCommonSettings(SettingsBaseModel):
         },
     )
 
+    self_consumption_interpolator: GridInterpolatorBackend = Field(
+        default=GridInterpolatorBackend.SCIPY,
+        json_schema_extra={
+            "description": (
+                "How the inverter's self-consumption probability table is evaluated "
+                "[scipy | numpy]. 'scipy' (default) uses scipy.interpolate. 'numpy' gives "
+                "the same values without importing SciPy, which saves about 35 MB RSS "
+                "after the first optimization run (useful on devices with little RAM); "
+                "a run takes about 2.5 % longer."
+            ),
+            "examples": ["scipy", "numpy"],
+        },
+    )
+
     # Computed fields
     @computed_field  # type: ignore[prop-decorator]
     @property

src/akkudoktoreos/prediction/interpolator.py

@@ -1,25 +1,28 @@
 #!/usr/bin/env python
-import pickle
 import threading
 from pathlib import Path
-from typing import TYPE_CHECKING, Optional
+from typing import Optional
 
 import numpy as np
 
-if TYPE_CHECKING:
-    # Unpickling the table imports scipy when it is first loaded.
-    from scipy.interpolate import RegularGridInterpolator
-
 from akkudoktoreos.core.cache import cache_energy_management
 from akkudoktoreos.core.coreabc import SingletonMixin
+from akkudoktoreos.utils.gridinterpolator import (
+    SwitchableGridInterpolator,
+    configured_grid_interpolator_backend,
+    load_grid_interpolator,
+)
 
 
 class SelfConsumptionProbabilityInterpolator:
     def __init__(self, filepath: str | Path):
         self.filepath = filepath
-        # Load the RegularGridInterpolator
-        with open(self.filepath, "rb") as file:
-            self.interpolator: RegularGridInterpolator = pickle.load(file)  # noqa: S301
+        # The table is a pickled SciPy RegularGridInterpolator. It is read
+        # without SciPy; optimization.self_consumption_interpolator decides
+        # whether SciPy (default) or NumPy evaluates it.
+        self.interpolator = SwitchableGridInterpolator(
+            load_grid_interpolator(self.filepath), configured_grid_interpolator_backend
+        )
         self.load_power_min_w = float(self.interpolator.grid[0][0])
         self.load_power_max_w = float(self.interpolator.grid[0][-1])
         self.minute_load_levels_w = np.asarray(self.interpolator.grid[1], dtype=float)
@@ -175,7 +178,7 @@ class EOSLoadInterpolator(SelfConsumptionProbabilityInterpolator, SingletonMixin
         super().__init__(filename)
 
 
-# Created on first use: loading the table imports scipy.
+# Created on first use.
 _eos_load_interpolator: Optional[EOSLoadInterpolator] = None
 _eos_load_interpolator_lock = threading.Lock()
 

src/akkudoktoreos/utils/gridinterpolator.py

@@ -0,0 +1,247 @@
+"""Linear interpolation on a 2-D rectilinear grid without SciPy.
+
+The self-consumption probability tables shipped in ``akkudoktoreos/data`` are
+pickled ``scipy.interpolate.RegularGridInterpolator`` objects. Importing SciPy
+only to evaluate them costs about 40 MB of resident memory, which matters on
+small devices. ``load_grid_interpolator()`` reads those pickles without
+importing SciPy and returns a ``LinearGridInterpolator`` that evaluates
+the same piecewise bilinear function.
+
+The arithmetic follows SciPy's ``method="linear"`` 2-D path step by step
+(interval search, normalised distances, corner weights summed in the same
+order), so the results are bit-identical to SciPy builds that do not fuse
+multiply-add operations (e.g. the Linux wheels). Builds that do fuse them
+(e.g. the macOS arm64 wheel) differ only by rounding (about 1e-16).
+
+``SwitchableGridInterpolator`` chooses the evaluation per call: SciPy (the
+default, unchanged behaviour) or the NumPy implementation. SciPy is imported
+only when it is actually used.
+"""
+
+import pickle
+import threading
+from enum import StrEnum
+from pathlib import Path
+from typing import Any, BinaryIO, Callable, Optional, Union
+
+import numpy as np
+
+
+class LinearGridInterpolator:
+    """Piecewise bilinear interpolation on a 2-D rectilinear grid.
+
+    Behaves like ``scipy.interpolate.RegularGridInterpolator(grid, values,
+    method="linear", bounds_error=False, fill_value=fill_value)``: points
+    outside the grid get ``fill_value``, points with a NaN coordinate get NaN.
+
+    Args:
+        grid: Two strictly ascending 1-D coordinate arrays, each with at least two points.
+        values: Values on the grid, shape ``(len(grid[0]), len(grid[1]))``.
+        fill_value: Value for points outside the grid. ``None`` extrapolates linearly.
+    """
+
+    def __init__(
+        self,
+        grid: tuple[Any, Any],
+        values: Any,
+        fill_value: Union[float, None] = np.nan,
+    ) -> None:
+        if len(grid) != 2:
+            raise ValueError(f"Expected a 2-D grid, got {len(grid)} dimensions.")
+        self._grid = tuple(np.asarray(axis, dtype=float) for axis in grid)
+        self._values = np.asarray(values, dtype=float)
+        for axis in self._grid:
+            if axis.ndim != 1 or axis.shape[0] < 2:
+                raise ValueError("Every grid axis needs at least two points.")
+            if not np.all(np.diff(axis) > 0):
+                raise ValueError("Grid axes must be strictly ascending.")
+        if self._values.shape != (self._grid[0].shape[0], self._grid[1].shape[0]):
+            raise ValueError(
+                f"Values shape {self._values.shape} does not match the grid "
+                f"({self._grid[0].shape[0]}, {self._grid[1].shape[0]})."
+            )
+        self.fill_value = fill_value
+
+    @property
+    def grid(self) -> tuple[np.ndarray, ...]:
+        """Grid coordinates per axis."""
+        return self._grid
+
+    @property
+    def values(self) -> np.ndarray:
+        """Values on the grid."""
+        return self._values
+
+    def __call__(self, xi: Any) -> np.ndarray:
+        """Interpolate at points ``xi`` of shape ``(..., 2)``."""
+        xi = np.asarray(xi, dtype=float)
+        if xi.shape[-1] != 2:
+            raise ValueError(f"Points must have 2 coordinates, got {xi.shape[-1]}.")
+        result_shape = xi.shape[:-1]
+        points = xi.reshape(-1, 2)
+
+        indices = []
+        distances = []
+        out_of_bounds = np.zeros(points.shape[0], dtype=bool)
+        for axis, coordinates in zip(self._grid, points.T):
+            # Interval i with axis[i] <= x < axis[i + 1]; the last interval also
+            # takes x == axis[-1] and points outside extrapolate from the edges.
+            index = np.searchsorted(axis, coordinates, side="right") - 1
+            index = np.clip(index, 0, axis.shape[0] - 2)
+            indices.append(index)
+            distances.append((coordinates - axis[index]) / (axis[index + 1] - axis[index]))
+            out_of_bounds |= coordinates < axis[0]
+            out_of_bounds |= coordinates > axis[-1]
+
+        i0, i1 = indices
+        y0, y1 = distances
+        values = self._values
+        result = values[i0, i1] * (1 - y0) * (1 - y1)
+        result = result + values[i0, i1 + 1] * (1 - y0) * y1
+        result = result + values[i0 + 1, i1] * y0 * (1 - y1)
+        result = result + values[i0 + 1, i1 + 1] * y0 * y1
+
+        if self.fill_value is not None:
+            result[out_of_bounds] = self.fill_value
+        result[np.any(np.isnan(points), axis=-1)] = np.nan
+        return result.reshape(result_shape)
+
+
+class _PickledRegularGridInterpolator:
+    """Stand-in for ``scipy.interpolate.RegularGridInterpolator`` while unpickling."""
+
+    def __setstate__(self, state: dict[str, Any]) -> None:
+        self.state = state
+
+
+class _GridInterpolatorUnpickler(pickle.Unpickler):
+    """Unpickler that only accepts a pickled SciPy RegularGridInterpolator.
+
+    SciPy classes are replaced by stand-ins, so SciPy is never imported; any
+    other global is rejected.
+    """
+
+    _NUMPY_GLOBALS = {
+        ("numpy._core.multiarray", "_reconstruct"),
+        ("numpy.core.multiarray", "_reconstruct"),
+        ("numpy", "ndarray"),
+        ("numpy", "dtype"),
+    }
+
+    def find_class(self, module: str, name: str) -> Any:
+        if module.startswith("scipy.interpolate") and name == "RegularGridInterpolator":
+            return _PickledRegularGridInterpolator
+        if module.startswith("scipy.") and name == "asarray":
+            # Array conversion helper stored on the interpolator; not needed here.
+            return np.asarray
+        if (module, name) in self._NUMPY_GLOBALS:
+            return super().find_class(module, name)
+        raise pickle.UnpicklingError(f"Unexpected object {module}.{name} in grid interpolator.")
+
+
+def load_grid_interpolator(source: Union[str, Path, BinaryIO]) -> LinearGridInterpolator:
+    """Load a pickled 2-D linear ``RegularGridInterpolator`` without SciPy.
+
+    Args:
+        source: Path to the pickle file or an open binary file.
+
+    Returns:
+        The interpolator as ``LinearGridInterpolator``.
+
+    Raises:
+        pickle.UnpicklingError: If the file holds anything else than a 2-D linear
+            ``RegularGridInterpolator``.
+    """
+    if isinstance(source, (str, Path)):
+        with open(source, "rb") as file:
+            pickled = _GridInterpolatorUnpickler(file).load()
+    else:
+        pickled = _GridInterpolatorUnpickler(source).load()
+    if not isinstance(pickled, _PickledRegularGridInterpolator):
+        raise pickle.UnpicklingError("File does not hold a RegularGridInterpolator.")
+    state = pickled.state
+    if state.get("method") != "linear":
+        raise pickle.UnpicklingError(f"Unsupported interpolation method {state.get('method')}.")
+    if state.get("bounds_error"):
+        raise pickle.UnpicklingError("Interpolators with bounds_error=True are not supported.")
+    return LinearGridInterpolator(
+        grid=state["_grid"], values=state["_values"], fill_value=state.get("fill_value")
+    )
+
+
+class GridInterpolatorBackend(StrEnum):
+    """How the self-consumption tables are evaluated."""
+
+    SCIPY = "scipy"
+    NUMPY = "numpy"
+
+
+class SwitchableGridInterpolator:
+    """Evaluate a grid table with SciPy or with ``LinearGridInterpolator``.
+
+    ``backend`` is called on every evaluation and returns the backend to use, so
+    a configuration change applies to the next evaluation. The SciPy
+    ``RegularGridInterpolator`` is built from the same grid, values and fill
+    value the pickle holds (``method="linear"``, ``bounds_error=False``) on its
+    first use; with the NumPy backend SciPy is never imported.
+
+    Args:
+        linear: The table, as loaded by ``load_grid_interpolator()``.
+        backend: Returns the backend for the next evaluation.
+    """
+
+    def __init__(
+        self,
+        linear: LinearGridInterpolator,
+        backend: Callable[[], Union[str, GridInterpolatorBackend]],
+    ) -> None:
+        self._linear = linear
+        self._backend = backend
+        self._scipy: Optional[Any] = None
+        self._lock = threading.Lock()
+
+    @property
+    def grid(self) -> tuple[np.ndarray, ...]:
+        """Grid coordinates per axis."""
+        return self._linear.grid
+
+    @property
+    def values(self) -> np.ndarray:
+        """Values on the grid."""
+        return self._linear.values
+
+    @property
+    def fill_value(self) -> Union[float, None]:
+        """Value for points outside the grid."""
+        return self._linear.fill_value
+
+    def _scipy_interpolator(self) -> Any:
+        if self._scipy is None:
+            with self._lock:
+                if self._scipy is None:
+                    from scipy.interpolate import RegularGridInterpolator
+
+                    self._scipy = RegularGridInterpolator(
+                        self._linear.grid,
+                        self._linear.values,
+                        method="linear",
+                        bounds_error=False,
+                        fill_value=self._linear.fill_value,
+                    )
+        return self._scipy
+
+    def __call__(self, xi: Any) -> np.ndarray:
+        """Interpolate at points ``xi`` of shape ``(..., 2)``."""
+        if self._backend() == GridInterpolatorBackend.NUMPY:
+            return self._linear(xi)
+        return self._scipy_interpolator()(xi)
+
+
+def configured_grid_interpolator_backend() -> GridInterpolatorBackend:
+    """Backend from ``optimization.self_consumption_interpolator``; SciPy if unavailable."""
+    try:
+        from akkudoktoreos.core.coreabc import get_config
+
+        return GridInterpolatorBackend(get_config().optimization.self_consumption_interpolator)
+    except Exception:
+        return GridInterpolatorBackend.SCIPY

tests/test_gridinterpolator.py

@@ -0,0 +1,206 @@
+"""Tests for the SciPy-free evaluation of the self-consumption probability tables."""
+
+import io
+import os
+import pickle
+import subprocess
+import sys
+from pathlib import Path
+
+import numpy as np
+import pytest
+from scipy.interpolate import RegularGridInterpolator
+
+from akkudoktoreos.utils.gridinterpolator import (
+    GridInterpolatorBackend,
+    LinearGridInterpolator,
+    SwitchableGridInterpolator,
+    load_grid_interpolator,
+)
+
+DATA_DIR = Path(__file__).parent.parent / "src" / "akkudoktoreos" / "data"
+TABLES = [
+    DATA_DIR / "regular_grid_interpolator.pkl",
+    DATA_DIR / "genetic0_load_interpolator.pkl",
+]
+
+
+def _reference(grid, values, point, fill_value):
+    """SciPy's linear 2-D evaluation, written out in plain (unfused) float arithmetic."""
+    x0, x1 = (float(c) for c in point)
+    if np.isnan(x0) or np.isnan(x1):
+        return float("nan")
+    if not (grid[0][0] <= x0 <= grid[0][-1] and grid[1][0] <= x1 <= grid[1][-1]):
+        return float(fill_value)
+    corners = []
+    for axis, x in zip(grid, (x0, x1)):
+        i = int(np.searchsorted(axis, x, side="right")) - 1
+        i = min(max(i, 0), len(axis) - 2)
+        corners.append((i, (x - float(axis[i])) / (float(axis[i + 1]) - float(axis[i]))))
+    (i0, y0), (i1, y1) = corners
+    v = values
+    result = 0.0
+    result = result + float(v[i0, i1]) * (1 - y0) * (1 - y1)
+    result = result + float(v[i0, i1 + 1]) * (1 - y0) * y1
+    result = result + float(v[i0 + 1, i1]) * y0 * (1 - y1)
+    result = result + float(v[i0 + 1, i1 + 1]) * y0 * y1
+    return result
+
+
+def _points(grid, n=20000, seed=0):
+    rng = np.random.default_rng(seed)
+    g0, g1 = grid
+    random = np.column_stack(
+        [
+            rng.uniform(g0[0] - 100, g0[-1] + 100, n),
+            rng.uniform(g1[0] - 100, g1[-1] + 100, n),
+        ]
+    )
+    nodes = np.array(np.meshgrid(g0, g1, indexing="ij")).reshape(2, -1).T
+    # Query pattern of the EOS interpolator: one mean load, PV steps of 50 W.
+    pattern = np.vstack(
+        [
+            np.column_stack(
+                [np.full(int(g1[-1] // 50) + 1, load), np.arange(0.0, g1[-1] + 1.0, 50.0)]
+            )
+            for load in rng.uniform(g0[0], g0[-1], 50)
+        ]
+    )
+    edges = np.array([[g0[0], g1[0]], [g0[-1], g1[-1]], [g0[-1], g1[0]], [np.nan, g1[0]]])
+    return np.vstack([random, nodes, pattern, edges])
+
+
+@pytest.mark.parametrize("table", TABLES, ids=lambda path: path.name)
+def test_loaded_table_matches_scipy(table: Path):
+    with table.open("rb") as file:
+        scipy_interpolator = pickle.load(file)  # noqa: S301 - shipped data file
+    interpolator = load_grid_interpolator(table)
+
+    for ours, theirs in zip(interpolator.grid, scipy_interpolator.grid):
+        np.testing.assert_array_equal(ours, theirs)
+    np.testing.assert_array_equal(interpolator.values, scipy_interpolator.values)
+    assert interpolator.fill_value == scipy_interpolator.fill_value
+
+    points = _points(interpolator.grid)
+    ours = interpolator(points)
+    theirs = scipy_interpolator(points)
+    # Bit-identical to SciPy builds without fused multiply-add (Linux wheels);
+    # builds that fuse it (macOS arm64) differ only by rounding (|diff| <= ~1e-16).
+    np.testing.assert_array_equal(np.isnan(ours), np.isnan(theirs))
+    finite = ~np.isnan(ours)
+    np.testing.assert_allclose(ours[finite], theirs[finite], rtol=0, atol=1e-15)
+
+    # Exactly SciPy's arithmetic, evaluated without any fusion.
+    sample = points[:: max(1, len(points) // 3000)]
+    expected = np.array(
+        [
+            _reference(interpolator.grid, interpolator.values, point, interpolator.fill_value)
+            for point in sample
+        ]
+    )
+    np.testing.assert_array_equal(interpolator(sample), expected)
+
+
+def test_interpolator_matches_scipy_semantics_on_small_grid():
+    grid = (np.array([0.0, 1.0, 3.0]), np.array([10.0, 20.0]))
+    values = np.array([[1.0, 2.0], [3.0, 5.0], [-1.0, 0.5]])
+    ours = LinearGridInterpolator(grid, values, fill_value=0)
+    theirs = RegularGridInterpolator(grid, values, bounds_error=False, fill_value=0)
+    points = np.array(
+        [[0.0, 10.0], [3.0, 20.0], [2.5, 12.5], [-0.1, 15.0], [1.0, 20.1], [np.nan, 10.0]]
+    )
+    np.testing.assert_allclose(ours(points)[:-1], theirs(points)[:-1], rtol=0, atol=1e-15)
+    assert np.isnan(ours(points)[-1])
+    assert ours(points)[3] == 0.0
+    assert ours(np.array([[2.5, 12.5]])).shape == (1,)
+    assert ours(np.zeros((3, 4, 2))).shape == (3, 4)
+
+
+def test_interpolator_rejects_invalid_grids():
+    with pytest.raises(ValueError):
+        LinearGridInterpolator((np.array([0.0]), np.array([0.0, 1.0])), np.zeros((1, 2)))
+    with pytest.raises(ValueError):
+        LinearGridInterpolator((np.array([1.0, 0.0]), np.array([0.0, 1.0])), np.zeros((2, 2)))
+    with pytest.raises(ValueError):
+        LinearGridInterpolator((np.array([0.0, 1.0]), np.array([0.0, 1.0])), np.zeros((3, 2)))
+
+
+def test_loader_rejects_foreign_pickles():
+    with pytest.raises(pickle.UnpicklingError):
+        load_grid_interpolator(io.BytesIO(pickle.dumps(os.system)))
+    with pytest.raises(pickle.UnpicklingError):
+        load_grid_interpolator(io.BytesIO(pickle.dumps({"a": 1})))
+    method_nearest = RegularGridInterpolator(
+        (np.array([0.0, 1.0]), np.array([0.0, 1.0])), np.zeros((2, 2)), method="nearest"
+    )
+    with pytest.raises(pickle.UnpicklingError):
+        load_grid_interpolator(io.BytesIO(pickle.dumps(method_nearest)))
+
+
+def test_loading_the_tables_does_not_import_scipy():
+    code = (
+        "import sys\n"
+        "from akkudoktoreos.utils.gridinterpolator import load_grid_interpolator\n"
+        + "".join(f"load_grid_interpolator({str(table)!r})([[500.0, 100.0]])\n" for table in TABLES)
+        + "print('scipy' in sys.modules)\n"
+    )
+    result = subprocess.run(  # noqa: S603
+        [sys.executable, "-c", code], capture_output=True, text=True, check=True, timeout=120
+    )
+    assert result.stdout.strip().splitlines()[-1] == "False"
+
+
+@pytest.mark.parametrize("table", TABLES, ids=lambda path: path.name)
+def test_switchable_interpolator_defaults_to_scipy(table: Path):
+    with table.open("rb") as file:
+        scipy_interpolator = pickle.load(file)  # noqa: S301 - shipped data file
+    backend: dict[str, str] = {"value": GridInterpolatorBackend.SCIPY}
+    switchable = SwitchableGridInterpolator(load_grid_interpolator(table), lambda: backend["value"])
+    points = _points(switchable.grid, n=5000)
+
+    # SciPy backend: the same RegularGridInterpolator as the pickle, so exactly its results.
+    np.testing.assert_array_equal(switchable(points), scipy_interpolator(points))
+
+    # NumPy backend: LinearGridInterpolator; the switch applies to the next call.
+    backend["value"] = GridInterpolatorBackend.NUMPY
+    np.testing.assert_array_equal(switchable(points), load_grid_interpolator(table)(points))
+    backend["value"] = "scipy"
+    np.testing.assert_array_equal(switchable(points), scipy_interpolator(points))
+
+
+def test_setting_selects_the_backend(config_eos):
+    from akkudoktoreos.prediction.interpolator import get_eos_load_interpolator
+
+    assert config_eos.optimization.self_consumption_interpolator == GridInterpolatorBackend.SCIPY
+    interpolator = get_eos_load_interpolator().interpolator
+    points = _points(interpolator.grid, n=2000)
+    with (DATA_DIR / "regular_grid_interpolator.pkl").open("rb") as file:
+        scipy_values = pickle.load(file)(points)  # noqa: S301 - shipped data file
+    np.testing.assert_array_equal(interpolator(points), scipy_values)
+
+    config_eos.merge_settings_from_dict({"optimization": {"self_consumption_interpolator": "numpy"}})
+    assert config_eos.optimization.self_consumption_interpolator == GridInterpolatorBackend.NUMPY
+    linear = load_grid_interpolator(DATA_DIR / "regular_grid_interpolator.pkl")
+    np.testing.assert_array_equal(interpolator(points), linear(points))
+
+    with pytest.raises(ValueError):
+        config_eos.merge_settings_from_dict({"optimization": {"self_consumption_interpolator": "cubic"}})
+
+
+@pytest.mark.parametrize(
+    ("backend", "imports_scipy"), [("numpy", "False"), ("scipy", "True")]
+)
+def test_numpy_backend_does_not_import_scipy(backend: str, imports_scipy: str):
+    code = (
+        "import sys\n"
+        "from akkudoktoreos.utils.gridinterpolator import SwitchableGridInterpolator, load_grid_interpolator\n"
+        + "".join(
+            f"SwitchableGridInterpolator(load_grid_interpolator({str(table)!r}), lambda: {backend!r})([[500.0, 100.0]])\n"
+            for table in TABLES
+        )
+        + "print('scipy' in sys.modules)\n"
+    )
+    result = subprocess.run(  # noqa: S603
+        [sys.executable, "-c", code], capture_output=True, text=True, check=True, timeout=120
+    )
+    assert result.stdout.strip().splitlines()[-1] == imports_scipy