MCKL
Monte Carlo Kernel Library
collision_test.hpp
Go to the documentation of this file.
1 //============================================================================
2 // MCKL/include/mckl/random/collision_test.hpp
3 //----------------------------------------------------------------------------
4 // MCKL: Monte Carlo Kernel Library
5 //----------------------------------------------------------------------------
6 // Copyright (c) 2013-2018, Yan Zhou
7 // All rights reserved.
8 //
9 // Redistribution and use in source and binary forms, with or without
10 // modification, are permitted provided that the following conditions are met:
11 //
12 // Redistributions of source code must retain the above copyright notice,
13 // this list of conditions and the following disclaimer.
14 //
15 // Redistributions in binary form must reproduce the above copyright notice,
16 // this list of conditions and the following disclaimer in the documentation
17 // and/or other materials provided with the distribution.
18 //
19 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS AS IS
20 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
23 // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 // POSSIBILITY OF SUCH DAMAGE.
30 //============================================================================
31 
32 #ifndef MCKL_RANDOM_COLLISION_TEST_HPP
33 #define MCKL_RANDOM_COLLISION_TEST_HPP
34 
37 #include <bitset>
38 
39 namespace mckl {
40 
52 template <std::size_t D, std::size_t T>
53 class CollisionTest : public PoissonTest<CollisionTest<D, T>>
54 {
55  static_assert(D > 1, "**CollisionTest** used with D less than two");
56 
57  static_assert(T > 0, "**CollisionTest** used with T equal to zero");
58 
59  public:
60  CollisionTest(std::size_t n) : n_(n), mean_(0)
61  {
62  double p = static_cast<double>(n);
63  mean_ = K_ * (p / K_ - 1 + std::pow(1 - 1.0 / K_, p));
64  }
65 
67 
68  template <typename RNGType, typename U01DistributionType>
69  std::size_t operator()(RNGType &rng, U01DistributionType &u01)
70  {
71  using result_type = typename U01DistributionType::result_type;
72 
73  const std::size_t k = internal::BufferSize<result_type, T>::value;
74  const std::size_t m = n_ / k;
75  const std::size_t l = n_ % k;
76  Vector<result_type> r(k * T);
77  std::size_t s = 0;
78  occurs_type occurs;
79  for (std::size_t i = 0; i != m; ++i) {
80  generate(rng, u01, k, r.data(), s, occurs);
81  }
82  generate(rng, u01, l, r.data(), s, occurs);
83 
84  return s;
85  }
86 
87  double mean() const { return mean_; }
88 
89  private:
90  static constexpr std::size_t K_ = internal::Pow<std::size_t, D, T>::value;
91 
92  using occurs_type = std::conditional_t<K_ <= (1U << 27), std::bitset<K_>,
93  std::set<std::size_t>>;
94 
95  std::size_t n_;
96  double mean_;
97 
98  template <typename RNGType, typename U01DistributionType>
99  void generate(RNGType &rng, U01DistributionType &u01, std::size_t n,
100  typename U01DistributionType::result_type *r, std::size_t &s,
101  occurs_type &occurs) const
102  {
103  rand(rng, u01, n * T, r);
104  mul(n * T, static_cast<typename U01DistributionType::result_type>(D),
105  r, r);
106  for (std::size_t i = 0; i != n; ++i, r += T) {
107  count(internal::serial_index<D, T>(r), s, occurs);
108  }
109  }
110 
111  void count(std::size_t u, std::size_t &s, std::bitset<K_> &occurs) const
112  {
113  if (occurs.test(u)) {
114  ++s;
115  } else {
116  occurs.set(u);
117  }
118  }
119 
120  void count(
121  std::size_t u, std::size_t &s, std::set<std::size_t> &occurs) const
122  {
123  if (occurs.count(u)) {
124  ++s;
125  } else {
126  occurs.insert(u);
127  }
128  }
129 }; // class CollisionTest
130 
131 } // namespace mckl
132 
133 #endif // MCKL_RANDOM_COLLISION_TEST_HPP
void pow(std::size_t n, const float *a, const float *b, float *y)
Definition: vmf.hpp:184
void mul(std::size_t n, const float *a, const float *b, float *y)
Definition: vmf.hpp:124
std::size_t operator()(RNGType &rng, U01DistributionType &u01)
std::vector< T, Alloc > Vector
std::vector with Allocator as the default allocator
Definition: memory.hpp:435
CollisionTest(std::size_t n)
RealType u01(UIntType u)
Convert uniform unsigned integers to floating points within [0, 1].
Definition: u01.hpp:88
Collision test.
#define MCKL_DEFINE_RANDOM_TEST_OPERATOR(ResultType)
Definition: mcmc.hpp:40
void rand(RNGType &rng, ArcsineDistribution< RealType > &distribution, std::size_t N, RealType *r)
Tests based on the Poisson distribution.
double mean() const