JQuantile.hh

Go to the documentation of this file.

#ifndef __JTOOLS__JQUANTILE__
#define __JTOOLS__JQUANTILE__
 
#include <ostream>
#include <iomanip>
#include <cmath>
#include <limits>
#include <algorithm>
#include <map>
 
#include "JLang/JException.hh"
#include "JLang/JTitle.hh"
#include "Jeep/JPrint.hh"
 
 
/**
 * \author mdejong
 */
 
namespace JTOOLS {}
namespace JPP { using namespace JTOOLS; }
 
namespace JTOOLS {
 
  using JLANG::JTitle;
  using JLANG::JDivisionByZero;
  using JLANG::JNoValue;
 
 
  /**
   * Quantile calculator.\n
   * This class acts as a zero-dimensional histogram.
   */
  class JQuantile :
    public JTitle
  {
  public:
    /**
     * Constructor.
     *
     * \param  title       title
     * \param  option      option
     */
    JQuantile(const JTitle& title  = "",
              const int     option = 0) :
      JTitle   (title),
      quantiles(false)
    {
      reset();
 
      quantiles = (option > 0);
    }
 
 
    /**
     * Constructor.\n
     * Set range according to minimal and maximal value in input data.
     *
     * \param  title            title
     * \param  __begin          begin of data
     * \param  __end            end   of data
     * \param  member           pointer to data member
     * \param  w                weight
     */
    template<class iterator_type, class value_type, class result_type>
    JQuantile(const JTitle& title,
              iterator_type __begin,
              iterator_type __end,
              result_type value_type::*member,
              const double w = 1.0) :
      JTitle   (title),
      quantiles(true)
    {
      put(__begin, __end, member, w);
    }
 
 
    /**
     * Constructor.\n
     * Set range according to minimal and maximal value in input data.
     *
     * \param  title            title
     * \param  __begin          begin of data
     * \param  __end            end   of data
     * \param  function         pointer to member method
     * \param  w                weight
     */
    template<class iterator_type, class value_type, class result_type>
    JQuantile(const JTitle& title,
              iterator_type __begin,
              iterator_type __end,
              result_type (value_type::*function)() const,
              const double w = 1.0) :
      JTitle   (title),
      quantiles(true)
    {
      put(__begin, __end, function, w);
    }
 
 
    /**
     * Reset.
     */
    void reset() 
    {
      sum = 0.0;
      rms = 0.0;
      tot = 0.0;
      num = 0;
      min = +std::numeric_limits<double>::max();
      max = -std::numeric_limits<double>::max();
 
      buffer.clear();
    }
    
    
    /**
     * Put value.
     *
     * \param  x           value
     * \param  w           weight
     */
    void put(const double x, const double w = 1.0)
    {
      sum += w*x;
      rms += w*x*x;
      tot += w;
      num += 1;
      min  = std::min(min, x);
      max  = std::max(max, x);
 
      if (quantiles) {
        buffer.insert(std::make_pair(x,w));
      }
    }
    
 
    /**
     * Put data.
     *
     * \param  __begin          begin of data
     * \param  __end            end   of data
     * \param  member           pointer to data member
     * \param  w                weight
     */
    template<class iterator_type, class value_type, class result_type>
    void put(iterator_type __begin,
             iterator_type __end,
             result_type value_type::*member,
             const double w = 1.0) 
    {
      for (iterator_type i = __begin; i != __end; ++i) {
        put((*i).*member, w);
      }
    }
 
 
    /**
     * Put data.
     *
     * \param  __begin          begin of data
     * \param  __end            end   of data
     * \param  function         pointer to member method
     * \param  w                weight
     */
    template<class iterator_type, class value_type, class result_type>
    void put(iterator_type __begin,
             iterator_type __end,
             result_type (value_type::*function)() const,
             const double w = 1.0)
    {
      for (iterator_type i = __begin; i != __end; ++i) {
        put(((*i).*function)(), w);
      }
    }
 
 
    /**
     * Get total count.
     *
     * \return             count
     */
    long long int getCount() const
    {
      return num;
    }
    
 
    /**
     * Get total weight.
     *
     * \return             weight
     */
    double getTotal() const
    {
      return tot;
    }
 
 
    
    /**
     * Get weighted sum.
     *
     * \return             weighted sum
     */
    double getSum() const
    {
      return sum;
    }
    
 
    /**
     * Get mean value.
     *
     * \return             mean value
     */
    double getMean() const
    {
      if (tot != 0.0)
        return sum / tot;
      else
        THROW(JDivisionByZero, "JQuantile::getMean()");
    }
 
 
    /**
     * Get RMS.
     *
     * \return             RMS
     */
    double getRMS() const
    {
      if (tot != 0.0)
        return sqrt(rms/tot);
      else
        THROW(JDivisionByZero, "JQuantile::getRMS()");
    }
 
 
    /**
     * Get standard deviation
     *
     * \return             standard deviation
     */
    double getSTDev() const
    {
      if (tot != 0.0)
        return sqrt(rms*tot - sum*sum)/tot;
      else
        THROW(JDivisionByZero, "JQuantile::getSTDev()");
    }
 
 
    /**
     * Get maximal deviation from average.
     *
     * \param  relative    relative to average or absolute
     * \return             deviation
     */
    double getDeviation(const bool relative = true) const
    {
      if (relative) 
        return std::max(max - getMean(), getMean() - min);
      else
        return std::max(fabs(max), fabs(min));
    }
 
 
    /**
     * Test relative accuracy.
     *
     * \param  precision   relative precision
     * \return             true if reached accuracy; else false
     */
    bool hasAccuracy(const double precision) const
    {
      return getCount() > 3 && getSTDev() < precision * getMean();
    }
    
 
    /**
     * Get quantile.
     *
     * \param  Q           quantile 
     * \param  reverse     reverse
     * \return             value
     */
    double getQuantile(const double Q, const bool reverse = false) const
    {
      if (quantiles) {
 
        double W = 0.0;
 
        for (std::map<double, double>::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
          W += i->second;
        }
 
        if (reverse)
          return getQuantile(buffer.rbegin(), buffer.rend(), Q*W);
        else
          return getQuantile(buffer. begin(), buffer. end(), Q*W);
      }
      
      THROW(JNoValue, "JQuantile::getQuantile()");
    }
 
    
    /**
     * Print quantile.
     *
     * \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(10) << left  << "    mean"      << ' '
            << setw(10) << left  << "    STD"       << ' '
            << setw(10) << left  << "  deviation"   << endl;
      }
      
      out << setw(nc) << left << getTitle()          << ' '
          << SCIENTIFIC(10,2) << getMean()           << ' ' 
          << SCIENTIFIC(10,2) << getSTDev()          << ' '
          << SCIENTIFIC(10,2) << getDeviation(false) << endl;
 
      return out;
    }
    
 
    /**
     * Print quantile.
     *
     * \param  out           output stream
     * \param  quantile      quantile
     * \return               output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JQuantile& quantile)
    {
      return quantile.print(out, getLongprint(out));
    }
      
  protected:
    double sum;
    double rms;
    double tot;
    double min;
    double max;
    long long int num;
 
    bool   quantiles;
    std::multimap<double, double> buffer;
 
 
    /**
     * Get quantile.
     *
     * \param  __begin     begin of data
     * \param  __end       end   of data
     * \param  W           weight
     * \return             value
     */
    template<class T>
    static double getQuantile(T __begin, T __end, const double W)
    {
      double w = 0.0;
      
      for (T i = __begin; i != __end; ++i) {
            
        w += i->second;
            
        if (w >= W) {
          return i->first;
        }
      }
 
      THROW(JNoValue, "JQuantile::getQuantile()");
    }
  };
}  
 
#endif