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] (if axis==0) 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

        """
        if self._y_axis == 0:
            # allow yi to be a ragged list
            for i in xrange(len(xi)):
                if orders is None or _isscalar(orders):
                    self.append(xi[i],yi[i],orders)
                else:
                    self.append(xi[i],yi[i],orders[i])
        else:
            preslice = (slice(None,None,None),) * self._y_axis
            for i in xrange(len(xi)):
                if orders is None or _isscalar(orders):
                    self.append(xi[i],yi[preslice + (i,)],orders)
                else:
                    self.append(xi[i],yi[preslice + (i,)],orders[i])

 def __add__(self, other):
     # First check if argument is a scalar
     if isscalarlike(other):
         new = dok_matrix(self.shape, dtype=self.dtype)
         # Add this scalar to every element.
         M, N = self.shape
         for i in xrange(M):
             for j in xrange(N):
                 aij = self.get((i, j), 0) + other
                 if aij != 0:
                     new[i, j] = aij
         # new.dtype.char = self.dtype.char
     elif isinstance(other, dok_matrix):
         if other.shape != self.shape:
             raise ValueError("matrix dimensions are not equal")
         # We could alternatively set the dimensions to the the largest of
         # the two matrices to be summed.  Would this be a good idea?
         new = dok_matrix(self.shape, dtype=self.dtype)
         new.update(self)
         for key in other.keys():
             new[key] += other[key]
     elif isspmatrix(other):
         csc = self.tocsc()
         new = csc + other
     elif isdense(other):
         new = self.todense() + other
     else:
         raise TypeError("data type not understood")
     return new

 def __radd__(self, other):
     # First check if argument is a scalar
     if isscalarlike(other):
         new = dok_matrix(self.shape, dtype=self.dtype)
         # Add this scalar to every element.
         M, N = self.shape
         for i in xrange(M):
             for j in xrange(N):
                 aij = self.get((i, j), 0) + other
                 if aij != 0:
                     new[i, j] = aij
     elif isinstance(other, dok_matrix):
         if other.shape != self.shape:
             raise ValueError("matrix dimensions are not equal")
         new = dok_matrix(self.shape, dtype=self.dtype)
         new.update(self)
         for key in other:
             new[key] += other[key]
     elif isspmatrix(other):
         csc = self.tocsc()
         new = csc + other
     elif isdense(other):
         new = other + self.todense()
     else:
         raise TypeError("data type not understood")
     return new

    def _evaluate_derivatives(self, x, der=None):
        n = self.n
        r = self.r

        if der is None:
            der = self.n
        pi = np.zeros((n, len(x)))
        w = np.zeros((n, len(x)))
        pi[0] = 1
        p = np.zeros((len(x), 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 += pi[k,:,np.newaxis]*self.c[k]

        cn = np.zeros((max(der,n+1), len(x), r), dtype=self.dtype)
        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
        return cn[:der]

    def __init__(self, xi, yi, axis=0):
        _Interpolator1DWithDerivatives.__init__(self, xi, yi, axis)

        self.xi = np.asarray(xi)
        self.yi = self._reshape_yi(yi)
        self.n, self.r = self.yi.shape

        c = np.zeros((self.n+1, self.r), dtype=self.dtype)
        c[0] = self.yi[0]
        Vk = np.zeros((self.n, self.r), dtype=self.dtype)
        for k in xrange(1,self.n):
            s = 0
            while s <= k and xi[k-s] == xi[k]:
                s += 1
            s -= 1
            Vk[0] = self.yi[k]/float(factorial(s))
            for i in xrange(k-s):
                if xi[i] == xi[k]:
                    raise ValueError("Elements if `xi` can't be equal.")
                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 lagrange(x, w):
    """
    Return a Lagrange interpolating polynomial.

    Given two 1-D arrays `x` and `w,` returns the Lagrange interpolating
    polynomial through the 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.

    Parameters
    ----------
    x : array_like
        `x` represents the x-coordinates of a set of datapoints.
    w : array_like
        `w` represents the y-coordinates of a set of datapoints, i.e. f(`x`).

    """
    M = len(x)
    p = poly1d(0.0)
    for j in xrange(M):
        pt = poly1d(w[j])
        for k in xrange(M):
            if k == j: continue
            fac = x[j]-x[k]
            pt *= poly1d([1.0,-x[k]])/fac
        p += pt
    return p

 def test_cascade(self):
     for J in xrange(1, 7):
         for i in xrange(1, 5):
             lpcoef = wavelets.daub(i)
             k = len(lpcoef)
             x, phi, psi = wavelets.cascade(lpcoef, J)
             assert_(len(x) == len(phi) == len(psi))
             assert_equal(len(x), (k - 1) * 2 ** J)

 def test_improvement(self):
     import time
     start = time.time()
     for i in xrange(100):
         quad(self.lib.sin, 0, 100)
     fast = time.time() - start
     start = time.time()
     for i in xrange(100):
         quad(math.sin, 0, 100)
     slow = time.time() - start
     assert_(fast < 0.5*slow, (fast, slow))

    def test_nd_simplex(self):
        # simple smoke test: triangulate a n-dimensional simplex
        for nd in xrange(2, 8):
            points = np.zeros((nd+1, nd))
            for j in xrange(nd):
                points[j,j] = 1.0
            points[-1,:] = 1.0

            tri = qhull.Delaunay(points)

            tri.vertices.sort()

            assert_equal(tri.vertices, np.arange(nd+1, dtype=np.int)[None,:])
            assert_equal(tri.neighbors, -1 + np.zeros((nd+1), dtype=np.int)[None,:])

 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_exponential(self):
     degree = 5
     p = approximate_taylor_polynomial(np.exp, 0, degree, 1, 15)
     for i in xrange(degree+1):
         assert_almost_equal(p(0),1)
         p = p.deriv()
     assert_almost_equal(p(0),0)

    def __getitem__(self, index):
        """Return the element(s) index=(i, j), where j may be a slice.
        This always returns a copy for consistency, since slices into
        Python lists return copies.
        """
        i, j = self._unpack_index(index)

        if isscalarlike(i) and isscalarlike(j):
            return self._get1(int(i), int(j))

        i, j = self._index_to_arrays(i, j)
        if i.size == 0:
            return lil_matrix((0,0), dtype=self.dtype)
        return self.__class__([[self._get1(int(i[ii, jj]), int(j[ii, jj])) for jj in
                                xrange(i.shape[1])] for ii in 
                               xrange(i.shape[0])])

    def test_concurrent_ok(self):
        f = lambda t, y: 1.0

        for k in xrange(3):
            for sol in ('vode', 'zvode', 'lsoda', 'dopri5', 'dop853'):
                r = ode(f).set_integrator(sol)
                r.set_initial_value(0, 0)

                r2 = ode(f).set_integrator(sol)
                r2.set_initial_value(0, 0)

                r.integrate(r.t + 0.1)
                r2.integrate(r2.t + 0.1)
                r2.integrate(r2.t + 0.1)

                assert_allclose(r.y, 0.1)
                assert_allclose(r2.y, 0.2)

            for sol in ('dopri5', 'dop853'):
                r = ode(f).set_integrator(sol)
                r.set_initial_value(0, 0)

                r2 = ode(f).set_integrator(sol)
                r2.set_initial_value(0, 0)

                r.integrate(r.t + 0.1)
                r.integrate(r.t + 0.1)
                r2.integrate(r2.t + 0.1)
                r.integrate(r.t + 0.1)
                r2.integrate(r2.t + 0.1)

                assert_allclose(r.y, 0.3)
                assert_allclose(r2.y, 0.2)

    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 _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.vector_valued:
                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

    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 _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.vector_valued:
                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 _setitem_setrow(self, row, data, j, xrow, xdata, xcols):
     if isinstance(j, slice):
         j = self._slicetoseq(j, self.shape[1])
     if issequence(j):
         if xcols == len(j):
             for jj, xi in zip(j, xrange(xcols)):
                 pos = bisect_left(xrow, xi)
                 if pos != len(xdata) and xrow[pos] == xi:
                     self._insertat2(row, data, jj, xdata[pos])
                 else:
                     self._insertat2(row, data, jj, 0)
         elif xcols == 1:           # OK, broadcast across row
             if len(xdata) > 0 and xrow[0] == 0:
                 val = xdata[0]
             else:
                 val = 0
             for jj in j:
                 self._insertat2(row, data, jj,val)
         else:
             raise IndexError('invalid index')
     elif np.isscalar(j):
         if not xcols == 1:
             raise ValueError('array dimensions are not compatible for copy')
         if len(xdata) > 0 and xrow[0] == 0:
             self._insertat2(row, data, j, xdata[0])
         else:
             self._insertat2(row, data, j, 0)
     else:
         raise ValueError('invalid column value: %s' % str(j))

    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 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 piecefuncgen(num):
        Mk = order // 2 - num
        if (Mk < 0):
            return 0  # final function is 0
        coeffs = [(1 - 2 * (k % 2)) * float(comb(order + 1, k, exact=1)) / fval
                  for k in xrange(Mk + 1)]
        shifts = [-bound - k for k in xrange(Mk + 1)]
        #print "Adding piece number %d with coeffs %s and shifts %s" \
        #      % (num, str(coeffs), str(shifts))

        def thefunc(x):
            res = 0.0
            for k in range(Mk + 1):
                res += coeffs[k] * (x + shifts[k]) ** order
            return res
        return thefunc

 def test_call(self):
     poly = []
     for n in xrange(5):
         poly.extend([x.strip() for x in
             ("""
             orth.jacobi(%(n)d,0.3,0.9)
             orth.sh_jacobi(%(n)d,0.3,0.9)
             orth.genlaguerre(%(n)d,0.3)
             orth.laguerre(%(n)d)
             orth.hermite(%(n)d)
             orth.hermitenorm(%(n)d)
             orth.gegenbauer(%(n)d,0.3)
             orth.chebyt(%(n)d)
             orth.chebyu(%(n)d)
             orth.chebyc(%(n)d)
             orth.chebys(%(n)d)
             orth.sh_chebyt(%(n)d)
             orth.sh_chebyu(%(n)d)
             orth.legendre(%(n)d)
             orth.sh_legendre(%(n)d)
             """ % dict(n=n)).split()
         ])
     olderr = np.seterr(all='ignore')
     try:
         for pstr in poly:
             p = eval(pstr)
             assert_almost_equal(p(0.315), np.poly1d(p)(0.315), err_msg=pstr)
     finally:
         np.seterr(**olderr)