VTK  9.2.6
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 
45 class VTKCOMMONMATH_EXPORT vtkReservoirSamplerBase
46 {
47 protected:
48  using SeedType = typename std::random_device::result_type;
49 
50  static SeedType RandomSeed();
51 };
52 
53 template <typename Integer, bool Monotonic = true>
55 {
56 public:
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  }
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 
85 protected:
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.
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).
const std::vector< Integer > & operator()(Integer kk, vtkAbstractArray *array) const
Choose kk items from a sequence of (0, array->GetNumberOfTuples() - 1).
vtkIdType GetNumberOfTuples() const
Get the number of complete tuples (a component group) in the array.
void GenerateSample(Integer kk, Integer nn, std::vector< Integer > &data) const
#define max(a, b)