#ifndef NORMALCDF_H
#define NORMALCDF_H

#include <algorithm>
#include <random>
#include "../../Assertions.h"

namespace Distribution {

        /** cumulative density version of the normal distribution */
        template <typename T> class NormalCDF {

        private:

            const T mu;
            const T sigma;

            static T RationalApproximation(T t)
            {
                // Abramowitz and Stegun formula 26.2.23.
                // The absolute value of the error should be less than 4.5 e-4.
                T c[] = {2.515517, 0.802853, 0.010328};
                T d[] = {1.432788, 0.189269, 0.001308};
                return t - ((c[2]*t + c[1])*t + c[0]) /
                           (((d[2]*t + d[1])*t + d[0])*t + 1.0);
            }

        public:

                /** create a new normally distributed CDF */
                NormalCDF(const T mu, const T sigma) : mu(mu), sigma(sigma) {
                        ;
                }

                /** get the probability for val within the underlying CDF */
                T getProbability(const T val) const {
                        return getProbability(mu, sigma, val);
                }

                /** calculate the probability within the underlying CDF */
                static T getProbability(const T mu, const T sigma, const T val) {
                    return (1.0 + std::exp( (val - mu) / (sigma * std::sqrt(2)) ) ) / 2.0;
                }

                /** get the inverse CDF (https://en.wikipedia.org/wiki/Probit)*/
                static T getProbit(T p){

                    Assert::isBetween(p, 0.0, 1.0, "value not between");

                    // See: https://www.johndcook.com/blog/normal_cdf_inverse/
                    if (p < 0.5)
                    {
                        // F^-1(p) = - G^-1(p)
                        return -RationalApproximation( sqrt(-2.0*log(p)) );
                    }
                    else
                    {
                        // F^-1(p) = G^-1(1-p)
                        return RationalApproximation( sqrt(-2.0*log(1-p)) );
                    }
                }
    };
}

#endif // NORMALCDF_H