Coverage for /usr/lib/python3/dist-packages/scipy/interpolate/_polyint.py: 26%

207 statements  

« prev     ^ index     » next       coverage.py v7.9.1, created at 2025-06-14 15:55 +0200

1import warnings 

2 

3import numpy as np 

4from scipy.special import factorial 

5from scipy._lib._util import _asarray_validated, float_factorial 

6 

7 

8__all__ = ["KroghInterpolator", "krogh_interpolate", "BarycentricInterpolator", 

9 "barycentric_interpolate", "approximate_taylor_polynomial"] 

10 

11 

12def _isscalar(x): 

13 """Check whether x is if a scalar type, or 0-dim""" 

14 return np.isscalar(x) or hasattr(x, 'shape') and x.shape == () 

15 

16 

17class _Interpolator1D: 

18 """ 

19 Common features in univariate interpolation 

20 

21 Deal with input data type and interpolation axis rolling. The 

22 actual interpolator can assume the y-data is of shape (n, r) where 

23 `n` is the number of x-points, and `r` the number of variables, 

24 and use self.dtype as the y-data type. 

25 

26 Attributes 

27 ---------- 

28 _y_axis 

29 Axis along which the interpolation goes in the original array 

30 _y_extra_shape 

31 Additional trailing shape of the input arrays, excluding 

32 the interpolation axis. 

33 dtype 

34 Dtype of the y-data arrays. Can be set via _set_dtype, which 

35 forces it to be float or complex. 

36 

37 Methods 

38 ------- 

39 __call__ 

40 _prepare_x 

41 _finish_y 

42 _reshape_yi 

43 _set_yi 

44 _set_dtype 

45 _evaluate 

46 

47 """ 

48 

49 __slots__ = ('_y_axis', '_y_extra_shape', 'dtype') 

50 

51 def __init__(self, xi=None, yi=None, axis=None): 

52 self._y_axis = axis 

53 self._y_extra_shape = None 

54 self.dtype = None 

55 if yi is not None: 

56 self._set_yi(yi, xi=xi, axis=axis) 

57 

58 def __call__(self, x): 

59 """ 

60 Evaluate the interpolant 

61 

62 Parameters 

63 ---------- 

64 x : array_like 

65 Points to evaluate the interpolant at. 

66 

67 Returns 

68 ------- 

69 y : array_like 

70 Interpolated values. Shape is determined by replacing 

71 the interpolation axis in the original array with the shape of x. 

72 

73 Notes 

74 ----- 

75 Input values `x` must be convertible to `float` values like `int` 

76 or `float`. 

77 

78 """ 

79 x, x_shape = self._prepare_x(x) 

80 y = self._evaluate(x) 

81 return self._finish_y(y, x_shape) 

82 

83 def _evaluate(self, x): 

84 """ 

85 Actually evaluate the value of the interpolator. 

86 """ 

87 raise NotImplementedError() 

88 

89 def _prepare_x(self, x): 

90 """Reshape input x array to 1-D""" 

91 x = _asarray_validated(x, check_finite=False, as_inexact=True) 

92 x_shape = x.shape 

93 return x.ravel(), x_shape 

94 

95 def _finish_y(self, y, x_shape): 

96 """Reshape interpolated y back to an N-D array similar to initial y""" 

97 y = y.reshape(x_shape + self._y_extra_shape) 

98 if self._y_axis != 0 and x_shape != (): 

99 nx = len(x_shape) 

100 ny = len(self._y_extra_shape) 

101 s = (list(range(nx, nx + self._y_axis)) 

102 + list(range(nx)) + list(range(nx+self._y_axis, nx+ny))) 

103 y = y.transpose(s) 

104 return y 

105 

106 def _reshape_yi(self, yi, check=False): 

107 yi = np.moveaxis(np.asarray(yi), self._y_axis, 0) 

108 if check and yi.shape[1:] != self._y_extra_shape: 

109 ok_shape = "{!r} + (N,) + {!r}".format(self._y_extra_shape[-self._y_axis:], 

110 self._y_extra_shape[:-self._y_axis]) 

111 raise ValueError("Data must be of shape %s" % ok_shape) 

112 return yi.reshape((yi.shape[0], -1)) 

113 

114 def _set_yi(self, yi, xi=None, axis=None): 

115 if axis is None: 

116 axis = self._y_axis 

117 if axis is None: 

118 raise ValueError("no interpolation axis specified") 

119 

120 yi = np.asarray(yi) 

121 

122 shape = yi.shape 

123 if shape == (): 

124 shape = (1,) 

125 if xi is not None and shape[axis] != len(xi): 

126 raise ValueError("x and y arrays must be equal in length along " 

127 "interpolation axis.") 

128 

129 self._y_axis = (axis % yi.ndim) 

130 self._y_extra_shape = yi.shape[:self._y_axis]+yi.shape[self._y_axis+1:] 

131 self.dtype = None 

132 self._set_dtype(yi.dtype) 

133 

134 def _set_dtype(self, dtype, union=False): 

135 if np.issubdtype(dtype, np.complexfloating) \ 

136 or np.issubdtype(self.dtype, np.complexfloating): 

137 self.dtype = np.complex_ 

138 else: 

139 if not union or self.dtype != np.complex_: 

140 self.dtype = np.float_ 

141 

142 

143class _Interpolator1DWithDerivatives(_Interpolator1D): 

144 def derivatives(self, x, der=None): 

145 """ 

146 Evaluate many derivatives of the polynomial at the point x 

147 

148 Produce an array of all derivative values at the point x. 

149 

150 Parameters 

151 ---------- 

152 x : array_like 

153 Point or points at which to evaluate the derivatives 

154 der : int or None, optional 

155 How many derivatives to extract; None for all potentially 

156 nonzero derivatives (that is a number equal to the number 

157 of points). This number includes the function value as 0th 

158 derivative. 

159 

160 Returns 

161 ------- 

162 d : ndarray 

163 Array with derivatives; d[j] contains the jth derivative. 

164 Shape of d[j] is determined by replacing the interpolation 

165 axis in the original array with the shape of x. 

166 

167 Examples 

168 -------- 

169 >>> from scipy.interpolate import KroghInterpolator 

170 >>> KroghInterpolator([0,0,0],[1,2,3]).derivatives(0) 

171 array([1.0,2.0,3.0]) 

172 >>> KroghInterpolator([0,0,0],[1,2,3]).derivatives([0,0]) 

173 array([[1.0,1.0], 

174 [2.0,2.0], 

175 [3.0,3.0]]) 

176 

177 """ 

178 x, x_shape = self._prepare_x(x) 

179 y = self._evaluate_derivatives(x, der) 

180 

181 y = y.reshape((y.shape[0],) + x_shape + self._y_extra_shape) 

182 if self._y_axis != 0 and x_shape != (): 

183 nx = len(x_shape) 

184 ny = len(self._y_extra_shape) 

185 s = ([0] + list(range(nx+1, nx + self._y_axis+1)) 

186 + list(range(1, nx+1)) + 

187 list(range(nx+1+self._y_axis, nx+ny+1))) 

188 y = y.transpose(s) 

189 return y 

190 

191 def derivative(self, x, der=1): 

192 """ 

193 Evaluate one derivative of the polynomial at the point x 

194 

195 Parameters 

196 ---------- 

197 x : array_like 

198 Point or points at which to evaluate the derivatives 

199 

200 der : integer, optional 

201 Which derivative to extract. This number includes the 

202 function value as 0th derivative. 

203 

204 Returns 

205 ------- 

206 d : ndarray 

207 Derivative interpolated at the x-points. Shape of d is 

208 determined by replacing the interpolation axis in the 

209 original array with the shape of x. 

210 

211 Notes 

212 ----- 

213 This is computed by evaluating all derivatives up to the desired 

214 one (using self.derivatives()) and then discarding the rest. 

215 

216 """ 

217 x, x_shape = self._prepare_x(x) 

218 y = self._evaluate_derivatives(x, der+1) 

219 return self._finish_y(y[der], x_shape) 

220 

221 

222class KroghInterpolator(_Interpolator1DWithDerivatives): 

223 """ 

224 Interpolating polynomial for a set of points. 

225 

226 The polynomial passes through all the pairs (xi,yi). One may 

227 additionally specify a number of derivatives at each point xi; 

228 this is done by repeating the value xi and specifying the 

229 derivatives as successive yi values. 

230 

231 Allows evaluation of the polynomial and all its derivatives. 

232 For reasons of numerical stability, this function does not compute 

233 the coefficients of the polynomial, although they can be obtained 

234 by evaluating all the derivatives. 

235 

236 Parameters 

237 ---------- 

238 xi : array_like, shape (npoints, ) 

239 Known x-coordinates. Must be sorted in increasing order. 

240 yi : array_like, shape (..., npoints, ...) 

241 Known y-coordinates. When an xi occurs two or more times in 

242 a row, the corresponding yi's represent derivative values. The length of `yi` 

243 along the interpolation axis must be equal to the length of `xi`. Use the 

244 `axis` parameter to select the correct axis. 

245 axis : int, optional 

246 Axis in the `yi` array corresponding to the x-coordinate values. Defaults to 

247 ``axis=0``. 

248 

249 Notes 

250 ----- 

251 Be aware that the algorithms implemented here are not necessarily 

252 the most numerically stable known. Moreover, even in a world of 

253 exact computation, unless the x coordinates are chosen very 

254 carefully - Chebyshev zeros (e.g., cos(i*pi/n)) are a good choice - 

255 polynomial interpolation itself is a very ill-conditioned process 

256 due to the Runge phenomenon. In general, even with well-chosen 

257 x values, degrees higher than about thirty cause problems with 

258 numerical instability in this code. 

259 

260 Based on [1]_. 

261 

262 References 

263 ---------- 

264 .. [1] Krogh, "Efficient Algorithms for Polynomial Interpolation 

265 and Numerical Differentiation", 1970. 

266 

267 Examples 

268 -------- 

269 To produce a polynomial that is zero at 0 and 1 and has 

270 derivative 2 at 0, call 

271 

272 >>> from scipy.interpolate import KroghInterpolator 

273 >>> KroghInterpolator([0,0,1],[0,2,0]) 

274 

275 This constructs the quadratic 2*X**2-2*X. The derivative condition 

276 is indicated by the repeated zero in the xi array; the corresponding 

277 yi values are 0, the function value, and 2, the derivative value. 

278 

279 For another example, given xi, yi, and a derivative ypi for each 

280 point, appropriate arrays can be constructed as: 

281 

282 >>> import numpy as np 

283 >>> rng = np.random.default_rng() 

284 >>> xi = np.linspace(0, 1, 5) 

285 >>> yi, ypi = rng.random((2, 5)) 

286 >>> xi_k, yi_k = np.repeat(xi, 2), np.ravel(np.dstack((yi,ypi))) 

287 >>> KroghInterpolator(xi_k, yi_k) 

288 

289 To produce a vector-valued polynomial, supply a higher-dimensional 

290 array for yi: 

291 

292 >>> KroghInterpolator([0,1],[[2,3],[4,5]]) 

293 

294 This constructs a linear polynomial giving (2,3) at 0 and (4,5) at 1. 

295 

296 """ 

297 

298 def __init__(self, xi, yi, axis=0): 

299 _Interpolator1DWithDerivatives.__init__(self, xi, yi, axis) 

300 

301 self.xi = np.asarray(xi) 

302 self.yi = self._reshape_yi(yi) 

303 self.n, self.r = self.yi.shape 

304 

305 if (deg := self.xi.size) > 30: 

306 warnings.warn(f"{deg} degrees provided, degrees higher than about" 

307 " thirty cause problems with numerical instability " 

308 "with 'KroghInterpolator'", stacklevel=2) 

309 

310 c = np.zeros((self.n+1, self.r), dtype=self.dtype) 

311 c[0] = self.yi[0] 

312 Vk = np.zeros((self.n, self.r), dtype=self.dtype) 

313 for k in range(1, self.n): 

314 s = 0 

315 while s <= k and xi[k-s] == xi[k]: 

316 s += 1 

317 s -= 1 

318 Vk[0] = self.yi[k]/float_factorial(s) 

319 for i in range(k-s): 

320 if xi[i] == xi[k]: 

321 raise ValueError("Elements if `xi` can't be equal.") 

322 if s == 0: 

323 Vk[i+1] = (c[i]-Vk[i])/(xi[i]-xi[k]) 

324 else: 

325 Vk[i+1] = (Vk[i+1]-Vk[i])/(xi[i]-xi[k]) 

326 c[k] = Vk[k-s] 

327 self.c = c 

328 

329 def _evaluate(self, x): 

330 pi = 1 

331 p = np.zeros((len(x), self.r), dtype=self.dtype) 

332 p += self.c[0,np.newaxis,:] 

333 for k in range(1, self.n): 

334 w = x - self.xi[k-1] 

335 pi = w*pi 

336 p += pi[:,np.newaxis] * self.c[k] 

337 return p 

338 

339 def _evaluate_derivatives(self, x, der=None): 

340 n = self.n 

341 r = self.r 

342 

343 if der is None: 

344 der = self.n 

345 pi = np.zeros((n, len(x))) 

346 w = np.zeros((n, len(x))) 

347 pi[0] = 1 

348 p = np.zeros((len(x), self.r), dtype=self.dtype) 

349 p += self.c[0, np.newaxis, :] 

350 

351 for k in range(1, n): 

352 w[k-1] = x - self.xi[k-1] 

353 pi[k] = w[k-1] * pi[k-1] 

354 p += pi[k, :, np.newaxis] * self.c[k] 

355 

356 cn = np.zeros((max(der, n+1), len(x), r), dtype=self.dtype) 

357 cn[:n+1, :, :] += self.c[:n+1, np.newaxis, :] 

358 cn[0] = p 

359 for k in range(1, n): 

360 for i in range(1, n-k+1): 

361 pi[i] = w[k+i-1]*pi[i-1] + pi[i] 

362 cn[k] = cn[k] + pi[i, :, np.newaxis]*cn[k+i] 

363 cn[k] *= float_factorial(k) 

364 

365 cn[n, :, :] = 0 

366 return cn[:der] 

367 

368 

369def krogh_interpolate(xi, yi, x, der=0, axis=0): 

370 """ 

371 Convenience function for polynomial interpolation. 

372 

373 See `KroghInterpolator` for more details. 

374 

375 Parameters 

376 ---------- 

377 xi : array_like 

378 Known x-coordinates. 

379 yi : array_like 

380 Known y-coordinates, of shape ``(xi.size, R)``. Interpreted as 

381 vectors of length R, or scalars if R=1. 

382 x : array_like 

383 Point or points at which to evaluate the derivatives. 

384 der : int or list, optional 

385 How many derivatives to extract; None for all potentially 

386 nonzero derivatives (that is a number equal to the number 

387 of points), or a list of derivatives to extract. This number 

388 includes the function value as 0th derivative. 

389 axis : int, optional 

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

391 

392 Returns 

393 ------- 

394 d : ndarray 

395 If the interpolator's values are R-D then the 

396 returned array will be the number of derivatives by N by R. 

397 If `x` is a scalar, the middle dimension will be dropped; if 

398 the `yi` are scalars then the last dimension will be dropped. 

399 

400 See Also 

401 -------- 

402 KroghInterpolator : Krogh interpolator 

403 

404 Notes 

405 ----- 

406 Construction of the interpolating polynomial is a relatively expensive 

407 process. If you want to evaluate it repeatedly consider using the class 

408 KroghInterpolator (which is what this function uses). 

409 

410 Examples 

411 -------- 

412 We can interpolate 2D observed data using krogh interpolation: 

413 

414 >>> import numpy as np 

415 >>> import matplotlib.pyplot as plt 

416 >>> from scipy.interpolate import krogh_interpolate 

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

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

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

420 >>> y = krogh_interpolate(x_observed, y_observed, x) 

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

422 >>> plt.plot(x, y, label="krogh interpolation") 

423 >>> plt.legend() 

424 >>> plt.show() 

425 

426 """ 

427 P = KroghInterpolator(xi, yi, axis=axis) 

428 if der == 0: 

429 return P(x) 

430 elif _isscalar(der): 

431 return P.derivative(x,der=der) 

432 else: 

433 return P.derivatives(x,der=np.amax(der)+1)[der] 

434 

435 

436def approximate_taylor_polynomial(f,x,degree,scale,order=None): 

437 """ 

438 Estimate the Taylor polynomial of f at x by polynomial fitting. 

439 

440 Parameters 

441 ---------- 

442 f : callable 

443 The function whose Taylor polynomial is sought. Should accept 

444 a vector of `x` values. 

445 x : scalar 

446 The point at which the polynomial is to be evaluated. 

447 degree : int 

448 The degree of the Taylor polynomial 

449 scale : scalar 

450 The width of the interval to use to evaluate the Taylor polynomial. 

451 Function values spread over a range this wide are used to fit the 

452 polynomial. Must be chosen carefully. 

453 order : int or None, optional 

454 The order of the polynomial to be used in the fitting; `f` will be 

455 evaluated ``order+1`` times. If None, use `degree`. 

456 

457 Returns 

458 ------- 

459 p : poly1d instance 

460 The Taylor polynomial (translated to the origin, so that 

461 for example p(0)=f(x)). 

462 

463 Notes 

464 ----- 

465 The appropriate choice of "scale" is a trade-off; too large and the 

466 function differs from its Taylor polynomial too much to get a good 

467 answer, too small and round-off errors overwhelm the higher-order terms. 

468 The algorithm used becomes numerically unstable around order 30 even 

469 under ideal circumstances. 

470 

471 Choosing order somewhat larger than degree may improve the higher-order 

472 terms. 

473 

474 Examples 

475 -------- 

476 We can calculate Taylor approximation polynomials of sin function with 

477 various degrees: 

478 

479 >>> import numpy as np 

480 >>> import matplotlib.pyplot as plt 

481 >>> from scipy.interpolate import approximate_taylor_polynomial 

482 >>> x = np.linspace(-10.0, 10.0, num=100) 

483 >>> plt.plot(x, np.sin(x), label="sin curve") 

484 >>> for degree in np.arange(1, 15, step=2): 

485 ... sin_taylor = approximate_taylor_polynomial(np.sin, 0, degree, 1, 

486 ... order=degree + 2) 

487 ... plt.plot(x, sin_taylor(x), label=f"degree={degree}") 

488 >>> plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', 

489 ... borderaxespad=0.0, shadow=True) 

490 >>> plt.tight_layout() 

491 >>> plt.axis([-10, 10, -10, 10]) 

492 >>> plt.show() 

493 

494 """ 

495 if order is None: 

496 order = degree 

497 

498 n = order+1 

499 # Choose n points that cluster near the endpoints of the interval in 

500 # a way that avoids the Runge phenomenon. Ensure, by including the 

501 # endpoint or not as appropriate, that one point always falls at x 

502 # exactly. 

503 xs = scale*np.cos(np.linspace(0,np.pi,n,endpoint=n % 1)) + x 

504 

505 P = KroghInterpolator(xs, f(xs)) 

506 d = P.derivatives(x,der=degree+1) 

507 

508 return np.poly1d((d/factorial(np.arange(degree+1)))[::-1]) 

509 

510 

511class BarycentricInterpolator(_Interpolator1D): 

512 """The interpolating polynomial for a set of points 

513 

514 Constructs a polynomial that passes through a given set of points. 

515 Allows evaluation of the polynomial, efficient changing of the y 

516 values to be interpolated, and updating by adding more x values. 

517 For reasons of numerical stability, this function does not compute 

518 the coefficients of the polynomial. 

519 

520 The values yi need to be provided before the function is 

521 evaluated, but none of the preprocessing depends on them, so rapid 

522 updates are possible. 

523 

524 Parameters 

525 ---------- 

526 xi : array_like, shape (npoints, ) 

527 1-D array of x coordinates of the points the polynomial 

528 should pass through 

529 yi : array_like, shape (..., npoints, ...), optional 

530 N-D array of y coordinates of the points the polynomial should pass through. 

531 If None, the y values will be supplied later via the `set_y` method. 

532 The length of `yi` along the interpolation axis must be equal to the length 

533 of `xi`. Use the ``axis`` parameter to select correct axis. 

534 axis : int, optional 

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

536 to ``axis=0``. 

537 

538 Notes 

539 ----- 

540 This class uses a "barycentric interpolation" method that treats 

541 the problem as a special case of rational function interpolation. 

542 This algorithm is quite stable, numerically, but even in a world of 

543 exact computation, unless the x coordinates are chosen very 

544 carefully - Chebyshev zeros (e.g., cos(i*pi/n)) are a good choice - 

545 polynomial interpolation itself is a very ill-conditioned process 

546 due to the Runge phenomenon. 

547 

548 Based on Berrut and Trefethen 2004, "Barycentric Lagrange Interpolation". 

549 

550 """ 

551 

552 def __init__(self, xi, yi=None, axis=0): 

553 _Interpolator1D.__init__(self, xi, yi, axis) 

554 

555 self.xi = np.asfarray(xi) 

556 self.set_yi(yi) 

557 self.n = len(self.xi) 

558 

559 # See page 510 of Berrut and Trefethen 2004 for an explanation of the 

560 # capacity scaling and the suggestion of using a random permutation of 

561 # the input factors. 

562 # At the moment, the permutation is not performed for xi that are 

563 # appended later through the add_xi interface. It's not clear to me how 

564 # to implement that and it seems that most situations that require 

565 # these numerical stability improvements will be able to provide all 

566 # the points to the constructor. 

567 self._inv_capacity = 4.0 / (np.max(self.xi) - np.min(self.xi)) 

568 permute = np.random.permutation(self.n) 

569 inv_permute = np.zeros(self.n, dtype=np.int32) 

570 inv_permute[permute] = np.arange(self.n) 

571 

572 self.wi = np.zeros(self.n) 

573 for i in range(self.n): 

574 dist = self._inv_capacity * (self.xi[i] - self.xi[permute]) 

575 dist[inv_permute[i]] = 1.0 

576 self.wi[i] = 1.0 / np.prod(dist) 

577 

578 def set_yi(self, yi, axis=None): 

579 """ 

580 Update the y values to be interpolated 

581 

582 The barycentric interpolation algorithm requires the calculation 

583 of weights, but these depend only on the xi. The yi can be changed 

584 at any time. 

585 

586 Parameters 

587 ---------- 

588 yi : array_like 

589 The y coordinates of the points the polynomial should pass through. 

590 If None, the y values will be supplied later. 

591 axis : int, optional 

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

593 

594 """ 

595 if yi is None: 

596 self.yi = None 

597 return 

598 self._set_yi(yi, xi=self.xi, axis=axis) 

599 self.yi = self._reshape_yi(yi) 

600 self.n, self.r = self.yi.shape 

601 

602 def add_xi(self, xi, yi=None): 

603 """ 

604 Add more x values to the set to be interpolated 

605 

606 The barycentric interpolation algorithm allows easy updating by 

607 adding more points for the polynomial to pass through. 

608 

609 Parameters 

610 ---------- 

611 xi : array_like 

612 The x coordinates of the points that the polynomial should pass 

613 through. 

614 yi : array_like, optional 

615 The y coordinates of the points the polynomial should pass through. 

616 Should have shape ``(xi.size, R)``; if R > 1 then the polynomial is 

617 vector-valued. 

618 If `yi` is not given, the y values will be supplied later. `yi` 

619 should be given if and only if the interpolator has y values 

620 specified. 

621 

622 """ 

623 if yi is not None: 

624 if self.yi is None: 

625 raise ValueError("No previous yi value to update!") 

626 yi = self._reshape_yi(yi, check=True) 

627 self.yi = np.vstack((self.yi,yi)) 

628 else: 

629 if self.yi is not None: 

630 raise ValueError("No update to yi provided!") 

631 old_n = self.n 

632 self.xi = np.concatenate((self.xi,xi)) 

633 self.n = len(self.xi) 

634 self.wi **= -1 

635 old_wi = self.wi 

636 self.wi = np.zeros(self.n) 

637 self.wi[:old_n] = old_wi 

638 for j in range(old_n, self.n): 

639 self.wi[:j] *= self._inv_capacity * (self.xi[j]-self.xi[:j]) 

640 self.wi[j] = np.multiply.reduce( 

641 self._inv_capacity * (self.xi[:j]-self.xi[j]) 

642 ) 

643 self.wi **= -1 

644 

645 def __call__(self, x): 

646 """Evaluate the interpolating polynomial at the points x 

647 

648 Parameters 

649 ---------- 

650 x : array_like 

651 Points to evaluate the interpolant at. 

652 

653 Returns 

654 ------- 

655 y : array_like 

656 Interpolated values. Shape is determined by replacing 

657 the interpolation axis in the original array with the shape of x. 

658 

659 Notes 

660 ----- 

661 Currently the code computes an outer product between x and the 

662 weights, that is, it constructs an intermediate array of size 

663 N by len(x), where N is the degree of the polynomial. 

664 """ 

665 return _Interpolator1D.__call__(self, x) 

666 

667 def _evaluate(self, x): 

668 if x.size == 0: 

669 p = np.zeros((0, self.r), dtype=self.dtype) 

670 else: 

671 c = x[..., np.newaxis] - self.xi 

672 z = c == 0 

673 c[z] = 1 

674 c = self.wi/c 

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

676 p = np.dot(c, self.yi) / np.sum(c, axis=-1)[..., np.newaxis] 

677 # Now fix where x==some xi 

678 r = np.nonzero(z) 

679 if len(r) == 1: # evaluation at a scalar 

680 if len(r[0]) > 0: # equals one of the points 

681 p = self.yi[r[0][0]] 

682 else: 

683 p[r[:-1]] = self.yi[r[-1]] 

684 return p 

685 

686 

687def barycentric_interpolate(xi, yi, x, axis=0): 

688 """ 

689 Convenience function for polynomial interpolation. 

690 

691 Constructs a polynomial that passes through a given set of points, 

692 then evaluates the polynomial. For reasons of numerical stability, 

693 this function does not compute the coefficients of the polynomial. 

694 

695 This function uses a "barycentric interpolation" method that treats 

696 the problem as a special case of rational function interpolation. 

697 This algorithm is quite stable, numerically, but even in a world of 

698 exact computation, unless the `x` coordinates are chosen very 

699 carefully - Chebyshev zeros (e.g., cos(i*pi/n)) are a good choice - 

700 polynomial interpolation itself is a very ill-conditioned process 

701 due to the Runge phenomenon. 

702 

703 Parameters 

704 ---------- 

705 xi : array_like 

706 1-D array of x coordinates of the points the polynomial should 

707 pass through 

708 yi : array_like 

709 The y coordinates of the points the polynomial should pass through. 

710 x : scalar or array_like 

711 Points to evaluate the interpolator at. 

712 axis : int, optional 

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

714 

715 Returns 

716 ------- 

717 y : scalar or array_like 

718 Interpolated values. Shape is determined by replacing 

719 the interpolation axis in the original array with the shape of x. 

720 

721 See Also 

722 -------- 

723 BarycentricInterpolator : Bary centric interpolator 

724 

725 Notes 

726 ----- 

727 Construction of the interpolation weights is a relatively slow process. 

728 If you want to call this many times with the same xi (but possibly 

729 varying yi or x) you should use the class `BarycentricInterpolator`. 

730 This is what this function uses internally. 

731 

732 Examples 

733 -------- 

734 We can interpolate 2D observed data using barycentric interpolation: 

735 

736 >>> import numpy as np 

737 >>> import matplotlib.pyplot as plt 

738 >>> from scipy.interpolate import barycentric_interpolate 

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

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

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

742 >>> y = barycentric_interpolate(x_observed, y_observed, x) 

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

744 >>> plt.plot(x, y, label="barycentric interpolation") 

745 >>> plt.legend() 

746 >>> plt.show() 

747 

748 """ 

749 return BarycentricInterpolator(xi, yi, axis=axis)(x)