JStats.hh

Go to the documentation of this file.

#ifndef __JTOOLS__JSTATS__
#define __JTOOLS__JSTATS__
 
#include <ostream>
#include <iomanip>
#include <cmath>
#include <limits>
#include <algorithm>
 
#include "JLang/JException.hh"
#include "JLang/JTitle.hh"
#include "JLang/JVectorize.hh"
#include "JLang/JManip.hh"
 
#include "JMath/JMath.hh"
 
 
/**
 * \author mdejong
 */
 
namespace JTOOLS {}
namespace JPP { using namespace JTOOLS; }
 
namespace JTOOLS {
 
  using JMATH::JMath;
  using JLANG::JTitle;
  using JLANG::JDivisionByZero;
  using JLANG::array_type;
  using JLANG::make_array;
 
 
  /**
   * Auxiliary data structure for running average and standard deviation.
   *
   * See Knuth TAOCP vol 2, 3rd edition, page 232.\n
   * This class acts as a zero-dimensional histogram.\n
   * Note that if a weight is used, it should strictly be positive.
   */
  struct JStats :
    public JTitle,
    public JMath<JStats>
  {
    /**
     * Constructor.
     *
     * \param  title            title
     */
    JStats(const JTitle& title     = "") :
      JTitle(title)
    {
      reset();
    }
 
 
    /**
     * Constructor.
     *
     * \param  title            title
     * \param  buffer           input data
     * \param  w                weight
     */
    template<class JElement_t, class JAllocator_t>
    JStats(const JTitle& title,
           const array_type<JElement_t, JAllocator_t>& buffer,
           const double  w         = 1.0) :
      JTitle(title)
    {
      reset();
 
      put(buffer, w);
    }
 
 
    /**
     * Reset.
     */
    void reset()
    {
      mean  = 0.0;
      sigma = 0.0;
      total = 0.0;
      count = 0;
      xmin  = std::numeric_limits<double>::max();
      xmax  = std::numeric_limits<double>::lowest();
      wmin  = std::numeric_limits<double>::max();
      wmax  = std::numeric_limits<double>::lowest();
    }
 
 
    /**
     * Add stats.
     *
     * \param  Q           stats
     * \return             this stats
     */
    JStats& add(const JStats& Q)
    {
      mean  += Q.mean;
      sigma += Q.sigma;
      total += Q.total;
      count += Q.count;
      xmin   = std::min(xmin, Q.xmin);
      xmax   = std::max(xmax, Q.xmax);
      wmin   = std::min(wmin, Q.wmin);
      wmax   = std::max(wmax, Q.wmax);
 
      return *this;
    }
 
 
    /**
     * Put value.
     *
     * \param  x           value
     * \param  w           weight
     */
    void put(const double x, const double w = 1.0)
    {
      total += w;
      count += 1;
 
      if (count == 1) {
 
        mean  = x;
        sigma = 0.0;
 
      } else {
 
        const double new_mean  = mean   +  w * (x - mean) / total;
        const double new_sigma = sigma  +  w * (x - mean) * (x - new_mean);
    
        // set up for next iteration
 
        mean  = new_mean; 
        sigma = new_sigma;
      }
 
      xmin = std::min(xmin, x);
      xmax = std::max(xmax, x);
      wmin = std::min(wmin, w);
      wmax = std::max(wmax, w);
    }
 
 
    /**
     * Put data.
     *
     * \param  buffer           input data
     * \param  w                weight
     */
    template<class JElement_t, class JAllocator_t>
    void put(const array_type<JElement_t, JAllocator_t>& buffer,
             const double w = 1.0)
    {
      for (typename array_type<JElement_t, JAllocator_t>::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
        put(*i, w);
      }
    }
 
 
    /**
     * Get total count.
     *
     * \return             count
     */
    long long int getCount() const
    {
      return count;
    }
 
 
    /**
     * Get total weight.
     *
     * \return             weight
     */
    double getTotal() const
    {
      return total;
    }
 
 
    /**
     * Get minimum value.
     *
     * \return             minimum value
     */
    double getXmin() const
    {
      return xmin;
    }
 
 
    /**
     * Get maximum value.
     *
     * \return             maximum value
     */
    double getXmax() const
    {
      return xmax;
    }
 
 
    /**
     * Get minimum weight.
     *
     * \return             minimum weight
     */
    double getWmin() const
    {
      return wmin;
    }
 
 
    /**
     * Get maximum weight.
     *
     * \return             maximum weight
     */
    double getWmax() const
    {
      return wmax;
    }
 
 
    /**
     * Get mean value.
     *
     * \return             mean value
     */
    double getMean() const
    {
      if (count != 0.0)
        return mean;
      else
        THROW(JDivisionByZero, "JStats::getMean() count " << count);
    }
 
 
    /**
     * Get mean value.
     *
     * \param  value       default value
     * \return             mean value
     */
    double getMean(const double value) const
    {
      if (count != 0.0)
        return getMean();
      else
        return value;
    }
 
 
    /**
     * Get standard deviation
     *
     * \return             standard deviation
     */
    double getSTDev() const
    {
      if (count > 1)
        return sqrt(count * sigma/(total * (count - 1)));
      else
        THROW(JDivisionByZero, "JStats::getSTDev() count " << count);
    }
 
 
    /**
     * Get standard deviation
     *
     * \param  value       default value
     * \return             standard deviation
     */
    double getSTDev(const double value) const
    {
      if (count > 1)
        return getSTDev();
      else
        return value;
    }
 
 
    /**
     * Get maximal deviation from average.
     *
     * \param  relative    if true, relative to average, else absolute
     * \return             deviation
     */
    double getDeviation(const bool relative = true) const
    {
      if (relative)
        return std::max(getXmax() - getMean(), getMean() - getXmin());
      else
        return getXmax() - getXmin();
    }
 
 
    /**
     * Test relative accuracy.
     *
     * \param  precision   relative precision
     * \return             true if reached accuracy; else false
     */
    bool hasAccuracy(const double precision) const
    {
      return getCount() > 1 && getSTDev() < precision * getMean();
    }
 
 
    /**
     * Print stats.
     *
     * \param  out         output stream
     * \param  lpr         long print
     */
    std::ostream& print(std::ostream& out, bool lpr = true) const
    {
      using namespace std;
 
      const int nc = getTitle().size();
 
      if (lpr) {
        out << setw(nc) << left  << " "             << ' '
            << setw(12) << left  << "    mean"      << ' '
            << setw(12) << left  << "    STD"       << ' '
            << setw(12) << left  << "  deviation"   << endl;
      }
 
      out << setw(nc) << left << getTitle()          << ' '
          << SCIENTIFIC(12,5) << getMean()           << ' '
          << SCIENTIFIC(12,5) << getSTDev()          << ' '
          << SCIENTIFIC(12,5) << getDeviation(false) << endl;
 
      return out;
    }
 
 
    /**
     * Print stats.
     *
     * \param  out           output stream
     * \param  stats         stats
     * \return               output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JStats& stats)
    {
      return stats.print(out, getLongprint(out));
    }
 
  protected:
    double        mean;
    double        sigma;
    double        total;
    long long int count;
    double        xmin;
    double        xmax;
    double        wmin;
    double        wmax;
  };
}
 
#endif