Coverage for /usr/lib/python3/dist-packages/scipy/interpolate/_fitpack_impl.py: 6%

1""" 

2fitpack (dierckx in netlib) --- A Python-C wrapper to FITPACK (by P. Dierckx). 

3 FITPACK is a collection of FORTRAN programs for curve and surface 

4 fitting with splines and tensor product splines. 

5 

6See 

7 https://web.archive.org/web/20010524124604/http://www.cs.kuleuven.ac.be:80/cwis/research/nalag/research/topics/fitpack.html 

8or 

9 http://www.netlib.org/dierckx/ 

10 

11Copyright 2002 Pearu Peterson all rights reserved, 

12Pearu Peterson <pearu@cens.ioc.ee> 

13Permission to use, modify, and distribute this software is given under the 

14terms of the SciPy (BSD style) license. See LICENSE.txt that came with 

15this distribution for specifics. 

16 

17NO WARRANTY IS EXPRESSED OR IMPLIED. USE AT YOUR OWN RISK. 

18 

19TODO: Make interfaces to the following fitpack functions: 

20 For univariate splines: cocosp, concon, fourco, insert 

21 For bivariate splines: profil, regrid, parsur, surev 

22""" 

23 

24__all__ = ['splrep', 'splprep', 'splev', 'splint', 'sproot', 'spalde', 

25 'bisplrep', 'bisplev', 'insert', 'splder', 'splantider'] 

26 

27import warnings 

28import numpy as np 

29from . import _fitpack 

30from numpy import (atleast_1d, array, ones, zeros, sqrt, ravel, transpose, 

31 empty, iinfo, asarray) 

32 

33# Try to replace _fitpack interface with 

34# f2py-generated version 

35from . import dfitpack 

36 

37 

38dfitpack_int = dfitpack.types.intvar.dtype 

39 

40 

41def _int_overflow(x, msg=None): 

42 """Cast the value to an dfitpack_int and raise an OverflowError if the value 

43 cannot fit. 

44 """ 

45 if x > iinfo(dfitpack_int).max: 

46 if msg is None: 

47 msg = f'{x!r} cannot fit into an {dfitpack_int!r}' 

48 raise OverflowError(msg) 

49 return dfitpack_int.type(x) 

50 

51 

52_iermess = { 

53 0: ["The spline has a residual sum of squares fp such that " 

54 "abs(fp-s)/s<=0.001", None], 

55 -1: ["The spline is an interpolating spline (fp=0)", None], 

56 -2: ["The spline is weighted least-squares polynomial of degree k.\n" 

57 "fp gives the upper bound fp0 for the smoothing factor s", None], 

58 1: ["The required storage space exceeds the available storage space.\n" 

59 "Probable causes: data (x,y) size is too small or smoothing parameter" 

60 "\ns is too small (fp>s).", ValueError], 

61 2: ["A theoretically impossible result when finding a smoothing spline\n" 

62 "with fp = s. Probable cause: s too small. (abs(fp-s)/s>0.001)", 

63 ValueError], 

64 3: ["The maximal number of iterations (20) allowed for finding smoothing\n" 

65 "spline with fp=s has been reached. Probable cause: s too small.\n" 

66 "(abs(fp-s)/s>0.001)", ValueError], 

67 10: ["Error on input data", ValueError], 

68 'unknown': ["An error occurred", TypeError] 

69} 

70 

71_iermess2 = { 

72 0: ["The spline has a residual sum of squares fp such that " 

73 "abs(fp-s)/s<=0.001", None], 

74 -1: ["The spline is an interpolating spline (fp=0)", None], 

75 -2: ["The spline is weighted least-squares polynomial of degree kx and ky." 

76 "\nfp gives the upper bound fp0 for the smoothing factor s", None], 

77 -3: ["Warning. The coefficients of the spline have been computed as the\n" 

78 "minimal norm least-squares solution of a rank deficient system.", 

79 None], 

80 1: ["The required storage space exceeds the available storage space.\n" 

81 "Probable causes: nxest or nyest too small or s is too small. (fp>s)", 

82 ValueError], 

83 2: ["A theoretically impossible result when finding a smoothing spline\n" 

84 "with fp = s. Probable causes: s too small or badly chosen eps.\n" 

85 "(abs(fp-s)/s>0.001)", ValueError], 

86 3: ["The maximal number of iterations (20) allowed for finding smoothing\n" 

87 "spline with fp=s has been reached. Probable cause: s too small.\n" 

88 "(abs(fp-s)/s>0.001)", ValueError], 

89 4: ["No more knots can be added because the number of B-spline\n" 

90 "coefficients already exceeds the number of data points m.\n" 

91 "Probable causes: either s or m too small. (fp>s)", ValueError], 

92 5: ["No more knots can be added because the additional knot would\n" 

93 "coincide with an old one. Probable cause: s too small or too large\n" 

94 "a weight to an inaccurate data point. (fp>s)", ValueError], 

95 10: ["Error on input data", ValueError], 

96 11: ["rwrk2 too small, i.e., there is not enough workspace for computing\n" 

97 "the minimal least-squares solution of a rank deficient system of\n" 

98 "linear equations.", ValueError], 

99 'unknown': ["An error occurred", TypeError] 

100} 

101 

102_parcur_cache = {'t': array([], float), 'wrk': array([], float), 

103 'iwrk': array([], dfitpack_int), 'u': array([], float), 

104 'ub': 0, 'ue': 1} 

105 

106 

107def splprep(x, w=None, u=None, ub=None, ue=None, k=3, task=0, s=None, t=None, 

108 full_output=0, nest=None, per=0, quiet=1): 

109 # see the docstring of `_fitpack_py/splprep` 

110 if task <= 0: 

111 _parcur_cache = {'t': array([], float), 'wrk': array([], float), 

112 'iwrk': array([], dfitpack_int), 'u': array([], float), 

113 'ub': 0, 'ue': 1} 

114 x = atleast_1d(x) 

115 idim, m = x.shape 

116 if per: 

117 for i in range(idim): 

118 if x[i][0] != x[i][-1]: 

119 if not quiet: 

120 warnings.warn(RuntimeWarning('Setting x[%d][%d]=x[%d][0]' % 

121 (i, m, i))) 

122 x[i][-1] = x[i][0] 

123 if not 0 < idim < 11: 

124 raise TypeError('0 < idim < 11 must hold') 

125 if w is None: 

126 w = ones(m, float) 

127 else: 

128 w = atleast_1d(w) 

129 ipar = (u is not None) 

130 if ipar: 

131 _parcur_cache['u'] = u 

132 if ub is None: 

133 _parcur_cache['ub'] = u[0] 

134 else: 

135 _parcur_cache['ub'] = ub 

136 if ue is None: 

137 _parcur_cache['ue'] = u[-1] 

138 else: 

139 _parcur_cache['ue'] = ue 

140 else: 

141 _parcur_cache['u'] = zeros(m, float) 

142 if not (1 <= k <= 5): 

143 raise TypeError('1 <= k= %d <=5 must hold' % k) 

144 if not (-1 <= task <= 1): 

145 raise TypeError('task must be -1, 0 or 1') 

146 if (not len(w) == m) or (ipar == 1 and (not len(u) == m)): 

147 raise TypeError('Mismatch of input dimensions') 

148 if s is None: 

149 s = m - sqrt(2*m) 

150 if t is None and task == -1: 

151 raise TypeError('Knots must be given for task=-1') 

152 if t is not None: 

153 _parcur_cache['t'] = atleast_1d(t) 

154 n = len(_parcur_cache['t']) 

155 if task == -1 and n < 2*k + 2: 

156 raise TypeError('There must be at least 2*k+2 knots for task=-1') 

157 if m <= k: 

158 raise TypeError('m > k must hold') 

159 if nest is None: 

160 nest = m + 2*k 

161 

162 if (task >= 0 and s == 0) or (nest < 0): 

163 if per: 

164 nest = m + 2*k 

165 else: 

166 nest = m + k + 1 

167 nest = max(nest, 2*k + 3) 

168 u = _parcur_cache['u'] 

169 ub = _parcur_cache['ub'] 

170 ue = _parcur_cache['ue'] 

171 t = _parcur_cache['t'] 

172 wrk = _parcur_cache['wrk'] 

173 iwrk = _parcur_cache['iwrk'] 

174 t, c, o = _fitpack._parcur(ravel(transpose(x)), w, u, ub, ue, k, 

175 task, ipar, s, t, nest, wrk, iwrk, per) 

176 _parcur_cache['u'] = o['u'] 

177 _parcur_cache['ub'] = o['ub'] 

178 _parcur_cache['ue'] = o['ue'] 

179 _parcur_cache['t'] = t 

180 _parcur_cache['wrk'] = o['wrk'] 

181 _parcur_cache['iwrk'] = o['iwrk'] 

182 ier = o['ier'] 

183 fp = o['fp'] 

184 n = len(t) 

185 u = o['u'] 

186 c.shape = idim, n - k - 1 

187 tcku = [t, list(c), k], u 

188 if ier <= 0 and not quiet: 

189 warnings.warn(RuntimeWarning(_iermess[ier][0] + 

190 "\tk=%d n=%d m=%d fp=%f s=%f" % 

191 (k, len(t), m, fp, s))) 

192 if ier > 0 and not full_output: 

193 if ier in [1, 2, 3]: 

194 warnings.warn(RuntimeWarning(_iermess[ier][0])) 

195 else: 

196 try: 

197 raise _iermess[ier][1](_iermess[ier][0]) 

198 except KeyError as e: 

199 raise _iermess['unknown'][1](_iermess['unknown'][0]) from e 

200 if full_output: 

201 try: 

202 return tcku, fp, ier, _iermess[ier][0] 

203 except KeyError: 

204 return tcku, fp, ier, _iermess['unknown'][0] 

205 else: 

206 return tcku 

207 

208 

209_curfit_cache = {'t': array([], float), 'wrk': array([], float), 

210 'iwrk': array([], dfitpack_int)} 

211 

212 

213def splrep(x, y, w=None, xb=None, xe=None, k=3, task=0, s=None, t=None, 

214 full_output=0, per=0, quiet=1): 

215 # see the docstring of `_fitpack_py/splrep` 

216 if task <= 0: 

217 _curfit_cache = {} 

218 x, y = map(atleast_1d, [x, y]) 

219 m = len(x) 

220 if w is None: 

221 w = ones(m, float) 

222 if s is None: 

223 s = 0.0 

224 else: 

225 w = atleast_1d(w) 

226 if s is None: 

227 s = m - sqrt(2*m) 

228 if not len(w) == m: 

229 raise TypeError('len(w)=%d is not equal to m=%d' % (len(w), m)) 

230 if (m != len(y)) or (m != len(w)): 

231 raise TypeError('Lengths of the first three arguments (x,y,w) must ' 

232 'be equal') 

233 if not (1 <= k <= 5): 

234 raise TypeError('Given degree of the spline (k=%d) is not supported. ' 

235 '(1<=k<=5)' % k) 

236 if m <= k: 

237 raise TypeError('m > k must hold') 

238 if xb is None: 

239 xb = x[0] 

240 if xe is None: 

241 xe = x[-1] 

242 if not (-1 <= task <= 1): 

243 raise TypeError('task must be -1, 0 or 1') 

244 if t is not None: 

245 task = -1 

246 if task == -1: 

247 if t is None: 

248 raise TypeError('Knots must be given for task=-1') 

249 numknots = len(t) 

250 _curfit_cache['t'] = empty((numknots + 2*k + 2,), float) 

251 _curfit_cache['t'][k+1:-k-1] = t 

252 nest = len(_curfit_cache['t']) 

253 elif task == 0: 

254 if per: 

255 nest = max(m + 2*k, 2*k + 3) 

256 else: 

257 nest = max(m + k + 1, 2*k + 3) 

258 t = empty((nest,), float) 

259 _curfit_cache['t'] = t 

260 if task <= 0: 

261 if per: 

262 _curfit_cache['wrk'] = empty((m*(k + 1) + nest*(8 + 5*k),), float) 

263 else: 

264 _curfit_cache['wrk'] = empty((m*(k + 1) + nest*(7 + 3*k),), float) 

265 _curfit_cache['iwrk'] = empty((nest,), dfitpack_int) 

266 try: 

267 t = _curfit_cache['t'] 

268 wrk = _curfit_cache['wrk'] 

269 iwrk = _curfit_cache['iwrk'] 

270 except KeyError as e: 

271 raise TypeError("must call with task=1 only after" 

272 " call with task=0,-1") from e 

273 if not per: 

274 n, c, fp, ier = dfitpack.curfit(task, x, y, w, t, wrk, iwrk, 

275 xb, xe, k, s) 

276 else: 

277 n, c, fp, ier = dfitpack.percur(task, x, y, w, t, wrk, iwrk, k, s) 

278 tck = (t[:n], c[:n], k) 

279 if ier <= 0 and not quiet: 

280 _mess = (_iermess[ier][0] + "\tk=%d n=%d m=%d fp=%f s=%f" % 

281 (k, len(t), m, fp, s)) 

282 warnings.warn(RuntimeWarning(_mess)) 

283 if ier > 0 and not full_output: 

284 if ier in [1, 2, 3]: 

285 warnings.warn(RuntimeWarning(_iermess[ier][0])) 

286 else: 

287 try: 

288 raise _iermess[ier][1](_iermess[ier][0]) 

289 except KeyError as e: 

290 raise _iermess['unknown'][1](_iermess['unknown'][0]) from e 

291 if full_output: 

292 try: 

293 return tck, fp, ier, _iermess[ier][0] 

294 except KeyError: 

295 return tck, fp, ier, _iermess['unknown'][0] 

296 else: 

297 return tck 

298 

299 

300def splev(x, tck, der=0, ext=0): 

301 # see the docstring of `_fitpack_py/splev` 

302 t, c, k = tck 

303 try: 

304 c[0][0] 

305 parametric = True 

306 except Exception: 

307 parametric = False 

308 if parametric: 

309 return list(map(lambda c, x=x, t=t, k=k, der=der: 

310 splev(x, [t, c, k], der, ext), c)) 

311 else: 

312 if not (0 <= der <= k): 

313 raise ValueError("0<=der=%d<=k=%d must hold" % (der, k)) 

314 if ext not in (0, 1, 2, 3): 

315 raise ValueError("ext = %s not in (0, 1, 2, 3) " % ext) 

316 

317 x = asarray(x) 

318 shape = x.shape 

319 x = atleast_1d(x).ravel() 

320 if der == 0: 

321 y, ier = dfitpack.splev(t, c, k, x, ext) 

322 else: 

323 y, ier = dfitpack.splder(t, c, k, x, der, ext) 

324 

325 if ier == 10: 

326 raise ValueError("Invalid input data") 

327 if ier == 1: 

328 raise ValueError("Found x value not in the domain") 

329 if ier: 

330 raise TypeError("An error occurred") 

331 

332 return y.reshape(shape) 

333 

334 

335def splint(a, b, tck, full_output=0): 

336 # see the docstring of `_fitpack_py/splint` 

337 t, c, k = tck 

338 try: 

339 c[0][0] 

340 parametric = True 

341 except Exception: 

342 parametric = False 

343 if parametric: 

344 return list(map(lambda c, a=a, b=b, t=t, k=k: 

345 splint(a, b, [t, c, k]), c)) 

346 else: 

347 aint, wrk = dfitpack.splint(t, c, k, a, b) 

348 if full_output: 

349 return aint, wrk 

350 else: 

351 return aint 

352 

353 

354def sproot(tck, mest=10): 

355 # see the docstring of `_fitpack_py/sproot` 

356 t, c, k = tck 

357 if k != 3: 

358 raise ValueError("sproot works only for cubic (k=3) splines") 

359 try: 

360 c[0][0] 

361 parametric = True 

362 except Exception: 

363 parametric = False 

364 if parametric: 

365 return list(map(lambda c, t=t, k=k, mest=mest: 

366 sproot([t, c, k], mest), c)) 

367 else: 

368 if len(t) < 8: 

369 raise TypeError("The number of knots %d>=8" % len(t)) 

370 z, m, ier = dfitpack.sproot(t, c, mest) 

371 if ier == 10: 

372 raise TypeError("Invalid input data. " 

373 "t1<=..<=t4<t5<..<tn-3<=..<=tn must hold.") 

374 if ier == 0: 

375 return z[:m] 

376 if ier == 1: 

377 warnings.warn(RuntimeWarning("The number of zeros exceeds mest")) 

378 return z[:m] 

379 raise TypeError("Unknown error") 

380 

381 

382def spalde(x, tck): 

383 # see the docstring of `_fitpack_py/spalde` 

384 t, c, k = tck 

385 try: 

386 c[0][0] 

387 parametric = True 

388 except Exception: 

389 parametric = False 

390 if parametric: 

391 return list(map(lambda c, x=x, t=t, k=k: 

392 spalde(x, [t, c, k]), c)) 

393 else: 

394 x = atleast_1d(x) 

395 if len(x) > 1: 

396 return list(map(lambda x, tck=tck: spalde(x, tck), x)) 

397 d, ier = dfitpack.spalde(t, c, k+1, x[0]) 

398 if ier == 0: 

399 return d 

400 if ier == 10: 

401 raise TypeError("Invalid input data. t(k)<=x<=t(n-k+1) must hold.") 

402 raise TypeError("Unknown error") 

403 

404# def _curfit(x,y,w=None,xb=None,xe=None,k=3,task=0,s=None,t=None, 

405# full_output=0,nest=None,per=0,quiet=1): 

406 

407 

408_surfit_cache = {'tx': array([], float), 'ty': array([], float), 

409 'wrk': array([], float), 'iwrk': array([], dfitpack_int)} 

410 

411 

412def bisplrep(x, y, z, w=None, xb=None, xe=None, yb=None, ye=None, 

413 kx=3, ky=3, task=0, s=None, eps=1e-16, tx=None, ty=None, 

414 full_output=0, nxest=None, nyest=None, quiet=1): 

415 """ 

416 Find a bivariate B-spline representation of a surface. 

417 

418 Given a set of data points (x[i], y[i], z[i]) representing a surface 

419 z=f(x,y), compute a B-spline representation of the surface. Based on 

420 the routine SURFIT from FITPACK. 

421 

422 Parameters 

423 ---------- 

424 x, y, z : ndarray 

425 Rank-1 arrays of data points. 

426 w : ndarray, optional 

427 Rank-1 array of weights. By default ``w=np.ones(len(x))``. 

428 xb, xe : float, optional 

429 End points of approximation interval in `x`. 

430 By default ``xb = x.min(), xe=x.max()``. 

431 yb, ye : float, optional 

432 End points of approximation interval in `y`. 

433 By default ``yb=y.min(), ye = y.max()``. 

434 kx, ky : int, optional 

435 The degrees of the spline (1 <= kx, ky <= 5). 

436 Third order (kx=ky=3) is recommended. 

437 task : int, optional 

438 If task=0, find knots in x and y and coefficients for a given 

439 smoothing factor, s. 

440 If task=1, find knots and coefficients for another value of the 

441 smoothing factor, s. bisplrep must have been previously called 

442 with task=0 or task=1. 

443 If task=-1, find coefficients for a given set of knots tx, ty. 

444 s : float, optional 

445 A non-negative smoothing factor. If weights correspond 

446 to the inverse of the standard-deviation of the errors in z, 

447 then a good s-value should be found in the range 

448 ``(m-sqrt(2*m),m+sqrt(2*m))`` where m=len(x). 

449 eps : float, optional 

450 A threshold for determining the effective rank of an 

451 over-determined linear system of equations (0 < eps < 1). 

452 `eps` is not likely to need changing. 

453 tx, ty : ndarray, optional 

454 Rank-1 arrays of the knots of the spline for task=-1 

455 full_output : int, optional 

456 Non-zero to return optional outputs. 

457 nxest, nyest : int, optional 

458 Over-estimates of the total number of knots. If None then 

459 ``nxest = max(kx+sqrt(m/2),2*kx+3)``, 

460 ``nyest = max(ky+sqrt(m/2),2*ky+3)``. 

461 quiet : int, optional 

462 Non-zero to suppress printing of messages. 

463 

464 Returns 

465 ------- 

466 tck : array_like 

467 A list [tx, ty, c, kx, ky] containing the knots (tx, ty) and 

468 coefficients (c) of the bivariate B-spline representation of the 

469 surface along with the degree of the spline. 

470 fp : ndarray 

471 The weighted sum of squared residuals of the spline approximation. 

472 ier : int 

473 An integer flag about splrep success. Success is indicated if 

474 ier<=0. If ier in [1,2,3] an error occurred but was not raised. 

475 Otherwise an error is raised. 

476 msg : str 

477 A message corresponding to the integer flag, ier. 

478 

479 See Also 

480 -------- 

481 splprep, splrep, splint, sproot, splev 

482 UnivariateSpline, BivariateSpline 

483 

484 Notes 

485 ----- 

486 See `bisplev` to evaluate the value of the B-spline given its tck 

487 representation. 

488 

489 If the input data is such that input dimensions have incommensurate 

490 units and differ by many orders of magnitude, the interpolant may have 

491 numerical artifacts. Consider rescaling the data before interpolation. 

492 

493 References 

494 ---------- 

495 .. [1] Dierckx P.:An algorithm for surface fitting with spline functions 

496 Ima J. Numer. Anal. 1 (1981) 267-283. 

497 .. [2] Dierckx P.:An algorithm for surface fitting with spline functions 

498 report tw50, Dept. Computer Science,K.U.Leuven, 1980. 

499 .. [3] Dierckx P.:Curve and surface fitting with splines, Monographs on 

500 Numerical Analysis, Oxford University Press, 1993. 

501 

502 Examples 

503 -------- 

504 Examples are given :ref:`in the tutorial <tutorial-interpolate_2d_spline>`. 

505 

506 """ 

507 x, y, z = map(ravel, [x, y, z]) # ensure 1-d arrays. 

508 m = len(x) 

509 if not (m == len(y) == len(z)): 

510 raise TypeError('len(x)==len(y)==len(z) must hold.') 

511 if w is None: 

512 w = ones(m, float) 

513 else: 

514 w = atleast_1d(w) 

515 if not len(w) == m: 

516 raise TypeError('len(w)=%d is not equal to m=%d' % (len(w), m)) 

517 if xb is None: 

518 xb = x.min() 

519 if xe is None: 

520 xe = x.max() 

521 if yb is None: 

522 yb = y.min() 

523 if ye is None: 

524 ye = y.max() 

525 if not (-1 <= task <= 1): 

526 raise TypeError('task must be -1, 0 or 1') 

527 if s is None: 

528 s = m - sqrt(2*m) 

529 if tx is None and task == -1: 

530 raise TypeError('Knots_x must be given for task=-1') 

531 if tx is not None: 

532 _surfit_cache['tx'] = atleast_1d(tx) 

533 nx = len(_surfit_cache['tx']) 

534 if ty is None and task == -1: 

535 raise TypeError('Knots_y must be given for task=-1') 

536 if ty is not None: 

537 _surfit_cache['ty'] = atleast_1d(ty) 

538 ny = len(_surfit_cache['ty']) 

539 if task == -1 and nx < 2*kx+2: 

540 raise TypeError('There must be at least 2*kx+2 knots_x for task=-1') 

541 if task == -1 and ny < 2*ky+2: 

542 raise TypeError('There must be at least 2*ky+2 knots_x for task=-1') 

543 if not ((1 <= kx <= 5) and (1 <= ky <= 5)): 

544 raise TypeError('Given degree of the spline (kx,ky=%d,%d) is not ' 

545 'supported. (1<=k<=5)' % (kx, ky)) 

546 if m < (kx + 1)*(ky + 1): 

547 raise TypeError('m >= (kx+1)(ky+1) must hold') 

548 if nxest is None: 

549 nxest = int(kx + sqrt(m/2)) 

550 if nyest is None: 

551 nyest = int(ky + sqrt(m/2)) 

552 nxest, nyest = max(nxest, 2*kx + 3), max(nyest, 2*ky + 3) 

553 if task >= 0 and s == 0: 

554 nxest = int(kx + sqrt(3*m)) 

555 nyest = int(ky + sqrt(3*m)) 

556 if task == -1: 

557 _surfit_cache['tx'] = atleast_1d(tx) 

558 _surfit_cache['ty'] = atleast_1d(ty) 

559 tx, ty = _surfit_cache['tx'], _surfit_cache['ty'] 

560 wrk = _surfit_cache['wrk'] 

561 u = nxest - kx - 1 

562 v = nyest - ky - 1 

563 km = max(kx, ky) + 1 

564 ne = max(nxest, nyest) 

565 bx, by = kx*v + ky + 1, ky*u + kx + 1 

566 b1, b2 = bx, bx + v - ky 

567 if bx > by: 

568 b1, b2 = by, by + u - kx 

569 msg = "Too many data points to interpolate" 

570 lwrk1 = _int_overflow(u*v*(2 + b1 + b2) + 

571 2*(u + v + km*(m + ne) + ne - kx - ky) + b2 + 1, 

572 msg=msg) 

573 lwrk2 = _int_overflow(u*v*(b2 + 1) + b2, msg=msg) 

574 tx, ty, c, o = _fitpack._surfit(x, y, z, w, xb, xe, yb, ye, kx, ky, 

575 task, s, eps, tx, ty, nxest, nyest, 

576 wrk, lwrk1, lwrk2) 

577 _curfit_cache['tx'] = tx 

578 _curfit_cache['ty'] = ty 

579 _curfit_cache['wrk'] = o['wrk'] 

580 ier, fp = o['ier'], o['fp'] 

581 tck = [tx, ty, c, kx, ky] 

582 

583 ierm = min(11, max(-3, ier)) 

584 if ierm <= 0 and not quiet: 

585 _mess = (_iermess2[ierm][0] + 

586 "\tkx,ky=%d,%d nx,ny=%d,%d m=%d fp=%f s=%f" % 

587 (kx, ky, len(tx), len(ty), m, fp, s)) 

588 warnings.warn(RuntimeWarning(_mess)) 

589 if ierm > 0 and not full_output: 

590 if ier in [1, 2, 3, 4, 5]: 

591 _mess = ("\n\tkx,ky=%d,%d nx,ny=%d,%d m=%d fp=%f s=%f" % 

592 (kx, ky, len(tx), len(ty), m, fp, s)) 

593 warnings.warn(RuntimeWarning(_iermess2[ierm][0] + _mess)) 

594 else: 

595 try: 

596 raise _iermess2[ierm][1](_iermess2[ierm][0]) 

597 except KeyError as e: 

598 raise _iermess2['unknown'][1](_iermess2['unknown'][0]) from e 

599 if full_output: 

600 try: 

601 return tck, fp, ier, _iermess2[ierm][0] 

602 except KeyError: 

603 return tck, fp, ier, _iermess2['unknown'][0] 

604 else: 

605 return tck 

606 

607 

608def bisplev(x, y, tck, dx=0, dy=0): 

609 """ 

610 Evaluate a bivariate B-spline and its derivatives. 

611 

612 Return a rank-2 array of spline function values (or spline derivative 

613 values) at points given by the cross-product of the rank-1 arrays `x` and 

614 `y`. In special cases, return an array or just a float if either `x` or 

615 `y` or both are floats. Based on BISPEV and PARDER from FITPACK. 

616 

617 Parameters 

618 ---------- 

619 x, y : ndarray 

620 Rank-1 arrays specifying the domain over which to evaluate the 

621 spline or its derivative. 

622 tck : tuple 

623 A sequence of length 5 returned by `bisplrep` containing the knot 

624 locations, the coefficients, and the degree of the spline: 

625 [tx, ty, c, kx, ky]. 

626 dx, dy : int, optional 

627 The orders of the partial derivatives in `x` and `y` respectively. 

628 

629 Returns 

630 ------- 

631 vals : ndarray 

632 The B-spline or its derivative evaluated over the set formed by 

633 the cross-product of `x` and `y`. 

634 

635 See Also 

636 -------- 

637 splprep, splrep, splint, sproot, splev 

638 UnivariateSpline, BivariateSpline 

639 

640 Notes 

641 ----- 

642 See `bisplrep` to generate the `tck` representation. 

643 

644 References 

645 ---------- 

646 .. [1] Dierckx P. : An algorithm for surface fitting 

647 with spline functions 

648 Ima J. Numer. Anal. 1 (1981) 267-283. 

649 .. [2] Dierckx P. : An algorithm for surface fitting 

650 with spline functions 

651 report tw50, Dept. Computer Science,K.U.Leuven, 1980. 

652 .. [3] Dierckx P. : Curve and surface fitting with splines, 

653 Monographs on Numerical Analysis, Oxford University Press, 1993. 

654 

655 Examples 

656 -------- 

657 Examples are given :ref:`in the tutorial <tutorial-interpolate_2d_spline>`. 

658 

659 """ 

660 tx, ty, c, kx, ky = tck 

661 if not (0 <= dx < kx): 

662 raise ValueError("0 <= dx = %d < kx = %d must hold" % (dx, kx)) 

663 if not (0 <= dy < ky): 

664 raise ValueError("0 <= dy = %d < ky = %d must hold" % (dy, ky)) 

665 x, y = map(atleast_1d, [x, y]) 

666 if (len(x.shape) != 1) or (len(y.shape) != 1): 

667 raise ValueError("First two entries should be rank-1 arrays.") 

668 z, ier = _fitpack._bispev(tx, ty, c, kx, ky, x, y, dx, dy) 

669 if ier == 10: 

670 raise ValueError("Invalid input data") 

671 if ier: 

672 raise TypeError("An error occurred") 

673 z.shape = len(x), len(y) 

674 if len(z) > 1: 

675 return z 

676 if len(z[0]) > 1: 

677 return z[0] 

678 return z[0][0] 

679 

680 

681def dblint(xa, xb, ya, yb, tck): 

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

683 

684 Parameters 

685 ---------- 

686 xa, xb : float 

687 The end-points of the x integration interval. 

688 ya, yb : float 

689 The end-points of the y integration interval. 

690 tck : list [tx, ty, c, kx, ky] 

691 A sequence of length 5 returned by bisplrep containing the knot 

692 locations tx, ty, the coefficients c, and the degrees kx, ky 

693 of the spline. 

694 

695 Returns 

696 ------- 

697 integ : float 

698 The value of the resulting integral. 

699 """ 

700 tx, ty, c, kx, ky = tck 

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

702 

703 

704def insert(x, tck, m=1, per=0): 

705 # see the docstring of `_fitpack_py/insert` 

706 t, c, k = tck 

707 try: 

708 c[0][0] 

709 parametric = True 

710 except Exception: 

711 parametric = False 

712 if parametric: 

713 cc = [] 

714 for c_vals in c: 

715 tt, cc_val, kk = insert(x, [t, c_vals, k], m) 

716 cc.append(cc_val) 

717 return (tt, cc, kk) 

718 else: 

719 tt, cc, ier = _fitpack._insert(per, t, c, k, x, m) 

720 if ier == 10: 

721 raise ValueError("Invalid input data") 

722 if ier: 

723 raise TypeError("An error occurred") 

724 return (tt, cc, k) 

725 

726 

727def splder(tck, n=1): 

728 # see the docstring of `_fitpack_py/splder` 

729 if n < 0: 

730 return splantider(tck, -n) 

731 

732 t, c, k = tck 

733 

734 if n > k: 

735 raise ValueError(("Order of derivative (n = {!r}) must be <= " 

736 "order of spline (k = {!r})").format(n, tck[2]))