BinnedStatistics.h

// -*- C++ -*-
//
// GeneralStatictis.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef Herwig_BinnedStatistics_H
#define Herwig_BinnedStatistics_H
//
// This is the declaration of the BinnedStatistics class.
//

#include "GeneralStatistics.h"
#include "ThePEG/Repository/UseRandom.h"

#include <boost/utility.hpp>
using boost::next;
using boost::prior;

namespace Herwig {

using namespace ThePEG;

class BinnedStatistics {

public:

  BinnedStatistics() 
    : lastPoint(0.), lastStatistics(0), theWeightThreshold(0.001) {}

  BinnedStatistics(unsigned int bins, double threshold = 0.001) 
    : lastPoint(0.), lastStatistics(0) {
    initialize(bins);
    theWeightThreshold = threshold;
  }

  virtual ~BinnedStatistics();

public:

  double sample(double& point) {
    const pair<double,double>& range =
      selectorMap.upper_bound(UseRandom::rnd())->second;
    lastPoint = UseRandom::rnd(range.first,range.second);
    point = lastPoint;
    lastStatistics = &(statisticsMap.upper_bound(lastPoint)->second);
    double weight = weightMap.upper_bound(lastPoint)->second;
    return 1./weight;
  }

  void bin(double point) {
    lastPoint = point;
    lastStatistics = &(statisticsMap.upper_bound(lastPoint)->second);
  }

  void select(double w) {
    lastStatistics->select(w);
  }

  void accept() {
    lastStatistics->accept();
  }

  void reject() {
    lastStatistics->reject();
  }

  void initialize(unsigned int bins);

  const map<double,GeneralStatistics>& statistics() const {
    return statisticsMap;
  }

  template<class Adaptor>
  void update(const Adaptor& ap) {
    double avgweight = 0.;
    size_t bins = 0;
    for ( map<double,GeneralStatistics>::const_iterator s =
            statisticsMap.begin(); s != statisticsMap.end(); ++s ) {
      avgweight += ap.importanceMeasure(s->second);
      ++bins;
    }
    avgweight /= bins;
    weightMap.clear();
    double norm = 0.;
    for ( map<double,GeneralStatistics>::const_iterator s =
            statisticsMap.begin(); s != statisticsMap.end(); ++s ) {
      double weight = ap.importanceMeasure(s->second);
      if ( weight < theWeightThreshold*avgweight )
        weight = theWeightThreshold*avgweight;
      weightMap[s->first] = weight;
      norm += 
        weight *
        (s != statisticsMap.begin() ? (s->first - prior(s)->first) : s->first);
    }
    selectorMap.clear();
    double current = 0.;
    for ( map<double,double>::iterator bw = weightMap.begin(); 
          bw != weightMap.end(); ++bw ) {
      bw->second /= norm;
      pair<double,double> range = 
        make_pair(bw != weightMap.begin() ? prior(bw)->first : 0.,
                  bw->first);
      current += bw->second*(range.second-range.first);
      selectorMap[current] = range;
    }
  }

  template<class Adaptor>
  void adapt(const Adaptor& ap) {
    update(ap);
    map<double,GeneralStatistics> newBins;
    for ( map<double,GeneralStatistics>::const_iterator b
            = statisticsMap.begin(); b != statisticsMap.end(); ++b ) {
      newBins[b->first] = GeneralStatistics();
      if ( ap.adapt(b->second) ) {
        double bound =
          b != statisticsMap.begin() ? (prior(b)->first + b->first)/2. : b->first/2.;
        newBins[bound] = GeneralStatistics();
      }
    }
    statisticsMap = newBins;
  }

public:

  void put(PersistentOStream & os) const;

  void get(PersistentIStream & is);

private:

  map<double,GeneralStatistics> statisticsMap;

  map<double,double> weightMap;

  map<double,pair<double,double> > selectorMap;

  double lastPoint;

  GeneralStatistics* lastStatistics;

  double theWeightThreshold;

};

inline PersistentOStream& operator<<(PersistentOStream& os, const BinnedStatistics& s) {
  s.put(os); return os;
}

inline PersistentIStream& operator>>(PersistentIStream& is, BinnedStatistics& s) {
  s.get(is); return is;
}

}

#endif /* Herwig_BinnedStatistics_H */