rejection_sampling.hpp

Go to the documentation of this file.

/* Copyright (c) 2025 Otto Link. Distributed under the terms of the GNU General
   Public License. The full license is in the file LICENSE, distributed with
   this software. */
#pragma once
#include <optional>
 
#include "point_sampler/function_rejection_filter.hpp"
 
namespace ps
{
 
template <typename T, size_t N, typename DensityFn>
std::vector<Point<T, N>> rejection_sampling(
    size_t                                count,
    const std::array<std::pair<T, T>, N> &axis_ranges,
    DensityFn                             density_fn,
    std::optional<unsigned int>           seed = std::nullopt)
{
  std::mt19937 gen(seed ? *seed : std::random_device{}());
 
  std::array<std::uniform_real_distribution<T>, N> coord_dists;
  for (size_t i = 0; i < N; ++i)
    coord_dists[i] = std::uniform_real_distribution<T>(axis_ranges[i].first,
                                                       axis_ranges[i].second);
 
  std::vector<Point<T, N>> candidates;
  candidates.reserve(count * 2); // overgenerate
 
  while (candidates.size() < count * 2)
  {
    Point<T, N> p;
    for (size_t i = 0; i < N; ++i)
      p[i] = coord_dists[i](gen);
 
    candidates.push_back(p);
  }
 
  return function_rejection_filter<T, N>(candidates, density_fn, seed);
}
 
} // namespace ps