herwig is hosted by Hepforge, IPPP Durham
Herwig  7.2.1
SimpleCellGrid.h
1 // -*- C++ -*-
2 //
3 // SimpleCellGrid.hpp is a part of ExSample
4 // Copyright (C) 2012-2019 Simon Platzer, The Herwig Collaboration
5 //
6 // ExSample is licenced under version 3 of the GPL, see COPYING for details.
7 //
8 
9 #ifndef EXSAMPLE_SimpleCellGrid_hpp_included
10 #define EXSAMPLE_SimpleCellGrid_hpp_included
11 
12 #include "CellGrid.h"
13 #include <cmath>
14 
15 namespace ExSample {
16 
17 
23  : public CellGrid {
24 
25  public:
26 
31  : CellGrid() {}
32 
36  SimpleCellGrid(const std::vector<double>& newLowerLeft,
37  const std::vector<double>& newUpperRight,
38  bool keepWeightInformation = true,
39  double newWeight = 0.0);
40 
44  virtual CellGrid* makeInstance() const;
45 
49  virtual CellGrid* makeInstance(const std::vector<double>& newLowerLeft,
50  const std::vector<double>& newUpperRight,
51  double newWeight = 0.0) const;
52 
57  virtual void split(std::size_t newSplitDimension, double newSplitCoordinate);
58 
59 
60  virtual void splitter(size_t dim, int rat);
61 
62  public:
63 
67  const SimpleCellGrid& firstChild() const {
68  return dynamic_cast<const SimpleCellGrid&>(CellGrid::firstChild());
69  }
70 
75  return dynamic_cast<SimpleCellGrid&>(CellGrid::firstChild());
76  }
77 
81  const SimpleCellGrid& secondChild() const {
82  return dynamic_cast<const SimpleCellGrid&>(CellGrid::secondChild());
83  }
84 
89  return dynamic_cast<SimpleCellGrid&>(CellGrid::secondChild());
90  }
91 
92  public:
93 
97  struct Counter {
98 
103  : nPoints(0.0), sumOfWeights(0.0),
104  sumOfSquaredWeights(0.0),
105  maxWeight(0.0) {}
106 
110  double nPoints;
111 
115  double sumOfWeights;
116 
121 
125  double maxWeight;
126 
130  void book(double weight) {
131  nPoints += 1.0;
132  sumOfWeights += std::abs(weight);
133  sumOfSquaredWeights += sqr(weight);
134  maxWeight = std::max(std::abs(weight),maxWeight);
135  }
136 
140  double averageWeight() const { return nPoints != 0.0 ? sumOfWeights/nPoints : 0.0; }
141 
145  double varianceOfAverage() const {
146  return
147  nPoints > 1.0 ?
148  fabs(sumOfSquaredWeights/nPoints - sqr(sumOfWeights/nPoints))/(nPoints-1) : 0.0;
149  }
150 
151  };
152 
156  const std::vector<std::pair<Counter,Counter> >& weightInformation() const { return theWeightInformation; }
157 
161  std::vector<std::pair<Counter,Counter> >& weightInformation() { return theWeightInformation; }
162 
166  virtual void updateWeightInformation(const std::vector<double>& p,
167  double w);
168 
172  void adjustReferenceWeight(double w) {
173  theReferenceWeight = std::max(theReferenceWeight,std::abs(w));
174  }
175 
179  double getReferenceWeight() const {
180  return theReferenceWeight;
181  }
182 
188  virtual void adapt(double gain, double epsilon,
189  std::set<SimpleCellGrid*>& newCells);
190 
195  virtual void setWeights();
196 
197  public:
198 
202  template<class RndGenerator>
203  void sampleFlatPoint(std::vector<double>& p,
204  RndGenerator& rnd) const {
205  assert(p.size() == lowerLeft().size());
206  for ( size_t k = 0; k < p.size(); ++k ) {
207  p[k] = lowerLeft()[k] + rnd.rnd()*(upperRight()[k]-lowerLeft()[k]);
208  }
209  }
210 
214  template<class RndGenerator>
215  void sampleFlatPoint(std::vector<double>& p,
216  const std::vector<bool>& parameterFlags,
217  RndGenerator& rnd) const {
218  assert(p.size() == lowerLeft().size());
219  for ( size_t k = 0; k < p.size(); ++k ) {
220  if ( parameterFlags[k] )
221  continue;
222  p[k] = lowerLeft()[k] + rnd.rnd()*(upperRight()[k]-lowerLeft()[k]);
223  }
224  }
225 
232  template<class RndGenerator, class Function>
233  void explore(std::size_t nPoints,
234  RndGenerator& rnd,
235  Function& f,
236  std::set<SimpleCellGrid*>& newCells,
237  std::ostream& warn) {
238  unsigned long nanPoints = 0;
239  if ( !isLeaf() ) {
240  firstChild().explore(nPoints,rnd,f,newCells,warn);
241  secondChild().explore(nPoints,rnd,f,newCells,warn);
242  return;
243  }
244  if ( !newCells.empty() ) {
245  if ( newCells.find(this) == newCells.end() )
246  return;
247  }
248  std::vector<double> point(lowerLeft().size());
249  for ( std::size_t k = 0; k < nPoints; ++k ) {
250  sampleFlatPoint(point,rnd);
251  double w = f.evaluate(point);
252  if ( ! std::isfinite(w) ) {
253  ++nanPoints;
254  continue;
255  }
256  updateWeightInformation(point,std::abs(w));
257  }
258  if ( nanPoints ) {
259  warn << "Warning: " << nanPoints << " out of "
260  << nPoints << " points with nan or inf weight encountered while "
261  << "exploring a cell.\n" << std::flush;
262  }
263  }
264 
268  template<class RndGenerator>
269  SimpleCellGrid* selectCell(RndGenerator& rnd) {
270  if ( isLeaf() )
271  return this;
272  if ( firstChild().active() &&
273  secondChild().active() ) {
274  double p = firstChild().integral()/integral();
275  if ( rnd.rnd() <= p )
276  return firstChild().selectCell(rnd);
277  else
278  return secondChild().selectCell(rnd);
279  }
280  if ( firstChild().active() &&
281  !secondChild().active() )
282  return firstChild().selectCell(rnd);
283  else
284  return secondChild().selectCell(rnd);
285  }
286 
290  template<class RndGenerator, class Function>
291  double sample(RndGenerator& rnd,
292  Function& f,
293  std::vector<double>& p,
294  bool unweight,
295  bool adjustReference) {
296  SimpleCellGrid* selected = selectCell(rnd);
297  selected->sampleFlatPoint(p,rnd);
298  double w = f.evaluate(p);
299  selected->updateWeightInformation(p,w);
300  double xw = integral()*w/selected->weight();
301  if ( adjustReference ) {
302  selected->adjustReferenceWeight(xw);
303  }
304  if ( unweight ) {
305  double r = selected->getReferenceWeight();
306  if ( r == 0. )
307  return xw;
308  double p = std::min(std::abs(xw),r)/r;
309  double sign = xw >= 0. ? 1. : -1.;
310  if ( p < 1 && rnd.rnd() > p )
311  xw = 0.;
312  else
313  xw = sign*std::max(std::abs(xw),r);
314  }
315  return xw;
316  }
317 
321  template<class RndGenerator, class Function>
322  std::pair<double,double> generate(RndGenerator& rnd,
323  Function& f,
324  std::vector<double>& p) {
325  SimpleCellGrid* selected = selectCell(rnd);
326  selected->sampleFlatPoint(p,rnd);
327  double w = f.evaluate(p);
328  selected->updateWeightInformation(p,w);
329  return std::make_pair(w,selected->weight());
330  }
331 
335  template<class RndGenerator, class Function>
336  std::pair<double,double> generate(RndGenerator& rnd,
337  Function& f,
338  std::vector<double>& p,
339  const std::vector<bool>& parameterFlags) {
340  SimpleCellGrid* selected = selectCell(rnd);
341  selected->sampleFlatPoint(p,parameterFlags,rnd);
342  double w = f.evaluate(p);
343  selected->updateWeightInformation(p,w);
344  return std::make_pair(w,selected->weight());
345  }
346 
347  public:
348 
352  virtual void fromXML(const XML::Element&);
353 
357  virtual XML::Element toXML() const;
358 
359  private:
360 
364  std::vector<std::pair<Counter,Counter> > theWeightInformation;
365 
370 
371  };
372 
373 }
374 
375 #endif // EXSAMPLE_SimpleCellGrid_hpp_included
376 
double sample(RndGenerator &rnd, Function &f, std::vector< double > &p, bool unweight, bool adjustReference)
Sample a point and return its weight.
const SimpleCellGrid & firstChild() const
Return the first child.
std::vector< std::pair< Counter, Counter > > & weightInformation()
Access weight information for adaption steps.
virtual XML::Element toXML() const
Return an XML element for the data of this CellGrid.
const CellGrid & firstChild() const
Return the first child.
A simple counter to store information used for adaption.
std::pair< double, double > generate(RndGenerator &rnd, Function &f, std::vector< double > &p, const std::vector< bool > &parameterFlags)
Sample a point and return its weight.
double averageWeight() const
Return the average weight.
A binary cell grid.
Definition: CellGrid.h:42
void adjustReferenceWeight(double w)
Adjust the reference weight.
A simple cell grid providing basic adaption and sampling.
virtual void adapt(double gain, double epsilon, std::set< SimpleCellGrid *> &newCells)
Perform a default adaption step, splitting along the dimension which shows up the largest difference ...
double weight() const
Return the weight.
Definition: CellGrid.h:178
constexpr auto sqr(const T &x) -> decltype(x *x)
const CellGrid & secondChild() const
Return the second child.
virtual void split(std::size_t newSplitDimension, double newSplitCoordinate)
Split this cell grid in the given dimension and coordinate, if it is a leaf.
SimpleCellGrid()
Default constructor.
virtual void fromXML(const XML::Element &)
Fill CellGrid data from an XML element.
double nPoints
The number of points.
void sampleFlatPoint(std::vector< double > &p, RndGenerator &rnd) const
Sample a point flat in this cell.
virtual void setWeights()
Update the weights of the cells from information accumulated so far.
std::pair< double, double > generate(RndGenerator &rnd, Function &f, std::vector< double > &p)
Sample a point and return its weight.
double integral() const
Return the integral relevant for the last parameter point set.
virtual CellGrid * makeInstance() const
Produce a new instance of a cell grid.
double sumOfSquaredWeights
The sum of squared weights.
void sampleFlatPoint(std::vector< double > &p, const std::vector< bool > &parameterFlags, RndGenerator &rnd) const
Sample a point flat in this cell, keeping parameters fixed.
const SimpleCellGrid & secondChild() const
Return the second child.
std::size_t size() const
Return the number of nodes contained in this grid.
double sumOfWeights
The sum of weights.
SimpleCellGrid & firstChild()
Access the first child.
Element represents a (tree of) XML elements.
Definition: Element.h:56
double theReferenceWeight
The reference weight to be used for unweighting.
double getReferenceWeight() const
Return the reference weight.
Counter()
Default constructor.
virtual void updateWeightInformation(const std::vector< double > &p, double w)
Update the weight information for the given point.
const std::vector< double > & lowerLeft() const
Return the lower left corner of the cell grid.
Definition: CellGrid.h:87
double sign(double x)
Small helper.
Definition: RandomHelpers.h:30
const std::vector< std::pair< Counter, Counter > > & weightInformation() const
Return weight information for adaption steps.
void explore(std::size_t nPoints, RndGenerator &rnd, Function &f, std::set< SimpleCellGrid *> &newCells, std::ostream &warn)
Explore the cell grid, given a number of points to be sampled in each cell; the weights of the cell w...
bool active() const
Return true, if this grid is active with respect to the last parameter point passed.
double varianceOfAverage() const
Return the variance of the weights.
SimpleCellGrid * selectCell(RndGenerator &rnd)
Select a cell.
const std::vector< double > & upperRight() const
Return the upper right corner of the cell grid.
Definition: CellGrid.h:92
void book(double weight)
Book a point.
double maxWeight
The maximum weight.
SimpleCellGrid & secondChild()
Access the second child.
std::vector< std::pair< Counter, Counter > > theWeightInformation
Weight information for adaption steps.
bool isLeaf() const
Return true, if this is a leaf in the tree.
Definition: CellGrid.h:108