Sampling/exsample/cell.h

0001 // -*- C++ -*-
0002 //
0003 // cell.h 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 #ifndef EXSAMPLE_cell_h_included
0012 #define EXSAMPLE_cell_h_included
0013
0014 #include "utility.h"
0015 #include "adaption_info.h"
0016 #include "statistics.h"
0017
0018 namespace exsample {
0019
0020   /// \brief Information contained in a leaf cell
0021   class cell_info {
0022
0023   public:
0024
0025     /// the default constructor
0026     cell_info();
0027
0028     /// construct from boundaries and adaption info
0029     cell_info(const std::vector<double>& ll,
0030           const std::vector<double>& ur,
0031           const adaption_info& ainfo);
0032
0033     /// construct from boundaries, flags for variables to be sampled,
0034     /// and adaption info
0035     cell_info(const std::vector<double>& ll,
0036           const std::vector<double>& ur,
0037           const std::vector<bool>& sampled_variables,
0038           const adaption_info& ainfo);
0039
0040   public:
0041
0042     /// generate a flat trial point in this cell
0043     template<class Random>
0044     void select(Random&,
0045         std::vector<double>&);
0046
0047     /// generate a flat trial point in this cell
0048     /// only for the variables falgged as true
0049     template<class Random>
0050     void select(Random&,
0051         std::vector<double>&,
0052         const std::vector<bool>&);
0053
0054     /// indicate a function value for the given point
0055     void selected(const std::vector<double>&,
0056           double,
0057           const adaption_info&);
0058
0059     /// indicate that a point has been
0060     /// accepted in this cell
0061     void accept() { ++accepted_; }
0062
0063     /// reject a previously accepted event
0064     void reject() { --accepted_; }
0065
0066   public:
0067
0068     /// return true, if below efficiency threshold
0069     bool bad(const adaption_info& ainfo) const {
0070       return ((static_cast<double>(accepted_)/static_cast<double>(attempted_)) <
0071           ainfo.efficiency_threshold);
0072     }
0073
0074     /// suggest a split and indicate wether it is worth
0075     /// to be performed
0076     std::pair<std::size_t,double> get_split(const adaption_info&,
0077                         bool&) const;
0078
0079     /// explore this cell performing a flat sampling,
0080     /// updating the given statistics object and pre-filling
0081     /// the efficiency histogram by a trial unweighting
0082     template<class Random, class Function, class SlaveStatistics>
0083     void explore(Random&, const adaption_info&, Function*, statistics*,
0084          SlaveStatistics& opt);
0085
0086     /// explore this cell in a more refined way, which
0087     /// is however not suited for already calculating integrals
0088     /// and stuff
0089     template<class Random, class Function>
0090     void explore(Random&, const adaption_info&, Function*);
0091
0092   public:
0093
0094     /// get the current overestimate
0095     double overestimate() const { return overestimate_; }
0096
0097     /// return the position of the last maximum
0098     const std::vector<double>& last_max_position() const { return last_max_position_; }
0099
0100     /// set the current overestimate and maximum position
0101     void overestimate(double v, const std::vector<double>& pos) {
0102       overestimate_ = v;
0103       last_max_position_ = pos;
0104     }
0105
0106     /// get the volume
0107     double volume() const { return volume_; }
0108
0109     /// get the lower left corner
0110     const std::vector<double>& lower_left() const { return lower_left_; }
0111
0112     /// get the upper right corner
0113     const std::vector<double>& upper_right() const { return upper_right_; }
0114
0115     /// get the number of attempted events
0116     unsigned long attempted() const { return attempted_; }
0117
0118     /// get the number of accepted events
0119     unsigned long accepted() const { return accepted_; }
0120
0121   public:
0122
0123     /// return the number of missing events
0124     /// for the given parameter bin id
0125     int parametric_missing(const bit_container<parameter_hash_bits>& id) const;
0126
0127     /// set the number of missing events
0128     /// for the given parameter bin id
0129     void parametric_missing(const bit_container<parameter_hash_bits>& id, int n);
0130
0131     /// increase to the number of missing events
0132     /// for the given parameter bin id
0133     void increase_parametric_missing(const bit_container<parameter_hash_bits>& id);
0134
0135     /// decrease to the number of missing events
0136     /// for the given parameter bin id
0137     void decrease_parametric_missing(const bit_container<parameter_hash_bits>& id);
0138
0139     /// return true, if the cell is compensating in
0140     /// at least one parameter bin
0141     bool parametric_compensating() const {
0142       return !parametric_missing_map_.empty();
0143     }
0144
0145     /// return true, if the cell contains the
0146     /// indicated parameter point
0147     bool contains_parameter(const std::vector<double>& point,
0148                 const std::vector<bool>& sampled) const;
0149
0150   public:
0151
0152     /// put to ostream
0153     template<class OStream>
0154     void put(OStream& os) const;
0155
0156     /// get from istream
0157     template<class IStream>
0158     void get(IStream& is);
0159
0160   private:
0161
0162     /// the value of the overestimate in this cell
0163     double overestimate_;
0164
0165     /// the volume of this cell
0166     double volume_;
0167
0168     /// the lower left corner of this cell
0169     std::vector<double> lower_left_;
0170
0171     /// the upper right corner of this cell
0172     std::vector<double> upper_right_;
0173
0174     /// midpoint of this cell
0175     std::vector<double> mid_point_;
0176
0177     /// the position of the last encountered
0178     /// maximum in this cell
0179     std::vector<double> last_max_position_;
0180
0181     /// left-right statistics of average weight
0182     std::vector<std::pair<double,double> > avg_weight_;
0183
0184     /// the number of attempts in this cell
0185     unsigned long attempted_;
0186
0187     /// the number of accepted events in this cell
0188     unsigned long accepted_;
0189
0190     /// an optional map of parameter bin ids
0191     /// to the number of missing events
0192     std::map<bit_container<parameter_hash_bits>,int> parametric_missing_map_;
0193
0194   };
0195
0196   /// \brief the general cell class
0197   class cell {
0198
0199   public:
0200
0201     /// default constructor
0202     cell();
0203
0204     /// construct from boundaries and adaption info
0205     cell(const std::vector<double>& ll,
0206      const std::vector<double>& ur,
0207      const adaption_info& ainfo);
0208
0209     /// construct from boundaries, flags for variables to be sampled,
0210     /// and adaption info
0211     cell(const std::vector<double>& ll,
0212      const std::vector<double>& ur,
0213      const std::vector<bool>& sampled_variables,
0214      const adaption_info& ainfo);
0215
0216     /// copy constructor
0217     cell(const cell& x);
0218
0219     /// assignment
0220     cell& operator=(const cell& x);
0221
0222   public:
0223
0224     /// split this cell, exploring the
0225     /// child not containing the current overestimate
0226     template<class Random, class Function>
0227     std::pair<cell,cell> split(std::pair<std::size_t,double> split_d,
0228                    Random& rnd_gen,
0229                    Function* f,
0230                    const adaption_info& ainfo,
0231                    const std::vector<bool>& sampled =
0232                    std::vector<bool>());
0233
0234   public:
0235
0236     /// return the split dimension
0237     std::size_t split_dimension() const { return split_dimension_; }
0238
0239     /// return the split value
0240     double split_point() const { return split_point_; }
0241
0242     /// return the integral
0243     double integral() const { return integral_; }
0244
0245     /// access the integral
0246     double& integral() { return integral_; }
0247
0248     /// set the integral
0249     void integral(double v) { integral_ = v; }
0250
0251     /// access the number of missing events
0252     int& missing_events() { return missing_events_; }
0253
0254     /// return the number of missing events
0255     int missing_events() const { return missing_events_; }
0256
0257     /// set the number of missing events
0258     void missing_events(int n) { missing_events_ = n; }
0259
0260     /// access the cell_info object
0261     cell_info& info() { assert(cell_info_); return *cell_info_; }
0262
0263     /// return the cell_info object
0264     const cell_info& info() const { assert(cell_info_); return *cell_info_; }
0265
0266   public:
0267
0268     /// put to ostream
0269     template<class OStream>
0270     void put(OStream& os) const;
0271
0272     /// get from istream
0273     template<class IStream>
0274     void get(IStream& is);
0275
0276   private:
0277
0278     /// the dimension along this cell
0279     /// was split
0280     std::size_t split_dimension_;
0281
0282     /// the value, where this cell was split
0283     double split_point_;
0284
0285     /// the integral of the absolute value
0286     /// of the overestimate over all the
0287     /// children cells
0288     double integral_;
0289
0290     /// the number of missing events in this cell
0291     int missing_events_;
0292
0293     /// a pointer to the cell info object,
0294     /// if this is a leaf cell
0295     std::unique_ptr<cell_info> cell_info_;
0296
0297
0298   };
0299
0300 }
0301
0302 #include "cell.icc"
0303
0304 #endif // EXSAMPLE_cell_h_included