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""" 

2fitpack --- curve and surface fitting with splines 

3 

4fitpack is based on a collection of Fortran routines DIERCKX 

5by P. Dierckx (see http://www.netlib.org/dierckx/) transformed 

6to double routines by Pearu Peterson. 

7""" 

8# Created by Pearu Peterson, June,August 2003 

9__all__ = [ 

10 'UnivariateSpline', 

11 'InterpolatedUnivariateSpline', 

12 'LSQUnivariateSpline', 

13 'BivariateSpline', 

14 'LSQBivariateSpline', 

15 'SmoothBivariateSpline', 

16 'LSQSphereBivariateSpline', 

17 'SmoothSphereBivariateSpline', 

18 'RectBivariateSpline', 

19 'RectSphereBivariateSpline'] 

20 

21 

22import warnings 

23 

24from numpy import zeros, concatenate, ravel, diff, array, ones 

25import numpy as np 

26 

27from . import fitpack 

28from . import dfitpack 

29 

30 

31dfitpack_int = dfitpack.types.intvar.dtype 

32 

33 

34# ############### Univariate spline #################### 

35 

36_curfit_messages = {1: """ 

37The required storage space exceeds the available storage space, as 

38specified by the parameter nest: nest too small. If nest is already 

39large (say nest > m/2), it may also indicate that s is too small. 

40The approximation returned is the weighted least-squares spline 

41according to the knots t[0],t[1],...,t[n-1]. (n=nest) the parameter fp 

42gives the corresponding weighted sum of squared residuals (fp>s). 

43""", 

44 2: """ 

45A theoretically impossible result was found during the iteration 

46process for finding a smoothing spline with fp = s: s too small. 

47There is an approximation returned but the corresponding weighted sum 

48of squared residuals does not satisfy the condition abs(fp-s)/s < tol.""", 

49 3: """ 

50The maximal number of iterations maxit (set to 20 by the program) 

51allowed for finding a smoothing spline with fp=s has been reached: s 

52too small. 

53There is an approximation returned but the corresponding weighted sum 

54of squared residuals does not satisfy the condition abs(fp-s)/s < tol.""", 

55 10: """ 

56Error on entry, no approximation returned. The following conditions 

57must hold: 

58xb<=x[0]<x[1]<...<x[m-1]<=xe, w[i]>0, i=0..m-1 

59if iopt=-1: 

60 xb<t[k+1]<t[k+2]<...<t[n-k-2]<xe""" 

61 } 

62 

63 

64# UnivariateSpline, ext parameter can be an int or a string 

65_extrap_modes = {0: 0, 'extrapolate': 0, 

66 1: 1, 'zeros': 1, 

67 2: 2, 'raise': 2, 

68 3: 3, 'const': 3} 

69 

70 

71class UnivariateSpline(object): 

72 """ 

73 1-D smoothing spline fit to a given set of data points. 

74 

75 Fits a spline y = spl(x) of degree `k` to the provided `x`, `y` data. `s` 

76 specifies the number of knots by specifying a smoothing condition. 

77 

78 Parameters 

79 ---------- 

80 x : (N,) array_like 

81 1-D array of independent input data. Must be increasing; 

82 must be strictly increasing if `s` is 0. 

83 y : (N,) array_like 

84 1-D array of dependent input data, of the same length as `x`. 

85 w : (N,) array_like, optional 

86 Weights for spline fitting. Must be positive. If None (default), 

87 weights are all equal. 

88 bbox : (2,) array_like, optional 

89 2-sequence specifying the boundary of the approximation interval. If 

90 None (default), ``bbox=[x[0], x[-1]]``. 

91 k : int, optional 

92 Degree of the smoothing spline. Must be 1 <= `k` <= 5. 

93 Default is `k` = 3, a cubic spline. 

94 s : float or None, optional 

95 Positive smoothing factor used to choose the number of knots. Number 

96 of knots will be increased until the smoothing condition is satisfied:: 

97 

98 sum((w[i] * (y[i]-spl(x[i])))**2, axis=0) <= s 

99 

100 If None (default), ``s = len(w)`` which should be a good value if 

101 ``1/w[i]`` is an estimate of the standard deviation of ``y[i]``. 

102 If 0, spline will interpolate through all data points. 

103 ext : int or str, optional 

104 Controls the extrapolation mode for elements 

105 not in the interval defined by the knot sequence. 

106 

107 * if ext=0 or 'extrapolate', return the extrapolated value. 

108 * if ext=1 or 'zeros', return 0 

109 * if ext=2 or 'raise', raise a ValueError 

110 * if ext=3 of 'const', return the boundary value. 

111 

112 The default value is 0. 

113 

114 check_finite : bool, optional 

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

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

117 (crashes, non-termination or non-sensical results) if the inputs 

118 do contain infinities or NaNs. 

119 Default is False. 

120 

121 See Also 

122 -------- 

123 InterpolatedUnivariateSpline : Subclass with smoothing forced to 0 

124 LSQUnivariateSpline : Subclass in which knots are user-selected instead of 

125 being set by smoothing condition 

126 splrep : An older, non object-oriented wrapping of FITPACK 

127 splev, sproot, splint, spalde 

128 BivariateSpline : A similar class for two-dimensional spline interpolation 

129 

130 Notes 

131 ----- 

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

133 

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

135 is not useful, since the underlying spline fitting routines cannot deal 

136 with ``nan``. A workaround is to use zero weights for not-a-number 

137 data points: 

138 

139 >>> from scipy.interpolate import UnivariateSpline 

140 >>> x, y = np.array([1, 2, 3, 4]), np.array([1, np.nan, 3, 4]) 

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

142 >>> y[w] = 0. 

143 >>> spl = UnivariateSpline(x, y, w=~w) 

144 

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

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

147 

148 Examples 

149 -------- 

150 >>> import matplotlib.pyplot as plt 

151 >>> from scipy.interpolate import UnivariateSpline 

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

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

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

155 

156 Use the default value for the smoothing parameter: 

157 

158 >>> spl = UnivariateSpline(x, y) 

159 >>> xs = np.linspace(-3, 3, 1000) 

160 >>> plt.plot(xs, spl(xs), 'g', lw=3) 

161 

162 Manually change the amount of smoothing: 

163 

164 >>> spl.set_smoothing_factor(0.5) 

165 >>> plt.plot(xs, spl(xs), 'b', lw=3) 

166 >>> plt.show() 

167 

168 """ 

169 def __init__(self, x, y, w=None, bbox=[None]*2, k=3, s=None, 

170 ext=0, check_finite=False): 

171 

172 x, y, w, bbox, self.ext = self.validate_input(x, y, w, bbox, k, s, ext, 

173 check_finite) 

174 

175 # _data == x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier 

176 data = dfitpack.fpcurf0(x, y, k, w=w, xb=bbox[0], 

177 xe=bbox[1], s=s) 

178 if data[-1] == 1: 

179 # nest too small, setting to maximum bound 

180 data = self._reset_nest(data) 

181 self._data = data 

182 self._reset_class() 

183 

184 @staticmethod 

185 def validate_input(x, y, w, bbox, k, s, ext, check_finite): 

186 x, y, bbox = np.asarray(x), np.asarray(y), np.asarray(bbox) 

187 if w is not None: 

188 w = np.asarray(w) 

189 if check_finite: 

190 w_finite = np.isfinite(w).all() if w is not None else True 

191 if (not np.isfinite(x).all() or not np.isfinite(y).all() or 

192 not w_finite): 

193 raise ValueError("x and y array must not contain " 

194 "NaNs or infs.") 

195 if s is None or s > 0: 

196 if not np.all(diff(x) >= 0.0): 

197 raise ValueError("x must be increasing if s > 0") 

198 else: 

199 if not np.all(diff(x) > 0.0): 

200 raise ValueError("x must be strictly increasing if s = 0") 

201 if x.size != y.size: 

202 raise ValueError("x and y should have a same length") 

203 elif w is not None and not x.size == y.size == w.size: 

204 raise ValueError("x, y, and w should have a same length") 

205 elif bbox.shape != (2,): 

206 raise ValueError("bbox shape should be (2,)") 

207 elif not (1 <= k <= 5): 

208 raise ValueError("k should be 1 <= k <= 5") 

209 elif s is not None and not s >= 0.0: 

210 raise ValueError("s should be s >= 0.0") 

211 

212 try: 

213 ext = _extrap_modes[ext] 

214 except KeyError: 

215 raise ValueError("Unknown extrapolation mode %s." % ext) 

216 

217 return x, y, w, bbox, ext 

218 

219 @classmethod 

220 def _from_tck(cls, tck, ext=0): 

221 """Construct a spline object from given tck""" 

222 self = cls.__new__(cls) 

223 t, c, k = tck 

224 self._eval_args = tck 

225 # _data == x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier 

226 self._data = (None, None, None, None, None, k, None, len(t), t, 

227 c, None, None, None, None) 

228 self.ext = ext 

229 return self 

230 

231 def _reset_class(self): 

232 data = self._data 

233 n, t, c, k, ier = data[7], data[8], data[9], data[5], data[-1] 

234 self._eval_args = t[:n], c[:n], k 

235 if ier == 0: 

236 # the spline returned has a residual sum of squares fp 

237 # such that abs(fp-s)/s <= tol with tol a relative 

238 # tolerance set to 0.001 by the program 

239 pass 

240 elif ier == -1: 

241 # the spline returned is an interpolating spline 

242 self._set_class(InterpolatedUnivariateSpline) 

243 elif ier == -2: 

244 # the spline returned is the weighted least-squares 

245 # polynomial of degree k. In this extreme case fp gives 

246 # the upper bound fp0 for the smoothing factor s. 

247 self._set_class(LSQUnivariateSpline) 

248 else: 

249 # error 

250 if ier == 1: 

251 self._set_class(LSQUnivariateSpline) 

252 message = _curfit_messages.get(ier, 'ier=%s' % (ier)) 

253 warnings.warn(message) 

254 

255 def _set_class(self, cls): 

256 self._spline_class = cls 

257 if self.__class__ in (UnivariateSpline, InterpolatedUnivariateSpline, 

258 LSQUnivariateSpline): 

259 self.__class__ = cls 

260 else: 

261 # It's an unknown subclass -- don't change class. cf. #731 

262 pass 

263 

264 def _reset_nest(self, data, nest=None): 

265 n = data[10] 

266 if nest is None: 

267 k, m = data[5], len(data[0]) 

268 nest = m+k+1 # this is the maximum bound for nest 

269 else: 

270 if not n <= nest: 

271 raise ValueError("`nest` can only be increased") 

272 t, c, fpint, nrdata = [np.resize(data[j], nest) for j in 

273 [8, 9, 11, 12]] 

274 

275 args = data[:8] + (t, c, n, fpint, nrdata, data[13]) 

276 data = dfitpack.fpcurf1(*args) 

277 return data 

278 

279 def set_smoothing_factor(self, s): 

280 """ Continue spline computation with the given smoothing 

281 factor s and with the knots found at the last call. 

282 

283 This routine modifies the spline in place. 

284 

285 """ 

286 data = self._data 

287 if data[6] == -1: 

288 warnings.warn('smoothing factor unchanged for' 

289 'LSQ spline with fixed knots') 

290 return 

291 args = data[:6] + (s,) + data[7:] 

292 data = dfitpack.fpcurf1(*args) 

293 if data[-1] == 1: 

294 # nest too small, setting to maximum bound 

295 data = self._reset_nest(data) 

296 self._data = data 

297 self._reset_class() 

298 

299 def __call__(self, x, nu=0, ext=None): 

300 """ 

301 Evaluate spline (or its nu-th derivative) at positions x. 

302 

303 Parameters 

304 ---------- 

305 x : array_like 

306 A 1-D array of points at which to return the value of the smoothed 

307 spline or its derivatives. Note: `x` can be unordered but the 

308 evaluation is more efficient if `x` is (partially) ordered. 

309 nu : int 

310 The order of derivative of the spline to compute. 

311 ext : int 

312 Controls the value returned for elements of `x` not in the 

313 interval defined by the knot sequence. 

314 

315 * if ext=0 or 'extrapolate', return the extrapolated value. 

316 * if ext=1 or 'zeros', return 0 

317 * if ext=2 or 'raise', raise a ValueError 

318 * if ext=3 or 'const', return the boundary value. 

319 

320 The default value is 0, passed from the initialization of 

321 UnivariateSpline. 

322 

323 """ 

324 x = np.asarray(x) 

325 # empty input yields empty output 

326 if x.size == 0: 

327 return array([]) 

328 if ext is None: 

329 ext = self.ext 

330 else: 

331 try: 

332 ext = _extrap_modes[ext] 

333 except KeyError: 

334 raise ValueError("Unknown extrapolation mode %s." % ext) 

335 return fitpack.splev(x, self._eval_args, der=nu, ext=ext) 

336 

337 def get_knots(self): 

338 """ Return positions of interior knots of the spline. 

339 

340 Internally, the knot vector contains ``2*k`` additional boundary knots. 

341 """ 

342 data = self._data 

343 k, n = data[5], data[7] 

344 return data[8][k:n-k] 

345 

346 def get_coeffs(self): 

347 """Return spline coefficients.""" 

348 data = self._data 

349 k, n = data[5], data[7] 

350 return data[9][:n-k-1] 

351 

352 def get_residual(self): 

353 """Return weighted sum of squared residuals of the spline approximation. 

354 

355 This is equivalent to:: 

356 

357 sum((w[i] * (y[i]-spl(x[i])))**2, axis=0) 

358 

359 """ 

360 return self._data[10] 

361 

362 def integral(self, a, b): 

363 """ Return definite integral of the spline between two given points. 

364 

365 Parameters 

366 ---------- 

367 a : float 

368 Lower limit of integration. 

369 b : float 

370 Upper limit of integration. 

371 

372 Returns 

373 ------- 

374 integral : float 

375 The value of the definite integral of the spline between limits. 

376 

377 Examples 

378 -------- 

379 >>> from scipy.interpolate import UnivariateSpline 

380 >>> x = np.linspace(0, 3, 11) 

381 >>> y = x**2 

382 >>> spl = UnivariateSpline(x, y) 

383 >>> spl.integral(0, 3) 

384 9.0 

385 

386 which agrees with :math:`\\int x^2 dx = x^3 / 3` between the limits 

387 of 0 and 3. 

388 

389 A caveat is that this routine assumes the spline to be zero outside of 

390 the data limits: 

391 

392 >>> spl.integral(-1, 4) 

393 9.0 

394 >>> spl.integral(-1, 0) 

395 0.0 

396 

397 """ 

398 return dfitpack.splint(*(self._eval_args+(a, b))) 

399 

400 def derivatives(self, x): 

401 """ Return all derivatives of the spline at the point x. 

402 

403 Parameters 

404 ---------- 

405 x : float 

406 The point to evaluate the derivatives at. 

407 

408 Returns 

409 ------- 

410 der : ndarray, shape(k+1,) 

411 Derivatives of the orders 0 to k. 

412 

413 Examples 

414 -------- 

415 >>> from scipy.interpolate import UnivariateSpline 

416 >>> x = np.linspace(0, 3, 11) 

417 >>> y = x**2 

418 >>> spl = UnivariateSpline(x, y) 

419 >>> spl.derivatives(1.5) 

420 array([2.25, 3.0, 2.0, 0]) 

421 

422 """ 

423 d, ier = dfitpack.spalde(*(self._eval_args+(x,))) 

424 if not ier == 0: 

425 raise ValueError("Error code returned by spalde: %s" % ier) 

426 return d 

427 

428 def roots(self): 

429 """ Return the zeros of the spline. 

430 

431 Restriction: only cubic splines are supported by fitpack. 

432 """ 

433 k = self._data[5] 

434 if k == 3: 

435 z, m, ier = dfitpack.sproot(*self._eval_args[:2]) 

436 if not ier == 0: 

437 raise ValueError("Error code returned by spalde: %s" % ier) 

438 return z[:m] 

439 raise NotImplementedError('finding roots unsupported for ' 

440 'non-cubic splines') 

441 

442 def derivative(self, n=1): 

443 """ 

444 Construct a new spline representing the derivative of this spline. 

445 

446 Parameters 

447 ---------- 

448 n : int, optional 

449 Order of derivative to evaluate. Default: 1 

450 

451 Returns 

452 ------- 

453 spline : UnivariateSpline 

454 Spline of order k2=k-n representing the derivative of this 

455 spline. 

456 

457 See Also 

458 -------- 

459 splder, antiderivative 

460 

461 Notes 

462 ----- 

463 

464 .. versionadded:: 0.13.0 

465 

466 Examples 

467 -------- 

468 This can be used for finding maxima of a curve: 

469 

470 >>> from scipy.interpolate import UnivariateSpline 

471 >>> x = np.linspace(0, 10, 70) 

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

473 >>> spl = UnivariateSpline(x, y, k=4, s=0) 

474 

475 Now, differentiate the spline and find the zeros of the 

476 derivative. (NB: `sproot` only works for order 3 splines, so we 

477 fit an order 4 spline): 

478 

479 >>> spl.derivative().roots() / np.pi 

480 array([ 0.50000001, 1.5 , 2.49999998]) 

481 

482 This agrees well with roots :math:`\\pi/2 + n\\pi` of 

483 :math:`\\cos(x) = \\sin'(x)`. 

484 

485 """ 

486 tck = fitpack.splder(self._eval_args, n) 

487 # if self.ext is 'const', derivative.ext will be 'zeros' 

488 ext = 1 if self.ext == 3 else self.ext 

489 return UnivariateSpline._from_tck(tck, ext=ext) 

490 

491 def antiderivative(self, n=1): 

492 """ 

493 Construct a new spline representing the antiderivative of this spline. 

494 

495 Parameters 

496 ---------- 

497 n : int, optional 

498 Order of antiderivative to evaluate. Default: 1 

499 

500 Returns 

501 ------- 

502 spline : UnivariateSpline 

503 Spline of order k2=k+n representing the antiderivative of this 

504 spline. 

505 

506 Notes 

507 ----- 

508 

509 .. versionadded:: 0.13.0 

510 

511 See Also 

512 -------- 

513 splantider, derivative 

514 

515 Examples 

516 -------- 

517 >>> from scipy.interpolate import UnivariateSpline 

518 >>> x = np.linspace(0, np.pi/2, 70) 

519 >>> y = 1 / np.sqrt(1 - 0.8*np.sin(x)**2) 

520 >>> spl = UnivariateSpline(x, y, s=0) 

521 

522 The derivative is the inverse operation of the antiderivative, 

523 although some floating point error accumulates: 

524 

525 >>> spl(1.7), spl.antiderivative().derivative()(1.7) 

526 (array(2.1565429877197317), array(2.1565429877201865)) 

527 

528 Antiderivative can be used to evaluate definite integrals: 

529 

530 >>> ispl = spl.antiderivative() 

531 >>> ispl(np.pi/2) - ispl(0) 

532 2.2572053588768486 

533 

534 This is indeed an approximation to the complete elliptic integral 

535 :math:`K(m) = \\int_0^{\\pi/2} [1 - m\\sin^2 x]^{-1/2} dx`: 

536 

537 >>> from scipy.special import ellipk 

538 >>> ellipk(0.8) 

539 2.2572053268208538 

540 

541 """ 

542 tck = fitpack.splantider(self._eval_args, n) 

543 return UnivariateSpline._from_tck(tck, self.ext) 

544 

545 

546class InterpolatedUnivariateSpline(UnivariateSpline): 

547 """ 

548 1-D interpolating spline for a given set of data points. 

549 

550 Fits a spline y = spl(x) of degree `k` to the provided `x`, `y` data. 

551 Spline function passes through all provided points. Equivalent to 

552 `UnivariateSpline` with s=0. 

553 

554 Parameters 

555 ---------- 

556 x : (N,) array_like 

557 Input dimension of data points -- must be strictly increasing 

558 y : (N,) array_like 

559 input dimension of data points 

560 w : (N,) array_like, optional 

561 Weights for spline fitting. Must be positive. If None (default), 

562 weights are all equal. 

563 bbox : (2,) array_like, optional 

564 2-sequence specifying the boundary of the approximation interval. If 

565 None (default), ``bbox=[x[0], x[-1]]``. 

566 k : int, optional 

567 Degree of the smoothing spline. Must be 1 <= `k` <= 5. 

568 ext : int or str, optional 

569 Controls the extrapolation mode for elements 

570 not in the interval defined by the knot sequence. 

571 

572 * if ext=0 or 'extrapolate', return the extrapolated value. 

573 * if ext=1 or 'zeros', return 0 

574 * if ext=2 or 'raise', raise a ValueError 

575 * if ext=3 of 'const', return the boundary value. 

576 

577 The default value is 0. 

578 

579 check_finite : bool, optional 

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

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

582 (crashes, non-termination or non-sensical results) if the inputs 

583 do contain infinities or NaNs. 

584 Default is False. 

585 

586 See Also 

587 -------- 

588 UnivariateSpline : Superclass -- allows knots to be selected by a 

589 smoothing condition 

590 LSQUnivariateSpline : spline for which knots are user-selected 

591 splrep : An older, non object-oriented wrapping of FITPACK 

592 splev, sproot, splint, spalde 

593 BivariateSpline : A similar class for two-dimensional spline interpolation 

594 

595 Notes 

596 ----- 

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

598 

599 Examples 

600 -------- 

601 >>> import matplotlib.pyplot as plt 

602 >>> from scipy.interpolate import InterpolatedUnivariateSpline 

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

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

605 >>> spl = InterpolatedUnivariateSpline(x, y) 

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

607 >>> xs = np.linspace(-3, 3, 1000) 

608 >>> plt.plot(xs, spl(xs), 'g', lw=3, alpha=0.7) 

609 >>> plt.show() 

610 

611 Notice that the ``spl(x)`` interpolates `y`: 

612 

613 >>> spl.get_residual() 

614 0.0 

615 

616 """ 

617 def __init__(self, x, y, w=None, bbox=[None]*2, k=3, 

618 ext=0, check_finite=False): 

619 

620 x, y, w, bbox, self.ext = self.validate_input(x, y, w, bbox, k, None, 

621 ext, check_finite) 

622 if not np.all(diff(x) > 0.0): 

623 raise ValueError('x must be strictly increasing') 

624 

625 # _data == x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier 

626 self._data = dfitpack.fpcurf0(x, y, k, w=w, xb=bbox[0], 

627 xe=bbox[1], s=0) 

628 self._reset_class() 

629 

630 

631_fpchec_error_string = """The input parameters have been rejected by fpchec. \ 

632This means that at least one of the following conditions is violated: 

633 

6341) k+1 <= n-k-1 <= m 

6352) t(1) <= t(2) <= ... <= t(k+1) 

636 t(n-k) <= t(n-k+1) <= ... <= t(n) 

6373) t(k+1) < t(k+2) < ... < t(n-k) 

6384) t(k+1) <= x(i) <= t(n-k) 

6395) The conditions specified by Schoenberg and Whitney must hold 

640 for at least one subset of data points, i.e., there must be a 

641 subset of data points y(j) such that 

642 t(j) < y(j) < t(j+k+1), j=1,2,...,n-k-1 

643""" 

644 

645 

646class LSQUnivariateSpline(UnivariateSpline): 

647 """ 

648 1-D spline with explicit internal knots. 

649 

650 Fits a spline y = spl(x) of degree `k` to the provided `x`, `y` data. `t` 

651 specifies the internal knots of the spline 

652 

653 Parameters 

654 ---------- 

655 x : (N,) array_like 

656 Input dimension of data points -- must be increasing 

657 y : (N,) array_like 

658 Input dimension of data points 

659 t : (M,) array_like 

660 interior knots of the spline. Must be in ascending order and:: 

661 

662 bbox[0] < t[0] < ... < t[-1] < bbox[-1] 

663 

664 w : (N,) array_like, optional 

665 weights for spline fitting. Must be positive. If None (default), 

666 weights are all equal. 

667 bbox : (2,) array_like, optional 

668 2-sequence specifying the boundary of the approximation interval. If 

669 None (default), ``bbox = [x[0], x[-1]]``. 

670 k : int, optional 

671 Degree of the smoothing spline. Must be 1 <= `k` <= 5. 

672 Default is `k` = 3, a cubic spline. 

673 ext : int or str, optional 

674 Controls the extrapolation mode for elements 

675 not in the interval defined by the knot sequence. 

676 

677 * if ext=0 or 'extrapolate', return the extrapolated value. 

678 * if ext=1 or 'zeros', return 0 

679 * if ext=2 or 'raise', raise a ValueError 

680 * if ext=3 of 'const', return the boundary value. 

681 

682 The default value is 0. 

683 

684 check_finite : bool, optional 

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

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

687 (crashes, non-termination or non-sensical results) if the inputs 

688 do contain infinities or NaNs. 

689 Default is False. 

690 

691 Raises 

692 ------ 

693 ValueError 

694 If the interior knots do not satisfy the Schoenberg-Whitney conditions 

695 

696 See Also 

697 -------- 

698 UnivariateSpline : Superclass -- knots are specified by setting a 

699 smoothing condition 

700 InterpolatedUnivariateSpline : spline passing through all points 

701 splrep : An older, non object-oriented wrapping of FITPACK 

702 splev, sproot, splint, spalde 

703 BivariateSpline : A similar class for two-dimensional spline interpolation 

704 

705 Notes 

706 ----- 

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

708 

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

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

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

712 

713 Examples 

714 -------- 

715 >>> from scipy.interpolate import LSQUnivariateSpline, UnivariateSpline 

716 >>> import matplotlib.pyplot as plt 

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

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

719 

720 Fit a smoothing spline with a pre-defined internal knots: 

721 

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

723 >>> spl = LSQUnivariateSpline(x, y, t) 

724 

725 >>> xs = np.linspace(-3, 3, 1000) 

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

727 >>> plt.plot(xs, spl(xs), 'g-', lw=3) 

728 >>> plt.show() 

729 

730 Check the knot vector: 

731 

732 >>> spl.get_knots() 

733 array([-3., -1., 0., 1., 3.]) 

734 

735 Constructing lsq spline using the knots from another spline: 

736 

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

738 >>> s = UnivariateSpline(x, x, s=0) 

739 >>> s.get_knots() 

740 array([ 0., 2., 3., 4., 5., 6., 7., 9.]) 

741 >>> knt = s.get_knots() 

742 >>> s1 = LSQUnivariateSpline(x, x, knt[1:-1]) # Chop 1st and last knot 

743 >>> s1.get_knots() 

744 array([ 0., 2., 3., 4., 5., 6., 7., 9.]) 

745 

746 """ 

747 

748 def __init__(self, x, y, t, w=None, bbox=[None]*2, k=3, 

749 ext=0, check_finite=False): 

750 

751 x, y, w, bbox, self.ext = self.validate_input(x, y, w, bbox, k, None, 

752 ext, check_finite) 

753 if not np.all(diff(x) >= 0.0): 

754 raise ValueError('x must be increasing') 

755 

756 # _data == x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier 

757 xb = bbox[0] 

758 xe = bbox[1] 

759 if xb is None: 

760 xb = x[0] 

761 if xe is None: 

762 xe = x[-1] 

763 t = concatenate(([xb]*(k+1), t, [xe]*(k+1))) 

764 n = len(t) 

765 if not np.all(t[k+1:n-k]-t[k:n-k-1] > 0, axis=0): 

766 raise ValueError('Interior knots t must satisfy ' 

767 'Schoenberg-Whitney conditions') 

768 if not dfitpack.fpchec(x, t, k) == 0: 

769 raise ValueError(_fpchec_error_string) 

770 data = dfitpack.fpcurfm1(x, y, k, t, w=w, xb=xb, xe=xe) 

771 self._data = data[:-3] + (None, None, data[-1]) 

772 self._reset_class() 

773 

774 

775# ############### Bivariate spline #################### 

776 

777class _BivariateSplineBase(object): 

778 """ Base class for Bivariate spline s(x,y) interpolation on the rectangle 

779 [xb,xe] x [yb, ye] calculated from a given set of data points 

780 (x,y,z). 

781 

782 See Also 

783 -------- 

784 bisplrep, bisplev : an older wrapping of FITPACK 

785 BivariateSpline : 

786 implementation of bivariate spline interpolation on a plane grid 

787 SphereBivariateSpline : 

788 implementation of bivariate spline interpolation on a spherical grid 

789 """ 

790 

791 def get_residual(self): 

792 """ Return weighted sum of squared residuals of the spline 

793 approximation: sum ((w[i]*(z[i]-s(x[i],y[i])))**2,axis=0) 

794 """ 

795 return self.fp 

796 

797 def get_knots(self): 

798 """ Return a tuple (tx,ty) where tx,ty contain knots positions 

799 of the spline with respect to x-, y-variable, respectively. 

800 The position of interior and additional knots are given as 

801 t[k+1:-k-1] and t[:k+1]=b, t[-k-1:]=e, respectively. 

802 """ 

803 return self.tck[:2] 

804 

805 def get_coeffs(self): 

806 """ Return spline coefficients.""" 

807 return self.tck[2] 

808 

809 def __call__(self, x, y, dx=0, dy=0, grid=True): 

810 """ 

811 Evaluate the spline or its derivatives at given positions. 

812 

813 Parameters 

814 ---------- 

815 x, y : array_like 

816 Input coordinates. 

817 

818 If `grid` is False, evaluate the spline at points ``(x[i], 

819 y[i]), i=0, ..., len(x)-1``. Standard Numpy broadcasting 

820 is obeyed. 

821 

822 If `grid` is True: evaluate spline at the grid points 

823 defined by the coordinate arrays x, y. The arrays must be 

824 sorted to increasing order. 

825 

826 Note that the axis ordering is inverted relative to 

827 the output of meshgrid. 

828 dx : int 

829 Order of x-derivative 

830 

831 .. versionadded:: 0.14.0 

832 dy : int 

833 Order of y-derivative 

834 

835 .. versionadded:: 0.14.0 

836 grid : bool 

837 Whether to evaluate the results on a grid spanned by the 

838 input arrays, or at points specified by the input arrays. 

839 

840 .. versionadded:: 0.14.0 

841 

842 """ 

843 x = np.asarray(x) 

844 y = np.asarray(y) 

845 

846 tx, ty, c = self.tck[:3] 

847 kx, ky = self.degrees 

848 if grid: 

849 if x.size == 0 or y.size == 0: 

850 return np.zeros((x.size, y.size), dtype=self.tck[2].dtype) 

851 

852 if dx or dy: 

853 z, ier = dfitpack.parder(tx, ty, c, kx, ky, dx, dy, x, y) 

854 if not ier == 0: 

855 raise ValueError("Error code returned by parder: %s" % ier) 

856 else: 

857 z, ier = dfitpack.bispev(tx, ty, c, kx, ky, x, y) 

858 if not ier == 0: 

859 raise ValueError("Error code returned by bispev: %s" % ier) 

860 else: 

861 # standard Numpy broadcasting 

862 if x.shape != y.shape: 

863 x, y = np.broadcast_arrays(x, y) 

864 

865 shape = x.shape 

866 x = x.ravel() 

867 y = y.ravel() 

868 

869 if x.size == 0 or y.size == 0: 

870 return np.zeros(shape, dtype=self.tck[2].dtype) 

871 

872 if dx or dy: 

873 z, ier = dfitpack.pardeu(tx, ty, c, kx, ky, dx, dy, x, y) 

874 if not ier == 0: 

875 raise ValueError("Error code returned by pardeu: %s" % ier) 

876 else: 

877 z, ier = dfitpack.bispeu(tx, ty, c, kx, ky, x, y) 

878 if not ier == 0: 

879 raise ValueError("Error code returned by bispeu: %s" % ier) 

880 

881 z = z.reshape(shape) 

882 return z 

883 

884 

885_surfit_messages = {1: """ 

886The required storage space exceeds the available storage space: nxest 

887or nyest too small, or s too small. 

888The weighted least-squares spline corresponds to the current set of 

889knots.""", 

890 2: """ 

891A theoretically impossible result was found during the iteration 

892process for finding a smoothing spline with fp = s: s too small or 

893badly chosen eps. 

894Weighted sum of squared residuals does not satisfy abs(fp-s)/s < tol.""", 

895 3: """ 

896the maximal number of iterations maxit (set to 20 by the program) 

897allowed for finding a smoothing spline with fp=s has been reached: 

898s too small. 

899Weighted sum of squared residuals does not satisfy abs(fp-s)/s < tol.""", 

900 4: """ 

901No more knots can be added because the number of b-spline coefficients 

902(nx-kx-1)*(ny-ky-1) already exceeds the number of data points m: 

903either s or m too small. 

904The weighted least-squares spline corresponds to the current set of 

905knots.""", 

906 5: """ 

907No more knots can be added because the additional knot would (quasi) 

908coincide with an old one: s too small or too large a weight to an 

909inaccurate data point. 

910The weighted least-squares spline corresponds to the current set of 

911knots.""", 

912 10: """ 

913Error on entry, no approximation returned. The following conditions 

914must hold: 

915xb<=x[i]<=xe, yb<=y[i]<=ye, w[i]>0, i=0..m-1 

916If iopt==-1, then 

917 xb<tx[kx+1]<tx[kx+2]<...<tx[nx-kx-2]<xe 

918 yb<ty[ky+1]<ty[ky+2]<...<ty[ny-ky-2]<ye""", 

919 -3: """ 

920The coefficients of the spline returned have been computed as the 

921minimal norm least-squares solution of a (numerically) rank deficient 

922system (deficiency=%i). If deficiency is large, the results may be 

923inaccurate. Deficiency may strongly depend on the value of eps.""" 

924 } 

925 

926 

927class BivariateSpline(_BivariateSplineBase): 

928 """ 

929 Base class for bivariate splines. 

930 

931 This describes a spline ``s(x, y)`` of degrees ``kx`` and ``ky`` on 

932 the rectangle ``[xb, xe] * [yb, ye]`` calculated from a given set 

933 of data points ``(x, y, z)``. 

934 

935 This class is meant to be subclassed, not instantiated directly. 

936 To construct these splines, call either `SmoothBivariateSpline` or 

937 `LSQBivariateSpline`. 

938 

939 See Also 

940 -------- 

941 UnivariateSpline : 

942 a similar class for univariate spline interpolation 

943 SmoothBivariateSpline : 

944 to create a bivariate spline through the given points 

945 LSQBivariateSpline : 

946 to create a bivariate spline using weighted least-squares fitting 

947 RectSphereBivariateSpline : 

948 to create a bivariate spline over a rectangular mesh on a sphere 

949 SmoothSphereBivariateSpline : 

950 to create a smooth bivariate spline in spherical coordinates 

951 LSQSphereBivariateSpline : 

952 to create a bivariate spline in spherical coordinates using 

953 weighted least-squares fitting 

954 bisplrep : older wrapping of FITPACK 

955 bisplev : older wrapping of FITPACK 

956 """ 

957 

958 @classmethod 

959 def _from_tck(cls, tck): 

960 """Construct a spline object from given tck and degree""" 

961 self = cls.__new__(cls) 

962 if len(tck) != 5: 

963 raise ValueError("tck should be a 5 element tuple of tx," 

964 " ty, c, kx, ky") 

965 self.tck = tck[:3] 

966 self.degrees = tck[3:] 

967 return self 

968 

969 def ev(self, xi, yi, dx=0, dy=0): 

970 """ 

971 Evaluate the spline at points 

972 

973 Returns the interpolated value at ``(xi[i], yi[i]), 

974 i=0,...,len(xi)-1``. 

975 

976 Parameters 

977 ---------- 

978 xi, yi : array_like 

979 Input coordinates. Standard Numpy broadcasting is obeyed. 

980 dx : int, optional 

981 Order of x-derivative 

982 

983 .. versionadded:: 0.14.0 

984 dy : int, optional 

985 Order of y-derivative 

986 

987 .. versionadded:: 0.14.0 

988 """ 

989 return self.__call__(xi, yi, dx=dx, dy=dy, grid=False) 

990 

991 def integral(self, xa, xb, ya, yb): 

992 """ 

993 Evaluate the integral of the spline over area [xa,xb] x [ya,yb]. 

994 

995 Parameters 

996 ---------- 

997 xa, xb : float 

998 The end-points of the x integration interval. 

999 ya, yb : float 

1000 The end-points of the y integration interval. 

1001 

1002 Returns 

1003 ------- 

1004 integ : float 

1005 The value of the resulting integral. 

1006 

1007 """ 

1008 tx, ty, c = self.tck[:3] 

1009 kx, ky = self.degrees 

1010 return dfitpack.dblint(tx, ty, c, kx, ky, xa, xb, ya, yb) 

1011 

1012 @staticmethod 

1013 def _validate_input(x, y, z, w, kx, ky, eps): 

1014 x, y, z = np.asarray(x), np.asarray(y), np.asarray(z) 

1015 if not x.size == y.size == z.size: 

1016 raise ValueError('x, y, and z should have a same length') 

1017 

1018 if w is not None: 

1019 w = np.asarray(w) 

1020 if x.size != w.size: 

1021 raise ValueError('x, y, z, and w should have a same length') 

1022 elif not np.all(w >= 0.0): 

1023 raise ValueError('w should be positive') 

1024 if (eps is not None) and (not 0.0 < eps < 1.0): 

1025 raise ValueError('eps should be between (0, 1)') 

1026 if not x.size >= (kx + 1) * (ky + 1): 

1027 raise ValueError('The length of x, y and z should be at least' 

1028 ' (kx+1) * (ky+1)') 

1029 return x, y, z, w 

1030 

1031 

1032class SmoothBivariateSpline(BivariateSpline): 

1033 """ 

1034 Smooth bivariate spline approximation. 

1035 

1036 Parameters 

1037 ---------- 

1038 x, y, z : array_like 

1039 1-D sequences of data points (order is not important). 

1040 w : array_like, optional 

1041 Positive 1-D sequence of weights, of same length as `x`, `y` and `z`. 

1042 bbox : array_like, optional 

1043 Sequence of length 4 specifying the boundary of the rectangular 

1044 approximation domain. By default, 

1045 ``bbox=[min(x), max(x), min(y), max(y)]``. 

1046 kx, ky : ints, optional 

1047 Degrees of the bivariate spline. Default is 3. 

1048 s : float, optional 

1049 Positive smoothing factor defined for estimation condition: 

1050 ``sum((w[i]*(z[i]-s(x[i], y[i])))**2, axis=0) <= s`` 

1051 Default ``s=len(w)`` which should be a good value if ``1/w[i]`` is an 

1052 estimate of the standard deviation of ``z[i]``. 

1053 eps : float, optional 

1054 A threshold for determining the effective rank of an over-determined 

1055 linear system of equations. `eps` should have a value within the open 

1056 interval ``(0, 1)``, the default is 1e-16. 

1057 

1058 See Also 

1059 -------- 

1060 bisplrep : an older wrapping of FITPACK 

1061 bisplev : an older wrapping of FITPACK 

1062 UnivariateSpline : a similar class for univariate spline interpolation 

1063 LSQBivariateSpline : to create a BivariateSpline using weighted least-squares fitting 

1064 

1065 Notes 

1066 ----- 

1067 The length of `x`, `y` and `z` should be at least ``(kx+1) * (ky+1)``. 

1068 

1069 """ 

1070 

1071 def __init__(self, x, y, z, w=None, bbox=[None] * 4, kx=3, ky=3, s=None, 

1072 eps=1e-16): 

1073 

1074 x, y, z, w = self._validate_input(x, y, z, w, kx, ky, eps) 

1075 bbox = ravel(bbox) 

1076 if not bbox.shape == (4,): 

1077 raise ValueError('bbox shape should be (4,)') 

1078 if s is not None and not s >= 0.0: 

1079 raise ValueError("s should be s >= 0.0") 

1080 

1081 xb, xe, yb, ye = bbox 

1082 nx, tx, ny, ty, c, fp, wrk1, ier = dfitpack.surfit_smth(x, y, z, w, 

1083 xb, xe, yb, 

1084 ye, kx, ky, 

1085 s=s, eps=eps, 

1086 lwrk2=1) 

1087 if ier > 10: # lwrk2 was to small, re-run 

1088 nx, tx, ny, ty, c, fp, wrk1, ier = dfitpack.surfit_smth(x, y, z, w, 

1089 xb, xe, yb, 

1090 ye, kx, ky, 

1091 s=s, 

1092 eps=eps, 

1093 lwrk2=ier) 

1094 if ier in [0, -1, -2]: # normal return 

1095 pass 

1096 else: 

1097 message = _surfit_messages.get(ier, 'ier=%s' % (ier)) 

1098 warnings.warn(message) 

1099 

1100 self.fp = fp 

1101 self.tck = tx[:nx], ty[:ny], c[:(nx-kx-1)*(ny-ky-1)] 

1102 self.degrees = kx, ky 

1103 

1104 

1105class LSQBivariateSpline(BivariateSpline): 

1106 """ 

1107 Weighted least-squares bivariate spline approximation. 

1108 

1109 Parameters 

1110 ---------- 

1111 x, y, z : array_like 

1112 1-D sequences of data points (order is not important). 

1113 tx, ty : array_like 

1114 Strictly ordered 1-D sequences of knots coordinates. 

1115 w : array_like, optional 

1116 Positive 1-D array of weights, of the same length as `x`, `y` and `z`. 

1117 bbox : (4,) array_like, optional 

1118 Sequence of length 4 specifying the boundary of the rectangular 

1119 approximation domain. By default, 

1120 ``bbox=[min(x,tx),max(x,tx), min(y,ty),max(y,ty)]``. 

1121 kx, ky : ints, optional 

1122 Degrees of the bivariate spline. Default is 3. 

1123 eps : float, optional 

1124 A threshold for determining the effective rank of an over-determined 

1125 linear system of equations. `eps` should have a value within the open 

1126 interval ``(0, 1)``, the default is 1e-16. 

1127 

1128 See Also 

1129 -------- 

1130 bisplrep : an older wrapping of FITPACK 

1131 bisplev : an older wrapping of FITPACK 

1132 UnivariateSpline : a similar class for univariate spline interpolation 

1133 SmoothBivariateSpline : create a smoothing BivariateSpline 

1134 

1135 Notes 

1136 ----- 

1137 The length of `x`, `y` and `z` should be at least ``(kx+1) * (ky+1)``. 

1138 

1139 """ 

1140 

1141 def __init__(self, x, y, z, tx, ty, w=None, bbox=[None]*4, kx=3, ky=3, 

1142 eps=None): 

1143 

1144 x, y, z, w = self._validate_input(x, y, z, w, kx, ky, eps) 

1145 bbox = ravel(bbox) 

1146 if not bbox.shape == (4,): 

1147 raise ValueError('bbox shape should be (4,)') 

1148 

1149 nx = 2*kx+2+len(tx) 

1150 ny = 2*ky+2+len(ty) 

1151 tx1 = zeros((nx,), float) 

1152 ty1 = zeros((ny,), float) 

1153 tx1[kx+1:nx-kx-1] = tx 

1154 ty1[ky+1:ny-ky-1] = ty 

1155 

1156 xb, xe, yb, ye = bbox 

1157 tx1, ty1, c, fp, ier = dfitpack.surfit_lsq(x, y, z, tx1, ty1, w, 

1158 xb, xe, yb, ye, 

1159 kx, ky, eps, lwrk2=1) 

1160 if ier > 10: 

1161 tx1, ty1, c, fp, ier = dfitpack.surfit_lsq(x, y, z, tx1, ty1, w, 

1162 xb, xe, yb, ye, 

1163 kx, ky, eps, lwrk2=ier) 

1164 if ier in [0, -1, -2]: # normal return 

1165 pass 

1166 else: 

1167 if ier < -2: 

1168 deficiency = (nx-kx-1)*(ny-ky-1)+ier 

1169 message = _surfit_messages.get(-3) % (deficiency) 

1170 else: 

1171 message = _surfit_messages.get(ier, 'ier=%s' % (ier)) 

1172 warnings.warn(message) 

1173 self.fp = fp 

1174 self.tck = tx1, ty1, c 

1175 self.degrees = kx, ky 

1176 

1177 

1178class RectBivariateSpline(BivariateSpline): 

1179 """ 

1180 Bivariate spline approximation over a rectangular mesh. 

1181 

1182 Can be used for both smoothing and interpolating data. 

1183 

1184 Parameters 

1185 ---------- 

1186 x,y : array_like 

1187 1-D arrays of coordinates in strictly ascending order. 

1188 z : array_like 

1189 2-D array of data with shape (x.size,y.size). 

1190 bbox : array_like, optional 

1191 Sequence of length 4 specifying the boundary of the rectangular 

1192 approximation domain. By default, 

1193 ``bbox=[min(x,tx),max(x,tx), min(y,ty),max(y,ty)]``. 

1194 kx, ky : ints, optional 

1195 Degrees of the bivariate spline. Default is 3. 

1196 s : float, optional 

1197 Positive smoothing factor defined for estimation condition: 

1198 ``sum((w[i]*(z[i]-s(x[i], y[i])))**2, axis=0) <= s`` 

1199 Default is ``s=0``, which is for interpolation. 

1200 

1201 See Also 

1202 -------- 

1203 SmoothBivariateSpline : a smoothing bivariate spline for scattered data 

1204 bisplrep : an older wrapping of FITPACK 

1205 bisplev : an older wrapping of FITPACK 

1206 UnivariateSpline : a similar class for univariate spline interpolation 

1207 

1208 """ 

1209 

1210 def __init__(self, x, y, z, bbox=[None] * 4, kx=3, ky=3, s=0): 

1211 x, y, bbox = ravel(x), ravel(y), ravel(bbox) 

1212 z = np.asarray(z) 

1213 if not np.all(diff(x) > 0.0): 

1214 raise ValueError('x must be strictly increasing') 

1215 if not np.all(diff(y) > 0.0): 

1216 raise ValueError('y must be strictly increasing') 

1217 if not x.size == z.shape[0]: 

1218 raise ValueError('x dimension of z must have same number of ' 

1219 'elements as x') 

1220 if not y.size == z.shape[1]: 

1221 raise ValueError('y dimension of z must have same number of ' 

1222 'elements as y') 

1223 if not bbox.shape == (4,): 

1224 raise ValueError('bbox shape should be (4,)') 

1225 if s is not None and not s >= 0.0: 

1226 raise ValueError("s should be s >= 0.0") 

1227 

1228 z = ravel(z) 

1229 xb, xe, yb, ye = bbox 

1230 nx, tx, ny, ty, c, fp, ier = dfitpack.regrid_smth(x, y, z, xb, xe, yb, 

1231 ye, kx, ky, s) 

1232 

1233 if ier not in [0, -1, -2]: 

1234 msg = _surfit_messages.get(ier, 'ier=%s' % (ier)) 

1235 raise ValueError(msg) 

1236 

1237 self.fp = fp 

1238 self.tck = tx[:nx], ty[:ny], c[:(nx - kx - 1) * (ny - ky - 1)] 

1239 self.degrees = kx, ky 

1240 

1241 

1242_spherefit_messages = _surfit_messages.copy() 

1243_spherefit_messages[10] = """ 

1244ERROR. On entry, the input data are controlled on validity. The following 

1245 restrictions must be satisfied: 

1246 -1<=iopt<=1, m>=2, ntest>=8 ,npest >=8, 0<eps<1, 

1247 0<=teta(i)<=pi, 0<=phi(i)<=2*pi, w(i)>0, i=1,...,m 

1248 lwrk1 >= 185+52*v+10*u+14*u*v+8*(u-1)*v**2+8*m 

1249 kwrk >= m+(ntest-7)*(npest-7) 

1250 if iopt=-1: 8<=nt<=ntest , 9<=np<=npest 

1251 0<tt(5)<tt(6)<...<tt(nt-4)<pi 

1252 0<tp(5)<tp(6)<...<tp(np-4)<2*pi 

1253 if iopt>=0: s>=0 

1254 if one of these conditions is found to be violated,control 

1255 is immediately repassed to the calling program. in that 

1256 case there is no approximation returned.""" 

1257_spherefit_messages[-3] = """ 

1258WARNING. The coefficients of the spline returned have been computed as the 

1259 minimal norm least-squares solution of a (numerically) rank 

1260 deficient system (deficiency=%i, rank=%i). Especially if the rank 

1261 deficiency, which is computed by 6+(nt-8)*(np-7)+ier, is large, 

1262 the results may be inaccurate. They could also seriously depend on 

1263 the value of eps.""" 

1264 

1265 

1266class SphereBivariateSpline(_BivariateSplineBase): 

1267 """ 

1268 Bivariate spline s(x,y) of degrees 3 on a sphere, calculated from a 

1269 given set of data points (theta,phi,r). 

1270 

1271 .. versionadded:: 0.11.0 

1272 

1273 See Also 

1274 -------- 

1275 bisplrep, bisplev : an older wrapping of FITPACK 

1276 UnivariateSpline : a similar class for univariate spline interpolation 

1277 SmoothUnivariateSpline : 

1278 to create a BivariateSpline through the given points 

1279 LSQUnivariateSpline : 

1280 to create a BivariateSpline using weighted least-squares fitting 

1281 """ 

1282 

1283 def __call__(self, theta, phi, dtheta=0, dphi=0, grid=True): 

1284 """ 

1285 Evaluate the spline or its derivatives at given positions. 

1286 

1287 Parameters 

1288 ---------- 

1289 theta, phi : array_like 

1290 Input coordinates. 

1291 

1292 If `grid` is False, evaluate the spline at points 

1293 ``(theta[i], phi[i]), i=0, ..., len(x)-1``. Standard 

1294 Numpy broadcasting is obeyed. 

1295 

1296 If `grid` is True: evaluate spline at the grid points 

1297 defined by the coordinate arrays theta, phi. The arrays 

1298 must be sorted to increasing order. 

1299 dtheta : int, optional 

1300 Order of theta-derivative 

1301 

1302 .. versionadded:: 0.14.0 

1303 dphi : int 

1304 Order of phi-derivative 

1305 

1306 .. versionadded:: 0.14.0 

1307 grid : bool 

1308 Whether to evaluate the results on a grid spanned by the 

1309 input arrays, or at points specified by the input arrays. 

1310 

1311 .. versionadded:: 0.14.0 

1312 

1313 """ 

1314 theta = np.asarray(theta) 

1315 phi = np.asarray(phi) 

1316 

1317 if theta.size > 0 and (theta.min() < 0. or theta.max() > np.pi): 

1318 raise ValueError("requested theta out of bounds.") 

1319 if phi.size > 0 and (phi.min() < 0. or phi.max() > 2. * np.pi): 

1320 raise ValueError("requested phi out of bounds.") 

1321 

1322 return _BivariateSplineBase.__call__(self, theta, phi, 

1323 dx=dtheta, dy=dphi, grid=grid) 

1324 

1325 def ev(self, theta, phi, dtheta=0, dphi=0): 

1326 """ 

1327 Evaluate the spline at points 

1328 

1329 Returns the interpolated value at ``(theta[i], phi[i]), 

1330 i=0,...,len(theta)-1``. 

1331 

1332 Parameters 

1333 ---------- 

1334 theta, phi : array_like 

1335 Input coordinates. Standard Numpy broadcasting is obeyed. 

1336 dtheta : int, optional 

1337 Order of theta-derivative 

1338 

1339 .. versionadded:: 0.14.0 

1340 dphi : int, optional 

1341 Order of phi-derivative 

1342 

1343 .. versionadded:: 0.14.0 

1344 """ 

1345 return self.__call__(theta, phi, dtheta=dtheta, dphi=dphi, grid=False) 

1346 

1347 

1348class SmoothSphereBivariateSpline(SphereBivariateSpline): 

1349 """ 

1350 Smooth bivariate spline approximation in spherical coordinates. 

1351 

1352 .. versionadded:: 0.11.0 

1353 

1354 Parameters 

1355 ---------- 

1356 theta, phi, r : array_like 

1357 1-D sequences of data points (order is not important). Coordinates 

1358 must be given in radians. Theta must lie within the interval 

1359 ``[0, pi]``, and phi must lie within the interval ``[0, 2pi]``. 

1360 w : array_like, optional 

1361 Positive 1-D sequence of weights. 

1362 s : float, optional 

1363 Positive smoothing factor defined for estimation condition: 

1364 ``sum((w(i)*(r(i) - s(theta(i), phi(i))))**2, axis=0) <= s`` 

1365 Default ``s=len(w)`` which should be a good value if ``1/w[i]`` is an 

1366 estimate of the standard deviation of ``r[i]``. 

1367 eps : float, optional 

1368 A threshold for determining the effective rank of an over-determined 

1369 linear system of equations. `eps` should have a value within the open 

1370 interval ``(0, 1)``, the default is 1e-16. 

1371 

1372 Notes 

1373 ----- 

1374 For more information, see the FITPACK_ site about this function. 

1375 

1376 .. _FITPACK: http://www.netlib.org/dierckx/sphere.f 

1377 

1378 Examples 

1379 -------- 

1380 Suppose we have global data on a coarse grid (the input data does not 

1381 have to be on a grid): 

1382 

1383 >>> theta = np.linspace(0., np.pi, 7) 

1384 >>> phi = np.linspace(0., 2*np.pi, 9) 

1385 >>> data = np.empty((theta.shape[0], phi.shape[0])) 

1386 >>> data[:,0], data[0,:], data[-1,:] = 0., 0., 0. 

1387 >>> data[1:-1,1], data[1:-1,-1] = 1., 1. 

1388 >>> data[1,1:-1], data[-2,1:-1] = 1., 1. 

1389 >>> data[2:-2,2], data[2:-2,-2] = 2., 2. 

1390 >>> data[2,2:-2], data[-3,2:-2] = 2., 2. 

1391 >>> data[3,3:-2] = 3. 

1392 >>> data = np.roll(data, 4, 1) 

1393 

1394 We need to set up the interpolator object 

1395 

1396 >>> lats, lons = np.meshgrid(theta, phi) 

1397 >>> from scipy.interpolate import SmoothSphereBivariateSpline 

1398 >>> lut = SmoothSphereBivariateSpline(lats.ravel(), lons.ravel(), 

1399 ... data.T.ravel(), s=3.5) 

1400 

1401 As a first test, we'll see what the algorithm returns when run on the 

1402 input coordinates 

1403 

1404 >>> data_orig = lut(theta, phi) 

1405 

1406 Finally we interpolate the data to a finer grid 

1407 

1408 >>> fine_lats = np.linspace(0., np.pi, 70) 

1409 >>> fine_lons = np.linspace(0., 2 * np.pi, 90) 

1410 

1411 >>> data_smth = lut(fine_lats, fine_lons) 

1412 

1413 >>> import matplotlib.pyplot as plt 

1414 >>> fig = plt.figure() 

1415 >>> ax1 = fig.add_subplot(131) 

1416 >>> ax1.imshow(data, interpolation='nearest') 

1417 >>> ax2 = fig.add_subplot(132) 

1418 >>> ax2.imshow(data_orig, interpolation='nearest') 

1419 >>> ax3 = fig.add_subplot(133) 

1420 >>> ax3.imshow(data_smth, interpolation='nearest') 

1421 >>> plt.show() 

1422 

1423 """ 

1424 

1425 def __init__(self, theta, phi, r, w=None, s=0., eps=1E-16): 

1426 

1427 # input validation 

1428 if not ((0.0 <= theta).all() and (theta <= np.pi).all()): 

1429 raise ValueError('theta should be between [0, pi]') 

1430 if not ((0.0 <= phi).all() and (phi <= 2.0 * np.pi).all()): 

1431 raise ValueError('phi should be between [0, 2pi]') 

1432 if w is not None and not (w >= 0.0).all(): 

1433 raise ValueError('w should be positive') 

1434 if not s >= 0.0: 

1435 raise ValueError('s should be positive') 

1436 if not 0.0 < eps < 1.0: 

1437 raise ValueError('eps should be between (0, 1)') 

1438 

1439 if np.issubclass_(w, float): 

1440 w = ones(len(theta)) * w 

1441 nt_, tt_, np_, tp_, c, fp, ier = dfitpack.spherfit_smth(theta, phi, 

1442 r, w=w, s=s, 

1443 eps=eps) 

1444 if ier not in [0, -1, -2]: 

1445 message = _spherefit_messages.get(ier, 'ier=%s' % (ier)) 

1446 raise ValueError(message) 

1447 

1448 self.fp = fp 

1449 self.tck = tt_[:nt_], tp_[:np_], c[:(nt_ - 4) * (np_ - 4)] 

1450 self.degrees = (3, 3) 

1451 

1452 

1453class LSQSphereBivariateSpline(SphereBivariateSpline): 

1454 """ 

1455 Weighted least-squares bivariate spline approximation in spherical 

1456 coordinates. 

1457 

1458 Determines a smooth bicubic spline according to a given 

1459 set of knots in the `theta` and `phi` directions. 

1460 

1461 .. versionadded:: 0.11.0 

1462 

1463 Parameters 

1464 ---------- 

1465 theta, phi, r : array_like 

1466 1-D sequences of data points (order is not important). Coordinates 

1467 must be given in radians. Theta must lie within the interval 

1468 ``[0, pi]``, and phi must lie within the interval ``[0, 2pi]``. 

1469 tt, tp : array_like 

1470 Strictly ordered 1-D sequences of knots coordinates. 

1471 Coordinates must satisfy ``0 < tt[i] < pi``, ``0 < tp[i] < 2*pi``. 

1472 w : array_like, optional 

1473 Positive 1-D sequence of weights, of the same length as `theta`, `phi` 

1474 and `r`. 

1475 eps : float, optional 

1476 A threshold for determining the effective rank of an over-determined 

1477 linear system of equations. `eps` should have a value within the 

1478 open interval ``(0, 1)``, the default is 1e-16. 

1479 

1480 Notes 

1481 ----- 

1482 For more information, see the FITPACK_ site about this function. 

1483 

1484 .. _FITPACK: http://www.netlib.org/dierckx/sphere.f 

1485 

1486 Examples 

1487 -------- 

1488 Suppose we have global data on a coarse grid (the input data does not 

1489 have to be on a grid): 

1490 

1491 >>> theta = np.linspace(0., np.pi, 7) 

1492 >>> phi = np.linspace(0., 2*np.pi, 9) 

1493 >>> data = np.empty((theta.shape[0], phi.shape[0])) 

1494 >>> data[:,0], data[0,:], data[-1,:] = 0., 0., 0. 

1495 >>> data[1:-1,1], data[1:-1,-1] = 1., 1. 

1496 >>> data[1,1:-1], data[-2,1:-1] = 1., 1. 

1497 >>> data[2:-2,2], data[2:-2,-2] = 2., 2. 

1498 >>> data[2,2:-2], data[-3,2:-2] = 2., 2. 

1499 >>> data[3,3:-2] = 3. 

1500 >>> data = np.roll(data, 4, 1) 

1501 

1502 We need to set up the interpolator object. Here, we must also specify the 

1503 coordinates of the knots to use. 

1504 

1505 >>> lats, lons = np.meshgrid(theta, phi) 

1506 >>> knotst, knotsp = theta.copy(), phi.copy() 

1507 >>> knotst[0] += .0001 

1508 >>> knotst[-1] -= .0001 

1509 >>> knotsp[0] += .0001 

1510 >>> knotsp[-1] -= .0001 

1511 >>> from scipy.interpolate import LSQSphereBivariateSpline 

1512 >>> lut = LSQSphereBivariateSpline(lats.ravel(), lons.ravel(), 

1513 ... data.T.ravel(), knotst, knotsp) 

1514 

1515 As a first test, we'll see what the algorithm returns when run on the 

1516 input coordinates 

1517 

1518 >>> data_orig = lut(theta, phi) 

1519 

1520 Finally we interpolate the data to a finer grid 

1521 

1522 >>> fine_lats = np.linspace(0., np.pi, 70) 

1523 >>> fine_lons = np.linspace(0., 2*np.pi, 90) 

1524 

1525 >>> data_lsq = lut(fine_lats, fine_lons) 

1526 

1527 >>> import matplotlib.pyplot as plt 

1528 >>> fig = plt.figure() 

1529 >>> ax1 = fig.add_subplot(131) 

1530 >>> ax1.imshow(data, interpolation='nearest') 

1531 >>> ax2 = fig.add_subplot(132) 

1532 >>> ax2.imshow(data_orig, interpolation='nearest') 

1533 >>> ax3 = fig.add_subplot(133) 

1534 >>> ax3.imshow(data_lsq, interpolation='nearest') 

1535 >>> plt.show() 

1536 

1537 """ 

1538 

1539 def __init__(self, theta, phi, r, tt, tp, w=None, eps=1E-16): 

1540 

1541 if not ((0.0 <= theta).all() and (theta <= np.pi).all()): 

1542 raise ValueError('theta should be between [0, pi]') 

1543 if not ((0.0 <= phi).all() and (phi <= 2*np.pi).all()): 

1544 raise ValueError('phi should be between [0, 2pi]') 

1545 if not ((0.0 < tt).all() and (tt < np.pi).all()): 

1546 raise ValueError('tt should be between (0, pi)') 

1547 if not ((0.0 < tp).all() and (tp < 2*np.pi).all()): 

1548 raise ValueError('tp should be between (0, 2pi)') 

1549 if w is not None and not (w >= 0.0).all(): 

1550 raise ValueError('w should be positive') 

1551 if not 0.0 < eps < 1.0: 

1552 raise ValueError('eps should be between (0, 1)') 

1553 

1554 if np.issubclass_(w, float): 

1555 w = ones(len(theta)) * w 

1556 nt_, np_ = 8 + len(tt), 8 + len(tp) 

1557 tt_, tp_ = zeros((nt_,), float), zeros((np_,), float) 

1558 tt_[4:-4], tp_[4:-4] = tt, tp 

1559 tt_[-4:], tp_[-4:] = np.pi, 2. * np.pi 

1560 tt_, tp_, c, fp, ier = dfitpack.spherfit_lsq(theta, phi, r, tt_, tp_, 

1561 w=w, eps=eps) 

1562 if ier < -2: 

1563 deficiency = 6 + (nt_ - 8) * (np_ - 7) + ier 

1564 message = _spherefit_messages.get(-3) % (deficiency, -ier) 

1565 warnings.warn(message, stacklevel=2) 

1566 elif ier not in [0, -1, -2]: 

1567 message = _spherefit_messages.get(ier, 'ier=%s' % (ier)) 

1568 raise ValueError(message) 

1569 

1570 self.fp = fp 

1571 self.tck = tt_, tp_, c 

1572 self.degrees = (3, 3) 

1573 

1574 

1575_spfit_messages = _surfit_messages.copy() 

1576_spfit_messages[10] = """ 

1577ERROR: on entry, the input data are controlled on validity 

1578 the following restrictions must be satisfied. 

1579 -1<=iopt(1)<=1, 0<=iopt(2)<=1, 0<=iopt(3)<=1, 

1580 -1<=ider(1)<=1, 0<=ider(2)<=1, ider(2)=0 if iopt(2)=0. 

1581 -1<=ider(3)<=1, 0<=ider(4)<=1, ider(4)=0 if iopt(3)=0. 

1582 mu >= mumin (see above), mv >= 4, nuest >=8, nvest >= 8, 

1583 kwrk>=5+mu+mv+nuest+nvest, 

1584 lwrk >= 12+nuest*(mv+nvest+3)+nvest*24+4*mu+8*mv+max(nuest,mv+nvest) 

1585 0< u(i-1)<u(i)< pi,i=2,..,mu, 

1586 -pi<=v(1)< pi, v(1)<v(i-1)<v(i)<v(1)+2*pi, i=3,...,mv 

1587 if iopt(1)=-1: 8<=nu<=min(nuest,mu+6+iopt(2)+iopt(3)) 

1588 0<tu(5)<tu(6)<...<tu(nu-4)< pi 

1589 8<=nv<=min(nvest,mv+7) 

1590 v(1)<tv(5)<tv(6)<...<tv(nv-4)<v(1)+2*pi 

1591 the schoenberg-whitney conditions, i.e. there must be 

1592 subset of grid co-ordinates uu(p) and vv(q) such that 

1593 tu(p) < uu(p) < tu(p+4) ,p=1,...,nu-4 

1594 (iopt(2)=1 and iopt(3)=1 also count for a uu-value 

1595 tv(q) < vv(q) < tv(q+4) ,q=1,...,nv-4 

1596 (vv(q) is either a value v(j) or v(j)+2*pi) 

1597 if iopt(1)>=0: s>=0 

1598 if s=0: nuest>=mu+6+iopt(2)+iopt(3), nvest>=mv+7 

1599 if one of these conditions is found to be violated,control is 

1600 immediately repassed to the calling program. in that case there is no 

1601 approximation returned.""" 

1602 

1603 

1604class RectSphereBivariateSpline(SphereBivariateSpline): 

1605 """ 

1606 Bivariate spline approximation over a rectangular mesh on a sphere. 

1607 

1608 Can be used for smoothing data. 

1609 

1610 .. versionadded:: 0.11.0 

1611 

1612 Parameters 

1613 ---------- 

1614 u : array_like 

1615 1-D array of colatitude coordinates in strictly ascending order. 

1616 Coordinates must be given in radians and lie within the interval 

1617 ``[0, pi]``. 

1618 v : array_like 

1619 1-D array of longitude coordinates in strictly ascending order. 

1620 Coordinates must be given in radians. First element (``v[0]``) must lie 

1621 within the interval ``[-pi, pi)``. Last element (``v[-1]``) must satisfy 

1622 ``v[-1] <= v[0] + 2*pi``. 

1623 r : array_like 

1624 2-D array of data with shape ``(u.size, v.size)``. 

1625 s : float, optional 

1626 Positive smoothing factor defined for estimation condition 

1627 (``s=0`` is for interpolation). 

1628 pole_continuity : bool or (bool, bool), optional 

1629 Order of continuity at the poles ``u=0`` (``pole_continuity[0]``) and 

1630 ``u=pi`` (``pole_continuity[1]``). The order of continuity at the pole 

1631 will be 1 or 0 when this is True or False, respectively. 

1632 Defaults to False. 

1633 pole_values : float or (float, float), optional 

1634 Data values at the poles ``u=0`` and ``u=pi``. Either the whole 

1635 parameter or each individual element can be None. Defaults to None. 

1636 pole_exact : bool or (bool, bool), optional 

1637 Data value exactness at the poles ``u=0`` and ``u=pi``. If True, the 

1638 value is considered to be the right function value, and it will be 

1639 fitted exactly. If False, the value will be considered to be a data 

1640 value just like the other data values. Defaults to False. 

1641 pole_flat : bool or (bool, bool), optional 

1642 For the poles at ``u=0`` and ``u=pi``, specify whether or not the 

1643 approximation has vanishing derivatives. Defaults to False. 

1644 

1645 See Also 

1646 -------- 

1647 RectBivariateSpline : bivariate spline approximation over a rectangular 

1648 mesh 

1649 

1650 Notes 

1651 ----- 

1652 Currently, only the smoothing spline approximation (``iopt[0] = 0`` and 

1653 ``iopt[0] = 1`` in the FITPACK routine) is supported. The exact 

1654 least-squares spline approximation is not implemented yet. 

1655 

1656 When actually performing the interpolation, the requested `v` values must 

1657 lie within the same length 2pi interval that the original `v` values were 

1658 chosen from. 

1659 

1660 For more information, see the FITPACK_ site about this function. 

1661 

1662 .. _FITPACK: http://www.netlib.org/dierckx/spgrid.f 

1663 

1664 Examples 

1665 -------- 

1666 Suppose we have global data on a coarse grid 

1667 

1668 >>> lats = np.linspace(10, 170, 9) * np.pi / 180. 

1669 >>> lons = np.linspace(0, 350, 18) * np.pi / 180. 

1670 >>> data = np.dot(np.atleast_2d(90. - np.linspace(-80., 80., 18)).T, 

1671 ... np.atleast_2d(180. - np.abs(np.linspace(0., 350., 9)))).T 

1672 

1673 We want to interpolate it to a global one-degree grid 

1674 

1675 >>> new_lats = np.linspace(1, 180, 180) * np.pi / 180 

1676 >>> new_lons = np.linspace(1, 360, 360) * np.pi / 180 

1677 >>> new_lats, new_lons = np.meshgrid(new_lats, new_lons) 

1678 

1679 We need to set up the interpolator object 

1680 

1681 >>> from scipy.interpolate import RectSphereBivariateSpline 

1682 >>> lut = RectSphereBivariateSpline(lats, lons, data) 

1683 

1684 Finally we interpolate the data. The `RectSphereBivariateSpline` object 

1685 only takes 1-D arrays as input, therefore we need to do some reshaping. 

1686 

1687 >>> data_interp = lut.ev(new_lats.ravel(), 

1688 ... new_lons.ravel()).reshape((360, 180)).T 

1689 

1690 Looking at the original and the interpolated data, one can see that the 

1691 interpolant reproduces the original data very well: 

1692 

1693 >>> import matplotlib.pyplot as plt 

1694 >>> fig = plt.figure() 

1695 >>> ax1 = fig.add_subplot(211) 

1696 >>> ax1.imshow(data, interpolation='nearest') 

1697 >>> ax2 = fig.add_subplot(212) 

1698 >>> ax2.imshow(data_interp, interpolation='nearest') 

1699 >>> plt.show() 

1700 

1701 Choosing the optimal value of ``s`` can be a delicate task. Recommended 

1702 values for ``s`` depend on the accuracy of the data values. If the user 

1703 has an idea of the statistical errors on the data, she can also find a 

1704 proper estimate for ``s``. By assuming that, if she specifies the 

1705 right ``s``, the interpolator will use a spline ``f(u,v)`` which exactly 

1706 reproduces the function underlying the data, she can evaluate 

1707 ``sum((r(i,j)-s(u(i),v(j)))**2)`` to find a good estimate for this ``s``. 

1708 For example, if she knows that the statistical errors on her 

1709 ``r(i,j)``-values are not greater than 0.1, she may expect that a good 

1710 ``s`` should have a value not larger than ``u.size * v.size * (0.1)**2``. 

1711 

1712 If nothing is known about the statistical error in ``r(i,j)``, ``s`` must 

1713 be determined by trial and error. The best is then to start with a very 

1714 large value of ``s`` (to determine the least-squares polynomial and the 

1715 corresponding upper bound ``fp0`` for ``s``) and then to progressively 

1716 decrease the value of ``s`` (say by a factor 10 in the beginning, i.e. 

1717 ``s = fp0 / 10, fp0 / 100, ...`` and more carefully as the approximation 

1718 shows more detail) to obtain closer fits. 

1719 

1720 The interpolation results for different values of ``s`` give some insight 

1721 into this process: 

1722 

1723 >>> fig2 = plt.figure() 

1724 >>> s = [3e9, 2e9, 1e9, 1e8] 

1725 >>> for ii in range(len(s)): 

1726 ... lut = RectSphereBivariateSpline(lats, lons, data, s=s[ii]) 

1727 ... data_interp = lut.ev(new_lats.ravel(), 

1728 ... new_lons.ravel()).reshape((360, 180)).T 

1729 ... ax = fig2.add_subplot(2, 2, ii+1) 

1730 ... ax.imshow(data_interp, interpolation='nearest') 

1731 ... ax.set_title("s = %g" % s[ii]) 

1732 >>> plt.show() 

1733 

1734 """ 

1735 

1736 def __init__(self, u, v, r, s=0., pole_continuity=False, pole_values=None, 

1737 pole_exact=False, pole_flat=False): 

1738 iopt = np.array([0, 0, 0], dtype=dfitpack_int) 

1739 ider = np.array([-1, 0, -1, 0], dtype=dfitpack_int) 

1740 if pole_values is None: 

1741 pole_values = (None, None) 

1742 elif isinstance(pole_values, (float, np.float32, np.float64)): 

1743 pole_values = (pole_values, pole_values) 

1744 if isinstance(pole_continuity, bool): 

1745 pole_continuity = (pole_continuity, pole_continuity) 

1746 if isinstance(pole_exact, bool): 

1747 pole_exact = (pole_exact, pole_exact) 

1748 if isinstance(pole_flat, bool): 

1749 pole_flat = (pole_flat, pole_flat) 

1750 

1751 r0, r1 = pole_values 

1752 iopt[1:] = pole_continuity 

1753 if r0 is None: 

1754 ider[0] = -1 

1755 else: 

1756 ider[0] = pole_exact[0] 

1757 

1758 if r1 is None: 

1759 ider[2] = -1 

1760 else: 

1761 ider[2] = pole_exact[1] 

1762 

1763 ider[1], ider[3] = pole_flat 

1764 

1765 u, v = np.ravel(u), np.ravel(v) 

1766 

1767 if not ((0.0 <= u).all() and (u <= np.pi).all()): 

1768 raise ValueError('u should be between [0, pi]') 

1769 if not -np.pi <= v[0] < np.pi: 

1770 raise ValueError('v[0] should be between [-pi, pi)') 

1771 if not v[-1] <= v[0] + 2*np.pi: 

1772 raise ValueError('v[-1] should be v[0] + 2pi or less ') 

1773 

1774 if not np.all(np.diff(u) > 0.0): 

1775 raise ValueError('u must be strictly increasing') 

1776 if not np.all(np.diff(v) > 0.0): 

1777 raise ValueError('v must be strictly increasing') 

1778 

1779 if not u.size == r.shape[0]: 

1780 raise ValueError('u dimension of r must have same number of ' 

1781 'elements as u') 

1782 if not v.size == r.shape[1]: 

1783 raise ValueError('v dimension of r must have same number of ' 

1784 'elements as v') 

1785 

1786 if pole_continuity[1] is False and pole_flat[1] is True: 

1787 raise ValueError('if pole_continuity is False, so must be ' 

1788 'pole_flat') 

1789 if pole_continuity[0] is False and pole_flat[0] is True: 

1790 raise ValueError('if pole_continuity is False, so must be ' 

1791 'pole_flat') 

1792 

1793 if not s >= 0.0: 

1794 raise ValueError('s should be positive') 

1795 

1796 r = np.ravel(r) 

1797 nu, tu, nv, tv, c, fp, ier = dfitpack.regrid_smth_spher(iopt, ider, 

1798 u.copy(), 

1799 v.copy(), 

1800 r.copy(), 

1801 r0, r1, s) 

1802 

1803 if ier not in [0, -1, -2]: 

1804 msg = _spfit_messages.get(ier, 'ier=%s' % (ier)) 

1805 raise ValueError(msg) 

1806 

1807 self.fp = fp 

1808 self.tck = tu[:nu], tv[:nv], c[:(nu - 4) * (nv-4)] 

1809 self.degrees = (3, 3)