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

1"""Interpolation algorithms using piecewise cubic polynomials.""" 

2 

3import numpy as np 

4 

5from . import PPoly 

6from .polyint import _isscalar 

7from scipy.linalg import solve_banded, solve 

8 

9 

10__all__ = ["CubicHermiteSpline", "PchipInterpolator", "pchip_interpolate", 

11 "Akima1DInterpolator", "CubicSpline"] 

12 

13 

14def prepare_input(x, y, axis, dydx=None): 

15 """Prepare input for cubic spline interpolators. 

16 

17 All data are converted to numpy arrays and checked for correctness. 

18 Axes equal to `axis` of arrays `y` and `dydx` are rolled to be the 0th 

19 axis. The value of `axis` is converted to lie in 

20 [0, number of dimensions of `y`). 

21 """ 

22 

23 x, y = map(np.asarray, (x, y)) 

24 if np.issubdtype(x.dtype, np.complexfloating): 

25 raise ValueError("`x` must contain real values.") 

26 x = x.astype(float) 

27 

28 if np.issubdtype(y.dtype, np.complexfloating): 

29 dtype = complex 

30 else: 

31 dtype = float 

32 

33 if dydx is not None: 

34 dydx = np.asarray(dydx) 

35 if y.shape != dydx.shape: 

36 raise ValueError("The shapes of `y` and `dydx` must be identical.") 

37 if np.issubdtype(dydx.dtype, np.complexfloating): 

38 dtype = complex 

39 dydx = dydx.astype(dtype, copy=False) 

40 

41 y = y.astype(dtype, copy=False) 

42 axis = axis % y.ndim 

43 if x.ndim != 1: 

44 raise ValueError("`x` must be 1-dimensional.") 

45 if x.shape[0] < 2: 

46 raise ValueError("`x` must contain at least 2 elements.") 

47 if x.shape[0] != y.shape[axis]: 

48 raise ValueError("The length of `y` along `axis`={0} doesn't " 

49 "match the length of `x`".format(axis)) 

50 

51 if not np.all(np.isfinite(x)): 

52 raise ValueError("`x` must contain only finite values.") 

53 if not np.all(np.isfinite(y)): 

54 raise ValueError("`y` must contain only finite values.") 

55 

56 if dydx is not None and not np.all(np.isfinite(dydx)): 

57 raise ValueError("`dydx` must contain only finite values.") 

58 

59 dx = np.diff(x) 

60 if np.any(dx <= 0): 

61 raise ValueError("`x` must be strictly increasing sequence.") 

62 

63 y = np.rollaxis(y, axis) 

64 if dydx is not None: 

65 dydx = np.rollaxis(dydx, axis) 

66 

67 return x, dx, y, axis, dydx 

68 

69 

70class CubicHermiteSpline(PPoly): 

71 """Piecewise-cubic interpolator matching values and first derivatives. 

72 

73 The result is represented as a `PPoly` instance. 

74 

75 Parameters 

76 ---------- 

77 x : array_like, shape (n,) 

78 1-D array containing values of the independent variable. 

79 Values must be real, finite and in strictly increasing order. 

80 y : array_like 

81 Array containing values of the dependent variable. It can have 

82 arbitrary number of dimensions, but the length along ``axis`` 

83 (see below) must match the length of ``x``. Values must be finite. 

84 dydx : array_like 

85 Array containing derivatives of the dependent variable. It can have 

86 arbitrary number of dimensions, but the length along ``axis`` 

87 (see below) must match the length of ``x``. Values must be finite. 

88 axis : int, optional 

89 Axis along which `y` is assumed to be varying. Meaning that for 

90 ``x[i]`` the corresponding values are ``np.take(y, i, axis=axis)``. 

91 Default is 0. 

92 extrapolate : {bool, 'periodic', None}, optional 

93 If bool, determines whether to extrapolate to out-of-bounds points 

94 based on first and last intervals, or to return NaNs. If 'periodic', 

95 periodic extrapolation is used. If None (default), it is set to True. 

96 

97 Attributes 

98 ---------- 

99 x : ndarray, shape (n,) 

100 Breakpoints. The same ``x`` which was passed to the constructor. 

101 c : ndarray, shape (4, n-1, ...) 

102 Coefficients of the polynomials on each segment. The trailing 

103 dimensions match the dimensions of `y`, excluding ``axis``. 

104 For example, if `y` is 1-D, then ``c[k, i]`` is a coefficient for 

105 ``(x-x[i])**(3-k)`` on the segment between ``x[i]`` and ``x[i+1]``. 

106 axis : int 

107 Interpolation axis. The same axis which was passed to the 

108 constructor. 

109 

110 Methods 

111 ------- 

112 __call__ 

113 derivative 

114 antiderivative 

115 integrate 

116 roots 

117 

118 See Also 

119 -------- 

120 Akima1DInterpolator : Akima 1D interpolator. 

121 PchipInterpolator : PCHIP 1-D monotonic cubic interpolator. 

122 CubicSpline : Cubic spline data interpolator. 

123 PPoly : Piecewise polynomial in terms of coefficients and breakpoints 

124 

125 Notes 

126 ----- 

127 If you want to create a higher-order spline matching higher-order 

128 derivatives, use `BPoly.from_derivatives`. 

129 

130 References 

131 ---------- 

132 .. [1] `Cubic Hermite spline 

133 <https://en.wikipedia.org/wiki/Cubic_Hermite_spline>`_ 

134 on Wikipedia. 

135 """ 

136 def __init__(self, x, y, dydx, axis=0, extrapolate=None): 

137 if extrapolate is None: 

138 extrapolate = True 

139 

140 x, dx, y, axis, dydx = prepare_input(x, y, axis, dydx) 

141 

142 dxr = dx.reshape([dx.shape[0]] + [1] * (y.ndim - 1)) 

143 slope = np.diff(y, axis=0) / dxr 

144 t = (dydx[:-1] + dydx[1:] - 2 * slope) / dxr 

145 

146 c = np.empty((4, len(x) - 1) + y.shape[1:], dtype=t.dtype) 

147 c[0] = t / dxr 

148 c[1] = (slope - dydx[:-1]) / dxr - t 

149 c[2] = dydx[:-1] 

150 c[3] = y[:-1] 

151 

152 super(CubicHermiteSpline, self).__init__(c, x, extrapolate=extrapolate) 

153 self.axis = axis 

154 

155 

156class PchipInterpolator(CubicHermiteSpline): 

157 r"""PCHIP 1-D monotonic cubic interpolation. 

158 

159 ``x`` and ``y`` are arrays of values used to approximate some function f, 

160 with ``y = f(x)``. The interpolant uses monotonic cubic splines 

161 to find the value of new points. (PCHIP stands for Piecewise Cubic 

162 Hermite Interpolating Polynomial). 

163 

164 Parameters 

165 ---------- 

166 x : ndarray 

167 A 1-D array of monotonically increasing real values. ``x`` cannot 

168 include duplicate values (otherwise f is overspecified) 

169 y : ndarray 

170 A 1-D array of real values. ``y``'s length along the interpolation 

171 axis must be equal to the length of ``x``. If N-D array, use ``axis`` 

172 parameter to select correct axis. 

173 axis : int, optional 

174 Axis in the y array corresponding to the x-coordinate values. 

175 extrapolate : bool, optional 

176 Whether to extrapolate to out-of-bounds points based on first 

177 and last intervals, or to return NaNs. 

178 

179 Methods 

180 ------- 

181 __call__ 

182 derivative 

183 antiderivative 

184 roots 

185 

186 See Also 

187 -------- 

188 CubicHermiteSpline : Piecewise-cubic interpolator. 

189 Akima1DInterpolator : Akima 1D interpolator. 

190 CubicSpline : Cubic spline data interpolator. 

191 PPoly : Piecewise polynomial in terms of coefficients and breakpoints. 

192 

193 Notes 

194 ----- 

195 The interpolator preserves monotonicity in the interpolation data and does 

196 not overshoot if the data is not smooth. 

197 

198 The first derivatives are guaranteed to be continuous, but the second 

199 derivatives may jump at :math:`x_k`. 

200 

201 Determines the derivatives at the points :math:`x_k`, :math:`f'_k`, 

202 by using PCHIP algorithm [1]_. 

203 

204 Let :math:`h_k = x_{k+1} - x_k`, and :math:`d_k = (y_{k+1} - y_k) / h_k` 

205 are the slopes at internal points :math:`x_k`. 

206 If the signs of :math:`d_k` and :math:`d_{k-1}` are different or either of 

207 them equals zero, then :math:`f'_k = 0`. Otherwise, it is given by the 

208 weighted harmonic mean 

209 

210 .. math:: 

211 

212 \frac{w_1 + w_2}{f'_k} = \frac{w_1}{d_{k-1}} + \frac{w_2}{d_k} 

213 

214 where :math:`w_1 = 2 h_k + h_{k-1}` and :math:`w_2 = h_k + 2 h_{k-1}`. 

215 

216 The end slopes are set using a one-sided scheme [2]_. 

217 

218 

219 References 

220 ---------- 

221 .. [1] F. N. Fritsch and R. E. Carlson, Monotone Piecewise Cubic Interpolation, 

222 SIAM J. Numer. Anal., 17(2), 238 (1980). 

223 :doi:`10.1137/0717021`. 

224 .. [2] see, e.g., C. Moler, Numerical Computing with Matlab, 2004. 

225 :doi:`10.1137/1.9780898717952` 

226 

227 

228 """ 

229 def __init__(self, x, y, axis=0, extrapolate=None): 

230 x, _, y, axis, _ = prepare_input(x, y, axis) 

231 xp = x.reshape((x.shape[0],) + (1,)*(y.ndim-1)) 

232 dk = self._find_derivatives(xp, y) 

233 super(PchipInterpolator, self).__init__(x, y, dk, axis=0, 

234 extrapolate=extrapolate) 

235 self.axis = axis 

236 

237 @staticmethod 

238 def _edge_case(h0, h1, m0, m1): 

239 # one-sided three-point estimate for the derivative 

240 d = ((2*h0 + h1)*m0 - h0*m1) / (h0 + h1) 

241 

242 # try to preserve shape 

243 mask = np.sign(d) != np.sign(m0) 

244 mask2 = (np.sign(m0) != np.sign(m1)) & (np.abs(d) > 3.*np.abs(m0)) 

245 mmm = (~mask) & mask2 

246 

247 d[mask] = 0. 

248 d[mmm] = 3.*m0[mmm] 

249 

250 return d 

251 

252 @staticmethod 

253 def _find_derivatives(x, y): 

254 # Determine the derivatives at the points y_k, d_k, by using 

255 # PCHIP algorithm is: 

256 # We choose the derivatives at the point x_k by 

257 # Let m_k be the slope of the kth segment (between k and k+1) 

258 # If m_k=0 or m_{k-1}=0 or sgn(m_k) != sgn(m_{k-1}) then d_k == 0 

259 # else use weighted harmonic mean: 

260 # w_1 = 2h_k + h_{k-1}, w_2 = h_k + 2h_{k-1} 

261 # 1/d_k = 1/(w_1 + w_2)*(w_1 / m_k + w_2 / m_{k-1}) 

262 # where h_k is the spacing between x_k and x_{k+1} 

263 y_shape = y.shape 

264 if y.ndim == 1: 

265 # So that _edge_case doesn't end up assigning to scalars 

266 x = x[:, None] 

267 y = y[:, None] 

268 

269 hk = x[1:] - x[:-1] 

270 mk = (y[1:] - y[:-1]) / hk 

271 

272 if y.shape[0] == 2: 

273 # edge case: only have two points, use linear interpolation 

274 dk = np.zeros_like(y) 

275 dk[0] = mk 

276 dk[1] = mk 

277 return dk.reshape(y_shape) 

278 

279 smk = np.sign(mk) 

280 condition = (smk[1:] != smk[:-1]) | (mk[1:] == 0) | (mk[:-1] == 0) 

281 

282 w1 = 2*hk[1:] + hk[:-1] 

283 w2 = hk[1:] + 2*hk[:-1] 

284 

285 # values where division by zero occurs will be excluded 

286 # by 'condition' afterwards 

287 with np.errstate(divide='ignore'): 

288 whmean = (w1/mk[:-1] + w2/mk[1:]) / (w1 + w2) 

289 

290 dk = np.zeros_like(y) 

291 dk[1:-1][condition] = 0.0 

292 dk[1:-1][~condition] = 1.0 / whmean[~condition] 

293 

294 # special case endpoints, as suggested in 

295 # Cleve Moler, Numerical Computing with MATLAB, Chap 3.4 

296 dk[0] = PchipInterpolator._edge_case(hk[0], hk[1], mk[0], mk[1]) 

297 dk[-1] = PchipInterpolator._edge_case(hk[-1], hk[-2], mk[-1], mk[-2]) 

298 

299 return dk.reshape(y_shape) 

300 

301 

302def pchip_interpolate(xi, yi, x, der=0, axis=0): 

303 """ 

304 Convenience function for pchip interpolation. 

305 

306 xi and yi are arrays of values used to approximate some function f, 

307 with ``yi = f(xi)``. The interpolant uses monotonic cubic splines 

308 to find the value of new points x and the derivatives there. 

309 

310 See `scipy.interpolate.PchipInterpolator` for details. 

311 

312 Parameters 

313 ---------- 

314 xi : array_like 

315 A sorted list of x-coordinates, of length N. 

316 yi : array_like 

317 A 1-D array of real values. `yi`'s length along the interpolation 

318 axis must be equal to the length of `xi`. If N-D array, use axis 

319 parameter to select correct axis. 

320 x : scalar or array_like 

321 Of length M. 

322 der : int or list, optional 

323 Derivatives to extract. The 0th derivative can be included to 

324 return the function value. 

325 axis : int, optional 

326 Axis in the yi array corresponding to the x-coordinate values. 

327 

328 See Also 

329 -------- 

330 PchipInterpolator : PCHIP 1-D monotonic cubic interpolator. 

331 

332 Returns 

333 ------- 

334 y : scalar or array_like 

335 The result, of length R or length M or M by R, 

336 

337 Examples 

338 -------- 

339 We can interpolate 2D observed data using pchip interpolation: 

340 

341 >>> import matplotlib.pyplot as plt 

342 >>> from scipy.interpolate import pchip_interpolate 

343 >>> x_observed = np.linspace(0.0, 10.0, 11) 

344 >>> y_observed = np.sin(x_observed) 

345 >>> x = np.linspace(min(x_observed), max(x_observed), num=100) 

346 >>> y = pchip_interpolate(x_observed, y_observed, x) 

347 >>> plt.plot(x_observed, y_observed, "o", label="observation") 

348 >>> plt.plot(x, y, label="pchip interpolation") 

349 >>> plt.legend() 

350 >>> plt.show() 

351 

352 """ 

353 P = PchipInterpolator(xi, yi, axis=axis) 

354 

355 if der == 0: 

356 return P(x) 

357 elif _isscalar(der): 

358 return P.derivative(der)(x) 

359 else: 

360 return [P.derivative(nu)(x) for nu in der] 

361 

362 

363class Akima1DInterpolator(CubicHermiteSpline): 

364 """ 

365 Akima interpolator 

366 

367 Fit piecewise cubic polynomials, given vectors x and y. The interpolation 

368 method by Akima uses a continuously differentiable sub-spline built from 

369 piecewise cubic polynomials. The resultant curve passes through the given 

370 data points and will appear smooth and natural. 

371 

372 Parameters 

373 ---------- 

374 x : ndarray, shape (m, ) 

375 1-D array of monotonically increasing real values. 

376 y : ndarray, shape (m, ...) 

377 N-D array of real values. The length of ``y`` along the first axis 

378 must be equal to the length of ``x``. 

379 axis : int, optional 

380 Specifies the axis of ``y`` along which to interpolate. Interpolation 

381 defaults to the first axis of ``y``. 

382 

383 Methods 

384 ------- 

385 __call__ 

386 derivative 

387 antiderivative 

388 roots 

389 

390 See Also 

391 -------- 

392 PchipInterpolator : PCHIP 1-D monotonic cubic interpolator. 

393 CubicSpline : Cubic spline data interpolator. 

394 PPoly : Piecewise polynomial in terms of coefficients and breakpoints 

395 

396 Notes 

397 ----- 

398 .. versionadded:: 0.14 

399 

400 Use only for precise data, as the fitted curve passes through the given 

401 points exactly. This routine is useful for plotting a pleasingly smooth 

402 curve through a few given points for purposes of plotting. 

403 

404 References 

405 ---------- 

406 [1] A new method of interpolation and smooth curve fitting based 

407 on local procedures. Hiroshi Akima, J. ACM, October 1970, 17(4), 

408 589-602. 

409 

410 """ 

411 

412 def __init__(self, x, y, axis=0): 

413 # Original implementation in MATLAB by N. Shamsundar (BSD licensed), see 

414 # https://www.mathworks.com/matlabcentral/fileexchange/1814-akima-interpolation 

415 x, dx, y, axis, _ = prepare_input(x, y, axis) 

416 # determine slopes between breakpoints 

417 m = np.empty((x.size + 3, ) + y.shape[1:]) 

418 dx = dx[(slice(None), ) + (None, ) * (y.ndim - 1)] 

419 m[2:-2] = np.diff(y, axis=0) / dx 

420 

421 # add two additional points on the left ... 

422 m[1] = 2. * m[2] - m[3] 

423 m[0] = 2. * m[1] - m[2] 

424 # ... and on the right 

425 m[-2] = 2. * m[-3] - m[-4] 

426 m[-1] = 2. * m[-2] - m[-3] 

427 

428 # if m1 == m2 != m3 == m4, the slope at the breakpoint is not defined. 

429 # This is the fill value: 

430 t = .5 * (m[3:] + m[:-3]) 

431 # get the denominator of the slope t 

432 dm = np.abs(np.diff(m, axis=0)) 

433 f1 = dm[2:] 

434 f2 = dm[:-2] 

435 f12 = f1 + f2 

436 # These are the mask of where the the slope at breakpoint is defined: 

437 ind = np.nonzero(f12 > 1e-9 * np.max(f12)) 

438 x_ind, y_ind = ind[0], ind[1:] 

439 # Set the slope at breakpoint 

440 t[ind] = (f1[ind] * m[(x_ind + 1,) + y_ind] + 

441 f2[ind] * m[(x_ind + 2,) + y_ind]) / f12[ind] 

442 

443 super(Akima1DInterpolator, self).__init__(x, y, t, axis=0, 

444 extrapolate=False) 

445 self.axis = axis 

446 

447 def extend(self, c, x, right=True): 

448 raise NotImplementedError("Extending a 1-D Akima interpolator is not " 

449 "yet implemented") 

450 

451 # These are inherited from PPoly, but they do not produce an Akima 

452 # interpolator. Hence stub them out. 

453 @classmethod 

454 def from_spline(cls, tck, extrapolate=None): 

455 raise NotImplementedError("This method does not make sense for " 

456 "an Akima interpolator.") 

457 

458 @classmethod 

459 def from_bernstein_basis(cls, bp, extrapolate=None): 

460 raise NotImplementedError("This method does not make sense for " 

461 "an Akima interpolator.") 

462 

463 

464class CubicSpline(CubicHermiteSpline): 

465 """Cubic spline data interpolator. 

466 

467 Interpolate data with a piecewise cubic polynomial which is twice 

468 continuously differentiable [1]_. The result is represented as a `PPoly` 

469 instance with breakpoints matching the given data. 

470 

471 Parameters 

472 ---------- 

473 x : array_like, shape (n,) 

474 1-D array containing values of the independent variable. 

475 Values must be real, finite and in strictly increasing order. 

476 y : array_like 

477 Array containing values of the dependent variable. It can have 

478 arbitrary number of dimensions, but the length along ``axis`` 

479 (see below) must match the length of ``x``. Values must be finite. 

480 axis : int, optional 

481 Axis along which `y` is assumed to be varying. Meaning that for 

482 ``x[i]`` the corresponding values are ``np.take(y, i, axis=axis)``. 

483 Default is 0. 

484 bc_type : string or 2-tuple, optional 

485 Boundary condition type. Two additional equations, given by the 

486 boundary conditions, are required to determine all coefficients of 

487 polynomials on each segment [2]_. 

488 

489 If `bc_type` is a string, then the specified condition will be applied 

490 at both ends of a spline. Available conditions are: 

491 

492 * 'not-a-knot' (default): The first and second segment at a curve end 

493 are the same polynomial. It is a good default when there is no 

494 information on boundary conditions. 

495 * 'periodic': The interpolated functions is assumed to be periodic 

496 of period ``x[-1] - x[0]``. The first and last value of `y` must be 

497 identical: ``y[0] == y[-1]``. This boundary condition will result in 

498 ``y'[0] == y'[-1]`` and ``y''[0] == y''[-1]``. 

499 * 'clamped': The first derivative at curves ends are zero. Assuming 

500 a 1D `y`, ``bc_type=((1, 0.0), (1, 0.0))`` is the same condition. 

501 * 'natural': The second derivative at curve ends are zero. Assuming 

502 a 1D `y`, ``bc_type=((2, 0.0), (2, 0.0))`` is the same condition. 

503 

504 If `bc_type` is a 2-tuple, the first and the second value will be 

505 applied at the curve start and end respectively. The tuple values can 

506 be one of the previously mentioned strings (except 'periodic') or a 

507 tuple `(order, deriv_values)` allowing to specify arbitrary 

508 derivatives at curve ends: 

509 

510 * `order`: the derivative order, 1 or 2. 

511 * `deriv_value`: array_like containing derivative values, shape must 

512 be the same as `y`, excluding ``axis`` dimension. For example, if 

513 `y` is 1-D, then `deriv_value` must be a scalar. If `y` is 3-D with 

514 the shape (n0, n1, n2) and axis=2, then `deriv_value` must be 2-D 

515 and have the shape (n0, n1). 

516 extrapolate : {bool, 'periodic', None}, optional 

517 If bool, determines whether to extrapolate to out-of-bounds points 

518 based on first and last intervals, or to return NaNs. If 'periodic', 

519 periodic extrapolation is used. If None (default), ``extrapolate`` is 

520 set to 'periodic' for ``bc_type='periodic'`` and to True otherwise. 

521 

522 Attributes 

523 ---------- 

524 x : ndarray, shape (n,) 

525 Breakpoints. The same ``x`` which was passed to the constructor. 

526 c : ndarray, shape (4, n-1, ...) 

527 Coefficients of the polynomials on each segment. The trailing 

528 dimensions match the dimensions of `y`, excluding ``axis``. 

529 For example, if `y` is 1-d, then ``c[k, i]`` is a coefficient for 

530 ``(x-x[i])**(3-k)`` on the segment between ``x[i]`` and ``x[i+1]``. 

531 axis : int 

532 Interpolation axis. The same axis which was passed to the 

533 constructor. 

534 

535 Methods 

536 ------- 

537 __call__ 

538 derivative 

539 antiderivative 

540 integrate 

541 roots 

542 

543 See Also 

544 -------- 

545 Akima1DInterpolator : Akima 1D interpolator. 

546 PchipInterpolator : PCHIP 1-D monotonic cubic interpolator. 

547 PPoly : Piecewise polynomial in terms of coefficients and breakpoints. 

548 

549 Notes 

550 ----- 

551 Parameters `bc_type` and ``interpolate`` work independently, i.e. the 

552 former controls only construction of a spline, and the latter only 

553 evaluation. 

554 

555 When a boundary condition is 'not-a-knot' and n = 2, it is replaced by 

556 a condition that the first derivative is equal to the linear interpolant 

557 slope. When both boundary conditions are 'not-a-knot' and n = 3, the 

558 solution is sought as a parabola passing through given points. 

559 

560 When 'not-a-knot' boundary conditions is applied to both ends, the 

561 resulting spline will be the same as returned by `splrep` (with ``s=0``) 

562 and `InterpolatedUnivariateSpline`, but these two methods use a 

563 representation in B-spline basis. 

564 

565 .. versionadded:: 0.18.0 

566 

567 Examples 

568 -------- 

569 In this example the cubic spline is used to interpolate a sampled sinusoid. 

570 You can see that the spline continuity property holds for the first and 

571 second derivatives and violates only for the third derivative. 

572 

573 >>> from scipy.interpolate import CubicSpline 

574 >>> import matplotlib.pyplot as plt 

575 >>> x = np.arange(10) 

576 >>> y = np.sin(x) 

577 >>> cs = CubicSpline(x, y) 

578 >>> xs = np.arange(-0.5, 9.6, 0.1) 

579 >>> fig, ax = plt.subplots(figsize=(6.5, 4)) 

580 >>> ax.plot(x, y, 'o', label='data') 

581 >>> ax.plot(xs, np.sin(xs), label='true') 

582 >>> ax.plot(xs, cs(xs), label="S") 

583 >>> ax.plot(xs, cs(xs, 1), label="S'") 

584 >>> ax.plot(xs, cs(xs, 2), label="S''") 

585 >>> ax.plot(xs, cs(xs, 3), label="S'''") 

586 >>> ax.set_xlim(-0.5, 9.5) 

587 >>> ax.legend(loc='lower left', ncol=2) 

588 >>> plt.show() 

589 

590 In the second example, the unit circle is interpolated with a spline. A 

591 periodic boundary condition is used. You can see that the first derivative 

592 values, ds/dx=0, ds/dy=1 at the periodic point (1, 0) are correctly 

593 computed. Note that a circle cannot be exactly represented by a cubic 

594 spline. To increase precision, more breakpoints would be required. 

595 

596 >>> theta = 2 * np.pi * np.linspace(0, 1, 5) 

597 >>> y = np.c_[np.cos(theta), np.sin(theta)] 

598 >>> cs = CubicSpline(theta, y, bc_type='periodic') 

599 >>> print("ds/dx={:.1f} ds/dy={:.1f}".format(cs(0, 1)[0], cs(0, 1)[1])) 

600 ds/dx=0.0 ds/dy=1.0 

601 >>> xs = 2 * np.pi * np.linspace(0, 1, 100) 

602 >>> fig, ax = plt.subplots(figsize=(6.5, 4)) 

603 >>> ax.plot(y[:, 0], y[:, 1], 'o', label='data') 

604 >>> ax.plot(np.cos(xs), np.sin(xs), label='true') 

605 >>> ax.plot(cs(xs)[:, 0], cs(xs)[:, 1], label='spline') 

606 >>> ax.axes.set_aspect('equal') 

607 >>> ax.legend(loc='center') 

608 >>> plt.show() 

609 

610 The third example is the interpolation of a polynomial y = x**3 on the 

611 interval 0 <= x<= 1. A cubic spline can represent this function exactly. 

612 To achieve that we need to specify values and first derivatives at 

613 endpoints of the interval. Note that y' = 3 * x**2 and thus y'(0) = 0 and 

614 y'(1) = 3. 

615 

616 >>> cs = CubicSpline([0, 1], [0, 1], bc_type=((1, 0), (1, 3))) 

617 >>> x = np.linspace(0, 1) 

618 >>> np.allclose(x**3, cs(x)) 

619 True 

620 

621 References 

622 ---------- 

623 .. [1] `Cubic Spline Interpolation 

624 <https://en.wikiversity.org/wiki/Cubic_Spline_Interpolation>`_ 

625 on Wikiversity. 

626 .. [2] Carl de Boor, "A Practical Guide to Splines", Springer-Verlag, 1978. 

627 """ 

628 def __init__(self, x, y, axis=0, bc_type='not-a-knot', extrapolate=None): 

629 x, dx, y, axis, _ = prepare_input(x, y, axis) 

630 n = len(x) 

631 

632 bc, y = self._validate_bc(bc_type, y, y.shape[1:], axis) 

633 

634 if extrapolate is None: 

635 if bc[0] == 'periodic': 

636 extrapolate = 'periodic' 

637 else: 

638 extrapolate = True 

639 

640 dxr = dx.reshape([dx.shape[0]] + [1] * (y.ndim - 1)) 

641 slope = np.diff(y, axis=0) / dxr 

642 

643 # If bc is 'not-a-knot' this change is just a convention. 

644 # If bc is 'periodic' then we already checked that y[0] == y[-1], 

645 # and the spline is just a constant, we handle this case in the same 

646 # way by setting the first derivatives to slope, which is 0. 

647 if n == 2: 

648 if bc[0] in ['not-a-knot', 'periodic']: 

649 bc[0] = (1, slope[0]) 

650 if bc[1] in ['not-a-knot', 'periodic']: 

651 bc[1] = (1, slope[0]) 

652 

653 # This is a very special case, when both conditions are 'not-a-knot' 

654 # and n == 3. In this case 'not-a-knot' can't be handled regularly 

655 # as the both conditions are identical. We handle this case by 

656 # constructing a parabola passing through given points. 

657 if n == 3 and bc[0] == 'not-a-knot' and bc[1] == 'not-a-knot': 

658 A = np.zeros((3, 3)) # This is a standard matrix. 

659 b = np.empty((3,) + y.shape[1:], dtype=y.dtype) 

660 

661 A[0, 0] = 1 

662 A[0, 1] = 1 

663 A[1, 0] = dx[1] 

664 A[1, 1] = 2 * (dx[0] + dx[1]) 

665 A[1, 2] = dx[0] 

666 A[2, 1] = 1 

667 A[2, 2] = 1 

668 

669 b[0] = 2 * slope[0] 

670 b[1] = 3 * (dxr[0] * slope[1] + dxr[1] * slope[0]) 

671 b[2] = 2 * slope[1] 

672 

673 s = solve(A, b, overwrite_a=True, overwrite_b=True, 

674 check_finite=False) 

675 else: 

676 # Find derivative values at each x[i] by solving a tridiagonal 

677 # system. 

678 A = np.zeros((3, n)) # This is a banded matrix representation. 

679 b = np.empty((n,) + y.shape[1:], dtype=y.dtype) 

680 

681 # Filling the system for i=1..n-2 

682 # (x[i-1] - x[i]) * s[i-1] +\ 

683 # 2 * ((x[i] - x[i-1]) + (x[i+1] - x[i])) * s[i] +\ 

684 # (x[i] - x[i-1]) * s[i+1] =\ 

685 # 3 * ((x[i+1] - x[i])*(y[i] - y[i-1])/(x[i] - x[i-1]) +\ 

686 # (x[i] - x[i-1])*(y[i+1] - y[i])/(x[i+1] - x[i])) 

687 

688 A[1, 1:-1] = 2 * (dx[:-1] + dx[1:]) # The diagonal 

689 A[0, 2:] = dx[:-1] # The upper diagonal 

690 A[-1, :-2] = dx[1:] # The lower diagonal 

691 

692 b[1:-1] = 3 * (dxr[1:] * slope[:-1] + dxr[:-1] * slope[1:]) 

693 

694 bc_start, bc_end = bc 

695 

696 if bc_start == 'periodic': 

697 # Due to the periodicity, and because y[-1] = y[0], the linear 

698 # system has (n-1) unknowns/equations instead of n: 

699 A = A[:, 0:-1] 

700 A[1, 0] = 2 * (dx[-1] + dx[0]) 

701 A[0, 1] = dx[-1] 

702 

703 b = b[:-1] 

704 

705 # Also, due to the periodicity, the system is not tri-diagonal. 

706 # We need to compute a "condensed" matrix of shape (n-2, n-2). 

707 # See https://web.archive.org/web/20151220180652/http://www.cfm.brown.edu/people/gk/chap6/node14.html 

708 # for more explanations. 

709 # The condensed matrix is obtained by removing the last column 

710 # and last row of the (n-1, n-1) system matrix. The removed 

711 # values are saved in scalar variables with the (n-1, n-1) 

712 # system matrix indices forming their names: 

713 a_m1_0 = dx[-2] # lower left corner value: A[-1, 0] 

714 a_m1_m2 = dx[-1] 

715 a_m1_m1 = 2 * (dx[-1] + dx[-2]) 

716 a_m2_m1 = dx[-2] 

717 a_0_m1 = dx[0] 

718 

719 b[0] = 3 * (dxr[0] * slope[-1] + dxr[-1] * slope[0]) 

720 b[-1] = 3 * (dxr[-1] * slope[-2] + dxr[-2] * slope[-1]) 

721 

722 Ac = A[:, :-1] 

723 b1 = b[:-1] 

724 b2 = np.zeros_like(b1) 

725 b2[0] = -a_0_m1 

726 b2[-1] = -a_m2_m1 

727 

728 # s1 and s2 are the solutions of (n-2, n-2) system 

729 s1 = solve_banded((1, 1), Ac, b1, overwrite_ab=False, 

730 overwrite_b=False, check_finite=False) 

731 

732 s2 = solve_banded((1, 1), Ac, b2, overwrite_ab=False, 

733 overwrite_b=False, check_finite=False) 

734 

735 # computing the s[n-2] solution: 

736 s_m1 = ((b[-1] - a_m1_0 * s1[0] - a_m1_m2 * s1[-1]) / 

737 (a_m1_m1 + a_m1_0 * s2[0] + a_m1_m2 * s2[-1])) 

738 

739 # s is the solution of the (n, n) system: 

740 s = np.empty((n,) + y.shape[1:], dtype=y.dtype) 

741 s[:-2] = s1 + s_m1 * s2 

742 s[-2] = s_m1 

743 s[-1] = s[0] 

744 else: 

745 if bc_start == 'not-a-knot': 

746 A[1, 0] = dx[1] 

747 A[0, 1] = x[2] - x[0] 

748 d = x[2] - x[0] 

749 b[0] = ((dxr[0] + 2*d) * dxr[1] * slope[0] + 

750 dxr[0]**2 * slope[1]) / d 

751 elif bc_start[0] == 1: 

752 A[1, 0] = 1 

753 A[0, 1] = 0 

754 b[0] = bc_start[1] 

755 elif bc_start[0] == 2: 

756 A[1, 0] = 2 * dx[0] 

757 A[0, 1] = dx[0] 

758 b[0] = -0.5 * bc_start[1] * dx[0]**2 + 3 * (y[1] - y[0]) 

759 

760 if bc_end == 'not-a-knot': 

761 A[1, -1] = dx[-2] 

762 A[-1, -2] = x[-1] - x[-3] 

763 d = x[-1] - x[-3] 

764 b[-1] = ((dxr[-1]**2*slope[-2] + 

765 (2*d + dxr[-1])*dxr[-2]*slope[-1]) / d) 

766 elif bc_end[0] == 1: 

767 A[1, -1] = 1 

768 A[-1, -2] = 0 

769 b[-1] = bc_end[1] 

770 elif bc_end[0] == 2: 

771 A[1, -1] = 2 * dx[-1] 

772 A[-1, -2] = dx[-1] 

773 b[-1] = 0.5 * bc_end[1] * dx[-1]**2 + 3 * (y[-1] - y[-2]) 

774 

775 s = solve_banded((1, 1), A, b, overwrite_ab=True, 

776 overwrite_b=True, check_finite=False) 

777 

778 super(CubicSpline, self).__init__(x, y, s, axis=0, 

779 extrapolate=extrapolate) 

780 self.axis = axis 

781 

782 @staticmethod 

783 def _validate_bc(bc_type, y, expected_deriv_shape, axis): 

784 """Validate and prepare boundary conditions. 

785 

786 Returns 

787 ------- 

788 validated_bc : 2-tuple 

789 Boundary conditions for a curve start and end. 

790 y : ndarray 

791 y casted to complex dtype if one of the boundary conditions has 

792 complex dtype. 

793 """ 

794 if isinstance(bc_type, str): 

795 if bc_type == 'periodic': 

796 if not np.allclose(y[0], y[-1], rtol=1e-15, atol=1e-15): 

797 raise ValueError( 

798 "The first and last `y` point along axis {} must " 

799 "be identical (within machine precision) when " 

800 "bc_type='periodic'.".format(axis)) 

801 

802 bc_type = (bc_type, bc_type) 

803 

804 else: 

805 if len(bc_type) != 2: 

806 raise ValueError("`bc_type` must contain 2 elements to " 

807 "specify start and end conditions.") 

808 

809 if 'periodic' in bc_type: 

810 raise ValueError("'periodic' `bc_type` is defined for both " 

811 "curve ends and cannot be used with other " 

812 "boundary conditions.") 

813 

814 validated_bc = [] 

815 for bc in bc_type: 

816 if isinstance(bc, str): 

817 if bc == 'clamped': 

818 validated_bc.append((1, np.zeros(expected_deriv_shape))) 

819 elif bc == 'natural': 

820 validated_bc.append((2, np.zeros(expected_deriv_shape))) 

821 elif bc in ['not-a-knot', 'periodic']: 

822 validated_bc.append(bc) 

823 else: 

824 raise ValueError("bc_type={} is not allowed.".format(bc)) 

825 else: 

826 try: 

827 deriv_order, deriv_value = bc 

828 except Exception: 

829 raise ValueError("A specified derivative value must be " 

830 "given in the form (order, value).") 

831 

832 if deriv_order not in [1, 2]: 

833 raise ValueError("The specified derivative order must " 

834 "be 1 or 2.") 

835 

836 deriv_value = np.asarray(deriv_value) 

837 if deriv_value.shape != expected_deriv_shape: 

838 raise ValueError( 

839 "`deriv_value` shape {} is not the expected one {}." 

840 .format(deriv_value.shape, expected_deriv_shape)) 

841 

842 if np.issubdtype(deriv_value.dtype, np.complexfloating): 

843 y = y.astype(complex, copy=False) 

844 

845 validated_bc.append((deriv_order, deriv_value)) 

846 

847 return validated_bc, y