Sampling/exsample/statistics.icc

0001 // -*- C++ -*-
0002 //
0003 // statistics.icc is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
0004 //
0005 // Copyright (C) 2008-2019 Simon Platzer -- simon.plaetzer@desy.de, The Herwig Collaboration
0006 //
0007 // ExSample is licenced under version 3 of the GPL, see COPYING for details.
0008 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0009 //
0010 //
0011 namespace exsample {
0012
0013   inline statistics::statistics()
0014     : average_weight_(0.), average_abs_weight_(0.), average_weight_variance_(-1.),
0015       iteration_points_(0),
0016       attempted_(0), accepted_(0), accepted_negative_(0),
0017       sum_weights_(0.), sum_abs_weights_(0.), sum_weights_squared_(0.),
0018       max_weight_(0.), n_iterations_(0) {
0019   }
0020
0021   inline void statistics::presampled(double weight) {
0022     ++iteration_points_;
0023     sum_weights_ += weight;
0024     sum_abs_weights_ += std::abs(weight);
0025     sum_weights_squared_ += sqr(weight);
0026     max_weight_ = std::max(max_weight_,std::abs(weight));
0027   }
0028
0029   inline void statistics::select(double weight,
0030                  bool calculate_integral) {
0031     ++attempted_;
0032     if (calculate_integral) {
0033       ++iteration_points_;
0034       sum_weights_ += weight;
0035       sum_abs_weights_ += std::abs(weight);
0036       sum_weights_squared_ += sqr(weight);
0037       max_weight_ = std::max(max_weight_,std::abs(weight));
0038     }
0039   }
0040
0041   inline void statistics::reset() {
0042     if (iteration_points_ == 0)
0043       return;
0044     double average = sum_weights_/iteration_points_;
0045     double average_abs = sum_abs_weights_/iteration_points_;
0046     double variance =
0047       std::abs(sum_weights_squared_/iteration_points_ - sqr(sum_weights_/iteration_points_)) /
0048       iteration_points_;
0049     if (n_iterations_ == 0 ||
0050     std::sqrt(variance)/average <= std::sqrt(average_weight_variance_)/average_weight_) {
0051
0052       if (n_iterations_ > 0) {
0053     average_weight_ = average_weight_ + average;
0054     average_abs_weight_ = average_abs_weight_ + average_abs;
0055     average_weight_variance_ = average_weight_variance_ + variance;
0056       } else {
0057     average_weight_ = average;
0058     average_abs_weight_ = average_abs;
0059     average_weight_variance_ = variance;
0060       }
0061       ++n_iterations_;
0062     }
0063     sum_weights_ = 0.;
0064     sum_abs_weights_ = 0.;
0065     sum_weights_squared_ = 0.;
0066     iteration_points_ = 0;
0067   }
0068
0069   inline std::pair<double,double> statistics::current() const {
0070     double average = sum_weights_/iteration_points_;
0071     double variance =
0072       std::abs(sum_weights_squared_/iteration_points_ - sqr(sum_weights_/iteration_points_)) /
0073       iteration_points_;
0074     std::pair<double,double> res;
0075     if (n_iterations_ > 0) {
0076       res.first = average_weight_ + average;
0077       res.second = average_weight_variance_ + variance;
0078     } else {
0079       res.first = average;
0080       res.second = variance;
0081     }
0082     res.first = (n_iterations_ == 0 ? 0. : res.first/n_iterations_);
0083     res.second = (n_iterations_ == 0 ? 0. : res.second/n_iterations_);
0084     return res;
0085   }
0086
0087
0088   template<class OStream>
0089   void statistics::put(OStream& os) const {
0090     os << average_weight_;
0091     ostream_traits<OStream>::separator(os);
0092     os << average_abs_weight_;
0093     ostream_traits<OStream>::separator(os);
0094     os << average_weight_variance_;
0095     ostream_traits<OStream>::separator(os);
0096     os << n_iterations_;
0097     ostream_traits<OStream>::separator(os);
0098   }
0099
0100   template<class IStream>
0101   void statistics::get(IStream& is) {
0102     is >> average_weight_
0103        >> average_abs_weight_
0104        >> average_weight_variance_
0105        >> n_iterations_;
0106   }
0107
0108 }