VTK  9.3.20240725
vtkReservoirSampler.h
Go to the documentation of this file.
21#ifndef vtkReservoirSampler_h
22#define vtkReservoirSampler_h
23
24#include "vtkAbstractArray.h"
25#include "vtkCommonMathModule.h"
26#include "vtkTypeTraits.h"
27
28#include <algorithm>
29#include <cmath>
30#include <limits>
31#include <random>
32#include <stdexcept>
33
34VTK_ABI_NAMESPACE_BEGIN
35
36class VTKCOMMONMATH_EXPORT vtkReservoirSamplerBase
37{
38protected:
39 using SeedType = typename std::random_device::result_type;
40
42};
43
44template <typename Integer, bool Monotonic = true>
46{
47public:
51 const std::vector<Integer>& operator()(Integer kk, Integer nn) const
52 {
54 this->GenerateSample(kk, nn, data);
55 return data;
56 }
57
61 const std::vector<Integer>& operator()(Integer kk, vtkAbstractArray* array) const
62 {
64 if (!array)
65 {
66 throw std::invalid_argument("Null arrays are disallowed.");
67 }
68 if (array->GetNumberOfTuples() > std::numeric_limits<Integer>::max())
69 {
70 throw std::invalid_argument("Array size would overflow integer type.");
71 }
72 this->GenerateSample(kk, array->GetNumberOfTuples(), data);
73 return data;
74 }
75
76protected:
77 void GenerateSample(Integer kk, Integer nn, std::vector<Integer>& data) const
78 {
79 if (nn < kk)
80 {
81 kk = nn;
82 }
83 if (kk < 0)
84 {
85 throw std::invalid_argument(
86 "You must choose a non-negative number of values from a proper sequence.");
87 }
88 data.resize(kk);
89 if (kk == 0)
90 {
91 return;
92 }
93 // I. Fill the output with the first kk values.
94 Integer ii;
95 for (ii = 0; ii < kk; ++ii)
96 {
97 data[ii] = ii;
98 }
99 if (kk == nn)
100 {
101 return;
102 }
103
104 std::mt19937 generator(vtkReservoirSampler::RandomSeed());
105 std::uniform_real_distribution<> unitUniform(0., 1.);
106 std::uniform_int_distribution<Integer> randomIndex(0, kk - 1);
107 double w = exp(log(unitUniform(generator)) / kk);
108
109 while (true)
110 {
111 double delta = floor(log(unitUniform(generator)) / log(1.0 - w)) + 1.0;
112 if (delta < 0.0 || delta > static_cast<double>(vtkTypeTraits<Integer>::Max()))
113 {
114 // If delta overflows the size of the integer, we are done.
115 break;
116 }
117 Integer intDelta = static_cast<Integer>(delta);
118 // Be careful here since delta may be large and nn may be
119 // at or near numeric_limits<Integer>::max().
120 if (nn - ii > intDelta)
121 {
122 Integer jj = randomIndex(generator);
123#if 0
124 std::cout << " i " << ii << " δ " << intDelta << " w " << w << " → j " << jj << "\n";
125#endif
126 ii += intDelta;
127 data[jj] = ii;
128 w *= exp(log(unitUniform(generator)) / kk);
129 }
130 else
131 {
132 // Adding delta to ii would step beyond our sequence size,
133 // so we are done.
134 break;
135 }
136 }
137 if (Monotonic)
138 {
139 std::sort(data.begin(), data.end());
140 }
141 }
142};
143
144VTK_ABI_NAMESPACE_END
145#endif // vtkReservoirSampler_h