random.h

#pragma once
#include <random>
#include <vector>
#include <cassert>
#include <stdexcept>
#include <Eigen/Core>
#include <nlohmann/json.hpp>

#ifdef ENABLE_PCG
#include <pcg_random.hpp>
#endif

namespace Faunus {

using Point = Eigen::Vector3d;

#ifdef ENABLE_PCG
using RandomNumberEngine = pcg32;
#else
using RandomNumberEngine = std::mt19937;
#endif

class Random
{
  private:
    std::uniform_real_distribution<double> dist01; 
  public:
    RandomNumberEngine engine; 
    Random();                  
    void seed();               
    double operator()();       

    auto range(std::integral auto min, std::integral auto max)
    {
        return std::uniform_int_distribution<>(min, max)(engine);
    }

    template <typename Iterator> Iterator sample(Iterator begin, Iterator end)
    {
        std::advance(begin, range<>(0, std::distance(begin, end) - 1));
        return begin;
    }
};

void to_json(nlohmann::json&, const Random&);   
void from_json(const nlohmann::json&, Random&); 

extern Random random; 

template <typename T> class WeightedDistribution
{
  private:
    std::discrete_distribution<> distribution;
    std::vector<double> weights; 
    size_t latest_index;         
  public:
    std::vector<T> data; 

    auto size() const { return data.size(); } 

    [[nodiscard]] bool empty() const { return data.empty(); } 

    [[nodiscard]] size_t getLastIndex() const
    {
        assert(!data.empty());
        return latest_index;
    } 

    void clear()
    {
        data.clear();
        weights.clear();
    } 

    template <typename Iterator> void setWeight(Iterator begin, Iterator end)
    {
        // Disable assert fow now
        // see
        // https://stackoverflow.com/questions/74496713/how-to-get-the-type-of-the-values-in-a-c20-stdranges-range
        // static_assert(std::is_convertible_v<std::ranges::range_value_t<Iterator>, double>);
        if (auto size = std::distance(begin, end); size == data.size()) {
            weights.resize(size);
            std::copy(begin, end, weights.begin());
            distribution = std::discrete_distribution(weights.begin(), weights.end());
            assert(size_t(distribution.max()) == data.size() - 1);
        }
        else {
            throw std::runtime_error("number of weights must match data");
        }
    }

    void push_back(const T& value, double weight = 1.0)
    {
        data.push_back(value);
        weights.push_back(weight);
        setWeight(weights.begin(), weights.end());
        latest_index = data.size() - 1;
    }

    template <typename RandomGenerator> const T& sample(RandomGenerator& engine)
    {
        assert(not empty());
        latest_index = distribution(engine);
        return data.at(latest_index);
    }
};

Point randomUnitVector(
    Random& rand,
    const Point& directions =
        Point::Ones()); 

Point randomUnitVectorPolar(
    Random& rand); 

} // namespace Faunus