Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1import functools 

2import operator 

3 

4import numpy as np 

5from numpy.core.multiarray import normalize_axis_index 

6from scipy.linalg import (get_lapack_funcs, LinAlgError, 

7 cholesky_banded, cho_solve_banded) 

8from . import _bspl 

9from . import _fitpack_impl 

10from . import _fitpack as _dierckx 

11from scipy._lib._util import prod 

12 

13__all__ = ["BSpline", "make_interp_spline", "make_lsq_spline"] 

14 

15 

16def _get_dtype(dtype): 

17 """Return np.complex128 for complex dtypes, np.float64 otherwise.""" 

18 if np.issubdtype(dtype, np.complexfloating): 

19 return np.complex_ 

20 else: 

21 return np.float_ 

22 

23 

24def _as_float_array(x, check_finite=False): 

25 """Convert the input into a C contiguous float array. 

26 

27 NB: Upcasts half- and single-precision floats to double precision. 

28 """ 

29 x = np.ascontiguousarray(x) 

30 dtyp = _get_dtype(x.dtype) 

31 x = x.astype(dtyp, copy=False) 

32 if check_finite and not np.isfinite(x).all(): 

33 raise ValueError("Array must not contain infs or nans.") 

34 return x 

35 

36 

37class BSpline(object): 

38 r"""Univariate spline in the B-spline basis. 

39 

40 .. math:: 

41 

42 S(x) = \sum_{j=0}^{n-1} c_j B_{j, k; t}(x) 

43 

44 where :math:`B_{j, k; t}` are B-spline basis functions of degree `k` 

45 and knots `t`. 

46 

47 Parameters 

48 ---------- 

49 t : ndarray, shape (n+k+1,) 

50 knots 

51 c : ndarray, shape (>=n, ...) 

52 spline coefficients 

53 k : int 

54 B-spline degree 

55 extrapolate : bool or 'periodic', optional 

56 whether to extrapolate beyond the base interval, ``t[k] .. t[n]``, 

57 or to return nans. 

58 If True, extrapolates the first and last polynomial pieces of b-spline 

59 functions active on the base interval. 

60 If 'periodic', periodic extrapolation is used. 

61 Default is True. 

62 axis : int, optional 

63 Interpolation axis. Default is zero. 

64 

65 Attributes 

66 ---------- 

67 t : ndarray 

68 knot vector 

69 c : ndarray 

70 spline coefficients 

71 k : int 

72 spline degree 

73 extrapolate : bool 

74 If True, extrapolates the first and last polynomial pieces of b-spline 

75 functions active on the base interval. 

76 axis : int 

77 Interpolation axis. 

78 tck : tuple 

79 A read-only equivalent of ``(self.t, self.c, self.k)`` 

80 

81 Methods 

82 ------- 

83 __call__ 

84 basis_element 

85 derivative 

86 antiderivative 

87 integrate 

88 construct_fast 

89 

90 Notes 

91 ----- 

92 B-spline basis elements are defined via 

93 

94 .. math:: 

95 

96 B_{i, 0}(x) = 1, \textrm{if $t_i \le x < t_{i+1}$, otherwise $0$,} 

97 

98 B_{i, k}(x) = \frac{x - t_i}{t_{i+k} - t_i} B_{i, k-1}(x) 

99 + \frac{t_{i+k+1} - x}{t_{i+k+1} - t_{i+1}} B_{i+1, k-1}(x) 

100 

101 **Implementation details** 

102 

103 - At least ``k+1`` coefficients are required for a spline of degree `k`, 

104 so that ``n >= k+1``. Additional coefficients, ``c[j]`` with 

105 ``j > n``, are ignored. 

106 

107 - B-spline basis elements of degree `k` form a partition of unity on the 

108 *base interval*, ``t[k] <= x <= t[n]``. 

109 

110 

111 Examples 

112 -------- 

113 

114 Translating the recursive definition of B-splines into Python code, we have: 

115 

116 >>> def B(x, k, i, t): 

117 ... if k == 0: 

118 ... return 1.0 if t[i] <= x < t[i+1] else 0.0 

119 ... if t[i+k] == t[i]: 

120 ... c1 = 0.0 

121 ... else: 

122 ... c1 = (x - t[i])/(t[i+k] - t[i]) * B(x, k-1, i, t) 

123 ... if t[i+k+1] == t[i+1]: 

124 ... c2 = 0.0 

125 ... else: 

126 ... c2 = (t[i+k+1] - x)/(t[i+k+1] - t[i+1]) * B(x, k-1, i+1, t) 

127 ... return c1 + c2 

128 

129 >>> def bspline(x, t, c, k): 

130 ... n = len(t) - k - 1 

131 ... assert (n >= k+1) and (len(c) >= n) 

132 ... return sum(c[i] * B(x, k, i, t) for i in range(n)) 

133 

134 Note that this is an inefficient (if straightforward) way to 

135 evaluate B-splines --- this spline class does it in an equivalent, 

136 but much more efficient way. 

137 

138 Here we construct a quadratic spline function on the base interval 

139 ``2 <= x <= 4`` and compare with the naive way of evaluating the spline: 

140 

141 >>> from scipy.interpolate import BSpline 

142 >>> k = 2 

143 >>> t = [0, 1, 2, 3, 4, 5, 6] 

144 >>> c = [-1, 2, 0, -1] 

145 >>> spl = BSpline(t, c, k) 

146 >>> spl(2.5) 

147 array(1.375) 

148 >>> bspline(2.5, t, c, k) 

149 1.375 

150 

151 Note that outside of the base interval results differ. This is because 

152 `BSpline` extrapolates the first and last polynomial pieces of B-spline 

153 functions active on the base interval. 

154 

155 >>> import matplotlib.pyplot as plt 

156 >>> fig, ax = plt.subplots() 

157 >>> xx = np.linspace(1.5, 4.5, 50) 

158 >>> ax.plot(xx, [bspline(x, t, c ,k) for x in xx], 'r-', lw=3, label='naive') 

159 >>> ax.plot(xx, spl(xx), 'b-', lw=4, alpha=0.7, label='BSpline') 

160 >>> ax.grid(True) 

161 >>> ax.legend(loc='best') 

162 >>> plt.show() 

163 

164 

165 References 

166 ---------- 

167 .. [1] Tom Lyche and Knut Morken, Spline methods, 

168 http://www.uio.no/studier/emner/matnat/ifi/INF-MAT5340/v05/undervisningsmateriale/ 

169 .. [2] Carl de Boor, A practical guide to splines, Springer, 2001. 

170 

171 """ 

172 def __init__(self, t, c, k, extrapolate=True, axis=0): 

173 super(BSpline, self).__init__() 

174 

175 self.k = operator.index(k) 

176 self.c = np.asarray(c) 

177 self.t = np.ascontiguousarray(t, dtype=np.float64) 

178 

179 if extrapolate == 'periodic': 

180 self.extrapolate = extrapolate 

181 else: 

182 self.extrapolate = bool(extrapolate) 

183 

184 n = self.t.shape[0] - self.k - 1 

185 

186 axis = normalize_axis_index(axis, self.c.ndim) 

187 

188 # Note that the normalized axis is stored in the object. 

189 self.axis = axis 

190 if axis != 0: 

191 # roll the interpolation axis to be the first one in self.c 

192 # More specifically, the target shape for self.c is (n, ...), 

193 # and axis !=0 means that we have c.shape (..., n, ...) 

194 # ^ 

195 # axis 

196 self.c = np.rollaxis(self.c, axis) 

197 

198 if k < 0: 

199 raise ValueError("Spline order cannot be negative.") 

200 if self.t.ndim != 1: 

201 raise ValueError("Knot vector must be one-dimensional.") 

202 if n < self.k + 1: 

203 raise ValueError("Need at least %d knots for degree %d" % 

204 (2*k + 2, k)) 

205 if (np.diff(self.t) < 0).any(): 

206 raise ValueError("Knots must be in a non-decreasing order.") 

207 if len(np.unique(self.t[k:n+1])) < 2: 

208 raise ValueError("Need at least two internal knots.") 

209 if not np.isfinite(self.t).all(): 

210 raise ValueError("Knots should not have nans or infs.") 

211 if self.c.ndim < 1: 

212 raise ValueError("Coefficients must be at least 1-dimensional.") 

213 if self.c.shape[0] < n: 

214 raise ValueError("Knots, coefficients and degree are inconsistent.") 

215 

216 dt = _get_dtype(self.c.dtype) 

217 self.c = np.ascontiguousarray(self.c, dtype=dt) 

218 

219 @classmethod 

220 def construct_fast(cls, t, c, k, extrapolate=True, axis=0): 

221 """Construct a spline without making checks. 

222 

223 Accepts same parameters as the regular constructor. Input arrays 

224 `t` and `c` must of correct shape and dtype. 

225 """ 

226 self = object.__new__(cls) 

227 self.t, self.c, self.k = t, c, k 

228 self.extrapolate = extrapolate 

229 self.axis = axis 

230 return self 

231 

232 @property 

233 def tck(self): 

234 """Equivalent to ``(self.t, self.c, self.k)`` (read-only). 

235 """ 

236 return self.t, self.c, self.k 

237 

238 @classmethod 

239 def basis_element(cls, t, extrapolate=True): 

240 """Return a B-spline basis element ``B(x | t[0], ..., t[k+1])``. 

241 

242 Parameters 

243 ---------- 

244 t : ndarray, shape (k+1,) 

245 internal knots 

246 extrapolate : bool or 'periodic', optional 

247 whether to extrapolate beyond the base interval, ``t[0] .. t[k+1]``, 

248 or to return nans. 

249 If 'periodic', periodic extrapolation is used. 

250 Default is True. 

251 

252 Returns 

253 ------- 

254 basis_element : callable 

255 A callable representing a B-spline basis element for the knot 

256 vector `t`. 

257 

258 Notes 

259 ----- 

260 The degree of the B-spline, `k`, is inferred from the length of `t` as 

261 ``len(t)-2``. The knot vector is constructed by appending and prepending 

262 ``k+1`` elements to internal knots `t`. 

263 

264 Examples 

265 -------- 

266 

267 Construct a cubic B-spline: 

268 

269 >>> from scipy.interpolate import BSpline 

270 >>> b = BSpline.basis_element([0, 1, 2, 3, 4]) 

271 >>> k = b.k 

272 >>> b.t[k:-k] 

273 array([ 0., 1., 2., 3., 4.]) 

274 >>> k 

275 3 

276 

277 Construct a quadratic B-spline on ``[0, 1, 1, 2]``, and compare 

278 to its explicit form: 

279 

280 >>> t = [-1, 0, 1, 1, 2] 

281 >>> b = BSpline.basis_element(t[1:]) 

282 >>> def f(x): 

283 ... return np.where(x < 1, x*x, (2. - x)**2) 

284 

285 >>> import matplotlib.pyplot as plt 

286 >>> fig, ax = plt.subplots() 

287 >>> x = np.linspace(0, 2, 51) 

288 >>> ax.plot(x, b(x), 'g', lw=3) 

289 >>> ax.plot(x, f(x), 'r', lw=8, alpha=0.4) 

290 >>> ax.grid(True) 

291 >>> plt.show() 

292 

293 """ 

294 k = len(t) - 2 

295 t = _as_float_array(t) 

296 t = np.r_[(t[0]-1,) * k, t, (t[-1]+1,) * k] 

297 c = np.zeros_like(t) 

298 c[k] = 1. 

299 return cls.construct_fast(t, c, k, extrapolate) 

300 

301 def __call__(self, x, nu=0, extrapolate=None): 

302 """ 

303 Evaluate a spline function. 

304 

305 Parameters 

306 ---------- 

307 x : array_like 

308 points to evaluate the spline at. 

309 nu: int, optional 

310 derivative to evaluate (default is 0). 

311 extrapolate : bool or 'periodic', optional 

312 whether to extrapolate based on the first and last intervals 

313 or return nans. If 'periodic', periodic extrapolation is used. 

314 Default is `self.extrapolate`. 

315 

316 Returns 

317 ------- 

318 y : array_like 

319 Shape is determined by replacing the interpolation axis 

320 in the coefficient array with the shape of `x`. 

321 

322 """ 

323 if extrapolate is None: 

324 extrapolate = self.extrapolate 

325 x = np.asarray(x) 

326 x_shape, x_ndim = x.shape, x.ndim 

327 x = np.ascontiguousarray(x.ravel(), dtype=np.float_) 

328 

329 # With periodic extrapolation we map x to the segment 

330 # [self.t[k], self.t[n]]. 

331 if extrapolate == 'periodic': 

332 n = self.t.size - self.k - 1 

333 x = self.t[self.k] + (x - self.t[self.k]) % (self.t[n] - 

334 self.t[self.k]) 

335 extrapolate = False 

336 

337 out = np.empty((len(x), prod(self.c.shape[1:])), dtype=self.c.dtype) 

338 self._ensure_c_contiguous() 

339 self._evaluate(x, nu, extrapolate, out) 

340 out = out.reshape(x_shape + self.c.shape[1:]) 

341 if self.axis != 0: 

342 # transpose to move the calculated values to the interpolation axis 

343 l = list(range(out.ndim)) 

344 l = l[x_ndim:x_ndim+self.axis] + l[:x_ndim] + l[x_ndim+self.axis:] 

345 out = out.transpose(l) 

346 return out 

347 

348 def _evaluate(self, xp, nu, extrapolate, out): 

349 _bspl.evaluate_spline(self.t, self.c.reshape(self.c.shape[0], -1), 

350 self.k, xp, nu, extrapolate, out) 

351 

352 def _ensure_c_contiguous(self): 

353 """ 

354 c and t may be modified by the user. The Cython code expects 

355 that they are C contiguous. 

356 

357 """ 

358 if not self.t.flags.c_contiguous: 

359 self.t = self.t.copy() 

360 if not self.c.flags.c_contiguous: 

361 self.c = self.c.copy() 

362 

363 def derivative(self, nu=1): 

364 """Return a B-spline representing the derivative. 

365 

366 Parameters 

367 ---------- 

368 nu : int, optional 

369 Derivative order. 

370 Default is 1. 

371 

372 Returns 

373 ------- 

374 b : BSpline object 

375 A new instance representing the derivative. 

376 

377 See Also 

378 -------- 

379 splder, splantider 

380 

381 """ 

382 c = self.c 

383 # pad the c array if needed 

384 ct = len(self.t) - len(c) 

385 if ct > 0: 

386 c = np.r_[c, np.zeros((ct,) + c.shape[1:])] 

387 tck = _fitpack_impl.splder((self.t, c, self.k), nu) 

388 return self.construct_fast(*tck, extrapolate=self.extrapolate, 

389 axis=self.axis) 

390 

391 def antiderivative(self, nu=1): 

392 """Return a B-spline representing the antiderivative. 

393 

394 Parameters 

395 ---------- 

396 nu : int, optional 

397 Antiderivative order. Default is 1. 

398 

399 Returns 

400 ------- 

401 b : BSpline object 

402 A new instance representing the antiderivative. 

403 

404 Notes 

405 ----- 

406 If antiderivative is computed and ``self.extrapolate='periodic'``, 

407 it will be set to False for the returned instance. This is done because 

408 the antiderivative is no longer periodic and its correct evaluation 

409 outside of the initially given x interval is difficult. 

410 

411 See Also 

412 -------- 

413 splder, splantider 

414 

415 """ 

416 c = self.c 

417 # pad the c array if needed 

418 ct = len(self.t) - len(c) 

419 if ct > 0: 

420 c = np.r_[c, np.zeros((ct,) + c.shape[1:])] 

421 tck = _fitpack_impl.splantider((self.t, c, self.k), nu) 

422 

423 if self.extrapolate == 'periodic': 

424 extrapolate = False 

425 else: 

426 extrapolate = self.extrapolate 

427 

428 return self.construct_fast(*tck, extrapolate=extrapolate, 

429 axis=self.axis) 

430 

431 def integrate(self, a, b, extrapolate=None): 

432 """Compute a definite integral of the spline. 

433 

434 Parameters 

435 ---------- 

436 a : float 

437 Lower limit of integration. 

438 b : float 

439 Upper limit of integration. 

440 extrapolate : bool or 'periodic', optional 

441 whether to extrapolate beyond the base interval, 

442 ``t[k] .. t[-k-1]``, or take the spline to be zero outside of the 

443 base interval. If 'periodic', periodic extrapolation is used. 

444 If None (default), use `self.extrapolate`. 

445 

446 Returns 

447 ------- 

448 I : array_like 

449 Definite integral of the spline over the interval ``[a, b]``. 

450 

451 Examples 

452 -------- 

453 Construct the linear spline ``x if x < 1 else 2 - x`` on the base 

454 interval :math:`[0, 2]`, and integrate it 

455 

456 >>> from scipy.interpolate import BSpline 

457 >>> b = BSpline.basis_element([0, 1, 2]) 

458 >>> b.integrate(0, 1) 

459 array(0.5) 

460 

461 If the integration limits are outside of the base interval, the result 

462 is controlled by the `extrapolate` parameter 

463 

464 >>> b.integrate(-1, 1) 

465 array(0.0) 

466 >>> b.integrate(-1, 1, extrapolate=False) 

467 array(0.5) 

468 

469 >>> import matplotlib.pyplot as plt 

470 >>> fig, ax = plt.subplots() 

471 >>> ax.grid(True) 

472 >>> ax.axvline(0, c='r', lw=5, alpha=0.5) # base interval 

473 >>> ax.axvline(2, c='r', lw=5, alpha=0.5) 

474 >>> xx = [-1, 1, 2] 

475 >>> ax.plot(xx, b(xx)) 

476 >>> plt.show() 

477 

478 """ 

479 if extrapolate is None: 

480 extrapolate = self.extrapolate 

481 

482 # Prepare self.t and self.c. 

483 self._ensure_c_contiguous() 

484 

485 # Swap integration bounds if needed. 

486 sign = 1 

487 if b < a: 

488 a, b = b, a 

489 sign = -1 

490 n = self.t.size - self.k - 1 

491 

492 if extrapolate != "periodic" and not extrapolate: 

493 # Shrink the integration interval, if needed. 

494 a = max(a, self.t[self.k]) 

495 b = min(b, self.t[n]) 

496 

497 if self.c.ndim == 1: 

498 # Fast path: use FITPACK's routine 

499 # (cf _fitpack_impl.splint). 

500 t, c, k = self.tck 

501 integral, wrk = _dierckx._splint(t, c, k, a, b) 

502 return integral * sign 

503 

504 out = np.empty((2, prod(self.c.shape[1:])), dtype=self.c.dtype) 

505 

506 # Compute the antiderivative. 

507 c = self.c 

508 ct = len(self.t) - len(c) 

509 if ct > 0: 

510 c = np.r_[c, np.zeros((ct,) + c.shape[1:])] 

511 ta, ca, ka = _fitpack_impl.splantider((self.t, c, self.k), 1) 

512 

513 if extrapolate == 'periodic': 

514 # Split the integral into the part over period (can be several 

515 # of them) and the remaining part. 

516 

517 ts, te = self.t[self.k], self.t[n] 

518 period = te - ts 

519 interval = b - a 

520 n_periods, left = divmod(interval, period) 

521 

522 if n_periods > 0: 

523 # Evaluate the difference of antiderivatives. 

524 x = np.asarray([ts, te], dtype=np.float_) 

525 _bspl.evaluate_spline(ta, ca.reshape(ca.shape[0], -1), 

526 ka, x, 0, False, out) 

527 integral = out[1] - out[0] 

528 integral *= n_periods 

529 else: 

530 integral = np.zeros((1, prod(self.c.shape[1:])), 

531 dtype=self.c.dtype) 

532 

533 # Map a to [ts, te], b is always a + left. 

534 a = ts + (a - ts) % period 

535 b = a + left 

536 

537 # If b <= te then we need to integrate over [a, b], otherwise 

538 # over [a, te] and from xs to what is remained. 

539 if b <= te: 

540 x = np.asarray([a, b], dtype=np.float_) 

541 _bspl.evaluate_spline(ta, ca.reshape(ca.shape[0], -1), 

542 ka, x, 0, False, out) 

543 integral += out[1] - out[0] 

544 else: 

545 x = np.asarray([a, te], dtype=np.float_) 

546 _bspl.evaluate_spline(ta, ca.reshape(ca.shape[0], -1), 

547 ka, x, 0, False, out) 

548 integral += out[1] - out[0] 

549 

550 x = np.asarray([ts, ts + b - te], dtype=np.float_) 

551 _bspl.evaluate_spline(ta, ca.reshape(ca.shape[0], -1), 

552 ka, x, 0, False, out) 

553 integral += out[1] - out[0] 

554 else: 

555 # Evaluate the difference of antiderivatives. 

556 x = np.asarray([a, b], dtype=np.float_) 

557 _bspl.evaluate_spline(ta, ca.reshape(ca.shape[0], -1), 

558 ka, x, 0, extrapolate, out) 

559 integral = out[1] - out[0] 

560 

561 integral *= sign 

562 return integral.reshape(ca.shape[1:]) 

563 

564 

565################################# 

566# Interpolating spline helpers # 

567################################# 

568 

569def _not_a_knot(x, k): 

570 """Given data x, construct the knot vector w/ not-a-knot BC. 

571 cf de Boor, XIII(12).""" 

572 x = np.asarray(x) 

573 if k % 2 != 1: 

574 raise ValueError("Odd degree for now only. Got %s." % k) 

575 

576 m = (k - 1) // 2 

577 t = x[m+1:-m-1] 

578 t = np.r_[(x[0],)*(k+1), t, (x[-1],)*(k+1)] 

579 return t 

580 

581 

582def _augknt(x, k): 

583 """Construct a knot vector appropriate for the order-k interpolation.""" 

584 return np.r_[(x[0],)*k, x, (x[-1],)*k] 

585 

586 

587def _convert_string_aliases(deriv, target_shape): 

588 if isinstance(deriv, str): 

589 if deriv == "clamped": 

590 deriv = [(1, np.zeros(target_shape))] 

591 elif deriv == "natural": 

592 deriv = [(2, np.zeros(target_shape))] 

593 else: 

594 raise ValueError("Unknown boundary condition : %s" % deriv) 

595 return deriv 

596 

597 

598def _process_deriv_spec(deriv): 

599 if deriv is not None: 

600 try: 

601 ords, vals = zip(*deriv) 

602 except TypeError: 

603 msg = ("Derivatives, `bc_type`, should be specified as a pair of " 

604 "iterables of pairs of (order, value).") 

605 raise ValueError(msg) 

606 else: 

607 ords, vals = [], [] 

608 return np.atleast_1d(ords, vals) 

609 

610 

611def make_interp_spline(x, y, k=3, t=None, bc_type=None, axis=0, 

612 check_finite=True): 

613 """Compute the (coefficients of) interpolating B-spline. 

614 

615 Parameters 

616 ---------- 

617 x : array_like, shape (n,) 

618 Abscissas. 

619 y : array_like, shape (n, ...) 

620 Ordinates. 

621 k : int, optional 

622 B-spline degree. Default is cubic, k=3. 

623 t : array_like, shape (nt + k + 1,), optional. 

624 Knots. 

625 The number of knots needs to agree with the number of datapoints and 

626 the number of derivatives at the edges. Specifically, ``nt - n`` must 

627 equal ``len(deriv_l) + len(deriv_r)``. 

628 bc_type : 2-tuple or None 

629 Boundary conditions. 

630 Default is None, which means choosing the boundary conditions 

631 automatically. Otherwise, it must be a length-two tuple where the first 

632 element sets the boundary conditions at ``x[0]`` and the second 

633 element sets the boundary conditions at ``x[-1]``. Each of these must 

634 be an iterable of pairs ``(order, value)`` which gives the values of 

635 derivatives of specified orders at the given edge of the interpolation 

636 interval. 

637 Alternatively, the following string aliases are recognized: 

638 

639 * ``"clamped"``: The first derivatives at the ends are zero. This is 

640 equivalent to ``bc_type=([(1, 0.0)], [(1, 0.0)])``. 

641 * ``"natural"``: The second derivatives at ends are zero. This is 

642 equivalent to ``bc_type=([(2, 0.0)], [(2, 0.0)])``. 

643 * ``"not-a-knot"`` (default): The first and second segments are the same 

644 polynomial. This is equivalent to having ``bc_type=None``. 

645 

646 axis : int, optional 

647 Interpolation axis. Default is 0. 

648 check_finite : bool, optional 

649 Whether to check that the input arrays contain only finite numbers. 

650 Disabling may give a performance gain, but may result in problems 

651 (crashes, non-termination) if the inputs do contain infinities or NaNs. 

652 Default is True. 

653 

654 Returns 

655 ------- 

656 b : a BSpline object of the degree ``k`` and with knots ``t``. 

657 

658 Examples 

659 -------- 

660 

661 Use cubic interpolation on Chebyshev nodes: 

662 

663 >>> def cheb_nodes(N): 

664 ... jj = 2.*np.arange(N) + 1 

665 ... x = np.cos(np.pi * jj / 2 / N)[::-1] 

666 ... return x 

667 

668 >>> x = cheb_nodes(20) 

669 >>> y = np.sqrt(1 - x**2) 

670 

671 >>> from scipy.interpolate import BSpline, make_interp_spline 

672 >>> b = make_interp_spline(x, y) 

673 >>> np.allclose(b(x), y) 

674 True 

675 

676 Note that the default is a cubic spline with a not-a-knot boundary condition 

677 

678 >>> b.k 

679 3 

680 

681 Here we use a 'natural' spline, with zero 2nd derivatives at edges: 

682 

683 >>> l, r = [(2, 0.0)], [(2, 0.0)] 

684 >>> b_n = make_interp_spline(x, y, bc_type=(l, r)) # or, bc_type="natural" 

685 >>> np.allclose(b_n(x), y) 

686 True 

687 >>> x0, x1 = x[0], x[-1] 

688 >>> np.allclose([b_n(x0, 2), b_n(x1, 2)], [0, 0]) 

689 True 

690 

691 Interpolation of parametric curves is also supported. As an example, we 

692 compute a discretization of a snail curve in polar coordinates 

693 

694 >>> phi = np.linspace(0, 2.*np.pi, 40) 

695 >>> r = 0.3 + np.cos(phi) 

696 >>> x, y = r*np.cos(phi), r*np.sin(phi) # convert to Cartesian coordinates 

697 

698 Build an interpolating curve, parameterizing it by the angle 

699 

700 >>> from scipy.interpolate import make_interp_spline 

701 >>> spl = make_interp_spline(phi, np.c_[x, y]) 

702 

703 Evaluate the interpolant on a finer grid (note that we transpose the result 

704 to unpack it into a pair of x- and y-arrays) 

705 

706 >>> phi_new = np.linspace(0, 2.*np.pi, 100) 

707 >>> x_new, y_new = spl(phi_new).T 

708 

709 Plot the result 

710 

711 >>> import matplotlib.pyplot as plt 

712 >>> plt.plot(x, y, 'o') 

713 >>> plt.plot(x_new, y_new, '-') 

714 >>> plt.show() 

715 

716 See Also 

717 -------- 

718 BSpline : base class representing the B-spline objects 

719 CubicSpline : a cubic spline in the polynomial basis 

720 make_lsq_spline : a similar factory function for spline fitting 

721 UnivariateSpline : a wrapper over FITPACK spline fitting routines 

722 splrep : a wrapper over FITPACK spline fitting routines 

723 

724 """ 

725 # convert string aliases for the boundary conditions 

726 if bc_type is None or bc_type == 'not-a-knot': 

727 deriv_l, deriv_r = None, None 

728 elif isinstance(bc_type, str): 

729 deriv_l, deriv_r = bc_type, bc_type 

730 else: 

731 try: 

732 deriv_l, deriv_r = bc_type 

733 except TypeError: 

734 raise ValueError("Unknown boundary condition: %s" % bc_type) 

735 

736 y = np.asarray(y) 

737 

738 axis = normalize_axis_index(axis, y.ndim) 

739 

740 # special-case k=0 right away 

741 if k == 0: 

742 if any(_ is not None for _ in (t, deriv_l, deriv_r)): 

743 raise ValueError("Too much info for k=0: t and bc_type can only " 

744 "be None.") 

745 x = _as_float_array(x, check_finite) 

746 t = np.r_[x, x[-1]] 

747 c = np.asarray(y) 

748 c = np.rollaxis(c, axis) 

749 c = np.ascontiguousarray(c, dtype=_get_dtype(c.dtype)) 

750 return BSpline.construct_fast(t, c, k, axis=axis) 

751 

752 # special-case k=1 (e.g., Lyche and Morken, Eq.(2.16)) 

753 if k == 1 and t is None: 

754 if not (deriv_l is None and deriv_r is None): 

755 raise ValueError("Too much info for k=1: bc_type can only be None.") 

756 x = _as_float_array(x, check_finite) 

757 t = np.r_[x[0], x, x[-1]] 

758 c = np.asarray(y) 

759 c = np.rollaxis(c, axis) 

760 c = np.ascontiguousarray(c, dtype=_get_dtype(c.dtype)) 

761 return BSpline.construct_fast(t, c, k, axis=axis) 

762 

763 x = _as_float_array(x, check_finite) 

764 y = _as_float_array(y, check_finite) 

765 k = operator.index(k) 

766 

767 # come up with a sensible knot vector, if needed 

768 if t is None: 

769 if deriv_l is None and deriv_r is None: 

770 if k == 2: 

771 # OK, it's a bit ad hoc: Greville sites + omit 

772 # 2nd and 2nd-to-last points, a la not-a-knot 

773 t = (x[1:] + x[:-1]) / 2. 

774 t = np.r_[(x[0],)*(k+1), 

775 t[1:-1], 

776 (x[-1],)*(k+1)] 

777 else: 

778 t = _not_a_knot(x, k) 

779 else: 

780 t = _augknt(x, k) 

781 

782 t = _as_float_array(t, check_finite) 

783 

784 y = np.rollaxis(y, axis) # now internally interp axis is zero 

785 

786 if x.ndim != 1 or np.any(x[1:] < x[:-1]): 

787 raise ValueError("Expect x to be a 1-D sorted array_like.") 

788 if np.any(x[1:] == x[:-1]): 

789 raise ValueError("Expect x to not have duplicates") 

790 if k < 0: 

791 raise ValueError("Expect non-negative k.") 

792 if t.ndim != 1 or np.any(t[1:] < t[:-1]): 

793 raise ValueError("Expect t to be a 1-D sorted array_like.") 

794 if x.size != y.shape[0]: 

795 raise ValueError('x and y are incompatible.') 

796 if t.size < x.size + k + 1: 

797 raise ValueError('Got %d knots, need at least %d.' % 

798 (t.size, x.size + k + 1)) 

799 if (x[0] < t[k]) or (x[-1] > t[-k]): 

800 raise ValueError('Out of bounds w/ x = %s.' % x) 

801 

802 # Here : deriv_l, r = [(nu, value), ...] 

803 deriv_l = _convert_string_aliases(deriv_l, y.shape[1:]) 

804 deriv_l_ords, deriv_l_vals = _process_deriv_spec(deriv_l) 

805 nleft = deriv_l_ords.shape[0] 

806 

807 deriv_r = _convert_string_aliases(deriv_r, y.shape[1:]) 

808 deriv_r_ords, deriv_r_vals = _process_deriv_spec(deriv_r) 

809 nright = deriv_r_ords.shape[0] 

810 

811 # have `n` conditions for `nt` coefficients; need nt-n derivatives 

812 n = x.size 

813 nt = t.size - k - 1 

814 

815 if nt - n != nleft + nright: 

816 raise ValueError("The number of derivatives at boundaries does not " 

817 "match: expected %s, got %s+%s" % (nt-n, nleft, nright)) 

818 

819 # set up the LHS: the collocation matrix + derivatives at boundaries 

820 kl = ku = k 

821 ab = np.zeros((2*kl + ku + 1, nt), dtype=np.float_, order='F') 

822 _bspl._colloc(x, t, k, ab, offset=nleft) 

823 if nleft > 0: 

824 _bspl._handle_lhs_derivatives(t, k, x[0], ab, kl, ku, deriv_l_ords) 

825 if nright > 0: 

826 _bspl._handle_lhs_derivatives(t, k, x[-1], ab, kl, ku, deriv_r_ords, 

827 offset=nt-nright) 

828 

829 # set up the RHS: values to interpolate (+ derivative values, if any) 

830 extradim = prod(y.shape[1:]) 

831 rhs = np.empty((nt, extradim), dtype=y.dtype) 

832 if nleft > 0: 

833 rhs[:nleft] = deriv_l_vals.reshape(-1, extradim) 

834 rhs[nleft:nt - nright] = y.reshape(-1, extradim) 

835 if nright > 0: 

836 rhs[nt - nright:] = deriv_r_vals.reshape(-1, extradim) 

837 

838 # solve Ab @ x = rhs; this is the relevant part of linalg.solve_banded 

839 if check_finite: 

840 ab, rhs = map(np.asarray_chkfinite, (ab, rhs)) 

841 gbsv, = get_lapack_funcs(('gbsv',), (ab, rhs)) 

842 lu, piv, c, info = gbsv(kl, ku, ab, rhs, 

843 overwrite_ab=True, overwrite_b=True) 

844 

845 if info > 0: 

846 raise LinAlgError("Collocation matix is singular.") 

847 elif info < 0: 

848 raise ValueError('illegal value in %d-th argument of internal gbsv' % -info) 

849 

850 c = np.ascontiguousarray(c.reshape((nt,) + y.shape[1:])) 

851 return BSpline.construct_fast(t, c, k, axis=axis) 

852 

853 

854def make_lsq_spline(x, y, t, k=3, w=None, axis=0, check_finite=True): 

855 r"""Compute the (coefficients of) an LSQ B-spline. 

856 

857 The result is a linear combination 

858 

859 .. math:: 

860 

861 S(x) = \sum_j c_j B_j(x; t) 

862 

863 of the B-spline basis elements, :math:`B_j(x; t)`, which minimizes 

864 

865 .. math:: 

866 

867 \sum_{j} \left( w_j \times (S(x_j) - y_j) \right)^2 

868 

869 Parameters 

870 ---------- 

871 x : array_like, shape (m,) 

872 Abscissas. 

873 y : array_like, shape (m, ...) 

874 Ordinates. 

875 t : array_like, shape (n + k + 1,). 

876 Knots. 

877 Knots and data points must satisfy Schoenberg-Whitney conditions. 

878 k : int, optional 

879 B-spline degree. Default is cubic, k=3. 

880 w : array_like, shape (n,), optional 

881 Weights for spline fitting. Must be positive. If ``None``, 

882 then weights are all equal. 

883 Default is ``None``. 

884 axis : int, optional 

885 Interpolation axis. Default is zero. 

886 check_finite : bool, optional 

887 Whether to check that the input arrays contain only finite numbers. 

888 Disabling may give a performance gain, but may result in problems 

889 (crashes, non-termination) if the inputs do contain infinities or NaNs. 

890 Default is True. 

891 

892 Returns 

893 ------- 

894 b : a BSpline object of the degree `k` with knots `t`. 

895 

896 Notes 

897 ----- 

898 

899 The number of data points must be larger than the spline degree `k`. 

900 

901 Knots `t` must satisfy the Schoenberg-Whitney conditions, 

902 i.e., there must be a subset of data points ``x[j]`` such that 

903 ``t[j] < x[j] < t[j+k+1]``, for ``j=0, 1,...,n-k-2``. 

904 

905 Examples 

906 -------- 

907 Generate some noisy data: 

908 

909 >>> x = np.linspace(-3, 3, 50) 

910 >>> y = np.exp(-x**2) + 0.1 * np.random.randn(50) 

911 

912 Now fit a smoothing cubic spline with a pre-defined internal knots. 

913 Here we make the knot vector (k+1)-regular by adding boundary knots: 

914 

915 >>> from scipy.interpolate import make_lsq_spline, BSpline 

916 >>> t = [-1, 0, 1] 

917 >>> k = 3 

918 >>> t = np.r_[(x[0],)*(k+1), 

919 ... t, 

920 ... (x[-1],)*(k+1)] 

921 >>> spl = make_lsq_spline(x, y, t, k) 

922 

923 For comparison, we also construct an interpolating spline for the same 

924 set of data: 

925 

926 >>> from scipy.interpolate import make_interp_spline 

927 >>> spl_i = make_interp_spline(x, y) 

928 

929 Plot both: 

930 

931 >>> import matplotlib.pyplot as plt 

932 >>> xs = np.linspace(-3, 3, 100) 

933 >>> plt.plot(x, y, 'ro', ms=5) 

934 >>> plt.plot(xs, spl(xs), 'g-', lw=3, label='LSQ spline') 

935 >>> plt.plot(xs, spl_i(xs), 'b-', lw=3, alpha=0.7, label='interp spline') 

936 >>> plt.legend(loc='best') 

937 >>> plt.show() 

938 

939 **NaN handling**: If the input arrays contain ``nan`` values, the result is 

940 not useful since the underlying spline fitting routines cannot deal with 

941 ``nan``. A workaround is to use zero weights for not-a-number data points: 

942 

943 >>> y[8] = np.nan 

944 >>> w = np.isnan(y) 

945 >>> y[w] = 0. 

946 >>> tck = make_lsq_spline(x, y, t, w=~w) 

947 

948 Notice the need to replace a ``nan`` by a numerical value (precise value 

949 does not matter as long as the corresponding weight is zero.) 

950 

951 See Also 

952 -------- 

953 BSpline : base class representing the B-spline objects 

954 make_interp_spline : a similar factory function for interpolating splines 

955 LSQUnivariateSpline : a FITPACK-based spline fitting routine 

956 splrep : a FITPACK-based fitting routine 

957 

958 """ 

959 x = _as_float_array(x, check_finite) 

960 y = _as_float_array(y, check_finite) 

961 t = _as_float_array(t, check_finite) 

962 if w is not None: 

963 w = _as_float_array(w, check_finite) 

964 else: 

965 w = np.ones_like(x) 

966 k = operator.index(k) 

967 

968 axis = normalize_axis_index(axis, y.ndim) 

969 

970 y = np.rollaxis(y, axis) # now internally interp axis is zero 

971 

972 if x.ndim != 1 or np.any(x[1:] - x[:-1] <= 0): 

973 raise ValueError("Expect x to be a 1-D sorted array_like.") 

974 if x.shape[0] < k+1: 

975 raise ValueError("Need more x points.") 

976 if k < 0: 

977 raise ValueError("Expect non-negative k.") 

978 if t.ndim != 1 or np.any(t[1:] - t[:-1] < 0): 

979 raise ValueError("Expect t to be a 1-D sorted array_like.") 

980 if x.size != y.shape[0]: 

981 raise ValueError('x & y are incompatible.') 

982 if k > 0 and np.any((x < t[k]) | (x > t[-k])): 

983 raise ValueError('Out of bounds w/ x = %s.' % x) 

984 if x.size != w.size: 

985 raise ValueError('Incompatible weights.') 

986 

987 # number of coefficients 

988 n = t.size - k - 1 

989 

990 # construct A.T @ A and rhs with A the collocation matrix, and 

991 # rhs = A.T @ y for solving the LSQ problem ``A.T @ A @ c = A.T @ y`` 

992 lower = True 

993 extradim = prod(y.shape[1:]) 

994 ab = np.zeros((k+1, n), dtype=np.float_, order='F') 

995 rhs = np.zeros((n, extradim), dtype=y.dtype, order='F') 

996 _bspl._norm_eq_lsq(x, t, k, 

997 y.reshape(-1, extradim), 

998 w, 

999 ab, rhs) 

1000 rhs = rhs.reshape((n,) + y.shape[1:]) 

1001 

1002 # have observation matrix & rhs, can solve the LSQ problem 

1003 cho_decomp = cholesky_banded(ab, overwrite_ab=True, lower=lower, 

1004 check_finite=check_finite) 

1005 c = cho_solve_banded((cho_decomp, lower), rhs, overwrite_b=True, 

1006 check_finite=check_finite) 

1007 

1008 c = np.ascontiguousarray(c) 

1009 return BSpline.construct_fast(t, c, k, axis=axis)