File indexing completed on 2024-04-06 12:26:28
0001 #ifndef ErrorPropogationTypes_h
0002 #define ErrorPropogationTypes_h
0003
0004 #include "boost/operators.hpp"
0005 #include <cmath>
0006
0007 class count_t
0008 : boost::addable<count_t,
0009 boost::incrementable<count_t, boost::totally_ordered<count_t, boost::equivalent<count_t> > > > {
0010 private:
0011 unsigned long count;
0012
0013 public:
0014 count_t() : count(0) {}
0015 count_t(unsigned long n) : count(n) {}
0016 bool operator<(const count_t& R) const { return count < R.count; }
0017 count_t operator++() {
0018 ++count;
0019 return *this;
0020 }
0021 count_t operator+=(const count_t& R) {
0022 count += R.count;
0023 return *this;
0024 }
0025 unsigned long operator()() const { return count; }
0026 unsigned long error2() const { return count; }
0027 double error() const { return sqrt(count); }
0028 double relative_error() const { return 1 / sqrt(count); }
0029 };
0030
0031 template <class charT, class traits>
0032 inline std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& strm, const count_t& f) {
0033 strm << f() << "(" << f.error() << ")";
0034 return strm;
0035 }
0036
0037 template <class T>
0038 class stats_t
0039 : boost::arithmetic1<
0040 stats_t<T>,
0041 boost::arithmetic2<
0042 stats_t<T>,
0043 T,
0044 boost::partially_ordered<
0045 stats_t<T>,
0046 boost::partially_ordered<
0047 stats_t<T>,
0048 count_t,
0049 boost::equality_comparable<stats_t<T>, boost::equality_comparable<stats_t<T>, count_t> > > > > > {
0050 private:
0051 T value, err2;
0052
0053 public:
0054 stats_t() : value(0), err2(1) {}
0055 stats_t(count_t c) : value(c()), err2(c.error2()) {}
0056 stats_t(T q, T e2) : value(q), err2(e2) {}
0057 static stats_t from_relative_uncertainty2(T q, T re2) { return stats_t(q, q * q * re2); }
0058
0059 bool operator<(const stats_t& R) const { return value < R.value; }
0060 bool operator<(const count_t& R) const { return value < R(); }
0061 bool operator==(const stats_t& R) const { return value == R.value && err2 == R.err2; }
0062 bool operator==(const count_t& R) const { return value == R() && err2 == R.error2(); }
0063
0064 stats_t operator+=(const stats_t& R) {
0065 value += R.value;
0066 err2 += R.err2;
0067 return *this;
0068 }
0069 stats_t operator-=(const stats_t& R) {
0070 value -= R.value;
0071 err2 += R.err2;
0072 return *this;
0073 }
0074 stats_t operator*=(const stats_t& R) {
0075 err2 = R.err2 * value * value + err2 * R.value * R.value;
0076 value *= R.value;
0077 return *this;
0078 }
0079 stats_t operator/=(const stats_t& R) {
0080 value /= R.value;
0081 err2 = (err2 + value * value * R.err2) / (R.value * R.value);
0082 return *this;
0083 }
0084
0085 stats_t operator+=(const T& r) {
0086 value += r;
0087 return *this;
0088 }
0089 stats_t operator-=(const T& r) {
0090 value -= r;
0091 return *this;
0092 }
0093 stats_t operator*=(const T& r) {
0094 value *= r;
0095 err2 *= r * r;
0096 return *this;
0097 }
0098 stats_t operator/=(const T& r) {
0099 value /= r;
0100 err2 /= r * r;
0101 return *this;
0102 }
0103
0104 stats_t inverse() const { return stats_t(1. / value, err2 / pow(value, 4)); }
0105
0106 T operator()() const { return value; }
0107 T error2() const { return err2; }
0108 T error() const { return sqrt(err2); }
0109 T relative_error() const { return sqrt(err2) / value; }
0110 T sigmaFrom(const T& x) const { return fabs(value - x) / error(); }
0111 };
0112
0113 template <class charT, class traits, class T>
0114 inline std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& strm, const stats_t<T>& f) {
0115 strm << f() << "(" << f.error() << ")";
0116 return strm;
0117 }
0118
0119 #endif
