generate_data.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef cMHN_UTILITY_GENERATE_DATA_H
#define cMHN_UTILITY_GENERATE_DATA_H
 
#include <cmhn/common/calculate_pTheta.hpp>
#include <cmhn/utility/utility.hpp>
 
#include <fstream>
#include <iomanip>
#include <string>
 
#include <prc.hpp>
 
namespace cMHN
{
    template<class T, pRC::Size D, class F>
    static inline auto generateData(pRC::RandomEngine<F> &rng,
        pRC::Tensor<T, D, D> const &smallThetaGT, pRC::Size const &size,
        std::string const &header, std::string const &filename)
    {
 
        using Subscripts =
            decltype(expand(pRC::makeConstantSequence<pRC::Size, D, 2>(),
                [](auto const... seq)
                {
                    return pRC::Subscripts<seq...>();
                }));
 
        // rng setup
        pRC::UniformDistribution<T> dist;
 
        pRC::Bool worked = true;
        pRC::Index at_try = 0;
        do
        {
            std::ofstream file(filename,
                std::ofstream::out | std::ofstream::trunc);
 
            if(!file.is_open())
            {
                pRC::Logging::error("Unable to open output file!");
            }
 
            // print header to file
            file << header << std::endl;
 
            ++at_try;
            Subscripts check{};
            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));
                    }
                }
                for(pRC::Index i = 0; i < D; ++i)
                {
                    file << sample[i];
                    if(i != D - 1)
                    {
                        file << " ";
                    }
                }
                file << std::endl;
            }
 
            // 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;
                }
            }
            file.close();
        }
        while(!worked && at_try < 10);
 
        return;
    }
 
    template<class T, pRC::Size D>
    static inline auto generateData(pRC::Tensor<T, D, D> const &smallThetaGT,
        pRC::Size const &size, std::string const &header,
        std::string const &filename)
    {
        // rng setup
        pRC::SeedSequence seq(8, 16);
        pRC::DefaultRandomEngine rng(seq);
 
        return generateData(rng, smallThetaGT, size, header, filename);
    }
}
 
#endif // cMHN_UTILITY_GENERATE_DATA_H