NormalDistribution.h

Go to the documentation of this file.

// Class NormalDistribution
//   This class can be used for sampling normal distribution function
//   on bounded domain, e.g. using Inverse Transform Sampling.
//
#ifndef IPPL_NORMAL_DISTRIBUTION_H
#define IPPL_NORMAL_DISTRIBUTION_H
 
#include "Random/Distribution.h"
#include "Random/Utility.h"
 
namespace ippl {
    namespace random {

        template <typename T>
        KOKKOS_FUNCTION T normal_cdf_func(T x, T mean, T stddev) {
            return 0.5 * (1 + Kokkos::erf((x - mean) / (stddev * Kokkos::sqrt(2.0))));
        }

        template <typename T>
        KOKKOS_FUNCTION T normal_pdf_func(T x, T mean, T stddev) {
            const T pi = Kokkos::numbers::pi_v<T>;
            return (1.0 / (stddev * Kokkos::sqrt(2 * pi)))
                   * Kokkos::exp(-(x - mean) * (x - mean) / (2 * stddev * stddev));
        }

        template <typename T>
        KOKKOS_FUNCTION T normal_estimate_func(T u, T mean, T stddev) {
            return mean + 0. * u * stddev;
        }

        template <typename T>
        struct normal_cdf {
            KOKKOS_INLINE_FUNCTION double operator()(T x, unsigned int d, const T* params_p) const {
                T mean   = params_p[2 * d + 0];
                T stddev = params_p[2 * d + 1];
                return ippl::random::normal_cdf_func<T>(x, mean, stddev);
            }
        };

        template <typename T>
        struct normal_pdf {
            KOKKOS_INLINE_FUNCTION double operator()(T x, unsigned int d, T const* params_p) const {
                T mean   = params_p[2 * d + 0];
                T stddev = params_p[2 * d + 1];
                return ippl::random::normal_pdf_func<T>(x, mean, stddev);
            }
        };

        template <typename T>
        struct normal_estimate {
            KOKKOS_INLINE_FUNCTION double operator()(T u, unsigned int d, T const* params_p) const {
                T mean   = params_p[2 * d + 0];
                T stddev = params_p[2 * d + 1];
                return ippl::random::normal_estimate_func<T>(u, mean, stddev);
            }
        };
 
        template <typename T>
        struct NormalDistributionFunctions {
            // Functor to calculate the probability density function (PDF) for a normal
            // distribution.
            struct PDF {
                KOKKOS_INLINE_FUNCTION double operator()(T x, unsigned int d,
                                                         const T* params_p) const {
                    T mean   = params_p[2 * d + 0];
                    T stddev = params_p[2 * d + 1];
                    return ippl::random::normal_pdf_func<T>(x, mean, stddev);
                }
            };
 
            // Functor to calculate the cumulative distribution function (CDF) for a normal
            // distribution.
            struct CDF {
                KOKKOS_INLINE_FUNCTION double operator()(T x, unsigned int d,
                                                         const T* params_p) const {
                    T mean   = params_p[2 * d + 0];
                    T stddev = params_p[2 * d + 1];
                    return ippl::random::normal_cdf_func<T>(x, mean, stddev);
                }
            };
 
            // Functor to estimate the initial guess for sampling a normal distribution.
            struct Estimate {
                KOKKOS_INLINE_FUNCTION double operator()(T u, unsigned int d,
                                                         T const* params_p) const {
                    T mean   = params_p[2 * d + 0];
                    T stddev = params_p[2 * d + 1];
                    return ippl::random::normal_estimate_func<T>(u, mean, stddev);
                }
            };
        };

        template <typename T, unsigned Dim>
        class NormalDistribution
            : public ippl::random::Distribution<T, Dim, 2 * Dim, NormalDistributionFunctions<T>> {
        public:
            KOKKOS_INLINE_FUNCTION NormalDistribution(const T* par_p)
                : ippl::random::Distribution<T, Dim, 2 * Dim, NormalDistributionFunctions<T>>(
                    par_p) {}
        };
 
    }  // namespace random
}  // namespace ippl
 
#endif