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

namespace Mathematics
{
    public class HoltWinters
    {
        /// <summary>
        /// returns best initial trend factor so that we can start forecasting.
        /// </summary>
        public static double InitialTrendFactor(double[] data, int seasonLength)
        {
            if (data.Length < seasonLength * 2)
            {
                throw new Exception("Not enough data to compute initial trend factor, need at least two seasons");
            }

            double tot = 0;
            for (int i = 0; i < seasonLength; i++)
            {
                tot += (data[i + seasonLength] - data[i]) / (double)seasonLength;

            }
            return tot / (double)seasonLength;
        }

        /// <summary>
        /// returns a list of best initial seasonal factor so that we can start forecasting.
        /// </summary>
        public static double[] InitialSeasonalFactors(double[] data, int seasonLength)
        {

            int nSeasons = data.Length / seasonLength;
            if (nSeasons == 0)
            {
                throw new Exception("Not enough data to compute initial seasonal factors, at least one full season required");
            }

            // take the same offset from each of the seasons and average it
            double[] seasonAverages = new double[nSeasons];

            for (int j = 0; j < nSeasons; j++)
            {
                double sum = 0;
                for (int i = 0; i < seasonLength; i++)
                {
                    sum += data[j * seasonLength + i];

                }
                seasonAverages[j] = sum / (double)seasonLength;
            }

            double[] result = new double[seasonLength];
            for (int i = 0; i < seasonLength; i++)
            {
                double sumOverAvg = 0;
                for (int j = 0; j < nSeasons; j++)
                {
                    sumOverAvg += data[j * seasonLength + i] / seasonAverages[j];
                }
                result[i] = sumOverAvg / (double)nSeasons;
            }

            return result;
        }

        /// <summary>
        /// performs triple exponential smoothing (multiplicative) on data
        /// You can obtain trend and seasonal from InitialTrendFactor and InitialSeasonalFactors.
        /// Best alpha, beta and gamma can be calculated by repeatedly calling
        /// TripleExponentialSmoothing to find the combination with the smallest SSE,
        /// you can use MinimiseSSE for this.
        /// </summary>
        /// <param name="seasonLength">number of data points per season</param>
        /// <param name="trend">initial trend</param>
        /// <param name="seasonal">slice of initial seasonal</param>
        /// <param name="alpha">smoothing factor</param>
        /// <param name="beta">smoothing factor</param>
        /// <param name="gamma">smoothing factor</param>
        /// <returns>
        /// A new slice of smoothed and forcasted data points, a slice of deviation
        /// such that y±y*d could be used as "confidence bands" as well as SSE
        /// (Sum of Squared Errors)
        /// </returns>
        public static double TripleExponentialSmoothing(double[] data, int seasonLength, double trend, double[] seasonal
            , int nPredictions, double alpha, double beta, double gamma
            , out double[] result, out double[] deviation)
        {
            double sse = 0;
            double level = 0, lastLevel = 0;
            double[] seasonalDev = new double[seasonLength];

            result = new double[data.Length + nPredictions];
            deviation = new double[data.Length + nPredictions];

            if ((alpha <= 0 || alpha >= 1) || (beta <= 0 || beta >= 1) || (gamma <= 0 || gamma >= 1))
            {
                return double.MaxValue;
            }

            for (int i = 0; i < data.Length + nPredictions; i++)
            {

                if (i == 0)
                {
                    level = data[i];
                    result[i] = data[i];
                    continue;
                }

                if (i >= data.Length)
                {
                    // we are forecasting
                    int m = i - data.Length + 1;

                    result[i] = (level + (double)(m) * trend) * seasonal[i % seasonLength];

                }
                else
                {
                    double val = data[i];

                    lastLevel = level;
                    level = alpha * val / seasonal[i % seasonLength] + (1 - alpha) * (level + trend);
                    trend = beta * (level - lastLevel) + (1 - beta) * trend;

                    seasonal[i % seasonLength] = gamma * val / level + (1 - gamma) * seasonal[i % seasonLength];

                    result[i] = (level + trend) * seasonal[i % seasonLength];

                    if (i > 2)
                    {
                        // smoothed deviation for confidence band
                        double factor = result[i] / val;
                        if (factor > 1)
                        {
                            factor = 1 / factor;
                        }
                        seasonalDev[i % seasonLength] = gamma * (1 - factor) + (1 - gamma) * seasonalDev[i % seasonLength];
                        // SSE
                        double diff = result[i] - val;
                        sse += diff * diff;
                    }
                }
                deviation[i] = seasonalDev[i % seasonLength];
            }

            return sse;
        }

        /// <summary>
        // HWMinimizeSSE repeatedly executes TripleExponentialSmoothing over
        // data, season, seasonLength, trend and nPred with varying alpha, beta and gamma using
        // the BOBYQA algorithm to minimize SSE.  It returns the final
        // (best) smooth and deviation returned by TripleExponentialSmoothing, as
        // well as the resulting alpha, beta, gamma (which can be reused later), k
        // (number of passes) and e (evocation count for
        // TripleExponentialSmoothing).
        /// </summary>
        public static double[] MinimizeSSE(double[] data, int seasonLength, int nPredictions,
            int maxIterations, out double alpha, out double beta, out double gamma)
        {
            double trend = InitialTrendFactor(data, seasonLength);
            double[] seasonal = InitialSeasonalFactors(data, seasonLength);

            BobyqaFunctionDelegate callFn = delegate (int l, double[] x)
            {
                double[] smooth, deviation;

                double[] sc = new double[seasonal.Length];
                Array.Copy(seasonal, sc, sc.Length);

                double sse = TripleExponentialSmoothing(data, seasonLength, trend, sc, nPredictions
                    , x[0], x[1], x[2], out smooth, out deviation);

                return sse;
            };

            double[] startPoints = new double[] { 0.5, 0.2, 0.2 };
            int n = startPoints.Length;

            var optimizer = new BOBYQA(n, callFn
                , new double[] { 0, 0, 0 }, new double[] { 1, 1, 1 })
            {
                TrustRegionRadiusEnd = 1.0e-8,
                MaximumFunctionCalls = maxIterations,
                Logger = null // Console.Out
            };

            optimizer.InterpolationConditions = 2 * n + 1;
            var r = optimizer.FindMinimum(startPoints);

            alpha = r.X[0];
            beta = r.X[1];
            gamma = r.X[2];

            double[] result, dev;

            TripleExponentialSmoothing(data, seasonLength, trend, seasonal, nPredictions
                    , alpha, beta, gamma, out result, out dev);

            var pred = new double[nPredictions];
            Array.Copy(result, result.Length - nPredictions, pred, 0, nPredictions);

            return pred;
        }
    }
}