Sampling/exsample/generator.icc

0001 // -*- C++ -*-
0002 //
0003 // generator.icc is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
0004 //
0005 // Copyright (C) 2008-2019 Simon Platzer -- simon.plaetzer@desy.de, The Herwig Collaboration
0006 //
0007 // ExSample is licenced under version 3 of the GPL, see COPYING for details.
0008 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0009 //
0010 //
0011 namespace exsample {
0012
0013   template<class Function, class Random>
0014   template<class SlaveStatistics>
0015   void generator<Function,Random>::initialize(SlaveStatistics& opt) {
0016     adaption_info_.dimension = function_->dimension();
0017     adaption_info_.lower_left = function_->support().first;
0018     adaption_info_.upper_right = function_->support().second;
0019     if (adaption_info_.adapt.empty())
0020       adaption_info_.adapt = std::vector<bool>(adaption_info_.dimension,true);
0021     last_point_.resize(adaption_info_.dimension);
0022     if (initialized_) return;
0023     root_cell_ =
0024       binary_tree<cell>(cell(adaption_info_.lower_left,
0025                  adaption_info_.upper_right,
0026                  adaption_info_));
0027     root_cell_.value().info().explore(rnd_gen_,adaption_info_,function_,&statistics_,opt);
0028     root_cell_.value().integral(root_cell_.value().info().volume() * root_cell_.value().info().overestimate());
0029     statistics_.reset();
0030     check_events_ = adaption_info_.presampling_points;
0031     initialized_ = true;
0032   }
0033
0034   template<class Function, class Random>
0035   bool generator<Function,Random>::split() {
0036     if (adaption_info_.freeze_grid < statistics_.accepted() &&
0037     adaption_info_.freeze_grid != 0)
0038       return false;
0039     if (compensating_) return false;
0040     if (!(*last_cell_).info().bad(adaption_info_)) return false;
0041     bool dosplit = false;
0042     std::pair<std::size_t,double> sp =
0043       (*last_cell_).info().get_split(adaption_info_,dosplit);
0044     if (!dosplit || !adaption_info_.adapt[sp.first]) {
0045       return false;
0046     }
0047     last_cell_.node().split((*last_cell_).split(sp,rnd_gen_,function_,adaption_info_));
0048     integral_accessor iacc;
0049     root_cell_.tree_accumulate(iacc,std::plus<double>());
0050     did_split_ = true;
0051     statistics_.reset();
0052     return true;
0053   }
0054
0055   template<class Function, class Random>
0056   void generator<Function,Random>::compensate() {
0057     if (!did_split_) {
0058       root_cell_.value().info().overestimate(std::abs(last_value_),last_point_);
0059       root_cell_.value().integral(root_cell_.value().info().volume() * root_cell_.value().info().overestimate());
0060       return;
0061     }
0062     double old_norm = root_cell_.value().integral();
0063     double new_norm = old_norm - (*last_cell_).integral() + std::abs(last_value_) * (*last_cell_).info().volume();
0064     compensating_ = false;
0065     last_cell_->missing_events() +=
0066       static_cast<int>(round(((std::abs(last_value_) * old_norm)/(last_cell_->info().overestimate() * new_norm) - 1.) *
0067                  (last_cell_->info().attempted())));
0068     if (last_cell_->missing_events() != 0)
0069       compensating_ = true;
0070     last_cell_->info().overestimate(std::abs(last_value_),last_point_);
0071     last_cell_->integral(last_cell_->info().volume() * last_cell_->info().overestimate());
0072     for (binary_tree<cell>::iterator it = root_cell_.begin();
0073      it != root_cell_.end(); ++it)
0074       if (it != last_cell_) {
0075     it->missing_events() += static_cast<int>(round((old_norm/new_norm - 1.) * (it->info().attempted())));
0076     if (it->missing_events() != 0) {
0077       compensating_ = true;
0078     }
0079       }
0080     integral_accessor iacc;
0081     root_cell_.tree_accumulate(iacc,std::plus<double>());
0082     statistics_.reset();
0083     root_cell_.tree_accumulate(missing_accessor(),std::plus<int>());
0084   }
0085
0086   template<class Function, class Random>
0087   template<class SlaveStatistics>
0088   double generator<Function,Random>::generate(SlaveStatistics& opt) {
0089     unsigned long n_hit_miss = 0;
0090     unsigned long n_select = 0;
0091     sampling_selector<rnd_generator<Random> > sampler (rnd_gen_,compensating_);
0092     missing_accessor macc;
0093     if (compensating_) {
0094       compensating_ = false;
0095       for (binary_tree<cell>::iterator it = root_cell_.begin();
0096        it != root_cell_.end(); ++it)
0097     if (it->missing_events() != 0) {
0098       compensating_ = true;
0099       break;
0100     }
0101     }
0102     while (true) {
0103       sampler.compensate = compensating_;
0104       n_select = 0;
0105       if (did_split_)
0106     while ((last_cell_ = root_cell_.select(sampler)) == root_cell_.end()) {
0107       root_cell_.tree_accumulate(macc,std::plus<int>());
0108       if(++n_select > adaption_info_.maxtry)
0109         throw selection_maxtry();
0110     }
0111       else
0112     last_cell_ = root_cell_.begin();
0113       last_cell_->info().select(rnd_gen_,last_point_);
0114       last_value_ = function_->evaluate(last_point_);
0115       last_cell_->info().selected(last_point_,std::abs(last_value_),adaption_info_);
0116       if (std::abs(last_value_) > last_cell_->info().overestimate()) {
0117     if ( std::abs(last_value_)/last_cell_->info().overestimate() > 2. ) {
0118       last_value_ =
0119         last_cell_->info().overestimate()*
0120         (1.+exp(2.*(2.-std::abs(last_value_)/last_cell_->info().overestimate())));
0121     }
0122     compensate();
0123     n_hit_miss = 0;
0124     continue;
0125       }
0126       if (last_cell_->info().attempted() > check_events_) {
0127     if (split()) {
0128       throw generator_update();
0129     }
0130       }
0131       if (did_split_) {
0132     statistics_.select(last_value_ * root_cell_.value().integral() /
0133                last_cell_->info().overestimate(), !compensating_);
0134     opt.select(last_value_ * root_cell_.value().integral() /
0135            last_cell_->info().overestimate(), !compensating_);
0136       } else {
0137     statistics_.select(last_value_, !compensating_);
0138     opt.select(last_value_, !compensating_);
0139       }
0140       if (std::abs(last_value_)/last_cell_->info().overestimate() > rnd_gen_())
0141     break;
0142       if(++n_hit_miss > adaption_info_.maxtry)
0143     throw hit_and_miss_maxtry();
0144     }
0145     last_cell_->info().accept();
0146     if (did_split_)
0147       last_value_ *= root_cell_.value().integral() / last_cell_->info().overestimate();
0148     statistics_.accept(last_value_);
0149     ++check_events_;
0150     if (last_value_ < 0.)
0151       return -1.;
0152     return 1.;
0153   }
0154
0155   template<class Function, class Random>
0156   template<class OStream>
0157   void generator<Function,Random>::put (OStream& os) const {
0158     adaption_info_.put(os);
0159     root_cell_.put(os);
0160     statistics_.put(os);
0161     os << check_events_;
0162     ostream_traits<OStream>::separator(os);
0163     os << did_split_;
0164     ostream_traits<OStream>::separator(os);
0165     os << initialized_;
0166     ostream_traits<OStream>::separator(os);
0167   }
0168
0169   template<class Function, class Random>
0170   template<class IStream>
0171   void generator<Function,Random>::get (IStream& is) {
0172     adaption_info_.get(is);
0173     root_cell_.get(is);
0174     statistics_.get(is);
0175     is >> check_events_ >> did_split_ >> initialized_;
0176   }
0177
0178 }