Coverage for /usr/lib/python3/dist-packages/sympy/functions/special/bsplines.py: 11%
107 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from sympy.core import S, sympify
2from sympy.core.symbol import (Dummy, symbols)
3from sympy.functions import Piecewise, piecewise_fold
4from sympy.logic.boolalg import And
5from sympy.sets.sets import Interval
7from functools import lru_cache
10def _ivl(cond, x):
11 """return the interval corresponding to the condition
13 Conditions in spline's Piecewise give the range over
14 which an expression is valid like (lo <= x) & (x <= hi).
15 This function returns (lo, hi).
16 """
17 if isinstance(cond, And) and len(cond.args) == 2:
18 a, b = cond.args
19 if a.lts == x:
20 a, b = b, a
21 return a.lts, b.gts
22 raise TypeError('unexpected cond type: %s' % cond)
25def _add_splines(c, b1, d, b2, x):
26 """Construct c*b1 + d*b2."""
28 if S.Zero in (b1, c):
29 rv = piecewise_fold(d * b2)
30 elif S.Zero in (b2, d):
31 rv = piecewise_fold(c * b1)
32 else:
33 new_args = []
34 # Just combining the Piecewise without any fancy optimization
35 p1 = piecewise_fold(c * b1)
36 p2 = piecewise_fold(d * b2)
38 # Search all Piecewise arguments except (0, True)
39 p2args = list(p2.args[:-1])
41 # This merging algorithm assumes the conditions in
42 # p1 and p2 are sorted
43 for arg in p1.args[:-1]:
44 expr = arg.expr
45 cond = arg.cond
47 lower = _ivl(cond, x)[0]
49 # Check p2 for matching conditions that can be merged
50 for i, arg2 in enumerate(p2args):
51 expr2 = arg2.expr
52 cond2 = arg2.cond
54 lower_2, upper_2 = _ivl(cond2, x)
55 if cond2 == cond:
56 # Conditions match, join expressions
57 expr += expr2
58 # Remove matching element
59 del p2args[i]
60 # No need to check the rest
61 break
62 elif lower_2 < lower and upper_2 <= lower:
63 # Check if arg2 condition smaller than arg1,
64 # add to new_args by itself (no match expected
65 # in p1)
66 new_args.append(arg2)
67 del p2args[i]
68 break
70 # Checked all, add expr and cond
71 new_args.append((expr, cond))
73 # Add remaining items from p2args
74 new_args.extend(p2args)
76 # Add final (0, True)
77 new_args.append((0, True))
79 rv = Piecewise(*new_args, evaluate=False)
81 return rv.expand()
84@lru_cache(maxsize=128)
85def bspline_basis(d, knots, n, x):
86 """
87 The $n$-th B-spline at $x$ of degree $d$ with knots.
89 Explanation
90 ===========
92 B-Splines are piecewise polynomials of degree $d$. They are defined on a
93 set of knots, which is a sequence of integers or floats.
95 Examples
96 ========
98 The 0th degree splines have a value of 1 on a single interval:
100 >>> from sympy import bspline_basis
101 >>> from sympy.abc import x
102 >>> d = 0
103 >>> knots = tuple(range(5))
104 >>> bspline_basis(d, knots, 0, x)
105 Piecewise((1, (x >= 0) & (x <= 1)), (0, True))
107 For a given ``(d, knots)`` there are ``len(knots)-d-1`` B-splines
108 defined, that are indexed by ``n`` (starting at 0).
110 Here is an example of a cubic B-spline:
112 >>> bspline_basis(3, tuple(range(5)), 0, x)
113 Piecewise((x**3/6, (x >= 0) & (x <= 1)),
114 (-x**3/2 + 2*x**2 - 2*x + 2/3,
115 (x >= 1) & (x <= 2)),
116 (x**3/2 - 4*x**2 + 10*x - 22/3,
117 (x >= 2) & (x <= 3)),
118 (-x**3/6 + 2*x**2 - 8*x + 32/3,
119 (x >= 3) & (x <= 4)),
120 (0, True))
122 By repeating knot points, you can introduce discontinuities in the
123 B-splines and their derivatives:
125 >>> d = 1
126 >>> knots = (0, 0, 2, 3, 4)
127 >>> bspline_basis(d, knots, 0, x)
128 Piecewise((1 - x/2, (x >= 0) & (x <= 2)), (0, True))
130 It is quite time consuming to construct and evaluate B-splines. If
131 you need to evaluate a B-spline many times, it is best to lambdify them
132 first:
134 >>> from sympy import lambdify
135 >>> d = 3
136 >>> knots = tuple(range(10))
137 >>> b0 = bspline_basis(d, knots, 0, x)
138 >>> f = lambdify(x, b0)
139 >>> y = f(0.5)
141 Parameters
142 ==========
144 d : integer
145 degree of bspline
147 knots : list of integer values
148 list of knots points of bspline
150 n : integer
151 $n$-th B-spline
153 x : symbol
155 See Also
156 ========
158 bspline_basis_set
160 References
161 ==========
163 .. [1] https://en.wikipedia.org/wiki/B-spline
165 """
166 # make sure x has no assumptions so conditions don't evaluate
167 xvar = x
168 x = Dummy()
170 knots = tuple(sympify(k) for k in knots)
171 d = int(d)
172 n = int(n)
173 n_knots = len(knots)
174 n_intervals = n_knots - 1
175 if n + d + 1 > n_intervals:
176 raise ValueError("n + d + 1 must not exceed len(knots) - 1")
177 if d == 0:
178 result = Piecewise(
179 (S.One, Interval(knots[n], knots[n + 1]).contains(x)), (0, True)
180 )
181 elif d > 0:
182 denom = knots[n + d + 1] - knots[n + 1]
183 if denom != S.Zero:
184 B = (knots[n + d + 1] - x) / denom
185 b2 = bspline_basis(d - 1, knots, n + 1, x)
186 else:
187 b2 = B = S.Zero
189 denom = knots[n + d] - knots[n]
190 if denom != S.Zero:
191 A = (x - knots[n]) / denom
192 b1 = bspline_basis(d - 1, knots, n, x)
193 else:
194 b1 = A = S.Zero
196 result = _add_splines(A, b1, B, b2, x)
197 else:
198 raise ValueError("degree must be non-negative: %r" % n)
200 # return result with user-given x
201 return result.xreplace({x: xvar})
204def bspline_basis_set(d, knots, x):
205 """
206 Return the ``len(knots)-d-1`` B-splines at *x* of degree *d*
207 with *knots*.
209 Explanation
210 ===========
212 This function returns a list of piecewise polynomials that are the
213 ``len(knots)-d-1`` B-splines of degree *d* for the given knots.
214 This function calls ``bspline_basis(d, knots, n, x)`` for different
215 values of *n*.
217 Examples
218 ========
220 >>> from sympy import bspline_basis_set
221 >>> from sympy.abc import x
222 >>> d = 2
223 >>> knots = range(5)
224 >>> splines = bspline_basis_set(d, knots, x)
225 >>> splines
226 [Piecewise((x**2/2, (x >= 0) & (x <= 1)),
227 (-x**2 + 3*x - 3/2, (x >= 1) & (x <= 2)),
228 (x**2/2 - 3*x + 9/2, (x >= 2) & (x <= 3)),
229 (0, True)),
230 Piecewise((x**2/2 - x + 1/2, (x >= 1) & (x <= 2)),
231 (-x**2 + 5*x - 11/2, (x >= 2) & (x <= 3)),
232 (x**2/2 - 4*x + 8, (x >= 3) & (x <= 4)),
233 (0, True))]
235 Parameters
236 ==========
238 d : integer
239 degree of bspline
241 knots : list of integers
242 list of knots points of bspline
244 x : symbol
246 See Also
247 ========
249 bspline_basis
251 """
252 n_splines = len(knots) - d - 1
253 return [bspline_basis(d, tuple(knots), i, x) for i in range(n_splines)]
256def interpolating_spline(d, x, X, Y):
257 """
258 Return spline of degree *d*, passing through the given *X*
259 and *Y* values.
261 Explanation
262 ===========
264 This function returns a piecewise function such that each part is
265 a polynomial of degree not greater than *d*. The value of *d*
266 must be 1 or greater and the values of *X* must be strictly
267 increasing.
269 Examples
270 ========
272 >>> from sympy import interpolating_spline
273 >>> from sympy.abc import x
274 >>> interpolating_spline(1, x, [1, 2, 4, 7], [3, 6, 5, 7])
275 Piecewise((3*x, (x >= 1) & (x <= 2)),
276 (7 - x/2, (x >= 2) & (x <= 4)),
277 (2*x/3 + 7/3, (x >= 4) & (x <= 7)))
278 >>> interpolating_spline(3, x, [-2, 0, 1, 3, 4], [4, 2, 1, 1, 3])
279 Piecewise((7*x**3/117 + 7*x**2/117 - 131*x/117 + 2, (x >= -2) & (x <= 1)),
280 (10*x**3/117 - 2*x**2/117 - 122*x/117 + 77/39, (x >= 1) & (x <= 4)))
282 Parameters
283 ==========
285 d : integer
286 Degree of Bspline strictly greater than equal to one
288 x : symbol
290 X : list of strictly increasing real values
291 list of X coordinates through which the spline passes
293 Y : list of real values
294 list of corresponding Y coordinates through which the spline passes
296 See Also
297 ========
299 bspline_basis_set, interpolating_poly
301 """
302 from sympy.solvers.solveset import linsolve
303 from sympy.matrices.dense import Matrix
305 # Input sanitization
306 d = sympify(d)
307 if not (d.is_Integer and d.is_positive):
308 raise ValueError("Spline degree must be a positive integer, not %s." % d)
309 if len(X) != len(Y):
310 raise ValueError("Number of X and Y coordinates must be the same.")
311 if len(X) < d + 1:
312 raise ValueError("Degree must be less than the number of control points.")
313 if not all(a < b for a, b in zip(X, X[1:])):
314 raise ValueError("The x-coordinates must be strictly increasing.")
315 X = [sympify(i) for i in X]
317 # Evaluating knots value
318 if d.is_odd:
319 j = (d + 1) // 2
320 interior_knots = X[j:-j]
321 else:
322 j = d // 2
323 interior_knots = [
324 (a + b)/2 for a, b in zip(X[j : -j - 1], X[j + 1 : -j])
325 ]
327 knots = [X[0]] * (d + 1) + list(interior_knots) + [X[-1]] * (d + 1)
329 basis = bspline_basis_set(d, knots, x)
331 A = [[b.subs(x, v) for b in basis] for v in X]
333 coeff = linsolve((Matrix(A), Matrix(Y)), symbols("c0:{}".format(len(X)), cls=Dummy))
334 coeff = list(coeff)[0]
335 intervals = {c for b in basis for (e, c) in b.args if c != True}
337 # Sorting the intervals
338 # ival contains the end-points of each interval
339 ival = [_ivl(c, x) for c in intervals]
340 com = zip(ival, intervals)
341 com = sorted(com, key=lambda x: x[0])
342 intervals = [y for x, y in com]
344 basis_dicts = [{c: e for (e, c) in b.args} for b in basis]
345 spline = []
346 for i in intervals:
347 piece = sum(
348 [c * d.get(i, S.Zero) for (c, d) in zip(coeff, basis_dicts)], S.Zero
349 )
350 spline.append((piece, i))
351 return Piecewise(*spline)