CompSelector.h

// -*- C++ -*-
//
// CompSelector.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2019 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef THEPEG_CompSelector_H
#define THEPEG_CompSelector_H
//
// This is the declaration of the CompSelector class.
//
 
#include "ThePEG/Utilities/Selector.h"
 
namespace ThePEG {
 
template <typename T, typename WeightType = double>
class CompSelector {
 
public:
 
  CompSelector(double newMargin = 1.1, double newTolerance = 1.0e-6)
    : N(0), last(), theMargin(newMargin), theTolerance(newTolerance) {}
 
public:
 
  WeightType insert(WeightType d, const T & t) {
    reset();
    return selector.insert(d, t);
  }
 
  template <typename RNDGEN>
  T & select(RNDGEN & rnd) {
    ++N;
    if ( !compensating() ) last = selector.select(rnd);
    return last;
  }
 
  WeightType reweight(double & weight) {
    if ( abs(weight) > 1.0 + tolerance() ) {
      // Retrieve the old overestimate of the object by seing how much
      // the summed weights are decreased when removing the object.
      WeightType oldtot = selector.sum();
      WeightType oldmax = oldtot - selector.erase(last);
      WeightType newmax = oldmax*abs(weight)*margin();
      WeightType newtot = selector.insert(newmax, last);
      double rat = newmax/oldmax;
      
      // Setup the new compensation level.
      Level level;
      level.weight = 1.0/rat;
      level.lastN = long(N*newtot/oldtot);
      
      // If we are already compensating, reweight the previous
      // compensation levels.
      for ( int i = 0, M = levels.size(); i < M; ++i ) {
        levels[i].lastN = long(levels[i].lastN*newtot/oldtot);
        levels[i].weight /= rat;
      }
      levels.push_back(level);
      weight /= rat;
      return newmax;
    }
    
    // If we are compensating we should only accept the selection if the
    // weight is above the previous overestimate.
    if ( compensating() ) if ( abs(weight) < levels.back().weight ) weight = 0.0;
    
    return WeightType();
  }
 
  void reset() {
    N = 0;
    levels.clear();
    last = T();
  }
 
  void clear() {
    selector.clear();
    reset();
  }
 
  void margin(double m) { theMargin = m; }
 
  void tolerance(double t) { theTolerance = t; }
 
 
  bool compensating() {
    // Leave all levels which has reached there 'expiry date'.
    while ( levels.size() && levels.back().lastN < N ) levels.pop_back();
    return !levels.empty();
  }
 
  long compleft() const { return levels.empty()? 0: levels.back().lastN - N; }
 
  WeightType sum() const { return selector.sum(); }
 
  double margin() const { return theMargin; }
 
  double tolerance() const { return theTolerance; }
 
  template <typename OStream>
  void output(OStream & os) const {
    os << selector << N << last << theMargin << theTolerance << levels.size();
    for ( int i = 0, M = levels.size(); i < M; ++i )
      os << levels[i].lastN << levels[i].weight;
  }
 
  template <typename IStream>
  void input(IStream & is) {
    long M;
    is >> selector >> N >> last >> theMargin >> theTolerance >> M;
    levels.resize(M);
    for ( int i = 0; i < M; ++i ) is >> levels[i].lastN >> levels[i].weight;
  }
 
private:
 
  struct Level {
 
    long lastN;
 
    double weight;
 
  };
 
private:
 
  Selector<T,WeightType> selector;
 
  long N;
 
  T last;
 
  double theMargin;
 
  double theTolerance;
 
  vector<Level> levels;
 
};
 
template <typename OStream, typename T, typename WeightType>
inline OStream & operator<<(OStream & os,
                            const CompSelector<T,WeightType> & s) {
  s.output(os);
  return os;
}
 
template <typename IStream, typename T, typename WeightType>
inline IStream & operator>>(IStream & is,
                            CompSelector<T,WeightType> & s) {
  s.input(is);
  return is;
}
 
}
 
#endif /* THEPEG_CompSelector_H */