herwig is hosted by Hepforge, IPPP Durham
Herwig  7.2.1
utility.h
1 // -*- C++ -*-
2 //
3 // utility.h is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
4 //
5 // Copyright (C) 2008-2019 Simon Platzer -- simon.plaetzer@desy.de, The Herwig Collaboration
6 //
7 // ExSample is licenced under version 3 of the GPL, see COPYING for details.
8 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
9 //
10 //
11 #ifndef EXSAMPLE_utility_h_included
12 #define EXSAMPLE_utility_h_included
13 
14 #include "config.h"
15 #include <memory>
16 
17 namespace exsample {
18 
20  template<class OStream>
21  struct ostream_traits {
22 
24  static void separator(OStream& os) { os << "\n"; }
25 
26  };
27 
28 #ifdef EXSAMPLE_has_ThePEG
29 
31  template<>
32  struct ostream_traits<ThePEG::PersistentOStream> {
33 
36 
37  };
38 
39 #endif // EXSAMPLE_has_ThePEG
40 
42  template<unsigned long>
43  struct static_binary {
44  enum { value = 1 };
45  };
46 
48  template<>
49  struct static_binary<0> {
50  enum { value = 0 };
51  };
52 
54  template<unsigned long bits>
55  struct bit_container {
56 
57  enum {
59  n_bits = bits,
61  uint_bits = CHAR_BIT * sizeof(unsigned long),
63  n_segments = bits / uint_bits + static_binary<bits % uint_bits>::value
64  };
65 
66 
69  for (std::size_t i = 0; i < n_segments; ++i)
70  segments[i] = 0;
71  }
72 
74  void reset() {
75  for (std::size_t i = 0; i < n_segments; ++i)
76  segments[i] = 0;
77  }
78 
80  bool operator==(const bit_container& x) const {
81  for (std::size_t i = 0; i < n_segments; ++i)
82  if(segments[i] != x.segments[i])
83  return false;
84  return true;
85  }
86 
88  bool operator<(const bit_container& x) const {
89  for (std::size_t i = 0; i < n_segments; ++i)
90  if(segments[i] != x.segments[i])
91  return (segments[i] < x.segments[i]);
92  return false;
93  }
94 
96  void bit(unsigned long k, bool value) {
97  assert(k<n_bits);
98  if (value)
99  segments[n_segments-k/uint_bits-1] |= (1ul << (k % uint_bits));
100  else
101  segments[n_segments-k/uint_bits-1] &= ~(1ul << (k % uint_bits));
102  }
103 
105  bool bit(unsigned long k) const {
106  assert(k<n_bits);
107  return (segments[n_segments-k/uint_bits-1] & (1ul << (k % uint_bits)));
108  }
109 
111  template<class OStream>
112  void dump(OStream& os) const {
113  for ( unsigned int k = 0; k < n_segments; ++k )
114  os << segments[k] << " ";
115  }
116 
118  template<class OStream>
119  void put(OStream& os) const {
120  for ( size_t k = 0; k < n_segments; ++k ) {
121  os << segments[k];
123  }
124  }
125 
127  template<class IStream>
128  void get(IStream& is) {
129  for ( size_t k = 0; k < n_segments; ++k ) {
130  is >> segments[k];
131  }
132  }
133 
134  private:
135 
137  unsigned long segments[n_segments];
138 
139  };
140 
142  template<class T>
143  T sqr(T x) {
144  return x*x;
145  }
146 
148  template<class T>
149  T cube(T x) {
150  return x*x*x;
151  }
152 
153 
156  template<class T>
157  T round(T x) {
158  T f = std::floor(x);
159  T c = std::ceil(x);
160  if (x < (f+c)/2.)
161  return f;
162  return c;
163  }
164 
166  inline std::size_t two_to(std::size_t n) {
167  assert(n <= sizeof(std::size_t)*CHAR_BIT);
168  return (1 << n);
169  }
170 
173  template<class Statistics>
175 
178  : depth(0), bins(nullptr) {}
179 
182  : depth(x.depth), bins(0) {
183  if (x.bins) {
184  bins.reset(new Statistics[two_to(depth)]);
185  for(std::size_t k = 0; k < two_to(depth); ++k)
186  bins[k] = x.bins[k];
187  }
188  }
189 
192  if (&x == this)
193  return *this;
194  depth = x.depth;
195  bins.reset(nullptr);
196  if (x.bins) {
197  bins.reset(new Statistics[two_to(depth)]);
198  for(std::size_t k = 0; k < two_to(depth); ++k)
199  bins[k] = x.bins[k];
200  }
201  return *this;
202  }
203 
205  explicit fast_small_histogram(std::size_t d)
206  : depth(d), bins(nullptr) {
207  bins.reset(new Statistics[two_to(d)]);
208  }
209 
211  Statistics& bin(double lower,
212  double upper,
213  double event) {
214  double thelower = lower;
215  double theupper = upper;
216  std::size_t bindex = 0;
217  std::size_t current_depth = 0;
218  while (true) {
219  double cut
220  = (thelower+theupper)/2.;
221  if (event < cut) {
222  theupper = cut;
223  } else {
224  thelower = cut;
225  bindex += two_to(depth-current_depth-1);
226  }
227  if(++current_depth == depth)
228  break;
229  }
230  return bins[bindex];
231  }
232 
234  std::size_t depth;
235 
237  std::unique_ptr<Statistics[]> bins;
238 
240  template<class OStream>
241  void put(OStream& os) const {
242  os << depth;
244  for (std::size_t k = 0; k < two_to(depth); ++k) {
245  bins[k].put(os);
246  }
247  }
248 
250  template<class IStream>
251  void get(IStream& is) {
252  is >> depth;
253  bins.reset(new Statistics[two_to(depth)]);
254  for(std::size_t k = 0; k < two_to(depth); ++k) {
255  bins[k].get(is);
256  }
257  }
258 
259  };
260 
263  template<class FirstInputIterator,
264  class SecondInputIterator,
265  class FlagIterator,
266  class OutputIterator,
267  class BinaryOperation>
268  OutputIterator conditional_transform(FirstInputIterator first1,
269  FirstInputIterator last1,
270  SecondInputIterator first2,
271  FlagIterator firstf,
272  OutputIterator result,
273  BinaryOperation binary_op) {
274  for (; first1 != last1; ++first1, ++first2, ++firstf, ++result)
275  if (*firstf)
276  *result = binary_op(*first1, *first2);
277  return result;
278  }
279 
282  inline double volume(const std::vector<double>& lower_left,
283  const std::vector<double>& upper_right) {
284  std::vector<double> delta;
285  std::transform(upper_right.begin(),upper_right.end(),
286  lower_left.begin(),std::back_inserter(delta),
287  std::minus<double>());
288  return
289  std::accumulate(delta.begin(),delta.end(),1.,std::multiplies<double>());
290  }
291 
296  inline double volume(const std::vector<double>& lower_left,
297  const std::vector<double>& upper_right,
298  const std::vector<bool>& flags) {
299  std::vector<double> delta;
300  conditional_transform(upper_right.begin(),upper_right.end(),
301  lower_left.begin(),flags.begin(),
302  std::back_inserter(delta),
303  std::minus<double>());
304  return
305  std::accumulate(delta.begin(),delta.end(),1.,std::multiplies<double>());
306  }
307 
308 
312 
316 
318  template<class Random>
319  struct rnd_generator {
320 
322  double operator()() const {
323  return Random::rnd();
324  }
325 
327  double operator()(double a) const {
328  return a*Random::rnd();
329  }
330 
332  double operator()(double a, double b) const {
333  return (a + (b-a)*Random::rnd());
334  }
335 
336  };
337 
338 }
339 
340 #endif // EXSAMPLE_utility_h_included
Fast, zero memory-overhead one-dimensional histogram with 2^n equally spaced bins.
Definition: utility.h:174
void bit(unsigned long k, bool value)
set the k&#39;th bit
Definition: utility.h:96
Compile time conversion of unsigned long to bool.
Definition: utility.h:43
static void separator(OStream &os)
put the separator to the ostream
Definition: utility.h:24
Fixed-size, packed vector of bools.
Definition: utility.h:55
bool bit(unsigned long k) const
get the k&#39;th bit
Definition: utility.h:105
fast_small_histogram(std::size_t d)
construct from depth d, creating 2^d bins
Definition: utility.h:205
constexpr auto sqr(const T &x) -> decltype(x *x)
Definition: Bin.h:16
Random generator traits.
Definition: utility.h:319
void put(OStream &os) const
put to ostream
Definition: utility.h:119
Exception thrown, if the maximum number of misses has been reached.
Definition: utility.h:315
fast_small_histogram(const fast_small_histogram &x)
copy constructor
Definition: utility.h:181
double operator()(double a, double b) const
Generate uniform random number on [a,b].
Definition: utility.h:332
bool operator<(const bit_container &x) const
compare for ordering
Definition: utility.h:88
std::unique_ptr< Statistics[]> bins
the contained statistics objects
Definition: utility.h:237
Statistics & bin(double lower, double upper, double event)
return the bin from event belongs to given outer boundaries
Definition: utility.h:211
void dump(OStream &os) const
print to ostream
Definition: utility.h:112
fast_small_histogram & operator=(const fast_small_histogram &x)
assignment
Definition: utility.h:191
fast_small_histogram()
default constructor
Definition: utility.h:177
double operator()(double a) const
Generate uniform random number on [0,a].
Definition: utility.h:327
void reset()
put all values to false
Definition: utility.h:74
unsigned long segments[n_segments]
segments needed to keep the hash value
Definition: utility.h:137
std::size_t depth
the depth, defining a histogram of 2^depth bins
Definition: utility.h:234
bool operator==(const bit_container &x) const
compare for equality
Definition: utility.h:80
bit_container()
the default constructor
Definition: utility.h:68
void put(OStream &os) const
put histogram to an ostream
Definition: utility.h:241
separate quantities written to an ostream
Definition: utility.h:21
Exception thrown if the maximum number of attempts to select a cell has been reached.
Definition: utility.h:311
double operator()() const
Generate uniform random number on [0,1].
Definition: utility.h:322
static void separator(ThePEG::PersistentOStream &)
put the separator to the ostream
Definition: utility.h:35