// <copyright file="NewtonCotesTrapeziumRule.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
//
// Copyright (c) 2009-2013 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>

using System;
using System.Collections.Generic;
using System.Numerics;

namespace IStation.Numerics.Integration
{
    /// <summary>
    /// Approximation algorithm for definite integrals by the Trapezium rule of the Newton-Cotes family.
    /// </summary>
    /// <remarks>
    /// <a href="http://en.wikipedia.org/wiki/Trapezium_rule">Wikipedia - Trapezium Rule</a>
    /// </remarks>
    public static class NewtonCotesTrapeziumRule
    {
        /// <summary>
        /// Direct 2-point approximation of the definite integral in the provided interval by the trapezium rule.
        /// </summary>
        /// <param name="f">The analytic smooth function to integrate.</param>
        /// <param name="intervalBegin">Where the interval starts, inclusive and finite.</param>
        /// <param name="intervalEnd">Where the interval stops, inclusive and finite.</param>
        /// <returns>Approximation of the finite integral in the given interval.</returns>
        public static double IntegrateTwoPoint(Func<double, double> f, double intervalBegin, double intervalEnd)
        {
            if (f == null)
            {
                throw new ArgumentNullException(nameof(f));
            }

            return (intervalEnd - intervalBegin)/2*(f(intervalBegin) + f(intervalEnd));
        }

        /// <summary>
        /// Direct 2-point approximation of the definite integral in the provided interval by the trapezium rule.
        /// </summary>
        /// <param name="f">The analytic smooth complex function to integrate, defined on real domain.</param>
        /// <param name="intervalBegin">Where the interval starts, inclusive and finite.</param>
        /// <param name="intervalEnd">Where the interval stops, inclusive and finite.</param>
        /// <returns>Approximation of the finite integral in the given interval.</returns>
        public static Complex ContourIntegrateTwoPoint(Func<double, Complex> f, double intervalBegin, double intervalEnd)
        {
            if (f == null)
            {
                throw new ArgumentNullException(nameof(f));
            }

            return (intervalEnd - intervalBegin) / 2 * (f(intervalBegin) + f(intervalEnd));
        }

        /// <summary>
        /// Composite N-point approximation of the definite integral in the provided interval by the trapezium rule.
        /// </summary>
        /// <param name="f">The analytic smooth function to integrate.</param>
        /// <param name="intervalBegin">Where the interval starts, inclusive and finite.</param>
        /// <param name="intervalEnd">Where the interval stops, inclusive and finite.</param>
        /// <param name="numberOfPartitions">Number of composite subdivision partitions.</param>
        /// <returns>Approximation of the finite integral in the given interval.</returns>
        public static double IntegrateComposite(Func<double, double> f, double intervalBegin, double intervalEnd, int numberOfPartitions)
        {
            if (f == null)
            {
                throw new ArgumentNullException(nameof(f));
            }

            if (numberOfPartitions <= 0)
            {
                throw new ArgumentOutOfRangeException(nameof(numberOfPartitions), "Value must be positive (and not zero).");
            }

            double step = (intervalEnd - intervalBegin)/numberOfPartitions;

            double offset = step;
            double sum = 0.5*(f(intervalBegin) + f(intervalEnd));
            for (int i = 0; i < numberOfPartitions - 1; i++)
            {
                // NOTE (ruegg, 2009-01-07): Do not combine intervalBegin and offset (numerical stability!)
                sum += f(intervalBegin + offset);
                offset += step;
            }

            return step*sum;
        }

        /// <summary>
        /// Composite N-point approximation of the definite integral in the provided interval by the trapezium rule.
        /// </summary>
        /// <param name="f">The analytic smooth complex function to integrate, defined on real domain.</param>
        /// <param name="intervalBegin">Where the interval starts, inclusive and finite.</param>
        /// <param name="intervalEnd">Where the interval stops, inclusive and finite.</param>
        /// <param name="numberOfPartitions">Number of composite subdivision partitions.</param>
        /// <returns>Approximation of the finite integral in the given interval.</returns>
        public static Complex ContourIntegrateComposite(Func<double, Complex> f, double intervalBegin, double intervalEnd, int numberOfPartitions)
        {
            if (f == null)
            {
                throw new ArgumentNullException(nameof(f));
            }

            if (numberOfPartitions <= 0)
            {
                throw new ArgumentOutOfRangeException(nameof(numberOfPartitions), "Value must be positive (and not zero).");
            }

            double step = (intervalEnd - intervalBegin) / numberOfPartitions;

            double offset = step;
            Complex sum = 0.5 * (f(intervalBegin) + f(intervalEnd));
            for (int i = 0; i < numberOfPartitions - 1; i++)
            {
                // NOTE (ruegg, 2009-01-07): Do not combine intervalBegin and offset (numerical stability!)
                sum += f(intervalBegin + offset);
                offset += step;
            }

            return step * sum;
        }

        /// <summary>
        /// Adaptive approximation of the definite integral in the provided interval by the trapezium rule.
        /// </summary>
        /// <param name="f">The analytic smooth function to integrate.</param>
        /// <param name="intervalBegin">Where the interval starts, inclusive and finite.</param>
        /// <param name="intervalEnd">Where the interval stops, inclusive and finite.</param>
        /// <param name="targetError">The expected accuracy of the approximation.</param>
        /// <returns>Approximation of the finite integral in the given interval.</returns>
        public static double IntegrateAdaptive(Func<double, double> f, double intervalBegin, double intervalEnd, double targetError)
        {
            if (f == null)
            {
                throw new ArgumentNullException(nameof(f));
            }

            int numberOfPartitions = 1;
            double step = intervalEnd - intervalBegin;
            double sum = 0.5*step*(f(intervalBegin) + f(intervalEnd));
            for (int k = 0; k < 20; k++)
            {
                double midpointsum = 0;
                for (int i = 0; i < numberOfPartitions; i++)
                {
                    midpointsum += f(intervalBegin + ((i + 0.5)*step));
                }

                midpointsum *= step;
                sum = 0.5*(sum + midpointsum);
                step *= 0.5;
                numberOfPartitions *= 2;

                if (sum.AlmostEqualRelative(midpointsum, targetError))
                {
                    break;
                }
            }

            return sum;
        }

        /// <summary>
        /// Adaptive approximation of the definite integral in the provided interval by the trapezium rule.
        /// </summary>
        /// <param name="f">The analytic smooth complex function to integrate, define don real domain.</param>
        /// <param name="intervalBegin">Where the interval starts, inclusive and finite.</param>
        /// <param name="intervalEnd">Where the interval stops, inclusive and finite.</param>
        /// <param name="targetError">The expected accuracy of the approximation.</param>
        /// <returns>Approximation of the finite integral in the given interval.</returns>
        public static Complex ContourIntegrateAdaptive(Func<double, Complex> f, double intervalBegin, double intervalEnd, double targetError)
        {
            if (f == null)
            {
                throw new ArgumentNullException(nameof(f));
            }

            int numberOfPartitions = 1;
            double step = intervalEnd - intervalBegin;
            Complex sum = 0.5 * step * (f(intervalBegin) + f(intervalEnd));
            for (int k = 0; k < 20; k++)
            {
                Complex midpointsum = 0;
                for (int i = 0; i < numberOfPartitions; i++)
                {
                    midpointsum += f(intervalBegin + ((i + 0.5) * step));
                }

                midpointsum *= step;
                sum = 0.5 * (sum + midpointsum);
                step *= 0.5;
                numberOfPartitions *= 2;

                if (sum.AlmostEqualRelative(midpointsum, targetError))
                {
                    break;
                }
            }

            return sum;
        }


        /// <summary>
        /// Adaptive approximation of the definite integral by the trapezium rule.
        /// </summary>
        /// <param name="f">The analytic smooth function to integrate.</param>
        /// <param name="intervalBegin">Where the interval starts, inclusive and finite.</param>
        /// <param name="intervalEnd">Where the interval stops, inclusive and finite.</param>
        /// <param name="levelAbscissas">Abscissa vector per level provider.</param>
        /// <param name="levelWeights">Weight vector per level provider.</param>
        /// <param name="levelOneStep">First Level Step</param>
        /// <param name="targetRelativeError">The expected relative accuracy of the approximation.</param>
        /// <returns>Approximation of the finite integral in the given interval.</returns>
        public static double IntegrateAdaptiveTransformedOdd(
            Func<double, double> f,
            double intervalBegin, double intervalEnd,
            IEnumerable<double[]> levelAbscissas, IEnumerable<double[]> levelWeights,
            double levelOneStep, double targetRelativeError)
        {
            if (f == null)
            {
                throw new ArgumentNullException(nameof(f));
            }

            if (levelAbscissas == null)
            {
                throw new ArgumentNullException(nameof(levelAbscissas));
            }

            if (levelWeights == null)
            {
                throw new ArgumentNullException(nameof(levelWeights));
            }

            double linearSlope = 0.5*(intervalEnd - intervalBegin);
            double linearOffset = 0.5*(intervalEnd + intervalBegin);
            targetRelativeError /= 5*linearSlope;

            using (var abcissasIterator = levelAbscissas.GetEnumerator())
            using (var weightsIterator = levelWeights.GetEnumerator())
            {
                double step = levelOneStep;

                // First Level
                abcissasIterator.MoveNext();
                weightsIterator.MoveNext();
                double[] abcissasL1 = abcissasIterator.Current;
                double[] weightsL1 = weightsIterator.Current;

                double sum = f(linearOffset)*weightsL1[0];
                for (int i = 1; i < abcissasL1.Length; i++)
                {
                    sum += weightsL1[i]*(f((linearSlope*abcissasL1[i]) + linearOffset) + f(-(linearSlope*abcissasL1[i]) + linearOffset));
                }

                sum *= step;

                // Additional Levels
                double previousDelta = double.MaxValue;
                for (int level = 1; abcissasIterator.MoveNext() && weightsIterator.MoveNext(); level++)
                {
                    double[] abcissas = abcissasIterator.Current;
                    double[] weights = weightsIterator.Current;

                    double midpointsum = 0;
                    for (int i = 0; i < abcissas.Length; i++)
                    {
                        midpointsum += weights[i]*(f((linearSlope*abcissas[i]) + linearOffset) + f(-(linearSlope*abcissas[i]) + linearOffset));
                    }

                    midpointsum *= step;
                    sum = 0.5*(sum + midpointsum);
                    step *= 0.5;

                    double delta = Math.Abs(sum - midpointsum);

                    if (level == 1)
                    {
                        previousDelta = delta;
                        continue;
                    }

                    double r = Math.Log(delta)/Math.Log(previousDelta);
                    previousDelta = delta;

                    if (r > 1.9 && r < 2.1)
                    {
                        // convergence region
                        delta = Math.Sqrt(delta);
                    }

                    if (sum.AlmostEqualNormRelative(midpointsum, delta, targetRelativeError))
                    {
                        break;
                    }
                }

                return sum*linearSlope;
            }
        }

        /// <summary>
        /// Adaptive approximation of the definite integral by the trapezium rule.
        /// </summary>
        /// <param name="f">The analytic smooth complex function to integrate, defined on the real domain.</param>
        /// <param name="intervalBegin">Where the interval starts, inclusive and finite.</param>
        /// <param name="intervalEnd">Where the interval stops, inclusive and finite.</param>
        /// <param name="levelAbscissas">Abscissa vector per level provider.</param>
        /// <param name="levelWeights">Weight vector per level provider.</param>
        /// <param name="levelOneStep">First Level Step</param>
        /// <param name="targetRelativeError">The expected relative accuracy of the approximation.</param>
        /// <returns>Approximation of the finite integral in the given interval.</returns>
        public static Complex ContourIntegrateAdaptiveTransformedOdd(
            Func<double, Complex> f,
            double intervalBegin, double intervalEnd,
            IEnumerable<double[]> levelAbscissas, IEnumerable<double[]> levelWeights,
            double levelOneStep, double targetRelativeError)
        {
            if (f == null)
            {
                throw new ArgumentNullException(nameof(f));
            }

            if (levelAbscissas == null)
            {
                throw new ArgumentNullException(nameof(levelAbscissas));
            }

            if (levelWeights == null)
            {
                throw new ArgumentNullException(nameof(levelWeights));
            }

            double linearSlope = 0.5 * (intervalEnd - intervalBegin);
            double linearOffset = 0.5 * (intervalEnd + intervalBegin);
            targetRelativeError /= 5 * linearSlope;

            using (var abcissasIterator = levelAbscissas.GetEnumerator())
            using (var weightsIterator = levelWeights.GetEnumerator())
            {
                double step = levelOneStep;

                // First Level
                abcissasIterator.MoveNext();
                weightsIterator.MoveNext();
                double[] abcissasL1 = abcissasIterator.Current;
                double[] weightsL1 = weightsIterator.Current;

                Complex sum = f(linearOffset) * weightsL1[0];
                for (int i = 1; i < abcissasL1.Length; i++)
                {
                    sum += weightsL1[i] * (f((linearSlope * abcissasL1[i]) + linearOffset) + f(-(linearSlope * abcissasL1[i]) + linearOffset));
                }

                sum *= step;

                // Additional Levels
                double previousDelta = double.MaxValue;
                for (int level = 1; abcissasIterator.MoveNext() && weightsIterator.MoveNext(); level++)
                {
                    double[] abcissas = abcissasIterator.Current;
                    double[] weights = weightsIterator.Current;

                    Complex midpointsum = 0;
                    for (int i = 0; i < abcissas.Length; i++)
                    {
                        midpointsum += weights[i] * (f((linearSlope * abcissas[i]) + linearOffset) + f(-(linearSlope * abcissas[i]) + linearOffset));
                    }

                    midpointsum *= step;
                    sum = 0.5 * (sum + midpointsum);
                    step *= 0.5;

                    double delta = Complex.Abs(sum - midpointsum);

                    if (level == 1)
                    {
                        previousDelta = delta;
                        continue;
                    }

                    double r = Math.Log(delta) / Math.Log(previousDelta);
                    previousDelta = delta;

                    if (r > 1.9 && r < 2.1)
                    {
                        // convergence region
                        delta = Math.Sqrt(delta);
                    }

                    if (sum.Real.AlmostEqualNormRelative(midpointsum.Real, delta, targetRelativeError)
                        && sum.Imaginary.AlmostEqualNormRelative(midpointsum.Imaginary, delta, targetRelativeError))
                    {
                        break;
                    }
                }

                return sum * linearSlope;
            }
        }
    }
}