Sampling/CellGrids/SimpleCellGrid.h

0001 // -*- C++ -*-
0002 //
0003 // SimpleCellGrid.hpp is a part of ExSample
0004 // Copyright (C) 2012-2019 Simon Platzer, The Herwig Collaboration
0005 //
0006 // ExSample is licenced under version 3 of the GPL, see COPYING for details.
0007 //
0008
0009 #ifndef EXSAMPLE_SimpleCellGrid_hpp_included
0010 #define EXSAMPLE_SimpleCellGrid_hpp_included
0011
0012 #include "CellGrid.h"
0013 #include <cmath>
0014
0015 namespace ExSample {
0016
0017
0018   /**
0019    * \brief A simple cell grid providing basic adaption and sampling
0020    * \author Simon Platzer
0021    */
0022   class SimpleCellGrid
0023     : public CellGrid {
0024
0025   public:
0026
0027     /**
0028      * Default constructor
0029      */
0030     SimpleCellGrid()
0031       : CellGrid() {}
0032
0033     /**
0034      * Construct given boundaries and a weight
0035      */
0036     SimpleCellGrid(const std::vector<double>& newLowerLeft,
0037            const std::vector<double>& newUpperRight,
0038            bool keepWeightInformation = true,
0039            double newWeight = 0.0);
0040
0041     /**
0042      * Produce a new instance of a cell grid
0043      */
0044     virtual CellGrid* makeInstance() const;
0045
0046     /**
0047      * Produce a new instance of a cell grid
0048      */
0049     virtual CellGrid* makeInstance(const std::vector<double>& newLowerLeft,
0050                    const std::vector<double>& newUpperRight,
0051                    double newWeight = 0.0) const;
0052
0053     /**
0054      * Split this cell grid in the given dimension and coordinate, if
0055      * it is a leaf
0056      */
0057     virtual void split(std::size_t newSplitDimension, double newSplitCoordinate);
0058
0059
0060     virtual void splitter(size_t dim, int rat);
0061
0062   public:
0063
0064     /**
0065      * Return the first child
0066      */
0067     const SimpleCellGrid& firstChild() const {
0068       return dynamic_cast<const SimpleCellGrid&>(CellGrid::firstChild());
0069     }
0070
0071     /**
0072      * Access the first child
0073      */
0074     SimpleCellGrid& firstChild() {
0075       return dynamic_cast<SimpleCellGrid&>(CellGrid::firstChild());
0076     }
0077
0078     /**
0079      * Return the second child
0080      */
0081     const SimpleCellGrid& secondChild() const {
0082       return dynamic_cast<const SimpleCellGrid&>(CellGrid::secondChild());
0083     }
0084
0085     /**
0086      * Access the second child
0087      */
0088     SimpleCellGrid& secondChild() {
0089       return dynamic_cast<SimpleCellGrid&>(CellGrid::secondChild());
0090     }
0091
0092   public:
0093
0094     /**
0095      * A simple counter to store information used for adaption
0096      */
0097     struct Counter {
0098
0099       /**
0100        * Default constructor
0101        */
0102       Counter()
0103     : nPoints(0.0), sumOfWeights(0.0),
0104       sumOfSquaredWeights(0.0),
0105       maxWeight(0.0) {}
0106
0107       /**
0108        * The number of points
0109        */
0110       double nPoints;
0111
0112       /**
0113        * The sum of weights
0114        */
0115       double sumOfWeights;
0116
0117       /**
0118        * The sum of squared weights
0119        */
0120       double sumOfSquaredWeights;
0121
0122       /**
0123        * The maximum weight
0124        */
0125       double maxWeight;
0126
0127       /**
0128        * Book a point
0129        */
0130       void book(double weight) {
0131     nPoints += 1.0;
0132     sumOfWeights += std::abs(weight);
0133     sumOfSquaredWeights += sqr(weight);
0134     maxWeight = std::max(std::abs(weight),maxWeight);
0135       }
0136
0137       /**
0138        * Return the average weight
0139        */
0140       double averageWeight() const { return nPoints != 0.0 ? sumOfWeights/nPoints : 0.0; }
0141
0142       /**
0143        * Return the variance of the weights
0144        */
0145       double varianceOfAverage() const {
0146     return
0147       nPoints > 1.0 ?
0148       fabs(sumOfSquaredWeights/nPoints - sqr(sumOfWeights/nPoints))/(nPoints-1) : 0.0;
0149       }
0150
0151     };
0152
0153     /**
0154      * Return weight information for adaption steps
0155      */
0156     const std::vector<std::pair<Counter,Counter> >& weightInformation() const { return theWeightInformation; }
0157
0158     /**
0159      * Access weight information for adaption steps
0160      */
0161     std::vector<std::pair<Counter,Counter> >& weightInformation() { return theWeightInformation; }
0162
0163     /**
0164      * Update the weight information for the given point
0165      */
0166     virtual void updateWeightInformation(const std::vector<double>& p,
0167                      double w);
0168
0169     /**
0170      * Adjust the reference weight
0171      */
0172     void adjustReferenceWeight(double w) {
0173       theReferenceWeight = std::max(theReferenceWeight,std::abs(w));
0174     }
0175
0176     /**
0177      * Return the reference weight
0178      */
0179     double getReferenceWeight() const {
0180       return theReferenceWeight;
0181     }
0182
0183     /**
0184      * Perform a default adaption step, splitting along the dimension
0185      * which shows up the largest difference in average weights; if
0186      * this exceeds gain, perform the split.
0187      */
0188     virtual void adapt(double gain, double epsilon,
0189                std::set<SimpleCellGrid*>& newCells);
0190
0191     /**
0192      * Update the weights of the cells from information accumulated so
0193      * far
0194      */
0195     virtual void setWeights();
0196
0197   public:
0198
0199     /**
0200      * Sample a point flat in this cell
0201      */
0202     template<class RndGenerator>
0203     void sampleFlatPoint(std::vector<double>& p,
0204              RndGenerator& rnd) const {
0205       assert(p.size() == lowerLeft().size());
0206       for ( size_t k = 0; k < p.size(); ++k ) {
0207     p[k] = lowerLeft()[k] + rnd.rnd()*(upperRight()[k]-lowerLeft()[k]);
0208       }
0209     }
0210
0211     /**
0212      * Sample a point flat in this cell, keeping parameters fixed
0213      */
0214     template<class RndGenerator>
0215     void sampleFlatPoint(std::vector<double>& p,
0216              const std::vector<bool>& parameterFlags,
0217              RndGenerator& rnd) const {
0218       assert(p.size() == lowerLeft().size());
0219       for ( size_t k = 0; k < p.size(); ++k ) {
0220     if ( parameterFlags[k] )
0221       continue;
0222     p[k] = lowerLeft()[k] + rnd.rnd()*(upperRight()[k]-lowerLeft()[k]);
0223       }
0224     }
0225
0226     /**
0227      * Explore the cell grid, given a number of points to be sampled
0228      * in each cell; the weights of the cell will contain the maximum
0229      * weight encountered. If newCells is non-empty explore only these
0230      * cells, otherwise explore all cells.
0231      */
0232     template<class RndGenerator, class Function>
0233     void explore(std::size_t nPoints,
0234          RndGenerator& rnd,
0235          Function& f,
0236          std::set<SimpleCellGrid*>& newCells,
0237          std::ostream& warn) {
0238       unsigned long nanPoints = 0;
0239       if ( !isLeaf() ) {
0240     firstChild().explore(nPoints,rnd,f,newCells,warn);
0241     secondChild().explore(nPoints,rnd,f,newCells,warn);
0242     return;
0243       }
0244       if ( !newCells.empty() ) {
0245     if ( newCells.find(this) == newCells.end() )
0246       return;
0247       }
0248       std::vector<double> point(lowerLeft().size());
0249       for ( std::size_t k = 0; k < nPoints; ++k ) {
0250     sampleFlatPoint(point,rnd);
0251     double w = f.evaluate(point);
0252     if ( ! std::isfinite(w) ) {
0253       ++nanPoints;
0254       continue;
0255     }
0256     updateWeightInformation(point,std::abs(w));
0257       }
0258       if ( nanPoints ) {
0259     warn << "Warning: " << nanPoints << " out of "
0260          << nPoints << " points with nan or inf weight encountered while "
0261          << "exploring a cell.\n" << std::flush;
0262       }
0263     }
0264
0265     /**
0266      * Select a cell
0267      */
0268     template<class RndGenerator>
0269     SimpleCellGrid* selectCell(RndGenerator& rnd) {
0270       if ( isLeaf() )
0271     return this;
0272       if ( firstChild().active() &&
0273        secondChild().active() ) {
0274     double p = firstChild().integral()/integral();
0275     if ( rnd.rnd() <= p )
0276       return firstChild().selectCell(rnd);
0277     else
0278       return secondChild().selectCell(rnd);
0279       }
0280       if ( firstChild().active() &&
0281        !secondChild().active() )
0282     return firstChild().selectCell(rnd);
0283       else
0284     return secondChild().selectCell(rnd);
0285     }
0286
0287     /**
0288      * Sample a point and return its weight
0289      */
0290     template<class RndGenerator, class Function>
0291     double sample(RndGenerator& rnd,
0292           Function& f,
0293           std::vector<double>& p,
0294           bool unweight,
0295           bool adjustReference) {
0296       SimpleCellGrid* selected = selectCell(rnd);
0297       selected->sampleFlatPoint(p,rnd);
0298       double w = f.evaluate(p);
0299       selected->updateWeightInformation(p,w);
0300       double xw = integral()*w/selected->weight();
0301       if ( adjustReference ) {
0302     selected->adjustReferenceWeight(xw);
0303       }
0304       if ( unweight ) {
0305     double r = selected->getReferenceWeight();
0306     if ( r == 0. )
0307       return xw;
0308     double p = std::min(std::abs(xw),r)/r;
0309     double sign = xw >= 0. ? 1. : -1.;
0310     if ( p < 1 && rnd.rnd() > p )
0311       xw = 0.;
0312     else
0313       xw = sign*std::max(std::abs(xw),r);
0314       }
0315       return xw;
0316     }
0317
0318     /**
0319      * Sample a point and return its weight
0320      */
0321     template<class RndGenerator, class Function>
0322     std::pair<double,double> generate(RndGenerator& rnd,
0323                       Function& f,
0324                       std::vector<double>& p) {
0325       SimpleCellGrid* selected = selectCell(rnd);
0326       selected->sampleFlatPoint(p,rnd);
0327       double w = f.evaluate(p);
0328       selected->updateWeightInformation(p,w);
0329       return std::make_pair(w,selected->weight());
0330     }
0331
0332     /**
0333      * Sample a point and return its weight
0334      */
0335     template<class RndGenerator, class Function>
0336     std::pair<double,double> generate(RndGenerator& rnd,
0337                       Function& f,
0338                       std::vector<double>& p,
0339                       const std::vector<bool>& parameterFlags) {
0340       SimpleCellGrid* selected = selectCell(rnd);
0341       selected->sampleFlatPoint(p,parameterFlags,rnd);
0342       double w = f.evaluate(p);
0343       selected->updateWeightInformation(p,w);
0344       return std::make_pair(w,selected->weight());
0345     }
0346
0347   public:
0348
0349     /**
0350      * Fill CellGrid data from an XML element
0351      */
0352     virtual void fromXML(const XML::Element&);
0353
0354     /**
0355      * Return an XML element for the data of this CellGrid
0356      */
0357     virtual XML::Element toXML() const;
0358
0359   private:
0360
0361     /**
0362      * Weight information for adaption steps
0363      */
0364     std::vector<std::pair<Counter,Counter> > theWeightInformation;
0365
0366     /**
0367      * The reference weight to be used for unweighting
0368      */
0369     double theReferenceWeight;
0370
0371   };
0372
0373 }
0374
0375 #endif // EXSAMPLE_SimpleCellGrid_hpp_included
0376