Estimate.h

//-*-C++-*-
/***************************************************************************
 *
 *   Copyright (C) 2003 by Willem van Straten
 *   Licensed under the Academic Free License version 2.1
 *
 ***************************************************************************/
 
// epsic/src/util/Estimate.h
 
#ifndef __Estimate_h
#define __Estimate_h
 
#include "PromoteTraits.h"
#include "Traits.h"
#include "true_math.h"
 
#include <iostream>
#include <math.h>
#include <limits>
 
// forward declarations
template <typename T, typename U=T> class MeanEstimate;
template <typename T, typename U=T> class MeanRadian;


template <typename T, typename U=T>
class Estimate
{
 
 public:
 
  typedef T val_type;
  typedef U var_type;

  T val;
  U var;

  Estimate (T _val=0, U _var=0) { val=_val; var=_var; }

  template <typename V, typename W>
  Estimate (const Estimate<V,W>& d) { val=d.val; var=d.var; }

  template <typename V, typename W>
  Estimate (const MeanEstimate<V,W>& m) { operator=( m.get_Estimate() ); }

  template <typename V, typename W>
  Estimate (const MeanRadian<V,W>& m) { operator=( m.get_Estimate() ); }

  const Estimate& operator= (const Estimate& d)
  { val=d.val; var=d.var; return *this; }

  void set_value (const T& t) { val = t; }

  T get_value () const { return val; }

  void set_variance (const U& u) { var = u; }

  U get_variance () const { return var; }

  void set_error (const U& u) { var = u*u; }

  U get_error () const { return sqrt(var); }

  T& operator [] (unsigned) { return val; }

  T operator [] (unsigned) const { return val; }

  const Estimate& operator+= (const Estimate& d)
  { val += d.val; var += d.var; return *this; }

  const Estimate& operator-= (const Estimate& d)
  { val -= d.val; var += d.var; return *this; }


  const Estimate& operator*= (const Estimate& d)
  { var=val*val*d.var+d.val*d.val*var; val*=d.val; return *this; }

  const Estimate& operator/= (const Estimate& d)
  { return operator *= (d.inverse()); }

  bool operator == (const Estimate& d) const
  { return val == d.val; }

  bool operator != (const Estimate& d) const
  { return ! operator == (d); }

  bool operator < (const Estimate& d) const
  { return val < d.val; }

  bool operator > (const Estimate& d) const
  { return val > d.val; }


  const Estimate inverse () const
  { T v=1.0/val; return Estimate (v,var*v*v*v*v); }
 
  friend const Estimate operator + (Estimate a, const Estimate& b)
  { return a+=b; }
 
  friend const Estimate operator - (Estimate a, const Estimate& b)
  { return a-=b; }
 
  friend const Estimate operator * (Estimate a, const Estimate& b)
  { return a*=b; }
  
  friend const Estimate operator / (Estimate a, const Estimate& b)
  { return a/=b; }

  friend const Estimate operator - (Estimate a)
  { return Estimate (-a.val, a.var); }
 
#if COMPLEX_SPECIALIZE
  friend const std::complex<Estimate> operator * (const std::complex<Estimate>& a, const std::complex<Estimate>& b)
  {
    return std::complex<Estimate> (
        a.real()*b.real() - a.imag()*b.imag(),
        a.real()*b.imag() + a.imag()*b.real()
        );
  }
 
  friend const Estimate norm (const std::complex<Estimate>& u)
  { return u.real()*u.real() + u.imag()*u.imag(); }
#endif

  friend const Estimate exp (const Estimate& u)
  { T val = ::exp (u.val); return Estimate (val, val*val*u.var); }

  friend const Estimate log (const Estimate& u)
  { return Estimate (::log (u.val), u.var/(u.val*u.val)); }

  friend const Estimate sqrt (const Estimate& u)
  { return Estimate (::sqrt (u.val), 0.25*u.var/fabs(u.val)); }

  friend const Estimate sin (const Estimate& u)
  { T val = ::sin (u.val); return Estimate (val, (1-val*val)*u.var); }

  friend const Estimate cos (const Estimate& u)
  { T val = ::cos (u.val); return Estimate (val, (1-val*val)*u.var); }

  friend const Estimate acos (const Estimate& u)
  { T val = ::acos (u.val); T del=-1.0/sqrt(1.0-u.val*u.val);
    return Estimate (val, del*del*u.var); }

  friend const Estimate atan (const Estimate& u)
  { T val = ::atan (u.val); T del=1/(1+u.val*u.val);
    return Estimate (val, del*del*u.var); }

  friend const Estimate atan2 (const Estimate& s, const Estimate& c)
  { T c2 = c.val*c.val;  T s2 = s.val*s.val;  T sc2 = c2+s2;
    return Estimate (::atan2 (s.val, c.val),(c2*s.var+s2*c.var)/(sc2*sc2)); }

  friend const Estimate sinh (const Estimate& u)
  { T val = ::sinh (u.val); return Estimate (val, (1+val*val)*u.var); }

  friend const Estimate cosh (const Estimate& u)
  { T val = ::cosh (u.val); return Estimate (val, (val*val-1)*u.var); }

  friend const Estimate atanh (const Estimate& u)
  { T val = ::atanh (u.val); T del=1/(1-u.val*u.val);
    return Estimate (val, del*del*u.var); }
 
  friend int finite (const Estimate& u)
  { return true_math::finite (u.val); }
 
  friend int isinf (const Estimate& u)
  { return isinf (u.val); }
 
  friend int isnan (const Estimate& u)
  { return isnan (u.val); }
 
  friend const T abs (const Estimate& u)
  { return std::abs (u.val); }
 
  friend const Estimate copysign (const Estimate& u, const Estimate& v)
  { return Estimate (::copysign (u.val, v.val), u.var); }
};
 

template <class T, class U> struct DatumTraits< Estimate<T,U> >
{
  ElementTraits<T> element_traits;
  static inline unsigned ndim () { return 1; }
  static inline T& element (Estimate<T,U>& t, unsigned) 
  { return t.val; }
  static inline const T& element (const Estimate<T,U>& t, unsigned)
  { return t.val; }
};
 
template <class T, class U, class V, class W>
class PromoteTraits< Estimate<T,U>, Estimate<V,W> >
{
  public:
    typedef Estimate< typename PromoteTraits<T,V>::promote_type, 
                      typename PromoteTraits<U,W>::promote_type > promote_type;
};
 
template <class T, class U, class V>
class PromoteTraits< Estimate<T,U>, V >
{
  public:
    typedef Estimate<typename PromoteTraits<T,V>::promote_type,U> promote_type;
};
 
template <class T, class U, class V>
class PromoteTraits< V, Estimate<T,U> >
{
  public:
    typedef Estimate<typename PromoteTraits<T,V>::promote_type,U> promote_type;
};
 
namespace std
{
  template <class T, class U>
  struct numeric_limits< Estimate<T,U> >
    : public numeric_limits<T> { };
}

template<typename T, typename U>
std::ostream& operator<< (std::ostream& ostr, const Estimate<T,U>& estimate)
{
  return ostr << "(" << estimate.val << "+-" << sqrt(estimate.var) << ")";
}
 
static inline bool expect (std::istream& is, char c)
{
  if (is.peek() != c) {
    is.setstate (std::ios::failbit);
    return false;
  }
  is.get();
  return true;
}
 
template<typename T, typename U>
std::istream& operator >> (std::istream& is, Estimate<T,U>& estimate)
{
  char open_brace;
  is >> open_brace;
 
  bool bracketed=true;
  if (open_brace != '(') {
    is.unget ();
    bracketed=false;
  }
 
  double value;
  is >> value;
 
  if (!expect(is, '+'))
    return is;
  if (!expect(is, '-'))
    return is;
 
  double error;
  is >> error;
 
  if (bracketed) 
  {
    if (!expect(is, ')'))
      return is;
  }
 
  estimate.set_value (value);
  estimate.set_error (error);
  return is;
}

std::string latex (const Estimate<double>&);

template <typename T, typename U>
class MeanEstimate
{
 public:
  T norm_val;
  U inv_var;

  MeanEstimate (T _val=0, U _var=0) { norm_val=_val; inv_var=_var; }

  MeanEstimate (const MeanEstimate& d)
  { norm_val=d.norm_val; inv_var=d.inv_var; }

  MeanEstimate (const Estimate<T,U>& d)
  { norm_val = 0; inv_var = 0; operator += (d); }

  const MeanEstimate& operator= (const MeanEstimate& d)
  { norm_val=d.norm_val; inv_var=d.inv_var; return *this; }

  const MeanEstimate& operator+= (const MeanEstimate& d)
  { norm_val += d.norm_val; inv_var += d.inv_var; return *this; }

  const MeanEstimate& operator+= (const Estimate<T,U>& d)
  { if (d.var){ U v=1.0/d.var; norm_val+=d.val*v; inv_var+=v; } return *this; }

  bool operator == (T _norm_val) const
  { return norm_val == _norm_val; }
 
  Estimate<T,U> get_Estimate () const
  { U var=0.0; if (inv_var != 0.0) var = 1.0/inv_var; 
    return Estimate<T,U> (norm_val*var, var); }
 
};

template<typename T, typename U>
std::ostream& operator<< (std::ostream& ostr, const MeanEstimate<T,U>& mean)
{
  return ostr << mean.get_Estimate();
}


template <typename T, typename U>
class MeanRadian
{
 
 public:

  MeanRadian () { }

  MeanRadian (const MeanRadian& mean)
  { operator = (mean); }

  const MeanRadian& operator= (const MeanRadian& mean)
  { cosine = mean.cosine; sine = mean.sine; return *this; }

  MeanRadian (const Estimate<T,U>& d)
  { operator = (d); }

  const MeanRadian& operator= (const Estimate<T,U>& d)
  { cosine = cos(d); sine = sin(d); return *this; }

  const MeanRadian& operator+= (const MeanRadian& d)
  { cosine += d.cosine; sine += d.sine; return *this; }

  const MeanRadian& operator+= (const Estimate<T,U>& d)
  { cosine += cos(d); sine += sin(d); return *this; }
 
  Estimate<T,U> get_Estimate () const
  { if (sine.norm_val==0 && cosine.norm_val==0) return Estimate<T,U>(0,0);
    return atan2 (sine.get_Estimate(), cosine.get_Estimate()); }
 
  Estimate<T,U> get_sin () const
  { return sine.get_Estimate(); }
 
  Estimate<T,U> get_cos () const
  { return cosine.get_Estimate(); }
 
 protected:
  MeanEstimate<T,U> cosine;

  MeanEstimate<T,U> sine;
 
};

template<typename T, typename U>
std::ostream& operator<< (std::ostream& ostr, const MeanRadian<T,U>& mean)
{
  return ostr << mean.get_Estimate();
}
 
#endif