// <copyright file="StudentT.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-2020 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 IStation.Numerics.Random;
using System;
using System.Collections.Generic;

namespace IStation.Numerics.Distributions
{
    /// <summary>
    /// Continuous Univariate Skewed Generalized Error Distribution (SGED).
    /// Implements the univariate SSkewed Generalized Error Distribution. For details about this
    /// distribution, see
    /// <a href="https://en.wikipedia.org/wiki/Generalized_normal_distribution">
    /// Wikipedia - Generalized Error Distribution</a>.
    /// It includes Laplace, Normal and Student-t distributions.
    /// This is the <see cref="SkewedGeneralizedT"/> distribution with q=Inf.
    /// </summary>
    /// <remarks><para>This implementation is based on the R package dsgt and corresponding viginette, see
    /// <a href="">https://cran.r-project.org/web/packages/sgt/vignettes/sgt.pdf</a>. Compared to that
    /// implementation, the options for mean adjustment and variance adjustment are always true.
    /// The location (μ) is the mean of the distribution.
    /// The scale (σ) squared is the variance of the distribution.
    /// </para>
    /// <para>The distribution will use the <see cref="System.Random"/> by
    /// default.  Users can get/set the random number generator by using the
    /// <see cref="RandomSource"/> property.</para>
    /// <para>The statistics classes will check all the incoming parameters
    /// whether they are in the allowed range.</para></remarks>
    public class SkewedGeneralizedError : IContinuousDistribution
    {
        System.Random _random;

        readonly double _skewness;

        /// <summary>
        /// Initializes a new instance of the SkewedGeneralizedError class. This is a generalized error distribution
        /// with location=0.0, scale=1.0, skew=0.0 and p=2.0 (a standard normal distribution).
        /// </summary>
        public SkewedGeneralizedError()
        {
            _random = SystemRandomSource.Default;
            Location = 0.0;
            Scale = 1.0;
            Skew = 0.0;
            P = 2.0;
        }

        /// <summary>
        /// Initializes a new instance of the SkewedGeneralizedT class with a particular location, scale, skew
        /// and kurtosis parameters. Different parameterizations result in different distributions.
        /// </summary>
        /// <param name="location">The location (μ) of the distribution.</param>
        /// <param name="scale">The scale (σ) of the distribution. Range: σ > 0.</param>
        /// <param name="skew">The skew, 1 > λ > -1</param>
        /// <param name="p">Parameter that controls kurtosis. Range: p > 0</param>
        public SkewedGeneralizedError(double location, double scale, double skew, double p)
        {
            if (!IsValidParameterSet(location, scale, skew, p))
            {
                throw new ArgumentException("Invalid parametrization for the distribution.");
            }

            _random = SystemRandomSource.Default;
            Location = location;
            Scale = scale;
            Skew = skew;
            P = p;

            _skewness = CalculateSkewness();
        }

        /// <summary>
        /// Gets or sets the random number generator which is used to draw random samples.
        /// </summary>
        public System.Random RandomSource
        {
            get => _random;
            set => _random = value ?? SystemRandomSource.Default;
        }

        /// <summary>
        /// A string representation of the distribution.
        /// </summary>
        /// <returns>a string representation of the distribution.</returns>
        public override string ToString()
        {
            return $"SkewedGeneralizedError(μ = {Location}, σ = {Scale}, λ = { Skew }, p = {P}";
        }

        /// <summary>
        /// Tests whether the provided values are valid parameters for this distribution.
        /// </summary>
        /// <param name="location">The location (μ) of the distribution.</param>
        /// <param name="scale">The scale (σ) of the distribution. Range: σ > 0.</param>
        /// <param name="skew">The skew, 1 > λ > -1</param>
        /// <param name="p">Parameter that controls kurtosis. Range: p > 0</param>
        public static bool IsValidParameterSet(double location, double scale, double skew, double p)
        {
            return scale > 0.0 && skew > -1.0 && skew < 1.0 && p > 0.0 && !double.IsNaN(location);
        }

        /// <summary>
        /// Gets the location (μ) of the Skewed Generalized t-distribution.
        /// </summary>
        public double Location { get; private set; }

        /// <summary>
        /// Gets the scale (σ) of the Skewed Generalized t-distribution. Range: σ > 0.
        /// </summary>
        public double Scale { get; private set; }

        /// <summary>
        /// Gets the skew (λ) of the Skewed Generalized t-distribution. Range: 1 > λ > -1.
        /// </summary>
        public double Skew { get; private set; }

        /// <summary>
        /// Gets the parameter that controls the kurtosis of the distribution. Range: p > 0.
        /// </summary>
        public double P { get; private set; }

        // No skew implies Median=Mode=Mean
        public double Mode =>
            Skew == 0 ? Mean : Mean - AdjustAddend(AdjustScale(Scale, Skew, P), Skew, P);

        public double Minimum => double.NegativeInfinity;

        public double Maximum => double.PositiveInfinity;

        // Mean=Location due to our adjustments made
        public double Mean => Location;

        // Variance=Scale*Scale due to our adjustments made
        public double Variance => Scale * Scale;

        public double StdDev => Scale;

        public double Entropy => throw new NotImplementedException();

        public double Skewness => _skewness;

        // No skew implies Median=Mode=Mean
        // Else find it via the point where CDF gives 0.5
        public double Median =>
            Skew == 0 ? Mean : InverseCumulativeDistribution(0.5);

        double CalculateSkewness()
        {
            if (Skew == 0)
            {
                return 0.0;
            }

            var piPow = Math.Pow(Constants.Pi, 3.0 / 2.0);
            var g1 = SpecialFunctions.Gamma(1.0 / P);
            var g2 = SpecialFunctions.Gamma(0.5 + 1.0 / P);
            var g3 = SpecialFunctions.Gamma(3.0 / P);
            var g4 = SpecialFunctions.Gamma(4.0 / P);

            var t1 = Skew * Scale * Scale * Scale / (piPow * g1);
            var t2 = Math.Pow(2.0, (6.0 + P) / P) * Skew * Skew * Math.Pow(g2, 3.0) * g1;
            var t3 = 3.0 * Math.Pow(4.0, 1.0 / P) * Constants.Pi * (1.0 + 3.0 * Skew * Skew) * g2 * g3;
            var t4 = 4.0 * piPow * (1.0 + Skew * Skew) * g4;

            return t1 * (t2 - t3 + t4);
        }

        static double AdjustScale(double scale, double skew, double p)
        {
            var g1 = SpecialFunctions.Gamma(3.0 / p);
            var g2 = SpecialFunctions.Gamma(0.5 + 1.0 / p);
            var g3 = SpecialFunctions.Gamma(1.0 / p);
            var g4 = SpecialFunctions.Gamma(1.0 / p);
            var n1 = Constants.Pi * (1.0 + 3.0 * skew * skew) * g1;
            var n2 = Math.Pow(16.0, 1.0 / p) * skew * skew * Math.Pow(g2, 2) * g3;
            var d = Constants.Pi * g4;
            return scale / Math.Sqrt((n1 - n2) / d);
        }

        static double AdjustX(double x, double scale, double skew, double p)
        {
            return x + AdjustAddend(scale, skew, p);
        }

        static double AdjustAddend(double scale, double skew, double p)
        {
            return (Math.Pow(2.0, 2.0 / p) * scale * skew * SpecialFunctions.Gamma(1.0 / 2.0 + 1.0 / p)) /
                Math.Sqrt(Constants.Pi);
        }

        public static double PDF(double location, double scale, double skew, double p, double x)
        {
            if (!IsValidParameterSet(location, scale, skew, p))
            {
                throw new ArgumentException("Invalid parametrization for the distribution.");
            }

            scale = AdjustScale(scale, skew, p);
            x = AdjustX(x, scale, skew, p);

           // p / (2 * sigma * gamma(1 / p) * exp((abs(x - mu) / (sigma * (1 + lambda * sgn(x - mu)))) ^ p))
            var d1 = Math.Abs(x - location);
            var d2 = scale * (1.0 + skew * Math.Sign(x - location));
            var d3 = 2.0 * scale * SpecialFunctions.Gamma(1.0 / p);
            return p / (Math.Exp(Math.Pow(d1 / d2, p)) * d3);
        }

        public static double PDFLn(double location, double scale, double skew, double p, double x)
        {
            if (!IsValidParameterSet(location, scale, skew, p))
            {
                throw new ArgumentException("Invalid parametrization for the distribution.");
            }

            scale = AdjustScale(scale, skew, p);
            x = AdjustX(x, scale, skew, p);

            return Math.Log(p) - Math.Log(2.0) - Math.Log(scale) - SpecialFunctions.GammaLn(1.0 / p) -
                Math.Pow(Math.Abs(x - location) / (scale * (1.0 + skew * Math.Sign(x - location))), p);
        }

        public static double CDF(double location, double scale, double skew, double p, double x)
        {
            if (!IsValidParameterSet(location, scale, skew, p))
            {
                throw new ArgumentException("Invalid parametrization for the distribution.");
            }

            scale = AdjustScale(scale, skew, p);
            x = AdjustX(x, scale, skew, p) - location;

            var flip = x < 0;
            if (flip)
            {
                skew = -skew;
                x = -x;
            }

            var res = (1.0 - skew) / 2.0 + (1.0 + skew) / 2.0 * Gamma.CDF(1.0 / p, 1.0, Math.Pow(x / (scale * (1.0 + skew)), p));
            return flip ? 1.0 - res : res;
        }

        public static double InvCDF(double location, double scale, double skew, double p, double pr)
        {
            if (!IsValidParameterSet(location, scale, skew, p))
            {
                throw new ArgumentException("Invalid parametrization for the distribution.");
            }

            scale = AdjustScale(scale, skew, p);

            var flip = pr < (1.0 - skew) / 2.0;
            var lambda = skew;
            if (flip)
            {
                pr = 1.0 - pr;
                lambda = -lambda;
            }

            var res = scale * (1.0 + lambda) * Math.Pow(Gamma.InvCDF(1.0 / p, 1.0, 2 * pr / (1.0 + lambda) + (lambda - 1.0) / (lambda + 1.0)), 1.0 / p);

            if (flip)
                res = -res;
            res += location;
            return res - AdjustAddend(scale, skew, p);
        }

        public double InverseCumulativeDistribution(double p)
        {
            return InvCDF(Location, Scale, Skew, P, p);
        }

        public double CumulativeDistribution(double x)
        {
            return CDF(Location, Scale, Skew, P, x);
        }

        public double Density(double x)
        {
            return PDF(Location, Scale, Skew, P, x);
        }

        public double DensityLn(double x)
        {
            return PDFLn(Location, Scale, Skew, P, x);
        }

        public double Sample()
        {
            return SampleUnchecked(_random, Location, Scale, Skew, P);
        }

        public void Samples(double[] values)
        {
            SamplesUnchecked(_random, values, Location, Scale, Skew, P);
        }

        public IEnumerable<double> Samples()
        {
            return SamplesUnchecked(_random, Location, Scale, Skew, P);
        }

        static double SampleUnchecked(System.Random rnd, double location, double scale, double skew, double p)
        {
            var u = ContinuousUniform.Sample(rnd, 0, 1);
            return InvCDF(location, scale, skew, p, u);
        }

        static void SamplesUnchecked(System.Random rnd, double[] values, double location, double scale, double skew, double p)
        {
            for (int i = 0; i < values.Length; i++)
            {
                values[i] = SampleUnchecked(rnd, location, scale, skew, p);
            }
        }

        static IEnumerable<double> SamplesUnchecked(System.Random rnd, double location, double scale, double skew, double p)
        {
            while (true)
            {
                yield return SampleUnchecked(rnd, location, scale, skew, p);
            }
        }

        /// <summary>
        /// Generates a sample from the Skew Generalized Error distribution.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="location">The location (μ) of the distribution.</param>
        /// <param name="scale">The scale (σ) of the distribution. Range: σ > 0.</param>
        /// <param name="skew">The skew, 1 > λ > -1</param>
        /// <param name="p">Parameter that controls kurtosis. Range: p > 0</param>
        /// <returns>a sample from the distribution.</returns>
        public static double Sample(System.Random rnd, double location, double scale, double skew, double p)
        {
            if (!IsValidParameterSet(location, scale, skew, p))
            {
                throw new ArgumentException("Invalid parametrization for the distribution.");
            }

            return SampleUnchecked(rnd, location, scale, skew, p);
        }

        /// <summary>
        /// Generates a sequence of samples from the Skew Generalized Error distribution using inverse transform.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="location">The location (μ) of the distribution.</param>
        /// <param name="scale">The scale (σ) of the distribution. Range: σ > 0.</param>
        /// <param name="skew">The skew, 1 > λ > -1</param>
        /// <param name="p">Parameter that controls kurtosis. Range: p > 0</param>
        /// <returns>a sequence of samples from the distribution.</returns>
        public static IEnumerable<double> Samples(System.Random rnd, double location, double scale, double skew, double p)
        {
            if (!IsValidParameterSet(location, scale, skew, p))
            {
                throw new ArgumentException("Invalid parametrization for the distribution.");
            }

            return SamplesUnchecked(rnd, location, scale, skew, p);
        }

        /// <summary>
        /// Fills an array with samples from the Skew Generalized Error distribution using inverse transform.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="values">The array to fill with the samples.</param>
        /// <param name="location">The location (μ) of the distribution.</param>
        /// <param name="scale">The scale (σ) of the distribution. Range: σ > 0.</param>
        /// <param name="skew">The skew, 1 > λ > -1</param>
        /// <param name="p">Parameter that controls kurtosis. Range: p > 0</param>
        /// <returns>a sequence of samples from the distribution.</returns>
        public static void Samples(System.Random rnd, double[] values, double location, double scale, double skew, double p)
        {
            if (!IsValidParameterSet(location, scale, skew, p))
            {
                throw new ArgumentException("Invalid parametrization for the distribution.");
            }

            SamplesUnchecked(rnd, values, location, scale, skew, p);
        }

        /// <summary>
        /// Generates a sample from the Skew Generalized Error distribution.
        /// </summary>
        /// <param name="location">The location (μ) of the distribution.</param>
        /// <param name="scale">The scale (σ) of the distribution. Range: σ > 0.</param>
        /// <param name="skew">The skew, 1 > λ > -1</param>
        /// <param name="p">Parameter that controls kurtosis. Range: p > 0</param>
        /// <returns>a sample from the distribution.</returns>
        public static double Sample(double location, double scale, double skew, double p)
        {
            if (!IsValidParameterSet(location, scale, skew, p))
            {
                throw new ArgumentException("Invalid parametrization for the distribution.");
            }

            return SampleUnchecked(SystemRandomSource.Default, location, scale, skew, p);
        }

        /// <summary>
        /// Generates a sequence of samples from the Skew Generalized Error distribution using inverse transform.
        /// </summary>
        /// <param name="location">The location (μ) of the distribution.</param>
        /// <param name="scale">The scale (σ) of the distribution. Range: σ > 0.</param>
        /// <param name="skew">The skew, 1 > λ > -1</param>
        /// <param name="p">Parameter that controls kurtosis. Range: p > 0</param>
        /// <returns>a sequence of samples from the distribution.</returns>
        public static IEnumerable<double> Samples(double location, double scale, double skew, double p)
        {
            if (!IsValidParameterSet(location, scale, skew, p))
            {
                throw new ArgumentException("Invalid parametrization for the distribution.");
            }

            return SamplesUnchecked(SystemRandomSource.Default, location, scale, skew, p);
        }

        /// <summary>
        /// Fills an array with samples from the Skew Generalized Error distribution using inverse transform.
        /// </summary>
        /// <param name="values">The array to fill with the samples.</param>
        /// <param name="location">The location (μ) of the distribution.</param>
        /// <param name="scale">The scale (σ) of the distribution. Range: σ > 0.</param>
        /// <param name="skew">The skew, 1 > λ > -1</param>
        /// <param name="p">Parameter that controls kurtosis. Range: p > 0</param>
        /// <returns>a sequence of samples from the distribution.</returns>
        public static void Samples(double[] values, double location, double scale, double skew, double p)
        {
            if (!IsValidParameterSet(location, scale, skew, p))
            {
                throw new ArgumentException("Invalid parametrization for the distribution.");
            }

            SamplesUnchecked(SystemRandomSource.Default, values, location, scale, skew, p);
        }
    }
}