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

namespace Azimuth.UnitTests.SparseLinearAlgebra
{
    public class TridiagonalTests : UnitTestSharp.TestFixture
    {
        public void Basic()
        {
            var tri = new Tridiagonal(
                new DenseVector(new Scalar[] { 1, 1, 1 }),
                new DenseVector(new Scalar[] { 1, 2, 2, 2 }),
                new DenseVector(new Scalar[] { 1, 1, 1 })
            );

            var lower = new CoreTransformation[3];
            var upper = new CoreTransformation[3];
            bool success = tri.Decompose(ref lower, ref upper);

            Check(success);

            CheckEqual(new CoreTransformation(new Matrix2x2(1, 0, 1, 1), 0), lower[0]);
            CheckEqual(new CoreTransformation(new Matrix2x2(1, 0, 1, 1), 1), lower[1]);
            CheckEqual(new CoreTransformation(new Matrix2x2(1, 0, 1, 1), 2), lower[2]);

            CheckEqual(new CoreTransformation(new Matrix2x2(1, 1, 0, 1), 2), upper[0]);
            CheckEqual(new CoreTransformation(new Matrix2x2(1, 1, 0, 1), 1), upper[1]);
            CheckEqual(new CoreTransformation(new Matrix2x2(1, 1, 0, 1), 0), upper[2]);
        }

        public void UpperBidiagonalAlready()
        {
            var tri = new Tridiagonal(
                new DenseVector(new Scalar[] { 1, 1, 1 }),
                new DenseVector(new Scalar[] { 2, 2, 2, 2 }),
                new DenseVector(new Scalar[] { 0, 0, 0 })
            );

            var lower = new CoreTransformation[3];
            var upper = new CoreTransformation[3];
            bool success = tri.Decompose(ref lower, ref upper);

            Check(success);

            CheckEqual(new CoreTransformation(Matrix2x2.Identity, 0), lower[0]);
            CheckEqual(new CoreTransformation(Matrix2x2.Identity, 1), lower[1]);
            CheckEqual(new CoreTransformation(Matrix2x2.Identity, 2), lower[2]);

            CheckEqual(new CoreTransformation(new Matrix2x2(2, 1, 0, 2), 2), upper[0]);
            CheckEqual(new CoreTransformation(new Matrix2x2(2, 1, 0, 1), 1), upper[1]);
            CheckEqual(new CoreTransformation(new Matrix2x2(2, 1, 0, 1), 0), upper[2]);
        }

        public void LowerBidiagonalAlready()
        {
            var tri = new Tridiagonal(
                new DenseVector(new Scalar[] { 0, 0, 0 }),
                new DenseVector(new Scalar[] { 2, 2, 2, 2 }),
                new DenseVector(new Scalar[] { 1, 1, 1 })
            );

            var lower = new CoreTransformation[3];
            var upper = new CoreTransformation[3];
            bool success = tri.Decompose(ref lower, ref upper);

            Check(success);

            CheckEqual(new CoreTransformation(new Matrix2x2(1, 0, 0.5, 1), 0), lower[0]);
            CheckEqual(new CoreTransformation(new Matrix2x2(1, 0, 0.5, 1), 1), lower[1]);
            CheckEqual(new CoreTransformation(new Matrix2x2(1, 0, 0.5, 1), 2), lower[2]);

            CheckEqual(new CoreTransformation(new Matrix2x2(2, 0, 0, 2), 2), upper[0]);
            CheckEqual(new CoreTransformation(new Matrix2x2(2, 0, 0, 1), 1), upper[1]);
            CheckEqual(new CoreTransformation(new Matrix2x2(2, 0, 0, 1), 0), upper[2]);
        }

        public void SingularPrincipalSubmatrix()
        {
            var tri = new Tridiagonal(
                new DenseVector(new Scalar[] { 1, 2, 3 }),
                new DenseVector(new Scalar[] { 4, 5, 6, 7 }),
                new DenseVector(new Scalar[] { 8, 9, 10 })
            );

            var lower = new CoreTransformation[3];
            var upper = new CoreTransformation[3];
            bool success = tri.Decompose(ref lower, ref upper);

            CheckFalse(success);
        }

        public void Diagonal()
        {
            var tri = new Tridiagonal(
                new DenseVector(new Scalar[] { 0, 0, 0 }),
                new DenseVector(new Scalar[] { 1, 2, 3, 4 }),
                new DenseVector(new Scalar[] { 0, 0, 0 })
            );

            var lower = new CoreTransformation[3];
            var upper = new CoreTransformation[3];
            bool success = tri.Decompose(ref lower, ref upper);

            Check(success);

            CheckEqual(new CoreTransformation(Matrix2x2.Identity, 0), lower[0]);
            CheckEqual(new CoreTransformation(Matrix2x2.Identity, 1), lower[1]);
            CheckEqual(new CoreTransformation(Matrix2x2.Identity, 2), lower[2]);

            CheckEqual(new CoreTransformation(new Matrix2x2(3, 0, 0, 4), 2), upper[0]);
            CheckEqual(new CoreTransformation(new Matrix2x2(2, 0, 0, 1), 1), upper[1]);
            CheckEqual(new CoreTransformation(new Matrix2x2(1, 0, 0, 1), 0), upper[2]);
        }

        public void Identity()
        {
            var tri = new Tridiagonal(
                new DenseVector(new Scalar[] { 0, 0, 0 }),
                new DenseVector(new Scalar[] { 1, 1, 1, 1 }),
                new DenseVector(new Scalar[] { 0, 0, 0 })
            );

            var lower = new CoreTransformation[3];
            var upper = new CoreTransformation[3];
            bool success = tri.Decompose(ref lower, ref upper);

            Check(success);

            CheckEqual(new CoreTransformation(Matrix2x2.Identity, 0), lower[0]);
            CheckEqual(new CoreTransformation(Matrix2x2.Identity, 1), lower[1]);
            CheckEqual(new CoreTransformation(Matrix2x2.Identity, 2), lower[2]);

            CheckEqual(new CoreTransformation(Matrix2x2.Identity, 2), upper[0]);
            CheckEqual(new CoreTransformation(Matrix2x2.Identity, 1), upper[1]);
            CheckEqual(new CoreTransformation(Matrix2x2.Identity, 0), upper[2]);
        }

        public void LotsOfNumbers()
        {
            var tri = new Tridiagonal(
                new DenseVector(new Scalar[] { -2, 4, -2 }),
                new DenseVector(new Scalar[] { 2, -3, 14, 24 }),
                new DenseVector(new Scalar[] { 6, 6, -42 })
            );

            var lower = new CoreTransformation[3];
            var upper = new CoreTransformation[3];
            bool success = tri.Decompose(ref lower, ref upper);

            Check(success);

            CheckEqual(new CoreTransformation(new Matrix2x2(1, 0,  3, 1), 0), lower[0]);
            CheckEqual(new CoreTransformation(new Matrix2x2(1, 0,  2, 1), 1), lower[1]);
            CheckEqual(new CoreTransformation(new Matrix2x2(1, 0, -7, 1), 2), lower[2]);

            CheckEqual(new CoreTransformation(new Matrix2x2(6, -2, 0, 10), 2), upper[0]); 
            CheckEqual(new CoreTransformation(new Matrix2x2(3,  4, 0,  1), 1), upper[1]);
            CheckEqual(new CoreTransformation(new Matrix2x2(2, -2, 0,  1), 0), upper[2]);
        }

        [IgnoreTest]
        public void FauxChebyshev()
        {
            // pre balanced matrix
            var tri = new Tridiagonal(
                new DenseVector(new Scalar[] { 1, 0.5, 0.5, }),
                new DenseVector(new Scalar[] { 1, 1, 1, 4 }),
                new DenseVector(new Scalar[] { 0.5, 0.5, 2 })
            );

            var lower = new CoreTransformation[4];
            var upper = new CoreTransformation[4];
            bool success = tri.Decompose(ref lower, ref upper);

            Check(success);

            CheckEqual(new CoreTransformation(new Matrix2x2(1, 0, 3, 1), 0), lower[0]);
            CheckEqual(new CoreTransformation(new Matrix2x2(1, 0, 2, 1), 1), lower[1]);
            CheckEqual(new CoreTransformation(new Matrix2x2(1, 0, -7, 1), 2), lower[2]);

            CheckEqual(new CoreTransformation(new Matrix2x2(6, -2, 0, 10), 2), upper[0]);
            CheckEqual(new CoreTransformation(new Matrix2x2(3, 4, 0, 1), 1), upper[1]);
            CheckEqual(new CoreTransformation(new Matrix2x2(2, -2, 0, 1), 0), upper[2]);
        }
    }
}