Python numpy.tril函数代码示例

def _chol_blocked_fwd(L, Adot, NB=256, inplace=False):
    """
    Forwards-mode differentiation through the Cholesky decomposition
    
    Obtain L_dot from Sigma_dot, where "_dot" means sensitivities in
    forwards-mode differentiation, and Sigma = L @ L.T.

    This version uses a blocked algorithm to update sensitivities Adot
    in place. tril(Adot) should start containing Sigma_dot, and will
    end containing the L_dot. Take tril() of the answer if
    triu(Adot,1) did not start out filled with zeros. Unlike the
    unblocked routine, if the upper triangular part of Adot started
    with non-zero values, some of these will be overwritten.

    If inplace=False, a copy of Adot is modified instead of the
    original. The Abar that was modified is returned.
    """
    if not inplace:
        Adot = Adot.copy()
    for j in range(0, L.shape[0], NB):
        k = min(N, j + NB)
        R, D, B, C = _level3partition(L, j, k)
        Rdot, Ddot, Bdot, Cdot = _level3partition(Adot, j, k)
        Ddot[:] = tril(Ddot) - tril(np.dot(Rdot, R.T) + np.dot(R, Rdot.T))
        #chol_unblocked_fwd(D, Ddot, inplace=True) # slow in Python
        Ddot[:] = _chol_symbolic_fwd(D, Ddot + tril(Ddot, -1).T)
        Cdot -= (np.dot(Bdot, R.T) + np.dot(B, Rdot.T))
        #Cdot[:] = (Cdot - [email protected]) @ inv(tril(D)).T
        Cdot[:] = _st(D, Cdot.T - np.dot(Ddot, C.T)).T
    return Adot

def plot_distance_matrix(distance_csv_file, outfile="similarity_matrix_plot.pdf"):
    """
    plotting distance matrix between organisms 

    the distance between the organisms are calculated based on the difference in their sequence composition 
    """

    distance = pandas.read_csv(distance_csv_file, header=0)
    C = numpy.tril(distance)
    sim = 1-distance
    C = numpy.tril(sim)
    N = sim.shape[1]
    C = numpy.ma.masked_array(C, C == 0)

    A = numpy.array([(y, x) for x in range(N, -1, -1) for y in range(N + 1)])
    t = numpy.array([[0.5, 1], [0.5, -1]])
    A = numpy.dot(A, t)
    X = A[:, 1].reshape(N + 1, N + 1)
    Y = A[:, 0].reshape(N + 1, N + 1)
    fig = pylab.figure(figsize=(20,20))
    ax = fig.add_subplot(121, frame_on=False, aspect=2.0)
    ax.set_xticks([])
    ax.set_yticks([])
    caxes = pylab.pcolormesh(X, Y, np.flipud(C), axes=ax)
    ax.set_xlim(right=0)
    fig.savefig(outfile, bbox_inches='tight')

def cmat_for_key(connectome, key, number_of_nodes=None,
                 force_symmetric=True):
    """Return a N x N connection matrix for given connectome and
    key. The connection matrix is returned as a numpy ndarray.

    """

    # create our new shiny connection matrix
    import numpy
    if number_of_nodes is None:
        n = max(connectome.nodes())
    else:
        n = number_of_nodes
    new_cmat = numpy.zeros((n,n))

    # extract the value for key for every edge in the given connectome
    for i,j in connectome.edges_iter():
        new_cmat[i-1][j-1] = connectome[i][j][key]
        
    # do we need to do anything regarding symmetry?
    if force_symmetric and (new_cmat - new_cmat.T != 0).any():
        #...if one-sided (no information below diagonal)
        if (numpy.tril(new_cmat,-1) == 0).all():
            # project above diagonal onto below diagonal
            new_cmat += numpy.tril(new_cmat.T, -1)
        #...else, we will assume two-sided unequal
        else:
            # our solution will be to take the mean of each pair of
            # reflected indices
            new_cmat = (new_cmat + new_cmat.T ) / 2.0

    # return the cmat
    return new_cmat

 def test_al_mohy_higham_2012_experiment_1(self):
     # Matrix square root of a tricky upper triangular matrix.
     A = _get_al_mohy_higham_2012_experiment_1()
     A_sqrtm, info = sqrtm(A, disp=False)
     A_round_trip = A_sqrtm.dot(A_sqrtm)
     assert_allclose(A_round_trip, A, rtol=1e-5)
     assert_allclose(np.tril(A_round_trip), np.tril(A))

 def autocor_two_time(self,print_=False,save_=True,filename=None):
     global buf,num,cts,cur,g12, countl        
     global Ndel,Npix
     global time_ind  #generate a time-frame for each level
     global g12x, g12y, g12z #for interpolate
     start_time = time.time()
     buf=zeros([nolev,nobuf,nopixels])  #// matrix of buffers, for store img
     cts=zeros(nolev)
     cur=ones(nolev) * nobuf
     countl = array(zeros(  nolev ),dtype='int')        
     g12 =  zeros( [ noframes,noframes, noqs] ) 
     g12x=[]
     g12y=[]
     g12z=[]        
     num= array(zeros(  nolev ),dtype='int')        
     time_ind ={key: [] for key in range(nolev)}         
     ttx=0        
     for n in range(1,noframes +1 ):   ##do the work here
         self.insertimg_twotime(begframe+n-1, print_=print_)
         if  n %(noframes/10) ==0:
             sys.stdout.write("#")
             sys.stdout.flush()                
     for q in range(noqs):            
         x0 =  g12[:,:,q]
         g12[:,:,q] = tril(x0) +  tril(x0).T - diag(diag(x0))            
     elapsed_time = time.time() - start_time
     print 'Total time: %.2f min' %(elapsed_time/60.)
     if save_:
         if filename==None:
             filename =  'g12_-%s-%s_ImgReadMethod_'%(
         begframe,begframe+noframes-1)+FOUT
             
         save(  RES_DIR + filename+FOUT, g12)
         print 'the %s was stored in %s'%(filename,RES_DIR)
     return g12, (elapsed_time/60.)

def q150():
    x = numpy.arange(2**20, dtype=numpy.int64)
    x = (615949*x + 797807) % 2**20
    s = numpy.empty(500500+1,dtype=numpy.int64)
    t = 0
    s[0] = 0
    for k in xrange(1,500500+1):
        t = x[t]
        s[k] = t - 2**19
    del x
    print s[1:4]
    r,c = numpy.mgrid[:1000,:1000]
    s = s[numpy.tril(r*(r+1)/2 + c + 1)]
    del r,c
    s0 = s
    best = s.min()
    t = s
    p = numpy.zeros((1001,1001), dtype=numpy.int64)
    for i in xrange(1,1000):
        n = s[:-1,:-1] + numpy.tril(t[1:,:-1]) + t[1:,1:] - numpy.tril(p[2:,1:-1])
        #n = s[:-1,:-1] +t[1:,:-1] + t[1:,1:] - p[2:,1:-1]
        print i,best,t[0,0]
        p,t,s = t,n,s[:-1,:-1]
        best = min(best, t.min())
    print best

def _chol_blocked_rev(L, Abar, NB=256, inplace=False):
    """
    Reverse-mode differentiation through the Cholesky decomposition
    
    Obtain tril(Sigma_bar) from L_bar, where "_bar" means sensitivities
    in reverse-mode differentiation, and Sigma = L @ L.T.

    This version uses a blocked algorithm to update sensitivities Abar
    in place. tril(Abar) should start containing L_bar, and will end
    containing the tril(Sigma_bar). Take tril(Abar) at the end if
    triu(Abar,1) did not start out filled with zeros. Alternatively,
    (tril(Abar) + tril(Abar).T) will give the symmetric, redundant
    matrix of sensitivities.
    
    Unlike the unblocked routine, if the upper triangular part of Abar
    started with non-zero values, some of these will be overwritten.

    If inplace=False, a copy of Abar is modified instead of the
    original. The Abar that was modified is returned.
    """
    if not inplace:
        Abar = Abar.copy()
    for k in range(L.shape[0], -1, -NB):
        j = max(0, k - NB)
        R, D, B, C = _level3partition(L, j, k)
        Rbar, Dbar, Bbar, Cbar = _level3partition(Abar, j, k)
        #Cbar[:] = Cbar @ inv(tril(D))
        Cbar[:] = _st(D, Cbar.T, trans=1).T
        Bbar -= np.dot(Cbar, R)
        Dbar[:] = tril(Dbar) - tril(np.dot(Cbar.T, C))
        #chol_unblocked_rev(D, Dbar, inplace=True) # slow in Python
        Dbar[:] = _chol_symbolic_rev(D, Dbar)
        Rbar -= (np.dot(Cbar.T, B) + np.dot(Dbar + Dbar.T, R))
    return Abar

	def connectionParameterMatrix(self, parameter):
		if utils.isCallable(parameter):
			matrix = parameter((self.noNodes, self.noNodes)) * self.connections
			matrix = np.tril(matrix) + np.tril(matrix).T # ensure symmetry
		else:
			matrix = parameter * np.ones((self.noNodes, self.noNodes)) * self.connections
		return matrix

 def getMechStiffStatistic(self, rangeK, minAA=0, AA='all'):
     """Return number of effective spring constant with set range of
     amino acids of protein structure.
     ``AA`` can be a list with a range of analysed amino acids as:
     [first_aa, last_aa, first_aa2, last_aa2],
     minAA - eliminate amino acids that are within 20aa and
     ``rangeK`` is a list [minK, maxK]"""
     
     model = self.getModel()
     if AA == 'all':
         sm = model.getStiffness()
     elif type(AA) == int:
         sm = model.getStiffness()[0: AA, (-1)*AA-1:-1]
     elif type(AA) == list and len(AA) == 1:
         sm = model.getStiffness()[0: AA, (-1)*AA-1:-1]
     elif type(AA) == list and len(AA) == 4:
         sm = model.getStiffness()[AA[0]:AA[1],AA[2]:AA[3]]
     if minAA > 0:
         sm2 = sm[minAA:-1,0:-1-minAA]  # matrix without close contacts
         sm3 = np.tril(sm2, k=-1)
         #sort_sm2 = np.sort((np.tril(sm2, k=-1)1).flatten())
         a = np.where(np.logical_and(sm3>rangeK[0], sm3<rangeK[1]))
     if minAA == 0:
         sm2 = np.tril(sm, k=-1)
         a = np.where(np.logical_and(sm2>rangeK[0], sm2<rangeK[1]))
     return len(a[0])

 def __init__(self, batch_size, mem_size, hidden_size):
     self.hidden_size = hidden_size
     self.mem_size = mem_size
     self.batch_size = batch_size
     N, M, d = batch_size, mem_size, hidden_size
     self.L = np.tril(np.ones([M, M], dtype='float32'))
     self.sL = np.tril(np.ones([M, M], dtype='float32'), k=-1)

def test_separate_independent_mok(session_tf):
    """
    We use different independent kernels for each of the output dimensions.
    We can achieve this in two ways:
        1) efficient: SeparateIndependentMok with Shared/SeparateIndependentMof
        2) inefficient: SeparateIndependentMok with InducingPoints
    However, both methods should return the same conditional,
    and after optimization return the same log likelihood.
    """
    # Model 1 (INefficient)
    q_mu_1 = np.random.randn(Data.M * Data.P, 1)
    q_sqrt_1 = np.tril(np.random.randn(Data.M * Data.P, Data.M * Data.P))[None, ...]  # 1 x MP x MP
    kern_list_1 = [RBF(Data.D, variance=0.5, lengthscales=1.2) for _ in range(Data.P)]
    kernel_1 = mk.SeparateIndependentMok(kern_list_1)
    feature_1 = InducingPoints(Data.X[:Data.M,...].copy())
    m1 = SVGP(Data.X, Data.Y, kernel_1, Gaussian(), feature_1, q_mu=q_mu_1, q_sqrt=q_sqrt_1)
    m1.set_trainable(False)
    m1.q_sqrt.set_trainable(True)
    m1.q_mu.set_trainable(True)
    gpflow.training.ScipyOptimizer().minimize(m1, maxiter=Data.MAXITER)

    # Model 2 (efficient)
    q_mu_2 = np.random.randn(Data.M, Data.P)
    q_sqrt_2 = np.array([np.tril(np.random.randn(Data.M, Data.M)) for _ in range(Data.P)])  # P x M x M
    kern_list_2 = [RBF(Data.D, variance=0.5, lengthscales=1.2) for _ in range(Data.P)]
    kernel_2 = mk.SeparateIndependentMok(kern_list_2)
    feature_2 = mf.SharedIndependentMof(InducingPoints(Data.X[:Data.M, ...].copy()))
    m2 = SVGP(Data.X, Data.Y, kernel_2, Gaussian(), feature_2, q_mu=q_mu_2, q_sqrt=q_sqrt_2)
    m2.set_trainable(False)
    m2.q_sqrt.set_trainable(True)
    m2.q_mu.set_trainable(True)
    gpflow.training.ScipyOptimizer().minimize(m2, maxiter=Data.MAXITER)

    check_equality_predictions(session_tf, [m1, m2])

    def hessian(self, x, lagrange, obj_factor):

        H = np.zeros((2*self._m, 2*self._m))
        H[:self._m, :self._m] = np.tril(np.tril(np.dot(self._A.T, self._A)))

        row, col = self.hessianstructure()

        return obj_factor*H[row, col]

def vec_to_sym(vec, shape):
    mask = np.tril(np.ones(shape)).astype(np.bool)
    sym = np.zeros(vec.shape[:-1] + mask.shape, vec.dtype)
    sym[..., mask] = vec
    sym -= (1 - np.sqrt(2))*np.diag(np.diag(sym))
    sym /= np.sqrt(2)
    sym += np.tril(sym, k=-1).T
    return sym

def absMDS(distmat, Z, weights, Vp):
	dZ = eucD(Z)
	dZ[dZ==0] = 1E-5
	bZ = Bcalc(weights, distmat, dZ)
	Xu = Vp.dot(bZ).dot(Z)
	dXu = eucD(Xu)
	stress = np.sqrt(np.tril(weights*(distmat-dXu)**2).sum() / np.tril(dXu**2).sum())
	return stress, Xu

def test_riemann():
    """Simple test of the Riemann matrix."""
    n = 10
    a = rogues.riemann(n)
    b = np.tril(-np.ones((n, n)), -1)
    c = np.tril(a, -1)
    # Kind of a goofy prop to check, but it's simple
    npt.assert_array_equal(b, c)