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,
        pRC::Tensor<T, D, D> const &smallThetaGT, pRC::Size const &size)
    {
        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;
 
        // rng setup
        pRC::UniformDistribution<T> dist;
 
        pRC::Bool worked = true;
        pRC::Index at_try = 0;
        do
        {
            ++at_try;
            Subscripts check{};
 
            pD = std::map<Subscripts, T>();
            pDSum = 0;
 
            for(pRC::Index i = 0; i < size; ++i)
            {
                Subscripts sample{};
 
                pRC::Tensor<T, D> transitionRates =
                    exp(extractDiagonal(smallThetaGT));
 
                // simulate until we terminate
                while(true)
                {
                    // simulate one step
                    T const rateSum =
                        reduce<pRC::Add>(transitionRates)() + pRC::unit<T>();
                    T const rand = pRC::random<T>(rng, dist);
                    T sumRejected = pRC::zero();
                    pRC::Index newEvent = 0;
                    while(sumRejected + transitionRates(newEvent) <
                        rand * rateSum)
                    {
                        sumRejected += transitionRates(newEvent);
                        ++newEvent;
                        if(newEvent == D)
                            break;
                    }
                    if(newEvent == D)
                        break;
                    sample[newEvent] = 1;
                    check[newEvent] = 1;
 
                    // update transitionRates
                    transitionRates[newEvent] = pRC::zero();
                    for(pRC::Index j = 0; j < D; ++j)
                    {
                        transitionRates[j] *= exp(smallThetaGT(j, newEvent));
                    }
                }
 
                pD.try_emplace(sample, pRC::zero<T>());
                pD[sample] += pRC::unit<T>();
 
                pDSum += pRC::unit<decltype(pDSum)>();
            }
 
            // check if there is a zero column
            worked = true;
            for(pRC::Index i = 0; i < D; ++i)
            {
                if(check[i] == 0)
                {
                    pRC::Logging::warning(
                        "Zero column detected in generated dataset, "
                        "regenerating!");
                    worked = false;
                }
            }
        }
        while(!worked && at_try < 10);
 
        // normalize pD
        for(auto &[k, v] : pD)
        {
            v /= pDSum;
        }
 
        return pD;
    }
 
    template<class T, pRC::Size D>
    static inline auto generatePD(pRC::Tensor<T, D, D> const &smallThetaGT,
        pRC::Size const &size)
    {
        // rng setup
        pRC::SeedSequence seq(8, 16);
        pRC::RandomEngine rng(seq);
 
        return generatePD(rng, smallThetaGT, size);
    }
}
 
#endif // cMHN_UTILITY_GENERATE_PD_H