Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/interpolate/fitpack2.py : 15%

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
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']
22import warnings
24from numpy import zeros, concatenate, ravel, diff, array, ones
25import numpy as np
27from . import fitpack
28from . import dfitpack
31dfitpack_int = dfitpack.types.intvar.dtype
34# ############### Univariate spline ####################
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 }
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}
71class UnivariateSpline(object):
72 """
73 1-D smoothing spline fit to a given set of data points.
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.
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::
98 sum((w[i] * (y[i]-spl(x[i])))**2, axis=0) <= s
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.
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.
112 The default value is 0.
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.
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
130 Notes
131 -----
132 The number of data points must be larger than the spline degree `k`.
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:
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)
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.)
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)
156 Use the default value for the smoothing parameter:
158 >>> spl = UnivariateSpline(x, y)
159 >>> xs = np.linspace(-3, 3, 1000)
160 >>> plt.plot(xs, spl(xs), 'g', lw=3)
162 Manually change the amount of smoothing:
164 >>> spl.set_smoothing_factor(0.5)
165 >>> plt.plot(xs, spl(xs), 'b', lw=3)
166 >>> plt.show()
168 """
169 def __init__(self, x, y, w=None, bbox=[None]*2, k=3, s=None,
170 ext=0, check_finite=False):
172 x, y, w, bbox, self.ext = self.validate_input(x, y, w, bbox, k, s, ext,
173 check_finite)
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()
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")
212 try:
213 ext = _extrap_modes[ext]
214 except KeyError:
215 raise ValueError("Unknown extrapolation mode %s." % ext)
217 return x, y, w, bbox, ext
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
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)
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
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]]
275 args = data[:8] + (t, c, n, fpint, nrdata, data[13])
276 data = dfitpack.fpcurf1(*args)
277 return data
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.
283 This routine modifies the spline in place.
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()
299 def __call__(self, x, nu=0, ext=None):
300 """
301 Evaluate spline (or its nu-th derivative) at positions x.
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.
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.
320 The default value is 0, passed from the initialization of
321 UnivariateSpline.
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)
337 def get_knots(self):
338 """ Return positions of interior knots of the spline.
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]
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]
352 def get_residual(self):
353 """Return weighted sum of squared residuals of the spline approximation.
355 This is equivalent to::
357 sum((w[i] * (y[i]-spl(x[i])))**2, axis=0)
359 """
360 return self._data[10]
362 def integral(self, a, b):
363 """ Return definite integral of the spline between two given points.
365 Parameters
366 ----------
367 a : float
368 Lower limit of integration.
369 b : float
370 Upper limit of integration.
372 Returns
373 -------
374 integral : float
375 The value of the definite integral of the spline between limits.
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
386 which agrees with :math:`\\int x^2 dx = x^3 / 3` between the limits
387 of 0 and 3.
389 A caveat is that this routine assumes the spline to be zero outside of
390 the data limits:
392 >>> spl.integral(-1, 4)
393 9.0
394 >>> spl.integral(-1, 0)
395 0.0
397 """
398 return dfitpack.splint(*(self._eval_args+(a, b)))
400 def derivatives(self, x):
401 """ Return all derivatives of the spline at the point x.
403 Parameters
404 ----------
405 x : float
406 The point to evaluate the derivatives at.
408 Returns
409 -------
410 der : ndarray, shape(k+1,)
411 Derivatives of the orders 0 to k.
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])
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
428 def roots(self):
429 """ Return the zeros of the spline.
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')
442 def derivative(self, n=1):
443 """
444 Construct a new spline representing the derivative of this spline.
446 Parameters
447 ----------
448 n : int, optional
449 Order of derivative to evaluate. Default: 1
451 Returns
452 -------
453 spline : UnivariateSpline
454 Spline of order k2=k-n representing the derivative of this
455 spline.
457 See Also
458 --------
459 splder, antiderivative
461 Notes
462 -----
464 .. versionadded:: 0.13.0
466 Examples
467 --------
468 This can be used for finding maxima of a curve:
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)
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):
479 >>> spl.derivative().roots() / np.pi
480 array([ 0.50000001, 1.5 , 2.49999998])
482 This agrees well with roots :math:`\\pi/2 + n\\pi` of
483 :math:`\\cos(x) = \\sin'(x)`.
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)
491 def antiderivative(self, n=1):
492 """
493 Construct a new spline representing the antiderivative of this spline.
495 Parameters
496 ----------
497 n : int, optional
498 Order of antiderivative to evaluate. Default: 1
500 Returns
501 -------
502 spline : UnivariateSpline
503 Spline of order k2=k+n representing the antiderivative of this
504 spline.
506 Notes
507 -----
509 .. versionadded:: 0.13.0
511 See Also
512 --------
513 splantider, derivative
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)
522 The derivative is the inverse operation of the antiderivative,
523 although some floating point error accumulates:
525 >>> spl(1.7), spl.antiderivative().derivative()(1.7)
526 (array(2.1565429877197317), array(2.1565429877201865))
528 Antiderivative can be used to evaluate definite integrals:
530 >>> ispl = spl.antiderivative()
531 >>> ispl(np.pi/2) - ispl(0)
532 2.2572053588768486
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`:
537 >>> from scipy.special import ellipk
538 >>> ellipk(0.8)
539 2.2572053268208538
541 """
542 tck = fitpack.splantider(self._eval_args, n)
543 return UnivariateSpline._from_tck(tck, self.ext)
546class InterpolatedUnivariateSpline(UnivariateSpline):
547 """
548 1-D interpolating spline for a given set of data points.
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.
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.
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.
577 The default value is 0.
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.
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
595 Notes
596 -----
597 The number of data points must be larger than the spline degree `k`.
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()
611 Notice that the ``spl(x)`` interpolates `y`:
613 >>> spl.get_residual()
614 0.0
616 """
617 def __init__(self, x, y, w=None, bbox=[None]*2, k=3,
618 ext=0, check_finite=False):
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')
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()
631_fpchec_error_string = """The input parameters have been rejected by fpchec. \
632This means that at least one of the following conditions is violated:
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"""
646class LSQUnivariateSpline(UnivariateSpline):
647 """
648 1-D spline with explicit internal knots.
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
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::
662 bbox[0] < t[0] < ... < t[-1] < bbox[-1]
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.
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.
682 The default value is 0.
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.
691 Raises
692 ------
693 ValueError
694 If the interior knots do not satisfy the Schoenberg-Whitney conditions
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
705 Notes
706 -----
707 The number of data points must be larger than the spline degree `k`.
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``.
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)
720 Fit a smoothing spline with a pre-defined internal knots:
722 >>> t = [-1, 0, 1]
723 >>> spl = LSQUnivariateSpline(x, y, t)
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()
730 Check the knot vector:
732 >>> spl.get_knots()
733 array([-3., -1., 0., 1., 3.])
735 Constructing lsq spline using the knots from another spline:
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.])
746 """
748 def __init__(self, x, y, t, w=None, bbox=[None]*2, k=3,
749 ext=0, check_finite=False):
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')
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()
775# ############### Bivariate spline ####################
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).
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 """
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
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]
805 def get_coeffs(self):
806 """ Return spline coefficients."""
807 return self.tck[2]
809 def __call__(self, x, y, dx=0, dy=0, grid=True):
810 """
811 Evaluate the spline or its derivatives at given positions.
813 Parameters
814 ----------
815 x, y : array_like
816 Input coordinates.
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.
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.
826 Note that the axis ordering is inverted relative to
827 the output of meshgrid.
828 dx : int
829 Order of x-derivative
831 .. versionadded:: 0.14.0
832 dy : int
833 Order of y-derivative
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.
840 .. versionadded:: 0.14.0
842 """
843 x = np.asarray(x)
844 y = np.asarray(y)
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)
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)
865 shape = x.shape
866 x = x.ravel()
867 y = y.ravel()
869 if x.size == 0 or y.size == 0:
870 return np.zeros(shape, dtype=self.tck[2].dtype)
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)
881 z = z.reshape(shape)
882 return z
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 }
927class BivariateSpline(_BivariateSplineBase):
928 """
929 Base class for bivariate splines.
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)``.
935 This class is meant to be subclassed, not instantiated directly.
936 To construct these splines, call either `SmoothBivariateSpline` or
937 `LSQBivariateSpline`.
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 """
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
969 def ev(self, xi, yi, dx=0, dy=0):
970 """
971 Evaluate the spline at points
973 Returns the interpolated value at ``(xi[i], yi[i]),
974 i=0,...,len(xi)-1``.
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
983 .. versionadded:: 0.14.0
984 dy : int, optional
985 Order of y-derivative
987 .. versionadded:: 0.14.0
988 """
989 return self.__call__(xi, yi, dx=dx, dy=dy, grid=False)
991 def integral(self, xa, xb, ya, yb):
992 """
993 Evaluate the integral of the spline over area [xa,xb] x [ya,yb].
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.
1002 Returns
1003 -------
1004 integ : float
1005 The value of the resulting integral.
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)
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')
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
1032class SmoothBivariateSpline(BivariateSpline):
1033 """
1034 Smooth bivariate spline approximation.
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.
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
1065 Notes
1066 -----
1067 The length of `x`, `y` and `z` should be at least ``(kx+1) * (ky+1)``.
1069 """
1071 def __init__(self, x, y, z, w=None, bbox=[None] * 4, kx=3, ky=3, s=None,
1072 eps=1e-16):
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")
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)
1100 self.fp = fp
1101 self.tck = tx[:nx], ty[:ny], c[:(nx-kx-1)*(ny-ky-1)]
1102 self.degrees = kx, ky
1105class LSQBivariateSpline(BivariateSpline):
1106 """
1107 Weighted least-squares bivariate spline approximation.
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.
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
1135 Notes
1136 -----
1137 The length of `x`, `y` and `z` should be at least ``(kx+1) * (ky+1)``.
1139 """
1141 def __init__(self, x, y, z, tx, ty, w=None, bbox=[None]*4, kx=3, ky=3,
1142 eps=None):
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,)')
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
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
1178class RectBivariateSpline(BivariateSpline):
1179 """
1180 Bivariate spline approximation over a rectangular mesh.
1182 Can be used for both smoothing and interpolating data.
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.
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
1208 """
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")
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)
1233 if ier not in [0, -1, -2]:
1234 msg = _surfit_messages.get(ier, 'ier=%s' % (ier))
1235 raise ValueError(msg)
1237 self.fp = fp
1238 self.tck = tx[:nx], ty[:ny], c[:(nx - kx - 1) * (ny - ky - 1)]
1239 self.degrees = kx, ky
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."""
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).
1271 .. versionadded:: 0.11.0
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 """
1283 def __call__(self, theta, phi, dtheta=0, dphi=0, grid=True):
1284 """
1285 Evaluate the spline or its derivatives at given positions.
1287 Parameters
1288 ----------
1289 theta, phi : array_like
1290 Input coordinates.
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.
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
1302 .. versionadded:: 0.14.0
1303 dphi : int
1304 Order of phi-derivative
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.
1311 .. versionadded:: 0.14.0
1313 """
1314 theta = np.asarray(theta)
1315 phi = np.asarray(phi)
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.")
1322 return _BivariateSplineBase.__call__(self, theta, phi,
1323 dx=dtheta, dy=dphi, grid=grid)
1325 def ev(self, theta, phi, dtheta=0, dphi=0):
1326 """
1327 Evaluate the spline at points
1329 Returns the interpolated value at ``(theta[i], phi[i]),
1330 i=0,...,len(theta)-1``.
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
1339 .. versionadded:: 0.14.0
1340 dphi : int, optional
1341 Order of phi-derivative
1343 .. versionadded:: 0.14.0
1344 """
1345 return self.__call__(theta, phi, dtheta=dtheta, dphi=dphi, grid=False)
1348class SmoothSphereBivariateSpline(SphereBivariateSpline):
1349 """
1350 Smooth bivariate spline approximation in spherical coordinates.
1352 .. versionadded:: 0.11.0
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.
1372 Notes
1373 -----
1374 For more information, see the FITPACK_ site about this function.
1376 .. _FITPACK: http://www.netlib.org/dierckx/sphere.f
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):
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)
1394 We need to set up the interpolator object
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)
1401 As a first test, we'll see what the algorithm returns when run on the
1402 input coordinates
1404 >>> data_orig = lut(theta, phi)
1406 Finally we interpolate the data to a finer grid
1408 >>> fine_lats = np.linspace(0., np.pi, 70)
1409 >>> fine_lons = np.linspace(0., 2 * np.pi, 90)
1411 >>> data_smth = lut(fine_lats, fine_lons)
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()
1423 """
1425 def __init__(self, theta, phi, r, w=None, s=0., eps=1E-16):
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)')
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)
1448 self.fp = fp
1449 self.tck = tt_[:nt_], tp_[:np_], c[:(nt_ - 4) * (np_ - 4)]
1450 self.degrees = (3, 3)
1453class LSQSphereBivariateSpline(SphereBivariateSpline):
1454 """
1455 Weighted least-squares bivariate spline approximation in spherical
1456 coordinates.
1458 Determines a smooth bicubic spline according to a given
1459 set of knots in the `theta` and `phi` directions.
1461 .. versionadded:: 0.11.0
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.
1480 Notes
1481 -----
1482 For more information, see the FITPACK_ site about this function.
1484 .. _FITPACK: http://www.netlib.org/dierckx/sphere.f
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):
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)
1502 We need to set up the interpolator object. Here, we must also specify the
1503 coordinates of the knots to use.
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)
1515 As a first test, we'll see what the algorithm returns when run on the
1516 input coordinates
1518 >>> data_orig = lut(theta, phi)
1520 Finally we interpolate the data to a finer grid
1522 >>> fine_lats = np.linspace(0., np.pi, 70)
1523 >>> fine_lons = np.linspace(0., 2*np.pi, 90)
1525 >>> data_lsq = lut(fine_lats, fine_lons)
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()
1537 """
1539 def __init__(self, theta, phi, r, tt, tp, w=None, eps=1E-16):
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)')
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)
1570 self.fp = fp
1571 self.tck = tt_, tp_, c
1572 self.degrees = (3, 3)
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."""
1604class RectSphereBivariateSpline(SphereBivariateSpline):
1605 """
1606 Bivariate spline approximation over a rectangular mesh on a sphere.
1608 Can be used for smoothing data.
1610 .. versionadded:: 0.11.0
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.
1645 See Also
1646 --------
1647 RectBivariateSpline : bivariate spline approximation over a rectangular
1648 mesh
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.
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.
1660 For more information, see the FITPACK_ site about this function.
1662 .. _FITPACK: http://www.netlib.org/dierckx/spgrid.f
1664 Examples
1665 --------
1666 Suppose we have global data on a coarse grid
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
1673 We want to interpolate it to a global one-degree grid
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)
1679 We need to set up the interpolator object
1681 >>> from scipy.interpolate import RectSphereBivariateSpline
1682 >>> lut = RectSphereBivariateSpline(lats, lons, data)
1684 Finally we interpolate the data. The `RectSphereBivariateSpline` object
1685 only takes 1-D arrays as input, therefore we need to do some reshaping.
1687 >>> data_interp = lut.ev(new_lats.ravel(),
1688 ... new_lons.ravel()).reshape((360, 180)).T
1690 Looking at the original and the interpolated data, one can see that the
1691 interpolant reproduces the original data very well:
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()
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``.
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.
1720 The interpolation results for different values of ``s`` give some insight
1721 into this process:
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()
1734 """
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)
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]
1758 if r1 is None:
1759 ider[2] = -1
1760 else:
1761 ider[2] = pole_exact[1]
1763 ider[1], ider[3] = pole_flat
1765 u, v = np.ravel(u), np.ravel(v)
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 ')
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')
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')
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')
1793 if not s >= 0.0:
1794 raise ValueError('s should be positive')
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)
1803 if ier not in [0, -1, -2]:
1804 msg = _spfit_messages.get(ier, 'ier=%s' % (ier))
1805 raise ValueError(msg)
1807 self.fp = fp
1808 self.tck = tu[:nu], tv[:nv], c[:(nu - 4) * (nv-4)]
1809 self.degrees = (3, 3)