fork(1) download
copy 
  1. r"""
  2. Show how to compute a derivative spline.
  3.  
  4. Scipy's splines are represented in terms of the standard B-spline basis
  5. functions.  In short, a spline of degree ``k`` is represented in terms of the
  6. knots ``t`` and coefficients ``c`` by:
  7.  
  8. .. math::
  9.  
  10.     s(x) = \sum_{j=-\infty}^\infty c_{j} B^k_{j}(x)
  11.     \\
  12.     B^0_i(x) = 1, \text{if $t_i \le x < t_{i+1}$, otherwise $0$,}
  13.     \\
  14.     B^k_i(x) = \frac{x - t_i}{t_{i+k} - t_i} B^{k-1}_i(x)
  15.              + \frac{t_{i+k+1} - x}{t_{i+k+1} - t_{i+1}} B^{k-1}_{i+1}(x)
  16.  
  17. where :math:`c_j` is nonzero only for ``0 <= j <= N`` and the first
  18. and last `k` knots are at the same position and only used to set up
  19. boundary conditions; terms with vanishing denominators are omitted.
  20.  
  21. One can follow standard spline textbooks here, and work out that by
  22. differentiating this, one obtains another B-spline, of one order
  23. lower:
  24.  
  25. .. math::
  26.  
  27.    s'(x) = \sum_{j=-\infty}^\infty d_j B^{k-1}_j(x)
  28.    \\
  29.    d_j = k \frac{c_{j+1} - c_{j}}{t_{j+k+1} - t_{j+1}}
  30.  
  31. Care needs to be paid at the endpoints: the first knot
  32. should be thrown away since the order is reduced by one.
  33.  
  34. """
  35. import numpy as np
  36. from scipy.interpolate import splev
  37. from scipy.interpolate import splrep
  38. from scipy.interpolate import UnivariateSpline
  39.  
  40. import matplotlib.pyplot as plt
  41.  
  42.  
  43. class MyUnivariateSpline(UnivariateSpline):
  44.     @classmethod
  45.     def _from_tck(cls, t, c, k):
  46.         self = cls.__new__(cls)
  47.         self._eval_args = t, c, k
  48.         #_data == x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier
  49.         self._data = [None,None,None,None,None,k,None,len(t),t,
  50.                       c,None,None,None,None]
  51.         return self
  52.  
  53.     def derivative_spline(self):
  54.         """
  55.         Compute the derivative spline of a spline in FITPACK tck
  56.         representation.
  57.         """
  58.         t, c, k = self._eval_args
  59.         if k <= 0:
  60.             raise ValueError("Cannot differentiate order 0 spline")
  61.         # Compute the denominator in the differentiation formula.
  62.         dt = t[k+1:-1] - t[1:-k-1]
  63.         # Compute the new coefficients
  64.         d = (c[1:-1-k] - c[:-2-k]) * k / dt
  65.         # Adjust knots
  66.         t2 = t[1:-1]
  67.         # Pad coefficient array to same size as knots (FITPACK convention)
  68.         d = np.r_[d, [0]*k]
  69.         # Done, return a new spline
  70.         return self._from_tck(t2, d, k-1)
  71.  
  72. def main():
  73.     y = [-1, 5, 6, 4.3, 2, 5.2, 8, 5.6, 1]
  74.     x = range(len(y))
  75.  
  76.     # Fit a spline (must use order 4, since roots() works only for order 3)
  77.     s0 = MyUnivariateSpline(x, y, s=0, k=5)
  78.  
  79.     # Compute the derivative splines of the fitted spline
  80.     s1 = s0.derivative_spline()
  81.     s2 = s1.derivative_spline()
  82.  
  83.     xs = np.linspace(0, len(y), 16*len(y))
  84.     ys0 = s0(xs)
  85.     ys1 = s1(xs)
  86.     ys2 = s2(xs)
  87.  
  88.     # Roots (works only with order 3 splines)
  89.     inflection_points = s2.roots()
  90.     print inflection_points
  91.     print "inflection points:  x =", zip(inflection_points, s0(inflection_points))
  92.  
  93.     plt.ylim(-10, 15)
  94.     plt.plot(x, y, 'o',
  95.             inflection_points, s0(inflection_points), '*',
  96.         	xs, ys0, '-',
  97.             xs, ys1, '-',
  98.             xs, ys2, '-')
  99.     plt.savefig('out.png', dpi=96)
  100.     plt.show()
  101.  
  102. if __name__ == "__main__":
  103.     main()
  104.  
Not running #stdin #stdout 0s 0KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/qKja7X
language:
Python (cpython 2.7.16)
created:
11 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!