Page MenuHomeHEPForge
Files F8724323

No OneTemporary
Actions

View Options

diff --git a/Exsample2/Makefile.am b/Exsample2/Makefile.am
new file mode 100644
--- /dev/null
+++ b/Exsample2/Makefile.am
@@ -0,0 +1,9 @@
+# SUBDIRS = exsample
+
+pkglib_LTLIBRARIES = HwExsample2.la
+HwExsample2_la_LDFLAGS = -module -version-info 11:1:0
+HwExsample2_la_SOURCES = \
+exsampler.h exsampler.cc
+
+# HwExsample2_la_LIBADD = \
+# exsample/libexsample.la
diff --git a/Exsample2/exsample/Makefile.am b/Exsample2/exsample/Makefile.am
new file mode 100644
--- /dev/null
+++ b/Exsample2/exsample/Makefile.am
@@ -0,0 +1,19 @@
+noinst_LTLIBRARIES = libexsample.la
+libexsample_la_SOURCES = \
+adaption_info.h \
+adaption_info.icc \
+binary_tree.h \
+cell.h \
+cell.icc \
+config.h \
+exponential_generator.h \
+exponential_generator.icc \
+generator.h \
+generator.icc \
+linear_interpolator.h \
+linear_interpolator.icc \
+selectors.h \
+selectors.icc \
+statistics.h \
+statistics.icc \
+utility.h
diff --git a/Exsample2/exsample/adaption_info.h b/Exsample2/exsample/adaption_info.h
new file mode 100644
--- /dev/null
+++ b/Exsample2/exsample/adaption_info.h
@@ -0,0 +1,72 @@
+// -*- C++ -*-
+//
+// adaption_info.h is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
+//
+// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
+//
+// ExSample is licenced under version 2 of the GPL, see COPYING for details.
+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
+//
+//
+#ifndef EXSAMPLE_adaption_info_h_included
+#define EXSAMPLE_adaption_info_h_included
+
+namespace exsample {
+
+ /// \brief adaption_info is a container for parameters relevant to
+ /// sampling and adaption.
+ struct adaption_info {
+
+ /// the default constructor
+ adaption_info()
+ : dimension(0), lower_left(), upper_right(),
+ presampling_points(100000),
+ histo_depth(2), adapt(), freeze_grid(0),
+ maxtry(200000), efficiency_threshold(.9),
+ gain_threshold(.1) {}
+
+ /// the phasespace dimension
+ std::size_t dimension;
+
+ /// the lower left corner of the function's support
+ std::vector<double> lower_left;
+
+ /// the upper right corner of the function's support
+ std::vector<double> upper_right;
+
+ /// the number of presampling points
+ unsigned long presampling_points;
+
+ /// use 2^histo_depth bins in efficiency histograms
+ std::size_t histo_depth;
+
+ /// indicate which dimensions should be adapted
+ std::vector<bool> adapt;
+
+ /// the number of accepted events after the grid is frozen
+ unsigned long freeze_grid;
+
+ /// the maximum number of misses allowed
+ unsigned long maxtry;
+
+ /// the efficiency threshold below which splits are considered
+ double efficiency_threshold;
+
+ /// a minimum gain for splits to be performed
+ double gain_threshold;
+
+ /// put to ostream
+ template<class OStream>
+ void put(OStream& os) const;
+
+ /// get from istream
+ template<class IStream>
+ void get(IStream& is);
+
+ };
+
+}
+
+#include "adaption_info.icc"
+
+#endif // EXSAMPLE_adaption_info_h_included
diff --git a/Exsample2/exsample/adaption_info.icc b/Exsample2/exsample/adaption_info.icc
new file mode 100644
--- /dev/null
+++ b/Exsample2/exsample/adaption_info.icc
@@ -0,0 +1,67 @@
+// -*- C++ -*-
+//
+// adaption_info.icc is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
+//
+// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
+//
+// ExSample is licenced under version 2 of the GPL, see COPYING for details.
+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
+//
+//
+namespace exsample {
+
+ template<class OStream>
+ void adaption_info::put(OStream& os) const {
+ os << dimension;
+ ostream_traits<OStream>::separator(os);
+ for (std::size_t k = 0; k < dimension; ++k) {
+ os << lower_left[k];
+ ostream_traits<OStream>::separator(os);
+ }
+ for (std::size_t k = 0; k < dimension; ++k) {
+ os << upper_right[k];
+ ostream_traits<OStream>::separator(os);
+ }
+ os << presampling_points;
+ ostream_traits<OStream>::separator(os);
+ os << histo_depth;
+ ostream_traits<OStream>::separator(os);
+ for (std::size_t k = 0; k < dimension; ++k) {
+ os << adapt[k];
+ ostream_traits<OStream>::separator(os);
+ }
+ os << freeze_grid;
+ ostream_traits<OStream>::separator(os);
+ os << maxtry;
+ ostream_traits<OStream>::separator(os);
+ os << efficiency_threshold;
+ ostream_traits<OStream>::separator(os);
+ os << gain_threshold;
+ ostream_traits<OStream>::separator(os);
+ }
+
+ template<class IStream>
+ void adaption_info::get(IStream& is) {
+ is >> dimension;
+ lower_left.resize(dimension);
+ upper_right.resize(dimension);
+ adapt.resize(dimension);
+ for (std::size_t k = 0; k < dimension; ++k) {
+ is >> lower_left[k];
+ }
+ for (std::size_t k = 0; k < dimension; ++k) {
+ is >> upper_right[k];
+ }
+ is >> presampling_points >> histo_depth;
+ for (std::size_t k = 0; k < dimension; ++k) {
+ // strange. is >> adapt[k]
+ // gives: warning: value computet is not used
+ // with gcc 4.2.1
+ bool v; is >> v;
+ adapt[k] = v;
+ }
+ is >> freeze_grid >> maxtry
+ >> efficiency_threshold >> gain_threshold;
+ }
+
+}
diff --git a/Exsample2/exsample/binary_tree.h b/Exsample2/exsample/binary_tree.h
new file mode 100644
--- /dev/null
+++ b/Exsample2/exsample/binary_tree.h
@@ -0,0 +1,887 @@
+// -*- C++ -*-
+//
+// binary_tree.h is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
+//
+// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
+//
+// ExSample is licenced under version 2 of the GPL, see COPYING for details.
+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
+//
+//
+#ifndef EXSAMPLE_binary_tree_h_included
+#define EXSAMPLE_binary_tree_h_included
+
+#include "utility.h"
+
+namespace exsample {
+
+ /// \brief binary_tree represents a binary tree with the ability to
+ /// `cascade' visitor objects down the tree
+ template<class Value>
+ class binary_tree {
+
+ public:
+
+ ///@name type definitions
+ //@{
+
+ /// define the object type
+ typedef Value value_type;
+
+ //@}
+
+ public:
+
+ ///@name constructors
+ //@{
+
+ /// default constructor
+ binary_tree()
+ : neighbours_(0,0),
+ parent_(0), value_(),
+ children_()
+ { }
+
+ /// construct giving key/cell and parent
+ binary_tree(const value_type& thevalue,
+ binary_tree * theparent = 0)
+ : neighbours_(0,0), parent_(theparent),
+ value_(new value_type(thevalue)),
+ children_()
+ { }
+
+ /// binary_tree has a strict ownership; on copying
+ /// binary trees ownership is transferred
+ binary_tree(const binary_tree& x)
+ : neighbours_(x.neighbours_),
+ parent_(x.parent_), value_(),
+ children_() {
+ assert(x.root());
+ binary_tree& nc_x = const_cast<binary_tree&>(x);
+ value_.swap(nc_x.value_);
+ children_.first.swap(nc_x.children_.first);
+ children_.second.swap(nc_x.children_.second);
+ nc_x.parent_ = 0;
+ nc_x.neighbours_.first = 0;
+ nc_x.neighbours_.second = 0;
+ }
+
+ /// binary_tree has a strict ownership; on copying
+ /// binary trees ownership is transferred
+ binary_tree& operator=(const binary_tree& x) {
+ if (this == &x)
+ return *this;
+ assert(x.root());
+ binary_tree& nc_x = const_cast<binary_tree&>(x);
+ value_.swap(nc_x.value_);
+ children_.first.swap(nc_x.children_.first);
+ children_.second.swap(nc_x.children_.second);
+ neighbours_ = x.neighbours_;
+ parent_ = x.parent_;
+ nc_x.parent_ = 0;
+ nc_x.neighbours_.first = 0;
+ nc_x.neighbours_.second = 0;
+ return *this;
+ }
+
+ //@}
+
+ ///@name standard-conforming leaf iterators
+ //@{
+
+ public:
+
+ class const_iterator;
+
+ /// iterator
+ class iterator {
+
+ public:
+
+ ///@name type definitions for iterator traits
+ //@{
+
+ /// define the iterator category
+ typedef std::bidirectional_iterator_tag iterator_category;
+
+ /// define the difference_type
+ typedef int difference_type;
+
+ /// define the value type
+ typedef Value value_type;
+
+ /// define the reference type
+ typedef value_type& reference;
+
+ /// define the pointer type
+ typedef value_type * pointer;
+
+ //@}
+
+ public:
+
+ ///@name constructors
+ //@{
+
+ /// default constructor
+ iterator() : pointee(0), post_end(0), pre_begin(0) { }
+
+ /// constructor taking pointee
+ iterator(binary_tree * p, std::size_t end = 0)
+ : pointee(p), post_end(end), pre_begin(0) { }
+
+ //@
+
+ public:
+
+ ///@name comparisons
+ //@{
+
+ /// comparison
+ bool operator==(const iterator& x) const {
+ return ((pointee == x.pointee) &&
+ (post_end == x.post_end) &&
+ (pre_begin == x.pre_begin));
+ }
+
+ /// comparison
+ bool operator!=(const iterator& x) const { return !(*this == x); }
+
+ //@}
+
+ public:
+
+ ///@name derefrence and indirection
+ //@{
+
+ /// dereference
+ reference operator*() { return pointee->value(); }
+
+ /// indirection
+ pointer operator->() { return &**this; }
+
+ /// return reference to the node
+ binary_tree& node() { return *pointee; }
+
+ //@}
+
+ /// return raw pointer to the element pointed to
+ binary_tree * get() const { return pointee; }
+
+ ///@name biderectional iterator increment/decrements
+ //@{
+
+ /// pre-increment
+ iterator& operator++() {
+ if (post_end) { ++post_end; return *this; }
+ if (pre_begin) { --pre_begin; return *this; }
+ if(!(pointee->right_neighbour())) { post_end = 1; return *this; }
+ pointee = pointee->right_neighbour();
+ return *this;
+ }
+
+ /// pre-decrement
+ iterator& operator--() {
+ if (post_end) { --post_end; return *this; }
+ if (pre_begin) { ++pre_begin; return *this; }
+ if(!(pointee->left_neighbour())) { pre_begin = 1; return *this; }
+ pointee = pointee->left_neighbour();
+ return *this;
+ }
+
+ /// post-increment
+ template<typename Dummy>
+ iterator operator++(Dummy) {
+ iterator tmp = *this;
+ ++(*this);
+ return tmp;
+ }
+
+ /// post-decrement
+ template<typename Dummy>
+ iterator operator--(Dummy) {
+ iterator tmp = *this;
+ --(*this);
+ return tmp;
+ }
+
+ //@}
+
+ private:
+
+ /// friend for conversion
+ friend class const_iterator;
+
+ /// the node pointed to
+ binary_tree * pointee;
+
+ /// the distance from --end() (if above --end())
+ std::size_t post_end;
+
+ /// the distance from begin() (if below begin())
+ std::size_t pre_begin;
+
+ };
+
+ /// return begin iterator
+ iterator begin() { return iterator(left_most()); }
+
+ /// return end iterator
+ iterator end() { return iterator(right_most(),1); }
+
+ /// return global begin iterator
+ iterator global_begin() {
+ if (!root())
+ return parent().global_begin();
+ return iterator(left_most());
+ }
+
+ /// return global end iterator
+ iterator global_end() {
+ if (!root())
+ return parent().global_end();
+ return iterator(right_most(),1);
+ }
+
+ /// const_iterator
+ class const_iterator {
+
+ public:
+
+ ///@name type definitions for iterator traits
+ //@{
+
+ /// define the iterator category
+ typedef std::bidirectional_iterator_tag iterator_category;
+
+ /// define the difference type
+ typedef int difference_type;
+
+ /// define the value type
+ typedef const Value value_type;
+
+ /// define the reference type
+ typedef const value_type& reference;
+
+ /// define the pointer type
+ typedef const value_type * pointer;
+
+ //@}
+
+ public:
+
+ ///@name constructors
+ //@{
+
+ /// default constructor
+ const_iterator() : pointee(0), post_end(0), pre_begin(0) { }
+
+ /// constructor taking pointee
+ const_iterator(binary_tree * p, std::size_t end = 0)
+ : pointee(p), post_end(end), pre_begin(0) { }
+
+ /// conversion from iterator
+ const_iterator(const iterator& x)
+ : pointee(x.pointee), post_end(x.post_end), pre_begin(x.pre_begin) { }
+
+ //@}
+
+ public:
+
+ ///@name comparisons
+ //@{
+
+ /// comparison
+ bool operator==(const const_iterator& x) const {
+ return ((pointee == x.pointee) &&
+ (post_end == x.post_end) &&
+ (pre_begin == x.pre_begin));
+ }
+
+ /// comparison
+ bool operator!=(const const_iterator& x) const { return !(*this == x); }
+
+ //@}
+
+ public:
+
+ ///@name dereference and indirection
+ //@{
+
+ /// dereference
+ reference operator*() const { return pointee->value(); }
+
+ /// indirection
+ pointer operator->() const { return &**this; }
+
+ /// return reference to the node
+ const binary_tree& node() const { return *pointee; }
+
+ //@}
+
+ ///@name biderectional iterator increment/decrements
+ //@{
+
+ /// pre-increment
+ const_iterator& operator++() {
+ if (post_end) { ++post_end; return *this; }
+ if (pre_begin) { --pre_begin; return *this; }
+ if(!(pointee->right_neighbour())) { post_end = 1; return *this; }
+ pointee = pointee->right_neighbour();
+ return *this;
+ }
+
+ /// pre-decrement
+ const_iterator& operator--() {
+ if (post_end) { --post_end; return *this; }
+ if (pre_begin) { ++pre_begin; return *this; }
+ if(!(pointee->left_neighbour())) { pre_begin = 1; return *this; }
+ pointee = pointee->left_neighbour();
+ return *this;
+ }
+
+ /// post-increment
+ template<typename Dummy>
+ const_iterator operator++(Dummy) {
+ const_iterator tmp = *this;
+ ++(*this);
+ return tmp;
+ }
+
+ /// post-decrement
+ template<typename Dummy>
+ const_iterator operator--(Dummy) {
+ const_iterator tmp = *this;
+ --(*this);
+ return tmp;
+ }
+
+ //@}
+
+ private:
+
+ /// the node pointed to
+ binary_tree * pointee;
+
+ /// the distance from --end() (if above --end())
+ std::size_t post_end;
+
+ /// the distance from begin() (if below begin())
+ std::size_t pre_begin;
+
+ };
+
+ /// return begin const_iterator
+ const_iterator begin() const { return const_iterator(left_most()); }
+
+ /// return end const_iterator
+ const_iterator end() const { return const_iterator(right_most(),1); }
+
+ /// return global begin iterator
+ const_iterator global_begin() const {
+ if (!root())
+ return parent().global_begin();
+ return iterator(left_most());
+ }
+
+ /// return global end iterator
+ const_iterator global_end() const {
+ if (!root())
+ return parent().global_end();
+ return iterator(right_most(),1);
+ }
+
+ private:
+
+ /// set the left neighbour
+ void left_neighbour(binary_tree * n) { neighbours_.first = n; }
+
+ /// set the right neighbour
+ void right_neighbour(binary_tree * n) { neighbours_.second = n; }
+
+ /// get the left neighbour
+ binary_tree * left_neighbour() const { return neighbours_.first; }
+
+ /// get the right neighbour
+ binary_tree * right_neighbour() const { return neighbours_.second; }
+
+ /// return the left-most leaf
+ binary_tree * left_most() {
+ if(leaf()) return this;
+ return left_child().left_most();
+ }
+
+ /// return the right-most leaf
+ binary_tree * right_most() {
+ if(leaf()) return this;
+ return right_child().right_most();
+ }
+
+ /// return the left-most leaf
+ const binary_tree * left_most() const {
+ if(leaf()) return this;
+ return left_child().left_most();
+ }
+
+ /// return the right-most leaf
+ const binary_tree * right_most() const {
+ if(leaf()) return this;
+ return right_child().right_most();
+ }
+
+ /// the iterator is a good friend
+ friend class binary_tree<value_type>::iterator;
+
+ /// the iterator is a good friend
+ friend class binary_tree<value_type>::const_iterator;
+
+ /// the left and right neighbours of this node
+ std::pair<binary_tree*,binary_tree*> neighbours_;
+
+ //@}
+
+ public:
+
+ /// return true, if this node is empty
+ bool empty() const { return root() && leaf() && !value_; }
+
+ /// clear this node
+ void clear() {
+ neighbours_ = std::make_pair<binary_tree*,binary_tree*>(0,0);
+ parent_ = 0;
+ value_.reset(0);
+ if (!leaf()) {
+ left_child().clear();
+ right_child().clear();
+ }
+ children_.first.reset(0);
+ children_.second.reset(0);
+ }
+
+ public:
+
+ /// split this node
+ std::pair<iterator,iterator> split(std::pair<value_type,value_type> children) {
+
+ assert(leaf());
+
+ children_.first.reset(new binary_tree(children.first,this));
+ children_.second.reset(new binary_tree(children.second,this));
+
+ children_.first->left_neighbour(neighbours_.first);
+ children_.first->right_neighbour(children_.second.get());
+ children_.second->left_neighbour(children_.first.get());
+ children_.second->right_neighbour(neighbours_.second);
+
+ // adjust original neighbours
+
+ if(neighbours_.first) {
+ neighbours_.first->right_neighbour(children_.first.get());
+ }
+
+ if (neighbours_.second) {
+ neighbours_.second->left_neighbour(children_.second.get());
+ }
+
+ neighbours_.first = 0; neighbours_.second = 0;
+
+ return std::make_pair(iterator(children_.first.get()),iterator(children_.second.get()));
+
+ }
+
+ public:
+
+ /// select using a selector
+ template<class Selector>
+ iterator select(const Selector& selector) {
+
+ if(leaf()) {
+ bool use = selector.use(value());
+ if (use) return iterator(this);
+ return global_end();
+ }
+
+ std::pair<bool,bool> which(selector.use(value(),left_child().value(),right_child().value()));
+
+ assert(!which.first || !which.second);
+
+ if (!which.first && !which.second) {
+ return global_end();
+ }
+
+ if (which.first) {
+ return left_child().select(selector);
+ }
+ else {
+ return right_child().select(selector);
+ }
+
+ return global_end();
+
+ }
+
+ /// generate a hash value for the sub-tree
+ /// selected by the given selector object
+ template<class Selector, unsigned long bits>
+ void subtree_hash(const Selector& selector, bit_container<bits>& bhash) {
+ bhash = bit_container<bits>();
+ unsigned long pos = 0;
+ do_subtree_hash<Selector,bits>(selector,bhash,pos);
+ }
+
+ /// accumulate values using a binary function
+ /// and accessor object
+ template<class Accessor, class BinaryOp>
+ typename BinaryOp::result_type accumulate(const Accessor& acc,
+ BinaryOp binary_op) const {
+
+ if (!leaf()) {
+ return
+ binary_op(left_child().accumulate(acc,binary_op),
+ right_child().accumulate(acc,binary_op));
+ }
+
+ return acc.get(value(),true);
+
+ }
+
+ /// accumulate values only from branches
+ /// matching a Selector
+ template<class Selector, class Accessor, class BinaryOp>
+ typename BinaryOp::result_type accumulate(const Selector& selector,
+ const Accessor& acc,
+ BinaryOp binary_op) const {
+
+ if (!leaf()) {
+ std::pair<bool,bool> which(selector.use(value(),left_child().value(),right_child().value()));
+ assert(which.first || which.second);
+ if (which.first && which.second) {
+ return
+ binary_op(left_child().accumulate(selector,acc,binary_op),
+ right_child().accumulate(selector,acc,binary_op));
+ } else if (which.first) {
+ return left_child().accumulate(selector,acc,binary_op);
+ } else if (which.second) {
+ return right_child().accumulate(selector,acc,binary_op);
+ }
+ }
+
+ return acc.get(value(),true);
+
+ }
+
+ /// accumulate values using a binary function
+ /// and accessor object, storing intermediate
+ /// values in nodes
+ template<class Accessor, class BinaryOp>
+ typename BinaryOp::result_type tree_accumulate(const Accessor& acc,
+ BinaryOp binary_op) {
+
+ if (!leaf()) {
+ acc.set(value()) =
+ binary_op(left_child().tree_accumulate(acc,binary_op),
+ right_child().tree_accumulate(acc,binary_op));
+ return acc.get(value(),false);
+ }
+
+ acc.set(value()) = acc.get(value(),true);
+ return acc.get(value(),true);
+
+ }
+
+ /// accumulate values only from branches
+ /// matching a Selector
+ template<class Selector, class Accessor, class BinaryOp>
+ typename BinaryOp::result_type tree_accumulate(const Selector& selector,
+ const Accessor& acc,
+ BinaryOp binary_op) {
+
+ if (!leaf()) {
+ std::pair<bool,bool> which(selector.use(value(),left_child().value(),right_child().value()));
+ assert(which.first || which.second);
+ if (which.first && which.second) {
+ acc.set(value()) =
+ binary_op(left_child().tree_accumulate(selector,acc,binary_op),
+ right_child().tree_accumulate(selector,acc,binary_op));
+ } else if (which.first) {
+ acc.set(value()) = left_child().tree_accumulate(selector,acc,binary_op);
+ } else if (which.second) {
+ acc.set(value()) = right_child().tree_accumulate(selector,acc,binary_op);
+ }
+ return acc.get(value(),false);
+ }
+
+ acc.set(value()) = acc.get(value(),true);
+ return acc.get(value(),true);
+
+ }
+
+ /// forward propagate a visitor to all children nodes
+ template<class Visitor>
+ void cascade(Visitor visitor) const {
+ if (leaf()) {
+ visitor.visit(value());
+ return;
+ } else visitor.visit(value(),left_child().value(),right_child().value());
+ left_child().cascade(visitor);
+ right_child().cascade(visitor);
+ }
+
+ /// succesively split using a generator
+ template<class Generator>
+ void generate(Generator generator) {
+ if (root())
+ value_.reset(new value_type(generator.root()));
+ if (generator.split()) {
+ std::pair<iterator,iterator> ch = split(generator.generate(value()));
+ ch.first.node().generate(generator);
+ ch.second.node().generate(generator);
+ }
+ }
+
+ public:
+
+ ///@name Public member access
+ //@{
+
+ /// return the value held by this node
+ value_type& value() { return *value_; }
+
+ /// return the value held by this node
+ const value_type& value() const { return *value_; }
+
+ /// return true, if this is the root node
+ bool root() const { return !parent_; }
+
+ /// return true, if this node has got children
+ bool leaf() const { return !(children_.first.get() && children_.second.get()); }
+
+ //@}
+
+ public:
+
+ ///@name put and get from streams
+ //@{
+
+ /// forward visitor writing out the tree to given ostream
+ template<class OStream>
+ struct ostream_visitor {
+
+ /// construct from ostream reference
+ explicit ostream_visitor(OStream& os) : os_(&os), first_time_(true) {}
+
+ /// visit a leaf node
+ void visit(const value_type&) {
+ (*os_) << "end_branch";
+ ostream_traits<OStream>::separator(*os_);
+ }
+
+ /// visit a branching
+ void visit(const value_type& parent,
+ const value_type& left, const value_type& right) {
+ if (first_time_) {
+ (*os_) << "root_node";
+ ostream_traits<OStream>::separator(*os_);
+ parent.put(*os_);
+ first_time_ = false;
+ }
+
+ (*os_) << "left_child";
+ ostream_traits<OStream>::separator(*os_);
+ left.put(*os_);
+
+ (*os_) << "right_child";
+ ostream_traits<OStream>::separator(*os_);
+ right.put(*os_);
+
+ }
+
+ private:
+
+ /// pointer to the ostream to write to
+ OStream* os_;
+
+ /// whether we are at the or not
+ bool first_time_;
+
+ };
+
+ /// generator reading binary tree from istream
+ template<class IStream>
+ struct istream_generator {
+
+ /// construct from istream reference
+ explicit istream_generator(IStream& is)
+ : is_(&is), children_(), tag_("") {}
+
+ /// copy constructor
+ istream_generator(const istream_generator& x)
+ : is_(x.is_), children_(), tag_("") {}
+
+ /// read the root node
+ value_type root() {
+
+ *is_ >> tag_;
+ assert(tag_ == "root_node");
+ value_type rnode;
+ rnode.get(*is_);
+
+ return rnode;
+ }
+
+ /// read children nodes
+ bool split() {
+
+ *is_ >> tag_;
+
+ if (tag_ == "end_branch") {
+ return false;
+ }
+
+ assert (tag_ == "left_child");
+ children_.first.get(*is_);
+
+ *is_ >> tag_;
+ assert(tag_ == "right_child");
+ children_.second.get(*is_);
+
+ return true;
+
+ }
+
+ /// return the children generated
+ std::pair<value_type,value_type> generate(const value_type&) {
+ return children_;
+ }
+
+ /// initialize a leaf
+ void initialize_leaf(const value_type&) {}
+
+ private:
+
+ /// pointer to the istream used
+ IStream* is_;
+
+ /// the children currently handled
+ std::pair<value_type,value_type> children_;
+
+ /// temporary storage for tags
+ std::string tag_;
+
+ };
+
+ /// put to ostream
+ template<class OStream>
+ void put(OStream& os) const {
+
+ if (empty()) {
+ os << "empty";
+ ostream_traits<OStream>::separator(os);
+ return;
+ } else if (root() && leaf()) {
+ os << "root_only";
+ ostream_traits<OStream>::separator(os);
+ value().put(os);
+ return;
+ } else {
+ os << "non_empty";
+ ostream_traits<OStream>::separator(os);
+ }
+
+ assert(root());
+ cascade(ostream_visitor<OStream>(os));
+
+ }
+
+ /// get from istream
+ template<class IStream>
+ void get(IStream& is) {
+
+ std::string state;
+ is >> state;
+
+ if (state == "empty") {
+ return;
+ }
+ if (state == "root_only") {
+ value_.reset(new value_type());
+ value().get(is);
+ return;
+ }
+
+ assert(empty());
+ generate(istream_generator<IStream>(is));
+
+ }
+
+ //@}
+
+
+ private:
+
+ /// calculate hash value
+ template<class Selector, unsigned long bits>
+ void do_subtree_hash(const Selector& selector,
+ bit_container<bits>& current,
+ unsigned long& position,
+ bool selected = true) const {
+
+ if (!leaf()) {
+
+ std::pair<bool,bool> which(false,false);
+ if (selected)
+ which = selector.use(value(),left_child().value(),right_child().value());
+
+ current.bit(position,which.first);
+ current.bit(position+1,which.second);
+
+ position += 2;
+
+ left_child().do_subtree_hash(selector,current,position,which.first && selected);
+ right_child().do_subtree_hash(selector,current,position,which.second && selected);
+ }
+
+ }
+
+ private:
+
+ ///@name private member access
+ //@{
+
+ /// return the parent of this node
+ binary_tree& parent() { assert(parent_); return *parent_; }
+
+ /// return the parent of this node
+ const binary_tree& parent() const { assert(parent_); return *parent_; }
+
+ /// return the left child of this node
+ binary_tree& left_child() { assert(children_.first.get()); return *children_.first; }
+
+ /// return the left child of this node
+ const binary_tree& left_child() const { assert(children_.first.get()); return *children_.first; }
+
+ /// return the right child of this node
+ binary_tree& right_child() { assert(children_.second.get()); return *children_.second; }
+
+ /// return the right child of this node
+ const binary_tree& right_child() const { assert(children_.second.get()); return *children_.second; }
+
+ //@}
+
+ private:
+
+ /// the parent of this node
+ binary_tree * parent_;
+
+ /// the cell held by this node
+ boost::scoped_ptr<value_type> value_;
+
+ /// the children of this node
+ std::pair<boost::scoped_ptr<binary_tree>,
+ boost::scoped_ptr<binary_tree> > children_;
+
+ };
+
+}
+
+#endif // EXSAMPLE_binary_tree_h_included
diff --git a/Exsample2/exsample/cell.h b/Exsample2/exsample/cell.h
new file mode 100644
--- /dev/null
+++ b/Exsample2/exsample/cell.h
@@ -0,0 +1,304 @@
+// -*- C++ -*-
+//
+// cell.h is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
+//
+// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
+//
+// ExSample is licenced under version 2 of the GPL, see COPYING for details.
+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
+//
+//
+#ifndef EXSAMPLE_cell_h_included
+#define EXSAMPLE_cell_h_included
+
+#include "utility.h"
+#include "adaption_info.h"
+#include "statistics.h"
+
+namespace exsample {
+
+ /// \brief Information contained in a leaf cell
+ class cell_info {
+
+ public:
+
+ /// the default constructor
+ cell_info();
+
+ /// construct from boundaries and adaption info
+ cell_info(const std::vector<double>& ll,
+ const std::vector<double>& ur,
+ const adaption_info& ainfo);
+
+ /// construct from boundaries, flags for variables to be sampled,
+ /// and adaption info
+ cell_info(const std::vector<double>& ll,
+ const std::vector<double>& ur,
+ const std::vector<bool>& sampled_variables,
+ const adaption_info& ainfo);
+
+ public:
+
+ /// generate a flat trial point in this cell
+ template<class Random>
+ void select(Random&,
+ std::vector<double>&);
+
+ /// generate a flat trial point in this cell
+ /// only for the variables falgged as true
+ template<class Random>
+ void select(Random&,
+ std::vector<double>&,
+ const std::vector<bool>&);
+
+ /// indicate a function value for the given point
+ void selected(const std::vector<double>&,
+ double,
+ const adaption_info&);
+
+ /// indicate that a point has been
+ /// accepted in this cell
+ void accept() { ++accepted_; }
+
+ /// reject a previously accepted event
+ void reject() { --accepted_; }
+
+ public:
+
+ /// return true, if below efficiency threshold
+ bool bad(const adaption_info& ainfo) const {
+ return ((static_cast<double>(accepted_)/static_cast<double>(attempted_)) <
+ ainfo.efficiency_threshold);
+ }
+
+ /// suggest a split and indicate wether it is worth
+ /// to be performed
+ std::pair<std::size_t,double> get_split(const adaption_info&,
+ bool&) const;
+
+ /// explore this cell performing a flat sampling,
+ /// updating the given statistics object and pre-filling
+ /// the efficiency histogram by a trial unweighting
+ template<class Random, class Function, class SlaveStatistics>
+ void explore(Random&, const adaption_info&, Function*, statistics*,
+ SlaveStatistics& opt);
+
+ /// explore this cell in a more refined way, which
+ /// is however not suited for already calculating integrals
+ /// and stuff
+ template<class Random, class Function>
+ void explore(Random&, const adaption_info&, Function*);
+
+ public:
+
+ /// get the current overestimate
+ double overestimate() const { return overestimate_; }
+
+ /// return the position of the last maximum
+ const std::vector<double>& last_max_position() const { return last_max_position_; }
+
+ /// set the current overestimate and maximum position
+ void overestimate(double v, const std::vector<double>& pos) {
+ overestimate_ = v;
+ last_max_position_ = pos;
+ }
+
+ /// get the volume
+ double volume() const { return volume_; }
+
+ /// get the lower left corner
+ const std::vector<double>& lower_left() const { return lower_left_; }
+
+ /// get the upper right corner
+ const std::vector<double>& upper_right() const { return upper_right_; }
+
+ /// get the number of attempted events
+ unsigned long attempted() const { return attempted_; }
+
+ /// get the number of accepted events
+ unsigned long accepted() const { return accepted_; }
+
+ public:
+
+ /// return the number of missing events
+ /// for the given parameter bin id
+ int parametric_missing(const bit_container<parameter_hash_bits>& id) const;
+
+ /// set the number of missing events
+ /// for the given parameter bin id
+ void parametric_missing(const bit_container<parameter_hash_bits>& id, int n);
+
+ /// increase to the number of missing events
+ /// for the given parameter bin id
+ void increase_parametric_missing(const bit_container<parameter_hash_bits>& id);
+
+ /// decrease to the number of missing events
+ /// for the given parameter bin id
+ void decrease_parametric_missing(const bit_container<parameter_hash_bits>& id);
+
+ /// return true, if the cell is compensating in
+ /// at least one parameter bin
+ bool parametric_compensating() const {
+ return !parametric_missing_map_.empty();
+ }
+
+ /// return true, if the cell contains the
+ /// indicated parameter point
+ bool contains_parameter(const std::vector<double>& point,
+ const std::vector<bool>& sampled) const;
+
+ public:
+
+ /// put to ostream
+ template<class OStream>
+ void put(OStream& os) const;
+
+ /// get from istream
+ template<class IStream>
+ void get(IStream& is);
+
+ private:
+
+ /// the value of the overestimate in this cell
+ double overestimate_;
+
+ /// the volume of this cell
+ double volume_;
+
+ /// the lower left corner of this cell
+ std::vector<double> lower_left_;
+
+ /// the upper right corner of this cell
+ std::vector<double> upper_right_;
+
+ /// midpoint of this cell
+ std::vector<double> mid_point_;
+
+ /// the position of the last encountered
+ /// maximum in this cell
+ std::vector<double> last_max_position_;
+
+ /// left-right statistics of average weight
+ std::vector<std::pair<double,double> > avg_weight_;
+
+ /// the number of attempts in this cell
+ unsigned long attempted_;
+
+ /// the number of accepted events in this cell
+ unsigned long accepted_;
+
+ /// an optional map of parameter bin ids
+ /// to the number of missing events
+ std::map<bit_container<parameter_hash_bits>,int> parametric_missing_map_;
+
+ };
+
+ /// \brief the general cell class
+ class cell {
+
+ public:
+
+ /// default constructor
+ cell();
+
+ /// construct from boundaries and adaption info
+ cell(const std::vector<double>& ll,
+ const std::vector<double>& ur,
+ const adaption_info& ainfo);
+
+ /// construct from boundaries, flags for variables to be sampled,
+ /// and adaption info
+ cell(const std::vector<double>& ll,
+ const std::vector<double>& ur,
+ const std::vector<bool>& sampled_variables,
+ const adaption_info& ainfo);
+
+ /// copy constructor
+ cell(const cell& x);
+
+ /// assignment
+ cell& operator=(const cell& x);
+
+ public:
+
+ /// split this cell, exploring the
+ /// child not containing the current overestimate
+ template<class Random, class Function>
+ std::pair<cell,cell> split(std::pair<std::size_t,double> split_d,
+ Random& rnd_gen,
+ Function* f,
+ const adaption_info& ainfo,
+ const std::vector<bool>& sampled =
+ std::vector<bool>());
+
+ public:
+
+ /// return the split dimension
+ std::size_t split_dimension() const { return split_dimension_; }
+
+ /// return the split value
+ double split_point() const { return split_point_; }
+
+ /// return the integral
+ double integral() const { return integral_; }
+
+ /// access the integral
+ double& integral() { return integral_; }
+
+ /// set the integral
+ void integral(double v) { integral_ = v; }
+
+ /// access the number of missing events
+ int& missing_events() { return missing_events_; }
+
+ /// return the number of missing events
+ int missing_events() const { return missing_events_; }
+
+ /// set the number of missing events
+ void missing_events(int n) { missing_events_ = n; }
+
+ /// access the cell_info object
+ cell_info& info() { assert(cell_info_); return *cell_info_; }
+
+ /// return the cell_info object
+ const cell_info& info() const { assert(cell_info_); return *cell_info_; }
+
+ public:
+
+ /// put to ostream
+ template<class OStream>
+ void put(OStream& os) const;
+
+ /// get from istream
+ template<class IStream>
+ void get(IStream& is);
+
+ private:
+
+ /// the dimension along this cell
+ /// was split
+ std::size_t split_dimension_;
+
+ /// the value, where this cell was split
+ double split_point_;
+
+ /// the integral of the absolute value
+ /// of the overestimate over all the
+ /// children cells
+ double integral_;
+
+ /// the number of missing events in this cell
+ int missing_events_;
+
+ /// a pointer to the cell info object,
+ /// if this is a leaf cell
+ boost::scoped_ptr<cell_info> cell_info_;
+
+
+ };
+
+}
+
+#include "cell.icc"
+
+#endif // EXSAMPLE_cell_h_included
diff --git a/Exsample2/exsample/cell.icc b/Exsample2/exsample/cell.icc
new file mode 100644
--- /dev/null
+++ b/Exsample2/exsample/cell.icc
@@ -0,0 +1,409 @@
+// -*- C++ -*-
+//
+// cell.icc is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
+//
+// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
+//
+// ExSample is licenced under version 2 of the GPL, see COPYING for details.
+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
+//
+//
+namespace exsample {
+
+
+ template<class Random>
+ void cell_info::select (Random& rnd_gen,
+ std::vector<double>& p) {
+ std::transform(lower_left_.begin(),lower_left_.end(),
+ upper_right_.begin(),p.begin(),
+ rnd_gen);
+ ++attempted_;
+ }
+
+
+ template<class Random>
+ void cell_info::select (Random& rnd_gen,
+ std::vector<double>& p,
+ const std::vector<bool>& sample) {
+ conditional_transform(lower_left_.begin(),lower_left_.end(),
+ upper_right_.begin(),sample.begin(),
+ p.begin(),rnd_gen);
+ ++attempted_;
+ }
+
+
+
+ template<class Random, class Function, class SlaveStatistics>
+ void cell_info::explore(Random& rnd_gen,
+ const adaption_info& ainfo,
+ Function* function, statistics* stats,
+ SlaveStatistics& opt) {
+ function->start_presampling();
+ unsigned long n_sampled = 0;
+ std::vector<double> ll = lower_left_;
+ std::vector<double> ur = upper_right_;
+ double val = 0.;
+ std::vector<double> pos (ll.size());
+ while (n_sampled < ainfo.presampling_points) {
+ std::transform(ll.begin(),ll.end(),
+ ur.begin(),pos.begin(),
+ rnd_gen);
+ val = function->evaluate(pos);
+ stats->presampled(val);
+ opt.select(val);
+ if (std::abs(val) > std::abs(overestimate_)) {
+ overestimate_ = std::abs(val);
+ last_max_position_ = pos;
+ }
+ selected(pos,std::abs(val),ainfo);
+ ++n_sampled;
+ }
+ function->stop_presampling();
+ }
+
+
+ template<class Random, class Function>
+ void cell_info::explore (Random& rnd_gen,
+ const adaption_info& ainfo, Function* function) {
+ function->start_presampling();
+ unsigned long n_sampled = 0;
+ std::vector<double> ll = lower_left_;
+ std::vector<double> ur = upper_right_;
+ double val = 0.;
+ std::vector<double> pos (ll.size());
+ std::vector<double> next_to_largest (ll.size());
+ while (n_sampled < ainfo.presampling_points) {
+ std::transform(ll.begin(),ll.end(),
+ ur.begin(),pos.begin(),
+ rnd_gen);
+ val = function->evaluate(pos);
+ if (std::abs(val) > std::abs(overestimate_)) {
+ next_to_largest = last_max_position_;
+ overestimate_ = std::abs(val);
+ last_max_position_ = pos;
+ }
+ ++n_sampled;
+ }
+ function->stop_presampling();
+ }
+
+
+ template<class OStream>
+ void cell_info::put (OStream& os) const {
+ os << overestimate_;
+ ostream_traits<OStream>::separator(os);
+ os << volume_;
+ ostream_traits<OStream>::separator(os);
+ os << lower_left_.size();
+ ostream_traits<OStream>::separator(os);
+ for (std::size_t k = 0; k < lower_left_.size(); ++k) {
+ os << lower_left_[k];
+ ostream_traits<OStream>::separator(os);
+ }
+ for (std::size_t k = 0; k < upper_right_.size(); ++k) {
+ os << upper_right_[k];
+ ostream_traits<OStream>::separator(os);
+ }
+ for (std::size_t k = 0; k < mid_point_.size(); ++k) {
+ os << mid_point_[k];
+ ostream_traits<OStream>::separator(os);
+ }
+ for (std::size_t k = 0; k < last_max_position_.size(); ++k) {
+ os << last_max_position_[k];
+ ostream_traits<OStream>::separator(os);
+ }
+ for (std::size_t k = 0; k < avg_weight_.size(); ++k) {
+ os << avg_weight_[k].first;
+ ostream_traits<OStream>::separator(os);
+ os << avg_weight_[k].second;
+ ostream_traits<OStream>::separator(os);
+ }
+ }
+
+
+ template<class IStream>
+ void cell_info::get (IStream& is) {
+ std::size_t dim;
+ is >> overestimate_ >> volume_ >> dim;
+ lower_left_.resize(dim);
+ for (std::size_t k = 0; k < lower_left_.size(); ++k) {
+ is >> lower_left_[k];
+ }
+ upper_right_.resize(dim);
+ for (std::size_t k = 0; k < upper_right_.size(); ++k) {
+ is >> upper_right_[k];
+ }
+ mid_point_.resize(dim);
+ for (std::size_t k = 0; k < mid_point_.size(); ++k) {
+ is >> mid_point_[k];
+ }
+ last_max_position_.resize(dim);
+ for (std::size_t k = 0; k < last_max_position_.size(); ++k) {
+ is >> last_max_position_[k];
+ }
+ avg_weight_.resize(dim);
+ for (std::size_t k = 0; k < avg_weight_.size(); ++k) {
+ is >> avg_weight_[k].first >> avg_weight_[k].second;
+ }
+ }
+
+ template<class Random, class Function>
+ std::pair<cell,cell >
+ cell::split (std::pair<std::size_t,double> split_d,
+ Random& rnd_gen,
+ Function* function,
+ const adaption_info& ainfo,
+ const std::vector<bool>& sampled) {
+ assert(!missing_events() && !info().parametric_compensating());
+ split_dimension_ = split_d.first;
+ split_point_ = split_d.second;
+ std::vector<double> lower_left1 = info().lower_left();
+ std::vector<double> upper_right1 = info().upper_right();
+ std::vector<double> lower_left2 = info().lower_left();
+ std::vector<double> upper_right2 = info().upper_right();
+ upper_right1[split_dimension_] = split_point_;
+ lower_left2[split_dimension_] = split_point_;
+ std::pair<cell,cell> children;
+ if (sampled.empty())
+ children = std::pair<cell,cell>(cell(lower_left1,upper_right1,ainfo),
+ cell(lower_left2,upper_right2,ainfo));
+ else
+ children = std::pair<cell,cell> (cell(lower_left1,upper_right1,sampled,ainfo),
+ cell(lower_left2,upper_right2,sampled,ainfo));
+ if (info().last_max_position()[split_dimension_] <= split_point_) {
+ children.first.info().overestimate(info().overestimate(),info().last_max_position());
+ children.second.info().explore(rnd_gen,ainfo,function);
+ } else {
+ children.second.info().overestimate(info().overestimate(),info().last_max_position());
+ children.first.info().explore(rnd_gen,ainfo,function);
+ }
+ cell_info_.reset(0);
+ children.first.integral(children.first.info().volume() * children.first.info().overestimate());
+ children.second.integral(children.second.info().volume() * children.second.info().overestimate());
+ return children;
+ }
+
+
+ template<class OStream>
+ void cell::put (OStream& os) const {
+ os << split_dimension_;
+ ostream_traits<OStream>::separator(os);
+ os << split_point_;
+ ostream_traits<OStream>::separator(os);
+ os << integral_;
+ ostream_traits<OStream>::separator(os);
+ if (cell_info_) {
+ os << "has_cell_info";
+ ostream_traits<OStream>::separator(os);
+ cell_info_->put(os);
+ } else {
+ os << "has_no_cell_info";
+ ostream_traits<OStream>::separator(os);
+ }
+ }
+
+
+ template<class IStream>
+ void cell::get (IStream& is) {
+ std::string info_tag;
+ is >> split_dimension_ >> split_point_
+ >> integral_ >> info_tag;
+ if (info_tag == "has_cell_info") {
+ cell_info_.reset(new cell_info());
+ cell_info_->get(is);
+ }
+ }
+
+ inline cell_info::cell_info()
+ : overestimate_(0.), volume_(0.),
+ lower_left_(), upper_right_(), mid_point_(),
+ last_max_position_(), avg_weight_(),
+ attempted_(0), accepted_(0) {}
+
+ inline cell_info::cell_info(const std::vector<double>& ll,
+ const std::vector<double>& ur,
+ const adaption_info& ainfo)
+ : overestimate_(0.), volume_(),
+ lower_left_(ll), upper_right_(ur), mid_point_(),
+ last_max_position_(),
+ avg_weight_(std::vector<std::pair<double,double> >
+ (ainfo.dimension,std::make_pair(0.,0.))),
+ attempted_(0), accepted_(0) {
+
+ std::vector<double> delta;
+ std::transform(ur.begin(),ur.end(),
+ ll.begin(),std::back_inserter(delta),
+ std::minus<double>());
+
+ volume_ =
+ std::accumulate(delta.begin(),delta.end(),1.,std::multiplies<double>());
+
+ std::transform(ur.begin(),ur.end(),
+ ll.begin(),std::back_inserter(mid_point_),
+ std::plus<double>());
+
+ for (std::size_t k = 0; k < ainfo.dimension; ++k)
+ mid_point_[k] /= 2.;
+
+ }
+
+ inline cell_info::cell_info(const std::vector<double>& ll,
+ const std::vector<double>& ur,
+ const std::vector<bool>& sampled_variables,
+ const adaption_info& ainfo)
+ : overestimate_(0.), volume_(),
+ lower_left_(ll), upper_right_(ur), mid_point_(),
+ last_max_position_(),
+ avg_weight_(std::vector<std::pair<double,double> >
+ (ainfo.dimension,std::make_pair(0.,0.))),
+ attempted_(0), accepted_(0) {
+
+ std::vector<double> delta;
+ conditional_transform(ur.begin(),ur.end(),
+ ll.begin(),sampled_variables.begin(),
+ std::back_inserter(delta),
+ std::minus<double>());
+
+ volume_ =
+ std::accumulate(delta.begin(),delta.end(),1.,std::multiplies<double>());
+
+ std::transform(ur.begin(),ur.end(),
+ ll.begin(),std::back_inserter(mid_point_),
+ std::plus<double>());
+
+ for (std::size_t k = 0; k < ainfo.dimension; ++k)
+ mid_point_[k] /= 2.;
+
+ }
+
+
+ inline int cell_info::parametric_missing(const bit_container<parameter_hash_bits>& id) const {
+ std::map<bit_container<parameter_hash_bits>,int>::const_iterator mit
+ = parametric_missing_map_.find(id);
+ if (mit == parametric_missing_map_.end())
+ return 0;
+ return mit->second;
+ }
+
+ inline void cell_info::parametric_missing(const bit_container<parameter_hash_bits>& id, int n) {
+ if (n == 0) {
+ std::map<bit_container<parameter_hash_bits>,int>::iterator mit
+ = parametric_missing_map_.find(id);
+ if (mit != parametric_missing_map_.end())
+ parametric_missing_map_.erase(mit);
+ return;
+ }
+ parametric_missing_map_[id] = n;
+ }
+
+ inline void cell_info::increase_parametric_missing(const bit_container<parameter_hash_bits>& id) {
+ std::map<bit_container<parameter_hash_bits>,int>::iterator mit
+ = parametric_missing_map_.find(id);
+ if (mit != parametric_missing_map_.end()) {
+ mit->second += 1;
+ if (mit->second == 0) parametric_missing_map_.erase(mit);
+ } else parametric_missing_map_[id] = 1;
+ }
+
+ inline void cell_info::decrease_parametric_missing(const bit_container<parameter_hash_bits>& id) {
+ std::map<bit_container<parameter_hash_bits>,int>::iterator mit
+ = parametric_missing_map_.find(id);
+ if (mit != parametric_missing_map_.end()) {
+ mit->second -= 1;
+ if (mit->second == 0) parametric_missing_map_.erase(mit);
+ } else assert(false);
+ }
+
+ inline void cell_info::selected(const std::vector<double>& p,
+ double weight,
+ const adaption_info& ainfo) {
+ for (std::size_t k = 0; k < p.size(); ++k) {
+ if (ainfo.adapt[k]) {
+ if (p[k] < mid_point_[k])
+ avg_weight_[k].first += weight;
+ else
+ avg_weight_[k].second += weight;
+ }
+ }
+ }
+
+ inline std::pair<std::size_t,double> cell_info::get_split (const adaption_info& ainfo,
+ bool& worth) const {
+ std::size_t split_d = 0;
+ double gain = 0.;
+ for (std::size_t k = 0; k < ainfo.dimension; ++k) {
+ double xgain = 0.;
+ double left = avg_weight_[k].first;
+ double right = avg_weight_[k].second;
+ if (left+right > 0.) {
+ xgain = std::abs(left-right)/(left+right);
+ }
+ if (xgain > gain) {
+ gain = xgain;
+ split_d = k;
+ }
+ }
+ worth = (gain >= ainfo.gain_threshold);
+ return std::make_pair(split_d,mid_point_[split_d]);
+ }
+
+ inline bool cell_info::contains_parameter (const std::vector<double>& point,
+ const std::vector<bool>& sampled) const {
+ std::vector<double>::const_iterator p = point.begin();
+ std::vector<double>::const_iterator l = lower_left_.begin();
+ std::vector<double>::const_iterator u = upper_right_.begin();
+ std::vector<bool>::const_iterator f = sampled.begin();
+ for (; p < point.end(); ++p, ++f, ++l, ++u)
+ if (!(*f)) {
+ if (((*l) > (*p)) ||
+ ((*u) < (*p)))
+ return false;
+ }
+ return true;
+ }
+
+
+ inline cell::cell()
+ : split_dimension_(0), split_point_(0.),
+ integral_(0.), missing_events_(0),
+ cell_info_(0) {}
+
+ inline cell::cell(const std::vector<double>& ll,
+ const std::vector<double>& ur,
+ const adaption_info& ainfo)
+ : split_dimension_(0), split_point_(0.),
+ integral_(0.), missing_events_(0),
+ cell_info_(new cell_info(ll,ur,ainfo)) {}
+
+ inline cell::cell(const std::vector<double>& ll,
+ const std::vector<double>& ur,
+ const std::vector<bool>& sampled_variables,
+ const adaption_info& ainfo)
+ : split_dimension_(0), split_point_(0.),
+ integral_(0.), missing_events_(0),
+ cell_info_(new cell_info(ll,ur,sampled_variables,ainfo)) {}
+
+ inline cell::cell(const cell& x)
+ : split_dimension_(x.split_dimension_),
+ split_point_(x.split_point_),
+ integral_(x.integral_),
+ missing_events_(x.missing_events_),
+ cell_info_(0) {
+ if (x.cell_info_)
+ cell_info_.reset(new cell_info(*x.cell_info_));
+ }
+
+ inline cell& cell::operator=(const cell& x) {
+ if (this == &x)
+ return *this;
+ split_dimension_ = x.split_dimension_;
+ split_point_ = x.split_point_;
+ integral_ = x.integral_;
+ missing_events_ = x.missing_events_;
+ if (x.cell_info_)
+ cell_info_.reset(new cell_info(*x.cell_info_));
+ return *this;
+ }
+
+
+}
diff --git a/Exsample2/exsample/config.h b/Exsample2/exsample/config.h
new file mode 100644
--- /dev/null
+++ b/Exsample2/exsample/config.h
@@ -0,0 +1,48 @@
+// -*- C++ -*-
+//
+// config.h is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
+//
+// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
+//
+// ExSample is licenced under version 2 of the GPL, see COPYING for details.
+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
+//
+//
+#ifndef EXSAMPLE_config_h_included
+#define EXSAMPLE_config_h_included
+
+#include <vector>
+#include <map>
+#include <set>
+#include <string>
+
+#include <cmath>
+#include <cassert>
+#include <climits>
+
+#include <algorithm>
+#include <numeric>
+#include <limits>
+
+#include <boost/utility.hpp>
+#include <boost/scoped_array.hpp>
+#include <boost/scoped_ptr.hpp>
+
+#define EXSAMPLE_has_ThePEG
+
+#ifdef EXSAMPLE_has_ThePEG
+
+#include "ThePEG/Persistency/PersistentOStream.h"
+
+#endif // EXSAMPLE_has_ThePEG
+
+namespace exsample {
+
+ static const unsigned long parameter_hash_bits = 512;
+
+ using boost::next;
+ using boost::prior;
+
+}
+
+#endif // EXSAMPLE_config_h_included
diff --git a/Exsample2/exsample/exponential_generator.h b/Exsample2/exsample/exponential_generator.h
new file mode 100644
--- /dev/null
+++ b/Exsample2/exsample/exponential_generator.h
@@ -0,0 +1,210 @@
+// -*- C++ -*-
+//
+// exponential_generator.h is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
+//
+// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
+//
+// ExSample is licenced under version 2 of the GPL, see COPYING for details.
+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
+//
+//
+#ifndef EXSAMPLE_exponential_generator_h_included
+#define EXSAMPLE_exponential_generator_h_included
+
+#include "cell.h"
+#include "selectors.h"
+#include "statistics.h"
+#include "linear_interpolator.h"
+#include "binary_tree.h"
+
+namespace exsample {
+
+ /// \brief Exception thrown, if the exponential_generator has just changed its
+ /// state. The attempt of generating an event should be repeated.
+ struct exponential_regenerate{};
+
+ /// \brief The generator for sudakov-type distributions.
+ template<class Function, class Random>
+ class exponential_generator {
+
+ public:
+
+ /// default constructor
+ exponential_generator()
+ : function_(0), check_events_(0), adaption_info_(), root_cell_(),
+ rnd_gen_(), did_split_(false), initialized_(false),
+ evolution_variable_(0), evolution_cutoff_(0.),
+ sample_variables_(), sample_other_variables_(),
+ parameter_splits_(),
+ last_cell_(), last_point_(), last_value_(0.),
+ last_parameter_bin_(), exponents_(),
+ last_exponent_integrand_(),
+ last_exponent_(), compensating_(false),
+ integral_accessor_(), missing_accessor_(),
+ parametric_selector_(), exponent_selector_(),
+ parametric_sampler_(), attempts_(0), accepts_(0),
+ splits_(0) {}
+
+ public:
+
+ /// initialize this generator
+ void initialize();
+
+ /// finalize this generator
+ void finalize() {}
+
+ /// generate an event, returning
+ /// the sign of the weight or zero
+ /// for an event below the evolution cutoff
+ double generate();
+
+ /// return the last sampled phase space point
+ const std::vector<double>& last_point() const { return last_point_; }
+
+ /// return the last evaluated function
+ double last_value() const { return last_value_; }
+
+ /// indicate that the last generated point has been rejected
+ void reject() {
+ last_cell_->info().reject();
+ }
+
+ public:
+
+ /// return true, if this generator has been initialized
+ bool initialized() const { return initialized_; }
+
+ /// return true, if at least one split has been performed
+ bool did_split() const { return did_split_; }
+
+ /// access the function
+ Function& function() { return *function_; }
+
+ /// set the function
+ void function(Function * f) { function_ = f; }
+
+ /// access the adaption_info object
+ adaption_info& sampling_parameters() { return adaption_info_; }
+
+ private:
+
+ /// check for and possibly split
+ /// the last selected cell
+ bool split();
+
+ /// get the projection of the density integrating over every
+ /// variable to be sampled, except the evolution variable for the
+ /// indicated parameter point. the k'th entry in
+ /// last_exponent_integrand_ is the value in the evolution
+ /// variable bin from evolution_splits_[k] to
+ /// evolution_splits_[k+1]
+ void get_exponent();
+
+ /// compensate
+ void compensate();
+
+ /// get all parameter points to build
+ /// all possible sub tree hashes
+ std::set<std::vector<double> > parameter_points();
+
+ /// get all parameter points to build
+ /// all possible sub tree hashes
+ void recursive_parameter_points(std::set<std::vector<double> >&,
+ std::vector<double>&,
+ size_t);
+
+ /// function to be sampled
+ Function * function_;
+
+ /// the number of events after which
+ /// a cell is checked for splits
+ unsigned long check_events_;
+
+ /// the adaption info object
+ adaption_info adaption_info_;
+
+ /// the root cell
+ binary_tree<cell> root_cell_;
+
+ /// the random number generator to be used
+ rnd_generator<Random> rnd_gen_;
+
+ /// wether a split has already been performed
+ bool did_split_;
+
+ /// wether this generator has been initialized
+ bool initialized_;
+
+ /// the position of the evolution variable
+ std::size_t evolution_variable_;
+
+ /// the cutoff on the evolution variable
+ double evolution_cutoff_;
+
+ /// flags of variables to be sampled
+ /// including the evolution variable
+ std::vector<bool> sample_variables_;
+
+ /// flags of variables to be sampled
+ /// excluding the evolution variable
+ std::vector<bool> sample_other_variables_;
+
+ /// the splits in any parameter done so far
+ /// (including the evolution variable)
+ std::map<std::size_t,std::vector<double> > parameter_splits_;
+
+ /// the last selected cell
+ binary_tree<cell>::iterator last_cell_;
+
+ /// the last sampled phasespace point
+ std::vector<double> last_point_;
+
+ /// the last function value
+ double last_value_;
+
+ /// the last parameter bin id
+ bit_container<parameter_hash_bits> last_parameter_bin_;
+
+ /// map parameter bin ids to exponent interpolations
+ std::map<bit_container<parameter_hash_bits>,linear_interpolator > exponents_;
+
+ /// the last exponent integrand
+ std::vector<double> last_exponent_integrand_;
+
+ /// the last exponent
+ std::map<bit_container<parameter_hash_bits>,linear_interpolator >::iterator last_exponent_;
+
+ /// wether or not we are compensating
+ bool compensating_;
+
+ /// the integral accessor to be used
+ integral_accessor integral_accessor_;
+
+ /// the missing events accessor to be used
+ parametric_missing_accessor missing_accessor_;
+
+ /// the parametric selector to be used
+ parametric_selector parametric_selector_;
+
+ /// the parametric selector to be used for parameter bins
+ parametric_selector exponent_selector_;
+
+ /// the parametric sampler to be used
+ parametric_sampling_selector<rnd_generator<Random> > parametric_sampler_;
+
+ /// the number of trials in the veto loo so far
+ unsigned long attempts_;
+
+ /// the number of accepted events so far
+ unsigned long accepts_;
+
+ /// number of splits done
+ unsigned long splits_;
+
+ };
+
+}
+
+#include "exponential_generator.icc"
+
+#endif // EXSAMPLE_exponential_generator_h_included
diff --git a/Exsample2/exsample/exponential_generator.icc b/Exsample2/exsample/exponential_generator.icc
new file mode 100644
--- /dev/null
+++ b/Exsample2/exsample/exponential_generator.icc
@@ -0,0 +1,323 @@
+// -*- C++ -*-
+//
+// exponential_generator.icc is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
+//
+// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
+//
+// ExSample is licenced under version 2 of the GPL, see COPYING for details.
+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
+//
+//
+namespace exsample {
+
+ template<class Function, class Random>
+ void exponential_generator<Function,Random>::initialize() {
+ adaption_info_.dimension = function_->dimension();
+ adaption_info_.lower_left = function_->support().first;
+ adaption_info_.upper_right = function_->support().second;
+ if (adaption_info_.adapt.empty())
+ adaption_info_.adapt = std::vector<bool>(adaption_info_.dimension,true);
+ evolution_variable_ = function_->evolution_variable();
+ evolution_cutoff_ = function_->evolution_cutoff();
+ sample_variables_ = function_->variable_flags();
+ sample_other_variables_ = sample_variables_;
+ sample_other_variables_[evolution_variable_] = false;
+ last_point_.resize(adaption_info_.dimension);
+ parametric_selector_ = parametric_selector(&last_point_,sample_other_variables_);
+ exponent_selector_ = parametric_selector(&last_point_,sample_variables_);
+ missing_accessor_ = parametric_missing_accessor(&last_parameter_bin_);
+ parametric_sampler_ = parametric_sampling_selector<rnd_generator<Random> >
+ (&last_point_,&last_parameter_bin_,sample_other_variables_,rnd_gen_);
+ if (initialized_) return;
+ splits_ = 0;
+ for ( std::size_t k = 0; k < adaption_info_.dimension; ++k ) {
+ if ( sample_other_variables_[k] )
+ continue;
+ parameter_splits_[k].push_back(adaption_info_.lower_left[k]);
+ parameter_splits_[k].push_back(adaption_info_.upper_right[k]);
+ }
+ root_cell_ =
+ binary_tree<cell>(cell(adaption_info_.lower_left,
+ adaption_info_.upper_right,
+ sample_other_variables_,
+ adaption_info_));
+ root_cell_.value().info().explore(rnd_gen_,adaption_info_,function_);
+ root_cell_.value().integral(root_cell_.value().info().volume() * root_cell_.value().info().overestimate());
+ last_exponent_integrand_.resize(1);
+ check_events_ = adaption_info_.presampling_points;
+ initialized_ = true;
+ }
+
+ template<class Function, class Random>
+ bool exponential_generator<Function,Random>::split () {
+ if (++splits_ > parameter_hash_bits/2)
+ return false;
+ if (adaption_info_.freeze_grid == accepts_)
+ return false;
+ if (compensating_)
+ return false;
+ if (!(*last_cell_).info().bad(adaption_info_)) return false;
+ bool dosplit = false;
+ std::pair<std::size_t,double> sp =
+ (*last_cell_).info().get_split(adaption_info_,dosplit);
+ if (!dosplit) return false;
+ if (!adaption_info_.adapt[sp.first]) return false;
+ last_cell_.node().split((*last_cell_).split(sp,rnd_gen_,function_,adaption_info_,sample_other_variables_));
+ if ( !sample_other_variables_[sp.first] ) {
+ if ( std::find(parameter_splits_[sp.first].begin(),parameter_splits_[sp.first].end(),sp.second)
+ == parameter_splits_[sp.first].end() ) {
+ parameter_splits_[sp.first].push_back(sp.second);
+ std::sort(parameter_splits_[sp.first].begin(),parameter_splits_[sp.first].end());
+ if ( sp.first == evolution_variable_ ) {
+ last_exponent_integrand_.push_back(0.);
+ }
+ }
+ }
+ did_split_ = true;
+ last_point_ = function_->parameter_point();
+ root_cell_.tree_accumulate(parametric_selector_,integral_accessor_,std::plus<double>());
+ exponents_.clear();
+ get_exponent();
+ return true;
+ }
+
+ template<class Function, class Random>
+ void exponential_generator<Function,Random>::get_exponent () {
+ last_parameter_bin_.reset();
+ root_cell_.subtree_hash (exponent_selector_,last_parameter_bin_);
+ last_exponent_ = exponents_.find(last_parameter_bin_);
+ if (last_exponent_ != exponents_.end())
+ return;
+ exponents_[last_parameter_bin_] = linear_interpolator();
+ last_exponent_ = exponents_.find(last_parameter_bin_);
+ double old_evo = last_point_[evolution_variable_];
+ std::vector<double>::iterator exp_it = last_exponent_integrand_.begin();
+ for (std::vector<double>::iterator esp = parameter_splits_[evolution_variable_].begin();
+ esp < prior(parameter_splits_[evolution_variable_].end()); ++esp, ++exp_it) {
+ last_point_[evolution_variable_] = (*esp + *next(esp))/2.;
+ *exp_it = root_cell_.accumulate(parametric_selector_,integral_accessor_,std::plus<double>());
+ }
+ exp_it = prior(last_exponent_integrand_.end());
+ double total = 0.;
+ for (std::vector<double>::iterator esp = prior(parameter_splits_[evolution_variable_].end());
+ esp > parameter_splits_[evolution_variable_].begin(); --esp, --exp_it) {
+ last_exponent_->second.set_interpolation(*esp,total);
+ total += (*exp_it) * ((*esp) - (*prior(esp)));
+ }
+ last_exponent_->second.set_interpolation(parameter_splits_[evolution_variable_].front(),total);
+ last_point_[evolution_variable_] = old_evo;
+ }
+
+ template<class Function, class Random>
+ std::set<std::vector<double> >
+ exponential_generator<Function,Random>::parameter_points() {
+ std::set<std::vector<double> > res;
+ std::vector<double> pt(adaption_info_.dimension,0.);
+ recursive_parameter_points(res,pt,0);
+ return res;
+ }
+
+ template<class Function, class Random>
+ void exponential_generator<Function,Random>::
+ recursive_parameter_points(std::set<std::vector<double> >& res,
+ std::vector<double>& pt,
+ size_t current) {
+ if ( current == adaption_info_.dimension ) {
+ res.insert(pt);
+ return;
+ }
+ if ( sample_variables_[current] ) {
+ recursive_parameter_points(res,pt,current+1);
+ return;
+ }
+ for ( std::vector<double>::const_iterator sp =
+ parameter_splits_[current].begin();
+ sp != prior(parameter_splits_[current].end()); ++sp ) {
+ pt[current] = (*sp + *next(sp))/2.;
+ recursive_parameter_points(res,pt,current+1);
+ }
+ }
+
+ template<class Function, class Random>
+ void exponential_generator<Function,Random>::compensate() {
+ if (!did_split_) {
+ exponents_.clear();
+ root_cell_.value().info().overestimate(last_value_,last_point_);
+ root_cell_.value().integral(root_cell_.value().info().volume() * root_cell_.value().info().overestimate());
+ last_point_ = function_->parameter_point();
+ get_exponent();
+ return;
+ }
+ std::vector<double> themaxpoint = last_point_;
+ std::set<std::vector<double> > id_points
+ = parameter_points();
+ for ( std::set<std::vector<double> >::const_iterator id =
+ id_points.begin(); id != id_points.end(); ++id ) {
+ last_point_ = *id;
+ get_exponent();
+ }
+ std::map<bit_container<parameter_hash_bits>,linear_interpolator >
+ old_exponents = exponents_;
+ double old_oe = last_cell_->info().overestimate();
+ last_cell_->info().overestimate(last_value_,themaxpoint);
+ last_cell_->integral(last_cell_->info().volume() * last_cell_->info().overestimate());
+ exponents_.clear();
+ for ( std::set<std::vector<double> >::const_iterator id =
+ id_points.begin(); id != id_points.end(); ++id ) {
+ last_point_ = *id;
+ get_exponent();
+ std::map<bit_container<parameter_hash_bits>,linear_interpolator >::iterator
+ old_exp = old_exponents.find(last_parameter_bin_);
+ std::map<bit_container<parameter_hash_bits>,linear_interpolator >::iterator
+ new_exp = exponents_.find(last_parameter_bin_);
+ assert(old_exp != old_exponents.end() && new_exp != exponents_.end());
+ double old_norm = 1. - std::exp(-(old_exp->second)(adaption_info_.lower_left[evolution_variable_]));
+ double new_norm = 1. - std::exp(-(new_exp->second)(adaption_info_.lower_left[evolution_variable_]));
+ for (binary_tree<cell>::iterator it = root_cell_.begin();
+ it != root_cell_.end(); ++it) {
+ if ( !it->info().contains_parameter(last_point_,sample_other_variables_) )
+ continue;
+ double old_int = 0.;
+ double new_int = 0.;
+ for ( std::vector<double>::const_iterator sp = parameter_splits_[evolution_variable_].begin();
+ sp != prior(parameter_splits_[evolution_variable_].end()); ++sp ) {
+ if ( *sp >= it->info().lower_left()[evolution_variable_] &&
+ *sp < it->info().upper_right()[evolution_variable_] ) {
+ double xl = *sp;
+ double xxl = *next(sp);
+ double old_al =
+ (old_exp->second.interpolation()[xxl] - old_exp->second.interpolation()[xl]) /
+ (xxl-xl);
+ double old_bl =
+ (xxl * old_exp->second.interpolation()[xl] -
+ xl * old_exp->second.interpolation()[xxl]) /
+ (xxl-xl);
+ double new_al =
+ (new_exp->second.interpolation()[xxl] - new_exp->second.interpolation()[xl]) /
+ (xxl-xl);
+ double new_bl =
+ (xxl * new_exp->second.interpolation()[xl] -
+ xl * new_exp->second.interpolation()[xxl]) /
+ (xxl-xl);
+ if ( std::abs(old_al) > std::numeric_limits<double>::epsilon() ) {
+ old_int += (exp(-(old_al*xl+old_bl)) - exp(-(old_al*xxl+old_bl)))/old_al;
+ } else {
+ old_int += (xxl-xl)*exp(-old_bl);
+ }
+ if ( std::abs(new_al) > std::numeric_limits<double>::epsilon() ) {
+ new_int += (exp(-(new_al*xl+new_bl)) - exp(-(new_al*xxl+new_bl)))/new_al;
+ } else {
+ new_int += (xxl-xl)*exp(-new_bl);
+ }
+ }
+ }
+ double scaling;
+ if (it != last_cell_) {
+ if (old_int > std::numeric_limits<double>::epsilon() &&
+ new_int > std::numeric_limits<double>::epsilon())
+ scaling = ((old_norm * new_int) /
+ (new_norm * old_int)) - 1.;
+ else
+ scaling = 0.;
+ } else {
+ if (old_int > std::numeric_limits<double>::epsilon() &&
+ new_int > std::numeric_limits<double>::epsilon())
+ scaling = ((last_value_ * old_norm * new_int) /
+ (old_oe * new_norm * old_int)) - 1.;
+ else
+ scaling = 0.;
+ }
+ it->info().parametric_missing(last_parameter_bin_,
+ it->info().parametric_missing(last_parameter_bin_) +
+ static_cast<int>(round(scaling * it->info().attempted())));
+ if (it->info().parametric_missing(last_parameter_bin_) != 0) {
+ compensating_ = true;
+ }
+ }
+ }
+ last_point_ = function_->parameter_point();
+ }
+
+ template<class Function, class Random>
+ double exponential_generator<Function,Random>::generate () {
+ if (compensating_) {
+ compensating_ = false;
+ for (binary_tree<cell>::iterator it = root_cell_.begin();
+ it != root_cell_.end(); ++it)
+ if (it->info().parametric_compensating()) {
+ compensating_ = true;
+ break;
+ }
+ parametric_sampler_.compensate(compensating_);
+ }
+ last_point_ = function_->parameter_point();
+ if (last_point_[evolution_variable_] < evolution_cutoff_) {
+ return 0.;
+ }
+ unsigned long n_hit_miss = 0;
+ unsigned long n_select = 0;
+ double minus_log_r;
+ root_cell_.tree_accumulate(parametric_selector_,integral_accessor_,std::plus<double>());
+ get_exponent();
+ while (true) {
+ n_select = 0;
+ minus_log_r = -std::log(rnd_gen_()) + last_exponent_->second(last_point_[evolution_variable_]);
+ if (!last_exponent_->second.invertible(minus_log_r)) {
+ return 0.;
+ }
+ try {
+ last_point_[evolution_variable_] = last_exponent_->second.unique_inverse(minus_log_r);
+ } catch (constant_interpolation& c) {
+ last_point_[evolution_variable_] = rnd_gen_(c.range.first,c.range.second);
+ }
+ assert(!std::isnan(last_point_[evolution_variable_]) &&
+ !std::isinf(last_point_[evolution_variable_]));
+ if (last_point_[evolution_variable_] < evolution_cutoff_) {
+ return 0.;
+ }
+ ++attempts_;
+ if (compensating_) {
+ root_cell_.tree_accumulate(missing_accessor_,std::plus<int>());
+ }
+ if (parameter_splits_[evolution_variable_].size() > 2)
+ root_cell_.tree_accumulate(parametric_selector_,integral_accessor_,std::plus<double>());
+ if (did_split_)
+ while ((last_cell_ = root_cell_.select(parametric_sampler_)) == root_cell_.end()) {
+ root_cell_.tree_accumulate(missing_accessor_,std::plus<int>());
+ if(++n_select > adaption_info_.maxtry)
+ throw selection_maxtry();
+ }
+ else
+ last_cell_ = root_cell_.begin();
+ last_cell_->info().select(rnd_gen_,last_point_,sample_other_variables_);
+ last_value_ = function_->evaluate(last_point_);
+ assert(last_value_ >= 0.);
+ last_cell_->info().selected(last_point_,last_value_,adaption_info_);
+ if (last_value_ > last_cell_->info().overestimate()) {
+ if ( std::abs(last_value_)/last_cell_->info().overestimate() > 2. ) {
+ last_value_ =
+ last_value_*(1.+exp(2.*(2.-std::abs(last_value_)/last_cell_->info().overestimate())));
+ }
+ compensate();
+ throw exponential_regenerate();
+ }
+ if (last_cell_->info().attempted() % check_events_ == 0) {
+ if (split()) {
+ throw exponential_regenerate();
+ }
+ }
+ if (last_value_/last_cell_->info().overestimate() > rnd_gen_())
+ break;
+ if(++n_hit_miss > adaption_info_.maxtry)
+ throw hit_and_miss_maxtry();
+ }
+ if (last_value_ == 0.)
+ return 0.;
+ ++accepts_;
+ ++check_events_;
+ last_cell_->info().accept();
+ return 1.;
+ }
+
+
+}
diff --git a/Exsample2/exsample/generator.h b/Exsample2/exsample/generator.h
new file mode 100644
--- /dev/null
+++ b/Exsample2/exsample/generator.h
@@ -0,0 +1,181 @@
+// -*- C++ -*-
+//
+// generator.h is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
+//
+// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
+//
+// ExSample is licenced under version 2 of the GPL, see COPYING for details.
+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
+//
+//
+#ifndef EXSAMPLE_generator_h_included
+#define EXSAMPLE_generator_h_included
+
+#include "cell.h"
+#include "selectors.h"
+#include "statistics.h"
+#include "binary_tree.h"
+
+namespace exsample {
+
+ /// \brief Exception thrown, if the generator has just changed its
+ /// state. The attempt of generating an event should be repeated.
+ struct generator_update{};
+
+ /// \brief A generator for plain sampling and integrating
+ template<class Function, class Random>
+ class generator {
+
+ public:
+
+ /// default constructor
+ generator()
+ : function_(0), statistics_(), check_events_(0),
+ adaption_info_(), root_cell_(),
+ rnd_gen_(), did_split_(false), initialized_(false),
+ last_cell_(), last_point_(), last_value_(0.),
+ compensating_(false) {}
+
+ public:
+
+ /// initialize this generator
+ template<class SlaveStatistics>
+ void initialize(SlaveStatistics&);
+
+ /// generate an event, returning
+ /// the sign of the weight
+ template<class SlaveStatistics>
+ double generate(SlaveStatistics&);
+
+ /// return the last sampled phase space point
+ const std::vector<double>& last_point() const { return last_point_; }
+
+ /// indicate that the last generated point has been rejected
+ void reject() {
+ statistics_.reject(last_value_);
+ last_cell_->info().reject();
+ }
+
+ /// finalize this generator
+ void finalize() {
+ statistics_.reset();
+ }
+
+ public:
+
+ /// return true, if this generator has been initialized
+ bool initialized() const { return initialized_; }
+
+ /// return true, if at least one split has been performed
+ bool did_split() const { return did_split_; }
+
+ /// access the function
+ Function& function() { return *function_; }
+
+ /// set the function
+ void function(Function * f) { function_ = f; }
+
+ /// return the statistics object
+ const statistics& stats() const { return statistics_; }
+
+ /// return the sampled volume
+ double volume() const {
+ return exsample::volume(adaption_info_.lower_left, adaption_info_.upper_right);
+ }
+
+ /// return the integral
+ double integral() const {
+ return volume() * statistics_.average_weight();
+ }
+
+ /// return the error on the integral
+ double integral_uncertainty() const {
+ return volume() * std::sqrt(statistics_.average_weight_variance());
+ }
+
+ /// return the integral
+ double current_integral() const {
+ return volume() * statistics_.current().first;
+ }
+
+ /// return the error on the integral
+ double current_integral_uncertainty() const {
+ return volume() * std::sqrt(statistics_.current().second);
+ }
+
+ /// return the variance of the integral estimate
+ double integral_variance() const {
+ return sqr(volume()) * statistics_.average_weight_variance();
+ }
+
+ /// access the adaption_info object
+ adaption_info& sampling_parameters() { return adaption_info_; }
+
+ /// return true, if still compensating
+ bool compensating() const { return compensating_; }
+
+ public:
+
+ /// put to ostream
+ template<class OStream>
+ void put(OStream& os) const;
+
+ /// get from istream
+ template<class IStream>
+ void get(IStream& is);
+
+ private:
+
+ /// check for and possibly split
+ /// the last selected cell
+ bool split();
+
+ /// compensate the last selected cell indicating the
+ /// value and position of the new overestimate
+ void compensate();
+
+ /// function to be sampled
+ Function * function_;
+
+ /// the global statistics object
+ statistics statistics_;
+
+ /// the number of events after which
+ /// a cell is checked for splits
+ unsigned long check_events_;
+
+ /// the adaption info object
+ adaption_info adaption_info_;
+
+ /// the root cell
+ binary_tree<cell> root_cell_;
+
+ /// the random number generator to be used
+ rnd_generator<Random> rnd_gen_;
+
+ /// wether a split has already been performed
+ bool did_split_;
+
+ /// wether this generator has been initialized
+ bool initialized_;
+
+ /// the last selected cell
+ binary_tree<cell>::iterator last_cell_;
+
+ /// the last sampled phasespace point
+ std::vector<double> last_point_;
+
+ /// the last function value
+ double last_value_;
+
+ /// wether or not we are compensating
+ bool compensating_;
+
+ };
+
+}
+
+#include "generator.icc"
+
+#endif // EXSAMPLE_generator_h_included
+
diff --git a/Exsample2/exsample/generator.icc b/Exsample2/exsample/generator.icc
new file mode 100644
--- /dev/null
+++ b/Exsample2/exsample/generator.icc
@@ -0,0 +1,178 @@
+// -*- C++ -*-
+//
+// generator.icc is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
+//
+// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
+//
+// ExSample is licenced under version 2 of the GPL, see COPYING for details.
+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
+//
+//
+namespace exsample {
+
+ template<class Function, class Random>
+ template<class SlaveStatistics>
+ void generator<Function,Random>::initialize(SlaveStatistics& opt) {
+ adaption_info_.dimension = function_->dimension();
+ adaption_info_.lower_left = function_->support().first;
+ adaption_info_.upper_right = function_->support().second;
+ if (adaption_info_.adapt.empty())
+ adaption_info_.adapt = std::vector<bool>(adaption_info_.dimension,true);
+ last_point_.resize(adaption_info_.dimension);
+ if (initialized_) return;
+ root_cell_ =
+ binary_tree<cell>(cell(adaption_info_.lower_left,
+ adaption_info_.upper_right,
+ adaption_info_));
+ root_cell_.value().info().explore(rnd_gen_,adaption_info_,function_,&statistics_,opt);
+ root_cell_.value().integral(root_cell_.value().info().volume() * root_cell_.value().info().overestimate());
+ statistics_.reset();
+ check_events_ = adaption_info_.presampling_points;
+ initialized_ = true;
+ }
+
+ template<class Function, class Random>
+ bool generator<Function,Random>::split() {
+ if (adaption_info_.freeze_grid < statistics_.accepted() &&
+ adaption_info_.freeze_grid != 0)
+ return false;
+ if (compensating_) return false;
+ if (!(*last_cell_).info().bad(adaption_info_)) return false;
+ bool dosplit = false;
+ std::pair<std::size_t,double> sp =
+ (*last_cell_).info().get_split(adaption_info_,dosplit);
+ if (!dosplit || !adaption_info_.adapt[sp.first]) {
+ return false;
+ }
+ std::pair<binary_tree<cell>::iterator,binary_tree<cell>::iterator> spit =
+ last_cell_.node().split((*last_cell_).split(sp,rnd_gen_,function_,adaption_info_));
+ integral_accessor iacc;
+ root_cell_.tree_accumulate(iacc,std::plus<double>());
+ did_split_ = true;
+ statistics_.reset();
+ return true;
+ }
+
+ template<class Function, class Random>
+ void generator<Function,Random>::compensate() {
+ if (!did_split_) {
+ root_cell_.value().info().overestimate(std::abs(last_value_),last_point_);
+ root_cell_.value().integral(root_cell_.value().info().volume() * root_cell_.value().info().overestimate());
+ return;
+ }
+ double old_norm = root_cell_.value().integral();
+ double new_norm = old_norm - (*last_cell_).integral() + std::abs(last_value_) * (*last_cell_).info().volume();
+ compensating_ = false;
+ last_cell_->missing_events() +=
+ static_cast<int>(round(((std::abs(last_value_) * old_norm)/(last_cell_->info().overestimate() * new_norm) - 1.) *
+ (last_cell_->info().attempted())));
+ if (last_cell_->missing_events() != 0)
+ compensating_ = true;
+ last_cell_->info().overestimate(std::abs(last_value_),last_point_);
+ last_cell_->integral(last_cell_->info().volume() * last_cell_->info().overestimate());
+ for (binary_tree<cell>::iterator it = root_cell_.begin();
+ it != root_cell_.end(); ++it)
+ if (it != last_cell_) {
+ it->missing_events() += static_cast<int>(round((old_norm/new_norm - 1.) * (it->info().attempted())));
+ if (it->missing_events() != 0) {
+ compensating_ = true;
+ }
+ }
+ integral_accessor iacc;
+ root_cell_.tree_accumulate(iacc,std::plus<double>());
+ statistics_.reset();
+ root_cell_.tree_accumulate(missing_accessor(),std::plus<int>());
+ }
+
+ template<class Function, class Random>
+ template<class SlaveStatistics>
+ double generator<Function,Random>::generate(SlaveStatistics& opt) {
+ unsigned long n_hit_miss = 0;
+ unsigned long n_select = 0;
+ sampling_selector<rnd_generator<Random> > sampler (rnd_gen_,compensating_);
+ missing_accessor macc;
+ if (compensating_) {
+ compensating_ = false;
+ for (binary_tree<cell>::iterator it = root_cell_.begin();
+ it != root_cell_.end(); ++it)
+ if (it->missing_events() != 0) {
+ compensating_ = true;
+ break;
+ }
+ }
+ while (true) {
+ sampler.compensate = compensating_;
+ n_select = 0;
+ if (did_split_)
+ while ((last_cell_ = root_cell_.select(sampler)) == root_cell_.end()) {
+ root_cell_.tree_accumulate(macc,std::plus<int>());
+ if(++n_select > adaption_info_.maxtry)
+ throw selection_maxtry();
+ }
+ else
+ last_cell_ = root_cell_.begin();
+ last_cell_->info().select(rnd_gen_,last_point_);
+ last_value_ = function_->evaluate(last_point_);
+ last_cell_->info().selected(last_point_,std::abs(last_value_),adaption_info_);
+ if (std::abs(last_value_) > last_cell_->info().overestimate()) {
+ if ( std::abs(last_value_)/last_cell_->info().overestimate() > 2. ) {
+ last_value_ =
+ last_value_*(1.+exp(2.*(2.-std::abs(last_value_)/last_cell_->info().overestimate())));
+ }
+ compensate();
+ n_hit_miss = 0;
+ continue;
+ }
+ if (last_cell_->info().attempted() % check_events_ == 0) {
+ if (split()) {
+ throw generator_update();
+ }
+ }
+ if (did_split_) {
+ statistics_.select(last_value_ * root_cell_.value().integral() /
+ last_cell_->info().overestimate(), !compensating_);
+ opt.select(last_value_ * root_cell_.value().integral() /
+ last_cell_->info().overestimate(), !compensating_);
+ } else {
+ statistics_.select(last_value_, !compensating_);
+ opt.select(last_value_, !compensating_);
+ }
+ if (std::abs(last_value_)/last_cell_->info().overestimate() > rnd_gen_())
+ break;
+ if(++n_hit_miss > adaption_info_.maxtry)
+ throw hit_and_miss_maxtry();
+ }
+ last_cell_->info().accept();
+ if (did_split_)
+ last_value_ *= root_cell_.value().integral() / last_cell_->info().overestimate();
+ statistics_.accept(last_value_);
+ ++check_events_;
+ if (last_value_ < 0.)
+ return -1.;
+ return 1.;
+ }
+
+ template<class Function, class Random>
+ template<class OStream>
+ void generator<Function,Random>::put (OStream& os) const {
+ adaption_info_.put(os);
+ root_cell_.put(os);
+ statistics_.put(os);
+ os << check_events_;
+ ostream_traits<OStream>::separator(os);
+ os << did_split_;
+ ostream_traits<OStream>::separator(os);
+ os << initialized_;
+ ostream_traits<OStream>::separator(os);
+ }
+
+ template<class Function, class Random>
+ template<class IStream>
+ void generator<Function,Random>::get (IStream& is) {
+ adaption_info_.get(is);
+ root_cell_.get(is);
+ statistics_.get(is);
+ is >> check_events_ >> did_split_ >> initialized_;
+ }
+
+}
diff --git a/Exsample2/exsample/linear_interpolator.h b/Exsample2/exsample/linear_interpolator.h
new file mode 100644
--- /dev/null
+++ b/Exsample2/exsample/linear_interpolator.h
@@ -0,0 +1,123 @@
+// -*- C++ -*-
+//
+// linear_interpolator.h is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
+//
+// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
+//
+// ExSample is licenced under version 2 of the GPL, see COPYING for details.
+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
+//
+//
+#ifndef EXSAMPLE_linear_interpolator_h_included
+#define EXSAMPLE_linear_interpolator_h_included
+
+#include "utility.h"
+
+namespace exsample {
+
+ /// \brief Exception thrown, if inversion of the interpolation has
+ /// no solution
+ struct inversion_has_no_solution { };
+
+ /// \brief Exception thrown, if a constant piece of the
+ /// interpolation has been hit
+ struct constant_interpolation {
+
+ /// standard constructor
+ constant_interpolation(double xlow,
+ double xhigh,
+ double h)
+ : range(xlow,xhigh), value(h) {}
+
+ /// The range in which the interpolation is constant
+ std::pair<double,double> range;
+
+ /// The value of the interpolation in the range
+ double value;
+
+ };
+
+ /// \brief A linear interpolation allowing for inversion of the
+ /// linear interpolation
+ class linear_interpolator {
+
+ public:
+
+ /// default constructor
+ linear_interpolator();
+
+ /// construct a linear interpolator given the
+ /// map of interpolation points and values
+ explicit linear_interpolator(const std::map<double,double>& points);
+
+ public:
+
+ /// return the interpolations value at the given point
+ double operator()(double x) const;
+
+ /// return the range of the interpolation
+ std::pair<double,double> range() const { return range_; }
+
+ public:
+
+ /// return true, if an inverse exists for
+ /// the given value
+ bool invertible(double f) const {
+ return (range_.first <= f &&
+ range_.second >= f);
+ }
+
+ /// return the inverse, assuming the inverse
+ /// is unique
+ double unique_inverse(double f) const;
+
+ /// access the interpolation map
+ std::map<double,double>& interpolation() {
+ return interpolation_;
+ }
+
+ /// set the value at the given point
+ void set_interpolation(double point, double value);
+
+ /// reset after interpolation points have been changed
+ void reset();
+
+ public:
+
+ /// put to ostream
+ template<class OStream>
+ void put(OStream& os) const;
+
+ /// get from istream
+ template<class IStream>
+ void get(IStream& is);
+
+ private:
+
+ /// map points to values
+ std::map<double, double> interpolation_;
+
+ /// the range of the interpolation
+ std::pair<double,double> range_;
+
+ };
+
+ /// ostream operator for linear_interpolator objects
+ template<class OStream>
+ inline OStream& operator<<(OStream& os, const linear_interpolator& ip) {
+ ip.put(os);
+ return os;
+ }
+
+ /// istream operator for linear_interpolator objects
+ template<class IStream>
+ inline IStream& operator>>(IStream& is, linear_interpolator& ip) {
+ ip.get(is);
+ return is;
+ }
+
+}
+
+#include "linear_interpolator.icc"
+
+#endif // EXSAMPLE_linear_interpolator_h_included
diff --git a/Exsample2/exsample/linear_interpolator.icc b/Exsample2/exsample/linear_interpolator.icc
new file mode 100644
--- /dev/null
+++ b/Exsample2/exsample/linear_interpolator.icc
@@ -0,0 +1,111 @@
+// -*- C++ -*-
+//
+// linear_interpolator.icc is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
+//
+// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
+//
+// ExSample is licenced under version 2 of the GPL, see COPYING for details.
+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
+//
+//
+namespace exsample {
+
+ inline linear_interpolator::linear_interpolator()
+ : interpolation_(), range_() {}
+
+ inline linear_interpolator::linear_interpolator(const std::map<double,double>& points)
+ : interpolation_(points) {
+ reset();
+ }
+
+
+ inline void linear_interpolator::set_interpolation(double point, double value) {
+ interpolation_[point] = value;
+ if (value > range_.second)
+ range_.second = value;
+ if(value < range_.first)
+ range_.first = value;
+ }
+
+
+ inline void linear_interpolator::reset() {
+ range_.first = interpolation_.begin()->second;
+ range_.second = interpolation_.begin()->second;
+ for (std::map<double,double>::const_iterator c = interpolation_.begin();
+ c != interpolation_.end(); ++c) {
+ if (c->second < range_.first)
+ range_.first = c->second;
+ if (c->second > range_.second)
+ range_.second = c->second;
+ }
+ }
+
+
+ inline double linear_interpolator::operator()(double x) const {
+ std::map<double, double>::const_iterator upper =
+ interpolation_.upper_bound(x);
+ if (upper == interpolation_.end()) {
+ upper = interpolation_.upper_bound(x-1e-10);
+ }
+ if (upper == interpolation_.end()) {
+ upper = interpolation_.upper_bound(x+1e-10);
+ }
+ assert(upper != interpolation_.begin() &&
+ upper != interpolation_.end());
+ return ((upper->second-prior(upper)->second)*x +
+ prior(upper)->second*upper->first - upper->second*prior(upper)->first)/
+ (upper->first - prior(upper)->first);
+ }
+
+
+ inline double linear_interpolator::unique_inverse(double f) const {
+ if(!invertible(f)) throw inversion_has_no_solution();
+ std::map<double, double>::const_iterator lower = interpolation_.begin();
+ bool gotone = false;
+ for (; lower != --interpolation_.end(); ++lower)
+ if ((lower->second >= f && next(lower)->second <= f) ||
+ (lower->second <= f && next(lower)->second >= f)) {
+ gotone = true;
+ break;
+ }
+ if(!gotone) throw inversion_has_no_solution();
+ if (lower->second == next(lower)->second) {
+ throw constant_interpolation(lower->first,next(lower)->first,lower->second);
+ }
+ double xdiff = next(lower)->first - lower->first;
+ double wdiff = next(lower)->second - lower->second;
+ double woffset = lower->second * next(lower)->first - next(lower)->second * lower->first;
+ return (xdiff/wdiff)*(f - woffset/xdiff);
+ }
+
+
+ template<class OStream>
+ void linear_interpolator::put(OStream& os) const {
+ os << interpolation_.size();
+ ostream_traits<OStream>::separator(os);
+ for (std::map<double, double>::const_iterator p
+ = interpolation_.begin(); p != interpolation_.end(); ++p) {
+ os << p->first;
+ ostream_traits<OStream>::separator(os);
+ os << p->second;
+ ostream_traits<OStream>::separator(os);
+ }
+ os << range_.first;
+ ostream_traits<OStream>::separator(os);
+ os << range_.second;
+ ostream_traits<OStream>::separator(os);
+ }
+
+ template<class IStream>
+ void linear_interpolator::get(IStream& is) {
+ std::size_t size;
+ is >> size;
+ std::pair<double,double> point;
+ for (std::size_t k = 0; k < size; ++k) {
+ is >> point.first >> point.second;
+ interpolation_.insert(point);
+ }
+ is >> range_.first >> range_.second;
+ }
+
+}
diff --git a/Exsample2/exsample/selectors.h b/Exsample2/exsample/selectors.h
new file mode 100644
--- /dev/null
+++ b/Exsample2/exsample/selectors.h
@@ -0,0 +1,201 @@
+// -*- C++ -*-
+//
+// selectors.h is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
+//
+// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
+//
+// ExSample is licenced under version 2 of the GPL, see COPYING for details.
+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
+//
+//
+#ifndef EXSAMPLE_selectors_h_included
+#define EXSAMPLE_selectors_h_included
+
+#include "cell.h"
+
+namespace exsample {
+
+ /// \brief flat sampling selector
+ template<class Random>
+ struct sampling_selector {
+
+ /// the default constructor
+ sampling_selector()
+ : rnd_gen(), compensate(false) {}
+
+ /// the standard constructor
+ explicit sampling_selector(const Random& r, bool comp = true)
+ : rnd_gen(r), compensate(comp) {}
+
+ /// return which of the children cells should be considered
+ std::pair<bool,bool> use(cell& parent,
+ const cell& first_child,
+ const cell& second_child) const;
+
+ /// return true, if the leaf cell should be considered
+ bool use(cell& leaf) const;
+
+ /// The random number generator to be used
+ Random rnd_gen;
+
+ /// Whether or not compensation is needed
+ bool compensate;
+
+ };
+
+ /// \brief selector selecting only bins
+ /// which contain the given parameter point
+ class parametric_selector {
+
+ public:
+
+ /// the default constructor
+ parametric_selector ()
+ : point_(), sampled_variables_() {}
+
+ /// construct from reference to point and flags to sample variables
+ parametric_selector(std::vector<double> * point,
+ const std::vector<bool>& sample)
+ : point_(point), sampled_variables_(sample) {}
+
+ public:
+
+ /// return which of the children cells should be considered
+ std::pair<bool,bool> use(const cell& parent,
+ const cell&,
+ const cell&) const;
+
+ /// return true, if the leaf cell should be considered
+ bool use(const cell&) const { return true; }
+
+ private:
+
+ /// the point chosen
+ std::vector<double> * point_;
+
+ /// flags for variables to be sampled
+ std::vector<bool> sampled_variables_;
+
+ };
+
+ /// \brief sampling selector selecting only bins
+ /// which contain the given parameter point
+ template<class Random>
+ class parametric_sampling_selector {
+
+ public:
+
+ /// the default constructor
+ parametric_sampling_selector()
+ : point_(), bin_id_(),
+ sampled_variables_(), rnd_gen_(),
+ compensate_(false) {}
+
+ /// construct from reference to point, subtree hash, flags of
+ /// variables to be sampled, and random number generator
+ parametric_sampling_selector(std::vector<double> * p,
+ bit_container<parameter_hash_bits> * bin_id,
+ const std::vector<bool>& sample,
+ const Random& rnd_gen)
+ : point_(p), bin_id_(bin_id),
+ sampled_variables_(sample), rnd_gen_(rnd_gen),
+ compensate_(false) {}
+
+ public:
+
+ /// return which of the children cells should be considered
+ std::pair<bool,bool> use(cell& parent,
+ const cell& first_child,
+ const cell& second_child) const;
+
+ /// return true, if the leaf cell should be considered
+ bool use(cell& leaf) const;
+
+ /// indicate that compensation is to take place
+ void compensate (bool doit = true) { compensate_ = doit; }
+
+ private:
+
+ /// the point chosen
+ std::vector<double> * point_;
+
+ /// the corresponding bin id
+ bit_container<parameter_hash_bits> * bin_id_;
+
+ /// flags for variables to be sampled
+ std::vector<bool> sampled_variables_;
+
+ /// the random number generator
+ Random rnd_gen_;
+
+ /// wether or not compensation is needed
+ bool compensate_;
+
+ };
+
+ /// \brief accessor returning the integral of a cell
+ struct integral_accessor {
+
+ /// update and return the value
+ double& set(cell& node) const {
+ return node.integral();
+ }
+
+ /// get the value
+ double get(const cell& node, bool) const {
+ return node.integral();
+ }
+
+ };
+
+ /// \brief accessor returning the number of missing events
+ struct missing_accessor {
+
+ /// update and return the value
+ int& set(cell& node) const {
+ return node.missing_events();
+ }
+
+ /// get the value
+ int get(const cell& node, bool) const {
+ return node.missing_events();
+ }
+
+ };
+
+ /// \brief accessor returning the number of missing events
+ /// for given parameter bin id
+ struct parametric_missing_accessor {
+
+ /// the default constructor
+ parametric_missing_accessor ()
+ : id_() {}
+
+ /// construct from subtree hash to consider
+ explicit parametric_missing_accessor (bit_container<parameter_hash_bits>* id)
+ : id_ (id) {}
+
+ /// update and return the value
+ int& set(cell& node) const {
+ return node.missing_events();
+ }
+
+ /// get the value
+ int get(const cell& node, bool isleaf) const {
+ if (isleaf)
+ return node.info().parametric_missing(*id_);
+ return node.missing_events();
+ }
+
+ private:
+
+ /// the subtree hash to consider
+ bit_container<parameter_hash_bits>* id_;
+
+ };
+
+}
+
+#include "selectors.icc"
+
+#endif // EXSAMPLE_selectors_h_included
diff --git a/Exsample2/exsample/selectors.icc b/Exsample2/exsample/selectors.icc
new file mode 100644
--- /dev/null
+++ b/Exsample2/exsample/selectors.icc
@@ -0,0 +1,133 @@
+// -*- C++ -*-
+//
+// selectors.icc is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
+//
+// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
+//
+// ExSample is licenced under version 2 of the GPL, see COPYING for details.
+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
+//
+//
+namespace exsample {
+
+ template<class Random>
+ std::pair<bool,bool> sampling_selector<Random>::use(cell& parent,
+ const cell& first_child,
+ const cell& second_child) const {
+ std::pair<bool,bool> selected (false,false);
+ if (compensate) {
+ if (first_child.missing_events() > 0 && second_child.missing_events() <= 0) {
+ selected.first = true;
+ --parent.missing_events();
+ return selected;
+ }
+ if (first_child.missing_events() <= 0 && second_child.missing_events() > 0) {
+ selected.second = true;
+ --parent.missing_events();
+ return selected;
+ }
+ if (first_child.missing_events() > 0 && second_child.missing_events() > 0) {
+ if (first_child.integral()/parent.integral() > rnd_gen())
+ selected.first = true;
+ else
+ selected.second = true;
+ --parent.missing_events();
+ return selected;
+ }
+ }
+ if (first_child.integral()/parent.integral() > rnd_gen())
+ selected.first = true;
+ else
+ selected.second = true;
+ return selected;
+ }
+
+ template<class Random>
+ bool sampling_selector<Random>::use(cell& leaf) const {
+ if (compensate) {
+ if (leaf.missing_events() < 0) {
+ ++leaf.missing_events();
+ return false;
+ }
+ if (leaf.missing_events() > 0)
+ --leaf.missing_events();
+ }
+ return true;
+ }
+
+ inline std::pair<bool,bool> parametric_selector::use(const cell& parent,
+ const cell&,
+ const cell&) const {
+ std::pair<bool,bool> match_para (true,true);
+ if (!sampled_variables_[parent.split_dimension()]) {
+ match_para.first = (parent.split_point() > (*point_)[parent.split_dimension()]);
+ match_para.second = (parent.split_point() <= (*point_)[parent.split_dimension()]);
+ }
+ return match_para;
+ }
+
+ template<class Random>
+ std::pair<bool,bool> parametric_sampling_selector<Random>::use(cell& parent,
+ const cell& first_child,
+ const cell& second_child) const {
+ std::pair<bool,bool> match_para (true,true);
+ std::pair<bool,bool> selected (false,false);
+ if (!sampled_variables_[parent.split_dimension()]) {
+ match_para.first = (parent.split_point() > (*point_)[parent.split_dimension()]);
+ match_para.second = (parent.split_point() <= (*point_)[parent.split_dimension()]);
+ }
+ if (match_para.first && match_para.second) {
+ if (compensate_) {
+ if (first_child.missing_events() > 0 && second_child.missing_events() <= 0) {
+ selected.first = true;
+ --parent.missing_events();
+ return selected;
+ }
+ if (first_child.missing_events() <= 0 && second_child.missing_events() > 0) {
+ selected.second = true;
+ --parent.missing_events();
+ return selected;
+ }
+
+ if (first_child.missing_events() > 0 && second_child.missing_events() > 0) {
+ if (first_child.integral()/parent.integral() > rnd_gen_())
+ selected.first = true;
+ else
+ selected.second = true;
+ --parent.missing_events();
+ return selected;
+ }
+ }
+ if (first_child.integral()/parent.integral() > rnd_gen_())
+ selected.first = true;
+ else
+ selected.second = true;
+ return selected;
+ }
+ assert(match_para.first || match_para.second);
+ if (compensate_)
+ if ((match_para.first && first_child.missing_events() > 0) ||
+ (match_para.second && second_child.missing_events() > 0))
+ --parent.missing_events();
+ return match_para;
+ }
+
+ template<class Random>
+ bool parametric_sampling_selector<Random>::use (cell& leaf) const {
+ if (compensate_) {
+ long pmissing = leaf.info().parametric_missing(*bin_id_);
+ assert(leaf.missing_events() == pmissing);
+ if (pmissing < 0) {
+ leaf.info().increase_parametric_missing(*bin_id_);
+ ++leaf.missing_events();
+ return false;
+ }
+ if (pmissing > 0) {
+ leaf.info().decrease_parametric_missing(*bin_id_);
+ --leaf.missing_events();
+ }
+ }
+ return true;
+ }
+
+}
diff --git a/Exsample2/exsample/statistics.h b/Exsample2/exsample/statistics.h
new file mode 100644
--- /dev/null
+++ b/Exsample2/exsample/statistics.h
@@ -0,0 +1,151 @@
+// -*- C++ -*-
+//
+// statistics.h is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
+//
+// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
+//
+// ExSample is licenced under version 2 of the GPL, see COPYING for details.
+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
+//
+//
+#ifndef EXSAMPLE_statistics_h_included
+#define EXSAMPLE_statistics_h_included
+
+#include "utility.h"
+
+namespace exsample {
+
+ /// \brief statistics is a helper class for keeping track of event
+ /// generation statistics.
+ class statistics {
+
+ public:
+
+ /// default constructor
+ statistics();
+
+ /// update the statistics for a weight encountered during
+ /// presampling
+ void presampled(double weight);
+
+ /// indicate that a weight has been selected; optionally preven
+ /// the weight from entering the caluclation of the integral
+ void select(double weight,
+ bool calculate_integral = true);
+
+ /// indicate that a point has been accepted
+ void accept(double weight) {
+ ++accepted_;
+ if (weight < 0) ++accepted_negative_;
+ }
+
+ /// reject a prviously accepted event
+ void reject(double weight) {
+ --accepted_;
+ if (weight < 0) --accepted_negative_;
+ }
+
+ /// reset the statistics object
+ void reset();
+
+ /// return the integral's estimate and its uncertainty at the
+ /// currently accumulated statistics
+ std::pair<double,double> current() const;
+
+ /// the average weight
+ double average_weight() const {
+ return (n_iterations_ == 0 ? 0. : average_weight_/n_iterations_);
+ }
+
+ /// the average absolute weight
+ double average_abs_weight() const {
+ return (n_iterations_ == 0 ? 0. : average_abs_weight_/n_iterations_);
+ }
+
+ /// the variance of the weight
+ double average_weight_variance() const {
+ return (n_iterations_ == 0 ? 0. : average_weight_variance_/n_iterations_);
+ }
+
+ /// the number of points in this iteration
+ unsigned long iteration_points() const { return iteration_points_; }
+
+ /// the number of iterations
+ unsigned long n_iterations() const { return n_iterations_; }
+
+ /// the total number of attempted in this bin
+ unsigned long attempted() const { return attempted_; }
+
+ /// the total number of finally accepted events in this bin
+ unsigned long accepted() const { return accepted_; }
+
+ /// the total number of acceptet events with negative weights
+ unsigned long accepted_negative() const { return accepted_negative_; }
+
+ /// the sum of weights
+ double sum_weights() const { return sum_weights_; }
+
+ /// the sum of absolute values of the weights
+ double sum_abs_weights() const { return sum_abs_weights_; }
+
+ /// the sum of weights squared
+ double sum_weights_squared() const { return sum_weights_squared_; }
+
+ /// the maximum weight
+ double max_weight() const { return max_weight_; }
+
+ public:
+
+ /// put ostream
+ template<class OStream>
+ void put(OStream& os) const;
+
+ /// get from istream
+ template<class IStream>
+ void get(IStream& is);
+
+ private:
+
+ /// the average weight
+ double average_weight_;
+
+ /// the average absolute weight
+ double average_abs_weight_;
+
+ /// the variance of the weight
+ double average_weight_variance_;
+
+ /// the number of points in this iteration
+ unsigned long iteration_points_;
+
+ /// the total number of attempted in this bin
+ unsigned long attempted_;
+
+ /// the total number of finally accepted events in this bin
+ unsigned long accepted_;
+
+ /// the total number of acceptet events with negative weights
+ unsigned long accepted_negative_;
+
+ /// the sum of weights
+ double sum_weights_;
+
+ /// the sum of absolute values of the weights
+ double sum_abs_weights_;
+
+ /// the sum of weights squared
+ double sum_weights_squared_;
+
+ /// the maximum weight
+ double max_weight_;
+
+ /// the number of iterations used to calculate the integral
+ unsigned long n_iterations_;
+
+ };
+
+}
+
+#include "statistics.icc"
+
+#endif // EXSAMPLE_statistics_h_included
diff --git a/Exsample2/exsample/statistics.icc b/Exsample2/exsample/statistics.icc
new file mode 100644
--- /dev/null
+++ b/Exsample2/exsample/statistics.icc
@@ -0,0 +1,108 @@
+// -*- C++ -*-
+//
+// statistics.icc is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
+//
+// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
+//
+// ExSample is licenced under version 2 of the GPL, see COPYING for details.
+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
+//
+//
+namespace exsample {
+
+ inline statistics::statistics()
+ : average_weight_(0.), average_abs_weight_(0.), average_weight_variance_(-1.),
+ iteration_points_(0),
+ attempted_(0), accepted_(0), accepted_negative_(0),
+ sum_weights_(0.), sum_abs_weights_(0.), sum_weights_squared_(0.),
+ max_weight_(0.), n_iterations_(0) {
+ }
+
+ inline void statistics::presampled(double weight) {
+ ++iteration_points_;
+ sum_weights_ += weight;
+ sum_abs_weights_ += std::abs(weight);
+ sum_weights_squared_ += sqr(weight);
+ max_weight_ = std::max(max_weight_,std::abs(weight));
+ }
+
+ inline void statistics::select(double weight,
+ bool calculate_integral) {
+ ++attempted_;
+ if (calculate_integral) {
+ ++iteration_points_;
+ sum_weights_ += weight;
+ sum_abs_weights_ += std::abs(weight);
+ sum_weights_squared_ += sqr(weight);
+ max_weight_ = std::max(max_weight_,std::abs(weight));
+ }
+ }
+
+ inline void statistics::reset() {
+ if (iteration_points_ == 0)
+ return;
+ double average = sum_weights_/iteration_points_;
+ double average_abs = sum_abs_weights_/iteration_points_;
+ double variance =
+ std::abs(sum_weights_squared_/iteration_points_ - sqr(sum_weights_/iteration_points_)) /
+ iteration_points_;
+ if (n_iterations_ == 0 ||
+ std::sqrt(variance)/average <= std::sqrt(average_weight_variance_)/average_weight_) {
+
+ if (n_iterations_ > 0) {
+ average_weight_ = average_weight_ + average;
+ average_abs_weight_ = average_abs_weight_ + average_abs;
+ average_weight_variance_ = average_weight_variance_ + variance;
+ } else {
+ average_weight_ = average;
+ average_abs_weight_ = average_abs;
+ average_weight_variance_ = variance;
+ }
+ ++n_iterations_;
+ }
+ sum_weights_ = 0.;
+ sum_abs_weights_ = 0.;
+ sum_weights_squared_ = 0.;
+ iteration_points_ = 0;
+ }
+
+ inline std::pair<double,double> statistics::current() const {
+ double average = sum_weights_/iteration_points_;
+ double variance =
+ std::abs(sum_weights_squared_/iteration_points_ - sqr(sum_weights_/iteration_points_)) /
+ iteration_points_;
+ std::pair<double,double> res;
+ if (n_iterations_ > 0) {
+ res.first = average_weight_ + average;
+ res.second = average_weight_variance_ + variance;
+ } else {
+ res.first = average;
+ res.second = variance;
+ }
+ res.first = (n_iterations_ == 0 ? 0. : res.first/n_iterations_);
+ res.second = (n_iterations_ == 0 ? 0. : res.second/n_iterations_);
+ return res;
+ }
+
+
+ template<class OStream>
+ void statistics::put(OStream& os) const {
+ os << average_weight_;
+ ostream_traits<OStream>::separator(os);
+ os << average_abs_weight_;
+ ostream_traits<OStream>::separator(os);
+ os << average_weight_variance_;
+ ostream_traits<OStream>::separator(os);
+ os << n_iterations_;
+ ostream_traits<OStream>::separator(os);
+ }
+
+ template<class IStream>
+ void statistics::get(IStream& is) {
+ is >> average_weight_
+ >> average_abs_weight_
+ >> average_weight_variance_
+ >> n_iterations_;
+ }
+
+}
diff --git a/Exsample2/exsample/utility.h b/Exsample2/exsample/utility.h
new file mode 100644
--- /dev/null
+++ b/Exsample2/exsample/utility.h
@@ -0,0 +1,322 @@
+// -*- C++ -*-
+//
+// utility.h is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
+//
+// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
+//
+// ExSample is licenced under version 2 of the GPL, see COPYING for details.
+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
+//
+//
+#ifndef EXSAMPLE_utility_h_included
+#define EXSAMPLE_utility_h_included
+
+#include "config.h"
+
+namespace exsample {
+
+ /// \brief Compile time conversion of unsigned long to bool
+ template<unsigned long>
+ struct static_binary {
+ enum { value = 1 };
+ };
+
+ /// \brief Compile time conversion of unsigned long to bool
+ template<>
+ struct static_binary<0> {
+ enum { value = 0 };
+ };
+
+ /// \brief Fixed-size, packed vector of bools
+ template<unsigned long bits>
+ struct bit_container {
+
+ enum {
+ /// the number of bits contained
+ n_bits = bits,
+ /// the number of bits fitting in a unsigned long
+ uint_bits = CHAR_BIT * sizeof(unsigned long),
+ /// the number of unsigned long segments needed
+ n_segments = bits / uint_bits + static_binary<bits % uint_bits>::value
+ };
+
+
+ /// the default constructor
+ bit_container() {
+ for (std::size_t i = 0; i < n_segments; ++i)
+ segments[i] = 0;
+ }
+
+ /// put all values to false
+ void reset() {
+ for (std::size_t i = 0; i < n_segments; ++i)
+ segments[i] = 0;
+ }
+
+ /// compare for equality
+ bool operator==(const bit_container& x) const {
+ for (std::size_t i = 0; i < n_segments; ++i)
+ if(segments[i] != x.segments[i])
+ return false;
+ return true;
+ }
+
+ /// compare for ordering
+ bool operator<(const bit_container& x) const {
+ for (std::size_t i = 0; i < n_segments; ++i)
+ if(segments[i] != x.segments[i])
+ return (segments[i] < x.segments[i]);
+ return false;
+ }
+
+ /// set the k'th bit
+ void bit(unsigned long k, bool value) {
+ assert(k<n_bits);
+ if (value)
+ segments[n_segments-k/uint_bits-1] |= (1ul << (k % uint_bits));
+ else
+ segments[n_segments-k/uint_bits-1] &= ~(1ul << (k % uint_bits));
+ }
+
+ /// get the k'th bit
+ bool bit(unsigned long k) {
+ assert(k<n_bits);
+ return (segments[n_segments-k/uint_bits-1] & (1ul << (k % uint_bits)));
+ }
+
+ /// print to ostream
+ template<class OStream>
+ void dump(OStream& os) const {
+ for ( unsigned int k = 0; k < n_segments; ++k )
+ os << segments[k] << " ";
+ }
+
+ private:
+
+ /// segments needed to keep the hash value
+ unsigned long segments[n_segments];
+
+ };
+
+ /// \brief square a number
+ template<class T>
+ T sqr(T x) {
+ return x*x;
+ }
+
+ /// \brief cube a number
+ template<class T>
+ T cube(T x) {
+ return x*x*x;
+ }
+
+
+ /// \brief Round a floating point value to an integer value of the
+ /// same type.
+ template<class T>
+ T round(T x) {
+ T f = std::floor(x);
+ T c = std::ceil(x);
+ if (x < (f+c)/2.)
+ return f;
+ return c;
+ }
+
+ /// \brief Calculate fast powers of two.
+ inline std::size_t two_to(std::size_t n) {
+ assert(n <= sizeof(std::size_t)*CHAR_BIT);
+ return (1 << n);
+ }
+
+ /// \brief separate quantities written to an ostream
+ template<class OStream>
+ struct ostream_traits {
+
+ /// put the separator to the ostream
+ static void separator(OStream& os) { os << "\n"; }
+
+ };
+
+#ifdef EXSAMPLE_has_ThePEG
+
+ /// \brief separate quantities written to a ThePEG::PersistentOStream
+ template<>
+ struct ostream_traits<ThePEG::PersistentOStream> {
+
+ /// put the separator to the ostream
+ static void separator(ThePEG::PersistentOStream&) { }
+
+ };
+
+#endif // EXSAMPLE_has_ThePEG
+
+ /// \brief Fast, zero memory-overhead one-dimensional
+ /// histogram with 2^n equally spaced bins
+ template<class Statistics>
+ struct fast_small_histogram {
+
+ /// default constructor
+ fast_small_histogram()
+ : depth(0), bins(0) {}
+
+ /// copy constructor
+ fast_small_histogram(const fast_small_histogram& x)
+ : depth(x.depth), bins(0) {
+ if (x.bins) {
+ bins.reset(new Statistics[two_to(depth)]);
+ for(std::size_t k = 0; k < two_to(depth); ++k)
+ bins[k] = x.bins[k];
+ }
+ }
+
+ /// assignment
+ fast_small_histogram& operator=(const fast_small_histogram& x) {
+ if (&x == this)
+ return *this;
+ depth = x.depth;
+ bins.reset(0);
+ if (x.bins) {
+ bins.reset(new Statistics[two_to(depth)]);
+ for(std::size_t k = 0; k < two_to(depth); ++k)
+ bins[k] = x.bins[k];
+ }
+ return *this;
+ }
+
+ /// construct from depth d, creating 2^d bins
+ explicit fast_small_histogram(std::size_t d)
+ : depth(d), bins(0) {
+ bins.reset(new Statistics[two_to(d)]);
+ }
+
+ /// return the bin from event belongs to given outer boundaries
+ Statistics& bin(double lower,
+ double upper,
+ double event) {
+ double thelower = lower;
+ double theupper = upper;
+ std::size_t bindex = 0;
+ std::size_t current_depth = 0;
+ while (true) {
+ double cut
+ = (thelower+theupper)/2.;
+ if (event < cut) {
+ theupper = cut;
+ } else {
+ thelower = cut;
+ bindex += two_to(depth-current_depth-1);
+ }
+ if(++current_depth == depth)
+ break;
+ }
+ return bins[bindex];
+ }
+
+ /// the depth, defining a histogram of 2^depth bins
+ std::size_t depth;
+
+ /// the contained statistics objects
+ boost::scoped_array<Statistics> bins;
+
+ /// put histogram to an ostream
+ template<class OStream>
+ void put(OStream& os) const {
+ os << depth;
+ ostream_traits<OStream>::separator(os);
+ for (std::size_t k = 0; k < two_to(depth); ++k) {
+ bins[k].put(os);
+ }
+ }
+
+ /// get histogram from an istream
+ template<class IStream>
+ void get(IStream& is) {
+ is >> depth;
+ bins.reset(new Statistics[two_to(depth)]);
+ for(std::size_t k = 0; k < two_to(depth); ++k) {
+ bins[k].get(is);
+ }
+ }
+
+ };
+
+ /// \brief Generalize the transform algorithm to only apply
+ /// depending on a range of flags accompanying the input range
+ template<class FirstInputIterator,
+ class SecondInputIterator,
+ class FlagIterator,
+ class OutputIterator,
+ class BinaryOperation>
+ OutputIterator conditional_transform(FirstInputIterator first1,
+ FirstInputIterator last1,
+ SecondInputIterator first2,
+ FlagIterator firstf,
+ OutputIterator result,
+ BinaryOperation binary_op) {
+ for (; first1 != last1; ++first1, ++first2, ++firstf, ++result)
+ if (*firstf)
+ *result = binary_op(*first1, *first2);
+ return result;
+ }
+
+ /// \brief calculate a volume given lower left and upper right
+ /// corners
+ inline double volume(const std::vector<double>& lower_left,
+ const std::vector<double>& upper_right) {
+ std::vector<double> delta;
+ std::transform(upper_right.begin(),upper_right.end(),
+ lower_left.begin(),std::back_inserter(delta),
+ std::minus<double>());
+ return
+ std::accumulate(delta.begin(),delta.end(),1.,std::multiplies<double>());
+ }
+
+ /// \brief calculate a volume given lower left and upper right
+ /// corners, taking into account only part of the dimensions, which
+ /// are flagged with true in the correspponding random access
+ /// container
+ inline double volume(const std::vector<double>& lower_left,
+ const std::vector<double>& upper_right,
+ const std::vector<bool>& flags) {
+ std::vector<double> delta;
+ conditional_transform(upper_right.begin(),upper_right.end(),
+ lower_left.begin(),flags.begin(),
+ std::back_inserter(delta),
+ std::minus<double>());
+ return
+ std::accumulate(delta.begin(),delta.end(),1.,std::multiplies<double>());
+ }
+
+
+ /// \brief Exception thrown if the maximum number of attempts to
+ /// select a cell has been reached.
+ struct selection_maxtry{};
+
+ /// \brief Exception thrown, if the maximum number of misses has
+ /// been reached.
+ struct hit_and_miss_maxtry{};
+
+ /// \brief Random generator traits.
+ template<class Random>
+ struct rnd_generator {
+
+ ///Generate uniform random number on [0,1]
+ double operator()() const {
+ return Random::rnd();
+ }
+
+ ///Generate uniform random number on [0,a]
+ double operator()(double a) const {
+ return a*Random::rnd();
+ }
+
+ ///Generate uniform random number on [a,b]
+ double operator()(double a, double b) const {
+ return (a + (b-a)*Random::rnd());
+ }
+
+ };
+
+}
+
+#endif // EXSAMPLE_utility_h_included
diff --git a/Exsample2/exsampler.cc b/Exsample2/exsampler.cc
new file mode 100644
--- /dev/null
+++ b/Exsample2/exsampler.cc
@@ -0,0 +1,312 @@
+// -*- C++ -*-
+
+// exsample2/exsampler.cc is part of matchbox
+// (C) 2008 Simon Platzer -- sp@particle.uni-karlsruhe.de
+
+#include "exsampler.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Interface/Parameter.h"
+
+#include "ThePEG/Handlers/StandardXComb.h"
+#include "ThePEG/PDT/ParticleData.h"
+
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+
+#include <sstream>
+
+using namespace matchbox;
+using namespace exsample;
+
+exsampler::exsampler()
+ : presampling_points_(2000),
+ freeze_grid_(0),
+ efficiency_threshold_(.95),
+ gain_threshold_(.1),
+ samplers_(), eg_wrappers_(),
+ bin_selector_(), integrals_(),
+ missing_events_(), oversampling_bins_(),
+ sum_weights_(0),
+ integral_(0.), integral_abs_(0.), variance_(0.),
+ max_integral_(0.), last_bin_(-1),
+ last_weight_(0.) {}
+
+exsampler::~exsampler() {}
+
+ThePEG::IBPtr exsampler::clone() const {
+ return ThePEG::new_ptr(*this);
+}
+
+ThePEG::IBPtr exsampler::fullclone() const {
+ return ThePEG::new_ptr(*this);
+}
+
+std::string exsampler::process(int bin) const {
+ std::ostringstream os("");
+ const ThePEG::StandardEventHandler& eh = *eventHandler();
+ const ThePEG::StandardXComb& xc = *eh.xCombs()[bin];
+ os << xc.mePartonData()[0]->PDGName() << " "
+ << xc.mePartonData()[1]->PDGName() << " -> ";
+ for ( ThePEG::cPDVector::const_iterator pid =
+ xc.mePartonData().begin() + 2;
+ pid != xc.mePartonData().end(); ++pid )
+ os << (**pid).PDGName() << " ";
+ return os.str();
+}
+
+void exsampler::dofinish() {
+ ThePEG::SamplerBase::dofinish();
+
+ integral_ = 0.;
+ integral_abs_ = 0.;
+ variance_ = 0.;
+ max_integral_ = 0.;
+
+ int b = 0;
+
+ bool compensating = false;
+
+ for (std::vector<sampler_type>::iterator s = samplers_.begin();
+ s != samplers_.end(); ++s, ++b) {
+
+ std::cout << "integrated cross section for "
+ << process(b) << "is ("
+ << samplers_[b].integral() << " +/- "
+ << samplers_[b].integral_uncertainty()
+ << ") nb"
+ << (s->compensating() ? " (*)" : "") << "\n";
+
+ if ( s->compensating() )
+ compensating = true;
+
+ s->finalize();
+ integral_ += s->integral();
+ integral_abs_ += std::abs(s->integral());
+ variance_ += s->integral_variance();
+ max_integral_ = std::max(max_integral_,s->integral());
+ }
+
+ std::cout << "integrated cross section is ("
+ << integral_ << " +/- "
+ << std::sqrt(variance_) << ") nb\n" << std::flush;
+
+ if ( compensating )
+ std::cout << "warning -- samplers marked with (*) are still in compensating mode\n" << std::flush;
+
+}
+
+void exsampler::doinitrun() {
+ ThePEG::SamplerBase::doinitrun();
+ initialize();
+}
+
+void exsampler::persistentOutput(ThePEG::PersistentOStream & os) const {
+ // *** ATTENTION *** os << ; // Add all member variable which should be written persistently here.
+ os << presampling_points_
+ << freeze_grid_ << efficiency_threshold_
+ << gain_threshold_ << samplers_.size();
+
+ for (std::vector<sampler_type>::const_iterator s = samplers_.begin();
+ s != samplers_.end(); ++s)
+ s->put(os);
+
+}
+
+void exsampler::persistentInput(ThePEG::PersistentIStream & is, int) {
+ // *** ATTENTION *** is >> ; // Add all member variable which should be read persistently here.
+ is >> presampling_points_
+ >> freeze_grid_ >> efficiency_threshold_
+ >> gain_threshold_;
+ std::size_t nsamplers;
+ is >> nsamplers;
+ samplers_.resize(nsamplers);
+
+ for (std::vector<sampler_type>::iterator s = samplers_.begin();
+ s != samplers_.end(); ++s)
+ s->get(is);
+
+}
+
+ThePEG::ClassDescription<exsampler> exsampler::initexsampler;
+// Definition of the static class description member.
+
+void exsampler::Init() {
+
+ static ThePEG::ClassDocumentation<exsampler> documentation
+ ("There is no documentation for the exsampler class");
+
+ static ThePEG::Parameter<exsampler,unsigned long> interfacepresampling_points
+ ("presampling_points",
+ "Set the number of presampling points per cell",
+ &exsampler::presampling_points_, 2000, 0, 0,
+ false, false, ThePEG::Interface::lowerlim);
+
+ static ThePEG::Parameter<exsampler,unsigned long> interfacefreeze_grid
+ ("freeze_grid",
+ "Set the number of events after which the grid should be frozen",
+ &exsampler::freeze_grid_, 0, 0, 0,
+ false, false, ThePEG::Interface::lowerlim);
+
+ static ThePEG::Parameter<exsampler,double> interfaceefficiency_threshold
+ ("efficiency_threshold",
+ "Set the efficiency threshold",
+ &exsampler::efficiency_threshold_, .95, 0., 1.,
+ false, false, ThePEG::Interface::limited);
+
+ static ThePEG::Parameter<exsampler,double> interfacegain_threshold
+ ("gain_threshold",
+ "Set the gain threshold",
+ &exsampler::gain_threshold_, .1, 0., 1.,
+ false, false, ThePEG::Interface::limited);
+
+}
+
+void exsampler::initialize() {
+
+ integral_ = 0.;
+ integral_abs_ = 0.;
+ variance_ = 0.;
+ max_integral_ = 0.;
+
+ if (samplers_.empty())
+ samplers_.resize(eventHandler()->nBins());
+
+ eg_wrappers_.resize(eventHandler()->nBins());
+ integrals_.resize(eventHandler()->nBins());
+ missing_events_.resize(eventHandler()->nBins(),0);
+
+ for (int b = 0; b < eventHandler()->nBins(); ++b) {
+
+ last_bin_ = b;
+
+ eg_wrappers_[b] = eg_exsample2_wrapper(b,eventHandler());
+ samplers_[b].function(&eg_wrappers_[b]);
+
+ samplers_[b].sampling_parameters().presampling_points = presampling_points_;
+ samplers_[b].sampling_parameters().freeze_grid = freeze_grid_;
+ samplers_[b].sampling_parameters().maxtry = eventHandler()->maxLoop();
+ samplers_[b].sampling_parameters().efficiency_threshold = efficiency_threshold_;
+ samplers_[b].sampling_parameters().gain_threshold = gain_threshold_;
+
+ samplers_[b].initialize(dummy_);
+ integrals_[b] = samplers_[b].integral();
+ integral_ += samplers_[b].integral();
+ integral_abs_ += std::abs(samplers_[b].integral());
+ variance_ += samplers_[b].integral_variance();
+ max_integral_ = std::max(max_integral_,samplers_[b].integral());
+
+ std::cout << "estimated cross section for " << process(b) << "is ("
+ << samplers_[b].integral() << " +/- "
+ << samplers_[b].integral_uncertainty()
+ << ") nb\n";
+
+ }
+
+ std::cout << "estimated cross section is ("
+ << integral_ << " +/- "
+ << std::sqrt(variance_) << ") nb\n" << std::flush;
+
+ double sum = 0.;
+
+ for (int b = 0; b < eventHandler()->nBins(); ++b) {
+ sum += std::abs(samplers_[b].integral());
+ bin_selector_.insert(std::make_pair(sum/integral_abs_,b));
+ }
+
+}
+
+void exsampler::update() {
+
+ double old_integral_abs = integral_abs_;
+ std::vector<double> old_integrals = integrals_;
+
+ integral_ = 0.;
+ integral_abs_ = 0.;
+ variance_ = 0.;
+ max_integral_ = 0.;
+
+ for (std::vector<sampler_type>::iterator s = samplers_.begin();
+ s != samplers_.end(); ++s) {
+
+ integral_ += s->integral();
+ integral_abs_ += std::abs(s->integral());
+ variance_ += s->integral_variance();
+ max_integral_ = std::max(max_integral_,s->integral());
+
+ }
+
+ bin_selector_.clear();
+
+ double sum = 0.;
+
+ for (int b = 0; b < eventHandler()->nBins(); ++b) {
+
+ sum += std::abs(samplers_[b].integral());
+ integrals_[b] = samplers_[b].integral();
+ bin_selector_.insert(std::make_pair(sum/integral_abs_,b));
+
+ long missing =
+ static_cast<long>(round(samplers_[b].stats().accepted()*
+ ((old_integral_abs * std::abs(integrals_[b])) /
+ (integral_abs_ * std::abs(old_integrals[b]))-1.)));
+
+ missing_events_[b] += missing;
+ if (missing_events_[b] > 0)
+ oversampling_bins_.insert(b);
+
+ }
+
+}
+
+double exsampler::generate() {
+
+ while (true) {
+
+ try {
+
+ if (oversampling_bins_.empty()) {
+ last_bin_ = bin_selector_.upper_bound(ThePEG::UseRandom::rnd())->second;
+ if (missing_events_[last_bin_] < 0) {
+ ++missing_events_[last_bin_];
+ continue;
+ }
+ } else {
+ last_bin_ = *oversampling_bins_.begin();
+ if (--missing_events_[last_bin_] == 0)
+ oversampling_bins_.erase(last_bin_);
+ }
+
+ last_weight_ = samplers_[last_bin_].generate(dummy_);
+ lastPoint() = samplers_[last_bin_].last_point();
+ break;
+
+ } catch (selection_maxtry&) {
+
+ throw maxtry_exception()
+ << "The maximum number of attempts to select a cell was exceeded\n"
+ << "for process " << process(last_bin_)
+ << ThePEG::Exception::eventerror;
+
+ } catch (hit_and_miss_maxtry&) {
+
+ throw maxtry_exception()
+ << "The maximum number of attempts to select an event was exceeded\n"
+ << "process " << process(last_bin_)
+ << ThePEG::Exception::eventerror;
+
+ } catch (generator_update&) {
+
+ update();
+ continue;
+
+ } catch (...) {
+ throw;
+ }
+
+ }
+
+ sum_weights_ += last_weight_;
+ return last_weight_;
+
+}
+
diff --git a/Exsample2/exsampler.h b/Exsample2/exsampler.h
new file mode 100644
--- /dev/null
+++ b/Exsample2/exsampler.h
@@ -0,0 +1,353 @@
+// -*- C++ -*-
+
+// exsample2/exsampler.h is part of matchbox
+// (C) 2008 Simon Platzer -- sp@particle.uni-karlsruhe.de
+
+#ifndef matchbox_exsample2_exsampler_h
+#define matchbox_exsample2_exsampler_h
+
+#include <map>
+
+#include "ThePEG/Handlers/SamplerBase.h"
+#include "ThePEG/Handlers/StandardEventHandler.h"
+#include "ThePEG/Utilities/Exception.h"
+#include "ThePEG/Repository/UseRandom.h"
+
+#include "exsample/generator.h"
+
+namespace matchbox {
+
+ /**
+ * exsample2 function interface conforming
+ * wrapper around standard event handler
+ */
+ struct eg_exsample2_wrapper {
+
+ inline eg_exsample2_wrapper ()
+ : bin_(-1), event_handler_() {}
+
+ inline eg_exsample2_wrapper (int b,
+ ThePEG::tStdEHPtr eh)
+ : bin_(b), event_handler_(eh) {}
+
+
+ /**
+ * Evaluate with given random numbers.
+ */
+ inline double evaluate (const std::vector<double>& p) const {
+ double ret;
+ try {
+ ret = event_handler_->dSigDR(p) / ThePEG::nanobarn;
+ } catch (ThePEG::Veto&) {
+ ret = 0.0;
+ } catch (...) {
+ throw;
+ }
+ return ret;
+ }
+
+ /**
+ * Return the dimensionality in random
+ * numbers needed to generate events
+ */
+ inline std::size_t dimension () const {
+ return event_handler_->nDim(bin_);
+ }
+
+ /**
+ * Return the lower left and upper right
+ * corners of the support of this function
+ */
+ inline std::pair<std::vector<double>,std::vector<double> > support () const {
+ std::vector<double> lower(dimension(),0.);
+ std::vector<double> upper(dimension(),1.);
+ return std::make_pair(lower,upper);
+ }
+
+ /**
+ * Indicate start of presampling
+ */
+ inline void start_presampling () { }
+
+ /**
+ * Indicate end of presampling
+ */
+ inline void stop_presampling () { }
+
+ private:
+
+ int bin_;
+ ThePEG::tStdEHPtr event_handler_;
+
+ };
+
+ /**
+ * Here is the documentation of the exsampler class.
+ *
+ * @see \ref exsamplerInterfaces "The interfaces"
+ * defined for exsampler.
+ */
+ class exsampler
+ : public ThePEG::SamplerBase {
+
+ public:
+
+ /** @name Standard constructors and destructors. */
+ //@{
+ /**
+ * The default constructor.
+ */
+ exsampler();
+
+ /**
+ * The destructor.
+ */
+ virtual ~exsampler();
+ //@}
+
+ public:
+
+ typedef exsample::generator<eg_exsample2_wrapper,ThePEG::UseRandom> sampler_type;
+
+ /** @name Virtual functions from SamplerBase. */
+ //@{
+ /**
+ * Initialize the the sampler, possibly doing presampling of the
+ * phase space.
+ */
+ virtual void initialize();
+
+ /**
+ * Generarate a new phase space point and return a weight associated
+ * with it. This weight should preferably be 1.
+ */
+ virtual double generate();
+
+ /**
+ * Reject the last chosen phase space point.
+ */
+ virtual inline void rejectLast() {
+ samplers_[last_bin_].reject();
+ sum_weights_ -= last_weight_;
+ }
+
+ /**
+ * If the sampler is able to sample several different functions
+ * separately, this function should return the last chosen
+ * function. This default version always returns 0.
+ */
+ virtual inline int lastBin() const { return last_bin_; }
+
+ /**
+ * Return the total integrated cross section determined from the
+ * Monte Carlo sampling so far.
+ */
+ virtual inline ThePEG::CrossSection integratedXSec() const {
+ return integral_ * ThePEG::nanobarn;
+ }
+
+ /**
+ * Return the error on the total integrated cross section determined
+ * from the Monte Carlo sampling so far.
+ */
+ virtual ThePEG::CrossSection integratedXSecErr() const {
+ return std::sqrt(variance_) * ThePEG::nanobarn;
+ }
+
+ /**
+ * Return the overestimated integrated cross section.
+ */
+ virtual ThePEG::CrossSection maxXSec() const {
+ return max_integral_ * ThePEG::nanobarn;
+ }
+
+ /**
+ * Return the sum of the weights returned by generate() so far (of
+ * the events that were not rejeted).
+ */
+ virtual inline double sumWeights() const { return sum_weights_; }
+ //@}
+
+ public:
+
+ /** @name Functions used by the persistent I/O system. */
+ //@{
+ /**
+ * Function used to write out object persistently.
+ * @name os the persistent output stream written to.
+ */
+ void persistentOutput(ThePEG::PersistentOStream & os) const;
+
+ /**
+ * Function used to read in object persistently.
+ * @name is the persistent input stream read from.
+ * @name version the version number of the object when written.
+ */
+ void persistentInput(ThePEG::PersistentIStream & is, int version);
+ //@}
+
+ /**
+ * The standard Init function used to initialize the interfaces.
+ * Called exactly once for each class by the class description system
+ * before the main function starts or
+ * when this class is dynamically loaded.
+ */
+ static void Init();
+
+ public:
+
+ struct maxtry_exception : public ThePEG::Exception {};
+
+ protected:
+
+ /** @name Standard Interfaced functions. */
+ //@{
+
+ /**
+ * Initialize this object. Called in the run phase just before
+ * a run begins.
+ */
+ virtual void doinitrun();
+
+ /**
+ * Finalize this object. Called in the run phase just after a
+ * run has ended. Used eg. to write out statistics.
+ */
+ virtual void dofinish();
+
+ //@}
+
+ protected:
+
+ /** @name Clone Methods. */
+ //@{
+ /**
+ * Make a simple clone of this object.
+ * @return a pointer to the new object.
+ */
+ virtual ThePEG::IBPtr clone() const;
+
+ /** Make a clone of this object, possibly modifying the cloned object
+ * to make it sane.
+ * @return a pointer to the new object.
+ */
+ virtual ThePEG::IBPtr fullclone() const;
+ //@}
+
+ private:
+
+ /// update integrals and stuff
+ void update ();
+
+ /// get process by bin id
+ std::string process(int bin) const;
+
+ ///@name interface helpers to adaption_info objects
+ //@{
+
+ unsigned long presampling_points_;
+
+ unsigned long freeze_grid_;
+
+ double efficiency_threshold_;
+
+ double gain_threshold_;
+
+ //@}
+
+ /// map bin ids to generators
+ std::vector<sampler_type> samplers_;
+
+ /// map bin ids to wrappers used
+ std::vector<eg_exsample2_wrapper> eg_wrappers_;
+
+ /// map used to select bins
+ std::map<double,int> bin_selector_;
+
+ /// current integrals of the bins
+ std::vector<double> integrals_;
+
+ /// map bins to number of missing events
+ std::vector<long> missing_events_;
+
+ /// bins which do need oversampling
+ std::set<int> oversampling_bins_;
+
+ /// the sum of accepted weights
+ double sum_weights_;
+
+ /// the total integrated corss section
+ double integral_;
+
+ /// sum of absoulte values of integrals
+ double integral_abs_;
+
+ /// the ucertainty on the integrated cross section
+ double variance_;
+
+ /// the maximum integrated cross section
+ double max_integral_;
+
+ /// the last selected bin
+ int last_bin_;
+
+ /// the last generated weight
+ double last_weight_;
+
+ /// dummy statistics to fix functionality required by rewrite
+ exsample::statistics dummy_;
+
+ private:
+
+ /**
+ * The static object used to initialize the description of this class.
+ * Indicates that this is an abstract class with persistent data.
+ */
+ static ThePEG::ClassDescription<exsampler> initexsampler;
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ exsampler & operator=(const exsampler &);
+
+ };
+
+}
+
+#include "ThePEG/Utilities/ClassTraits.h"
+
+namespace ThePEG {
+
+ /** @cond TRAITSPECIALIZATIONS */
+
+ /** This template specialization informs ThePEG about the
+ * base classes of exsampler. */
+ template <>
+ struct BaseClassTrait<matchbox::exsampler,1> {
+ /** Typedef of the first base class of exsampler. */
+ typedef ThePEG::SamplerBase NthBase;
+ };
+
+ /** This template specialization informs ThePEG about the name of
+ * the exsampler class and the shared object where it is defined. */
+ template <>
+ struct ClassTraits<matchbox::exsampler>
+ : public ClassTraitsBase<matchbox::exsampler> {
+ /** Return a platform-independent class name */
+ static string className() { return "matchbox::exsampler"; }
+ /**
+ * The name of a file containing the dynamic library where the class
+ * exsampler is implemented. It may also include several, space-separated,
+ * libraries if the class exsampler depends on other classes (base classes
+ * excepted). In this case the listed libraries will be dynamically
+ * linked in the order they are specified.
+ */
+ static string library() { return "HwExsample2.so"; }
+ };
+
+ /** @endcond */
+
+}
+
+#endif /* matchbox_exsample2_exsampler_h */

File Metadata

Mime Type
text/x-diff
Expires
Mon, Jan 20, 11:01 PM (1 d, 17 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242788
Default Alt Text
(133 KB)

Event Timeline