def test_type_preservation_and_conversion(self):
        # The fractional_matrix_power matrix function should preserve
        # the type of a matrix whose eigenvalues
        # are positive 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]]):

            # 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 various positive and negative powers
            # with absolute values bigger and smaller than 1.
            for p in (-2.4, -0.9, 0.2, 3.3):

                # check float type preservation
                A = np.array(matrix_as_list, dtype=float)
                A_power = fractional_matrix_power(A, p)
                assert_(A_power.dtype.char not in complex_dtype_chars)

                # check complex type preservation
                A = np.array(matrix_as_list, dtype=complex)
                A_power = fractional_matrix_power(A, p)
                assert_(A_power.dtype.char in complex_dtype_chars)

                # check float->complex for the matrix negation
                A = -np.array(matrix_as_list, dtype=float)
                A_power = fractional_matrix_power(A, p)
                assert_(A_power.dtype.char in complex_dtype_chars)

	def fidelity(self, m, n):
		"""Compute the fidelity between the density matrces m and n.

		:param numpy_array m: Density matrix
		:param numpy_array n: Density matrix

		:return: The fideltiy between m and n (:math:`\mathrm{Tr}(\sqrt{\sqrt{m}n\sqrt{m}})^2`).
		:rtype: complex
		"""
		sqrt_m = fractional_matrix_power(m, 0.5)
		F = np.trace(fractional_matrix_power(np.dot(sqrt_m,np.dot(n, sqrt_m)), 0.5))**2.0

		return F

    def get_mean_and_covariance(self, w2v, num_of_occurences):
        """ get mean and covariance of words vectors over the training set of word2vec model
            w2v -- word2vec model (in matrix form)
            num_of_occurences -- array that specifies weights for averaging over words
        """
        weights = num_of_occurences/np.sum(num_of_occurences)
        try:
            w2v_temp = np.multiply(w2v, weights)
        except MemoryError:
            w2v_temp = np.copy(w2v)
            for wn in range(np.shape(w2v)[1]):
                w2v_temp[:, wn] *= weights[wn]
        self.Mean = np.sum(w2v_temp, 1)
        try:
            w2v_except0 = w2v - np.reshape(self.Mean, (len(self.Mean), 1))
        except MemoryError:
            w2v_except0 = w2v_temp  # just to set the right shape (to avoid memoryError)
            for wn in range(np.shape(w2v)[1]):
                w2v_except0[:, wn] = w2v[:, wn] - self.Mean

        try:
            w2v_normalized = np.multiply(w2v_except0, np.power(weights, 0.5))
        except MemoryError:
            w2v_normalized = w2v_except0
            for wn in range(np.shape(w2v_except0)[1]):
                w2v_normalized[:, wn] *= weights[wn]**0.5

        self.Cov = np.dot(w2v_normalized, np.transpose(w2v_normalized))
        self.Cov = self.Cov/np.shape(w2v)[1]
        self.SqrtCov = fractional_matrix_power(self.Cov, -0.5)

 def test_larger_abs_fractional_matrix_powers(self):
     np.random.seed(1234)
     for n in (2, 3, 5):
         for i in range(10):
             M = np.random.randn(n, n) + 1j * np.random.randn(n, n)
             M_one_fifth = fractional_matrix_power(M, 0.2)
             # Test the round trip.
             M_round_trip = np.linalg.matrix_power(M_one_fifth, 5)
             assert_allclose(M, M_round_trip)
             # Test a large abs fractional power.
             X = fractional_matrix_power(M, -5.4)
             Y = np.linalg.matrix_power(M_one_fifth, -27)
             assert_allclose(X, Y)
             # Test another large abs fractional power.
             X = fractional_matrix_power(M, 3.8)
             Y = np.linalg.matrix_power(M_one_fifth, 19)
             assert_allclose(X, Y)

    def test_singular(self):
        # Negative fractional powers do not work with singular matrices.
        for matrix_as_list in (
                [[0, 0], [0, 0]],
                [[1, 1], [1, 1]],
                [[1, 2], [3, 6]],
                [[0, 0, 0], [0, 1, 1], [0, -1, 1]]):

            # Check fractional powers both for float and for complex types.
            for newtype in (float, complex):
                A = np.array(matrix_as_list, dtype=newtype)
                for p in (-0.7, -0.9, -2.4, -1.3):
                    A_power = fractional_matrix_power(A, p)
                    assert_(np.isnan(A_power).all())
                for p in (0.2, 1.43):
                    A_power = fractional_matrix_power(A, p)
                    A_round_trip = fractional_matrix_power(A_power, 1/p)
                    assert_allclose(A_round_trip, A)

 def test_round_trip_random_complex(self):
     np.random.seed(1234)
     for p in range(1, 5):
         for n in range(1, 5):
             M_unscaled = np.random.randn(n, n) + 1j * np.random.randn(n, n)
             for scale in np.logspace(-4, 4, 9):
                 M = M_unscaled * scale
                 M_root = fractional_matrix_power(M, 1/p)
                 M_round_trip = np.linalg.matrix_power(M_root, p)
                 assert_allclose(M_round_trip, M)

    def test_al_mohy_higham_2012_experiment_1(self):
        # Fractional powers of a tricky upper triangular matrix.
        A = _get_al_mohy_higham_2012_experiment_1()

        # Test remainder matrix power.
        A_funm_sqrt, info = funm(A, np.sqrt, disp=False)
        A_sqrtm, info = sqrtm(A, disp=False)
        A_rem_power = _matfuncs_inv_ssq._remainder_matrix_power(A, 0.5)
        A_power = fractional_matrix_power(A, 0.5)
        assert_array_equal(A_rem_power, A_power)
        assert_allclose(A_sqrtm, A_power)
        assert_allclose(A_sqrtm, A_funm_sqrt)

        # Test more fractional powers.
        for p in (1/2, 5/3):
            A_power = fractional_matrix_power(A, p)
            A_round_trip = fractional_matrix_power(A_power, 1/p)
            assert_allclose(A_round_trip, A, rtol=1e-2)
            assert_allclose(np.tril(A_round_trip, 1), np.tril(A, 1))

    def test_type_conversion_mixed_sign_or_complex_spectrum(self):
        complex_dtype_chars = ("F", "D", "G")
        for matrix_as_list in ([[1, 0], [0, -1]], [[0, 1], [1, 0]], [[0, 1, 0], [0, 0, 1], [1, 0, 0]]):

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

            # Check various positive and negative powers
            # with absolute values bigger and smaller than 1.
            for p in (-2.4, -0.9, 0.2, 3.3):

                # check complex->complex
                A = np.array(matrix_as_list, dtype=complex)
                A_power = fractional_matrix_power(A, p)
                assert_(A_power.dtype.char in complex_dtype_chars)

                # check float->complex
                A = np.array(matrix_as_list, dtype=float)
                A_power = fractional_matrix_power(A, p)
                assert_(A_power.dtype.char in complex_dtype_chars)

 def fit(self,X,Y):
     self._compute_covariance(X,Y)
     S_11_ = fractional_matrix_power(self.S_11,-0.5)
     
     S_22_ = fractional_matrix_power(self.S_22,-0.5)
     
     T = np.dot(np.dot(S_11_,self.S_12),S_22_)
     
     U ,S ,V = np.linalg.svd(T)
     
     self.U = U[:,:self.k]
     self.S = S[:self.k]
     self.V = V[:self.k,:]
     
     self.A = np.dot(S_11_,U)
     self.B = np.dot(S_22_,V.T)
     self.A = self.A[:,:self.k]
     self.B = self.B[:,:self.k]
     
     self.coeff_ = np.dot(self.A,self.V)
     
     return self

 def test_round_trip_random_float(self):
     # This test is more annoying because it can hit the branch cut;
     # this happens when the matrix has an eigenvalue
     # with no imaginary component and with a real negative component,
     # and it means that the principal branch does not exist.
     np.random.seed(1234)
     for p in range(1, 5):
         for n in range(1, 5):
             M_unscaled = np.random.randn(n, n)
             for scale in np.logspace(-4, 4, 9):
                 M = M_unscaled * scale
                 M_root = fractional_matrix_power(M, 1/p)
                 M_round_trip = np.linalg.matrix_power(M_root, p)
                 assert_allclose(M_round_trip, M)

    def test_singular(self):
        # Negative fractional powers do not work with singular matrices.
        # Neither do non-integer fractional powers,
        # because the scaling and squaring cannot deal with it.
        for matrix_as_list in (
                [[0, 0], [0, 0]],
                [[1, 1], [1, 1]],
                [[1, 2], [3, 6]],
                [[0, 0, 0], [0, 1, 1], [0, -1, 1]]):

            # check that the spectrum has the expected properties
            W = scipy.linalg.eigvals(matrix_as_list)
            assert_(np.count_nonzero(W) < len(W))

            # check fractional powers both for float and for complex types
            for newtype in (float, complex):
                A = np.array(matrix_as_list, dtype=newtype)
                for p in (-0.7, -0.9, -2.4, -1.3):
                    A_power = fractional_matrix_power(A, p)
                    assert_(np.isnan(A_power).all())
                for p in (0.2, 1.43):
                    A_power = fractional_matrix_power(A, p)
                    A_round_trip = fractional_matrix_power(A_power, 1/p)
                    assert_allclose(A_round_trip, A)

    def test_random_matrices_and_powers(self):
        # Each independent iteration of this fuzz test picks random parameters.
        # It tries to hit some edge cases.
        np.random.seed(1234)
        nsamples = 20
        for i in range(nsamples):
            # Sample a matrix size and a random real power.
            n = random.randrange(1, 5)
            p = np.random.randn()

            # Sample a random real or complex matrix.
            matrix_scale = np.exp(random.randrange(-4, 5))
            A = np.random.randn(n, n)
            if random.choice((True, False)):
                A = A + 1j * np.random.randn(n, n)
            A = A * matrix_scale

            # Check a couple of analytically equivalent ways
            # to compute the fractional matrix power.
            # These can be compared because they both use the principal branch.
            A_power = fractional_matrix_power(A, p)
            A_logm, info = logm(A, disp=False)
            A_power_expm_logm = expm(A_logm * p)
            assert_allclose(A_power, A_power_expm_logm)

temp = np.array(xrange(1829))+1
adDup = np.hstack((temp[:,np.newaxis],adMat))
temp = np.array(xrange(1830))
adDup = np.vstack((temp,adDup))
np.savetxt("aff.csv", adDup, delimiter=';')

# social similarity matrix
socialsim = sklearn.metrics.pairwise.pairwise_kernels(social,metric='rbf',gamma=0.8)

adMat = adMat+socialsim

# Make D matrix
D = np.diag(adMat.sum(1))

# Make laplacian
Dinv = la.fractional_matrix_power(D,-0.5)
L = np.eye(numProj) - np.dot(np.dot(Dinv,adMat),Dinv)

#L = D - adMat

# Cal eigenvector of laplacian
eigval,eigvec = np.linalg.eigh(L)

# Sort from smallest to largest
idx = np.argsort(eigval)
eigval = eigval[idx]
eigvec = eigvec[:,idx]

 Perform PCA
uXPca,sXPca,wXPca = la.svd(eigvec, full_matrices=False)
K=500

def sq_IG_distance(cov_1, cov_2) :
    
    cov_1_pow = fractional_matrix_power(cov_1, -0.5)
    return norm(logm(np.linalg.multi_dot([cov_1_pow, cov_2, cov_1_pow])), ord='fro') ** 2

#adjMat[2,9] = 1
#adjMat[5,0] = 1
#adjMat[0,5] = 1
#adjMat[1,8] = 1
#adjMat[8,1] = 1

#print adjMat
#print(np.sum(adjMat))
#show_graph(adjMat)

# Create diagonal matrix
diagMat = np.diag(adjMat.sum(axis=0))
#print diagMat

# Create Laplacian matrix
diagMatinv = la.fractional_matrix_power(diagMat,-0.5)

lapMat = np.eye(2*N+N2) - np.dot(np.dot(diagMatinv,adjMat),diagMatinv)

# Unnormalized
#lapMat = diagMat - adjMat

#print lapMat

# Cal eigenvector of laplacian
eigval,eigvec = np.linalg.eig(lapMat)

# Sort from smallest to largest
idx = np.argsort(eigval)
eigval = eigval[idx]
eigvec = eigvec[:,idx]

def build_Hk_5(QE_xml_data_file,shift,shift_type,Hk_space,Hk_outfile,nbnds_norm=0,nbnds_in=0):
    """
    returns Hk:
    build_Hk_2: includes all the bands that lay under the 'shift' energy.
    build_Hk_3: Optionally one can inclue a fixed number of bands with 'nbnds_in';
                this capability is similar to WanT's 'atmproj_nbnd'.
    shift_type: 0 = regular shifting. 1 = new shifting

    build_Hk_4: -a bug for nonortho shifting is corrected: Sks was needed for that case.
                -changed the name of the output variable from Hks to Hk.
    build_Hk_5: reads nawf,nkpnts,nspin,shift,eigsmat, from QE_xml_data_file
    """
    nproc = 1

    if not os.path.isfile(QE_xml_data_file):
        sys.exit('File not found: {0:s}'.format(QE_xml_data_file))
    data      = np.load(QE_xml_data_file)
    nawf      = int(data['nawf'])
    nkpnts    = int(data['nkpnts'])
    nspin     = int(data['nspin'])
    eigsmat   = data['eigsmat']
    Sks       = data['Sk']
    U         = data['U']


    if Hk_space.lower()=='ortho':
       del Sks 
    elif Hk_space.lower()=='nonortho':
       if len(Sks.shape) != 3: sys.exit('Need Sks[nawf,nawf,nkpnts]  for nonortho calculations')
    else:
       sys.exit('wrong Hk_space option. Only ortho and nonortho are accepted')

    tic = time.time()
    Hks = np.zeros((nawf,nawf,nkpnts,nspin),dtype=complex)
    for ispin in range(nspin):
        for ik in range(nkpnts):
            my_eigs=eigsmat[:,ik,ispin]
            E = np.diag(my_eigs)
            UU    = U[:,:,ik,ispin] #transpose of U. Now the columns of UU are the eigenvector of length nawf
            if nbnds_norm > 0:
                norms = 1/np.sqrt(np.real(np.sum(np.conj(UU)*UU,axis=0)))
                UU[:,:nbnds_norm] = UU[:,:nbnds_norm]*norms[:nbnds_norm]
            kappa = shift
            if nbnds_in == 0:
               iselect   = np.where(my_eigs <= shift)[0]
            elif nbnds_in > 0:
               iselect   = range(nbnds_in)
            else:
               sys.exit('build_Hk_4: wrong nbnd variable')

            ac    = UU[:,iselect]
            ee1   = E[np.ix_(iselect,iselect)]

            if shift_type ==0:
                Hks_aux = ac.dot(ee1).dot(np.conj(ac).T) + kappa*( -ac.dot(np.conj(ac).T))
            elif shift_type==1:
                aux_p=la.inv(np.dot(np.conj(ac).T,ac))
                Hks_aux = ac.dot(ee1).dot(np.conj(ac).T) + kappa*( -ac.dot(aux_p).dot(np.conj(ac).T))
            else:
                sys.exit('shift_type not recognized')


                Hks_aux = np.triu(Hks_aux,1)+np.diag(np.diag(Hks_aux))+np.conj(np.triu(Hks_aux,1)).T

            if Hk_space.lower()=='ortho':
                Hks[:,:,ik,ispin] = Hks_aux  + kappa*np.identity(nawf)
            elif Hk_space.lower()=='nonortho':
                Sk_half = sla.fractional_matrix_power(Sks[:,:,ik],0.5)
                if 'multi_dot' in dir (la):
                    Hks[:,:,ik,ispin] =la.multi_dot([Sk_half,Hks_aux,Sk_half])+kappa*Sks[:,:,ik]
                else:
                    Hks[:,:,ik,ispin] =np.dot(np.dot(Sk_half,Hks_aux),Sk_half)+kappa*Sks[:,:,ik]

            else:
                sys.exit('wrong Hk_space option. Only ortho and nonortho are accepted')


    
    np.savez(Hk_outfile,Hk=Hks,nbnds_norm=nbnds_norm,nbnds_in=nbnds_in,shift_type=shift_type,shift=shift)
            
    toc = time.time()
    hours, rem = divmod(toc-tic, 3600)
    minutes, seconds = divmod(rem, 60)
    print("Parallel calculation of H[k] with {0:d} processors".format(nproc))
    print("Elapsed time {:0>2}:{:0>2}:{:05.2f}".format(int(hours),int(minutes),seconds))
    return Hks

 def test_opposite_sign_complex_eigenvalues(self):
     M = [[2j, 4], [0, -2j]]
     R = [[1+1j, 2], [0, 1-1j]]
     assert_allclose(np.dot(R, R), M, atol=1e-14)
     assert_allclose(fractional_matrix_power(M, 0.5), R, atol=1e-14)