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 "JLang/JVectorize.hh"
#include "JLang/JManip.hh"
 
#include "JMath/JMath.hh"
 
#include "JTools/JStats.hh"
 
 
/**
 * \author mdejong
 */
 
namespace JTOOLS {}
namespace JPP { using namespace JTOOLS; }
 
namespace JTOOLS {
 
  using JLANG::JNoValue;
  using JLANG::array_type;
  using JLANG::make_array;
 
 
  /**
   * Auxiliary data structure for running average, standard deviation and quantiles.
   */
  struct JQuantile :
    public JStats,
    public JMath<JQuantile>
  {
    /**
     * Constructor.
     *
     * \param  title            title
     */
    JQuantile(const JTitle& title     = "") :
      JStats(title)
    {
      reset();
    }
 
 
    /**
     * Constructor.
     *
     * \param  title            title
     * \param  buffer           input data
     * \param  w                weight
     */
    template<class JElement_t, class JAllocator_t>
    JQuantile(const JTitle& title,
              const array_type<JElement_t, JAllocator_t>& buffer,
              const double  w         = 1.0) :
      JStats(title)
    {
      reset();
 
      put(buffer, w);
    }
 
 
    /**
     * Reset.
     */
    void reset()
    {
      JStats::reset();
 
      buffer.clear();
    }
 
 
    /**
     * Add quantile.
     *
     * \param  Q           quantile
     * \return             this quantile
     */
    JQuantile& add(const JQuantile& Q)
    {
      JStats::add(Q);
 
      std::copy(Q.buffer.begin(), Q.buffer.end(), std::inserter(buffer, buffer.end()));
 
      return *this;
    }
 
 
    /**
     * Put value.
     *
     * \param  x           value
     * \param  w           weight
     */
    void put(const double x, const double w = 1.0)
    {
      JStats::put(x, w);
 
      buffer.insert(std::make_pair(x,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);
      }
    }
 
 
    /**
     * Options for evaluation of quantile.
     */
    enum Quantile_t {
      forward_t    = +1,     //!< forward
      symmetric_t  =  0,     //!< symmatric
      backward_t   = -1      //!< backward
    };
 
    
    /**
     * Get quantile.
     *
     * \param  Q           quantile 
     * \param  option      option
     * \return             value
     */
    double getQuantile(const double Q, const Quantile_t option = forward_t) const
    {
      double W = 0.0;
 
      for (std::map<double, double>::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
        W += i->second;
      }
 
      switch (option) {
      case forward_t:
        return (getQuantile(buffer. begin(), buffer. end(), Q*W));
      case symmetric_t:
        return (getQuantile(buffer.rbegin(), buffer.rend(), 0.5*(1.0 - Q)*W)  -
                getQuantile(buffer. begin(), buffer. end(), 0.5*(1.0 - Q)*W));
      case backward_t:
        return (getQuantile(buffer.rbegin(), buffer.rend(), Q*W));
      default:
        THROW(JNoValue, "Invalid option " << option);
      }
    }
 
  protected:
    /**
     * 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, "Invalid weight " << W);
    }
 
    std::multimap<double, double> buffer;
  };
}
 
#endif