def FBFMEB(engine,app):
    '''
    This method calculates the energy spectrums of the spin excitations.
    '''
    path,ne=app.path,min(app.ne or engine.nmatrix,engine.nmatrix)
    if path is not None:
        bz,reciprocals=engine.basis.BZ,engine.lattice.reciprocals
        if not isinstance(path,HP.BaseSpace): path=bz.path(HP.KMap(reciprocals,path) if isinstance(path,str) else path,mode='Q')
        result=np.zeros((path.rank(0),ne+1))
        result[:,0]=path.mesh(0) if path.mesh(0).ndim==1 else np.array(range(path.rank(0)))
        engine.log<<'%s: '%path.rank(0)
        for i,paras in enumerate(path('+')):
            engine.log<<'%s%s'%(i,'..' if i<path.rank(0)-1 else '')
            m=engine.matrix(scalefree=app.scalefree,scaleint=app.scaleint,**paras)
            result[i,1:]=nl.eigvalsh(m)[:ne] if app.method=='eigvalsh' else HM.eigsh(m,k=ne,evon=False)
        engine.log<<'\n'
    else:
        result=np.zeros((2,ne+1))
        result[:,0]=np.array(range(2))
        m=engine.matrix(scalefree=app.scalefree,scaleint=app.scaleint)
        result[0,1:]=nl.eigvalsh(m)[:ne] if app.method=='eigvalsh' else HM.eigsh(m,k=ne,evon=False)
        result[1,1:]=result[0,1:]
    name='%s_%s%s'%(engine.tostr(mask=path.tags),app.name,app.suffix)
    if app.savedata: np.savetxt('%s/%s.dat'%(engine.dout,name),result)
    if app.plot: app.figure('L',result,'%s/%s'%(engine.dout,name))
    if app.returndata: return result

def is_psd(M: np.ndarray, strict=False, tol=1e-8) -> bool:
    """
    Tests if matrix is positive semi-definite

    Parameters
    ----------
    M: (N, N) ndarray
        Matrix to be tested for positive semi-definiteness

    strict: bool
        If True, tests for posotive definiteness

    tol: float
        Numeric tolerance to check for equality

    Returns
    -------
    bool
        True if matrix is positive (semi-)definite, else False
    """
    if isinstance(M, (int, float)):
        return M >= 0

    M = np.asarray(M)
    if not is_symmetric(M, tol):
        return False

    if strict:
        return (la.eigvalsh(M) > 0).all()
    else:
        return (la.eigvalsh(M) >= 0).all()

def Eigenvalues(Hgesamt):
	eigenwerte=[]
	eigenwerte2=[]
	for l in range(N/2):
		block=Hgesamt[l*8:8+l*8:,l*8:8+l*8:]
		eigenwerte.append(LA.eigvalsh(block))
		eigenwerte2.extend(LA.eigvalsh(block))
	#print(eigenwerte2)
	return(eigenwerte)

def linear_algebra():
    """ Use the `numpy.linalg` library to do Linear Algebra 
        For a reference on math, see 'Linear Algebra explained in four pages'
        http://minireference.com/static/tutorials/linear_algebra_in_4_pages.pdf
    """

    ### Setup two vectors
    x = np.array([1, 2, 3, 4])
    y = np.array([5, 6, 7, 8])

    ### Vector Operations include addition, subtraction, scaling, norm (length),
    # dot product, and cross product
    print np.vdot(x, y)  # Dot product of two vectors


    ### Setup two arrays / matrices
    a = np.array([[1, 2],
                  [3, 9]])
    b = np.array([[2, 4],
                  [5, 6]])


    ### Dot Product of two arrays
    print np.dot(a, b)


    ### Solving system of equations (i.e. 2 different equations with x and y)
    print LA.solve(a, b)


    ### Inverse of a matrix undoes the effects of the Matrix
    # The matrix multipled by the inverse matrix returns the 
    # 'identity matrix' (ones on the diagonal and zeroes everywhere else); 
    # identity matrix is useful for getting rid of the matrix in some equation
    print LA.inv(a)  # return inverse of the matrix
    print "\n"


    ### Determinant of a matrix is a special way to combine the entries of a
    # matrix that serves to check if matrix is invertible (!=0) or not (=0)
    print LA.det(a)  # returns the determinant of the array
    print "\n"  # e.g. 3, means that it is invertible


    ### Eigenvectors is a special set of input vectors for which the action of
    # the matrix is described as simple 'scaling'.  When a matrix is multiplied
    # by one of its eigenvectors, the output is the same eigenvector multipled
    # by a constant (that constant is the 'eigenvalue' of the matrix)
    print LA.eigvals(a)  # comput the eigenvalues of a general matrix
    print "\n"
    print LA.eigvalsh(a)  # Comput the eigenvalues of a Hermitian or real symmetric matrix
    print "\n"
    print LA.eig(a)  # return the eigenvalues for a square matrix
    print "\n"
    print LA.eigh(a)  # return the eigenvalues or eigenvectors of a Hermitian or symmetric matrix
    print "\n"

    def do(self, a, b):
        # note that eigenvalue arrays returned by eig must be sorted since
        # their order isn't guaranteed.
        ev = linalg.eigvalsh(a, 'L')
        evalues, evectors = linalg.eig(a)
        evalues.sort(axis=-1)
        assert_allclose(ev, evalues, rtol=get_rtol(ev.dtype))

        ev2 = linalg.eigvalsh(a, 'U')
        assert_allclose(ev2, evalues, rtol=get_rtol(ev.dtype))

def info_rate_ild_fast( H_TX, H_CH, H_RX, P_vec, SCHEME, SIM ):
    N0 = SIM['N0'];
    dt = SIM['dt'];
    REAL_DIM_PER_SYM = SIM['REAL_DIM_PER_SYM'];
    T_TRANSMISSION = SIM['T_TRANSMISSION'];
    
    # Power allocation (uniform)
    Sigma_X_NORMALIZED, layer = power_alloc(H_TX, SCHEME, SIM);
    K_prime = len(layer); #K_prime = SHCEME['K_prime']
    
    # Compute mutual informal layer by layer
    R_vec_per_layer = zeros((K_prime,len(P_vec)));
    for k in range(K_prime):
        sigma2_N = N0/2 * REAL_DIM_PER_SYM * (1/dt);  # DEPENDS ON WHETHER NOISE IS REAL (PASSBAND) OR COMPLEX (BASEBAND)
        #Sigma_N = sigma2_N * eye(length(layer{k}));
        
        H_RX_k = H_RX[:,layer[k]];
        H_RX_k, s, V = svd(H_RX_k); # ensure that H_RX_k is orthonormal
        H_RXCHTX = H_RX_k.transpose().dot(H_CH).dot(H_TX);
        
        
        Sigma_X = Sigma_X_NORMALIZED;     #XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
        
        Sigma_S = H_RXCHTX.dot(Sigma_X).dot(H_RXCHTX.transpose());
        
        #Sigma_Z = H_RX_k.transpose().dot(sigma2_N).dot(H_RX_k);
        
        Sigma_X_bar = Sigma_X;
        Sigma_X_bar[:,layer[k]] = 0;
        Sigma_X_bar[layer[k],:] = 0;
        
        Sigma_S_bar = H_RXCHTX.dot(Sigma_X_bar).dot(H_RXCHTX.transpose());
    
        EIG_S = eigvalsh(Sigma_S);
        EIG_S_bar = eigvalsh(Sigma_S_bar);
        
        I_vec = zeros(len(P_vec));
        loop_index = 0;
        for P in P_vec:
            #I = sum(log(eig( P*Sigma_S+Sigma_Z ))) - sum(log(eig( P*Sigma_S_bar+Sigma_Z )));
            I = sum(log( P * EIG_S + sigma2_N )) - sum(log( P * EIG_S_bar + sigma2_N ));
            I_vec[loop_index] = I;
            loop_index = loop_index + 1;
        
        R_vec_per_layer[k,:] = I_vec / 2 * REAL_DIM_PER_SYM / T_TRANSMISSION;

    
    R_vec = sum(R_vec_per_layer,0); # sum the rows (sum along the columns)
    
    return R_vec

def info_rate(G_PRE,H,Sigma_Z): 
    #G_PRE is a list for precoding matrices for all users
    #H is a list for channel matrices
    G1, G2 = G_PRE
    
    R_users = list()
    for u in [0,1]: # 2 users
        m, n = H[u].shape
        R_total = array(0)
        for i in range(n):
            #G = G_RX[u]
            
            #R_v1 = identity(n)
            #R_v2 = identity(n)
            R_v1 = zeros((n,n))
            R_v1[0,0] = 1
            R_v2 = zeros((n,n))
            R_v2[1,1] = 1

            Sigma_X = \
            G1.dot(R_v1).dot(G1.conj().transpose()) + \
            G2.dot(R_v2).dot(G2.conj().transpose())
            
            Sigma_Y = ( 
            H[u].dot(Sigma_X).dot(H[u].conj().transpose()) + 
            Sigma_Z 
            )
            
            if   u==0: # if user 1
                R_v1[i,i] = 0
            elif u==1: # if user 2
                R_v2[i,i] = 0

            Sigma_X_V = \
            G1.dot(R_v1).dot(G1.conj().transpose()) + \
            G2.dot(R_v2).dot(G2.conj().transpose())
            
            Sigma_Y_V = ( 
            H[u].dot(Sigma_X_V).dot(H[u].conj().transpose()) + 
            Sigma_Z 
            )
            
            H_Y = sum(log(eigvalsh(Sigma_Y)))
            H_Y_V = sum(log(eigvalsh(Sigma_Y_V)))
            R = H_Y - H_Y_V    
    
            R_total = R_total + R
            
        R_users.append(R_total)          
    return R_users

    def _fa_paf(self, smc = True):
        cor = self.cor.copy()
        if smc:
            np.fill_diagonal(cor, utils.SMC(cor))
        eig_val1 = la.eigvalsh(self.cor)
        eig_val1.sort()
        comm = np.sum(np.diag(cor))
        err = comm
        comm_list = []
        i = 0
        while err > self.lower:
            eig_val, eig_vec = utils.eigenh_sorted(cor)
            if (self.n_factors > 1):
                loadings = eig_vec[:, :self.n_factors].dot(
                    np.diag(np.sqrt(eig_val[:self.n_factors])))
            else:
                loadings = eig_vec[:, 0] * np.sqrt(eig_val[0])

            model = loadings.dot(loadings.T)
            new = np.diag(model)
            comm1 = np.sum(new)
            np.fill_diagonal(cor, new)
            err = np.abs(comm - comm1)
            if np.iscomplex(err):
                print('Imaginary eigenvalue condition'
                      ' occurred!')
                break
            comm = comm1
            comm_list.append(comm1)
            i += 1
            if i > 1000:
                print('maximum iteration exceeded')
                err = 0

        self.loadings = loadings

def eigenvalues_hadamard(matrix1, matrix2):
    """Computes the Hadamard product of 2 matrices. See
    https://www.johndcook.com/blog/2018/10/10/hadamard-product/ for details

    :param matrix1: first matrix
    :param matrix2: second matrix
    :return: lower and upper
    """

    matrix1 = array(matrix1)  # as arrays
    matrix2 = array(matrix2)

    eig_a = eigvalsh(matrix1)  # eigenvalues (optimized for Hermitian matrices)
    eig_b = eigvalsh(matrix2)

    return eig_a, eig_b

def normal_eval(mu, P, x, dP=None):
    """
    Probability of x under normal(mu,inv(P))

    Parameters
    ----------
    mu: array of shape (n): the mean parameter
    P: array of shape (n,n): the precision matrix 
    x: array of shape (n): the data to be evaluated

    Returns
    -------
    (float) the density
    """
    p = np.size(mu)
    if dP==None:
        dP = np.prod(eigvalsh(P))
    mu = np.reshape(mu,(1,p))
    w0 = np.log(dP)-p*np.log(2*np.pi)
    w0 /= 2               
    x = np.reshape(x,(1,p))
    q = np.dot(np.dot(mu-x,P),(mu-x).T)
    w = w0 - q/2
    L = np.exp(w)
    return np.squeeze(L)

def compactness(im):
    from numpy import array, meshgrid, arange, shape, mean, zeros
    from numpy import outer, sum, max, linalg
    from numpy import sqrt
    from math import exp

    (h, w) = shape(im)
    (X, Y) = meshgrid(arange(w), arange(h))
    x = X.flatten()
    y = Y.flatten()
    wgts = im[y, x]
    sw = sum(wgts)
    if sw == 0:
        return 1
    wgts /= sw
    wpts = array([wgts * x, wgts * y])
    wmean = sum(wpts, 1)
    N = len(x)
    s = array([x, y])
    P = zeros((2, 2))
    for i in range(0, N):
        P += wgts[i] * outer(s[:, i], s[:, i])
    P = P - outer(wmean, wmean)

    det = abs(linalg.det(P))
    if det > 0:
        v = linalg.eigvalsh(P)
        v = abs(v)
        r = min(v) / max(v)
        return 100.0 * sqrt(r / det)
    else:
        return 0.0

def calculate_gyration_tensor_parameters(points):
    """
    Calculates the gyration tensor parameters R_g^2, η, c, κ from a list of
    all points inside a cavity.
     - R_g^2 is the squared gyration radius
     - η is the asphericity
     - c is the acylindricity
     - κ is the anisotropy
    """

    points = np.array(points, dtype=np.float)
    mean = np.mean(points, axis=0)
    points -= mean

    gyration_tensor = np.zeros((3, 3))
    for i in range(3):
        for j in range(i, 3):
            gyration_tensor[i, j] = np.dot(points[:, i], points[:, j])
            gyration_tensor[j, i] = gyration_tensor[i, j]
    # cell volume is constant, cavity volume is proportional to len(points)
    gyration_tensor /= len(points)

    eigvals = list(sorted(la.eigvalsh(gyration_tensor), reverse=True))

    squared_gyration_radius = sum(eigvals)
    if squared_gyration_radius > 0:
        asphericity = (eigvals[0] - 0.5 * (eigvals[1] + eigvals[2])) / squared_gyration_radius
        acylindricity = (eigvals[1] - eigvals[2]) / squared_gyration_radius
        anisotropy = (asphericity ** 2 + 0.75 * acylindricity ** 2) ** 0.5
    else:
        asphericity = 0
        acylindricity = 0
        anisotropy = 0
    return mean, squared_gyration_radius, asphericity, acylindricity, anisotropy

def array_compactness(im):
        '''
        calculate the compactness of a 2D array. Each element of the 2D array
        should be proportional to the score of that pixel in the overall scoring scheme
        . '''
        from numpy import array,meshgrid,arange,shape,mean,zeros
        from numpy import outer,sum,max,linalg
        from numpy import sqrt
        from math import exp
        (h,w) = shape(im)
        (X,Y) = meshgrid(arange(w),arange(h))
        x = X.flatten()
        y = Y.flatten()
        wgts = im[y,x]
        sw = sum(wgts)
        if sw == 0:
                return 1
        wgts /= sw
        wpts = array([wgts*x, wgts*y])
        wmean = sum(wpts, 1)
        N = len(x)
        s = array([x,y])
        P = zeros((2,2))
        for i in range(0,N):
                P += wgts[i]*outer(s[:,i],s[:,i])
        P = P - outer(wmean,wmean);

        det = abs(linalg.det(P))
        if (det <= 0):
                return 0.0
        v = linalg.eigvalsh(P)
        v = abs(v)
        r = min(v)/max(v)
        return 100.0*sqrt(r/det)

def get_vn_entropy(psi, spin, N=None, mode='1spin', base='e'):
    """
    Compute the von Neumann entropy for a given state "psi".
    "psi" must be a column vector. Sparsity is optional.
    Available modes: "1spin" for the entanglement entropy for only the first
    spin, "eqsplit" for the entanglement entropy for an evenly split system.
    Available bases: "e" for natural log, "2" for log2 and "10" for log10.
    """
    if N is None and mode == 'eqsplit':
        raise Exception("N cannot be 'None' for mode='eqsplit'.")

    if mode == '1spin':
        red_rho_A = red_rho_A_1spin(psi, spin)
    elif mode == 'eqsplit':
        red_rho_A = red_rho_eqsplit(psi, spin, N)

    lamb = eigvalsh(red_rho_A.todense())  # Eigenvalues of the reduced matrix.
    S_AB_terms = []
    for i in range(np.shape(red_rho_A)[0]):
        if abs(lamb[i]) < 1e-6:
            # lim a->0 (alog(a)) = 0. It also removes some minuscule negative
            #  lambda values resulting from rounding errors.
            S_AB_terms.append(0)
        else:
            if base == 'e':
                S_AB_terms.append(-lamb[i] * np.log(lamb[i]))
            elif base == '2' or base == 2:
                S_AB_terms.append(-lamb[i] * np.log2(lamb[i]))
            elif base == '10' or base == 10:
                S_AB_terms.append(-lamb[i] * np.log10(lamb[i]))
            else:
                raise Exception('Available bases are "e", "2" and "10"')

    return float(np.sum(S_AB_terms))

def main():

    L = 42
    t = 1
    h_0 = t*np.eye(L, k=1)                                      ##Upper triagonal of hopping in y-direction
    h_0[0,L-1] = np.conjugate(t)                                ##Periodic boundary condition in -y-direction
    h_0 = h_0 + h_0.conj().T                                    ##Full matrix for one slice in y
    h_0 = np.kron(np.eye(L), h_0)                               ##Put all slices together

    diag = np.eye(L, k=1)
    period = np.eye(L, k=L-1)

##Now only the hopping between 1d slices is missing the for loop will add these

#    ps = np.array([1, 2, 6, 8, 12, 14, 16, 21])
    ps = np.linspace(0, 100, 50)                               ##All fluxes
    spectrum = []

    for p in np.nditer(ps):
        t_mag = np.exp(2j*np.pi*np.arange(L)*p/L)               ##Peierls phase for all fluxes
        X = t*np.diag(t_mag)
        X = np.kron(diag, X) + np.kron(period, X.conj().T)
        X = X + X.conj().T                                      ##Construct the hopping matrix in x-direction between 1d slices
        H = h_0 + X                                             ##Add the two hamiltonian parts
        spectrum.append(LA.eigvalsh(H))                         ##Get eigenvalues for each flux

    ax = plt.subplot(111)
    ps = ps/L
    ax.plot(ps, spectrum, marker='.', markersize=0.1, linestyle='None', color='k')

    plt.show()

  def __call__(self, function, point, state):
    """
    Computes Goldfeld step 
    """
    g = function.gradient(point)
    state['gradient'] = g
    G = function.hessian(point)
    state['hessian'] = G
    c = 1e-8 # is this one best?
    
    d0 = None

    try:
        L = cholesky(G)
        # reach here => isPositiveDefinite = True
        step = n_solve(L.T, n_solve(L, -g))
    except:
        # isPositiveDefinite = False
        G_eigvals = eigvalsh(G)
        minEig = min(G_eigvals)
        if minEig < 0:
            shift = -minEig + c
            
            #avoiding sparse case with big nVars
            for i in xrange(point):  G[i,i] += shift
                
        step = n_solve(G, -g)

    state['direction'] = step
    return step

    def unweighted_likelihood_(self, x):
        """
        return the likelihood of each data for each component
        the values are not weighted by the component weights

        Parameters
        ----------
        x: array of shape (n_samples,self.dim)
           the data used in the estimation process

        Returns
        -------
        like, array of shape(n_samples,self.k)
          unweighted component-wise likelihood
        """
        n = x.shape[0]
        like = np.zeros((n, self.k))

        for k in range(self.k):
            # compute the data-independent factor first
            w = -np.log(2 * np.pi) * self.dim
            m = np.reshape(self.means[k], (1, self.dim))
            b = self.precisions[k]
            if self.prec_type == "full":
                w += np.log(eigvalsh(b)).sum()
                dx = m - x
                q = np.sum(np.dot(dx, b) * dx, 1)
            else:
                w += np.sum(np.log(b))
                q = np.dot((m - x) ** 2, b)
            w -= q
            w /= 2
            like[:, k] = np.exp(w)
        return like

  def f(k):
    kp = kpgen(k) # get kpoint
    hk = hk_gen(kp) # generate hamiltonian
#    es,ew = lgs.eigsh(csc_matrix(hk),k=4,which="LM",sigma=0.0)
    es = lg.eigvalsh(hk) # get eigenvalues
    g = np.min(es[es>0.]) - np.max(es[es<0.])
    return g  # return gap

def genMulCov(size, numberOfCov, low, upper, mode, portion = 0.05):
    S_set = []   
    Cov_set = []
    minEVal_set = [] 
#    low = abs(low)
#    upper = abs(upper)
    m = size/3
    mm = m/2
#    print m, mm
    S_init = np.zeros((size,size))
    for k in range(numberOfCov):
        S = np.zeros((size,size))
        if k == 0:
            S = genInvCov(size, low, upper, portion)
            if mode == 5:      
                ind_zero = np.where(spy.sparse.rand(m, size-m, 0.5).todense() == 0)
                value = np.multiply((np.random.randint(2, size = (m, size -m)) - 0.5)*2,(low + (upper - low)*np.random.rand(m,size -m)))
                value[ind_zero] = 0
                hub = value
                S[:m, m:size] = hub
                S[m:size, :m] = hub.T        
                minEVal_set.append(alg.eigvalsh(S)[0])
            S_init = S
        elif mode == 3: #'laplacian'
            ind1 = range(m)
            ind2 = np.random.permutation(m)
            S = np.copy(S_init)
            S[ind1, :] = S[ind2, :]            
            S[:, ind1] = S[:, ind2]
        elif mode  == 5: #'perturbation'
            S = np.copy(S_init)
            ind_zero = np.where(spy.sparse.rand(mm, size-mm, 0.5).todense() == 0)
            pert = np.multiply((np.random.randint(2, size = (mm, size -mm)) - 0.5)*2,(low + (upper - low)*np.random.rand(mm,size -mm)))
            pert[ind_zero] = 0 
            S[:mm, mm:size] = pert
            S[mm:size, :mm] = pert.T
            minEVal_set.append(alg.eigvalsh(S)[0])
        else:
#            print 'Activate normal mode'
            S = genInvCov(size, low, upper, portion)
        S_set.append(S)
    
    for k in range(numberOfCov):
        if mode == 5:
            S_set[k] = S_set[k] + (0.1 - min(minEVal_set))*np.identity(size)
        Cov_set.append(alg.inv(S_set[k]))
    return S_set, Cov_set   

 def do(self, a, b):
     # note that eigenvalue arrays must be sorted since
     # their order isn't guaranteed.
     ev = linalg.eigvalsh(a)
     evalues, evectors = linalg.eig(a)
     ev.sort(axis=-1)
     evalues.sort(axis=-1)
     assert_almost_equal(ev, evalues)