[Scipy-svn] r4175 - in trunk/scipy/interpolate: . tests

scipy-svn@scip... scipy-svn@scip...
Thu Apr 24 19:11:36 CDT 2008


Author: peridot
Date: 2008-04-24 19:11:33 -0500 (Thu, 24 Apr 2008)
New Revision: 4175

Added:
   trunk/scipy/interpolate/polyint.py
   trunk/scipy/interpolate/tests/test_polyint.py
Modified:
   trunk/scipy/interpolate/__init__.py
   trunk/scipy/interpolate/info.py
   trunk/scipy/interpolate/interpolate.py
Log:
Polynomial interpolation classes. (Python-only implementation.)


Modified: trunk/scipy/interpolate/__init__.py
===================================================================
--- trunk/scipy/interpolate/__init__.py	2008-04-24 22:53:44 UTC (rev 4174)
+++ trunk/scipy/interpolate/__init__.py	2008-04-25 00:11:33 UTC (rev 4175)
@@ -12,6 +12,8 @@
 
 from rbf import Rbf
 
+from polyint import *
+
 __all__ = filter(lambda s:not s.startswith('_'),dir())
 from scipy.testing.pkgtester import Tester
 test = Tester().test

Modified: trunk/scipy/interpolate/info.py
===================================================================
--- trunk/scipy/interpolate/info.py	2008-04-24 22:53:44 UTC (rev 4174)
+++ trunk/scipy/interpolate/info.py	2008-04-25 00:11:33 UTC (rev 4175)
@@ -24,18 +24,27 @@
 
   SmoothBivariateSpline
 
-Interpolation Class
+Interpolation Classes (univariate)
 
   interp1d -- Create a class whose instances can linearly interpolate
                to compute unknown values of a univariate function.
+  BarycentricInterpolator -- Compute with a numerically-stable version 
+               of the Lagrange interpolating polynomial.
+  KroghInterpolator -- Compute with the Hermite interpolating polynomial
+               (allows the specification of derivatives at some points).
+  PiecewisePolynomial -- Spline that is specified by giving positions and
+               derivatives at every knot; allows high orders and 
+               efficient appending.
+
+Interpolation Classes (multivariate)
+
   interp2d -- Create a class whose instances can interpolate
                to compute unknown values of a bivariate function.
   Rbf -- Apply Radial Basis Functions to interpolate scattered N-D data.
 
 Additional tools
 
-  lagrange -- Compute the Lagrange interpolating polynomial
-
+  lagrange -- Compute the Lagrange interpolating polynomial.
 """
 
 postpone_import = 1

Modified: trunk/scipy/interpolate/interpolate.py
===================================================================
--- trunk/scipy/interpolate/interpolate.py	2008-04-24 22:53:44 UTC (rev 4174)
+++ trunk/scipy/interpolate/interpolate.py	2008-04-25 00:11:33 UTC (rev 4175)
@@ -25,6 +25,9 @@
 
 def lagrange(x, w):
     """Return the Lagrange interpolating polynomial of the data-points (x,w)
+
+    Warning: This implementation is numerically unstable; do not expect to
+    be able to use more than about 20 points even if they are chosen optimally.
     """
     M = len(x)
     p = poly1d(0.0)

Added: trunk/scipy/interpolate/polyint.py
===================================================================
--- trunk/scipy/interpolate/polyint.py	2008-04-24 22:53:44 UTC (rev 4174)
+++ trunk/scipy/interpolate/polyint.py	2008-04-25 00:11:33 UTC (rev 4175)
@@ -0,0 +1,618 @@
+import numpy as np
+from scipy import factorial
+from numpy import poly1d
+
+__all__ = ["KroghInterpolator", "BarycentricInterpolator", "PiecewisePolynomial"]
+
+class KroghInterpolator:
+    """The interpolating polynomial for a set of points
+
+    Constructs a polynomial that passes through a given set of points,         
+    optionally with specified derivatives at those points.
+    Allows evaluation of the polynomial and all its derivatives.
+    For reasons of numerical stability, this function does not compute
+    the coefficients of the polynomial, although they can be obtained
+    by evaluating all the derivatives.
+
+    Be aware that the algorithms implemented here are not necessarily 
+    the most numerically stable known. Moreover, even in a world of 
+    exact computation, unless the x coordinates are chosen very 
+    carefully - Chebyshev zeros (e.g. cos(i*pi/n)) are a good choice - 
+    polynomial interpolation itself is a very ill-conditioned process 
+    due to the Runge phenomenon. In general, even with well-chosen 
+    x values, degrees higher than about thirty cause problems with
+    numerical instability in this code.
+
+    Based on Krogh 1970, "Efficient Algorithms for Polynomial Interpolation 
+    and Numerical Differentiation"
+    """
+    def __init__(self, xi, yi):
+        """Construct an interpolator passing through the specified points
+
+        The polynomial passes through all the pairs (xi,yi). One may additionally
+        specify a number of derivatives at each point xi; this is done by 
+        repeating the value xi and specifying the derivatives as successive 
+        yi values.  
+
+        Parameters
+        ----------
+        xi : array-like, length N
+            known x-coordinates 
+        yi : array-like, N by R
+            known y-coordinates, interpreted as vectors of length R,
+            or scalars if R=1
+
+        Example
+        -------
+        To produce a polynomial that is zero at 0 and 1 and has 
+        derivative 2 at 0, call
+
+        >>> KroghInterpolator([0,0,1],[0,2,0])
+        """
+        self.xi = np.asarray(xi)
+        self.yi = np.asarray(yi)
+        if len(self.yi.shape)==1:
+            self.vector_valued = False
+            self.yi = self.yi[:,np.newaxis]
+        elif len(self.yi.shape)>2:
+            raise ValueError, "y coordinates must be either scalars or vectors"
+        else:
+            self.vector_valued = True
+
+        n = len(xi)
+        self.n = n
+        nn, r = self.yi.shape
+        if nn!=n:
+            raise ValueError, "%d x values provided and %d y values; must be equal" % (n, nn)
+        self.r = r
+
+        c = np.zeros((n+1,r))
+        c[0] = yi[0]
+        Vk = np.zeros((n,r))
+        for k in xrange(1,n):
+            s = 0
+            while s<=k and xi[k-s]==xi[k]:
+                s += 1
+            s -= 1
+            Vk[0] = yi[k]/float(factorial(s))
+            for i in xrange(k-s):
+                assert xi[i]!=xi[k]
+                if s==0:
+                    Vk[i+1] = (c[i]-Vk[i])/(xi[i]-xi[k])
+                else:
+                    Vk[i+1] = (Vk[i+1]-Vk[i])/(xi[i]-xi[k])
+            c[k] = Vk[k-s]
+        self.c = c
+
+    def __call__(self,x):
+        """Evaluate the polynomial at the point x
+        
+        Parameters
+        ----------
+        x : scalar or array-like of length N
+
+        Returns
+        -------
+        y : scalar, array of length R, array of length N, or array of length N by R
+            If x is a scalar, returns either a vector or a scalar depending on
+            whether the interpolator is vector-valued or scalar-valued.
+            If x is a vector, returns a vector of values.
+        """
+        if np.isscalar(x):
+            scalar = True
+            m = 1
+        else: 
+            scalar = False
+            m = len(x)
+        x = np.asarray(x)
+
+        n = self.n
+        pi = 1
+        p = np.zeros((m,self.r))
+        p += self.c[0,np.newaxis,:]
+        for k in xrange(1,n):
+            w = x - self.xi[k-1]
+            pi = w*pi
+            p = p + np.multiply.outer(pi,self.c[k])
+        if not self.vector_valued:
+            if scalar:
+                return p[0,0]
+            else:
+                return p[:,0]
+        else:
+            if scalar:
+                return p[0]
+            else:
+                return p
+
+    def derivatives(self,x,der=None): 
+        """Evaluate many derivatives of the polynomial at the point x
+
+        Produce an array of all derivative values at the point x.
+        
+        Parameters
+        ----------
+        x : scalar or array-like of length N
+            Point or points at which to evaluate the derivatives
+        der : None or integer
+            How many derivatives to extract; None for all potentially
+            nonzero derivatives (that is a number equal to the number
+            of points). This number includes the function value as 0th
+            derivative.
+        Returns
+        -------
+        d : array
+            If the interpolator's values are R-dimensional then the 
+            returned array will be der by N by R. If x is a scalar, 
+            the middle dimension will be dropped; if R is 1 then the 
+            last dimension will be dropped.
+
+        Example
+        -------
+        >>> KroghInterpolator([0,0,0],[1,2,3]).derivatives(0)
+        array([1.0,2.0,3.0])
+        >>> KroghInterpolator([0,0,0],[1,2,3]).derivatives([0,0])
+        array([[1.0,1.0],
+               [2.0,2.0],
+               [3.0,3.0]])
+        """
+        if np.isscalar(x):
+            scalar = True
+            m = 1
+        else: 
+            scalar = False
+            m = len(x)
+        x = np.asarray(x)
+
+        n = self.n
+        r = self.r
+
+        if der is None:
+            der = self.n
+        dern = min(self.n,der)
+        pi = np.zeros((n,m))
+        w = np.zeros((n,m))
+        pi[0] = 1
+        p = np.zeros((m,self.r))
+        p += self.c[0,np.newaxis,:]
+
+        for k in xrange(1,n):
+            w[k-1] = x - self.xi[k-1]
+            pi[k] = w[k-1]*pi[k-1]
+            p += np.multiply.outer(pi[k],self.c[k])
+
+        cn = np.zeros((max(der,n+1),m,r))
+        cn[:n+1,...] += self.c[:n+1,np.newaxis,:]
+        cn[0] = p
+        for k in xrange(1,n):
+            for i in xrange(1,n-k+1):
+                pi[i] = w[k+i-1]*pi[i-1]+pi[i]
+                cn[k] = cn[k]+pi[i,:,np.newaxis]*cn[k+i]
+            cn[k]*=factorial(k)
+
+        cn[n,...] = 0
+        if not self.vector_valued:
+            if scalar:
+                return cn[:der,0,0]
+            else:
+                return cn[:der,:,0]
+        else:
+            if scalar:
+                return cn[:der,0]
+            else:
+                return cn[:der]
+    def derivative(self,x,der): 
+        """Evaluate one derivative of the polynomial at the point x
+
+        Parameters
+        ----------
+        x : scalar or array-like of length N
+            Point or points at which to evaluate the derivatives
+        der : None or integer
+            Which derivative to extract. This number includes the 
+            function value as 0th derivative.
+        Returns
+        -------
+        d : array
+            If the interpolator's values are R-dimensional then the 
+            returned array will be N by R. If x is a scalar, 
+            the middle dimension will be dropped; if R is 1 then the 
+            last dimension will be dropped.
+
+        Notes
+        -----
+        This is computed by evaluating all derivatives up to the desired 
+        one and then discarding the rest.
+        """
+        return self.derivatives(x,der=der+1)[der]
+
+
+class BarycentricInterpolator:
+    """The interpolating polynomial for a set of points
+
+    Constructs a polynomial that passes through a given set of points.
+    Allows evaluation of the polynomial, efficient changing of the y 
+    values to be interpolated, and updating by adding more x values.  
+    For reasons of numerical stability, this function does not compute 
+    the coefficients of the polynomial.
+
+    This class uses a "barycentric interpolation" method that treats
+    the problem as a special case of rational function interpolation.
+    This algorithm is quite stable, numerically, but even in a world of 
+    exact computation, unless the x coordinates are chosen very 
+    carefully - Chebyshev zeros (e.g. cos(i*pi/n)) are a good choice - 
+    polynomial interpolation itself is a very ill-conditioned process 
+    due to the Runge phenomenon.
+
+    Based on Berrut and Trefethen 2004, "Barycentric Lagrange Interpolation".
+    """
+    def __init__(self, xi, yi=None):
+        """Construct an object capable of interpolating functions sampled at xi
+
+        The values yi need to be provided before the function is evaluated,
+        but none of the preprocessing depends on them, so rapid updates
+        are possible.
+
+        Parameters
+        ----------
+        xi : array-like of length N
+            The x coordinates of the points the polynomial should pass through
+        yi : array-like N by R or None
+            The y coordinates of the points the polynomial should pass through;
+            if R>1 the polynomial is vector-valued. If None the y values 
+            will be supplied later.
+        """
+        self.n = len(xi)
+        self.xi = np.asarray(xi)
+        if yi is not None and len(yi)!=len(self.xi):
+            raise ValueError, "yi dimensions do not match xi dimensions"
+        self.set_yi(yi)
+        self.wi = np.zeros(self.n)
+        self.wi[0] = 1
+        for j in xrange(1,self.n):
+            self.wi[:j]*=(self.xi[j]-self.xi[:j])
+            self.wi[j] = np.multiply.reduce(self.xi[:j]-self.xi[j])
+        self.wi**=-1
+            
+    def set_yi(self, yi):
+        """Update the y values to be interpolated
+
+        The barycentric interpolation algorithm requires the calculation 
+        of weights, but these depend only on the xi. The yi can be changed 
+        at any time.
+
+        Parameters
+        ----------
+        yi : array-like N by R
+            The y coordinates of the points the polynomial should pass through;
+            if R>1 the polynomial is vector-valued. If None the y values 
+            will be supplied later.
+        """
+        if yi is None:
+            self.yi = None
+            return
+        yi = np.asarray(yi)
+        if len(yi.shape)==1:
+            self.vector_valued = False
+            yi = yi[:,np.newaxis]
+        elif len(yi.shape)>2:
+            raise ValueError, "y coordinates must be either scalars or vectors"
+        else:
+            self.vector_valued = True
+
+        n, r = yi.shape
+        if n!=len(self.xi):
+            raise ValueError, "yi dimensions do not match xi dimensions"
+        self.yi = yi
+        self.r = r
+        
+
+    def add_xi(self, xi, yi=None):
+        """Add more x values to the set to be interpolated
+
+        The barycentric interpolation algorithm allows easy updating by 
+        adding more points for the polynomial to pass through.
+
+        Parameters
+        ----------
+        xi : array-like of length N1
+            The x coordinates of the points the polynomial should pass through
+        yi : array-like N1 by R or None
+            The y coordinates of the points the polynomial should pass through;
+            if R>1 the polynomial is vector-valued. If None the y values 
+            will be supplied later. The yi should be specified if and only if 
+            the interpolator has y values specified.
+        """
+        if yi is not None:
+            if self.yi is None:
+                raise ValueError, "No previous yi value to update!"
+            yi = np.asarray(yi)
+            if len(yi.shape)==1:
+                if self.vector_valued:
+                    raise ValueError, "Cannot extend dimension %d y vectors with scalars" % self.r
+                yi = yi[:,np.newaxis]
+            elif len(yi.shape)>2:
+                raise ValueError, "y coordinates must be either scalars or vectors"
+            else:
+                n, r = yi.shape
+                if r!=self.r:
+                    raise ValueError, "Cannot extend dimension %d y vectors with dimension %d y vectors" % (self.r, r)
+
+            self.yi = np.vstack((self.yi,yi))
+        else:
+            if self.yi is not None:
+                raise ValueError, "No update to yi provided!"
+        old_n = self.n
+        self.xi = np.concatenate((self.xi,xi))
+        self.n = len(self.xi)
+        self.wi**=-1
+        old_wi = self.wi
+        self.wi = np.zeros(self.n)
+        self.wi[:old_n] = old_wi
+        for j in xrange(old_n,self.n):
+            self.wi[:j]*=(self.xi[j]-self.xi[:j])
+            self.wi[j] = np.multiply.reduce(self.xi[:j]-self.xi[j])
+        self.wi**=-1
+
+    def __call__(self, x):
+        """Evaluate the interpolating polynomial at the points x
+
+        Parameters
+        ----------
+        x : scalar or array-like of length M
+
+        Returns
+        -------
+        y : scalar or array-like of length R or length M or M by R
+            The shape of y depends on the shape of x and whether the
+            interpolator is vector-valued or scalar-valued.
+
+        Notes
+        -----
+        Currently the code computes an outer product between x and the 
+        weights, that is, it constructs an intermediate array of size 
+        N by M, where N is the degree of the polynomial.
+        """
+        scalar = np.isscalar(x)
+        x = np.atleast_1d(x)
+        c = np.subtract.outer(x,self.xi)
+        z = c==0
+        c[z] = 1
+        c = self.wi/c
+        p = np.dot(c,self.yi)/np.sum(c,axis=-1)[:,np.newaxis]
+        i, j = np.nonzero(z)
+        p[i] = self.yi[j]
+        if not self.vector_valued:
+            if scalar:
+                return p[0,0]
+            else:
+                return p[:,0]
+        else:
+            if scalar:
+                return p[0]
+            else:
+                return p
+
+
+class PiecewisePolynomial:
+    """Piecewise polynomial curve specified by points and derivatives
+
+    This class represents a curve that is a piecewise polynomial. It 
+    passes through a list of points and has specified derivatives at 
+    each point. The degree of the polynomial may very from segment to 
+    segment, as may the number of derivatives available. The degree 
+    should not exceed about thirty.
+
+    Appending points to the end of the curve is efficient.
+    """
+    def __init__(self, xi, yi, orders=None, direction=None):
+        """Construct a piecewise polynomial
+
+        Parameters
+        ----------
+        xi : array-like of length N
+            a sorted list of x-coordinates
+        yi : list of lists of length N
+            yi[i] is the list of derivatives known at xi[i]
+        orders : list of integers, or integer
+            a list of polynomial orders, or a single universal order
+        direction : {None, 1, -1}
+            indicates whether the xi are increasing or decreasing
+            +1 indicates increasing
+            -1 indicates decreasing
+            None indicates that it should be deduced from the first two xi
+
+        Notes
+        -----
+        If orders is None, or orders[i] is None, then the degree of the 
+        polynomial segment is exactly the degree required to match all i
+        available derivatives at both endpoints. If orders[i] is not None, 
+        then some derivatives will be ignored. The code will try to use an     
+        equal number of derivatives from each end; if the total number of 
+        derivatives needed is odd, it will prefer the rightmost endpoint. If 
+        not enough derivatives are available, an exception is raised.
+        """
+        self.xi = [xi[0]]
+        self.yi = [yi[0]]
+        self.n = 1
+        
+        try:
+            self.r = len(yi[0][0])
+        except TypeError:
+            self.r = 1
+
+        self.n = 1
+        self.direction = direction
+        self.orders = []
+        self.polynomials = []
+        self.extend(xi[1:],yi[1:],orders)
+
+    def _make_polynomial(self,x1,y1,x2,y2,order,direction):
+        """Construct the interpolating polynomial object
+        
+        Deduces the number of derivatives to match at each end
+        from order and the number of derivatives available. If
+        possible it uses the same number of derivatives from 
+        each end; if the number is odd it tries to take the
+        extra one from y2. In any case if not enough derivatives 
+        are available at one end or another it draws enough to
+        make up the total from the other end.
+        """
+        n = order+1
+        n1 = min(n//2,len(y1))
+        n2 = min(n-n1,len(y2))
+        n1 = min(n-n2,len(y1))
+        if n1+n2!=n:
+            raise ValueError, "Point %g has %d derivatives, point %g has %d derivatives, but order %d requested" % (x1, len(y1), x2, len(y2), order)
+        assert n1<=len(y1)
+        assert n2<=len(y2)
+
+        xi = np.zeros(n)
+        if self.r==1:
+            yi = np.zeros(n)
+        else:
+            yi = np.zeros((n,self.r))
+        xi[:n1] = x1
+        yi[:n1] = y1[:n1]
+        xi[n1:] = x2
+        yi[n1:] = y2[:n2]
+
+        return KroghInterpolator(xi,yi)
+
+    def append(self, xi, yi, order=None):
+        """Append a single point with derivatives to the PiecewisePolynomial
+        
+        Parameters
+        ----------
+        xi : float
+        yi : array-like
+            yi is the list of derivatives known at xi
+        order : integer or None
+            a polynomial order, or instructions to use the highest 
+            possible order
+        """
+        if self.direction is None:
+            self.direction = np.sign(xi-self.xi[-1])
+        elif (xi-self.xi[-1])*self.direction < 0: 
+            raise ValueError, "x coordinates must be in the %d direction: %s" % (self.direction, self.xi)
+        self.xi.append(xi)
+        self.yi.append(yi)
+
+        if order is None:
+            n1 = len(self.yi[-2])
+            n2 = len(self.yi[-1])
+            n = n1+n2
+            order = n-1
+
+        self.orders.append(order)
+        self.polynomials.append(self._make_polynomial(
+            self.xi[-2], self.yi[-2],
+            self.xi[-1], self.yi[-1],
+            order, self.direction))
+        self.n += 1
+
+
+    def extend(self, xi, yi, orders=None):
+        """Extend the PiecewisePolynomial by a list of points
+
+        Parameters
+        ----------
+        xi : array-like of length N1
+            a sorted list of x-coordinates
+        yi : list of lists of length N1
+            yi[i] is the list of derivatives known at xi[i]
+        orders : list of integers, or integer
+            a list of polynomial orders, or a single universal order
+        direction : {None, 1, -1}
+            indicates whether the xi are increasing or decreasing
+            +1 indicates increasing
+            -1 indicates decreasing
+            None indicates that it should be deduced from the first two xi
+        """
+
+        for i in xrange(len(xi)):
+            if orders is None or np.isscalar(orders):
+                self.append(xi[i],yi[i],orders)
+            else:
+                self.append(xi[i],yi[i],orders[i])
+
+    def __call__(self, x):
+        """Evaluate the piecewise polynomial
+
+        Parameters
+        ----------
+        x : scalar or array-like of length N
+        
+        Returns
+        -------
+        y : scalar or array-like of length R or length N or N by R
+        """
+        if np.isscalar(x):
+            pos = np.clip(np.searchsorted(self.xi, x) - 1, 0, self.n-2)
+            y = self.polynomials[pos](x)
+        else:
+            x = np.asarray(x)
+            m = len(x)
+            pos = np.clip(np.searchsorted(self.xi, x) - 1, 0, self.n-2)
+            if self.r>1:
+                y = np.zeros((m,self.r))
+            else:
+                y = np.zeros(m)
+            for i in xrange(self.n-1):
+                c = pos==i
+                y[c] = self.polynomials[i](x[c])
+        return y
+
+    def derivative(self, x, der):
+        """Evaluate a derivative of the piecewise polynomial
+
+        Parameters
+        ----------
+        x : scalar or array-like of length N
+        der : integer
+            which single derivative to extract
+        
+        Returns
+        -------
+        y : scalar or array-like of length R or length N or N by R
+
+        Notes
+        -----
+        This currently computes all derivatives of the curve segment 
+        containing each x but returns only one. This is because the 
+        number of nonzero derivatives that a segment can have depends
+        on the degree of the segment, which may vary.
+        """
+        return self.derivatives(x,der=der+1)[der]
+
+    def derivatives(self, x, der):
+        """Evaluate a derivative of the piecewise polynomial
+
+        Parameters
+        ----------
+        x : scalar or array-like of length N
+        der : integer
+            how many derivatives (including the function value as 
+            0th derivative) to extract
+        
+        Returns
+        -------
+        y : array-like of shape der by R or der by N or der by N by R
+        
+        """
+        if np.isscalar(x):
+            pos = np.clip(np.searchsorted(self.xi, x) - 1, 0, self.n-2)
+            y = self.polynomials[pos].derivatives(x,der=der)
+        else:
+            x = np.asarray(x)
+            m = len(x)
+            pos = np.clip(np.searchsorted(self.xi, x) - 1, 0, self.n-2)
+            if self.r>1:
+                y = np.zeros((der,m,self.r))
+            else:
+                y = np.zeros((der,m))
+            for i in xrange(self.n-1):
+                c = pos==i
+                y[:,c] = self.polynomials[i].derivatives(x[c],der=der)
+        return y
+    # FIXME: provide multiderivative finder

Added: trunk/scipy/interpolate/tests/test_polyint.py
===================================================================
--- trunk/scipy/interpolate/tests/test_polyint.py	2008-04-24 22:53:44 UTC (rev 4174)
+++ trunk/scipy/interpolate/tests/test_polyint.py	2008-04-25 00:11:33 UTC (rev 4175)
@@ -0,0 +1,231 @@
+
+from scipy.testing import *
+from scipy.interpolate import KroghInterpolator, \
+        BarycentricInterpolator, PiecewisePolynomial
+import scipy
+import numpy as np
+from scipy.interpolate import splrep, splev
+
+class CheckKrogh(TestCase):
+    def setUp(self):
+        self.true_poly = scipy.poly1d([-2,3,1,5,-4])
+        self.test_xs = np.linspace(-1,1,100) 
+        self.xs = np.linspace(-1,1,5)
+        self.ys = self.true_poly(self.xs)
+
+    def test_lagrange(self):
+        P = KroghInterpolator(self.xs,self.ys)
+        assert_almost_equal(self.true_poly(self.test_xs),P(self.test_xs))
+    def test_scalar(self):
+        P = KroghInterpolator(self.xs,self.ys)
+        assert_almost_equal(self.true_poly(7),P(7))
+
+    def test_derivatives(self):
+        P = KroghInterpolator(self.xs,self.ys)
+        D = P.derivatives(self.test_xs)
+        for i in xrange(D.shape[0]):
+            assert_almost_equal(self.true_poly.deriv(i)(self.test_xs),
+                                D[i])
+    def test_low_derivatives(self):
+        P = KroghInterpolator(self.xs,self.ys)
+        D = P.derivatives(self.test_xs,len(self.xs)+2)
+        for i in xrange(D.shape[0]):
+            assert_almost_equal(self.true_poly.deriv(i)(self.test_xs),
+                                D[i])
+    def test_derivative(self):
+        P = KroghInterpolator(self.xs,self.ys)
+        m = 10
+        r = P.derivatives(self.test_xs,m)
+        for i in xrange(m):
+            assert_almost_equal(P.derivative(self.test_xs,i),r[i])
+    def test_high_derivative(self):
+        P = KroghInterpolator(self.xs,self.ys)
+        for i in xrange(len(self.xs),2*len(self.xs)):
+            assert_almost_equal(P.derivative(self.test_xs,i),
+                                np.zeros(len(self.test_xs)))
+    def test_hermite(self):
+        xs = [0,0,0,1,1,1,2]
+        ys = [self.true_poly(0),
+              self.true_poly.deriv(1)(0),
+              self.true_poly.deriv(2)(0),
+              self.true_poly(1),
+              self.true_poly.deriv(1)(1),
+              self.true_poly.deriv(2)(1),
+              self.true_poly(2)]
+        P = KroghInterpolator(self.xs,self.ys)
+        assert_almost_equal(self.true_poly(self.test_xs),P(self.test_xs))
+
+    def test_vector(self):
+        xs = [0, 1, 2]
+        ys = np.array([[0,1],[1,0],[2,1]])
+        P = KroghInterpolator(xs,ys)
+        Pi = [KroghInterpolator(xs,ys[:,i]) for i in xrange(ys.shape[1])]
+        test_xs = np.linspace(-1,3,100)
+        assert_almost_equal(P(test_xs),
+                np.rollaxis(np.asarray([p(test_xs) for p in Pi]),-1))
+        assert_almost_equal(P.derivatives(test_xs),
+                np.transpose(np.asarray([p.derivatives(test_xs) for p in Pi]),
+                    (1,2,0)))
+
+    def test_empty(self):
+        P = KroghInterpolator(self.xs,self.ys)
+        assert_array_equal(P([]), [])
+    def test_shapes_scalarvalue(self):
+        P = KroghInterpolator(self.xs,self.ys)
+        assert_array_equal(np.shape(P(0)), ())
+        assert_array_equal(np.shape(P([0])), (1,))
+        assert_array_equal(np.shape(P([0,1])), (2,))
+
+    def test_shapes_scalarvalue_derivative(self):
+        P = KroghInterpolator(self.xs,self.ys)
+        n = P.n
+        assert_array_equal(np.shape(P.derivatives(0)), (n,))
+        assert_array_equal(np.shape(P.derivatives([0])), (n,1))
+        assert_array_equal(np.shape(P.derivatives([0,1])), (n,2))
+
+    def test_shapes_vectorvalue(self):
+        P = KroghInterpolator(self.xs,np.outer(self.ys,np.arange(3)))
+        assert_array_equal(np.shape(P(0)), (3,))
+        assert_array_equal(np.shape(P([0])), (1,3))
+        assert_array_equal(np.shape(P([0,1])), (2,3))
+
+    def test_shapes_1d_vectorvalue(self):
+        P = KroghInterpolator(self.xs,np.outer(self.ys,[1]))
+        assert_array_equal(np.shape(P(0)), (1,))
+        assert_array_equal(np.shape(P([0])), (1,1))
+        assert_array_equal(np.shape(P([0,1])), (2,1))
+
+    def test_shapes_vectorvalue_derivative(self):
+        P = KroghInterpolator(self.xs,np.outer(self.ys,np.arange(3)))
+        n = P.n
+        assert_array_equal(np.shape(P.derivatives(0)), (n,3))
+        assert_array_equal(np.shape(P.derivatives([0])), (n,1,3))
+        assert_array_equal(np.shape(P.derivatives([0,1])), (n,2,3))
+
+
+
+class CheckBarycentric(TestCase):
+    def setUp(self):
+        self.true_poly = scipy.poly1d([-2,3,1,5,-4])
+        self.test_xs = np.linspace(-1,1,100) 
+        self.xs = np.linspace(-1,1,5)
+        self.ys = self.true_poly(self.xs)
+
+    def test_lagrange(self):
+        P = BarycentricInterpolator(self.xs,self.ys)
+        assert_almost_equal(self.true_poly(self.test_xs),P(self.test_xs))
+    def test_scalar(self):
+        P = BarycentricInterpolator(self.xs,self.ys)
+        assert_almost_equal(self.true_poly(7),P(7))
+
+    def test_delayed(self):
+        P = BarycentricInterpolator(self.xs)
+        P.set_yi(self.ys)
+        assert_almost_equal(self.true_poly(self.test_xs),P(self.test_xs))
+
+    def test_append(self):
+        P = BarycentricInterpolator(self.xs[:3],self.ys[:3])
+        P.add_xi(self.xs[3:],self.ys[3:])
+        assert_almost_equal(self.true_poly(self.test_xs),P(self.test_xs))
+
+    def test_vector(self):
+        xs = [0, 1, 2]
+        ys = np.array([[0,1],[1,0],[2,1]])
+        P = BarycentricInterpolator(xs,ys)
+        Pi = [BarycentricInterpolator(xs,ys[:,i]) for i in xrange(ys.shape[1])]
+        test_xs = np.linspace(-1,3,100)
+        assert_almost_equal(P(test_xs),
+                np.rollaxis(np.asarray([p(test_xs) for p in Pi]),-1))
+
+    def test_shapes_scalarvalue(self):
+        P = BarycentricInterpolator(self.xs,self.ys)
+        assert_array_equal(np.shape(P(0)), ())
+        assert_array_equal(np.shape(P([0])), (1,))
+        assert_array_equal(np.shape(P([0,1])), (2,))
+
+    def test_shapes_vectorvalue(self):
+        P = BarycentricInterpolator(self.xs,np.outer(self.ys,np.arange(3)))
+        assert_array_equal(np.shape(P(0)), (3,))
+        assert_array_equal(np.shape(P([0])), (1,3))
+        assert_array_equal(np.shape(P([0,1])), (2,3))
+    def test_shapes_1d_vectorvalue(self):
+        P = BarycentricInterpolator(self.xs,np.outer(self.ys,[1]))
+        assert_array_equal(np.shape(P(0)), (1,))
+        assert_array_equal(np.shape(P([0])), (1,1))
+        assert_array_equal(np.shape(P([0,1])), (2,1))
+
+
+class CheckPiecewise(TestCase):
+    def setUp(self):
+        self.tck = splrep([0,1,2,3,4,5],[0,10,-1,3,7,2],s=0)
+        self.test_xs = np.linspace(-1,6,100)
+        self.spline_ys = splev(self.test_xs, self.tck)
+        self.spline_yps = splev(self.test_xs, self.tck, der=1)
+        self.xi = np.unique(self.tck[0])
+        self.yi = [[splev(x,self.tck,der=j) for j in xrange(3)] for x in self.xi]
+
+    def test_construction(self):
+        P = PiecewisePolynomial(self.xi,self.yi,3)
+        assert_almost_equal(P(self.test_xs),self.spline_ys)
+    def test_scalar(self):
+        P = PiecewisePolynomial(self.xi,self.yi,3)
+        assert_almost_equal(P(self.test_xs[0]),self.spline_ys[0])
+        assert_almost_equal(P.derivative(self.test_xs[0],1),self.spline_yps[0])
+    def test_derivative(self):
+        P = PiecewisePolynomial(self.xi,self.yi,3)
+        assert_almost_equal(P.derivative(self.test_xs,1),self.spline_yps)
+    def test_derivatives(self):
+        P = PiecewisePolynomial(self.xi,self.yi,3)
+        m = 4
+        r = P.derivatives(self.test_xs,m)
+        #print r.shape, r
+        for i in xrange(m):
+            assert_almost_equal(P.derivative(self.test_xs,i),r[i])
+    def test_vector(self):
+        xs = [0, 1, 2]
+        ys = [[[0,1]],[[1,0],[-1,-1]],[[2,1]]]
+        P = PiecewisePolynomial(xs,ys)
+        Pi = [PiecewisePolynomial(xs,[[yd[i] for yd in y] for y in ys]) 
+            for i in xrange(len(ys[0][0]))]
+        test_xs = np.linspace(-1,3,100)
+        assert_almost_equal(P(test_xs),
+                np.rollaxis(np.asarray([p(test_xs) for p in Pi]),-1))
+        assert_almost_equal(P.derivative(test_xs,1),
+                np.transpose(np.asarray([p.derivative(test_xs,1) for p in Pi]),
+                    (1,0)))
+    def test_incremental(self):
+        P = PiecewisePolynomial([self.xi[0]], [self.yi[0]], 3)
+        for i in xrange(1,len(self.xi)):
+            P.append(self.xi[i],self.yi[i],3)
+        assert_almost_equal(P(self.test_xs),self.spline_ys)
+
+    def test_shapes_scalarvalue(self):
+        P = PiecewisePolynomial(self.xi,self.yi,4)
+        assert_array_equal(np.shape(P(0)), ())
+        assert_array_equal(np.shape(P([0])), (1,))
+        assert_array_equal(np.shape(P([0,1])), (2,))
+
+    def test_shapes_scalarvalue_derivative(self):
+        P = PiecewisePolynomial(self.xi,self.yi,4)
+        n = 4
+        assert_array_equal(np.shape(P.derivative(0,1)), ())
+        assert_array_equal(np.shape(P.derivative([0],1)), (1,))
+        assert_array_equal(np.shape(P.derivative([0,1],1)), (2,))
+
+    def test_shapes_vectorvalue(self):
+        yi = np.multiply.outer(np.asarray(self.yi),np.arange(3))
+        P = PiecewisePolynomial(self.xi,yi,4)
+        assert_array_equal(np.shape(P(0)), (3,))
+        assert_array_equal(np.shape(P([0])), (1,3))
+        assert_array_equal(np.shape(P([0,1])), (2,3))
+
+    def test_shapes_vectorvalue_derivative(self):
+        P = PiecewisePolynomial(self.xi,np.multiply.outer(self.yi,np.arange(3)),4)
+        n = 4
+        assert_array_equal(np.shape(P.derivative(0,1)), (3,))
+        assert_array_equal(np.shape(P.derivative([0],1)), (1,3))
+        assert_array_equal(np.shape(P.derivative([0,1],1)), (2,3))
+
+
+if __name__=='__main__':
+    nose.run(argv=['', __file__])



More information about the Scipy-svn mailing list