def interpolation(interP,knotP=None):
    """
        Interpolates the given points and returns an object of the Spline class 
            Arguments:
                interP: interpolation points, (L x 2) matrix
                knotP: knot points, (L+4 x 1) matrix
                    default: equidistant on [0,1]
    """
    nip=len(interP)
    ctrlP=zeros((nip,2))
    knotP=fixKnotPoints(nip, knotP)
    xi=(knotP[:-2]+knotP[1:-1]+knotP[2:])/3
    nMatrix=zeros((nip,nip))
    for i in xrange(nip):
        fun=basisFunction(i,knotP)
        for k in xrange(nip):
            nMatrix[k,i]=fun(xi[k],3)
    print nMatrix
    print knotP

    ctrlP[:,0]=sl.solve_banded((nip-4,nip-4),nMatrix,interP[:,0])
    ctrlP[:,1]=sl.solve_banded((nip-4,nip-4),nMatrix,interP[:,1])

    
    return Spline(ctrlP,knotP)

	def calculate_vd_vL_vR(self):
		for index in xrange(self.neuron.BPcount) :
			if self.test_list[index] :
				i = self.i_list[index]
				j = self.j_list[index]
				self.vL[i:j+1] = solve_banded((1,1),self.ab[:,i:j+1],self.bL[i:j+1],overwrite_ab=False,overwrite_b=False)
				self.vR[i:j+1] = solve_banded((1,1),self.ab[:,i:j+1],self.bR[i:j+1],overwrite_ab=False,overwrite_b=False)
				self.d[i:j+1] = solve_banded((1,1),self.ab[:,i:j+1],self.bd[i:j+1],overwrite_ab=False,overwrite_b=False)

 def solve(self,f): # Solves u'=f for u(-1)=0
     fhat = ifcgt(f*(1-self.x)) # Compute (psi_j,f)
     fhat[0] *= 2
     y = solve_banded((1,0),self.Lbands,fhat,\
         overwrite_ab=False,overwrite_b=True)
     uhat = solve_banded((0,1),self.Ubands,y,\
         overwrite_ab=False,overwrite_b=True)
     u = fcgt(uhat)*(1+self.x) # u=phi*uhat
     return u

    def fit_grad_descent(self,r_init=0.001,learning_rate=0.1):
        ''' Calculate fit using gradient decent, as opposed to iteratively
        computing exact solutions
        '''
        delta_r = 1
        r = r_init

        r_list = []
        error_list = []

        F_C = solve_banded((1,1),self.ab,self.F_M - r*self.F_N)

        itmax = 10000
        it = 0

        exceed_bounds = False
        while (delta_r > 0.0001 and it < itmax):

            F_C = solve_banded((1,1),self.ab,self.F_M - r*self.F_N)
            r_new = r + learning_rate*np.mean((self.F_M - F_C - r*self.F_N)*self.F_N)

            '''compute error on cross-validation set'''
            F_C_crossval = solve_banded((1,1),self.ab,self.F_M_crossval - r*self.F_N_crossval)
#            error_it = error_calc_outlier(self.F_M_crossval, self.F_N_crossval, F_C_crossval, r)
            error_it = abs(error_calc(self.F_M_crossval, self.F_N_crossval, F_C_crossval, r))
            
            delta_r = np.abs(r_new - r)/r
            r = r_new

            r_list.append(r)
            error_list.append(error_it)
            it+=1
            
            if error_it > error_list[it-1]: # early stopping
                break

        # if r or error_it go out of acceptable bounds, break 
        if r < 0.0 or r > 1.0:
            exceed_bounds = True
        if error_it > 0.2:
            exceed_bounds = True

        F_C = solve_banded((1,1),self.ab,self.F_M - r*self.F_N)

        r_list = np.array(r_list)
        error_list = np.array(error_list)
        self.r_vals = r_list
        self.error_vals = error_list
        self.r = self.r_vals[-1]
        self.F_C = F_C
        self.F_C_crossval = F_C_crossval
        self.it = it

        return exceed_bounds

def impl(V, L1, R1x, L2, R2x, dt, n, initial=None, callback=None):
    V = V.copy()

    # L1i = flatten_tensor(L1)
    L1i = L1.copy()
    R1 = np.array(R1x).T

    # L2i = flatten_tensor(L2)
    L2i = L2.copy()
    R2 = np.array(R2x)

    m = 2

    # L  = (As + Ass - H.interest_rate*np.eye(nspots))*-dt + np.eye(nspots)
    L1i.data *= -dt
    L1i.data[m, :] += 1
    R1 *= -dt

    L2i.data *= -dt
    L2i.data[m, :] += 1
    R2 *= -dt

    offsets1 = (abs(min(L1i.offsets)), abs(max(L1i.offsets)))
    offsets2 = (abs(min(L2i.offsets)), abs(max(L2i.offsets)))

    dx = np.gradient(spots)[:, np.newaxis]
    dy = np.gradient(vars)
    Y, X = np.meshgrid(vars, spots)
    gradgrid = dt * coeffs[(0, 1)](0, X, Y) / (dx * dy)
    gradgrid[:, 0] = 0
    gradgrid[:, -1] = 0
    gradgrid[0, :] = 0
    gradgrid[-1, :] = 0

    print_step = max(1, int(n / 10))
    to_percent = 100.0 / n
    utils.tic("Impl:\t")
    for k in xrange(n):
        if not k % print_step:
            if np.isnan(V).any():
                print "Impl fail @ t = %f (%i steps)" % (dt * k, k)
                return crumbs
            print int(k * to_percent),
        if callback is not None:
            callback(V, ((n - k) * dt))
        Vsv = np.gradient(np.gradient(V)[0])[1] * gradgrid
        # Vsv = 0.0
        V = spl.solve_banded(offsets2, L2i.data, (V + Vsv - R2).flat, overwrite_b=True).reshape(V.shape)
        V = spl.solve_banded(offsets1, L1i.data, (V - R1).T.flat, overwrite_b=True).reshape(V.shape[::-1]).T
    crumbs.append(V.copy())
    utils.toc(":  \t")
    return crumbs

def crank(V, L1, R1x, L2, R2x, dt, n, crumbs=[], callback=None):
    V = V.copy()
    dt *= 0.5

    L1e = flatten_tensor(L1)
    L1i = L1e.copy()
    R1 = np.array(R1x).T

    L2e = flatten_tensor(L2)
    L2i = L2e.copy()
    R2 = np.array(R2x)

    m = 2

    # L  = (As + Ass - r*np.eye(nspots))*-dt + np.eye(nspots)
    L1e.data *= dt
    L1e.data[m, :] += 1
    L1i.data *= -dt
    L1i.data[m, :] += 1
    R1 *= dt

    L2e.data *= dt
    L2e.data[m, :] += 1
    L2i.data *= -dt
    L2i.data[m, :] += 1
    R2 *= dt

    offsets1 = (abs(min(L1i.offsets)), abs(max(L1i.offsets)))
    offsets2 = (abs(min(L2i.offsets)), abs(max(L2i.offsets)))

    print_step = max(1, int(n / 10))
    to_percent = 100.0 / n
    utils.tic("Crank:")
    R = R1 + R2
    normal_shape = V.shape
    transposed_shape = normal_shape[::-1]
    for k in xrange(n):
        if not k % print_step:
            if isnan(V).any():
                print "Crank fail @ t = %f (%i steps)" % (dt * k, k)
                return crumbs
            print int(k * to_percent),
        if callback is not None:
            callback(V, ((n - k) * dt))
        V = (L2e.dot(V.flat).reshape(normal_shape) + R).T
        V = spl.solve_banded(offsets1, L1i.data, V.flat, overwrite_b=True)
        V = (L1e.dot(V).reshape(transposed_shape).T) + R
        V = spl.solve_banded(offsets2, L2i.data, V.flat, overwrite_b=True).reshape(normal_shape)
        crumbs.append(V.copy())
    utils.toc()
    return crumbs

 def test_check_finite(self):
     a = array([[1.0, 20, 0, 0], [-30, 4, 6, 0], [2, 1, 20, 2], [0, -1, 7, 14]])
     ab = array([[0.0, 20, 6, 2], [1, 4, 20, 14], [-30, 1, 7, 0], [2, -1, 0, 0]])
     l, u = 2, 1
     b4 = array([10.0, 0.0, 2.0, 14.0])
     x = solve_banded((l, u), ab, b4, check_finite=False)
     assert_array_almost_equal(dot(a, x), b4)

    def _solve_for_angular_rates(self, dt, angular_rates, rotvecs):
        angular_rate_first = angular_rates[0].copy()

        A = _angular_rate_to_rotvec_dot_matrix(rotvecs)
        A_inv = _rotvec_dot_to_angular_rate_matrix(rotvecs)
        M = _create_block_3_diagonal_matrix(
            2 * A_inv[1:-1] / dt[1:-1, None, None],
            2 * A[1:-1] / dt[1:-1, None, None],
            4 * (1 / dt[:-1] + 1 / dt[1:]))

        b0 = 6 * (rotvecs[:-1] * dt[:-1, None] ** -2 +
                  rotvecs[1:] * dt[1:, None] ** -2)
        b0[0] -= 2 / dt[0] * A_inv[0].dot(angular_rate_first)
        b0[-1] -= 2 / dt[-1] * A[-1].dot(angular_rates[-1])

        for iteration in range(self.MAX_ITER):
            rotvecs_dot = _matrix_vector_product_of_stacks(A, angular_rates)
            delta_beta = _angular_acceleration_nonlinear_term(
                rotvecs[:-1], rotvecs_dot[:-1])
            b = b0 - delta_beta
            angular_rates_new = solve_banded((5, 5), M, b.ravel())
            angular_rates_new = angular_rates_new.reshape((-1, 3))

            delta = np.abs(angular_rates_new - angular_rates[:-1])
            angular_rates[:-1] = angular_rates_new
            if np.all(delta < self.TOL * (1 + np.abs(angular_rates_new))):
                break

        rotvecs_dot = _matrix_vector_product_of_stacks(A, angular_rates)
        angular_rates = np.vstack((angular_rate_first, angular_rates[:-1]))

        return angular_rates, rotvecs_dot

	def solve_banded(self, rhs):
		rhsex = np.zeros(self.M)
		rhsex[::self.nBlockSize] = rhs
		if not hasattr(self, 'Aex_diag'):
			Aex_diag = sp.dia_matrix(self.Aex)
		solex = solve_banded((self.m + 1,self.m + 1), Aex_diag.data[::-1,:], rhsex)
		return solex[::self.nBlockSize]

def beamwarming(u, nt, dt, dx):
    ##Tridiagonal setup##
    a = numpy.zeros_like(u)
    b = numpy.ones_like(u)
    c = numpy.zeros_like(u)
    d = numpy.zeros_like(u)

    un = numpy.zeros((nt,len(u)))
    un[:]=u.copy()

    for n in range(1,nt):
        u[0] = 1
        E = utoE(u)
        au = utoA(u)

        a[0] = -dt/(4*dx)*u[0]
        a[1:] = -dt/(4*dx)*u[0:-1]
        a[-1] = -dt/(4*dx)*u[-1]

        #b is all ones

        c[:-1] = dt/(4*dx)*u[1:]

        print str(-.5*dt/dx*(E[2:]-E[0:-2])+dt/(4*dx)*(au[2:]-au[:-2]))
        #print str(dt/(4*dx)*(au[2:]-au[:-2]))
        d[1:-1] = u[1:-1]-.5*dt/dx*(E[2:]-E[0:-2])+dt/(4*dx)*(au[2:]-au[:-2])

        ###subtract a[0]*LHS B.C to 'fix' thomas algorithm
        d[0] = u[0] - .5*dt/dx*(E[1]-E[0])+dt/(4*dx)*(au[1]-au[0]) - a[0]

        ab = numpy.matrix([c,b,a])
        u = linalg.solve_banded((1,1), ab, d)
        u[0]=1
        un[n] = u.copy()
    return un

def solve(fk):
    
    N = len(fk)+2
    k = ST.wavenumbers(N)
    if solver == "banded":
        A = np.zeros((N-2, N-2))
        A[-1, :] = -2*np.pi*(k+1)*(k+2)
        for i in range(2, N-2, 2):
            A[-i-1, i:] = -4*np.pi*(k[:-i]+1)
        uk_hat = solve_banded((0, N-3), A, fk)
        
    elif solver == "sparse":
        aij = [-2*np.pi*(k+1)*(k+2)]
        for i in range(2, N-2, 2):
            aij.append(np.array(-4*np.pi*(k[:-i]+1)))    
        A = diags(aij, range(0, N-2, 2))
        uk_hat = la.spsolve(A, fk)
        
    elif solver == "bs":
        fc = fk.copy()
        uk_hat = np.zeros(N-2)        
        uk_hat = SFTc.BackSubstitution_1D(uk_hat, fc)
                
        #for i in range(N-3, -1, -1):
            #for l in range(i+2, N-2, 2):
                #fc[i] += (4*np.pi*(i+1)uk_hat[l])
            #uk_hat[i] = -fc[i] / (2*np.pi*(i+1)*(i+2))
            
            
    return uk_hat

    def _forward_control_points(self, dps):
        num_cps = dps.shape[0] - 1;

        ab = np.matrix([
          np.tile(1, num_cps),
          np.tile(4, num_cps),
          np.tile(1, num_cps)
        ])
        
        # Fixup first and last row
        ab[0,0] = 0
        ab[1,0] = 2
        ab[1, num_cps - 1] = 7
        ab[2, num_cps - 1] = 0
        
        b = np.empty([num_cps, dps.shape[1]])
        
        for i in range(0, num_cps - 1):
            b[i] = 4 * dps[i] + 2 * dps[i + 1]
        
        # Fixup first and last element
        b[0] = dps[0] + 2 * dps[1]
        b[num_cps - 1] = 8 * dps[num_cps - 1] + dps[num_cps]
        
        return solve_banded((1, 1), ab, b)

    def fit_block_coordinate_desc(self, r_init=5.0, min_delta_r=0.00000001):
        F_M = np.concatenate(self.F_M)
        F_N = np.concatenate(self.F_N)

        r_vals = []
        error_vals = []
        r = r_init

        delta_r = None
        it = 0

        ab = ab_from_T(self.T, self.lam, self.dt)
        while delta_r is None or delta_r > min_delta_r:
            F_C = solve_banded((1, 1), ab, F_M - r * F_N)
            new_r = - np.sum((F_C - F_M) * F_N) / np.sum(np.square(F_N))
            error = self.estimate_error(new_r)

            error_vals.append(error)
            r_vals.append(new_r)

            if r is not None:
                delta_r = np.abs(r - new_r) / r

            r = new_r
            it += 1

        self.r_vals = r_vals
        self.error_vals = error_vals
        self.r = r_vals[-1]
        self.error = error_vals.min()

def diffuseImplicit(gr, phi, k, dt):
    """ diffuse phi implicitly through timestep dt """

    phinew = gr.scratchArray()
    
    alpha = k*dt/gr.dx**2

    # create the RHS of the matrix
    R = phi[gr.ilo:gr.ihi+1] 
    
    # create the diagonal, d+1 and d-1 parts of the matrix
    d = (1.0 + 2.0*alpha)*numpy.ones(gr.nx)
    u = -alpha*numpy.ones(gr.nx)
    u[0] = 0.0

    l = -alpha*numpy.ones(gr.nx)
    l[gr.nx-1] = 0.0

    # set the boundary conditions by changing the matrix elements

    # homogeneous neumann
    d[0] = 1.0 + alpha
    d[gr.nx-1] = 1.0 + alpha


    # solve
    A = numpy.matrix([u,d,l])
    phinew[gr.ilo:gr.ihi+1] = linalg.solve_banded((1,1), A, R)

    return phinew

 def __init__(self, pixbound, flux):
     
     #- If pixbound and flux have same number of elements, assume they
     #- are centers rather than edges
     if len(pixbound) == len(flux):
         pixbound = cen2bound(pixbound)
     
     npix = len(flux)
     # Test for correct argument dimensions:
     if (len(pixbound) - npix) != 1:
         raise PixSplineError('Need one more element in pixbound than in flux!')
     # The array of "delta-x" values:
     dxpix = pixbound[1:] - pixbound[:-1]
     # Test for monotonic increase:
     if dxpix.min() <= 0.:
         raise PixSplineError('Pixel boundaries not monotonically increasing!')
     self.npix = npix
     self.pixbound = pixbound.copy()
     self.dxpix = dxpix.copy()
     self.xcen = 0.5 * (pixbound[1:] + pixbound[:-1]).copy()
     self.flux = flux.copy()
     maindiag = (dxpix[:-1] + dxpix[1:]) / 3.
     offdiag = dxpix[1:-1] / 6.
     upperdiag = np.append(0., offdiag)
     lowerdiag = np.append(offdiag, 0.)
     band_matrix = np.vstack((upperdiag, maindiag, lowerdiag))
     # The right-hand side:
     rhs = flux[1:] - flux[:-1]
     # Solve the banded matrix for the slopes at the ducks:
     acoeff = la.solve_banded((1,1), band_matrix, rhs)
     self.duckslopes = np.append(np.append(0., acoeff), 0.)

def _get_natural_f(knots):
    """Returns mapping of natural cubic spline values to 2nd derivatives.

    .. note:: See 'Generalized Additive Models', Simon N. Wood, 2006, pp 145-146

    :param knots: The 1-d array knots used for cubic spline parametrization,
     must be sorted in ascending order.
    :return: A 2-d array mapping natural cubic spline values at
     knots to second derivatives.

    :raise ImportError: if scipy is not found, required for
     ``linalg.solve_banded()``
    """
    try:
        from scipy import linalg
    except ImportError: # pragma: no cover
        raise ImportError("Cubic spline functionality requires scipy.")

    h = knots[1:] - knots[:-1]
    diag = (h[:-1] + h[1:]) / 3.
    ul_diag = h[1:-1] / 6.
    banded_b = np.array([np.r_[0., ul_diag], diag, np.r_[ul_diag, 0.]])
    d = np.zeros((knots.size - 2, knots.size))
    for i in range(knots.size - 2):
        d[i, i] = 1. / h[i]
        d[i, i + 2] = 1. / h[i + 1]
        d[i, i + 1] = - d[i, i] - d[i, i + 2]

    fm = linalg.solve_banded((1, 1), banded_b, d)

    return np.vstack([np.zeros(knots.size), fm, np.zeros(knots.size)])

    def __init__(self, a, b, y, alpha=0, beta=0):
        y = np.asarray(y)
        n = y.shape[0] - 1
        h = (b - a)/n

        coeff = np.zeros(n + 3, dtype=y.dtype)
        # Solutions to boundary coeffcients of spline
        coeff[1] = 1/6. * (y[0] - (alpha * h**2)/6) #C2 in paper
        coeff[n + 1] = 1/6. * (y[n] - (beta * h**2)/6) #cn+2 in paper

        # Compressed tridiagonal matrix 
        ab = np.ones((3, n - 1), dtype=float)
        ab[0,0] = 0 # Because top row is upper diag with one less elem
        ab[1, :] = 4
        ab[-1,-1] = 0 # Because bottom row is lower diag with one less elem
        
        B = y[1:-1].copy() #grabs elements y[1] - > y[n-2] for reduced array
        B[0] -= coeff[1]
        B[-1] -=  coeff[n + 1]

        coeff[2:-2] = la.solve_banded((1, 1), ab, B, overwrite_ab=True, 
                        overwrite_b=True, check_finite=False)

        coeff[0] = alpha * h**2/6. + 2 * coeff[1] - coeff[2]
        coeff[-1] = beta * h**2/6. + 2 * coeff[-2] - coeff[-3]

        self.a = a          # Lower-bound of domain
        self.b = b          # Uppser-bound of domain
        self.coeffs = coeff # Spline coefficients
        self.is_complex = (y.dtype == complex) #Tells which dtype solver to use

def phi_np1(phi_n, x, dx, t, dt, M, eps, a, Df, verbose=False):
    """ returns phi(x, t=n+1) given phi(x, t=n) using finite difference method.

    Uses *backward* finite difference in time (implicit), central finite
    difference in x
    """

    mu = M * dt / (dx * dx)
    eps2 = eps * eps

    # Set up matrix from 1 to n-2, and handle x = 0 and x = n-1 via boundary
    # conditions
    alpha, beta = boundary(t)
    boundary_vec = zeros_like(phi_n[1:-1])
    boundary_vec[0] = -mu * eps2 * alpha
    boundary_vec[-1] = -mu * eps2 * beta

    fp = fprime(phi_n[1:-1], a, Df)

    b = phi_n[1:-1] - boundary_vec + M * dt * fp

    d = ones_like(phi_n[1:-1]) * (1 + 2 * mu * eps2)
    u = ones_like(phi_n[1:-1]) * (-mu * eps2)
    u[0] = 0
    l = ones_like(phi_n[1:-1]) * (-mu * eps2)
    l[-1] = 0

    ab = np.matrix([u, d, l])

    res = zeros_like(phi_n)
    res[1:-1] = linalg.solve_banded((1, 1), ab, b)
    res[0] = alpha
    res[-1] = beta
    return res

    def __init__(self, xp, yp):

        if np.any(np.diff(xp) < 0):
            raise TypeError("xp must be in ascending order")

        # n = len(x)
        self.m = len(xp)

        xk = xp[1:-1]
        yk = yp[1:-1]
        xkp = xp[2:]
        ykp = yp[2:]
        xkm = xp[:-2]
        ykm = yp[:-2]

        b = (ykp - yk) / (xkp - xk) - (yk - ykm) / (xk - xkm)
        l = (xk - xkm) / 6.0
        d = (xkp - xkm) / 3.0
        u = (xkp - xk) / 6.0
        # u[0] = 0.0  # non-existent entries
        # l[-1] = 0.0

        # solve for second derivatives
        fpp = solve_banded((1, 1), np.matrix([u, d, l]), b)
        self.fpp = np.concatenate([[0.0], fpp, [0.0]])  # natural spline
        self.xp = xp
        self.yp = yp

def radau(alpha,beta,xr):
    """
        Compute the Radau nodes and weights with the preassigned node xr
        
        Inputs: 
        alpha - recursion coefficients
        beta - recursion coefficients
        xr - assigned node location

        Outputs: 
        x - quadrature nodes		
        w - quadrature weights
   
        Based on the section 7 of the paper "Some modified matrix eigenvalue
        problems" by Gene Golub, SIAM Review Vol 15, No. 2, April 1973, pp.318--334
    """
    from scipy.linalg import solve_banded
    n = len(alpha)-1
    f = np.zeros(n)
    f[-1] = beta[-1]
    A = np.vstack((np.sqrt(beta),alpha-xr))
    J = np.vstack((A[:,0:-1],A[0,1:]))
    delta = solve_banded((1,1),J,f)
    alphar = alpha
    alphar[-1] = xr+delta[-1]  
    x,w = gauss(alphar,beta)
    return x,w