"""Test Krylov methods."""
import numpy as np
from numpy.testing import (TestCase, assert_array_almost_equal,
                           assert_equal, assert_almost_equal)
from scipy.linalg import solve
from scipy import sparse

import pyamg
from pyamg.util.linalg import norm
from pyamg.krylov import bicgstab, cg, cgne, cgnr, cr, fgmres, gmres, steepest_descent
from pyamg.krylov._gmres_householder import gmres_householder
from pyamg.krylov._gmres_mgs import gmres_mgs

import pytest


class TestStoppingCriteria(TestCase):
    def setUp(self):
        self.cases = []

        np.random.seed(9062883)
        n = 10
        A = np.random.rand(n, n)
        b = np.random.rand(n)
        x0 = np.random.rand(n)
        A = 0.5 * (A + A.T) + n*np.eye(n, n)
        self.cases.append({'A': A, 'b': b, 'x0': x0, 'tol': 1e-8})
        self.cases.append({'A': sparse.csr_array(A), 'b': b, 'x0': x0, 'tol': 1e-8})

    def test_stoppingcriteria(self):
        for method, crits in [
                (cg,       ('rr', 'rr+', 'MrMr', 'rMr')),
                (bicgstab, ('rr', 'rr+')),
                (cgne,     ('rr', 'rr+', 'MrMr', 'rMr')),
                (cgnr,     ('rr', 'rr+', 'MrMr', 'rMr')),
                (cr,     ('rr', 'rr+', 'MrMr')),
                (steepest_descent,     ('rr', 'rr+', 'MrMr', 'rMr'))]:
            for criteria in crits:
                for case in self.cases:
                    A = case['A']
                    b = case['b']
                    maxiter = None
                    if method.__name__ == 'steepest_descent':
                        maxiter = 100
                    x1, info = method(A, b, x0=case['x0'], tol=case['tol'],
                                      criteria=criteria, maxiter=maxiter)
                    assert_equal(info, 0, err_msg=f'Problem in {method.__name__}.')
                    assert_almost_equal(np.linalg.norm(b - A @ x1), 0.0,
                                        decimal=6,
                                        err_msg=f'Problem in {method.__name__}.')


class TestKrylov(TestCase):
    def setUp(self):
        self.cases = []
        self.spd_cases = []
        self.symm_cases = []

        # self.oblique = [gmres, fgmres, cgnr,
        #                 krylov._gmres_householder.gmres_householder,
        #                 krylov._gmres_mgs.gmres_mgs]
        self.oblique = [gmres_householder, gmres_mgs, gmres, fgmres, cgnr]
        self.symm_oblique = [cr]
        self.orth = [cgne]
        self.inexact = [bicgstab]
        self.spd_orth = [cg]

        # 1x1
        A = np.array([[1.2]])
        b = np.array([3.9]).reshape(-1, 1)
        x0 = np.zeros((1, 1))
        self.cases.append({'A': A, 'b': b, 'x0': x0, 'tol': 1e-16,
                           'maxiter': 1, 'reduction_factor': 1e-10})
        self.spd_cases.append({'A': A, 'b': b, 'x0': x0, 'tol': 1e-16,
                               'maxiter': 1, 'reduction_factor': 1e-10})
        self.symm_cases.append({'A': A, 'b': b, 'x0': x0, 'tol': 1e-16,
                                'maxiter': 1, 'reduction_factor': 1e-10})

        # 4x4
        A = np.array([[1.2, 0., 0., 0.],
                      [0., 4., 2., 6.],
                      [0., 0., 9.3, -2.31],
                      [-4., 0., 0., -11.]])
        b = np.array([1., 3.9, 0., -1.23]).reshape(-1, 1)
        x0 = np.zeros((4, 1))
        self.cases.append({'A': A, 'b': b, 'x0': x0, 'tol': 1e-16,
                           'maxiter': 4, 'reduction_factor': 1e-10})
        self.spd_cases.append({'A': A.T.dot(A), 'b': b, 'x0': x0, 'tol': 1e-16,
                               'maxiter': 4, 'reduction_factor': 1e-10})
        self.symm_cases.append({'A': A.T + A, 'b': b, 'x0': x0, 'tol': 1e-16,
                                'maxiter': 4, 'reduction_factor': 1e-10})

        # 4x4 Imaginary
        A = np.array(A, dtype=complex)
        A[0, 0] += 3.1j
        A[3, 3] -= 1.34j
        A[1, 3] *= 1.0j
        A[1, 2] += 1.0j
        b = np.array([1. - 1.0j, 2.0 - 3.9j, 0., -1.23]).reshape(-1, 1)
        x0 = np.ones((4, 1))
        self.cases.append({'A': A, 'b': b, 'x0': x0, 'tol': 1e-16,
                           'maxiter': 4, 'reduction_factor': 1e-10})
        self.spd_cases.append({'A': A.conj().T.dot(A), 'b': b, 'x0': x0, 'tol': 1e-16,
                               'maxiter': 4, 'reduction_factor': 1e-10})
        self.symm_cases.append({'A': A.conj().T + A, 'b': b, 'x0': x0, 'tol': 1e-16,
                                'maxiter': 4, 'reduction_factor': 1e-10})

        # 10x10
        A = np.array([[-1.1, 0., 0., 0., 3.9, 0., 0., 11., -1., 0.],
                      [0., 4., 2.9, 0., 0., 6.8, 0., 0., 0., 0.],
                      [0., 0., 9.0, 0., 0., 0.8, 1., -2.2, 0., 9.],
                      [-4., 0., 0.0, 0., 0., 0.0, 2., 2.2, 0., 0.],
                      [0., 0., 0.0, 21., 0., 0.1, 0., 0., 0., 0.],
                      [0., 0., 0.0, 0., -4.7, 0.0, 0., 0., 0., 0.],
                      [2.1, 7., 22.0, 0., 0., 0.0, 0., 0., 0., 0.],
                      [0., 0., 0.0, 34., 0., 0.0, 0., 0., -12.3, 0.],
                      [0., 3.4, 0.0, 0., 0., -0.3, 0., 0., 0., 0.],
                      [9., 0., 0.0, 0., 87., 0.0, 0., 0., 0., -11.2]])
        b = np.array([1., 0., 0.2, 8., 0., -1.9,
                     11.3, 0.0, 0.1, 0.0]).reshape(-1, 1)
        x0 = np.zeros((10, 1))
        x0[4] = 11.1
        x0[7] = -2.
        self.cases.append({'A': A, 'b': b, 'x0': x0, 'tol': 1e-16,
                           'maxiter': 2, 'reduction_factor': 0.98})
        self.symm_cases.append({'A': A + A.T, 'b': b, 'x0': x0, 'tol': 1e-16,
                                'maxiter': 2, 'reduction_factor': 0.98})
        self.spd_cases.append({'A':
                               pyamg.gallery.poisson((10,)).toarray(),
                               'b': b, 'x0': x0, 'tol': 1e-16, 'maxiter': 2,
                               'reduction_factor': 0.98})

    def test_gmres(self):
        # Ensure repeatability
        np.random.seed(0)

        #  For these small matrices, Householder and MGS GMRES should give the
        #  same result, and for symmetric (but possibly indefinite) matrices CR
        #  and GMRES should give same result
        for maxiter in [1, 2, 3]:
            for case, symm_case in zip(self.cases, self.symm_cases):
                A = case['A']
                b = case['b']
                x0 = case['x0']
                A_symm = symm_case['A']
                b_symm = symm_case['b']
                x0_symm = symm_case['x0']

                # Test agreement between Householder and GMRES
                (x, flag) = gmres_householder(A, b, x0=x0,
                                              maxiter=min(A.shape[0], maxiter))
                (x2, flag2) = gmres_mgs(A, b, x0=x0, maxiter=min(A.shape[0],
                                                                 maxiter))
                err_msg = ('Householder GMRES and MGS GMRES gave '
                           'different results for small matrix')
                assert_array_almost_equal(x/norm(x), x2/norm(x2),
                                          err_msg=err_msg)

                err_msg = ('Householder GMRES and MGS GMRES returned '
                           'different convergence flags for small matrix')
                assert_equal(flag, flag2, err_msg=err_msg)

                # Test agreement between GMRES and CR
                if A_symm.shape[0] > 1:
                    residuals2 = []
                    (x2, flag2) = gmres_mgs(A_symm, b_symm, x0=x0_symm,
                                            maxiter=maxiter,
                                            restart=None,
                                            residuals=residuals2)
                    residuals3 = []
                    (x3, flag2) = cr(A_symm, b_symm, x0=x0_symm,
                                     maxiter=min(A.shape[0], maxiter),
                                     residuals=residuals3)
                    residuals2 = np.array(residuals2)
                    residuals3 = np.array(residuals3)

                    err_msg = 'CR and GMRES yield different residual vectors'
                    assert_array_almost_equal(residuals3/norm(residuals3),
                                              residuals2/norm(residuals2),
                                              err_msg=err_msg)

                    err_msg = 'CR and GMRES yield different answers'
                    assert_array_almost_equal(x2/norm(x2), x3/norm(x3),
                                              err_msg=err_msg)

    def test_krylov(self):
        # Oblique projectors reduce the residual
        for method in self.oblique:
            for case in self.cases:
                A = case['A']
                b = case['b']
                x0 = case['x0']
                factor = case['reduction_factor']
                xNew, _ = method(A, b, x0=x0, tol=case['tol'],
                                 maxiter=case['maxiter'])
                xNew = xNew.reshape(-1, 1)
                assert_equal((norm(b - A.dot(xNew)) / norm(b - A.dot(x0))) < factor,
                             True, err_msg='Oblique Krylov Method Failed Test')

        # Oblique projectors reduce the residual, here we consider oblique
        # projectors for symmetric matrices
        for method in self.symm_oblique:
            for case in self.symm_cases:
                A = case['A']
                b = case['b']
                x0 = case['x0']
                factor = case['reduction_factor']
                xNew, _ = method(A, b, x0=x0, tol=case['tol'],
                                 maxiter=case['maxiter'])
                xNew = xNew.reshape(-1, 1)
                assert_equal((norm(b - A.dot(xNew)) / norm(b - A.dot(x0))) < factor,
                             True, err_msg='Symmetric oblique Krylov Method Failed')

        # Orthogonal projectors reduce the error
        for method in self.orth:
            for case in self.cases:
                A = case['A']
                b = case['b']
                x0 = case['x0']
                factor = case['reduction_factor']
                xNew, _ = method(A, b, x0=x0, tol=case['tol'],
                                 maxiter=case['maxiter'])
                xNew = xNew.reshape(-1, 1)
                soln = solve(A, b)
                assert_equal((norm(soln - xNew) / norm(soln - x0)) < factor,
                             True, err_msg='Orthogonal Krylov Method Failed Test')

        # SPD Orthogonal projectors reduce the error
        for method in self.spd_orth:
            for case in self.spd_cases:
                A = case['A']
                b = case['b']
                x0 = case['x0']
                factor = case['reduction_factor']
                xNew, _ = method(A, b, x0=x0, tol=case['tol'],
                                 maxiter=case['maxiter'])
                xNew = xNew.reshape(-1, 1)
                soln = solve(A, b)
                assert_equal((norm(soln - xNew) / norm(soln - x0)) < factor,
                             True, err_msg='Orthogonal Krylov Method Failed Test')

        # Assume that Inexact Methods reduce the residual for these examples
        for method in self.inexact:
            for case in self.cases:
                A = case['A']
                b = case['b']
                x0 = case['x0']
                xNew, _ = method(A, b, x0=x0, tol=case['tol'],
                                 maxiter=A.shape[0])
                xNew = xNew.reshape(-1, 1)
                assert_equal((norm(b - A.dot(xNew)) / norm(b - A.dot(x0))) < 0.35,
                             True, err_msg='Inexact Krylov Method Failed Test')


np.random.seed(751537155)
A = pyamg.gallery.poisson((10,), format='csr')
b = np.random.rand(A.shape[0])


@pytest.mark.parametrize('method', [fgmres, gmres_mgs, gmres_householder, gmres,
                                    bicgstab, cg, cgne, cgnr, cr])
def test_defaults(method):
    x, info = method(A, b)
    assert info == 0
    assert np.linalg.norm(b - A @ x) < np.linalg.norm(b)
