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 */