using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Azimuth.UnitTests.RootFinding
{
    public class BisectionTests : UnitTestSharp.TestFixture
    {
        public class Polynomials : UnitTestSharp.TestFixture
        {
            public void Basic()
            {
                Scalar r1 = Azimuth.RootFinding.Bisection.FindRoot(x => x*(x+18)+17, new Interval(-400, -2));
                Scalar r2 = Azimuth.RootFinding.Bisection.FindRoot(x => x * (x + 18) + 17, new Interval(-2, 7000));

                CheckEqual(-17, r1);
                CheckEqual(-1, r2);
            }

            public void DoesNotStraddle()
            {
                CheckThrow(typeof(Exception));
                Scalar r1 = Azimuth.RootFinding.Bisection.FindRoot(x => x * (x + 18) + 17, new Interval(-400, -401));
            }

            public void MultipleRoots()
            {
                // Should give us the minimum root
                Scalar r1 = Azimuth.RootFinding.Bisection.FindRoot(
                    x => (x + 3) * (x - 7) * (x - 4) * (x - 4) * (x - 4) * (x - 4) * (x + 2) * (x - 1) * (x - 3),
                    new Interval(-400, 400));

                var possibleRoots = new List<Scalar>(new Scalar[] { -3, 7, 4, -2, 1, 3 });

                // There's no way to be sure which root the bisection will converge to.
                Check(possibleRoots.Contains(r1));
            }

            public void InfinityBounds()
            {
                CheckThrow(typeof(Exception));
                Scalar r1 = Azimuth.RootFinding.Bisection.FindRoot(
                    x => (x + 3) * (x - 7) * (x - 4) * (x - 4) * (x - 4) * (x - 4) * (x + 2) * (x - 1) * (x - 3),
                    new Interval(Scalar.NegativeInfinity, Scalar.PositiveInfinity));
            }

            public void NonPolynomial()
            {
                var func = new Func<Scalar, Scalar>(x => Math.Exp(x) - 0.5);

                Scalar r1 = Azimuth.RootFinding.Bisection.FindRoot(func, new Interval(-10, 10));

                CheckEqual(Math.Log(0.5), r1);
            }

            public void Singularity()
            {
                var func = new Func<Scalar, Scalar>(x => 1.0 / (x - 1.1));

                Scalar r1 = Azimuth.RootFinding.Bisection.FindRoot(func, new Interval(-10, 10));

                CheckEqual(1.1, r1);
            }

            public void ComplexEval()
            {
                // Root is at 2
                var initialLineOrientation = Math.PI / 4;
                var lineAngularVelocity = Math.PI / 8;
                var ptLeverArm = new Vector(10, 10);
                var ptLinearAcceleration = new Vector(0, -5);
                var ptLinearVelocity = new Vector(0, -5);
                var ptPos = new Vector(-10, 0);
                var dw = -Math.PI / 4;

                Func<Scalar, Scalar> derivative = (t) =>
                {
                    Scalar dAngle;
                    {
                        var dAngle_accum = new RobustAccumulator();
                        dAngle_accum.Add(dw);
                        dAngle_accum.Mul(t);
                        dAngle = dAngle_accum.ComputedSum;
                    }

                    Scalar angle1;
                    {
                        var angle1_accum = new RobustAccumulator();
                        angle1_accum.AddQuadratic(0, lineAngularVelocity, initialLineOrientation, t);
                        angle1 = angle1_accum.ComputedSum;
                    }

                    var w1 = lineAngularVelocity;

                    var displacementAccumulator_X = new RobustAccumulator();
                    displacementAccumulator_X.AddQuadratic(
                        -0.5 * ptLinearAcceleration.X,
                        -ptLinearVelocity.X,
                        -ptPos.X,
                        t);
                    displacementAccumulator_X.Mul(w1);
                    displacementAccumulator_X.AddQuadratic(
                        0,
                        -ptLinearAcceleration.X,
                        -ptLinearVelocity.X,
                        t);

                    var displacementAccumulator_Y = new RobustAccumulator();
                    displacementAccumulator_Y.AddQuadratic(
                        -0.5 * ptLinearAcceleration.Y,
                        -ptLinearVelocity.Y,
                        -ptPos.Y,
                        t);
                    displacementAccumulator_Y.Mul(w1);
                    displacementAccumulator_X.AddQuadratic(
                        0,
                        ptLinearAcceleration.Y,
                        ptLinearVelocity.Y,
                        t);

                    var accum = new RobustAccumulator();

                    accum.AddHarmonic(ptLeverArm.X * dw, -ptLeverArm.Y * dw, dAngle);
                    accum.AddHarmonic(displacementAccumulator_X, displacementAccumulator_Y, angle1);

                    return accum.ComputedSum;
                };

                int iterations;
                Scalar r1 = Azimuth.RootFinding.Bisection.FindRoot(derivative, new Interval(1.9, 2.1), out iterations);

                CheckEqual(2.0, r1);
                Check(64 >= iterations);
            }
        }
    }
}