Python linalg.solveh_banded函数代码示例

def minimize_Z_EL_cython(A, l1, l2, rho):
    """"""
    
    # build banded matrix:
    n = len(A)
    Bmat = numpy.zeros((2,n))
    Bmat[0,:] = -2*l2/rho
    Bmat[1,1:n-1] = 1 +4*l2/rho
    Bmat[1,0] = Bmat[1,n-1] = 1 + 2*l2/rho
    
    # convert A into an array:
    A_ = numpy.zeros((len(A), A[0].shape[0], A[0].shape[0]))
    for i in range(len(A)):
	A_[i,:,:] = A[i]
    
    sudoZ = A_[:]
    for i in range(A[0].shape[0]):
	for j in range(i, A[0].shape[0]):
	    resp = A_[:,i,j]
	    # get LS solution:
	    beta_hat = solveh_banded(Bmat, resp, overwrite_ab=True, overwrite_b=True)
	    
	    # shooting algorithm:
	    beta_hat = Z_shooting.Z_shooting(B=beta_hat, y=resp, l1=l1, l2=l2, tol=0.1, max_iter=100)
	    
	    sudoZ[:,i,j] = beta_hat
	    sudoZ[:,j,i] = beta_hat

    # return to a list (terribly inefficient! I have to change this!)
    Z_ = [None] * len(A)
    for i in range(len(A)):
	Z_[i] = sudoZ[i,:,:]
	
    return Z_

    def fit(self, y, x=None, weights=None, pen=0.):
        banded = True

        if x is None:
            x = self.tau[(self.M-1):-(self.M-1)] # internal knots

        if pen == 0.: # can't use cholesky for singular matrices
            banded = False
            
        if x.shape != y.shape:
            raise ValueError, 'x and y shape do not agree, by default x are the Bspline\'s internal knots'
        
        bt = self.basis(x)
        if pen >= self.penmax:
            pen = self.penmax

        if weights is None:
            weights = N.array(1.)

        wmean = weights.mean()
        _w = N.sqrt(weights / wmean)
        bt *= _w

        # throw out rows with zeros (this happens at boundary points!)

        mask = N.flatnonzero(1 - N.alltrue(N.equal(bt, 0), axis=0))

        bt = bt[:, mask]
        y = y[mask]

        self.df_total = y.shape[0]

        if bt.shape[1] != y.shape[0]:
            raise ValueError, "some x values are outside range of B-spline knots"
        bty = N.dot(bt, _w * y)
        self.N = y.shape[0]
        if not banded:
            self.btb = N.dot(bt, bt.T)
            _g = _band2array(self.g, lower=1, symmetric=True)
            self.coef, _, self.rank = L.lstsq(self.btb + pen*_g, bty)[0:3]
            self.rank = min(self.rank, self.btb.shape[0])
        else:
            self.btb = N.zeros(self.g.shape, N.float64)
            nband, nbasis = self.g.shape
            for i in range(nbasis):
                for k in range(min(nband, nbasis-i)):
                    self.btb[k, i] = (bt[i] * bt[i+k]).sum()

            bty.shape = (1, bty.shape[0])
            self.chol, self.coef = solveh_banded(self.btb + 
                                                 pen*self.g,
                                                 bty, lower=1)

        self.coef = N.squeeze(self.coef)
        self.resid = N.sqrt(wmean) * (y * _w - N.dot(self.coef, bt))
        self.pen = pen

 def test_01_float32(self):
     # Solve
     # [ 4 1 0]     [1]
     # [ 1 4 1] X = [4]
     # [ 0 1 4]     [1]
     #
     ab = array([[-99, 1.0, 1.0], [4.0, 4.0, 4.0]], dtype=float32)
     b = array([1.0, 4.0, 1.0], dtype=float32)
     x = solveh_banded(ab, b)
     assert_array_almost_equal(x, [0.0, 1.0, 0.0])

 def test_check_finite(self):
     # Solve
     # [ 4 1 0]     [1]
     # [ 1 4 1] X = [4]
     # [ 0 1 4]     [1]
     # with the RHS as a 1D array.
     ab = array([[-99, 1.0, 1.0], [4.0, 4.0, 4.0]])
     b = array([1.0, 4.0, 1.0])
     x = solveh_banded(ab, b, check_finite=False)
     assert_array_almost_equal(x, [0.0, 1.0, 0.0])

 def test_01_complex(self):
     # Solve
     # [ 4 -j 0]     [ -j]
     # [ j 4 -j] X = [4-j]
     # [ 0 j  4]     [4+j]
     #
     ab = array([[-99, -1.0j, -1.0j], [4.0, 4.0, 4.0]])
     b = array([-1.0j, 4.0-1j, 4+1j])
     x = solveh_banded(ab, b)
     assert_array_almost_equal(x, [0.0, 1.0, 1.0])

 def test_01_upper(self):
     # Solve
     # [ 4 1 0]     [1]
     # [ 1 4 1] X = [4]
     # [ 0 1 4]     [1]
     # with the RHS as a 1D array.
     ab = array([[-99, 1.0, 1.0], [4.0, 4.0, 4.0]])
     b = array([1.0, 4.0, 1.0])
     x = solveh_banded(ab, b)
     assert_array_almost_equal(x, [0.0, 1.0, 0.0])

 def test_tridiag_01_lower(self):
     # Solve
     # [ 4 1 0]     [1]
     # [ 1 4 1] X = [4]
     # [ 0 1 4]     [1]
     #
     ab = array([[4.0, 4.0, 4.0], [1.0, 1.0, -99]])
     b = array([1.0, 4.0, 1.0])
     x = solveh_banded(ab, b, lower=True)
     assert_array_almost_equal(x, [0.0, 1.0, 0.0])

 def test_03_upper(self):
     # Solve
     # [ 4 1 0]     [1]
     # [ 1 4 1] X = [4]
     # [ 0 1 4]     [1]
     # with the RHS as a 2D array with shape (3,1).
     ab = array([[-99, 1.0, 1.0], [4.0, 4.0, 4.0]])
     b = array([1.0, 4.0, 1.0]).reshape(-1,1)
     x = solveh_banded(ab, b)
     assert_array_almost_equal(x, array([0.0, 1.0, 0.0]).reshape(-1,1))

 def timestep(self, input_signal, output_condition=0,
              curr_noise=0, sub_noise=0, cond_noise=0):
     """
     Compute the time evolution for one timestep using increments of
     the driving noise (current, subunit or conductance or a mixture)
     and the corresponding input signal at the soma. The output as an
     outflux current can also be specified.
     
     Note that the noise has to be initialized outside the class in
     order to be able to compare models for given trajectories of the
     noise. Shape of the noise has to be
     curr_noise.shape = (N_axon+1, M)
     sub_noise.shape = (3, N_axon+1, M)
     cond_noise.shape = (11, N_axon+1, M)
     """
     [v, m, h, n] = self.state.copy()
     # Explicit Euler(-Maruyama) step for gating variables
     self.state[1] += self.dt * self.t_m(v) * (self.m_inf(v)-m)
     self.state[2] += self.dt * self.t_h(v) * (self.h_inf(v)-h)
     self.state[3] += self.dt * self.t_n(v) * (self.n_inf(v)-n)
     # If subunit noise is switched on
     if np.max(self.sigma_sub):
         N = self.N_axon + 1
         self.state[1][:N, :] += (np.sqrt(self.alpha_m(v[:N, :]) *
                 (1-m[:N, :]) + self.beta_m(v[:N, :]) * m[:N, :]) *
                 self.sigma_sub[0] * sub_noise[0])
         self.state[2][:N, :] += (np.sqrt(self.alpha_h(v[:N, :]) *
                 (1-h[:N, :]) + self.beta_h(v[:N, :]) * h[:N, :]) *
                 self.sigma_sub[1] * sub_noise[1])
         self.state[3][:N, :] += (np.sqrt(self.alpha_n(v[:N, :]) *
                 (1-n[:N, :]) + self.beta_n(v[:N, :]) * n[:N, :]) *
                 self.sigma_sub[2] * sub_noise[2])
     # Respect boundaries!
     for i in xrange(3):
         self.state[i][self.state[i]>1] = 1
         self.state[i][self.state[i]<0] = 0
     # Semi-implicit Euler step for voltage variable
     rhs = v - self.dt*self.fv(v, m, h, n)
     rhs[:self.N_axon+1, :] += self.sigma_curr*curr_noise
     # If conductance noise is switched on
     if self.sigma_cond:
         rhs[:self.N_axon+1,:] -= self.compute_cond_noise(cond_noise)
     # Input signal as Neumann boundary condition for left axon endpoint
     rhs[0] += self.boundary * input_signal
     # Output condition as Neumann boundary condition for right axon endpoint
     rhs[-1] += self.boundary * output_condition
     # Modification of RHS to make discrete Neumann Laplacian symmetric
     rhs[0] *= 0.5
     rhs[-1] *= 0.5
     # Solve the linear problem
     self.state[0] = lng.solveh_banded(self.bandedmatrix, rhs,
                                       check_finite=False)
     # Update state
     #self.state = [v_neu, m_neu, h_neu, n_neu]
     self.time_elapsed += self.dt

 def test_01_complex(self):
     # Solve
     # [ 4 -j  2  0]     [2-j]
     # [ j  4 -j  2] X = [4-j]
     # [ 2  j  4 -j]     [4+j]
     # [ 0  2  j  4]     [2+j]
     #
     ab = array([[0.0, 0.0, 2.0, 2.0], [-99, -1.0j, -1.0j, -1.0j], [4.0, 4.0, 4.0, 4.0]])
     b = array([2 - 1.0j, 4.0 - 1j, 4 + 1j, 2 + 1j])
     x = solveh_banded(ab, b)
     assert_array_almost_equal(x, [0.0, 1.0, 1.0, 0.0])

 def test_tridiag_02_lower(self):
     # Solve
     # [ 4 1 0]     [1 4]
     # [ 1 4 1] X = [4 2]
     # [ 0 1 4]     [1 4]
     #
     ab = array([[4.0, 4.0, 4.0], [1.0, 1.0, -99]])
     b = array([[1.0, 4.0], [4.0, 2.0], [1.0, 4.0]])
     x = solveh_banded(ab, b, lower=True)
     expected = array([[0.0, 1.0], [1.0, 0.0], [0.0, 1.0]])
     assert_array_almost_equal(x, expected)

 def test_01_lower(self):
     # Solve
     # [ 4 1 2 0]     [1]
     # [ 1 4 1 2] X = [4]
     # [ 2 1 4 1]     [1]
     # [ 0 2 1 4]     [2]
     #
     ab = array([[4.0, 4.0, 4.0, 4.0], [1.0, 1.0, 1.0, -99], [2.0, 2.0, 0.0, 0.0]])
     b = array([1.0, 4.0, 1.0, 2.0])
     x = solveh_banded(ab, b, lower=True)
     assert_array_almost_equal(x, [0.0, 1.0, 0.0, 0.0])

 def test_01_float32(self):
     warnings.simplefilter("ignore", category=DeprecationWarning)
     # Solve
     # [ 4 1 0]     [1]
     # [ 1 4 1] X = [4]
     # [ 0 1 4]     [1]
     #
     ab = array([[-99, 1.0, 1.0], [4.0, 4.0, 4.0]], dtype=float32)
     b = array([1.0, 4.0, 1.0], dtype=float32)
     c, x = solveh_banded(ab, b)
     assert_array_almost_equal(x, [0.0, 1.0, 0.0])

 def test_01_complex(self):
     warnings.simplefilter("ignore", category=DeprecationWarning)
     # Solve
     # [ 4 -j 0]     [ -j]
     # [ j 4 -j] X = [4-j]
     # [ 0 j  4]     [4+j]
     #
     ab = array([[-99, -1.0j, -1.0j], [4.0, 4.0, 4.0]])
     b = array([-1.0j, 4.0 - 1j, 4 + 1j])
     c, x = solveh_banded(ab, b)
     assert_array_almost_equal(x, [0.0, 1.0, 1.0])

 def test_03_upper(self):
     warnings.simplefilter("ignore", category=DeprecationWarning)
     # Solve
     # [ 4 1 0]     [1]
     # [ 1 4 1] X = [4]
     # [ 0 1 4]     [1]
     # with the RHS as a 2D array with shape (3,1).
     ab = array([[-99, 1.0, 1.0], [4.0, 4.0, 4.0]])
     b = array([1.0, 4.0, 1.0]).reshape(-1, 1)
     c, x = solveh_banded(ab, b)
     assert_array_almost_equal(x, array([0.0, 1.0, 0.0]).reshape(-1, 1))