Interpolation API

The simses.interpolation module provides two fast scalar interpolation helpers used on the inner-loop hot path of the simulation. They are drop-in replacements for numpy.interp (1-D) and scipy.interpolate.RegularGridInterpolator (2-D bilinear) that avoid the ~2–80 µs of per-call argument sanitisation those libraries pay for a single scalar query.

Both helpers operate on plain Python lists and use bisect for the index search. Models should convert lookup-table arrays to lists once at construction time and pass them in on every step.

Out-of-bounds inputs raise rather than clip — silent clipping would mask integration overshoot, initialisation bugs, and upstream control errors.

Functions

simses.interpolation.interp1d_scalar(x, xp, fp)

Linear interpolation of a scalar against a sorted axis.

Parameters:

Name	Type	Description	Default
x	float	Query value.	required
xp	list[float]	Strictly ascending axis (Python list).	required
fp	list[float]	Function values at xp (Python list, same length).	required

Returns:

Type	Description
float	f(x) by linear interpolation.

Raises:

Type	Description
ValueError	If x lies outside [xp[0], xp[-1]].

Source code in src/simses/interpolation.py

def interp1d_scalar(x: float, xp: list[float], fp: list[float]) -> float:
    """Linear interpolation of a scalar against a sorted axis.

    Args:
        x:  Query value.
        xp: Strictly ascending axis (Python list).
        fp: Function values at ``xp`` (Python list, same length).

    Returns:
        ``f(x)`` by linear interpolation.

    Raises:
        ValueError: If ``x`` lies outside ``[xp[0], xp[-1]]``.
    """
    if x < xp[0] or x > xp[-1]:
        raise ValueError(f"x={x} out of bounds [{xp[0]}, {xp[-1]}]")
    if x == xp[-1]:
        return fp[-1]
    i = bisect.bisect_right(xp, x) - 1
    x0 = xp[i]
    return fp[i] + (fp[i + 1] - fp[i]) * (x - x0) / (xp[i + 1] - x0)

simses.interpolation.interp2d_scalar(x, y, xp, yp, mat)

Bilinear interpolation of a scalar (x, y) against a 2-D LUT.

mat is indexed as mat[i][j] with xp[i] and yp[j] as the grid coordinates. Both axes must be strictly ascending. The result is bit-for-bit equivalent to RegularGridInterpolator((xp, yp), mat)((x, y)) at every grid node and at all interior points (machine epsilon).

Parameters:

Name	Type	Description	Default
x	float	Query value along the first axis.	required
y	float	Query value along the second axis.	required
xp	list[float]	Strictly ascending first-axis values.	required
yp	list[float]	Strictly ascending second-axis values.	required
mat	list[list[float]]	2-D table; mat[i][j] is the value at (xp[i], yp[j]).	required

Returns:

Type	Description
float	f(x, y) by bilinear interpolation.

Raises:

Type	Description
ValueError	If x or y lies outside its respective LUT range.

Source code in src/simses/interpolation.py

def interp2d_scalar(
    x: float,
    y: float,
    xp: list[float],
    yp: list[float],
    mat: list[list[float]],
) -> float:
    """Bilinear interpolation of a scalar (x, y) against a 2-D LUT.

    ``mat`` is indexed as ``mat[i][j]`` with ``xp[i]`` and ``yp[j]`` as the
    grid coordinates. Both axes must be strictly ascending. The result is
    bit-for-bit equivalent to ``RegularGridInterpolator((xp, yp), mat)((x, y))``
    at every grid node and at all interior points (machine epsilon).

    Args:
        x:   Query value along the first axis.
        y:   Query value along the second axis.
        xp:  Strictly ascending first-axis values.
        yp:  Strictly ascending second-axis values.
        mat: 2-D table; ``mat[i][j]`` is the value at ``(xp[i], yp[j])``.

    Returns:
        ``f(x, y)`` by bilinear interpolation.

    Raises:
        ValueError: If ``x`` or ``y`` lies outside its respective LUT range.
    """
    if x < xp[0] or x > xp[-1]:
        raise ValueError(f"x={x} out of bounds [{xp[0]}, {xp[-1]}]")
    if y < yp[0] or y > yp[-1]:
        raise ValueError(f"y={y} out of bounds [{yp[0]}, {yp[-1]}]")

    if x == xp[-1]:
        i = len(xp) - 2
        u = 1.0
    else:
        i = bisect.bisect_right(xp, x) - 1
        x0 = xp[i]
        u = (x - x0) / (xp[i + 1] - x0)

    if y == yp[-1]:
        j = len(yp) - 2
        v = 1.0
    else:
        j = bisect.bisect_right(yp, y) - 1
        y0 = yp[j]
        v = (y - y0) / (yp[j + 1] - y0)

    row_i = mat[i]
    row_i1 = mat[i + 1]
    f00 = row_i[j]
    f01 = row_i[j + 1]
    f10 = row_i1[j]
    f11 = row_i1[j + 1]

    one_minus_u = 1.0 - u
    one_minus_v = 1.0 - v
    return one_minus_u * one_minus_v * f00 + u * one_minus_v * f10 + one_minus_u * v * f01 + u * v * f11