using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using Azimuth.SparseLinearAlgebra;
using Azimuth.DenseLinearAlgebra;

namespace Azimuth.UnitTests.SparseLinearAlgebra
{
    public class SymmetricTridiagonalTests
    {
        public class Ctor : UnitTestSharp.TestFixture
        {
            public void Basic()
            {
                DenseVector diagonal = new Scalar[] { 1, 2, 3, 4, };
                DenseVector subDiagonal = new Scalar[] { 1, 1, 3, };

                var tridiagonal = new SymmetricTridiagonal(ref diagonal, ref subDiagonal);

                DenseVector expectedDiagonal = new Scalar[] { 1, 2, 3, 4, };
                DenseVector expectedSubDiagonal = new Scalar[] { 1, 1, 3, };

                CheckEqual(expectedDiagonal, tridiagonal.Diagonal);
                CheckEqual(expectedSubDiagonal, tridiagonal.SubDiagonal);
            }

            public void LengthMismatch()
            {
                DenseVector diagonal = new Scalar[] { 1, 2, 3, 4, };
                DenseVector subDiagonal = new Scalar[] { 1, 1, 3, 3, };

                CheckThrow(typeof(Exception));
                var tridiagonal = new SymmetricTridiagonal(ref diagonal, ref subDiagonal);
            }
        }

        public class EigenvaluesAndVectors : UnitTestSharp.TestFixture
        {
            public void Identity2()
            {
                DenseVector diagonal = new Scalar[] { 1, 1, };
                DenseVector subDiagonal = new Scalar[] { 0, };

                var tridiagonal = new SymmetricTridiagonal(ref diagonal, ref subDiagonal);

                var eigenvectors = new SquareMatrix(2);
                var eigenvalues = new DenseVector(2);

                eigenvectors.SetIdentity();
                tridiagonal.ComputeEigenValuesAndVectors(ref eigenvalues, ref eigenvectors);

                DenseVector expectedEigenvalues = new Scalar[] { 1, 1, };
                SquareMatrix expectedEigenvectors = new Scalar[,] 
                {
                    { 1, 0, },
                    { 0, 1, },
                };

                CheckEqual(expectedEigenvalues, eigenvalues);
                CheckEqual(expectedEigenvectors, eigenvectors);
            }
            
            public void Identity4()
            {
                DenseVector diagonal = new Scalar[] { 1, 1, 1, 1, };
                DenseVector subDiagonal = new Scalar[] { 0, 0, 0, };

                var tridiagonal = new SymmetricTridiagonal(ref diagonal, ref subDiagonal);

                var eigenvectors = new SquareMatrix(4);
                var eigenvalues = new DenseVector(4);

                eigenvectors.SetIdentity();
                tridiagonal.ComputeEigenValuesAndVectors(ref eigenvalues, ref eigenvectors);

                DenseVector expectedEigenvalues = new Scalar[] { 1, 1, 1, 1, };
                SquareMatrix expectedEigenvectors = new Scalar[,] 
                {
                    { 1, 0, 0, 0, },
                    { 0, 1, 0, 0, },
                    { 0, 0, 1, 0, },
                    { 0, 0, 0, 1, },
                };

                CheckEqual(expectedEigenvalues, eigenvalues);
                CheckEqual(expectedEigenvectors, eigenvectors);
            }

            public void Diagonal_RandomOrder()
            {
                DenseVector diagonal = new Scalar[] { 4, 2, 1, 3, };
                DenseVector subDiagonal = new Scalar[] { 0, 0, 0, };

                var tridiagonal = new SymmetricTridiagonal(ref diagonal, ref subDiagonal);

                var eigenvectors = new SquareMatrix(4);
                var eigenvalues = new DenseVector(4);

                eigenvectors.SetIdentity();
                tridiagonal.ComputeEigenValuesAndVectors(ref eigenvalues, ref eigenvectors);

                DenseVector expectedEigenvalues = new Scalar[] { 4, 2, 1, 3, };
                SquareMatrix expectedEigenvectors = new Scalar[,] 
                {
                    { 1, 0, 0, 0, },
                    { 0, 1, 0, 0, },
                    { 0, 0, 1, 0, },
                    { 0, 0, 0, 1, },
                };

                CheckEqual(expectedEigenvalues, eigenvalues);
                CheckEqual(expectedEigenvectors, eigenvectors);
            }

            public void DiagonalDegenerate()
            {
                DenseVector diagonal = new Scalar[] { 0, 0, 0, 0, };
                DenseVector subDiagonal = new Scalar[] { 0, 0, 0, };

                var tridiagonal = new SymmetricTridiagonal(ref diagonal, ref subDiagonal);

                var eigenvectors = new SquareMatrix(4);
                var eigenvalues = new DenseVector(4);

                eigenvectors.SetIdentity();
                tridiagonal.ComputeEigenValuesAndVectors(ref eigenvalues, ref eigenvectors);

                DenseVector expectedEigenvalues = new Scalar[] { 0, 0, 0, 0, };
                SquareMatrix expectedEigenvectors = new Scalar[,] 
                {
                    { 1, 0, 0, 0, },
                    { 0, 1, 0, 0, },
                    { 0, 0, 1, 0, },
                    { 0, 0, 0, 1, },
                };

                CheckEqual(expectedEigenvalues, eigenvalues);
                CheckEqual(expectedEigenvectors, eigenvectors);
            }

            public void DiagonalMultiplicities()
            {
                DenseVector diagonal = new Scalar[] { 1, 2, 3, 2, };
                DenseVector subDiagonal = new Scalar[] { 0, 0, 0, };

                var tridiagonal = new SymmetricTridiagonal(ref diagonal, ref subDiagonal);

                var eigenvectors = new SquareMatrix(4);
                var eigenvalues = new DenseVector(4);

                eigenvectors.SetIdentity();
                tridiagonal.ComputeEigenValuesAndVectors(ref eigenvalues, ref eigenvectors);

                DenseVector expectedEigenvalues = new Scalar[] { 1, 2, 3, 2, };
                SquareMatrix expectedEigenvectors = new Scalar[,] 
                {
                    { 1, 0, 0, 0, },
                    { 0, 1, 0, 0, },
                    { 0, 0, 1, 0, },
                    { 0, 0, 0, 1, },
                };

                CheckEqual(expectedEigenvalues, eigenvalues);
                CheckEqual(expectedEigenvectors, eigenvectors);
            }

            public void FullSubdiagonal2()
            {
                DenseVector diagonal = new Scalar[] { 1, 2, };
                DenseVector subDiagonal = new Scalar[] { 3 };

                var tridiagonal = new SymmetricTridiagonal(ref diagonal, ref subDiagonal);

                var eigenvectors = new SquareMatrix(2);
                var eigenvalues = new DenseVector(2);

                eigenvectors.SetIdentity();
                tridiagonal.ComputeEigenValuesAndVectors(ref eigenvalues, ref eigenvectors);

                DenseVector expectedEigenvalues = new Scalar[]
                {
                    (3 - Math.Sqrt(37)) / 2.0,
                    (3 + Math.Sqrt(37)) / 2.0,                    
                };

                SquareMatrix expectedEigenvectors = new Scalar[,] 
                {
                    { 0.763019982472726, 0.646374896130196, },                    
                    { -0.646374896130196, 0.763019982472726, },                    
                };

                CheckEqual(expectedEigenvalues, eigenvalues);
                CheckEqual(expectedEigenvectors, eigenvectors);
            }

            public void FullSubdiagonal3()
            {
                DenseVector diagonal = new Scalar[] { 1, 2, 3, };
                DenseVector subDiagonal = new Scalar[] { 3, 2 };

                var tridiagonal = new SymmetricTridiagonal(ref diagonal, ref subDiagonal);

                var eigenvectors = new SquareMatrix(3);
                var eigenvalues = new DenseVector(3);

                eigenvectors.SetIdentity();
                tridiagonal.ComputeEigenValuesAndVectors(ref eigenvalues, ref eigenvectors);

                DenseVector expectedEigenvalues = new Scalar[]
                {
                    -1.908854944928920,
                    2.360489028035530,
                    5.548365916893390,
                };

                SquareMatrix expectedEigenvectors = new Scalar[,] 
                {
                    { 0.690686942305700, 0.557527622245240, 0.460558897603740 },
                    { -0.669702709174591, 0.252836737630463, 0.698263464194281 },
                    { 0.272854959736150, -0.790718998467821, 0.548008792273838 },
                };

                CheckEqual(expectedEigenvalues, eigenvalues);
                CheckEqual(expectedEigenvectors, eigenvectors);
            }

            public void FullSubdiagonal4()
            {
                DenseVector diagonal = new Scalar[] { 1, 2, 3, 4, };
                DenseVector subDiagonal = new Scalar[] { 2, 4, 8, };

                var tridiagonal = new SymmetricTridiagonal(ref diagonal, ref subDiagonal);

                var eigenvectors = new SquareMatrix(4);
                var eigenvalues = new DenseVector(4);

                eigenvectors.SetIdentity();
                tridiagonal.ComputeEigenValuesAndVectors(ref eigenvalues, ref eigenvectors);

                DenseVector expectedEigenvalues = new Scalar[]
                {
                    -0.155418156176762,
                    3.620747630960274, 
                    -5.774652139314513,
                    12.309322664530997,                   
                };

                SquareMatrix expectedEigenvectors = new Scalar[,] 
                {
                    { 0.8161350555, 0.5639586241, 0.1160284037, 0.0491085634 },
                    { -0.4714886305, 0.738996614, -0.3930260366, 0.2776922943 },
                    { -0.1540037391, 0.0174524408, 0.7058959772, 0.6911505842 },
                    { 0.2964875896, -0.3681441118, -0.5777359375, 0.6654218276 },
                };

                CheckEqual(expectedEigenvalues, eigenvalues);
                CheckEqual(expectedEigenvectors, eigenvectors);
            }
        }
    }
}