VTK  9.2.5
vtkReservoirSampler.h
Go to the documentation of this file.
1/*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkReservoirSampler.h
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13
14=========================================================================*/
33#ifndef vtkReservoirSampler_h
34#define vtkReservoirSampler_h
35
36#include "vtkAbstractArray.h"
37#include "vtkCommonMathModule.h"
38
39#include <algorithm>
40#include <cmath>
41#include <limits>
42#include <random>
43#include <stdexcept>
44
45class VTKCOMMONMATH_EXPORT vtkReservoirSamplerBase
46{
47protected:
48 using SeedType = typename std::random_device::result_type;
49
51};
52
53template <typename Integer, bool Monotonic = true>
55{
56public:
60 const std::vector<Integer>& operator()(Integer kk, Integer nn) const
61 {
62 VTK_THREAD_LOCAL static std::vector<Integer> data;
63 this->GenerateSample(kk, nn, data);
64 return data;
65 }
66
70 const std::vector<Integer>& operator()(Integer kk, vtkAbstractArray* array) const
71 {
72 VTK_THREAD_LOCAL static std::vector<Integer> data;
73 if (!array)
74 {
75 throw std::invalid_argument("Null arrays are disallowed.");
76 }
77 if (array->GetNumberOfTuples() > std::numeric_limits<Integer>::max())
78 {
79 throw std::invalid_argument("Array size would overflow integer type.");
80 }
81 this->GenerateSample(kk, array->GetNumberOfTuples(), data);
82 return data;
83 }
84
85protected:
86 void GenerateSample(Integer kk, Integer nn, std::vector<Integer>& data) const
87 {
88 if (nn < kk)
89 {
90 kk = nn;
91 }
92 if (kk < 0)
93 {
94 throw std::invalid_argument(
95 "You must choose a non-negative number of values from a proper sequence.");
96 }
97 data.resize(kk);
98 if (kk == 0)
99 {
100 return;
101 }
102 // I. Fill the output with the first kk values.
103 Integer ii;
104 for (ii = 0; ii < kk; ++ii)
105 {
106 data[ii] = ii;
107 }
108 if (kk == nn)
109 {
110 return;
111 }
112
113 std::mt19937 generator(vtkReservoirSampler::RandomSeed());
114 std::uniform_real_distribution<> unitUniform(0., 1.);
115 std::uniform_int_distribution<Integer> randomIndex(0, kk - 1);
116 double w = exp(log(unitUniform(generator)) / kk);
117
118 while (true)
119 {
120 Integer delta = static_cast<Integer>(floor(log(unitUniform(generator)) / log(1.0 - w)) + 1);
121 if (delta < 0)
122 {
123 // If delta overflows the size of the integer, we are done.
124 break;
125 }
126 // Be careful here since delta may be large and nn may be
127 // at or near numeric_limits<Integer>::max().
128 if (nn - ii > delta)
129 {
130 Integer jj = randomIndex(generator);
131#if 0
132 std::cout << " i " << ii << " δ " << delta << " w " << w << " → j " << jj << "\n";
133#endif
134 ii += delta;
135 data[jj] = ii;
136 w *= exp(log(unitUniform(generator)) / kk);
137 }
138 else
139 {
140 // Adding delta to ii would step beyond our sequence size,
141 // so we are done.
142 break;
143 }
144 }
145 if (Monotonic)
146 {
147 std::sort(data.begin(), data.end());
148 }
149 }
150};
151
152#endif // vtkReservoirSampler_h
153// VTK-HeaderTest-Exclude: vtkReservoirSampler.h
Abstract superclass for all arrays.
vtkIdType GetNumberOfTuples() const
Get the number of complete tuples (a component group) in the array.
static SeedType RandomSeed()
typename std::random_device::result_type SeedType
void GenerateSample(Integer kk, Integer nn, std::vector< Integer > &data) const
const std::vector< Integer > & operator()(Integer kk, Integer nn) const
Choose kk items from a sequence of (0, nn - 1).
const std::vector< Integer > & operator()(Integer kk, vtkAbstractArray *array) const
Choose kk items from a sequence of (0, array->GetNumberOfTuples() - 1).