rng.h

Go to the documentation of this file.

#pragma once
 
#include "artstime.h"
#include "debug.h"
#include "matpack_concepts.h"
 
#include <algorithm>
#include <limits>
#include <mutex>
#include <ostream>
#include <random>
#include <vector>
 
template <std::uniform_random_bit_generator Generator = std::mt19937_64>
class RandomNumberGenerator {
  using int_t = typename Generator::result_type;
 
  static constexpr int_t min_v = std::numeric_limits<int_t>::lowest();
  static constexpr int_t max_v = std::numeric_limits<int_t>::max();
 
  // A mutable generator that
  mutable Generator v;
 
  // The
  mutable std::mutex g;
 
  // Keep the seed so we can "copy" the RandomNumberGenerator
  int_t current_seed;
 
public:
  RandomNumberGenerator(const Time &t = {}) { seed(t); }
 
  RandomNumberGenerator(const int_t &s) { seed(s); }
 
  RandomNumberGenerator(const RandomNumberGenerator &rng) {
    std::lock_guard lock{rng.g};
    seed(rng.current_seed);
  }
 
  Generator new_generator() const {
    const std::lock_guard my_lock(this->g);
    return Generator{std::uniform_int_distribution<int_t>(min_v, max_v)(v)};
  }
 
  template <template <typename>
            class random_distribution = std::uniform_real_distribution,
            typename... Ts>
  auto get(Ts &&...x) const {
    return [v = new_generator(),
            draw = random_distribution(std::forward<Ts>(x)...)]() mutable {
      return draw(v);
    };
  }
 
  bool seed(int_t n) {
    static std::mutex seed_mtx;
    static std::vector<int_t> seeds{};
 
    std::lock_guard class_lock{g};
    std::lock_guard method_lock{seed_mtx};
    bool seed_as_is = true;
 
    constexpr auto add_one = [](int_t i) {
      return (i >= max_v - 1) ? min_v : i + 1;
    };
 
    if (std::find(seeds.begin(), seeds.end(), n) != seeds.end()) {
      seed_as_is = false;
 
      auto minmax = std::minmax_element(seeds.begin(), seeds.end());
      int_t min = *minmax.first;
      int_t max = *minmax.second;
 
      if (min not_eq min_v) {
        n = add_one(max);
      } else if (static_cast<std::size_t>(max - min) <= seeds.size()) {
        n = add_one(max);
        if (n == min)
          seeds.clear();
      } else {
        do {
          n = add_one(n);
        } while (std::find(seeds.begin(), seeds.end(), n) != seeds.end());
      }
    }
    seeds.push_back(n);
 
    force_seed(n);
 
    return seed_as_is;
  }
 
  bool seed(const Time &t = {}) {
    return seed(static_cast<int_t>(t.EpochTime().count() % max_v));
  }
 
  void force_seed(int_t n) {
    v.seed(n);
    current_seed = n;
  }
 
  friend std::ostream &operator<<(std::ostream &os,
                                  const RandomNumberGenerator &rng) {
    return os << rng.current_seed;
  }
};
 
template <std::size_t N,
          template <typename>
          class random_distribution = std::uniform_real_distribution,
          std::uniform_random_bit_generator Generator = std::mt19937_64,
          typename... param_types, class RNG = RandomNumberGenerator<Generator>,
          class T = decltype(RNG{}.template get<random_distribution>(
              param_types{}...)())>
matpack::matpack_data<T, N> random_numbers(std::array<Index, N> sz,
                                           param_types &&...x) {
  matpack::matpack_data<T, N> out(sz);
  std::generate(out.elem_begin(), out.elem_end(),
                RNG{}.template get<random_distribution>(
                    std::forward<param_types>(x)...));
  return out;
}
 
template <template <typename>
          class random_distribution = std::uniform_real_distribution,
          std::uniform_random_bit_generator Generator = std::mt19937_64,
          typename... param_types>
auto random_numbers(Index n, param_types &&...x) {
  return random_numbers<1, random_distribution, Generator>(
      std::array<Index, 1>{n}, std::forward<param_types>(x)...);
}