gaussian.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef pRC_CORE_RANDOM_GAUSSIAN_H
#define pRC_CORE_RANDOM_GAUSSIAN_H
 
#include <prc/core/basic/functions/identity.hpp>
#include <prc/core/basic/functions/zero.hpp>
#include <prc/core/log/log.hpp>
#include <prc/core/random/type_traits.hpp>
#include <prc/core/random/uniform.hpp>
#include <prc/core/value/functions/sqrt.hpp>
 
namespace pRC
{
    template<class T>
    class GaussianDistribution<T, If<IsFloat<T>>>
    {
    public:
        constexpr explicit GaussianDistribution(T const mean = zero(),
            T const stdDev = identity())
            : mMean(mean)
            , mStdDev(stdDev)
        {
            if(!(mStdDev > zero()))
            {
                Logging::error(
                    "Parameter standard deviation of Gaussian distribution "
                    "must be positive.");
            }
        }
 
        constexpr auto reset()
        {
            mHot = false;
        }
 
        constexpr auto &mean() const
        {
            return mMean;
        }
 
        constexpr auto &stdDev() const
        {
            return mStdDev;
        }
 
        constexpr auto min() const
        {
            return NumericLimits<T>::lowest();
        }
 
        constexpr auto max() const
        {
            return NumericLimits<T>::max();
        }
 
        template<class URNG>
        constexpr auto operator()(URNG &rng)
        {
            auto basic = [&]() -> T
            {
                if(mHot)
                {
                    mHot = false;
                    return mNextValue;
                }
                else
                {
                    T u, v, s;
                    UniformDistribution<T> uniform(T(-1), T(1));
 
                    do
                    {
                        u = uniform(rng);
                        v = uniform(rng);
                        s = u * u + v * v;
                    }
                    while(s > T(1) || s == zero());
 
                    auto m = sqrt(T(-2) * log(s) / s);
 
                    mNextValue = v * m;
                    mHot = true;
 
                    return u * m;
                }
            }();
            return basic * stdDev() + mean();
        }
 
    private:
        T mMean;
        T mStdDev;
        T mNextValue = zero();
        Bool mHot = false;
    };
 
    template<class T>
    constexpr auto operator==(GaussianDistribution<T> const &lhs,
        GaussianDistribution<T> const &rhs)
    {
        if(lhs.mean() == rhs.mean() && lhs.stdDev() == rhs.stdDev())
        {
            if(lhs.mHot == rhs.mHot)
            {
                if(!lhs.mHot || lhs.mNextValue == rhs.mNextValue)
                {
                    return true;
                }
            }
        }
 
        return false;
    }
 
    template<class T>
    constexpr auto operator!=(GaussianDistribution<T> const &lhs,
        GaussianDistribution<T> const &rhs)
    {
        return !(lhs == rhs);
    }
}
#endif // pRC_CORE_RANDOM_GAUSSIAN_H