herwig is hosted by Hepforge, IPPP Durham
Herwig 7.3.0
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
15namespace 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),
105 maxWeight(0.0) {}
106
110 double nPoints;
111
116
121
125 double maxWeight;
126
130 void book(double weight) {
131 nPoints += 1.0;
132 sumOfWeights += std::abs(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 ?
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
A binary cell grid.
Definition: CellGrid.h:42
const std::vector< double > & lowerLeft() const
Return the lower left corner of the cell grid.
Definition: CellGrid.h:87
const CellGrid & firstChild() const
Return the first child.
double integral() const
Return the integral relevant for the last parameter point set.
bool isLeaf() const
Return true, if this is a leaf in the tree.
Definition: CellGrid.h:108
bool active() const
Return true, if this grid is active with respect to the last parameter point passed.
std::size_t size() const
Return the number of nodes contained in this grid.
double weight() const
Return the weight.
Definition: CellGrid.h:178
const std::vector< double > & upperRight() const
Return the upper right corner of the cell grid.
Definition: CellGrid.h:92
const CellGrid & secondChild() const
Return the second child.
A simple cell grid providing basic adaption and sampling.
double getReferenceWeight() const
Return the reference weight.
SimpleCellGrid & secondChild()
Access the second child.
void sampleFlatPoint(std::vector< double > &p, const std::vector< bool > &parameterFlags, RndGenerator &rnd) const
Sample a point flat in this cell, keeping parameters fixed.
virtual CellGrid * makeInstance(const std::vector< double > &newLowerLeft, const std::vector< double > &newUpperRight, double newWeight=0.0) const
Produce a new instance of a cell grid.
std::pair< double, double > generate(RndGenerator &rnd, Function &f, std::vector< double > &p)
Sample a point and return its weight.
SimpleCellGrid()
Default constructor.
SimpleCellGrid & firstChild()
Access the first child.
const std::vector< std::pair< Counter, Counter > > & weightInformation() const
Return weight information for adaption steps.
const SimpleCellGrid & firstChild() const
Return the first child.
void sampleFlatPoint(std::vector< double > &p, RndGenerator &rnd) const
Sample a point flat in this cell.
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...
void adjustReferenceWeight(double w)
Adjust the reference weight.
std::vector< std::pair< Counter, Counter > > & weightInformation()
Access weight information for adaption steps.
virtual void fromXML(const XML::Element &)
Fill CellGrid data from an XML element.
virtual CellGrid * makeInstance() const
Produce a new instance of a cell grid.
double sample(RndGenerator &rnd, Function &f, std::vector< double > &p, bool unweight, bool adjustReference)
Sample a point and return its weight.
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 theReferenceWeight
The reference weight to be used for unweighting.
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 ...
virtual void setWeights()
Update the weights of the cells from information accumulated so far.
SimpleCellGrid(const std::vector< double > &newLowerLeft, const std::vector< double > &newUpperRight, bool keepWeightInformation=true, double newWeight=0.0)
Construct given boundaries and a weight.
virtual XML::Element toXML() const
Return an XML element for the data of this CellGrid.
virtual void split(std::size_t newSplitDimension, double newSplitCoordinate)
Split this cell grid in the given dimension and coordinate, if it is a leaf.
virtual void updateWeightInformation(const std::vector< double > &p, double w)
Update the weight information for the given point.
const SimpleCellGrid & secondChild() const
Return the second child.
std::vector< std::pair< Counter, Counter > > theWeightInformation
Weight information for adaption steps.
SimpleCellGrid * selectCell(RndGenerator &rnd)
Select a cell.
Element represents a (tree of) XML elements.
Definition: Element.h:56
constexpr auto sqr(const T &x) -> decltype(x *x)
A simple counter to store information used for adaption.
double averageWeight() const
Return the average weight.
void book(double weight)
Book a point.
double maxWeight
The maximum weight.
double sumOfWeights
The sum of weights.
double nPoints
The number of points.
double sumOfSquaredWeights
The sum of squared weights.
double varianceOfAverage() const
Return the variance of the weights.