def test_solver_parameters(self):
        A = poisson((50, 50), format='csr')

        for method in methods:
            # method = ('richardson', {'omega':4.0/3.0})
            ml = smoothed_aggregation_solver(A, presmoother=method,
                                             postsmoother=method,
                                             max_coarse=10)

            residuals = profile_solver(ml)
            assert((residuals[-1]/residuals[0])**(1.0/len(residuals)) < 0.95)
            assert(ml.symmetric_smoothing)

        for method in methods2:
            ml = smoothed_aggregation_solver(A, max_coarse=10)
            change_smoothers(ml, presmoother=method[0], postsmoother=method[1])

            residuals = profile_solver(ml)
            assert((residuals[-1]/residuals[0])**(1.0/len(residuals)) < 0.95)
            assert(not ml.symmetric_smoothing)

        for method in methods3:
            ml = smoothed_aggregation_solver(A, max_coarse=10)
            change_smoothers(ml, presmoother=method[0], postsmoother=method[1])
            assert(ml.symmetric_smoothing)

        for method in methods4:
            ml = smoothed_aggregation_solver(A, max_coarse=10)
            change_smoothers(ml, presmoother=method[0], postsmoother=method[1])
            assert(not ml.symmetric_smoothing)

def solveMG(A, b):
    ml = pyamg.smoothed_aggregation_solver(A)

    x = np.zeros(b.shape)
    for bi in range(3):
        x[:, bi] = ml.solve(b[:, bi], tol=1e-10)
    return x

def random_walker(mask, prior, gamma):
    """
    Assume prior is given on the mask, of shape (NPTS, K).
    Return random walker probability map.
    """
    gc.enable()

    # Assembling graph Laplacian
    L = make_laplacian(mask)
    n = L.shape[0]
    L = L + sparse.coo_matrix((gamma * prior.sum(axis=1),
                               (range(n), range(n))))

    # Creating sparse solver
    mls = smoothed_aggregation_solver(L.tocsr())
    del L
    gc.collect()

    # Loop over classes
    X = []
    for k in range(prior.shape[-1]):
        X += [mls.solve(gamma * prior[:, k])]

    del mls
    gc.collect()
    return np.array(X).T

def test():
    class EnergyEdgeKernel(object):
        def __init__(self):
            self.subdomains = [None]
            return

        def eval(self, mesh, cell_mask):
            edge_ce_ratio = mesh.ce_ratios[..., cell_mask]
            beta = 1.0
            return numpy.array(
                [
                    [edge_ce_ratio, -edge_ce_ratio * numpy.exp(1j * beta)],
                    [-edge_ce_ratio * numpy.exp(-1j * beta), edge_ce_ratio],
                ]
            )

    vertices, cells = meshzoo.rectangle(0.0, 2.0, 0.0, 1.0, 101, 51)
    mesh = meshplex.MeshTri(vertices, cells)

    matrix = pyfvm.get_fvm_matrix(mesh, [EnergyEdgeKernel()], [], [], [])
    rhs = mesh.control_volumes.copy()

    sa = pyamg.smoothed_aggregation_solver(matrix, smooth="energy")
    u = sa.solve(rhs, tol=1e-10)

    # Cannot write complex data ot VTU; split real and imaginary parts first.
    # <http://stackoverflow.com/a/38902227/353337>
    mesh.write("out.vtk", point_data={"u": u.view("(2,)float")})
    return

def edgeAMG(Anode, Acurl, D):
    nodalAMG = pyamg.smoothed_aggregation_solver(Anode, max_coarse=10, keep=True)

    # construct multilevel structure
    levels = []
    levels.append(pyamg.multilevel_solver.level())
    levels[-1].A = Acurl
    levels[-1].D = D
    for i in range(1, len(nodalAMG.levels)):
        A = levels[-1].A
        Pnode = nodalAMG.levels[i - 1].AggOp
        P = findPEdge(D, Pnode)
        R = P.T
        levels[-1].P = P
        levels[-1].R = R
        levels.append(pyamg.multilevel_solver.level())
        A = R * A * P
        M = sparse.dia_matrix((1.0 / ((P.T * P).diagonal()), 0),
                              shape=(P.shape[1], P.shape[1]))
        D = M * (P.T * D * Pnode)
        D = D.tocsr()
        levels[-1].A = A
        levels[-1].D = D

    edgeML = pyamg.multilevel_solver(levels)
    for i in range(0, len(edgeML.levels)):
        edgeML.levels[i].presmoother = setup_hiptmair(levels[i])
        edgeML.levels[i].postsmoother = setup_hiptmair(levels[i])
    return edgeML

def random_walker(mask, prior, gamma):
    """
    Assume prior is given on the mask, of shape (NPTS, K)
    """
    gc.enable()

    print('Assembling graph Laplacian...')
    L = make_laplacian(mask)
    n = L.shape[0]
    L = L + sparse.coo_matrix((gamma * prior.sum(axis=1),
                               (range(n), range(n))))

    print('Creating sparse solver...')
    mls = smoothed_aggregation_solver(L.tocsr())
    del L
    gc.collect()

    print('Loop over classes...')
    X = []
    for k in range(prior.shape[-1]):
        print('  Doing class %d...' % k)
        X += [mls.solve(gamma * prior[:, k])]

    del mls
    gc.collect()
    return np.array(X).T

def laplacian_basis(W, k, largest = False, method = "arpack"):
    """Build laplacian basis matrix with k bases from weighted adjacency matrix W."""

    logger.info(
        "solving for %i %s eigenvector(s) of the Laplacian (using %s)",
        k,
        "largest" if largest else "smallest",
        method,
        )

    L = laplacian_operator(W) 

    assert isinstance(L, scipy.sparse.csr_matrix)
    assert k > 0
    assert k < L.shape[0]

    if method == "amg":
        solver = pyamg.smoothed_aggregation_solver(L)
        pre = solver.aspreconditioner()
        initial = scipy.rand(L.shape[0], k)

        (evals, basis) = scipy.sparse.linalg.lobpcg(L, initial, M = pre, tol = 1e-10, largest = largest)
        logger.info('amg eigen values: '+str(evals))
    elif method == "arpack":
        logger.info('using arpack')
        if largest:
            which = "LM"
        else:
            which = "SM"

        if hasattr(scipy.sparse.linalg, "eigsh"): # check scipy version
            which = "LM" # use sigma=0 and ask for the large eigenvalues (shift trick, see arpack doc)
            (evals, basis) = scipy.sparse.linalg.eigsh(L, k, which = which, tol=1e-10, sigma = 0, maxiter=15000)
            try:
                for i in xrange(len(basis)):
                    b = basis[i]
                    print 'basis vector shape: ', b.shape
                    residual = np.linalg.norm((np.dot(L,b)).todense()-evals[i]*b)
                    perp_test = np.linalg.norm(np.dot(basis[i],basis[1]))
                    logger.info('eigenvalue residual: %f',residual)
                    logger.info('dot of %ith eigenvector with first: %f',i,perp_test)

            except:
                print 'error in eigensolver test code'

            logger.info('arpack eigen values: '+str(evals))
        else: 
            (evals, basis) = scipy.sparse.linalg.eigen_symmetric(L, k, which = which)
            logger.info('arpack (old) eigen values: '+str(evals))
    elif method == "dense":
        (evals, full_basis) = np.linalg.eigh(L.todense())

        basis = full_basis[:, :k]
    else:
        raise ValueError("unrecognized eigenvector method name")

    assert basis.shape[1] == k

    return basis

def solver(A, config):
    """
    Given a matrix A and a solver configuration dictionary, generate a
    smoothed_aggregation_solver

    Parameters
    ----------
    A : {array, matrix, csr_matrix, bsr_matrix}
        Matrix to invert, CSR or BSR format preferred for efficiency
    config : {dict}
        A dictionary of solver configuration parameters that is used to
        generate a smoothed aggregation solver

    Returns
    -------
    ml : {smoothed_aggregation_solver}
        smoothed aggregation hierarchy

    Notes
    -----
    config must contain the following parameter entries for
    smoothed_aggregation_solver:
        symmetry, smooth, presmoother, postsmoother, B, strength,
        max_levels, max_coarse, coarse_solver, aggregate, keep

    Examples
    --------
    >>> from pyamg.gallery import poisson
    >>> from pyamg import solver_configuration,solver
    >>> A = poisson((40,40),format='csr')
    >>> config = solver_configuration(A,verb=False)
    >>> ml = solver(A,config)

    """

    # Convert A to acceptable format
    A = make_csr(A)

    # Generate smoothed aggregation solver
    try:
        return smoothed_aggregation_solver(
            A,
            B=config["B"],
            BH=config["BH"],
            smooth=config["smooth"],
            strength=config["strength"],
            max_levels=config["max_levels"],
            max_coarse=config["max_coarse"],
            coarse_solver=config["coarse_solver"],
            symmetry=config["symmetry"],
            aggregate=config["aggregate"],
            presmoother=config["presmoother"],
            postsmoother=config["postsmoother"],
            keep=config["keep"],
        )
    except:
        raise TypeError("Failed generating smoothed_aggregation_solver")

def MGsetup(nx):
    import numpy as np
    import scipy as sp
    import scipy.sparse
    import pyamg
    import scipy.io
    # scipy.io.savemat('A.mat', {'A': A})
    A = scipy.io.loadmat('A.mat')['A'].tocsr()
    ml = pyamg.smoothed_aggregation_solver(A, max_coarse=10)
    b = np.random.rand(A.shape[0])

def solve(A, b, tol=1e-8):
    ml = pyamg.smoothed_aggregation_solver(A)

    if len(b.shape) == 1:
        x = ml.solve(b, tol=tol)
        return x

    x = np.zeros(b.shape)
    for bi in xrange(b.shape[1]):
        x[:, bi] = ml.solve(b[:, bi], tol=tol)
    return x

    def _linear_analysis(self, method, **kwargs):
        """Performs the linear analysis, in which the pressure and flow fields
        are computed.
        INPUT: method: This can be either 'direct' or 'iterative'
               **kwargs
               precision: The accuracy to which the ls is to be solved. If not
                          supplied, machine accuracy will be used. (This only
                          applies to the iterative solver)
        OUTPUT: The maximum, mean, and median pressure change. Moreover,
                pressure and flow are modified in-place.
        """

        G = self._G
        A = self._A.tocsr()
        if method == 'direct':
            linalg.use_solver(useUmfpack=True)
            x = linalg.spsolve(A, self._b)
        elif method == 'iterative':
            if kwargs.has_key('precision'):
                eps = kwargs['precision']
            else:
                eps = self._eps
            AA = smoothed_aggregation_solver(A, max_levels=10, max_coarse=500)
            x = abs(AA.solve(self._b, x0=None, tol=eps, accel='cg', cycle='V', maxiter=150))
            # abs required, as (small) negative pressures may arise
        elif method == 'iterative2':
         # Set linear solver
             ml = rootnode_solver(A, smooth=('energy', {'degree':2}), strength='evolution' )
             M = ml.aspreconditioner(cycle='V')
             # Solve pressure system
             #x,info = gmres(A, self._b, tol=self._eps, maxiter=50, M=M, x0=self._x)
             #x,info = gmres(A, self._b, tol=self._eps/10000000000000, maxiter=50, M=M)
             x,info = gmres(A, self._b, tol=self._eps/10000, maxiter=50, M=M)
             if info != 0:
                 print('SOLVEERROR in Solving the Matrix')

        pdiff = map(abs, [(p - xx) / p if p > 0 else 0.0
                          for p, xx in zip(G.vs['pressure'], x)])
        maxPDiff = max(pdiff)
        meanPDiff = np.mean(pdiff)
        medianPDiff = np.median(pdiff)
        log.debug(np.nonzero(np.array(pdiff) == maxPDiff)[0])

        G.vs['pressure'] = x
        G.es['flow'] = [abs(G.vs[edge.source]['pressure'] -   \
                            G.vs[edge.target]['pressure']) *  \
                        edge['conductance'] for edge in G.es]

	self._maxPDiff=maxPDiff
        self._meanPDiff=meanPDiff
        self._medianPDiff=medianPDiff

        return maxPDiff, meanPDiff, medianPDiff

    def keo_amg(self, psi):
        """Algebraic multigrid solve.
        """
        import pyamg

        if self._keo_amg_solver is None:
            if self._modeleval._keo is None:
                self._modeleval._assemble_keo()
            self._keo_amg_solver = pyamg.smoothed_aggregation_solver(
                self._modeleval._keo
            )
        return self._keo_amg_solver.solve(psi, tol=1e-12, accel=None)

    def test_solver_parameters(self):
        A = poisson((50,50), format='csr')

        for method in methods:
            #method = ('richardson', {'omega':4.0/3.0})
            ml = smoothed_aggregation_solver(A, presmoother=method, postsmoother=method, max_coarse=10)

            residuals = profile_solver(ml)
            #print "method",method
            #print "residuals",residuals
            #print "convergence rate:",(residuals[-1]/residuals[0])**(1.0/len(residuals))
            assert( (residuals[-1]/residuals[0])**(1.0/len(residuals)) < 0.95 )

        for method in methods2:
            ml = smoothed_aggregation_solver(A, max_coarse=10)
            change_smoothers(ml, presmoother=method[0], postsmoother=method[1])

            residuals = profile_solver(ml)
            #print "method",method
            #print "residuals",residuals
            #print "convergence rate:",(residuals[-1]/residuals[0])**(1.0/len(residuals))
            assert( (residuals[-1]/residuals[0])**(1.0/len(residuals)) < 0.95 )

  def solve_on_coarse_level(self):

    if comm.rank == 0:
      if self.verbosity >= 2:
        print pid+"  Solving on coarse level"

      timer = Timer("Coarse level solution")

      if self.problem.switch_matrices_on_coarse_level:
        A = self.B_coarse
        B = self.A_coarse
        largest = True
        which = 'LM'
      else:
        A = self.A_coarse
        B = self.B_coarse
        largest = False
        which = 'SM'

      # Set initial approximation
      self.v_coarse.fill(0.0)
      self.v_coarse[0] = 1.0

      if self.use_lobpcg_on_coarse_level:
        if self.precond_lobpcg_by_ml:
          if self.update_lobpcg_prec or self.M is None:
            if self.verbosity >= 3:
              print0(pid+"    Creating coarse level preconditioner")

            ml = smoothed_aggregation_solver(A)
            self.M = ml.aspreconditioner()

        w, v, h = lobpcg(A, self.v_coarse, B, self.M, tol=self.coarse_level_tol, maxiter=self.coarse_level_maxit,
                         largest=largest, verbosityLevel=self.lobpcg_verb, retResidualNormsHistory=True)
      else:
        if self.problem.sym:
          w, v = eigsh(A, 1, B, which=which, v0=self.v_coarse,
                       ncv=self.coarse_level_num_ritz_vec, maxiter=self.coarse_level_maxit, tol=self.coarse_level_tol)
        else:
          w, v = eigs(A, 1, B, which=which, v0=self.v_coarse,
                      ncv=self.coarse_level_num_ritz_vec, maxiter=self.coarse_level_maxit, tol=self.coarse_level_tol)

      self.lam = w[0]
      self.v_coarse = v[0]

      try:
        self.num_it_coarse += len(h)
      except NameError:
        pass  # There seems to be no way to obtain number of iterations for eigs/eigsh

def fiedler(adj_list,plot=False,fn="FiedlerPlots",n_fied=2):
    """calculate the first fiedler vector of a graph adjascancy list and optionally write associated plots to file.

    Takes:
    adj_list:
    An Nx2 nested list of ints of the form:
    [[node1,node2],
    ...]
    Representing the adjascancy list.

    plot=False: make plots or not.
    fn="FiedlerPlots": filename to prepend to the plot png file names
    n_fied=2: the number of fiedler vectors to calculate (values above 2 will not be output)

    Returns a Dictionary of the form:



    {"f1": the first fiedler vector,
    "f2": (if caclulated) the second fideler vector
    "d": the node degrees,
    "r1": the rank of each node in the first fiedler vector
    "r2": the rank of each node in the second fiedler vector}


    """
    

    A = graph_laplacian(adj_list)

    # construct preconditioner
    ml = smoothed_aggregation_solver(A, coarse_solver='pinv2',max_coarse=10)
    M = ml.aspreconditioner()

    # solve for lowest two modes: constant vector and Fiedler vector
    X = scipy.rand(A.shape[0], n_fied+1)
    (eval,evec,res) = lobpcg(A, X, M=None, tol=1e-12, largest=False, \
            verbosityLevel=0, retResidualNormsHistory=True)

    if plot:
        doPlots(evec[:,1],evec[:,2],A.diagonal(),adj_list,fn)
        
        
    out = {"f1":list(evec[:,1]),"d":list(A.diagonal()),"r1":[int(i) for i in list(numpy.argsort(numpy.argsort(evec[:,1])))]}
    if n_fied > 1:
        out["f2"]=list(evec[:,2])
        out["r2"]=[int(i) for i in list(numpy.argsort(numpy.argsort(evec[:,2])))]
    return out

	def __init__(self, p):
		"""Parses the file _infoSolver to initialize the solver.
		
		Arguments
		---------
		p -- field variable (Variable object) of the Poisson equation.
		"""
		with open(Case.path+'/_infoSolver.yaml', 'r') as infile:
			info = yaml.load(infile)['poisson']
		self.solver = Poisson.solvers[info['solver']]
		self.tol = info['tol']
		self.maxiter = info['maxiter']
		self.M = None
		if 'precond' in info:
			ml = pyamg.smoothed_aggregation_solver(p.laplacian.mat)
			self.M = ml.aspreconditioner(cycle=info['precond']['cycle'])
		self.iterations, self.residuals = [], []

def locally_linear_embedding(X, n_neighbors, out_dim, tol=1e-6, max_iter=200):
    W = neighbors.kneighbors_graph(X, n_neighbors=n_neighbors, mode="barycenter")

    # M = (I-W)' (I-W)
    A = eye(*W.shape, format=W.format) - W
    A = (A.T).dot(A).tocsr()

    # initial approximation to the eigenvectors
    X = np.random.rand(W.shape[0], out_dim)
    ml = smoothed_aggregation_solver(A, symmetry="symmetric")
    prec = ml.aspreconditioner()

    # compute eigenvalues and eigenvectors with LOBPCG
    eigen_values, eigen_vectors = linalg.lobpcg(A, X, M=prec, largest=False, tol=tol, maxiter=max_iter)

    index = np.argsort(eigen_values)
    return eigen_vectors[:, index], np.sum(eigen_values)

def solver_diagnostic(A):
    ##
    # Generate B
    B = ones((A.shape[0],1), dtype=A.dtype); BH = B.copy()

    ##
    # Random initial guess, zero right-hand side
    random.seed(0)
    b = zeros((A.shape[0],1))
    x0 = rand(A.shape[0],1)

    ##
    # Create solver
    ml = smoothed_aggregation_solver(A, B=B, BH=BH,
        strength=('symmetric', {'theta': 0.0}),
        smooth=('energy', {'weighting': 'local', 'krylov': 'gmres', 'degree': 1, 'maxiter': 2}),
        improve_candidates=[('gauss_seidel_nr', {'sweep': 'symmetric', 'iterations': 4}), None],
        aggregate="standard",
        presmoother=('gauss_seidel_nr', {'sweep': 'symmetric', 'iterations': 2}),
        postsmoother=('gauss_seidel_nr', {'sweep': 'symmetric', 'iterations': 2}),
        max_levels=15,
        max_coarse=300,
        coarse_solver="pinv")

    ##
    # Solve system
    res = []
    x = ml.solve(b, x0=x0, tol=1e-08, residuals=res, accel="gmres", maxiter=300, cycle="V")
    res_rate = (res[-1]/res[0])**(1.0/(len(res)-1.))
    normr0 = norm(ravel(b) - ravel(A*x0))
    print " "
    print ml
    print "System size:                " + str(A.shape)
    print "Avg. Resid Reduction:       %1.2f"%res_rate
    print "Iterations:                 %d"%len(res)
    print "Operator Complexity:        %1.2f"%ml.operator_complexity()
    print "Work per DOA:               %1.2f"%(ml.cycle_complexity()/abs(log10(res_rate)))
    print "Relative residual norm:     %1.2e"%(norm(ravel(b) - ravel(A*x))/normr0)

    ##
    # Plot residual history
    pylab.semilogy(array(res)/normr0)
    pylab.title('Residual Histories')
    pylab.xlabel('Iteration')
    pylab.ylabel('Relative Residual Norm')
    pylab.show()

def spectral_partition(A):
    ml = smoothed_aggregation_solver(A,
            coarse_solver='pinv2',max_coarse=100,smooth=None, strength=None)
    print ml

    M = ml.aspreconditioner()

    X = sp.rand(A.shape[0], 2) 
    (evals,evecs,res) = lobpcg(A, X, M=M, tol=1e-12, largest=False, \
        verbosityLevel=0, retResidualNormsHistory=True, maxiter=200)

    fiedler = evecs[:,1]
    vmed = np.median(fiedler)
    v = np.zeros((A.shape[0],))
    K = np.where(fiedler<=vmed)[0]
    v[K]=-1
    K = np.where(fiedler>vmed)[0]
    v[K]=1
    return v, res

def isoperimetric(A, ground=None, residuals=None) :

  #from pyamg.graph import pseudo_peripheral_node
  #ground = pseudo_peripheral_node(A)[0]

  # select random ground 'dead' node
  if ground is None :
    seed()
    ground = randint(0,A.shape[0]-1)

  coarse = numpy.arange(0,A.shape[0])
  coarse = numpy.delete(coarse,ground,0)

  # remove ground node row and column
  L = A[coarse,:][:,coarse]
  r = numpy.ones((L.shape[0],))

  if residuals is None :
    res = []

  # construct preconditioner
  ml = smoothed_aggregation_solver(L,coarse_solver='pinv2')
  M = ml.aspreconditioner()

  # solve system using cg
  (x,flag) = cg(L,r,residuals=res,tol=1e-12,M=M)

  # use the median of solution, x, as the separator
  vmed = numpy.median(x)
  vmin = numpy.min(x)
  P1 = coarse[numpy.where(x<=vmed)[0]]
  P2 = coarse[numpy.where(x>vmed)[0]]

  weights = numpy.zeros((A.shape[0],))
  weights[P1] = x[numpy.where(x<=vmed)[0]]
  weights[P2] = x[numpy.where(x>vmed)[0]]
  weights[ground] = vmin-1

  P1 = numpy.append(P1,ground)

  return P1,P2,weights