Skip to content

image/svg+xml

Reference🔗

berhoel.splines 🔗

Various Splines.

Modules:

  • akima –

    Calculate interpolating spline with the akima subspline algorithm.

  • line –

    Class for representing lines between two points.

  • renner –

    Class for calculating Renner subsplines for given nodes.

  • tests –

Classes:

  • AkimaSubspline –

    Akima subsplines interpolator acc. to [^l1], chapter 13.

Classes🔗

berhoel.splines.AkimaSubspline 🔗

AkimaSubspline(
    x: ToFloat1D,
    y: ToComplexND,
    axis: int = 0,
    extrapolate: bool | None | Literal["periodic"] = None,
)

Akima subsplines interpolator acc. to 1, chapter 13.

The original definition is found in 2.

Fit piecewise cubic polynomials, given vectors \(x\) and \(y\). The interpolation method by Akima uses a continuously differentiable sub-spline built from piecewise cubic polynomials. The resultant curve passes through the given data points and will appear smooth and natural.

Notes🔗

Use only for precise data, as the fitted curve passes through the given points exactly. This routine is useful for plotting a pleasingly smooth curve through a few given points for purposes of plotting.

Given: \(\left(x_i,y_i\right)\), \(i=0(1)n, n>=4, x_i<x_{i+1}\) for \(i=0(1)n-1\), not periodic or (with \(y_n = y_0\)) periodic.

Let \(m_i = (y_{i+1} - y_i) / (x_{i+1} - x_i)\) be the slopes of the interval \(\left[x_i, x_{i+1}\right)\).

Further slopes are for the not periodic variant:

\[ \begin{aligned} m_{-2} &:= 3 m_0 - 2 m_1 \\ m_{-1} &:= 2 m_0 - m_1 \\ m_n &:= 2 m_{n-1} - m_{n-2} \\ m_{n+1} &:= 3 m_{n-1} - 2 m_{n-2} \end{aligned} \]

For the periodic variant:

\[ \begin{aligned} m_{-2} &:= m_{n-2} \\ m_{-1} &:= m_{n-1} \\ m_n &:= m_0 \\ m_{n+1} &:= m_1 \end{aligned} \]
Examples🔗

If the 3 or more bases one side of the subspline are colinear, a corner will appear in the interpoltion spline.

# mkdocs: render
import numpy as np
from berhoel.splines import AkimaSubspline
import matplotlib.pyplot as plt
x = np.linspace(1, 7, 7)
y = np.array([-1, -1, -1, 0, 1, 1, 1])
xs = np.linspace(min(x), max(x), num=100)
y_akima = AkimaSubspline(x, y)(xs)
fig, ax = plt.subplots()
ax.plot(x, y, "o", label="bases")
ax.plot(xs, y_akima, label="akima")
ax.legend()
fig.show()

The AkimaSubspline implementation allows for periodic splines, in contrast to scipy.interpolate.Akima1DInterpolator.

# mkdocs: render
from berhoel.splines import AkimaSubspline
x = np.linspace(-np.pi, np.pi, 6)
y = np.sin(x)
xs = np.linspace(-np.pi, np.pi, 100) * 1.7
y_akima = AkimaSubspline(x, y, extrapolate="periodic")(xs)
fig, ax = plt.subplots()
ax.plot(x, y, "o", label="bases")
ax.plot(xs, y_akima, label="akima")
ax.plot(xs, np.sin(xs), label="sin")
ax.legend()
fig.show()
References🔗

  1. Engeln-Müllges, Gisela, und Fritz Reutter. Numerik-Algorithmen: Entscheidungshilfe zur Auswahl und Nutzung. 8., Neubearb. und erw. Aufl. Düsseldorf: VDI-Verl, 1996., ISBN 3-18-401539-4 ↩

  2. Akima, Hiroshi. „A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures“. Journal of the ACM 17, Nr. 4 (Oktober 1970): 589-602. https://doi.org/10.1145/321607.321609. ↩

Parameters:

  • x 🔗
    (ToFloat1D) –

    1-D array of monotonically increasing real values.

  • y 🔗
    (ToComplexND) –

    N-D array of real values. The length of y along the interpolation axis must be equal to the length of x. Use the axis parameter to select the interpolation axis.

  • axis 🔗
    (int, default: 0 ) –

    Axis in the y array corresponding to the x-coordinate values. Defaults to axis=0.

  • extrapolate 🔗
    (bool | None | Literal['periodic'], default: None ) –

    If bool, determines whether to extrapolate to out-of-bounds points based on first and last intervals, or to return NaNs. If None, extrapolate is set to False. If “periodic” the starting and finishing slopes are calculate for allowing interpolating periodic functions.

Methods:

  • __call__ –

    Evaluate the piecewise polynomial or its derivative.

Source code in src/berhoel/splines/akima.py 
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
def __init__(
    self,
    x: onp.ToFloat1D,
    y: onp.ToComplexND,
    axis: int = 0,
    extrapolate: bool | None | Literal["periodic"] = None,  # noqa: FBT001
):
    """Initialize ``AkimaSubspline``.

    Args:
      x: 1-D array of monotonically increasing real values.
      y: N-D array of real values. The length of ``y`` along the interpolation axis
          must be equal to the length of ``x``. Use the ``axis`` parameter to
          select the interpolation axis.
      axis: Axis in the ``y`` array corresponding to the x-coordinate values.
          Defaults to ``axis=0``.
      extrapolate: If bool, determines whether to extrapolate to out-of-bounds
        points based on first and last intervals, or to return NaNs. If None,
        ``extrapolate`` is set to False. If "periodic" the starting and finishing
        slopes are calculate for allowing interpolating periodic functions.
    """
    self.x, dx, self.y, axis, _ = prepare_input(x, y, axis)

    if np.iscomplexobj(self.y):
        msg = "`AkimaSubspline` only works with real values for `y`."
        raise NotImplementedError(msg)

    # Akima extrapolation historically False, this implementation also allows
    # "periodic".
    extrapolate = False if extrapolate is None else extrapolate

    n = self.x.shape[0] - 1

    m_shape = np.array(self.y.shape)
    y_multidim = self.y.ndim == 2  # noqa:PLR2004
    # To avoid shifiting indices this implementation works with indices -1 and -2
    # as described in [1]_.
    m = np.zeros([self.y.shape[0] + 4, 2] if y_multidim else (self.y.shape[0] + 4))

    m[:n] = np.diff(self.y, axis=0) / (dx if (len(m_shape) == 1) else dx[:, None])
    self.periodic = extrapolate == "periodic"
    if not self.periodic:
        m[-2:] = 3 * m[0] - 2 * m[1], 2 * m[0] - m[1]
        m[n : n + 2] = 2 * m[n - 1] - m[n - 2], 3 * m[n - 1] - 2 * m[n - 2]
    else:
        if not np.allclose(self.y[0], self.y[-1]):
            msg = "For periodic curves y[0] hast to be equal to y[n]."
            raise ValueError(msg)
        m[-2:] = m[n - 2], m[n - 1]
        m[n : n + 2] = m[0:2]

    i = np.linspace(0, n, n + 1, dtype=int)
    NE = np.abs(m[i + 1] - m[i]) + np.abs(m[i - 1] - m[i - 2])  # noqa:N806
    h = np.diff(self.x)
    self.b, self.c, self.d = map(np.empty, ([n, self.y.ndim],) * 3)
    t_L, t_R = np.empty(n + 1), np.empty(n + 1)  # noqa:N806
    for j in range(self.y.ndim):
        line_splice = (np.s_[:, j]) if (y_multidim) else np.s_[:]
        # where NE>0
        i_s = i[NE[line_splice] > 0]
        i_s_splice, i_s_splice_1, i_s_splice_2 = (
            (np.s_[i_s, j], np.s_[i_s - 1, j], np.s_[i_s - 2, j])
            if y_multidim
            else (i_s, i_s - 1, i_s - 2)
        )
        alpha = np.abs(m[i_s_splice_1] - m[i_s_splice_2]) / NE[i_s_splice]
        t_L[i_s] = m[i_s_splice_1] + alpha * (m[i_s_splice] - m[i_s_splice_1])
        t_R[i_s] = t_L[i_s]
        # where NE==0
        i_s = i[NE[line_splice] == 0]
        i_s_splice, i_s_splice_1, m_splice = (
            (np.s_[i_s, j], np.s_[i_s - 1, j], np.s_[:n, j])
            if y_multidim
            else (i_s, i_s - 1, np.s_[:n])
        )
        t_L[i_s], t_R[i_s] = m[i_s_splice_1], m[i_s_splice]
        self.c[:, j], self.d[:, j] = (
            1 / h * (3 * m[m_splice] - 2 * t_R[:n] - t_L[1 : n + 1]),
            1 / h / h * (t_R[:n] + t_L[1 : n + 1] - 2 * m[m_splice]),
        )
        self.b[:, j] = t_R[:n]
Functions🔗
__call__ 🔗
__call__(x: NDArray, nu: int = 0, extrapolate: str | None = None) -> NDArray

Evaluate the piecewise polynomial or its derivative.

Parameters:

  • x 🔗 (NDArray) –

    Points to evaluate the interpolant at.

  • nu 🔗 (int, default: 0 ) –

    Order of derivative to evaluate. Must be non-negative.

  • extrapolate 🔗 (str | None, default: None ) –

    Extrapolation level

Returns:

  • y ( NDArray ) –

    Interpolated values. Shape is determined by replacing the interpolation axis in the original array with the shape of x.

Notes🔗

Derivatives are evaluated piecewise for each polynomial segment, even if the polynomial is not differentiable at the breakpoints. The polynomial intervals are considered half-open, \([a, b)\), except for the last interval which is closed \([a, b]\).

Source code in src/berhoel/splines/akima.py 
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
def __call__(
    self, x: npt.NDArray, nu: int = 0, extrapolate: str | None = None
) -> npt.NDArray:
    """Evaluate the piecewise polynomial or its derivative.

    Args:
      x: Points to evaluate the interpolant at.
      nu: Order of derivative to evaluate. Must be non-negative.
      extrapolate: Extrapolation level

    Returns:
      y : Interpolated values. Shape is determined by replacing
        the interpolation axis in the original array with the shape of x.

    # Notes

    Derivatives are evaluated piecewise for each polynomial
    segment, even if the polynomial is not differentiable at the
    breakpoints. The polynomial intervals are considered half-open,
    $[a, b)$, except for the last interval which is closed
    $[a, b]$.
    """
    # import here to avoid cyclic imports
    from . import intervall  # noqa: PLC0415

    periodic = self.periodic

    if extrapolate is not None:
        msg = "extrapolate must be set on class level."
        raise ValueError(msg)

    _x = np.array(x)

    if periodic:
        min_val, max_val = np.min(self.x), np.max(self.x)
        diff = max_val - min_val

        while max(_x) > max_val:
            _x[_x > max_val] -= diff
        while min(_x) < min_val:
            _x[_x < min_val] += diff

    indices = intervall(_x, self.x)
    _x -= self.x[indices]
    y_multidim = self.y.ndim == 2  # noqa:PLR2004
    param_slice, x_slice = (
        (indices, np.s_[:, None]) if y_multidim else (np.s_[indices, 0], np.s_[:])
    )
    return {
        0: (
            (
                ((self.d[param_slice] * _x[x_slice]) + self.c[param_slice])
                * _x[x_slice]
                + self.b[param_slice]
            )
            * _x[x_slice]
        )
        + self.y[indices],
        1: ((3 * self.d[param_slice] * _x[x_slice]) + 2 * self.c[param_slice])
        * _x[x_slice]
        + self.b[param_slice],
        2: (6 * self.d[param_slice] * _x[x_slice]) + 2 * self.c[param_slice],
        3: 6 * self.d[param_slice],
    }.get(
        nu,
        (np.zeros((len(x), self.y.shape[1])) if y_multidim else np.zeros(len(x))),
    )

Functions🔗

berhoel.splines.intervall 🔗

intervall(vals: float | NDArray, boundaries: NDArray) -> NDArray

Find the interval number of a decomposition.

This is a generally usable procedure for finding an interval that contains a given value xvalue. Given is a strictly monotonically increasing sequence of values \(x_i, i=0(1)n\). The index \(i_x\) is sought such that \(x_{i_x} <= x_\text{value} < x_{i_x+1}\). If \(x_\text{value} < x_0\) or \(x_\text{value} >= x_{n-1}\), \(i_x\) gets the value \(0\) or \(n-1\).

Parameters:

Returns:

  • NDArray –

    index ix to be searched for

Source code in src/berhoel/splines/__init__.py 
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
def intervall(vals: float | npt.NDArray, boundaries: npt.NDArray) -> npt.NDArray:
    r"""Find the interval number of a decomposition.

    This is a generally usable procedure for finding an interval that
    contains a given value xvalue.  Given is a strictly monotonically
    increasing sequence of values $x_i, i=0(1)n$.  The index $i_x$ is
    sought such that $x_{i_x} <= x_\text{value} < x_{i_x+1}$.  If
    $x_\text{value} < x_0$ or $x_\text{value} >= x_{n-1}$, $i_x$ gets
    the value $0$ or $n-1$.

    Args:
      vals: value whose interval is being searched
      boundaries: field with the interval limits

    Returns:
      index ix to be searched for
    """
    if isinstance(vals, float):
        vals = np.array([vals], float)
    boundaries = np.array(boundaries, float)
    boundaries[-1] = 1e32
    res = (boundaries > vals[:, np.newaxis]).argmax(axis=1) - 1
    res[res < 0] = 0
    return res

berhoel.splines.line 🔗

Class for representing lines between two points.

Requires NumPy and, for the exampes, Michael Haggerty’s Gnuplot module.

Classes:

  • Line –

    Represent a line.

Classes🔗

berhoel.splines.line.Line 🔗

Line(points: object)

Represent a line.

Methods:

  • __call__ –

    Calculate steps points for line.

  • intersect –

    Calculate intersection of line with area.

Source code in src/berhoel/splines/line.py 
34
35
36
37
38
39
def __init__(self, points: object) -> None:
    """Initialize class instance."""
    self.points = np.asarray(points, float)
    self.x = self.points[0]
    self.y = self.points[1]
    self.p = self.y - self.x
Functions🔗
__call__ 🔗
__call__(steps: int = 10) -> NDArray

Calculate steps points for line.

Source code in src/berhoel/splines/line.py 
41
42
43
def __call__(self, steps: int = 10) -> npt.NDArray:
    """Calculate `steps` points for line."""
    return np.outer(np.arange(steps + 1, dtype=float) / float(steps), self.p)
intersect 🔗
intersect(norm: tuple[float, ...], dist: float) -> NDArray

Calculate intersection of line with area.

Parameters:

Source code in src/berhoel/splines/line.py 
45
46
47
48
49
50
51
52
53
54
55
def intersect(self, norm: tuple[float, ...], dist: float) -> npt.NDArray:
    """Calculate intersection of line with area.

    Args:
       norm: normal of area
       dist: distance of area of zero point
    """
    out = (dist - np.dot(norm, self.x)) / np.dot(norm, self.p)
    if out >= 0.0 and out <= 1.0:
        return self.x + out * self.p
    return np.empty(0)

berhoel.splines.renner 🔗

Class for calculating Renner subsplines for given nodes.

Requires ‘NumPy’ and, for the exampes, Michael Haggerty’s ‘Gnuplot’ module and Konrad Hinsens ‘VRML’ module.

References:🔗

.. [1] Engeln-Müllges, Gisela, und Fritz Reutter. Numerik-Algorithmen:
       Entscheidungshilfe zur Auswahl und Nutzung. 8., Neubearb. und erw. Aufl.
       Düsseldorf: VDI-Verl, 1996., ISBN 3-18-401539-4

Classes:

Classes🔗

berhoel.splines.renner.Renner 🔗

Renner(points: NDArray, closed: bool = False)

Renner splines.

Methods:

  • __call__ –

    Calculate ‘steps’ interpolated points between each two input points.

  • intersect –

    Return a list of intersections.

Source code in src/berhoel/splines/renner.py 
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
def __init__(
    self,
    points: npt.NDArray,
    closed: bool = False,  # noqa:FBT001,FBT002
) -> None:
    """Intialize class instance."""
    self.points = np.asarray(points)
    self.oshape = self.points.shape
    ## We need at least 5 points
    if self.oshape[0] < 5:  # noqa:PLR2004
        msg = "Renner requires at least fife points"
        raise TypeError(msg)
    ## 1. calculate chord vectors
    s = np.zeros(np.array((3, 0)) + self.oshape, float)
    s[2:-2] = self.points[1:] - self.points[:-1]
    ## 2. calculate additional chord vectors
    if closed:
        ## 2.1 closed curve (P_n == P_0)
        s[0:2] = s[-4:-2]
        s[-2:] = s[2:4]
    else:
        ## 2.2 not closed curve (P_n != P_0)
        s[0] = 3.0 * s[2] - 2.0 * s[3]
        s[1] = 2.0 * s[2] - s[3]
        s[-2] = 2.0 * s[-3] - s[-4]
        s[-1] = 3.0 * s[-3] - 2.0 * s[-4]
    ## 3. calculate the standard vectors for the chord vectors
    s0 = np.zeros(s.shape, float)
    ## chordlength can be 0.
    for i, chord in enumerate(s):
        length = math.sqrt(np.sum(np.power(chord, 2)))
        if length > eps:
            s0[i, :] = chord / length
        else:
            s0[i, :] = np.zeros(self.oshape[1], float)
    ## 4. calculate the lefthand- and righthand-tangent unitiy vectors
    t_l = np.zeros(np.array(self.oshape), float)
    t_r = np.zeros(np.array(self.oshape), float)
    for n in range(self.oshape[0]):
        ## n+2 := i for s0
        ne = math.sqrt(max(0.0, 1.0 - np.dot(s0[n], s0[n + 1]) ** 2)) + math.sqrt(
            max(0.0, 1.0 - np.dot(s0[n + 2], s0[n + 3]) ** 2)
        )
        if self.NE_eps < ne:
            alp = math.sqrt(max(0.0, 1.0 - np.dot(s0[n], s0[n + 1]) ** 2)) / ne
            t_l[n, :] = (s[n + 1] + alp * (s[n + 2] - s[n + 1])) / math.sqrt(
                np.sum(np.power(s[n + 1] + alp * (s[n + 2] - s[n + 1]), 2))
            )
            t_r[n, :] = t_l[n, :]
        else:
            t_l[n, :] = s0[n + 1]
            t_r[n, :] = s0[n + 2]
    ## 5. calculate the length t_i of the parameterintervalls
    a = np.ones(self.oshape[0] - 1) * 16.0 - np.sum(pow(t_r[:-1] + t_l[1:], 2), 1)
    b = np.zeros(self.oshape[0] - 1, float)
    for i, (s_i, t_r_i, t_l_i) in enumerate(
        zip(s[2:-2], t_r[:-1], t_l[1:], strict=False)
    ):
        b[i] = np.dot(6.0 * s_i, (t_r_i + t_l_i))
    c = 36.0 * np.sum(pow(s[2:-2], 2), 1)
    self.T = (-b + math.sqrt(pow(b, 2) + a * c)) / a
    ## 6. calculate the coefficients for the Renner-subsplines
    self.a = self.points[:-1]
    self.b = t_r[:-1]
    self.c = np.zeros(np.array((-1, 0)) + self.oshape, float)
    for i, (s_i, t_i, t_r_i, t_l_i) in enumerate(
        zip(s[2:-2], self.T, t_r[:-1], t_l[1:], strict=False)
    ):
        self.c[i, :] = 3.0 * s_i / t_i**2 - (2.0 * t_r_i + t_l_i) / t_i
    self.d = np.zeros(np.array((-1, 0)) + self.oshape, float)
    for i, (s_i, t_i, t_r_i, t_l_i) in enumerate(
        zip(s[2:-2], self.T, t_r[:-1], t_l[1:], strict=False)
    ):
        self.d[i, :] = (t_r_i + t_l_i) / t_i**2 - 2.0 * s_i / t_i**3
    self.polynoms = zip(self.a, self.b, self.c, self.d, strict=False)
Functions🔗
__call__ 🔗
__call__(steps: int = 10) -> NDArray

Calculate ‘steps’ interpolated points between each two input points.

Source code in src/berhoel/splines/renner.py 
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
def __call__(self, steps: int = 10) -> npt.NDArray:
    """Calculate 'steps' interpolated points between each two input points."""
    out = np.zeros(
        ((np.array(self.oshape) - (1, 0)) * np.array((steps, 1))) + ((1, 0)), float
    )
    points = np.arange(steps, dtype=float) / float(steps)
    irange = range(steps)
    for i, (t_i, coeff_) in enumerate(zip(self.T, self.polynoms, strict=False)):
        coeff = np.asarray(coeff_)
        for j, t in zip(irange, points * t_i, strict=False):
            out[i * steps + j, :] = np.sum(
                coeff * (t[:, np.newaxis] ** np.arange(coeff.shape[0])), axis=1
            )
    out[-1, :] = np.sum(
        coeff * (self.T[-1][:, np.newaxis] ** np.arange(coeff.shape[0])), axis=1
    )
    return out
intersect 🔗
intersect(norm: float, dist: float) -> list[float]

Return a list of intersections.

Source code in src/berhoel/splines/renner.py 
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
def intersect(self, norm: float, dist: float) -> list[float]:
    """Return a list of intersections."""
    res_t = []
    res_c = []
    tsum = 0
    n = 0
    for t, coeff in zip(self.T, self.polynoms, strict=False):
        n = n + 1
        values = _normalize(np.dot(np.array(coeff), np.asarray(norm)).tolist())
        values[0] = values[0] - dist
        values = values[::-1]
        num = values.shape[0]
        if num < 2:  # noqa:PLR2004
            x = np.array([0.0, t] if num == 0 else [])
        else:
            # Build the matrix, whos "characteristic polynomial"
            # is the polynominal we152've build.
            matrix = np.zeros((num - 1, num - 1), complex)
            if num > 2:  # noqa:PLR2004
                matrix[1:, :-1] = np.identity(num - 2)
            matrix[0, :] = -values[1:] / values[0]
            ## The Eigenvalues of the matrix are the roots of the polynom.
            x = np.array(la.eigvals(matrix))
            ## We only want real values and values in [0., t]
            r = np.compress(np.less_equal(abs(x.imag), self.EPS), x.real)
            x = np.compress(
                np.logical_and(np.greater_equal(r, 0.0), np.less_equal(r, t)), r
            )
        for val in x:
            res_t.append(tsum + val)
            res_c.append((val, np.asarray(coeff)))
        tsum = tsum + t
    res_t.sort()
    delete = [
        i for i in range(len(res_t) - 1) if (res_t[i + 1] - res_t[i]) < self.EPS
    ]
    delete.reverse()
    for i in delete:
        del res_c[i]
    return [np.sum(coe * (val ** np.arange(coe.shape[0]))) for val, coe in res_c]