using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using UnitTestSharp;
using Azimuth.DenseLinearAlgebra;

namespace Azimuth.UnitTests.DenseLinearAlgebra
{
    /// <summary>
    /// Notation:
    /// PD = Positive Definite
    /// SPD = Semi positive definite
    /// H = Hermitian = Symmetric
    /// </summary>
    public class CholeskyDecompositionTests : TestFixture
    {
        public class FullRankPDH : TestFixture
        {
            public void x11()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    {7},
                });

                var expected_L = new SquareMatrix(new Scalar[,] {
                    {1},
                });

                var expected_D = new DenseVector(new Scalar[] {
                    7,
                });

                CholeskyDecomposition decomp = new CholeskyDecomposition(matrix);

                CheckEqual(1, decomp.Rank);
                CheckEqual(1, decomp.Size);
                CheckEqual(decomp.Size, decomp.ColumnCount);
                CheckEqual(decomp.Size, decomp.RowCount);
                
                CheckEqual(expected_L, decomp.L);
                CheckEqual(expected_D, decomp.D);                
            }

            public void x22()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    {5, 2,},
                    {2, 6},
                });

                Scalar sqr5 = Math.Sqrt(5);

                var expected_L = new SquareMatrix(new Scalar[,] {
                    {1,          0, },
                    {2.0 / 5.0,  1 },
                });

                var expected_D = new DenseVector(new Scalar[] {
                    5,
                    5.2,
                });

                CholeskyDecomposition decomp = new CholeskyDecomposition(matrix);

                CheckEqual(2, decomp.Rank);
                CheckEqual(2, decomp.Size);
                CheckEqual(decomp.Size, decomp.ColumnCount);
                CheckEqual(decomp.Size, decomp.RowCount);

                CheckEqual(expected_L, decomp.L);
                CheckEqual(expected_D, decomp.D);            
            }

            public void x33()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    {2, 1, 1,},
                    {1, 2, 1,},
                    {1, 1, 2,},
                });

                var expected_L = new SquareMatrix(new Scalar[,] {
                    { 1,         0, 0, },
                    { 0.5,       1, 0, },
                    { 0.5, 1.0/3.0, 1, },
                });

                var expected_D = new DenseVector(new Scalar[] {
                    2,
                    1.5,
                    4.0/3.0,
                });

                CholeskyDecomposition decomp = new CholeskyDecomposition(matrix);

                CheckEqual(3, decomp.Rank);
                CheckEqual(3, decomp.Size);
                CheckEqual(decomp.Size, decomp.ColumnCount);
                CheckEqual(decomp.Size, decomp.RowCount);

                CheckEqual(expected_L, decomp.L);
                CheckEqual(expected_D, decomp.D);            
            }

            public void PremultiplyByInverse_x11()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    {7},
                });

                CholeskyDecomposition decomp = new CholeskyDecomposition(matrix);

                var output = new DenseVector(new Scalar[] { 7.0 });
                decomp.PremultiplyByInverse(ref output);

                CheckEqual(1, decomp.Rank);
                CheckEqual(new DenseVector(new Scalar[] { 1 }), output);
            }

            public void PremultiplyByInverse_x22()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    {5, 2,},
                    {2, 6},
                });

                var output = new DenseVector(new Scalar[] { 2.0, 3.0 });

                CholeskyDecomposition decomp = new CholeskyDecomposition(matrix);                
                decomp.PremultiplyByInverse(ref output);

                CheckEqual(new DenseVector(new Scalar[] { 3.0/13, 11.0/26 }), output);
            }

            public void PremultiplyByInverse_x33()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    {2, 1, 1,},
                    {1, 2, 1,},
                    {1, 1, 2,},
                });

                var output = new DenseVector(new Scalar[] { 2.0, 3.0, 5.0 });

                CholeskyDecomposition decomp = new CholeskyDecomposition(matrix);
                decomp.PremultiplyByInverse(ref output);

                CheckEqual(new DenseVector(new Scalar[] { -.5, .5, 2.5 }), output);
            }
        }

        public class RankDeficientPDH : TestFixture
        {
            public void x22()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    {1, 2, },
                    {2, 4, },
                });

                var expected_L = new SquareMatrix(new Scalar[,] {
                    {1, 0, },
                    {2, 1, },
                });

                var expected_D = new DenseVector(new Scalar[] {
                    1, 
                    0,
                });

                CholeskyDecomposition decomp = new CholeskyDecomposition(matrix);

                CheckEqual(1, decomp.Rank);
                CheckEqual(2, decomp.Size);
                CheckEqual(decomp.Size, decomp.ColumnCount);
                CheckEqual(decomp.Size, decomp.RowCount);

                CheckEqual(expected_L, decomp.L);
                CheckEqual(expected_D, decomp.D);
            }

            public void x33()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    {1, 0, 1,},
                    {0, 1, 1,},
                    {1, 1, 2,},
                });

                var expected_L = new SquareMatrix(new Scalar[,] {
                    {1, 0, 0,},
                    {0, 1, 0,},
                    {1, 1, 1,},
                });

                var expected_D = new DenseVector(new Scalar[] {
                    1,
                    1,
                    0,
                });

                CholeskyDecomposition decomp = new CholeskyDecomposition(matrix);

                CheckEqual(2, decomp.Rank);
                CheckEqual(3, decomp.Size);
                CheckEqual(decomp.Size, decomp.ColumnCount);
                CheckEqual(decomp.Size, decomp.RowCount);

                CheckEqual(expected_L, decomp.L);
                CheckEqual(expected_D, decomp.D);
            }

            public void x33_Permutation()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    {1, 1, 0, },
                    {1, 2, 1, },
                    {0, 1, 1, },
                });

                var expected_L = new SquareMatrix(new Scalar[,] {
                    {1, 0, 0,},
                    {1, 1, 0,},
                    {0, 1, 1,},
                });

                var expected_D = new DenseVector(new Scalar[] {
                    1,
                    1,
                    0,
                });

                CholeskyDecomposition decomp = new CholeskyDecomposition(matrix);

                CheckEqual(2, decomp.Rank);
                CheckEqual(3, decomp.Size);
                CheckEqual(decomp.Size, decomp.ColumnCount);
                CheckEqual(decomp.Size, decomp.RowCount);

                CheckEqual(expected_L, decomp.L);
                CheckEqual(expected_D, decomp.D);
            }

            //[IgnoreTest]
            public void x88()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                       { 93,    62,    67,    71,    57,   107,   100,    48, },
                       { 62,    49,    50,    52,    33,    79,    74,    42, },
                       { 67,    50,    58,    65,    42,    89,    91,    42, },
                       { 71,    52,    65,   109,    67,    99,   120,    52, },
                       { 57,    33,    42,    67,    52,    66,    75,    28, },
                       {107,    79,    89,    99,    66,   138,   138,    66, },
                       {100,    74,    91,   120,    75,   138,   153,    66, },
                       { 48,    42,    42,    52,    28,    66,    66,    40, },
                });

                var expected_L = new SquareMatrix(new Scalar[,] {
                    { 1, 0, 0, 0, 0, 0, 0, 0 },
                    { 0.6666666667, 1, 0, 0, 0, 0, 0, 0 },
                    { 0.7204301075, 0.6956521739, 1, 0, 0, 0, 0, 0 },
                    { 0.7634408602, 0.6086956522, 1.7610062893, 1, 0, 0, 0, 0 },
                    { 0.6129032258, -0.652173913, 0.7330538085, 0.5634920635, 1, 0, 0, 0 },
                    { 1.1505376344, 1, 1.0929419986, 0.0317460317, 66, 1, 0, 0 },
                    { 1.0752688172, 0.9565217391, 2.3011879804, 0.4444444444, 75, 138, 1, 0 },
                    { 0.5161290323, 1.3043478261, 0.0768693222, 0.253968254, 28, 66, 66, 1 },
                });

                var expected_D = new DenseVector(new Scalar[] {
                    93, 7.66666666666667, 6.02103786816269, 33.2830188679245, 
                    0, 0, 0, 0
                });

                CholeskyDecomposition decomp = new CholeskyDecomposition(matrix);

                CheckEqual(4, decomp.Rank);
                CheckEqual(8, decomp.Size);
                CheckEqual(decomp.Size, decomp.ColumnCount);
                CheckEqual(decomp.Size, decomp.RowCount);

                CheckEqual(expected_L, decomp.L);
                CheckEqual(expected_D, decomp.D);
            }
        }

        public class FullRankPDNonH : TestFixture
        {
            public void x22()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    {1, 1, },
                    {-1, 1, },
                });

                CheckThrow(typeof(Exception));
                CholeskyDecomposition decomp = new CholeskyDecomposition(matrix);
            }
        }

        public class NotPSD : TestFixture
        {
            public void x22()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    {0, 1, },
                    {1, 0, },
                });

                var expected_L = new SquareMatrix(new Scalar[,] {
                    {1, 0,},
                    {1, 1,},
                });

                var expected_D = new DenseVector(new Scalar[] {
                    0,
                    0,
                });

                CholeskyDecomposition decomp = new CholeskyDecomposition(matrix);

                CheckEqual(0, decomp.Rank);
                CheckEqual(expected_L, decomp.L);
                CheckEqual(expected_D, decomp.D);
            }

            public void x33()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { 1, 0, 0, },
                    { 0, 0, 1, },
                    { 0, 1, 0, },
                });

                var expected_L = new SquareMatrix(new Scalar[,] {
                    { 1, 0, 0, },
                    { 0, 1, 0, },
                    { 0, 1, 1, },
                });

                var expected_D = new DenseVector(new Scalar[] {
                    1,
                    0,
                    0,
                });

                CholeskyDecomposition decomp = new CholeskyDecomposition(matrix);

                CheckEqual(1, decomp.Rank);
                CheckEqual(expected_L, decomp.L);
                CheckEqual(expected_D, decomp.D);
            }
        }

        public class SizeEtAl : TestFixture
        {
            SquareMatrix matrix = new Scalar[,] {
                    {1, 2, 3},
                    {2, 5, 6},
                    {3, 6, 9},
                };

            CholeskyDecomposition cholesky = null;

            public override void FixtureSetup()
            {
 	             cholesky = new CholeskyDecomposition(matrix);
            }

            public void Size()
            {
                CheckEqual(3, cholesky.Size);
            }

            public void RowCount()
            {
                CheckEqual(3, cholesky.RowCount);
            }

            public void ColumnCount()
            {
                CheckEqual(3, cholesky.ColumnCount);
            }
        }

        public class PremultiplyByInverse : TestFixture
        {
            public void x22()
            {
                SquareMatrix matrix = new Scalar[,] {
                    { 75, 15, },
                    { 15, 9, },
                };

                var input = new DenseVector(new Scalar[] {
                    420,
                    102,
                });

                var expected = new DenseVector(new Scalar[] {
                    5,
                    3,
                });

                var cholseky = new CholeskyDecomposition(matrix);

                var actual = input.Clone();

                cholseky.PremultiplyByInverse(ref actual);

                CheckEqual(expected, actual);
            }

            // PSD matrices don't have well defined inverses from a Cholesky decomposition
            [IgnoreTest]
            public void xPSD()
            {
                SquareMatrix matrix = new Scalar[,] {
                    { 10, 8, 6, },
                    {  8, 8, 8, },
                    {  6, 8, 10, },
                };

                var input = new DenseVector(new Scalar[] {
                    3,
                    4,
                    5,
                });

                var expected = new DenseVector(new Scalar[] {
                    1.0/12,
                    1.0/6,
                    1.0/2.4,
                });

                var cholseky = new CholeskyDecomposition(matrix);

                var actual = input.Clone();

                cholseky.PremultiplyByInverse(ref actual);

                CheckEqual(expected, actual);
            }
        }

        public class ExplicitlyCalculateInverseTests : TestFixture
        {
            public void x11()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { 4, },
                });

                var inverse = new SquareMatrix(new Scalar[,] {
                    { 0.25, },
                });

                CholeskyDecomposition decomp = new CholeskyDecomposition(matrix);

                CheckEqual(inverse, decomp.ExplicitlyCalculateInverse());
            }
            
            public void x22()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    {2, 1, },
                    {1, 2, },
                });

                var inverse  = new SquareMatrix(new Scalar[,] {
                    {  2.0/3.0, -1.0/3.0, },
                    { -1.0/3.0,  2.0/3.0, },
                });

                CholeskyDecomposition decomp = new CholeskyDecomposition(matrix);

                CheckEqual(inverse, decomp.ExplicitlyCalculateInverse());
            }

            public void x33()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { 1, 1, 1, },
                    { 1, 2, 2, },
                    { 1, 2, 3, },
                });

                var inverse = new SquareMatrix(new Scalar[,] {
                    {  2, -1,  0, },
                    { -1,  2, -1, },
                    {  0, -1,  1, },
                });

                CholeskyDecomposition decomp = new CholeskyDecomposition(matrix);

                CheckEqual(inverse, decomp.ExplicitlyCalculateInverse());
            }
        }
    }
}