using System;

using System.Collections.Generic;

using System.Linq;

using System.Text;

using System.Threading.Tasks;

using MathNet.Numerics;

​

namespace TF.Trading.Domain.Services

{

    public class VegasMath

    {

        public static double MuHat(double gl, double fp, double mu)

        {

            return fp + gl * (mu * gl + fp * (1 - gl));

        }

​

        public static double SigHat(double gl, double fp, double muHat)

        {

            return (0.0033 * Math.Pow(20 * gl + 1, 2) - 0.0311 * (20 * gl + 1) + 0.4039) * gl * (muHat - fp);

        }

​

        public static double SigHatConditional(double gl, double fp, double muHat)

        {

            if (gl == 0)

                return 0;

​

            if (muHat >= 0 && 10 >= muHat)

                return (-0.132 * Math.Pow(21 - 20 * gl, 2) - 0.0673 * (21 - 20 * gl) + 8.3644);

            if (muHat > 10 && 20 >= muHat)

                return (-0.0181 * Math.Pow(21 - 20 * gl, 2) + 0.0346 * (21 - 20 * gl) + 8.1582);

            if (muHat > 20 && 30 >= muHat)

                return (-0.0227 * Math.Pow(21 - 20 * gl, 2) + 0.0792 * (21 - 20 * gl) + 9.3414);

            if (muHat > 30 && 40 >= muHat)

                return (-0.0266 * Math.Pow(21 - 20 * gl, 2) + 0.1257 * (21 - 20 * gl) + 10.262);

            if (muHat > 40 && 50 >= muHat)

                return (-0.0293 * Math.Pow(21 - 20 * gl, 2) + 0.1678 * (21 - 20 * gl) + 11.04);

            if (muHat > 50 && 60 >= muHat)

                return (-0.0323 * Math.Pow(21 - 20 * gl, 2) + 0.1867 * (21 - 20 * gl) + 12.364);

            if (muHat > 60 && 70 >= muHat)

                return (-0.0527 * Math.Pow(21 - 20 * gl, 2) + 0.6006 * (21 - 20 * gl) + 12.688);

            if (muHat > 60)

                return -0.6603 * (21 - 20 * gl) + 16.568;

​

            return 0;

        }

​

​

        public static double Gaussian(double x, double mu, double sig)

        {

            if (Math.Abs(sig) < 0.01)

            {

                if (Math.Abs(x - mu) < 0.01)

                    return 1;

                else

                    return 0;

            }

            return Math.Exp(-Math.Pow(x - mu, 2.0) / (2.0 * Math.Pow(sig, 2.0))) /

                (sig * Math.Sqrt(2.0 * Math.PI));

        }

​

        public static double Phi(double x)

        {

            return 0.5 * (1.0 + SpecialFunctions.Erf(x / Math.Sqrt(2.0)));

        }

​

        public static double Integrand(double xi, int i, List<double> mus, List<double> sigs)

        {

            var phiArgs = mus.Zip(sigs, (mu, sig) =>

            (sig <= 0.01 && Math.Abs(xi - mu) <= 0.01) ?

            1 :

            (xi - mu) / sig)

            .ToList();

            phiArgs.RemoveAt(i);

​

            var phiProd = phiArgs

                .Select(phiArg => Phi(phiArg))

                .Aggregate((acc, next) => acc * next);

​

            return Gaussian(xi, mus[i], sigs[i]) * phiProd;

        }

​

        public static double Integral(int i, List<double> mus, List<double> sigs, double xMin, double xMax, int n)

        {

            if (Math.Abs(sigs[i]) < 0.01)

            {

                return Integrand(mus[i], i, mus, sigs);

            }

​

            var interval = (xMax - xMin) / (double)n;

​

            double acc = 0.0;

            for (double xi = xMin; xi <= xMax; xi += interval)

            {

                acc += Integrand(xi, i, mus, sigs);

            }

            return acc * interval;

        }

​

​

        //public static double Integral(int i, List<double> mus, List<double> sigs, double xMin, double xMax, int n)

        //{

        //    if (Math.Abs(sigs[i]) < 0.01)

        //    {

        //        return Integrand(mus[i], i, mus, sigs);

        //    }

​

        //    var interval = (xMax - xMin) / (double)n;

        //    double musCount = mus.Count;

​

        //    var result = Enumerable.Range(0, (int)musCount)

        //        .AsParallel()

        //        .Select(c => xMin + (c / musCount) * (xMax - xMin))

        //        .Select(xi => Integrand(xi, i, mus, sigs) * interval)

        //        .Sum();

        //    return result;

        //}

​

        public static List<double> VegasPrices(List<double> mus, List<double> sigs, double tol = 0.001, double maxIter = 4)

        {

            if (mus.Count() == 1)

            {

                return new List<double>() { 1.0 };

            }

            if (!mus.Any())

            {

                return new List<double>();

            }

​

            var musSigs = mus.Zip(sigs, (m, s) => new {Mu = m, Sig = s}).ToList();

            var xMax = musSigs.Max(ms => ms.Mu + 4 * ms.Sig);

            var xMin = musSigs.Min(ms => ms.Mu - 4 * ms.Sig); 

            int n = 50;

            List<double> vegasPrices = new List<double>();

            int currIter = 0;        

            while ((currIter < maxIter) &&

                ((vegasPrices.Sum() < 1.0 - tol) ||

                (vegasPrices.Sum() > 1.0 + tol)))

            {                

                vegasPrices = Enumerable.Range(0, mus.Count()).AsParallel()

                    .Select(i => new {

                        Index = i,

                        Result = Integral(i, mus, sigs, xMin, xMax, n)})

                    .OrderBy(item => item.Index)

                    .Select(item => item.Result)

                    .ToList();

                n *= 10;

                currIter++;

            }

            return vegasPrices;

        }

    }

}

​