def BipartiteGraphClusteringK(dist_mat,prefix1,prefix2,K,thresh):
    # build the graph
    G,Local_Seq,Long_Seq = GraphBuild(dist_mat)
    # compute graph laplacian spectrum
    U,S,V,D1,D2 = GraphSVD(G,Local_Seq,Long_Seq)
    # update K if necessary
    #if ((S>0.7).sum()<K):
    #    K = (S>0.7).sum()
    #
    # data matrix
    Z0 = np.round(la.sqrtm(la.inv(D1))*U[:,np.arange(K)],8)
    Z1 = np.round(la.sqrtm(la.inv(D2))*V[:,np.arange(K)],8)
    # compute the centroid
    estimator = KMeans(init='k-means++', n_clusters=K, n_init=10)
    estimator.fit(Z0)
    B = estimator.cluster_centers_
    # compute the distance between long seq and centroid
    Y = cdist(Z1,B)
    # compute the assignment of long seq
    C = Y.argmin(axis=1)
    # output the pairs of long-seqs
    for k in np.arange(K):
        print "#%d"%k + "\t" + ":".join(['%.2f']*len(S)) % tuple(S)
        I = np.where(C==k)[0]
        for a in Long_Seq[I]:
            print a

def solveBlockGlasso(signal):
    start = int(signal[0]) # include
    S_Matrix  = S_Matrix_bc.value
    W_matrix = W_Matrix_bc.value
    old_W = np.copy(W_matrix)
    end   = min(int(signal[1]),S_Matrix.shape[0]) # non-inclusive
    deltamatrix = np.zeros(S_Matrix.shape)
    NN = S_Matrix.shape[0]
    for n in range(start,end):
        W11 = np.delete(W_matrix,n,0)
        W11 = np.delete(W11,n,1)
        Z   = linalg.sqrtm(W11)

        s11 = S_Matrix[:,n]
        s11 = np.delete(s11,n)
        Y   = np.dot(nplinalg.inv(linalg.sqrtm(W11)),s11)
	Y = np.real(Y)
	Z = np.real(Z)
	B = lasso(Z,Y,beta_value)

    updated_column = np.dot(W11,B)

    matrix_ind = np.array(range(0,NN))
    matrix_ind = np.delete(matrix_ind,n)
    column_ind = 0
    for k in matrix_ind:
        deltamatrix[k,n]=updated_column[column_ind] - W_matrix[k,n]
        deltamatrix[n,k]=updated_column[column_ind] - W_matrix[k,n]
	    W_matrix[k,n] = updated_column[column_ind]
	    W_matrix[n,k] = updated_column[column_ind]
        column_ind = column_ind+1

def gibbs1_swr(S0, P0, P1, T):
    """
    function SA = GIBBS1_SWR(S0,P0,P1,T);

    This file executes the Carter-Kohn backward sampler for the
    Stock-Watson-Romer model.

    S0, P0, P1 are outputs of the forward Kalman filter

    VERIFIED (1x) SL (8-9-13)
    """
    A = np.array([[0., 1.], [0., 1.]])

    # initialize arrays for Gibbs sampler
    SA = zeros((2, T))  # artificial states
    SM = zeros((2, 1))  # backward update for conditional mean of state vector
    PM = zeros((2, 2))  # backward update for projection matrix
    P = zeros((2, 2))  # backward update for conditional variance matrix
    wa = np.random.randn(2, T)  # draws for state innovations

    # Backward recursions and sampling
    # Terminal state
    SA[:, T-1] = S0[:, T-1] + np.real(sqrtm(P0[:, :, T-1])).dot(wa[:, T-1])

    # iterating back through the rest of the sample
    for i in range(1, T):
        PM = np.dot(P0[:, :, T-i].dot(A.T), inv(P1[:, :, T-i]))
        P = P0[:, :, T-i] - np.dot(PM.dot(A), P0[:, :, T-i])
        SM = S0[:, T-i-1] + PM.dot(SA[:, T-i] - A.dot(S0[:, T-i-1]))
        SA[:, T-i-1] = SM + np.real(sqrtm(P)).dot(wa[:, T-i-1])

    return SA

def rcca(X, Y, reg):   
    X = np.array(X)
    Y = np.array(Y)
    n, p = X.shape
    n, q = Y.shape
        
    # zero mean
    X = X - X.mean(axis=0)
    Y = Y - Y.mean(axis=0)
        
    # covariances
    S = np.cov(X.T, Y.T, bias=1)
    
    SXX = S[:p,:p]
    SYY = S[p:,p:]
    SXY = S[:p,p:]

    ux, lx, vxt = SLA.svd(SXX)
    #uy, ly, vyt = SLA.svd(SYY)
    
    #print lx
    #正則化
    Rg = np.diag(np.ones(p)*reg)
    SXX = SXX + Rg
    SYY = SYY + Rg

    sqx = SLA.sqrtm(SLA.inv(SXX)) # SXX^(-1/2)
    sqy = SLA.sqrtm(SLA.inv(SYY)) # SYY^(-1/2)
    M = np.dot(np.dot(sqx, SXY), sqy.T) # SXX^(-1/2) * SXY * SYY^(-T/2)
    A, r, Bh = SLA.svd(M, full_matrices=False)
    B = Bh.T      
    #r = self.reg*r
    return r[0], lx[0], lx[0] 

    def stdNorm(self, U1, U2):
        print "U1"
        print U1
        print "U2"
        print U2

        mat1 = np.matrix(U1).T
        print mat1
        print mat1.mean(axis=1)
        mat1 = mat1 - mat1.mean(axis=1)
        print mat1
        mat1cov = np.cov(mat1)
        print mat1cov
        p1,l1,p1t = NLA.svd(mat1cov)
        print p1
        print l1
        print p1t
        l1sq = SLA.sqrtm(SLA.inv(np.diag(l1))) 
        snU1 =  np.dot(np.dot(l1sq, p1.T), mat1)

        mat2 = np.matrix(U2).T
        mat2 = mat2 - mat2.mean(axis=1)
        mat2cov = np.cov(mat2)
        p2,l2,p2t = NLA.svd(mat2cov)
        l2sq = SLA.sqrtm(SLA.inv(np.diag(l2))) 
        snU2 =  np.dot(np.dot(l2sq, p2.T), mat2)

        print "cov:"
        print np.cov(snU1)
        print np.cov(snU2)

        return snU1, snU2

 def cca1(self, X, Y):
     '''
     正準相関分析
     http://en.wikipedia.org/wiki/Canonical_correlation
     '''    
     X = np.array(X)
     Y = np.array(Y)
     n, p = X.shape
     n, q = Y.shape
     
     # zero mean
     X = X - X.mean(axis=0)
     Y = Y - Y.mean(axis=0)
     
     # covariances
     S = np.cov(X.T, Y.T)
     
     # S = np.corrcoef(X.T, Y.T)
     SXX = S[:p,:p]
     SYY = S[p:,p:]
     SXY = S[:p,p:]
     SYX = S[p:,:p]
     
     # 
     sqx = SLA.sqrtm(SLA.inv(SXX)) # SXX^(-1/2)
     sqy = SLA.sqrtm(SLA.inv(SYY)) # SYY^(-1/2)
     M = np.dot(np.dot(sqx, SXY), sqy.T) # SXX^(-1/2) * SXY * SYY^(-T/2)
     A, s, Bh = SLA.svd(M, full_matrices=False)
     B = Bh.T      
     #print np.dot(np.dot(A[:,0].T,SXX),A[:,0])
     return s, A, B

    def __init__(self, F, H, R_v, R_n, x_hat_0, P_0, alpha=1):
        """
        x_hat_0 and P_0 are initial estimates of the state mean and covariance.
        alpha determines the spread of sigma points, and should be small...
        """
        self.F = F
        self.H = H
        self.R_v = R_v
        self.R_n = R_n
        self.x_hat = x_hat_0
        self.P = P_0
        self.S = spla.sqrtm(P_0) #np.linalg.cholesky(P_0)
        self.sqrt_R_v = spla.sqrtm(R_v) #np.linalg.cholesky(R_v)
        self.sqrt_R_n = spla.sqrtm(R_n) #np.linalg.cholesky(R_n)
        self.L = len(x_hat_0)
        self.M = len(H(x_hat_0))
        self.alpha = alpha
        self.beta = 2
        kappa = 3 - self.L
        self.kappa = kappa
        self.Lambda = alpha**2 * (self.L + kappa) - self.L

        # NB: these don't change while we don't augment the sigma points
        self.weights_m = self.weights_m()
        self.weights_c = self.weights_c()

def gsvd(a, m, w):
	"""
	:param a: Matrix to GSVD
	:param m: 1st Constraint, (u.T * m * u) = I
	:param w: 2nd Constraint, (v.T * w * v) = I
	:return: (u ,s, v)
	"""

	(aHeight, aWidth) = a.shape
	(mHeight, mWidth) = m.shape
	(wHeight, mWidth) = w.shape

	assert(aHeight == mHeight)
	assert(aWidth == mWidth)

	mSqrt = sqrtm(m)
	wSqrt = sqrtm(w)


	mSqrtInv = np.linalg.inv(mSqrt)
	wSqrtInv = np.linalg.inv(wSqrt)

	_a = np.dot(np.dot(mSqrt, a), wSqrt)

	(_u, _s, _v) = np.linalg.svd(_a)

	u = np.dot(mSqrtInv, _u)
	v = np.dot(wSqrtInv, _v.T).T
	s = _s

	return (u, s, v)

def kcca(X, Y, kernel_x=gaussian_kernel, kernel_y=gaussian_kernel, eta=1.0):
    '''
    カーネル正準相関分析
    http://staff.aist.go.jp/s.akaho/papers/ibis00.pdf
    '''
    n, p = X.shape
    n, q = Y.shape

    Kx = DIST.squareform(DIST.pdist(X, kernel_x))
    Ky = DIST.squareform(DIST.pdist(Y, kernel_y))
    J = np.eye(n) - np.ones((n, n)) / n
    M = np.dot(np.dot(Kx.T, J), Ky) / n
    L = np.dot(np.dot(Kx.T, J), Kx) / n + eta * Kx
    N = np.dot(np.dot(Ky.T, J), Ky) / n + eta * Ky

    sqx = LA.sqrtm(LA.inv(L))
    sqy = LA.sqrtm(LA.inv(N))

    a = np.dot(np.dot(sqx, M), sqy.T)
    A, s, Bh = LA.svd(a, full_matrices=False)
    B = Bh.T

    # U = np.dot(np.dot(A.T, sqx), X).T
    # V = np.dot(np.dot(B.T, sqy), Y).T

    return s, A, B

    def cca(self, X, Y):
        '''
        正準相関分析
        http://en.wikipedia.org/wiki/Canonical_correlation
        '''    
        #X = np.array(X)
        #Y = np.array(Y)
        n, p = X.shape
        n, q = Y.shape
        
        # zero mean
        X = X - X.mean(axis=0)
        Y = Y - Y.mean(axis=0)
        
        # covariances
        S = np.cov(X.T, Y.T, bias=1)
        
        SXX = S[:p,:p]
        SYY = S[p:,p:]
        SXY = S[:p,p:]

        #正則化
        SXX = self.add_reg(SXX, self.reg) 
        SYY = self.add_reg(SYY, self.reg)

        sqx = SLA.sqrtm(SLA.inv(SXX)) # SXX^(-1/2)
        sqy = SLA.sqrtm(SLA.inv(SYY)) # SYY^(-1/2)
        M = np.dot(np.dot(sqx, SXY), sqy.T) # SXX^(-1/2) * SXY * SYY^(-T/2)
        A, r, Bh = SLA.svd(M, full_matrices=False)
        B = Bh.T      
        #r = self.reg*r
        return r, A, B

def quad_form(x, P):
    x,P = map(Expression.cast_to_const, (x,P))
    # Check dimensions.
    n = P.size[0]
    if P.size[1] != n or x.size != (n,1):
        raise Exception("Invalid dimensions for arguments.")
    if x.curvature.is_constant():
        return x.T*P*x
    elif P.curvature.is_constant():
        np_intf = intf.get_matrix_interface(np.ndarray)
        P = np_intf.const_to_matrix(P.value)
        # Replace P with symmetric version.
        P = (P + P.T)/2
        # Check if P is PSD.
        eigvals = LA.eigvalsh(P)
        if min(eigvals) > 0:
            P_sqrt = Constant(LA.sqrtm(P).real)
            return square(norm2(P_sqrt*x))
        elif max(eigvals) < 0:
            P_sqrt = Constant(LA.sqrtm(-P).real)
            return -square(norm2(P_sqrt*x))
        else:
            raise Exception("P has both positive and negative eigenvalues.")
    else:
        raise Exception("At least one argument to quad_form must be constant.")

    def dataNorm(self):
        SXX = np.cov(self.X)
        U, l, Ut = LA.svd(SXX, full_matrices=True) 
        H = np.dot(LA.sqrtm(LA.inv(np.diag(l))),Ut)
        self.nX = np.dot(H,self.X)

        #print np.cov(self.nX)
        #print "mean:"
        #print np.mean(self.nX)

        SYY = np.cov(self.Y)
        U, l, Ut = LA.svd(SYY, full_matrices=True) 
        H = np.dot(LA.sqrtm(LA.inv(np.diag(l))),Ut)
        #print "H"
        #print H
        self.nY = np.dot(H,self.Y)
        #print np.cov(self.nY)

        print "dataNorm_X:"
        for i in range(len(self.nX)):
            print(self.nX[i])
        print("---")

        print "dataNorm_Y:"
        for i in range(len(self.nY)):
            print(self.nY[i])
        print("---")

def cca(H1, H2):
    H1bar = copy.deepcopy(H1)
    H1bar = H1bar-H1bar.mean(axis=0)
    H2bar = copy.deepcopy(H2)
    H2bar = H2bar-H2bar.mean(axis=0)
    H1bar = H1bar.T
    H2bar = H2bar.T
    H1bar += np.random.random(H1bar.shape)*0.00001
    H2bar += np.random.random(H2bar.shape)*0.00001
    r1 = 0.00000001
    m = H1.shape[0]
    #H1bar = H1 - (1.0/m)*np.dot(H1, np.ones((m,m), dtype=np.float32))
    #H2bar = H2 - (1.0/m)*np.dot(H2, np.ones((m,m), dtype=np.float32))
    SigmaHat12 = (1.0/(m-1))*np.dot(H1bar, H2bar.T)
    SigmaHat11 = (1.0/(m-1))*np.dot(H1bar, H1bar.T)
    SigmaHat11 = SigmaHat11 + r1*np.identity(SigmaHat11.shape[0], dtype=np.float32)
    SigmaHat22 = (1.0/(m-1))*np.dot(H2bar, H2bar.T)
    SigmaHat22 = SigmaHat22 + r1*np.identity(SigmaHat22.shape[0], dtype=np.float32)
    SigmaHat11_2=mat_pow(SigmaHat11).real.astype(np.float32)
    SigmaHat22_2=mat_pow(SigmaHat22).real.astype(np.float32)
    ##TMP = np.dot(SigmaHat12, SigmaHat22_2) #unstable
    TMP2 = stable_inverse_A_dot_Bneg1(SigmaHat12, sqrtm(SigmaHat22))#np.dot(SigmaHat12, SigmaHat22_2)
    TMP3 = stable_inverse_A_dot_Bneg1_cholesky(SigmaHat12, sqrtm(SigmaHat22))#np.dot(SigmaHat12, SigmaHat22_2)

    ##Tval = np.dot(SigmaHat11_2, TMP) #unstable
    Tval = stable_inverse_Aneg1_dot_B(sqrtm(SigmaHat11), TMP2)
    Tval3 = stable_inverse_Aneg1_dot_B_cholesky(sqrtm(SigmaHat11), TMP3)

    ##U, D, V, = np.linalg.svd(Tval)

    ## corr =  np.trace(np.dot(Tval.T, Tval))**(0.5) #wrong
    corr =  np.trace(sqrtm(np.dot(Tval.T, Tval)))
    return corr

def pca_similarity(covar_a, covar_b):
    """
        Calculates the similarity between the two covariance matrices

        **Arguments:**

        covar_a
            The first covariance matrix.

        covar_b
            The second covariance matrix.
    """
    # Take the square root of the symmetric matrices
    a_sq = spla.sqrtm(covar_a)
    b_sq = spla.sqrtm(covar_b)

    # Check for imaginary entries
    for mat in [a_sq, b_sq]:
        max_imag = np.amax(np.abs(np.imag(mat)))
        mean_real = np.mean(np.abs(np.real(mat)))
        if(max_imag/mean_real > 1e-6):
            Warning('Covariance matrix is not diagonally dominant')

    # Return the PCA similarity (1 - PCA distance)
    return 1 - np.sqrt(np.trace(np.dot(a_sq-b_sq, a_sq-b_sq))/(np.trace(covar_a+covar_b)))

    def test_sqrtm_type_preservation_and_conversion(self):
        # The sqrtm matrix function should preserve the type of a matrix
        # whose eigenvalues are nonnegative with zero imaginary part.
        # Test this preservation for variously structured matrices.
        complex_dtype_chars = ('F', 'D', 'G')
        for matrix_as_list in (
                [[1, 0], [0, 1]],
                [[1, 0], [1, 1]],
                [[2, 1], [1, 1]],
                [[2, 3], [1, 2]],
                [[1, 1], [1, 1]]):

            # check that the spectrum has the expected properties
            W = scipy.linalg.eigvals(matrix_as_list)
            assert_(not any(w.imag or w.real < 0 for w in W))

            # check float type preservation
            A = np.array(matrix_as_list, dtype=float)
            A_sqrtm, info = sqrtm(A, disp=False)
            assert_(A_sqrtm.dtype.char not in complex_dtype_chars)

            # check complex type preservation
            A = np.array(matrix_as_list, dtype=complex)
            A_sqrtm, info = sqrtm(A, disp=False)
            assert_(A_sqrtm.dtype.char in complex_dtype_chars)

            # check float->complex type conversion for the matrix negation
            A = -np.array(matrix_as_list, dtype=float)
            A_sqrtm, info = sqrtm(A, disp=False)
            assert_(A_sqrtm.dtype.char in complex_dtype_chars)

def MSI(X, Y):
    try: 
        n, p = X.shape
        n, q = Y.shape
    except ValueError as v:
        return []

    X = X.astype('float32', copy=False)
    X -= X.mean(axis=0)
    X /= np.max(np.abs(X))
    Y = Y.astype('float32', copy=False)
    Y -= Y.mean(axis=0)
    Y /= np.max(np.abs(Y))

    C = np.cov(X.T, Y.T, bias=1)
    CXX = C[:p,:p]
    CYY = C[p:,p:]

    sqx,_ = LA.sqrtm(LA.inv(CXX),False) # SXX^(-1/2)
    sqy,_ = LA.sqrtm(LA.inv(CYY),False) # SYY^(-1/2)

    # build square matrix
    u1 = np.vstack((sqx, np.zeros((sqy.shape[0], sqx.shape[1]))))
    u2 = np.vstack((np.zeros((sqx.shape[0], sqy.shape[1])), sqy))
    U = np.hstack((u1, u2))
    
    R = np.dot(np.dot(U, C), U.T)

    eigvals = LA.eigh(R)[0]
    eigvals /= np.sum(eigvals)
    # Compute index
    return 1 + np.sum(eigvals * np.log(eigvals)) / np.log(eigvals.shape[0])

    def kcca(self, X, Y, kernel_x=gaussian_kernel, kernel_y=gaussian_kernel, eta=1.0):
        n, p = X.shape
        n, q = Y.shape
        
        Kx = DIST.squareform(DIST.pdist(X, kernel_x))
        Ky = DIST.squareform(DIST.pdist(Y, kernel_y))
        J = np.eye(n) - np.ones((n, n)) / n
        M = np.dot(np.dot(Kx.T, J), Ky) / n
        L = np.dot(np.dot(Kx.T, J), Kx) / n + eta * Kx
        N = np.dot(np.dot(Ky.T, J), Ky) / n + eta * Ky


        sqx = SLA.sqrtm(SLA.inv(L))
        sqy = SLA.sqrtm(SLA.inv(N))
        
        a = np.dot(np.dot(sqx, M), sqy.T)
        A, s, Bh = SLA.svd(a, full_matrices=False)
        B = Bh.T
        
        # U = np.dot(np.dot(A.T, sqx), X).T
        # V = np.dot(np.dot(B.T, sqy), Y).T
        print s.shape
        print A.shape
        print B.shape
        return s, A, B

def calc_vecp(l,C_l_hat,C_fl, C_l):

    C_fl_12 = sqrtm(C_fl[l])
    C_l_inv = LA.inv(C_l[l])
    C_l_inv_12= sqrtm(C_l_inv)
    # the order is inverted compared to matlab hamimeche_lewis_likelihood.m line 19

    # line 20 of hamimeche_lewis_likelihood.m
    res = np.dot(C_l_inv_12, np.dot(C_l_hat[l], C_l_inv_12))
    [d, u] = LA.eigh(res)
    d = np.diag(d)  # noticed that python returns the eigenvalues as a vector, not a matrix
    #np. dot( u, np.dot( np.diag(d), LA.inv(u))) should be equals to res
    # real symmetric matrices are diagnalized by orthogonal matrices (M^t M = 1)

    # this makes a diagonal matrix by applying g(x) to the eigenvalues, equation 10 in Barkats et al
    gd = np.sign(np.diag(d) - 1) * np.sqrt(2 * (np.diag(d) - np.log(np.diag(d)) - 1))
    gd = np.diag(gd)
    # Argument of vecp in equation 8; multiplying from right to left
    X = np.dot(np.transpose(u), C_fl_12)
    X = np.dot(gd, X)
    X = np.dot(u, X)
    X = np.dot(C_fl_12, X)
    # This is the vector of equation 7
    X = vecp(X)

    return X

 def compute_induced_kernel_matrix_on_data(self,data_x,data_y):
     '''Z follows the same distribution as X; W follows that of Y.
     The current data generating methods we use 
     generate X and Y at the same time. '''
     size_induced_set = max(self.num_inducex,self.num_inducey)
     #print "size_induce_set", size_induced_set
     if self.data_generator is None:
         subsample_idx = np.random.randint(self.num_samples, size=size_induced_set)
         self.data_z = data_x[subsample_idx,:]
         self.data_w = data_y[subsample_idx,:]
     else:
         self.data_z, self.data_w = self.data_generator(size_induced_set)
         self.data_z[[range(self.num_inducex)],:]
         self.data_w[[range(self.num_inducey)],:]
     #print 'Induce Set'
     if self.kernelX_use_median:
         sigmax = self.kernelX.get_sigma_median_heuristic(data_x)
         self.kernelX.set_width(float(sigmax))
     if self.kernelY_use_median:
         sigmay = self.kernelY.get_sigma_median_heuristic(data_y)
         self.kernelY.set_width(float(sigmay))
     Kxz = self.kernelX.kernel(data_x,self.data_z)
     Kzz = self.kernelX.kernel(self.data_z)
     #R = inv(sqrtm(Kzz))
     R = inv(sqrtm(Kzz + np.eye(np.shape(Kzz)[0])*10**(-6)))
     phix = Kxz.dot(R)
     Kyw = self.kernelY.kernel(data_y,self.data_w)
     Kww = self.kernelY.kernel(self.data_w)
     #S = inv(sqrtm(Kww))
     S = inv(sqrtm(Kww + np.eye(np.shape(Kww)[0])*10**(-6)))
     phiy = Kyw.dot(S)
     return phix, phiy

def GraphSVD(G,row_nodes,column_nodes):
    A = bipartite.biadjacency_matrix(G,row_order=row_nodes,column_order=column_nodes)
    D1 = np.diag(np.squeeze(np.asarray(A.sum(1))))
    D2 = np.diag(np.squeeze(np.asarray(A.sum(0))))
    An = la.sqrtm(la.inv(D1))*A*la.sqrtm(la.inv(D2))
    U,S,V = np.linalg.svd(An)
    V = V.T
    return (U,S,V,D1,D2)