Python linalg.solve_banded函数代码示例

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