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 _setdiag(self, values, k):
     M, N = self.shape
     if k < 0:
         if values.ndim == 0:
             # broadcast
             max_index = min(M+k, N)
             for i in xrange(max_index):
                 self[i - k, i] = values
         else:
             max_index = min(M+k, N, len(values))
             if max_index <= 0:
                 return
             for i,v in enumerate(values[:max_index]):
                 self[i - k, i] = v
     else:
         if values.ndim == 0:
             # broadcast
             max_index = min(M, N-k)
             for i in xrange(max_index):
                 self[i, i + k] = values
         else:
             max_index = min(M, N-k, len(values))
             if max_index <= 0:
                 return
             for i,v in enumerate(values[:max_index]):
                 self[i, i + k] = v

 def __add__(self, other):
     # First check if argument is a scalar
     if isscalarlike(other):
         res_dtype = upcast_scalar(self.dtype, other)
         new = dok_matrix(self.shape, dtype=res_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 largest of
         # the two matrices to be summed.  Would this be a good idea?
         res_dtype = upcast(self.dtype, other.dtype)
         new = dok_matrix(self.shape, dtype=res_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:
         return NotImplemented
     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 __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:
         return NotImplemented
     return new

 def __add__(self, other):
     if isscalarlike(other):
         res_dtype = upcast_scalar(self.dtype, other)
         new = dok_matrix(self.shape, dtype=res_dtype)
         # Add this scalar to every element.
         M, N = self.shape
         for key in itertools.product(xrange(M), xrange(N)):
             aij = dict.get(self, (key), 0) + other
             if aij:
                 new[key] = aij
         # new.dtype.char = self.dtype.char
     elif isspmatrix_dok(other):
         if other.shape != self.shape:
             raise ValueError("Matrix dimensions are not equal.")
         # We could alternatively set the dimensions to the largest of
         # the two matrices to be summed.  Would this be a good idea?
         res_dtype = upcast(self.dtype, other.dtype)
         new = dok_matrix(self.shape, dtype=res_dtype)
         dict.update(new, self)
         with np.errstate(over='ignore'):
             dict.update(new,
                        ((k, new[k] + other[k]) for k in iterkeys(other)))
     elif isspmatrix(other):
         csc = self.tocsc()
         new = csc + other
     elif isdense(other):
         new = self.todense() + other
     else:
         return NotImplemented
     return new

    def extend(self, xi, yi, orders=None):
        """
        Extend the PiecewisePolynomial by a list of points

        Parameters
        ----------
        xi : array_like
            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 : int or list of ints, optional
            A list of polynomial orders, or a single universal order.

        """
        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 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 _get_arrayXarray(self, row, col):
        # inner indexing
        i, j = map(np.atleast_2d, np.broadcast_arrays(row, col))
        newdok = dok_matrix(i.shape, dtype=self.dtype)

        for key in itertools.product(xrange(i.shape[0]), xrange(i.shape[1])):
            v = dict.get(self, (i[key], j[key]), 0)
            if v:
                dict.__setitem__(newdok, key, v)
        return newdok

 def test_correspond_2_and_up(self):
     # Tests correspond(Z, y) on linkage and CDMs over observation sets of
     # different sizes.
     for i in xrange(2, 4):
         y = np.random.rand(i*(i-1)//2)
         Z = linkage(y)
         assert_(correspond(Z, y))
     for i in xrange(4, 15, 3):
         y = np.random.rand(i*(i-1)//2)
         Z = linkage(y)
         assert_(correspond(Z, y))

 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 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)]

        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_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=int)[None,:])
            assert_equal(tri.neighbors, -1 + np.zeros((nd+1), dtype=int)[None,:])

def _printresmat(function, interval, resmat):
    # Print the Romberg result matrix.
    i = j = 0
    print('Romberg integration of', repr(function), end=' ')
    print('from', interval)
    print('')
    print('%6s %9s %9s' % ('Steps', 'StepSize', 'Results'))
    for i in xrange(len(resmat)):
        print('%6d %9f' % (2**i, (interval[1]-interval[0])/(2.**i)), end=' ')
        for j in xrange(i+1):
            print('%9f' % (resmat[i][j]), end=' ')
        print('')
    print('')
    print('The final result is', resmat[i][j], end=' ')
    print('after', 2**(len(resmat)-1)+1, 'function evaluations.')

 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_is_valid_linkage_4_and_up(self):
     # Tests is_valid_linkage(Z) on linkage on observation sets between
     # sizes 4 and 15 (step size 3).
     for i in xrange(4, 15, 3):
         y = np.random.rand(i*(i-1)//2)
         Z = linkage(y)
         assert_(is_valid_linkage(Z) == True)

 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_num_obs_linkage_4_and_up(self):
     # Tests num_obs_linkage(Z) on linkage on observation sets between sizes
     # 4 and 15 (step size 3).
     for i in xrange(4, 15, 3):
         y = np.random.rand(i*(i-1)//2)
         Z = linkage(y)
         assert_equal(num_obs_linkage(Z), i)

 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 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)