ErrorPropogationTypes.h

Go to the documentation of this file.

#ifndef ErrorPropogationTypes_h
#define ErrorPropogationTypes_h

#include "boost/operators.hpp"
#include <cmath>              

class count_t :
boost::addable<count_t,    
boost::incrementable<count_t, 
boost::totally_ordered<count_t, 
boost::equivalent<count_t> > > > {
  
 private: 
  unsigned long count;
 public:               
  count_t() : count(0) {}
  count_t(unsigned long n) : count(n) {}
  bool operator<(const count_t& R) const { return count < R.count;}
  count_t operator++() {++count; return *this;}                    
  count_t operator+=(const count_t& R) {count += R.count; return *this;}
  unsigned long operator()() const {return count;}                 
  unsigned long error2() const {return count;}                     
  double error() const {return sqrt(count);}                       
  double relative_error() const {return 1/sqrt(count);}            
};                                          
  
template <class charT, class traits> 
  inline
  std::basic_ostream<charT,traits>& operator<<(std::basic_ostream<charT,traits>& strm, const count_t& f)
{ strm << f() << "("<< f.error()<< ")"; return strm;}

template<class T>
class stats_t :  
    boost::arithmetic1<stats_t<T>, 
    boost::arithmetic2<stats_t<T>, T,
    boost::partially_ordered<stats_t<T>,
    boost::partially_ordered<stats_t<T>,count_t,
    boost::equality_comparable<stats_t<T>,      
    boost::equality_comparable<stats_t<T>,count_t> > > > > > {
  private:
  T value, err2;
  public:
  stats_t() : value(0), err2(1) {}
  stats_t(count_t c) : value(c()), err2(c.error2()) {}
  stats_t(T q, T e2) : value(q), err2(e2) {}
  static stats_t from_relative_uncertainty2(T q, T re2) { return stats_t(q,q*q*re2);}
  
  bool operator<(const stats_t& R) const { return value < R.value ;}
  bool operator<(const count_t& R) const { return value < R()     ;}
  bool operator==(const stats_t& R) const { return value == R.value && err2 == R.err2;}
  bool operator==(const count_t& R) const { return value == R()     && err2 == R.error2();}
  
  stats_t operator+=(const stats_t& R) { value += R.value;   err2 += R.err2; return *this;}
  stats_t operator-=(const stats_t& R) { value -= R.value;   err2 += R.err2; return *this;}
  stats_t operator*=(const stats_t& R) { err2 = R.err2*value*value + err2*R.value*R.value;        value *= R.value;  return *this;}
  stats_t operator/=(const stats_t& R) { value /= R.value;  err2 = (err2 + value*value*R.err2) / (R.value*R.value);  return *this;}
  
  stats_t operator+=(const T& r) { value += r; return *this;}
  stats_t operator-=(const T& r) { value -= r; return *this;}
  stats_t operator*=(const T& r) { value *= r; err2 *= r*r; return *this;}
  stats_t operator/=(const T& r) { value /= r; err2 /= r*r; return *this;}
  
  stats_t inverse() const { return stats_t(1./value, err2/pow(value,4));}
  
  T operator()() const {return value;}
  T error2() const {return err2;}
  T error() const {return sqrt(err2);}
  T relative_error() const {return sqrt(err2)/value;}
  T sigmaFrom(const T& x) const {return fabs(value-x)/error();}
};

template <class charT, class traits, class T> inline
  std::basic_ostream<charT,traits>& operator<<(std::basic_ostream<charT,traits>& strm, const stats_t<T>& f)
{ strm << f() << "("<< f.error()<< ")"; return strm; }

#endif