VTK  9.3.20240304
vtkReservoirSampler.h
Go to the documentation of this file.
1 // SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
2 // SPDX-License-Identifier: BSD-3-Clause
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 
34 VTK_ABI_NAMESPACE_BEGIN
35 
36 class VTKCOMMONMATH_EXPORT vtkReservoirSamplerBase
37 {
38 protected:
39  using SeedType = typename std::random_device::result_type;
40 
41  static SeedType RandomSeed();
42 };
43 
44 template <typename Integer, bool Monotonic = true>
46 {
47 public:
51  const std::vector<Integer>& operator()(Integer kk, Integer nn) const
52  {
53  VTK_THREAD_LOCAL static std::vector<Integer> data;
54  this->GenerateSample(kk, nn, data);
55  return data;
56  }
57 
61  const std::vector<Integer>& operator()(Integer kk, vtkAbstractArray* array) const
62  {
63  VTK_THREAD_LOCAL static std::vector<Integer> data;
64  if (!array)
65  {
66  throw std::invalid_argument("Null arrays are disallowed.");
67  }
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 
76 protected:
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 
144 VTK_ABI_NAMESPACE_END
145 #endif // vtkReservoirSampler_h
146 // 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
const std::vector< Integer > & operator()(Integer kk, Integer nn) const
Choose kk items from a sequence of (0, nn - 1).
void GenerateSample(Integer kk, Integer nn, std::vector< Integer > &data) const
const std::vector< Integer > & operator()(Integer kk, vtkAbstractArray *array) const
Choose kk items from a sequence of (0, array->GetNumberOfTuples() - 1).
@ data
Definition: vtkX3D.h:315
Template defining traits of native types used by VTK.
Definition: vtkTypeTraits.h:23
#define max(a, b)