herwig is hosted by Hepforge, IPPP Durham
Herwig++  2.7.0
BinnedStatistics.h
00001 // -*- C++ -*-
00002 //
00003 // GeneralStatictis.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
00004 // Copyright (C) 2002-2012 The Herwig Collaboration
00005 //
00006 // Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
00007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
00008 //
00009 #ifndef Herwig_BinnedStatistics_H
00010 #define Herwig_BinnedStatistics_H
00011 //
00012 // This is the declaration of the BinnedStatistics class.
00013 //
00014 
00015 #include "GeneralStatistics.h"
00016 #include "ThePEG/Repository/UseRandom.h"
00017 
00018 #include <boost/utility.hpp>
00019 using boost::next;
00020 using boost::prior;
00021 
00022 namespace Herwig {
00023 
00024 using namespace ThePEG;
00025 
00033 class BinnedStatistics {
00034 
00035 public:
00036 
00042   BinnedStatistics() 
00043     : lastPoint(0.), lastStatistics(0), theWeightThreshold(0.001) {}
00044 
00048   BinnedStatistics(unsigned int bins, double threshold = 0.001) 
00049     : lastPoint(0.), lastStatistics(0) {
00050     initialize(bins);
00051     theWeightThreshold = threshold;
00052   }
00053 
00057   virtual ~BinnedStatistics();
00059 
00060 public:
00061 
00065   double sample(double& point) {
00066     const pair<double,double>& range =
00067       selectorMap.upper_bound(UseRandom::rnd())->second;
00068     lastPoint = UseRandom::rnd(range.first,range.second);
00069     point = lastPoint;
00070     lastStatistics = &(statisticsMap.upper_bound(lastPoint)->second);
00071     double weight = weightMap.upper_bound(lastPoint)->second;
00072     return 1./weight;
00073   }
00074 
00078   void bin(double point) {
00079     lastPoint = point;
00080     lastStatistics = &(statisticsMap.upper_bound(lastPoint)->second);
00081   }
00082 
00086   void select(double w) {
00087     lastStatistics->select(w);
00088   }
00089 
00093   void accept() {
00094     lastStatistics->accept();
00095   }
00096 
00100   void reject() {
00101     lastStatistics->reject();
00102   }
00103 
00108   void initialize(unsigned int bins);
00109 
00113   const map<double,GeneralStatistics>& statistics() const {
00114     return statisticsMap;
00115   }
00116 
00121   template<class Adaptor>
00122   void update(const Adaptor& ap) {
00123     double avgweight = 0.;
00124     size_t bins = 0;
00125     for ( map<double,GeneralStatistics>::const_iterator s =
00126             statisticsMap.begin(); s != statisticsMap.end(); ++s ) {
00127       avgweight += ap.importanceMeasure(s->second);
00128       ++bins;
00129     }
00130     avgweight /= bins;
00131     weightMap.clear();
00132     double norm = 0.;
00133     for ( map<double,GeneralStatistics>::const_iterator s =
00134             statisticsMap.begin(); s != statisticsMap.end(); ++s ) {
00135       double weight = ap.importanceMeasure(s->second);
00136       if ( weight < theWeightThreshold*avgweight )
00137         weight = theWeightThreshold*avgweight;
00138       weightMap[s->first] = weight;
00139       norm += 
00140         weight *
00141         (s != statisticsMap.begin() ? (s->first - prior(s)->first) : s->first);
00142     }
00143     selectorMap.clear();
00144     double current = 0.;
00145     for ( map<double,double>::iterator bw = weightMap.begin(); 
00146           bw != weightMap.end(); ++bw ) {
00147       bw->second /= norm;
00148       pair<double,double> range = 
00149         make_pair(bw != weightMap.begin() ? prior(bw)->first : 0.,
00150                   bw->first);
00151       current += bw->second*(range.second-range.first);
00152       selectorMap[current] = range;
00153     }
00154   }
00155 
00160   template<class Adaptor>
00161   void adapt(const Adaptor& ap) {
00162     update(ap);
00163     map<double,GeneralStatistics> newBins;
00164     for ( map<double,GeneralStatistics>::const_iterator b
00165             = statisticsMap.begin(); b != statisticsMap.end(); ++b ) {
00166       newBins[b->first] = GeneralStatistics();
00167       if ( ap.adapt(b->second) ) {
00168         double bound =
00169           b != statisticsMap.begin() ? (prior(b)->first + b->first)/2. : b->first/2.;
00170         newBins[bound] = GeneralStatistics();
00171       }
00172     }
00173     statisticsMap = newBins;
00174   }
00175 
00176 public:
00177 
00184   void put(PersistentOStream & os) const;
00185 
00191   void get(PersistentIStream & is);
00193 
00194 private:
00195 
00200   map<double,GeneralStatistics> statisticsMap;
00201 
00205   map<double,double> weightMap;
00206 
00210   map<double,pair<double,double> > selectorMap;
00211 
00215   double lastPoint;
00216 
00220   GeneralStatistics* lastStatistics;
00221 
00225   double theWeightThreshold;
00226 
00227 };
00228 
00229 inline PersistentOStream& operator<<(PersistentOStream& os, const BinnedStatistics& s) {
00230   s.put(os); return os;
00231 }
00232 
00233 inline PersistentIStream& operator>>(PersistentIStream& is, BinnedStatistics& s) {
00234   s.get(is); return is;
00235 }
00236 
00237 }
00238 
00239 #endif /* Herwig_BinnedStatistics_H */