Python gallery.load_example函数代码示例

    def setUp(self):
        self.cases = []

        # Poisson problems in 1D and 2D
        N = 20
        self.cases.append((poisson((2*N,), format='csr'), rand(2*N,)))     # 0
        self.cases.append((poisson((N, N), format='csr'), rand(N*N,)))  # 1
        # Boxed examples
        A = load_example('recirc_flow')['A'].tocsr()                      # 2
        self.cases.append((A, rand(A.shape[0],)))
        A = load_example('bar')['A'].tobsr(blocksize=(3, 3))               # 3
        self.cases.append((A, rand(A.shape[0],)))

    def test_nonhermitian(self):
        # problem data
        data = load_example('helmholtz_2D')
        A = data['A'].tocsr()
        B = data['B']
        numpy.random.seed(625)
        x0 = scipy.rand(A.shape[0]) + 1.0j * scipy.rand(A.shape[0])
        b = A * scipy.rand(A.shape[0]) + 1.0j * (A * scipy.rand(A.shape[0]))
        # solver parameters
        smooth = ('energy', {'krylov': 'gmres'})
        SA_build_args = {'max_coarse': 25, 'coarse_solver': 'pinv2',
                         'symmetry': 'symmetric'}
        SA_solve_args = {'cycle': 'V', 'maxiter': 20, 'tol': 1e-8}
        strength = [('evolution', {'k': 2, 'epsilon': 2.0})]
        smoother = ('gauss_seidel_nr', {'sweep': 'symmetric', 'iterations': 1})
        # Construct solver with nonsymmetric parameters
        sa = smoothed_aggregation_solver(A, B=B, smooth=smooth,
                                         strength=strength,
                                         presmoother=smoother,
                                         postsmoother=smoother,
                                         **SA_build_args)
        residuals = []
        # stand-alone solve
        x = sa.solve(b, x0=x0, residuals=residuals, **SA_solve_args)
        residuals = array(residuals)
        avg_convergence_ratio =\
            (residuals[-1] / residuals[0]) ** (1.0 / len(residuals))
        assert(avg_convergence_ratio < 0.85)
        # accelerated solve
        residuals = []
        x = sa.solve(b, x0=x0, residuals=residuals, accel='gmres',
                     **SA_solve_args)
        del x
        residuals = array(residuals)
        avg_convergence_ratio =\
            (residuals[-1] / residuals[0]) ** (1.0 / len(residuals))
        assert(avg_convergence_ratio < 0.6)

        # test that nonsymmetric parameters give the same result as symmetric
        # parameters for the complex-symmetric matrix A
        strength = 'symmetric'
        SA_build_args['symmetry'] = 'nonsymmetric'
        sa_nonsymm = smoothed_aggregation_solver(A, B=ones((A.shape[0], 1)),
                                                 smooth=smooth,
                                                 strength=strength,
                                                 presmoother=smoother,
                                                 postsmoother=smoother,
                                                 improve_candidates=None,
                                                 **SA_build_args)
        SA_build_args['symmetry'] = 'symmetric'
        sa_symm = smoothed_aggregation_solver(A, B=ones((A.shape[0], 1)),
                                              smooth=smooth,
                                              strength=strength,
                                              presmoother=smoother,
                                              postsmoother=smoother,
                                              improve_candidates=None,
                                              **SA_build_args)
        for (symm_lvl, nonsymm_lvl) in zip(sa_nonsymm.levels, sa_symm.levels):
            assert_array_almost_equal(symm_lvl.A.todense(),
                                      nonsymm_lvl.A.todense())

    def test_nonsymmetric(self):
        # problem data
        data = load_example('recirc_flow')
        A = data['A'].tocsr()
        B = data['B']
        numpy.random.seed(625)
        x0 = scipy.rand(A.shape[0])
        b = A * scipy.rand(A.shape[0])
        # solver parameters
        smooth = ('energy', {'krylov': 'gmres'})
        SA_build_args = {'max_coarse': 25, 'coarse_solver': 'pinv2',
                         'symmetry': 'nonsymmetric'}
        SA_solve_args = {'cycle': 'V', 'maxiter': 20, 'tol': 1e-8}
        strength = [('evolution', {'k': 2, 'epsilon': 8.0})]
        smoother = ('gauss_seidel_nr', {'sweep': 'symmetric', 'iterations': 1})
        improve_candidates = [('gauss_seidel_nr', {'sweep': 'symmetric',
                                                   'iterations': 4}), None]
        # Construct solver with nonsymmetric parameters
        sa = rootnode_solver(A, B=B, smooth=smooth,
                             improve_candidates=improve_candidates,
                             strength=strength,
                             presmoother=smoother,
                             postsmoother=smoother, **SA_build_args)
        residuals = []
        # stand-alone solve
        x = sa.solve(b, x0=x0, residuals=residuals, **SA_solve_args)
        residuals = array(residuals)
        avg_convergence_ratio =\
            (residuals[-1] / residuals[0]) ** (1.0 / len(residuals))
        # print "Test 1  %1.3e,  %1.3e" % (avg_convergence_ratio, 0.7)
        assert(avg_convergence_ratio < 0.7)
        # accelerated solve
        residuals = []
        x = sa.solve(b, x0=x0, residuals=residuals, accel='gmres',
                     **SA_solve_args)
        residuals = array(residuals)
        avg_convergence_ratio =\
            (residuals[-1] / residuals[0]) ** (1.0 / len(residuals))
        # print "Test 2  %1.3e,  %1.3e" % (avg_convergence_ratio, 0.45)
        assert(avg_convergence_ratio < 0.45)

        # test that nonsymmetric parameters give the same result as symmetric
        # parameters for Poisson problem
        A = poisson((15, 15), format='csr')
        strength = 'symmetric'
        SA_build_args['symmetry'] = 'nonsymmetric'
        sa_nonsymm = rootnode_solver(A, B=ones((A.shape[0], 1)), smooth=smooth,
                                     strength=strength,
                                     presmoother=smoother,
                                     postsmoother=smoother,
                                     improve_candidates=None, **SA_build_args)
        SA_build_args['symmetry'] = 'symmetric'
        sa_symm = rootnode_solver(A, B=ones((A.shape[0], 1)), smooth=smooth,
                                  strength=strength, presmoother=smoother,
                                  postsmoother=smoother,
                                  improve_candidates=None, **SA_build_args)
        for (symm_lvl, nonsymm_lvl) in zip(sa_nonsymm.levels, sa_symm.levels):
            assert_array_almost_equal(symm_lvl.A.todense(),
                                      nonsymm_lvl.A.todense())

    def test_nonhermitian(self):
        # problem data
        data = load_example("helmholtz_2D")
        A = data["A"].tocsr()
        B = data["B"]
        numpy.random.seed(625)
        x0 = scipy.rand(A.shape[0]) + 1.0j * scipy.rand(A.shape[0])
        b = A * scipy.rand(A.shape[0]) + 1.0j * (A * scipy.rand(A.shape[0]))
        # solver parameters
        smooth = ("energy", {"krylov": "gmres"})
        SA_build_args = {"max_coarse": 25, "coarse_solver": "pinv2", "symmetry": "symmetric"}
        SA_solve_args = {"cycle": "V", "maxiter": 20, "tol": 1e-8}
        strength = [("evolution", {"k": 2, "epsilon": 2.0})]
        smoother = ("gauss_seidel_nr", {"sweep": "symmetric", "iterations": 1})
        # Construct solver with nonsymmetric parameters
        sa = smoothed_aggregation_solver(
            A, B=B, smooth=smooth, strength=strength, presmoother=smoother, postsmoother=smoother, **SA_build_args
        )
        residuals = []
        # stand-alone solve
        x = sa.solve(b, x0=x0, residuals=residuals, **SA_solve_args)
        residuals = array(residuals)
        avg_convergence_ratio = (residuals[-1] / residuals[0]) ** (1.0 / len(residuals))
        assert avg_convergence_ratio < 0.85
        # accelerated solve
        residuals = []
        x = sa.solve(b, x0=x0, residuals=residuals, accel="gmres", **SA_solve_args)
        residuals = array(residuals)
        avg_convergence_ratio = (residuals[-1] / residuals[0]) ** (1.0 / len(residuals))
        assert avg_convergence_ratio < 0.6

        # test that nonsymmetric parameters give the same result as symmetric
        # parameters for the complex-symmetric matrix A
        strength = "symmetric"
        SA_build_args["symmetry"] = "nonsymmetric"
        sa_nonsymm = smoothed_aggregation_solver(
            A,
            B=ones((A.shape[0], 1)),
            smooth=smooth,
            strength=strength,
            presmoother=smoother,
            postsmoother=smoother,
            improve_candidates=None,
            **SA_build_args
        )
        SA_build_args["symmetry"] = "symmetric"
        sa_symm = smoothed_aggregation_solver(
            A,
            B=ones((A.shape[0], 1)),
            smooth=smooth,
            strength=strength,
            presmoother=smoother,
            postsmoother=smoother,
            improve_candidates=None,
            **SA_build_args
        )
        for (symm_lvl, nonsymm_lvl) in zip(sa_nonsymm.levels, sa_symm.levels):
            assert_array_almost_equal(symm_lvl.A.todense(), nonsymm_lvl.A.todense())

    def setUp(self):
        self.cases = []

        # Poisson problems in 1D and 2D
        for N in [2, 3, 5, 7, 10, 11, 19]:
            self.cases.append(poisson((N,), format='csr'))
        for N in [2, 3, 7, 9]:
            self.cases.append(poisson((N, N), format='csr'))

        for name in ['knot', 'airfoil', 'bar']:
            ex = load_example(name)
            self.cases.append(ex['A'].tocsr())

    def test_prefilter(self):
        """Check that using prefilter reduces NNZ in P"""
        np.random.seed(0)  # make tests repeatable
        cases = []

        # Simple, real-valued diffusion problems
        X = load_example('airfoil')
        A = X['A'].tocsr()
        B = X['B']

        cases.append((A, B,
                      ('energy', {'krylov': 'cg', 'degree': 2, 'maxiter': 3}),
                      {'theta': 0.05}))

        cases.append((A, B,
                      ('energy', {'krylov': 'gmres', 'degree': 2,
                                  'maxiter': 3}),
                      {'k': 3}))

        cases.append((A.tobsr(blocksize=(2, 2)),
                      np.hstack((B, np.random.rand(B.shape[0], 1))),
                      ('energy', {'krylov': 'cg', 'degree': 2,
                                  'maxiter': 3}),
                      {'theta': 0.1}))

        # Simple, imaginary-valued problems
        iA = 1.0j * A
        iB = 1.0 + np.random.rand(iA.shape[0], 2)\
                 + 1.0j * (1.0 + np.random.rand(iA.shape[0], 2))

        cases.append((iA, B,
                      ('energy', {'krylov': 'cg', 'degree': 2, 'maxiter': 3}),
                      {'theta': 0.05}))

        cases.append((iA, iB,
                      ('energy', {'krylov': 'gmres', 'degree': 2,
                                  'maxiter': 3}),
                      {'k': 3}))

        cases.append((A.tobsr(blocksize=(2, 2)),
                      np.hstack((B, np.random.rand(B.shape[0], 1))),
                      ('energy', {'krylov': 'cg', 'degree': 2, 'maxiter': 3}),
                      {'theta': 0.1}))

        for A, B, smooth, prefilter in cases:
            ml_nofilter = rootnode_solver(A, B=B, max_coarse=1, max_levels=2,
                                          smooth=smooth, keep=True)
            smooth[1]['prefilter'] = prefilter
            ml_filter = rootnode_solver(A, B=B, max_coarse=1, max_levels=2,
                                        smooth=smooth, keep=True)
            assert_equal(ml_nofilter.levels[0].P.nnz >
                         ml_filter.levels[0].P.nnz, True)

    def setUp(self):
        self.cases = []

        # random matrices
        np.random.seed(0)
        for N in [2, 3, 5]:
            self.cases.append(csr_matrix(np.random.rand(N, N)))

        # Poisson problems in 1D and 2D
        for N in [2, 3, 5, 7, 10, 11, 19]:
            self.cases.append(poisson((N,), format='csr'))
        for N in [2, 3, 5, 7, 8]:
            self.cases.append(poisson((N, N), format='csr'))

        for name in ['knot', 'airfoil', 'bar']:
            ex = load_example(name)
            self.cases.append(ex['A'].tocsr())

    def test_distance_strength_of_connection(self):
        data = load_example('airfoil')
        cases = []
        cases.append((data['A'].tocsr(), data['vertices']))

        for (A, V) in cases:
            for theta in [1.5, 2.0, 2.5]:
                result = distance_soc(A, V, theta=theta)
                expected = reference_distance_soc(A, V, theta=theta)
                assert_equal(result.nnz,       expected.nnz)
                assert_array_almost_equal(result.todense(), expected.todense())

        for (A, V) in cases:
            for theta in [0.5, 1.0, 1.5]:
                result = distance_soc(A, V, theta=theta, relative_drop=False)
                expected = reference_distance_soc(A, V, theta=theta,
                                                  relative_drop=False)
                assert_equal(result.nnz,       expected.nnz)
                assert_array_almost_equal(result.todense(), expected.todense())

    def setUp(self):
        self.cases = []

        # random matrices
        numpy.random.seed(0)
        for N in [2,3,5]:
            self.cases.append( csr_matrix(rand(N,N)) + csr_matrix(1.0j*rand(N,N)))

        # Poisson problems in 1D and 2D
        for N in [2,3,5,7,10,11,19]:
            A = poisson( (N,), format='csr'); A.data = A.data + 1.0j*A.data;
            self.cases.append(A)
        for N in [2,3,7,9]:
            A = poisson( (N,N), format='csr'); A.data = A.data + 1.0j*rand(A.data.shape[0],);
            self.cases.append(A)

        for name in ['knot','airfoil','bar']:
            ex = load_example(name)
            A = ex['A'].tocsr(); A.data = A.data + 0.5j*rand(A.data.shape[0],);
            self.cases.append(A)

    def setUp(self):
        self.cases = []

        #

        # Random matrices, cases 0-2
        np.random.seed(0)
        for N in [2, 3, 5]:
            self.cases.append(csr_matrix(np.random.rand(N, N)))

        # Poisson problems in 1D, cases 3-9
        for N in [2, 3, 5, 7, 10, 11, 19]:
            self.cases.append(poisson((N,), format="csr"))

        # Poisson problems in 2D, cases 10-15
        for N in [2, 3, 5, 7, 10, 11]:
            self.cases.append(poisson((N, N), format="csr"))

        for name in ["knot", "airfoil", "bar"]:
            ex = load_example(name)
            self.cases.append(ex["A"].tocsr())

    def setUp(self):
        cases = []
        seed(0)

        for i in range(5):
            A = rand(8, 8) > 0.5
            cases.append(canonical_graph(A + A.T).astype(float))

        cases.append(zeros((1, 1)))
        cases.append(zeros((2, 2)))
        cases.append(zeros((8, 8)))
        cases.append(ones((2, 2)) - eye(2))
        cases.append(poisson((5,)))
        cases.append(poisson((5, 5)))
        cases.append(poisson((11, 11)))
        cases.append(poisson((5, 5, 5)))
        for name in ['airfoil', 'bar', 'knot']:
            cases.append(load_example(name)['A'])

        cases = [canonical_graph(G) for G in cases]

        self.cases = cases

    def test_distance(self):
        data = load_example('airfoil')
        cases = []
        cases.append((data['A'].tocsr(), data['vertices']))
       
        for (A, V) in cases:
            dim = V.shape[1]
            for theta in [1.5, 2.0, 2.5]:
                cost = [0]
                lower_bound = 3*dim + float(A.shape[0]) / A.nnz
                upper_bound = 3*dim + 3
                distance_strength_of_connection(A, V, theta=theta, relative_drop=True, cost=cost)
                assert(cost[0] >= lower_bound)
                assert(cost[0] <= upper_bound)

        for (A, V) in cases:
            for theta in [0.5, 1.0, 1.5]:
                cost = [0]
                lower_bound = 3*dim + float(A.shape[0]) / A.nnz
                upper_bound = 3*dim + 3
                distance_strength_of_connection(A, V, theta=theta, relative_drop=False, cost=cost)
                assert(cost[0] >= lower_bound)
                assert(cost[0] <= upper_bound)

    def test_nonsymmetric(self):
        # problem data
        data = load_example("recirc_flow")
        A = data["A"].tocsr()
        B = data["B"]
        numpy.random.seed(625)
        x0 = scipy.rand(A.shape[0])
        b = A * scipy.rand(A.shape[0])
        # solver parameters
        smooth = ("energy", {"krylov": "gmres"})
        SA_build_args = {"max_coarse": 25, "coarse_solver": "pinv2", "symmetry": "nonsymmetric"}
        SA_solve_args = {"cycle": "V", "maxiter": 20, "tol": 1e-8}
        strength = [("evolution", {"k": 2, "epsilon": 8.0})]
        smoother = ("gauss_seidel_nr", {"sweep": "symmetric", "iterations": 1})
        improve_candidates = [("gauss_seidel_nr", {"sweep": "symmetric", "iterations": 4}), None]
        # Construct solver with nonsymmetric parameters
        sa = rootnode_solver(
            A,
            B=B,
            smooth=smooth,
            improve_candidates=improve_candidates,
            strength=strength,
            presmoother=smoother,
            postsmoother=smoother,
            **SA_build_args
        )
        residuals = []
        # stand-alone solve
        x = sa.solve(b, x0=x0, residuals=residuals, **SA_solve_args)
        residuals = array(residuals)
        avg_convergence_ratio = (residuals[-1] / residuals[0]) ** (1.0 / len(residuals))
        # print "Test 1  %1.3e,  %1.3e" % (avg_convergence_ratio, 0.7)
        assert avg_convergence_ratio < 0.7
        # accelerated solve
        residuals = []
        x = sa.solve(b, x0=x0, residuals=residuals, accel="gmres", **SA_solve_args)
        del x
        residuals = array(residuals)
        avg_convergence_ratio = (residuals[-1] / residuals[0]) ** (1.0 / len(residuals))
        # print "Test 2  %1.3e,  %1.3e" % (avg_convergence_ratio, 0.45)
        assert avg_convergence_ratio < 0.45

        # test that nonsymmetric parameters give the same result as symmetric
        # parameters for Poisson problem
        A = poisson((15, 15), format="csr")
        strength = "symmetric"
        SA_build_args["symmetry"] = "nonsymmetric"
        sa_nonsymm = rootnode_solver(
            A,
            B=ones((A.shape[0], 1)),
            smooth=smooth,
            strength=strength,
            presmoother=smoother,
            postsmoother=smoother,
            improve_candidates=None,
            **SA_build_args
        )
        SA_build_args["symmetry"] = "symmetric"
        sa_symm = rootnode_solver(
            A,
            B=ones((A.shape[0], 1)),
            smooth=smooth,
            strength=strength,
            presmoother=smoother,
            postsmoother=smoother,
            improve_candidates=None,
            **SA_build_args
        )
        for (symm_lvl, nonsymm_lvl) in zip(sa_nonsymm.levels, sa_symm.levels):
            assert_array_almost_equal(symm_lvl.A.todense(), nonsymm_lvl.A.todense())

    def test_range(self):
        """Check that P*R=B"""
        np.random.seed(0)  # make tests repeatable

        cases = []

        # Simple, real-valued diffusion problems
        X = load_example("airfoil")
        A = X["A"].tocsr()
        B = X["B"]
        cases.append((A, B, ("jacobi", {"filter": True, "weighting": "local"})))
        cases.append((A, B, ("jacobi", {"filter": True, "weighting": "block"})))

        cases.append((A, B, ("energy", {"maxiter": 3})))
        cases.append((A, B, ("energy", {"krylov": "cgnr"})))
        cases.append((A, B, ("energy", {"krylov": "gmres", "degree": 2})))

        A = poisson((10, 10), format="csr")
        B = np.ones((A.shape[0], 1))
        cases.append((A, B, ("jacobi", {"filter": True, "weighting": "diagonal"})))
        cases.append((A, B, ("jacobi", {"filter": True, "weighting": "local"})))

        cases.append((A, B, "energy"))
        cases.append((A, B, ("energy", {"degree": 2})))
        cases.append((A, B, ("energy", {"krylov": "cgnr", "degree": 2})))
        cases.append((A, B, ("energy", {"krylov": "gmres"})))

        # Simple, imaginary-valued problems
        iA = 1.0j * A
        iB = 1.0 + np.random.rand(iA.shape[0], 2) + 1.0j * (1.0 + np.random.rand(iA.shape[0], 2))

        cases.append((iA, B, ("jacobi", {"filter": True, "weighting": "diagonal"})))
        cases.append((iA, B, ("jacobi", {"filter": True, "weighting": "block"})))
        cases.append((iA, iB, ("jacobi", {"filter": True, "weighting": "local"})))
        cases.append((iA, iB, ("jacobi", {"filter": True, "weighting": "block"})))

        cases.append((iA.tobsr(blocksize=(5, 5)), B, ("jacobi", {"filter": True, "weighting": "block"})))
        cases.append((iA.tobsr(blocksize=(5, 5)), iB, ("jacobi", {"filter": True, "weighting": "block"})))

        cases.append((iA, B, ("energy", {"krylov": "cgnr", "degree": 2})))
        cases.append((iA, iB, ("energy", {"krylov": "cgnr"})))
        cases.append(
            (
                iA.tobsr(blocksize=(5, 5)),
                B,
                ("energy", {"krylov": "cgnr", "degree": 2, "maxiter": 3, "postfilter": {"theta": 0.05}}),
            )
        )
        cases.append(
            (
                iA.tobsr(blocksize=(5, 5)),
                B,
                ("energy", {"krylov": "cgnr", "degree": 2, "maxiter": 3, "prefilter": {"theta": 0.05}}),
            )
        )
        cases.append((iA.tobsr(blocksize=(5, 5)), B, ("energy", {"krylov": "cgnr", "degree": 2, "maxiter": 3})))
        cases.append((iA.tobsr(blocksize=(5, 5)), iB, ("energy", {"krylov": "cgnr"})))

        cases.append((iA, B, ("energy", {"krylov": "gmres"})))
        cases.append((iA, iB, ("energy", {"krylov": "gmres", "degree": 2})))
        cases.append((iA.tobsr(blocksize=(5, 5)), B, ("energy", {"krylov": "gmres", "degree": 2, "maxiter": 3})))
        cases.append((iA.tobsr(blocksize=(5, 5)), iB, ("energy", {"krylov": "gmres"})))

        # Simple, imaginary-valued problems
        iA = A + 1.0j * scipy.sparse.eye(A.shape[0], A.shape[1])

        cases.append((iA, B, ("jacobi", {"filter": True, "weighting": "local"})))
        cases.append((iA, B, ("jacobi", {"filter": True, "weighting": "block"})))
        cases.append((iA, iB, ("jacobi", {"filter": True, "weighting": "diagonal"})))
        cases.append((iA, iB, ("jacobi", {"filter": True, "weighting": "block"})))
        cases.append((iA.tobsr(blocksize=(4, 4)), iB, ("jacobi", {"filter": True, "weighting": "block"})))

        cases.append((iA, B, ("energy", {"krylov": "cgnr"})))
        cases.append((iA.tobsr(blocksize=(4, 4)), iB, ("energy", {"krylov": "cgnr"})))

        cases.append((iA, B, ("energy", {"krylov": "gmres"})))
        cases.append((iA.tobsr(blocksize=(4, 4)), iB, ("energy", {"krylov": "gmres", "degree": 2, "maxiter": 3})))

        cases.append(
            (
                iA.tobsr(blocksize=(4, 4)),
                iB,
                ("energy", {"krylov": "gmres", "degree": 2, "maxiter": 3, "postfilter": {"theta": 0.05}}),
            )
        )
        cases.append(
            (
                iA.tobsr(blocksize=(4, 4)),
                iB,
                ("energy", {"krylov": "gmres", "degree": 2, "maxiter": 3, "prefilter": {"theta": 0.05}}),
            )
        )

        A = gauge_laplacian(10, spacing=1.0, beta=0.21)
        B = np.ones((A.shape[0], 1))
        cases.append((A, iB, ("jacobi", {"filter": True, "weighting": "diagonal"})))
        cases.append((A, iB, ("jacobi", {"filter": True, "weighting": "local"})))

        cases.append((A, B, ("energy", {"krylov": "cg"})))
        cases.append((A, iB, ("energy", {"krylov": "cgnr"})))
#.........这里部分代码省略.........

    def test_evolution_strength_of_connection(self):
        # Params:  A, B, epsilon=4.0, k=2, proj_type="l2"
        cases = []

        # Ensure that isotropic diffusion results in isotropic strength stencil
        for N in [3, 5, 7, 10]:
            A = poisson((N,), format='csr')
            B = np.ones((A.shape[0], 1))
            cases.append({'A': A.copy(), 'B': B.copy(), 'epsilon': 4.0,
                          'k': 2, 'proj': 'l2'})

        # Ensure that anisotropic diffusion results in an anisotropic
        # strength stencil
        for N in [3, 6, 7]:
            u = np.ones(N*N)
            A = spdiags([-u, -0.001*u, 2.002*u, -0.001*u, -u],
                        [-N, -1, 0, 1, N], N*N, N*N, format='csr')
            B = np.ones((A.shape[0], 1))
            cases.append({'A': A.copy(), 'B': B.copy(), 'epsilon': 4.0,
                          'k': 2, 'proj': 'l2'})

        # Ensure that isotropic elasticity results in an isotropic stencil
        for N in [3, 6, 7]:
            (A, B) = linear_elasticity((N, N), format='bsr')
            cases.append({'A': A.copy(), 'B': B.copy(), 'epsilon': 32.0,
                          'k': 8, 'proj': 'D_A'})

        # Run an example with a non-uniform stencil
        ex = load_example('airfoil')
        A = ex['A'].tocsr()
        B = np.ones((A.shape[0], 1))
        cases.append({'A': A.copy(), 'B': B.copy(), 'epsilon': 8.0, 'k': 4,
                      'proj': 'D_A'})
        Absr = A.tobsr(blocksize=(5, 5))
        cases.append({'A': Absr.copy(), 'B': B.copy(), 'epsilon': 8.0, 'k': 4,
                      'proj': 'D_A'})
        # Different B
        B = arange(1, 2*A.shape[0]+1, dtype=float).reshape(-1, 2)
        cases.append({'A': A.copy(), 'B': B.copy(), 'epsilon': 4.0, 'k': 2,
                      'proj': 'l2'})
        cases.append({'A': Absr.copy(), 'B': B.copy(), 'epsilon': 4.0, 'k': 2,
                      'proj': 'l2'})

        # Zero row and column
        A.data[A.indptr[4]:A.indptr[5]] = 0.0
        A = A.tocsc()
        A.data[A.indptr[4]:A.indptr[5]] = 0.0
        A.eliminate_zeros()
        A = A.tocsr()
        A.sort_indices()
        cases.append({'A': A.copy(), 'B': B.copy(), 'epsilon': 4.0, 'k': 2,
                      'proj': 'l2'})
        Absr = A.tobsr(blocksize=(5, 5))
        cases.append({'A': Absr.copy(), 'B': B.copy(), 'epsilon': 4.0, 'k': 2,
                      'proj': 'l2'})

        for ca in cases:
            scipy.random.seed(0)  # make results deterministic
            result = evolution_soc(ca['A'], ca['B'], epsilon=ca['epsilon'],
                                   k=ca['k'], proj_type=ca['proj'],
                                   symmetrize_measure=False)
            scipy.random.seed(0)  # make results deterministic
            expected = reference_evolution_soc(ca['A'], ca['B'],
                                               epsilon=ca['epsilon'],
                                               k=ca['k'], proj_type=ca['proj'])
            assert_array_almost_equal(result.todense(), expected.todense(),
                                      decimal=4)

        # Test Scale Invariance for multiple near nullspace candidates
        (A, B) = linear_elasticity((5, 5), format='bsr')
        scipy.random.seed(0)  # make results deterministic
        result_unscaled = evolution_soc(A, B, epsilon=4.0,
                                        k=2, proj_type="D_A",
                                        symmetrize_measure=False)
        # create scaled A
        D = spdiags([arange(A.shape[0], 2*A.shape[0], dtype=float)],
                    [0], A.shape[0], A.shape[0], format='csr')
        Dinv = spdiags([1.0/arange(A.shape[0], 2*A.shape[0], dtype=float)],
                       [0], A.shape[0], A.shape[0], format='csr')
        scipy.random.seed(0)  # make results deterministic
        result_scaled = evolution_soc((D*A*D).tobsr(blocksize=(2, 2)),
                                      Dinv*B, epsilon=4.0, k=2,
                                      proj_type="D_A",
                                      symmetrize_measure=False)
        assert_array_almost_equal(result_scaled.todense(),
                                  result_unscaled.todense(), decimal=2)