utility.h

// -*- C++ -*-
//
// utility.h is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
//
// Copyright (C) 2008-2019 Simon Platzer -- simon.plaetzer@desy.de, The Herwig Collaboration
//
// ExSample is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
#ifndef EXSAMPLE_utility_h_included
#define EXSAMPLE_utility_h_included
 
#include "config.h"
#include <memory>
 
namespace exsample {
 
  template<class OStream>
  struct ostream_traits {
 
    static void separator(OStream& os) { os << "\n"; }
 
  };
 
#ifdef EXSAMPLE_has_ThePEG
 
  template<>
  struct ostream_traits<ThePEG::PersistentOStream> {
 
    static void separator(ThePEG::PersistentOStream&) { }
 
  };  
 
#endif // EXSAMPLE_has_ThePEG
 
  template<unsigned long>
  struct static_binary {
    enum { value = 1 };
  };
 
  template<>
  struct static_binary<0> {
    enum { value = 0 };
  };
 
  template<unsigned long bits>
  struct bit_container {
 
    enum {
      n_bits = bits,
      uint_bits = CHAR_BIT * sizeof(unsigned long),
      n_segments = bits / uint_bits + static_binary<bits % uint_bits>::value
    };
 
 
    bit_container() {
      for (std::size_t i = 0; i < n_segments; ++i)
        segments[i] = 0;
    }
 
    void reset() {
      for (std::size_t i = 0; i < n_segments; ++i)
        segments[i] = 0;
    }
 
    bool operator==(const bit_container& x) const {
      for (std::size_t i = 0; i < n_segments; ++i)
        if(segments[i] != x.segments[i])
          return false;
      return true;
    }
 
    bool operator<(const bit_container& x) const {
      for (std::size_t i = 0; i < n_segments; ++i)
        if(segments[i] != x.segments[i])
          return (segments[i] < x.segments[i]);
      return false;
    }
 
    void bit(unsigned long k, bool value) {
      assert(k<n_bits);
      if (value)
        segments[n_segments-k/uint_bits-1] |= (1ul << (k % uint_bits));
      else
        segments[n_segments-k/uint_bits-1] &= ~(1ul << (k % uint_bits));
    }
 
    bool bit(unsigned long k) const {
      assert(k<n_bits);
      return (segments[n_segments-k/uint_bits-1] & (1ul << (k % uint_bits)));
    }
 
    template<class OStream>
    void dump(OStream& os) const {
      for ( unsigned int k = 0; k < n_segments; ++k )
        os << segments[k] << " ";
    }
 
    template<class OStream>
    void put(OStream& os) const {
      for ( size_t k = 0; k < n_segments; ++k ) {
        os << segments[k];
        ostream_traits<OStream>::separator(os);
      }
    }
 
    template<class IStream>
    void get(IStream& is) {
      for ( size_t k = 0; k < n_segments; ++k ) {
        is >> segments[k];
      }
    }
 
  private:
 
    unsigned long segments[n_segments];
 
  };
 
  template<class T>
  T sqr(T x) {
    return x*x;
  }
 
  template<class T>
  T cube(T x) {
    return x*x*x;
  }
 
 
  template<class T>
  T round(T x) {
    T f = std::floor(x);
    T c = std::ceil(x);
    if (x < (f+c)/2.)
      return f;
    return c;
  }
 
  inline std::size_t two_to(std::size_t n) {
    assert(n <= sizeof(std::size_t)*CHAR_BIT);
    return (1 << n);
  }
 
  template<class Statistics>
  struct fast_small_histogram {
 
    fast_small_histogram()
      : depth(0), bins(nullptr) {}
 
    fast_small_histogram(const fast_small_histogram& x)
      : depth(x.depth), bins(0) {
      if (x.bins) {
        bins.reset(new Statistics[two_to(depth)]);
        for(std::size_t k = 0; k < two_to(depth); ++k)
          bins[k] = x.bins[k];
      }
    }
 
    fast_small_histogram& operator=(const fast_small_histogram& x) {
      if (&x == this)
        return *this;
      depth = x.depth;
      bins.reset(nullptr);
      if (x.bins) {
        bins.reset(new Statistics[two_to(depth)]);
        for(std::size_t k = 0; k < two_to(depth); ++k)
          bins[k] = x.bins[k];
      }
      return *this;
    }
 
    explicit fast_small_histogram(std::size_t d)
      : depth(d), bins(nullptr) {
      bins.reset(new Statistics[two_to(d)]);
    }
 
    Statistics& bin(double lower,
                    double upper,
                    double event) {
      double thelower = lower;
      double theupper = upper;
      std::size_t bindex = 0;
      std::size_t current_depth = 0;
      while (true) {
        double cut
          = (thelower+theupper)/2.;
        if (event < cut) {
          theupper = cut;
        } else {
          thelower = cut;
          bindex += two_to(depth-current_depth-1);
        }
        if(++current_depth == depth)
          break;
      }
      return bins[bindex];
    }
 
    std::size_t depth;
 
    std::unique_ptr<Statistics[]> bins;
 
    template<class OStream>
    void put(OStream& os) const {
      os << depth;
      ostream_traits<OStream>::separator(os);
      for (std::size_t k = 0; k < two_to(depth); ++k) {
        bins[k].put(os);
      }
    }
 
    template<class IStream>
    void get(IStream& is) {
      is >> depth;
      bins.reset(new Statistics[two_to(depth)]);
      for(std::size_t k = 0; k < two_to(depth); ++k) {
        bins[k].get(is);
      }
    }
     
  };
 
  template<class FirstInputIterator, 
           class SecondInputIterator, 
           class FlagIterator,
           class OutputIterator,
           class BinaryOperation>
  OutputIterator conditional_transform(FirstInputIterator first1,
                                       FirstInputIterator last1,
                                       SecondInputIterator first2,
                                       FlagIterator firstf,
                                       OutputIterator result,
                                       BinaryOperation binary_op) {
    for (; first1 != last1; ++first1, ++first2, ++firstf, ++result)
      if (*firstf)
        *result = binary_op(*first1, *first2);
    return result;
  }
 
  inline double volume(const std::vector<double>& lower_left,
                       const std::vector<double>& upper_right) {
    std::vector<double> delta;
    std::transform(upper_right.begin(),upper_right.end(),
                   lower_left.begin(),std::back_inserter(delta),
                   std::minus<double>());
    return
      std::accumulate(delta.begin(),delta.end(),1.,std::multiplies<double>());
  }
 
  inline double volume(const std::vector<double>& lower_left,
                       const std::vector<double>& upper_right,
                       const std::vector<bool>& flags) {
    std::vector<double> delta;
    conditional_transform(upper_right.begin(),upper_right.end(),
                          lower_left.begin(),flags.begin(),
                          std::back_inserter(delta),
                          std::minus<double>());
    return
      std::accumulate(delta.begin(),delta.end(),1.,std::multiplies<double>());
  }
 
 
  struct selection_maxtry{};
 
  struct hit_and_miss_maxtry{};
 
  template<class Random>
  struct rnd_generator {
 
    double operator()() const {
      return Random::rnd();
    }
 
    double operator()(double a) const {
      return a*Random::rnd();
    }
 
    double operator()(double a, double b) const {
      return (a + (b-a)*Random::rnd());
    }
 
  };
 
}
 
#endif // EXSAMPLE_utility_h_included