generate_pD.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef cMHN_UTILITY_GENERATE_PD_H
#define cMHN_UTILITY_GENERATE_PD_H
 
#include <cmhn/common/calculate_pTheta.hpp>
#include <cmhn/utility/utility.hpp>
 
#include <prc.hpp>
 
namespace cMHN
{
    template<class T, pRC::Size D>
    static inline auto generatePD(pRC::RandomEngine &rng,
        nonTT::MHNOperator<T, D> const &op, pRC::Size const &size,
        T const &toleranceSolverP)
    {
        using Subscripts =
            decltype(expand(pRC::makeConstantSequence<pRC::Size, D, 2>(),
                [](auto const... seq)
                {
                    return pRC::Subscripts<seq...>();
                }));
 
        std::map<Subscripts, T> pD;
        pRC::UnsignedInteger<64> pDSum = 0;
 
        auto const pTheta = calculatePTheta(op, toleranceSolverP);
 
        // rng setup
        pRC::UniformDistribution<T> dist;
 
        for(pRC::Index i = 0; i < size; ++i)
        {
            T sum = pRC::zero();
            T const r = pRC::random<T>(rng, dist);
            pRC::Index index = 0;
            while(index < Subscripts::size()-1)
            {
                auto pThetaE = pTheta(Subscripts(index));
                sum += pThetaE;
                if(sum >= r)
                {
                    break;
                }
                ++index;
            }
            pD.try_emplace(Subscripts(index), pRC::zero<T>());
            pD[Subscripts(index)] += pRC::unit<T>();
 
            pDSum += pRC::unit<decltype(pDSum)>();
        }
 
        // normalize pD
        for(auto &[k, v] : pD)
        {
            v /= pDSum;
        }
 
        return pD;
    }
 
    template<class T, pRC::Size D>
    static inline auto generatePD(nonTT::MHNOperator<T, D> const &op,
        pRC::Size const &size, T const &toleranceSolverP)
    {
        // rng setup
        pRC::SeedSequence seq(8, 16);
        pRC::RandomEngine rng(seq);
 
        return generatePD(rng, op, size, toleranceSolverP);
    }
 
    template<pRC::Size RP, class T, pRC::Size D>
    static inline auto generatePD(pRC::RandomEngine &rng,
        TT::MHNOperator<T, D> const &op, pRC::Size const &size,
        T const &toleranceSolverP)
    {
        using Subscripts =
            decltype(expand(pRC::makeConstantSequence<pRC::Size, D, 2>(),
                [](auto const... seq)
                {
                    return pRC::Subscripts<seq...>();
                }));
        
        auto const pTheta = calculatePTheta<RP>(op, toleranceSolverP);
 
        // get heuristic for maximum entry of pTheta;
        T maxHeuristic = pRC::zero();
        for(pRC::Index i = 0; i < 1000; ++i)
        {
            auto const subscripts = getRandomSubscripts<Subscripts>();
            auto const value = pTheta(subscripts);
            if(value > maxHeuristic)
            {
                maxHeuristic = value;
            }
        }
 
        // choose prefactor according to heuristic
        T prefactor = T(1) / maxHeuristic;
 
        std::map<Subscripts, T> pD;
        pRC::UnsignedInteger<64> pDSum = 0;
 
        // rng setup
        pRC::UniformDistribution<T> dist;
 
        pRC::Index i = 0;
        while(i < size)
        {
            auto const subscripts = getRandomSubscripts<Subscripts>();
            T const r = pRC::random<T>(rng, dist);
            if(r < prefactor * pTheta(subscripts))
            {
                pD.try_emplace(subscripts, pRC::zero<T>());
                pD[subscripts] += pRC::unit<T>();
                pDSum += pRC::unit<decltype(pDSum)>();
                ++i;
            }
        }
 
        // normalize pD
        for(auto &[k, v] : pD)
        {
            v /= pDSum;
        }
 
        return pD;
    }
 
    template<pRC::Size RP, class T, pRC::Size D>
    static inline auto generatePD(TT::MHNOperator<T, D> const &op,
        pRC::Size const &size, T const &toleranceSolverP)
    {
        // rng setup
        pRC::SeedSequence seq(8, 16);
        pRC::RandomEngine rng(seq);
 
        return generatePD<RP>(rng, op, size, toleranceSolverP);
    }
}
 
#endif // cMHN_UTILITY_GENERATE_PD_H