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 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 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);
            }
        }

        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());
            }
        }
    }
}