//******************************************************************************************************
//  GaussianDistribution.cs - Gbtc
//
//  Copyright © 2014, Grid Protection Alliance.  All Rights Reserved.
//
//  Licensed to the Grid Protection Alliance (GPA) under one or more contributor license agreements. See
//  the NOTICE file distributed with this work for additional information regarding copyright ownership.
//  The GPA licenses this file to you under the MIT License (MIT), the "License"; you may
//  not use this file except in compliance with the License. You may obtain a copy of the License at:
//
//      http://www.opensource.org/licenses/MIT
//
//  Unless agreed to in writing, the subject software distributed under the License is distributed on an
//  "AS-IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. Refer to the
//  License for the specific language governing permissions and limitations.
//
//  Code Modification History:
//  ----------------------------------------------------------------------------------------------------
//  04/30/2014 - Stephen E. Chisholm
//       Generated original version of source code.
//
//******************************************************************************************************

//------------------------------------------------------------------------------------------------------
// Code base on Wikipedia article http://en.wikipedia.org/wiki/Box–Muller_transform
//------------------------------------------------------------------------------------------------------

using System;

namespace GSF.NumericalAnalysis
{
    /// <summary>
    /// Implements a BoxMuller method for generating statistically normal random numbers.
    /// </summary>
    public class GaussianDistribution
    {
        #region [ Members ]

        // Fields
        private readonly Random m_random;
        private readonly double m_mean;
        private readonly double m_standardDeviation;
        private readonly double m_min;
        private readonly double m_max;
        private bool m_z1IsValid;
        private double m_z1;

        #endregion

        #region [ Constructors ]

        /// <summary>
        /// Creates a <see cref="GaussianDistribution"/>
        /// </summary>
        /// <param name="mean">the mean of the distribution</param>
        /// <param name="standardDeviation">the standard deviation</param>
        /// <param name="min">a clipping boundary</param>
        /// <param name="max">a clipping boundary</param>
        public GaussianDistribution(double mean, double standardDeviation, double min, double max)
        {
            standardDeviation = Math.Abs(standardDeviation);

            //These limits are set to prevent excessive looping when a value is calculated outside this range.
            if (min > mean - 0.25 * standardDeviation)
                throw new ArgumentOutOfRangeException(nameof(min), "must be less than 1/4 standard deviations away from the mean");

            if (max < mean + 0.25 * standardDeviation)
                throw new ArgumentOutOfRangeException(nameof(max), "must be greater than 1/4 standard deviations away from the mean");

            m_random = new Random(Guid.NewGuid().GetHashCode());
            m_mean = mean;
            m_standardDeviation = standardDeviation;
            m_min = min;
            m_max = max;
        }

        #endregion

        #region [ Methods ]

        /// <summary>
        /// Gets the next random value.
        /// </summary>
        /// <returns></returns>
        public double Next()
        {
            double value;

        TryAgain:

            if (m_z1IsValid)
            {
                m_z1IsValid = false;
                value = m_mean + m_standardDeviation * m_z1;
            }
            else
            {
                double u1 = m_random.NextDouble();
                double u2 = m_random.NextDouble();

                if (u1 < 1e-100)
                    u1 = 1e-100;

                if (u2 < 1e-100)
                    u1 = 1e-100;

                double sqrt = Math.Sqrt(-2 * Math.Log(u1));
                double z0 = sqrt * Math.Sin(2 * Math.PI * u2);

                m_z1 = sqrt * Math.Cos(2 * Math.PI * u2);
                m_z1IsValid = true;

                value = m_mean + m_standardDeviation * z0;
            }

            if (value < m_min || value > m_max)
                goto TryAgain;

            return value;
        }

        #endregion
    }
}