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)

def run_tril(dtype, shape, order, inplace):
    ac, ag = gen_gpuarray(shape, dtype, order=order, ctx=context)
    result = tril(ag, inplace=inplace)
    assert numpy.all(numpy.tril(ac) == result)
    if inplace:
        assert numpy.all(numpy.tril(ac) == ag)
    else:
        assert numpy.all(ac == ag)

def test_tfttr_trttf():
    """
    Test conversion routines between the Rectengular Full Packed (RFP) format
    and Standard Triangular Array (TR)
    """
    seed(1234)
    for ind, dtype in enumerate(DTYPES):
        n = 20
        if ind > 1:
            A_full = (rand(n, n) + rand(n, n)*1j).astype(dtype)
            transr = 'C'
        else:
            A_full = (rand(n, n)).astype(dtype)
            transr = 'T'

        trttf, tfttr = get_lapack_funcs(('trttf', 'tfttr'), dtype=dtype)
        A_tf_U, info = trttf(A_full)
        assert_(info == 0)
        A_tf_L, info = trttf(A_full, uplo='L')
        assert_(info == 0)
        A_tf_U_T, info = trttf(A_full, transr=transr, uplo='U')
        assert_(info == 0)
        A_tf_L_T, info = trttf(A_full, transr=transr, uplo='L')
        assert_(info == 0)

        # Create the RFP array manually (n is even!)
        A_tf_U_m = zeros((n+1, n//2), dtype=dtype)
        A_tf_U_m[:-1, :] = triu(A_full)[:, n//2:]
        A_tf_U_m[n//2+1:, :] += triu(A_full)[:n//2, :n//2].conj().T

        A_tf_L_m = zeros((n+1, n//2), dtype=dtype)
        A_tf_L_m[1:, :] = tril(A_full)[:, :n//2]
        A_tf_L_m[:n//2, :] += tril(A_full)[n//2:, n//2:].conj().T

        assert_array_almost_equal(A_tf_U, A_tf_U_m.reshape(-1, order='F'))
        assert_array_almost_equal(A_tf_U_T,
                                  A_tf_U_m.conj().T.reshape(-1, order='F'))

        assert_array_almost_equal(A_tf_L, A_tf_L_m.reshape(-1, order='F'))
        assert_array_almost_equal(A_tf_L_T,
                                  A_tf_L_m.conj().T.reshape(-1, order='F'))

        # Get the original array from RFP
        A_tr_U, info = tfttr(n, A_tf_U)
        assert_(info == 0)
        A_tr_L, info = tfttr(n, A_tf_L, uplo='L')
        assert_(info == 0)
        A_tr_U_T, info = tfttr(n, A_tf_U_T, transr=transr, uplo='U')
        assert_(info == 0)
        A_tr_L_T, info = tfttr(n, A_tf_L_T, transr=transr, uplo='L')
        assert_(info == 0)

        assert_array_almost_equal(A_tr_U, triu(A_full))
        assert_array_almost_equal(A_tr_U_T, triu(A_full))
        assert_array_almost_equal(A_tr_L, tril(A_full))
        assert_array_almost_equal(A_tr_L_T, tril(A_full))

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

  def __init__(self, s):
    """Initialize the elastic tensor from a string"""

    if not s:
      raise ValueError("no matrix was provided")

    # Remove braces and pipes
    s = s.replace("|", " ").replace("(", " ").replace(")", " ")

    # Remove empty lines
    lines = [line for line in s.split('\n') if line.strip()]
    if len(lines) != 6:
      raise ValueError("should have six rows")

    # Convert to float
    try:
      mat = [map(float, line.split()) for line in lines]
    except:
      raise ValueError("not all entries are numbers")

    # Make it into a square matrix
    mat = np.array(mat)
    if mat.shape != (6,6):
      # Is it upper triangular?
      if map(len, mat) == [6,5,4,3,2,1]:
	mat = [ [0]*i + mat[i] for i in range(6) ]
        mat = np.array(mat)

      # Is it lower triangular?
      if map(len, mat) == [1,2,3,4,5,6]:
	mat = [ mat[i] + [0]*(5-i) for i in range(6) ]
        mat = np.array(mat)

    if mat.shape != (6,6):
      raise ValueError("should be a square matrix")

    # Check that is is symmetric, or make it symmetric
    if la.norm(np.tril(mat, -1)) == 0:
      mat = mat + np.triu(mat, 1).transpose()
    if la.norm(np.triu(mat, 1)) == 0:
      mat = mat + np.tril(mat, -1).transpose()
    if la.norm(mat - mat.transpose()) > 0:
      raise ValueError("should be symmetric, or triangular")

    # Store it
    self.CVoigt = mat

    # Put it in a more useful representation
    self.SVoigt = la.inv(self.CVoigt)
    VoigtMat = [[0, 5, 4], [5, 1, 3], [4, 3, 2]]
    def SVoigtCoeff(p,q): return 1. / ((1+p/3)*(1+q/3))

    self.Smat = [[[[ SVoigtCoeff(VoigtMat[i][j], VoigtMat[k][l]) * self.SVoigt[VoigtMat[i][j]][VoigtMat[k][l]]
                     for i in range(3) ] for j in range(3) ] for k in range(3) ] for l in range(3) ]
    return

def passed_test(dtype, as_matrix, provide_C, uplo, trans):
    """
    Run one symmetric rank-2k update test.

    Arguments:
        dtype:        either 'float64' or 'float32', the NumPy dtype to test
        as_matrix:    True to test a NumPy matrix, False to test a NumPy ndarray
        provide_C:    True if C is to be provided to the BLASpy function, False otherwise
        uplo:         BLASpy 'uplo' parameter to test
        trans:        BLASpy 'trans' parameter to test

    Returns:
        True if the expected result is within the margin of error of the actual result,
        False otherwise.
    """

    transpose_a = trans == 't' or trans == 'T'
    upper = uplo == 'u' or uplo == 'U'

    # generate random sizes for matrix dimensions
    m_A = randint(N_MIN, N_MAX)
    n_A = randint(N_MIN, N_MAX)
    n = m_A if not transpose_a else n_A

    # create random scalars and matrices to test
    alpha = uniform(SCAL_MIN, SCAL_MAX)
    beta = uniform(SCAL_MIN, SCAL_MAX)
    A = random_matrix(m_A, n_A, dtype, as_matrix)
    B = random_matrix(m_A, n_A, dtype, as_matrix)
    C = random_symmetric_matrix(n, dtype, as_matrix) if provide_C else None

    # create a copy of  C that can be used to calculate the expected result
    C_2 = copy(C) if C is not None else zeros((n, n))

    # compute the expected result
    if not transpose_a:
        C_2 = (beta * C_2) + (alpha * dot(A, B.T)) + (alpha * dot(B, A.T))
    else:
        C_2 = (beta * C_2) + (alpha * dot(A.T, B)) + (alpha * dot(B.T, A))

    # ensure C and C_2 are upper or lower triangular representations of symmetric matrices
    if upper:
        C_2 = triu(C_2)
        if provide_C:
            C = triu(C)
    else:
        C_2 = tril(C_2)
        if provide_C:
            C = tril(C)

    # get the actual result
    C = syr2k(A, B, C, uplo, trans, alpha, beta)

    # compare the actual result to the expected result and return result of the test
    return allclose(C, C_2, RTOL, ATOL)