JGarage.hh

Go to the documentation of this file.

#ifndef __JASTRONOMY__JGARAGE__
#define __JASTRONOMY__JGARAGE__
 
#include <vector>
#include <numeric>
#include <algorithm>
#include <stdexcept>
 
 
/**
 * \author mdejong
 */
 
namespace JASTRONOMY {}
namespace JPP { using namespace JASTRONOMY; }
 
namespace JASTRONOMY {
 
  /**
   * Auxiliary container for statistical analysis of a large number of positive (or zero) values.
   *
   * The data are organised in groups according to a numerical hash key 
   * so that the computation involving quantiles can be made faster.\n
   * The numerical hash key is obtained by multiplying the input value with the given factor and then taking the unsigned integer value of the result.\n
   * The computed quantiles correspond to a distribution with ascending values.
   */
  template<class T, const size_t factor = 0>
  struct JGarage :
    public std::vector< std::vector<T> >
  {
    typedef std::vector< std::vector<T> >                map_type;
    typedef typename map_type::value_type                value_type;
    typedef typename map_type::const_iterator            const_iterator;
    typedef typename map_type::const_reverse_iterator    const_reverse_iterator;
 
 
    /**
     * Get total number of values.
     *
     * \return            number of values
     */
    size_t getN() const
    {
      using namespace std;
 
      return accumulate(this->begin(), this->end(), 0, [](const size_t n, const value_type& buffer) { return n + buffer.size(); });
    }
 
 
    /**
     * Get minimal value.
     *
     * \return            value
     */
    T getMin() const
    {
      using namespace std;
 
      for (const_iterator i = this->begin(); i != this->end(); ++i) {
        if (!i->empty()) {
          return *min_element(i->begin(), i->end());
        }
      }
 
      throw length_error("No data.");
    }
 
 
    /**
     * Get maximal value.
     *
     * \return            value
     */
    T getMax() const
    {
      using namespace std;
 
      for (const_reverse_iterator i = this->rbegin(); i != this->rend(); ++i) {
        if (!i->empty()) {
          return *max_element(i->begin(), i->end());
        }
      }
 
      throw length_error("No data.");
    }
 
 
    /**
     * Get hash key corresponding to given value.
     *
     * \param  value      value
     * \return            key
     */
    static size_t getKey(const T value)
    {
      if (value > 0.0)
        return value * factor;
      else
        return 0;
    }
 
 
    /**
     * Add value to container.
     *
     * \param  value      value
     */
    void put(const T value)
    {
      const size_t index = getKey(value);
 
      if (index >= this->size()) {
        this->resize(index + 1);
      }
 
      (*this)[index].push_back(value);
    }
 
 
    /**
     * Add values to container.
     *
     * \param  N          number of values
     * \param  object     function object returning value
     */
    template<class JFunction_t>
    void put(const size_t N, const JFunction_t& object)
    {
      for (size_t i = 0; i != N; ++i) {
        this->put(object());
      }
    }
 
 
    /**
     * For-each method.
     *
     * \param  object     function object
     */
    template<class JFunction_t>
    inline void operator()(const JFunction_t& object) const
    {
      for (const_iterator i = this->begin(); i != this->end(); ++i) {
        for (const T x : *i) {
          object(x);
        }
      }
    }
 
 
    /**
     * Get minimal value corresponding given maximal probability.
     *
     * \param  P          probability [0,1>
     * \return            value
     */
    T getValue(const double P) const
    {
      using namespace std;
 
      if (P < 0.0 || P > 1.0) {
        throw out_of_range("Invalid probability.");
      }
 
      const size_t N = this->getN();
      const size_t M = N - (P * N);
 
      if (N == 0) {
        throw length_error("No data.");
      }
 
      if (M == 0) { return getMin(); }
      if (M == N) { return getMax(); }
 
      size_t n = N;
 
      const_reverse_iterator q = this->rbegin();
 
      while ((n -= q->size()) > M) {
        ++q;
      }
 
      value_type buffer(q->begin(), q->end());
 
      typename value_type::iterator p = next(buffer.begin(), M - n);
 
      nth_element(buffer.begin(), p, buffer.end());
 
      return *p;
    }
   
 
    /**
     * Get maximal probability corresponding given minimal value.
     *
     * \param  value      value
     * \return            probability
     */
    double getProbability(const T value) const
    {
      using namespace std;
 
      const size_t N = this->getN();
      const int    M = getKey(value);
 
      if (N == 0) {
        throw length_error("No data.");
      }
 
      size_t n = 0;
      int    i = this->size() - 1;
 
      for ( ; i != -1 && i > M; --i) {
        n += (*this)[i].size();
      }
        
      if (i != -1) {
        for (const T x : (*this)[i]) {
          if (x >= value) {
            n += 1;
          }
        }
      }
 
      return (double) n / (double) N;
    }
  };
}
 
#endif