Sampling/CellGrids/CellGrid.h

0001 // -*- C++ -*-
0002 //
0003 // CellGrid.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_CellGrid_hpp_included
0010 #define EXSAMPLE_CellGrid_hpp_included
0011
0012 #include <vector>
0013 #include <utility>
0014 #include <set>
0015 #include <map>
0016 #include <cstdlib>
0017 #include <cassert>
0018 #include <cmath>
0019
0020 #define units
0021
0022 #ifdef units
0023 #include <iostream>
0024 #include <string>
0025 #endif
0026
0027 #include "Herwig/Utilities/XML/Element.h"
0028
0029 namespace ExSample {
0030
0031   /**
0032    * Simple helper
0033    */
0034   inline double sqr(double x) {
0035     return x*x;
0036   }
0037
0038   /**
0039    * \brief A binary cell grid
0040    * \author Simon Platzer
0041    */
0042   class CellGrid {
0043
0044   public:
0045
0046     /**
0047      * Default constructor
0048      */
0049     CellGrid()
0050       : theVolumeOrIntegral(0.0), theWeight(0.0) {}
0051
0052     /**
0053      * Construct given boundaries and a weight
0054      */
0055     CellGrid(const std::vector<double>& newLowerLeft,
0056          const std::vector<double>& newUpperRight,
0057          double newWeight = 0.0);
0058
0059     /**
0060      * Produce a new instance of a cell grid
0061      */
0062     virtual CellGrid* makeInstance() const;
0063
0064     /**
0065      * Produce a new instance of a cell grid
0066      */
0067     virtual CellGrid* makeInstance(const std::vector<double>& newLowerLeft,
0068                    const std::vector<double>& newUpperRight,
0069                    double newWeight = 0.0) const;
0070
0071     /**
0072      * The destructor
0073      */
0074     virtual ~CellGrid();
0075
0076   public:
0077
0078     /**
0079      * Set the boundaries
0080      */
0081     void boundaries(const std::vector<double>& newLowerLeft,
0082             const std::vector<double>& newUpperRight);
0083
0084     /**
0085      * Return the lower left corner of the cell grid
0086      */
0087     const std::vector<double>& lowerLeft() const { return theLowerLeft; }
0088
0089     /**
0090      * Return the upper right corner of the cell grid
0091      */
0092     const std::vector<double>& upperRight() const { return theUpperRight; }
0093
0094     /**
0095      * Flag in which dimension this cell is upper bound inclusive
0096      */
0097     const std::vector<bool>& upperBoundInclusive() const { return theUpperBoundInclusive; }
0098
0099     /**
0100      * Calculate a volume given upper and lower bound
0101      */
0102     double volume(const std::vector<double>& lowerLeft,
0103           const std::vector<double>& upperRight) const;
0104
0105     /**
0106      * Return true, if this is a leaf in the tree
0107      */
0108     bool isLeaf() const { return theChildren.empty(); }
0109
0110     /**
0111      * Return the depth of this grid
0112      */
0113     std::size_t depth() const;
0114
0115     /**
0116      * Return the number of nodes contained in this grid
0117      */
0118     std::size_t size() const;
0119
0120     /**
0121      * Split this cell grid in the given dimension and coordinate, if
0122      * it is a leaf
0123      */
0124     virtual void split(std::size_t newSplitDimension, double newSplitCoordinate);
0125
0126     /**
0127      * Return the dimension and coordinate along which the first split
0128      * of this grid occurs
0129      */
0130     std::pair<std::size_t,double> splitPoint() const;
0131
0132     /**
0133      * Return the first child
0134      */
0135     const CellGrid& firstChild() const;
0136
0137     /**
0138      * Access the first child
0139      */
0140     CellGrid& firstChild();
0141
0142     /**
0143      * Return the second child
0144      */
0145     const CellGrid& secondChild() const;
0146
0147     /**
0148      * Access the second child
0149      */
0150     CellGrid& secondChild();
0151
0152     /**
0153      * Fill split coordinates along a given dimension
0154      */
0155     void splitCoordinates(std::size_t, std::set<double>&) const;
0156
0157   public:
0158
0159     /**
0160      * Return true, if this grid is active with respect to the last
0161      * parameter point passed.
0162      */
0163     bool active() const;
0164
0165     /**
0166      * Return the volume relevant for the last parameter point set.
0167      */
0168     double volume() const;
0169
0170     /**
0171      * Return the integral relevant for the last parameter point set.
0172      */
0173     double integral() const;
0174
0175     /**
0176      * Return the weight
0177      */
0178     double weight() const { return theWeight; }
0179
0180     /**
0181      * Set the weight
0182      */
0183     void weight(double w);
0184
0185     /**
0186      * Update the integrals
0187      */
0188     void updateIntegral();
0189
0190     /**
0191      * Ensure a minimum cell selection probability
0192      */
0193     void minimumSelection(double p = 0.1);
0194
0195     /**
0196      * Return true, if this grid contains the given parameter point
0197      */
0198     bool contains(const std::vector<double>& point,
0199           const std::vector<bool>& parameterFlags) const;
0200
0201     /**
0202      * Get the volume relevant for the dimensions not considered
0203      * parameters
0204      */
0205     double nonParametricVolume(const std::vector<double>& point,
0206                    const std::vector<bool>& parameterFlags) const;
0207
0208     /**
0209      * Set a parameter point and flag which dimensions are considered
0210      * parameters. Calculate a hash for the relevant subgrid.
0211      */
0212     void updateIntegral(const std::vector<double>& point,
0213             const std::vector<bool>& parameterFlags,
0214             std::vector<bool>::iterator hashPosition);
0215
0216     /**
0217      * Return the projection to the given interval and dimension,
0218      * provided the interval is not overlapping with more than one
0219      * cell in the given dimension. Use splitCoordinates() to
0220      * determine appropriate intervals.
0221      */
0222     double projectInterval(const std::pair<double,double>& interval,
0223                std::size_t dimension) const;
0224
0225     /**
0226      * Return the projection to the given interval and dimension
0227      */
0228     std::map<std::pair<double,double>,double>
0229     project(std::pair<double,double> interval,
0230         std::size_t dimension) const;
0231
0232   public:
0233
0234     /**
0235      * Fill CellGrid data from an XML element
0236      */
0237     virtual void fromXML(const XML::Element&);
0238
0239     /**
0240      * Return an XML element for the data of this CellGrid
0241      */
0242     virtual XML::Element toXML() const;
0243
0244   private:
0245
0246     /**
0247      * Ensure a minimum cell selection probability
0248      */
0249     void doMinimumSelection(double r,
0250                 double ref);
0251
0252     /**
0253      * The lower left corner of the cell grid
0254      */
0255     std::vector<double> theLowerLeft;
0256
0257     /**
0258      * The upper right corner of the cell grid
0259      */
0260     std::vector<double> theUpperRight;
0261
0262     /**
0263      * Flag in which dimension this cell is upper bound inclusive
0264      */
0265     std::vector<bool> theUpperBoundInclusive;
0266
0267     /**
0268      * The volume (leafs) or integral (nodes) relevant for the last
0269      * parameter point set.
0270      */
0271     double theVolumeOrIntegral;
0272
0273     /**
0274      * The weight
0275      */
0276     double theWeight;
0277
0278     /**
0279      * The children cell grids
0280      */
0281     std::vector<CellGrid*> theChildren;
0282
0283 #ifdef units
0284
0285   public:
0286
0287     /**
0288      * Generate a random grid of given maximum depth
0289      */
0290     template<class RndGenerator>
0291     void randomGrid(RndGenerator& rnd,
0292             std::size_t maxDepth,
0293             bool forceSplit = true) {
0294       if ( maxDepth > 0 ) {
0295     if ( !forceSplit )
0296       if ( rnd.rnd() < 0.5 )
0297         return;
0298     std::size_t dimension = (std::size_t)(std::floor(rnd.rnd()*lowerLeft().size()));
0299     double point = 0.5*(lowerLeft()[dimension]+upperRight()[dimension]);
0300     split(dimension,point);
0301     firstChild().weight(rnd.rnd());
0302     secondChild().weight(rnd.rnd());
0303     firstChild().randomGrid(rnd,maxDepth-1,false);
0304     secondChild().randomGrid(rnd,maxDepth-1,false);
0305       }
0306     }
0307
0308     /**
0309      * Write out C code corresponding to the cell grid function
0310      */
0311     void dumpToC(std::ostream& os,
0312          const std::string& name) const;
0313
0314     /**
0315      * Write out C code corresponding to the cell grid function
0316      */
0317     void dumpPartToC(std::ostream& os,
0318              std::string prefix = "") const;
0319
0320 #endif
0321
0322   };
0323
0324 }
0325
0326 #endif // EXSAMPLE_CellGrid_hpp_included
0327