Loading...
Searching...
No Matches
WeightedSampler.h
Go to the documentation of this file.
/* ----------------------------------------------------------------------------
* GTSAM Copyright 2010, Georgia Tech Research Corporation,
* Atlanta, Georgia 30332-0415
* All Rights Reserved
* Authors: Frank Dellaert, et al. (see THANKS for the full author list)
* See LICENSE for the license information
* -------------------------------------------------------------------------- */
#pragma once
#include <cmath>
#include <queue>
#include <random>
#include <stdexcept>
#include <utility>
#include <vector>
namespace gtsam {
/*
* Fast sampling without replacement.
* Example usage:
* std::mt19937 rng(42);
* WeightedSampler<std::mt19937> sampler(&rng);
* auto samples = sampler.sampleWithoutReplacement(5, weights);
*/
template <class Engine = std::mt19937>
class WeightedSampler {
private:
Engine* engine_; // random number generation engine
public:
explicit WeightedSampler(Engine* engine) : engine_(engine) {}
std::vector<size_t> sampleWithoutReplacement(
size_t numSamples, const std::vector<double>& weights) {
// Implementation adapted from code accompanying paper at
// https://www.ethz.ch/content/dam/ethz/special-interest/baug/ivt/ivt-dam/vpl/reports/1101-1200/ab1141.pdf
const size_t n = weights.size();
if (n < numSamples) {
throw std::runtime_error(
"numSamples must be smaller than weights.size()");
}
// Return empty array if numSamples==0
std::vector<size_t> result(numSamples);
if (numSamples == 0) return result;
// Step 1: The first m items of V are inserted into reservoir
// Step 2: For each item v_i ∈ reservoir: Calculate a key k_i = u_i^(1/w),
// where u_i = random(0, 1)
// (Modification: Calculate and store -log k_i = e_i / w where e_i = exp(1),
// reservoir is a priority queue that pops the *maximum* elements)
std::priority_queue<std::pair<double, size_t> > reservoir;
static const double kexp1 = std::exp(1.0);
for (auto it = weights.begin(); it != weights.begin() + numSamples; ++it) {
const double k_i = kexp1 / *it;
reservoir.push(std::make_pair(k_i, it - weights.begin() + 1));
}
// Step 4: Repeat Steps 5–10 until the population is exhausted
{
// Step 3: The threshold T_w is the minimum key of reservoir
// (Modification: This is now the logarithm)
// Step 10: The new threshold T w is the new minimum key of reservoir
const std::pair<double, size_t>& T_w = reservoir.top();
// Incrementing it is part of Step 7
for (auto it = weights.begin() + numSamples; it != weights.end(); ++it) {
// Step 5: Let r = random(0, 1) and X_w = log(r) / log(T_w)
// (Modification: Use e = -exp(1) instead of log(r))
const double X_w = kexp1 / T_w.first;
// Step 6: From the current item v_c skip items until item v_i, such
// that:
double w = 0.0;
// Step 7: w_c + w_{c+1} + ··· + w_{i−1} < X_w <= w_c + w_{c+1} + ··· +
// w_{i−1} + w_i
for (; it != weights.end(); ++it) {
w += *it;
if (X_w <= w) break;
}
// Step 7: No such item, terminate
if (it == weights.end()) break;
// Step 9: Let t_w = T_w^{w_i}, r_2 = random(t_w, 1) and v_i’s key: k_i
// = (r_2)^{1/w_i} (Mod: Let t_w = log(T_w) * {w_i}, e_2 =
// log(random(e^{t_w}, 1)) and v_i’s key: k_i = -e_2 / w_i)
const double t_w = -T_w.first * *it;
std::uniform_real_distribution<double> randomAngle(std::exp(t_w), 1.0);
const double e_2 = std::log(randomAngle(*engine_));
const double k_i = -e_2 / *it;
// Step 8: The item in reservoir with the minimum key is replaced by
// item v_i
reservoir.pop();
reservoir.push(std::make_pair(k_i, it - weights.begin() + 1));
}
}
for (auto iret = result.end(); iret != result.begin();) {
--iret;
if (reservoir.empty()) {
throw std::runtime_error(
"Reservoir empty before all elements have been filled");
}
*iret = reservoir.top().second - 1;
reservoir.pop();
}
if (!reservoir.empty()) {
throw std::runtime_error(
"Reservoir not empty after all elements have been filled");
}
return result;
}
}; // namespace gtsam
} // namespace gtsam
gtsam
Global functions in a separate testing namespace.
Definition chartTesting.h:28
gtsam::WeightedSampler::WeightedSampler
WeightedSampler(Engine *engine)
Construct from random number generation engine We only store a pointer to it.
Definition WeightedSampler.h:46