function_rejection_filter.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 <random>
 
#include "point_sampler/point.hpp"
 
namespace ps
{
 
template <typename T, size_t N, typename DensityFn>
std::vector<Point<T, N>> function_rejection_filter(
    const std::vector<Point<T, N>> &points,
    DensityFn                       density_fn,
    std::optional<unsigned int>     seed = std::nullopt)
{
  std::mt19937                      gen(seed ? *seed : std::random_device{}());
  std::uniform_real_distribution<T> accept_dist(0.0, 1.0);
 
  std::vector<Point<T, N>> result;
  result.reserve(points.size());
 
  for (const auto &p : points)
  {
    T prob = density_fn(p);
    T threshold = accept_dist(gen);
 
    if (prob >= threshold)
      result.push_back(p);
  }
 
  return result;
}
 
} // namespace ps