DataGenerator.h

Go to the documentation of this file.

// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
// This software is distributed under the terms of the GNU General Public
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
 
#ifndef DATAGENERATOR_H
#define DATAGENERATOR_H
 
#include <stdexcept> // exeptions, range_error
#include <utility>   // std::forward
#include <random>    // random distribution
#include <cmath>     // exp
//#include <iostream> // lets see if needed, or keep class fre from output
//#include <iomanip>  // lets see if needed, or keep class fre from output
 
namespace o2
{
namespace test
{
 
template <typename ValueT, typename ModelT>
class DataGenerator
{
 public:
  using size_type = std::size_t;
  using value_type = ValueT;
  using result_type = value_type;
  using self_type = DataGenerator;
 
  template <typename... Args>
  DataGenerator(result_type _min, result_type _max, result_type _step, Args&&... args)
    : min(_min), max(_max), step(_step), nbins((max - min) / step), mGenerator(), mModel(std::forward<Args>(args)...)
  {
  }
  ~DataGenerator() = default;
  DataGenerator(const DataGenerator&) = default;
  DataGenerator& operator=(const DataGenerator&) = default;
 
  const result_type min;
  const result_type max;
  const result_type step;
  const size_type nbins;
 
  using random_engine = std::default_random_engine;
 
  // TODO: can it be const?
  value_type operator()()
  {
    value_type v;
    int trials = 0;
    while ((v = mModel(mGenerator)) < min || v >= max) {
      if (trials++ > 1000) {
        // this is a protection, just picked a reasonable threshold for number of trials
        throw std::range_error("random value outside configured range for too many trials");
      }
    }
    int bin = (v - min) / step;
    return min + bin * step + 0.5 * step;
  }
 
  value_type getRandom() const { return (*this)(); }
 
  value_type getMin() const { return ModelT::min; }
 
  value_type getMax() const { return ModelT::max; }
 
  double getProbability(value_type v) const { return mModel.getProbability(v); }
 
  template <class ContainerT>
  class iterator
  {
   public:
    iterator(const ContainerT& parent, size_type count = 0) : mParent(parent), mCount(count) {}
    ~iterator() = default;
 
    using self_type = iterator;
    using reference = typename ContainerT::value_type&;
    using pointer = typename ContainerT::value_type*;
    using difference_type = typename std::iterator_traits<pointer>::difference_type;
    using iterator_category = std::forward_iterator_tag;
 
    // prefix increment
    self_type& operator++()
    {
      if (mCount < mParent.nbins) {
        mCount++;
      }
      return *this;
    }
 
    // postfix increment
    self_type operator++(int /*unused*/)
    {
      self_type copy(*this);
      ++*this;
      return copy;
    }
 
    // addition
    self_type operator+(size_type n) const
    {
      self_type copy(*this);
      if (copy.mCount + n < mParent.nbins) {
        copy.mCount += n;
      } else {
        copy.mCount = mParent.nbins;
      }
      return copy;
    }
 
    value_type operator*() { return mParent.min + mCount * mParent.step + .5 * mParent.step; }
    // pointer operator->() const {return &mValue;}
    // reference operator[](size_type n) const;
 
    bool operator==(const self_type& other) const { return mCount == other.mCount; }
    bool operator!=(const self_type& other) const { return not(*this == other); }
 
   private:
    const ContainerT& mParent;
    size_type mCount;
  };
 
  iterator<self_type> begin() { return iterator<self_type>(*this); }
 
  iterator<self_type> end() { return iterator<self_type>(*this, nbins); }
 
 private:
  random_engine mGenerator;
  ModelT mModel;
};
 
template <class RealType = double, class _BASE = std::normal_distribution<RealType>>
class normal_distribution : public _BASE
{
 public:
  using result_type = typename _BASE::result_type;
 
  normal_distribution(result_type _mean, result_type _stddev) : _BASE(_mean, _stddev), mean(_mean), stddev(_stddev) {}
 
  const double sqrt2pi = 2.5066283;
  const result_type mean;
  const result_type stddev;
 
  // if value_type is an integral type we want to have the probability
  // that the result value is in the range [v, v+1) whereas the step
  // can be something else than 1
  // also the values outside the specified range should be excluded
  // and the probability for intervals in the range has to be scaled
  template <typename value_type>
  double getProbability(value_type v) const
  {
    return (exp(-(v - mean) * (v - mean) / (2 * stddev * stddev))) / (stddev * sqrt2pi);
  }
};
 
template <class IntType = int, class _BASE = std::poisson_distribution<IntType>>
class poisson_distribution : public _BASE
{
 public:
  using result_type = typename _BASE::result_type;
 
  poisson_distribution(result_type _mean) : _BASE(_mean), mean(_mean) {}
  ~poisson_distribution() = default;
 
  const result_type mean;
 
  int factorial(unsigned int n) const { return (n <= 1) ? 1 : factorial(n - 1) * n; }
 
  template <typename value_type>
  double getProbability(value_type v) const
  {
    if (v < 0) {
      return 0.;
    }
    return pow(mean, v) * exp(-mean) / factorial(v);
  }
};
 
template <class IntType = int, class _BASE = std::geometric_distribution<IntType>>
class geometric_distribution : public _BASE
{
 public:
  geometric_distribution(float _parameter) : _BASE(_parameter), parameter(_parameter) {}
 
  const float parameter;
 
  template <typename value_type>
  double getProbability(value_type v) const
  {
    if (v < 0) {
      return 0.;
    }
    return parameter * pow((1 - parameter), v);
  }
};
 
}; // namespace test
}; // namespace o2
#endif