JAstronomyToolkit.hh

Go to the documentation of this file.

#ifndef __JASTRONOMY__JASTRONOMYTOOLKIT__
#define __JASTRONOMY__JASTRONOMYTOOLKIT__
 
#include <vector>
#include <map>
#include <cmath>
#include <numeric>
#include <limits>
#include <algorithm>
 
#include "JLang/JException.hh"
 
 
/**
 * \author mdejong
 */
 
namespace JASTRONOMY {}
namespace JPP { using namespace JASTRONOMY; }
 
namespace JASTRONOMY {
 
  using JLANG::JValueOutOfRange;
  using JLANG::JEmptyCollection;
 
 
  /**
   * Auxiliary container for statistical analysis of a large number of values.
   *
   * The data are organised in groups according to a numerical key 
   * so that the computation involving quantiles can be speed up.\n
   * The numerical key is obtained by scaling the value with a given factor and then taking the floor of the result.\n
   * The quantiles correspond to a distribution with ascending values.
   */
  template<class T>
  struct JQuantile :
    public std::map<int, std::vector<T> >
  {
    typedef std::map<int, 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;
 
 
    /**
     * Constructor.
     *
     * \param  factor     factor
     */
    JQuantile(const double factor) :
      factor(factor)
    {}
    
 
    /**
     * Get total number of values.
     *
     * \return            number of values
     */
    size_t getN() const
    {
      return std::accumulate(this->begin(), this->end(), 0, [](const size_t n, const value_type& i) { return n + i.second.size(); });
    }
 
 
    /**
     * Get minimal value.
     *
     * \return            value
     */
    T getMin() const
    {
      for (const_iterator i = this->begin(); i != this->end(); ++i) {
        if (!i->second.empty()) {
          return *std::min_element(i->second.begin(), i->second.end());
        }
      }
 
      THROW(JEmptyCollection, "Missing data.");
    }
 
 
    /**
     * Get maximal value.
     *
     * \return            value
     */
    T getMax() const
    {
      for (const_reverse_iterator i = this->rbegin(); i != this->rend(); ++i) {
        if (!i->second.empty()) {
          return *std::max_element(i->second.begin(), i->second.end());
        }
      }
 
      THROW(JEmptyCollection, "Missing data.");
    }
 
 
    /**
     * Get key of given value.
     *
     * \param  value      value
     * \return            key
     */
    int getKey(const T value) const
    {
      return floor(value * factor);
    }
 
 
    /**
     * Add value to container.
     *
     * \param  value      value
     */
    void put(const T value)
    {
      (*this)[getKey(value)].push_back(value);
    }
 
 
    /**
     * 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(JValueOutOfRange, "Invalid probability " << P);
      }
 
      const size_t N = this->getN();
      const size_t M = N - (size_t) (P * N);
 
      if (N == 0) {
        THROW(JEmptyCollection, "Missing data.");
      }
 
      if (M == N) { return this->getMax(); }
      if (M == 0) { return this->getMin(); }
 
      size_t n = N;
        
      for (const_reverse_iterator p = this->rbegin(); p != this->rend(); ++p) {
 
        n -= p->second.size();
 
        if (n <= M) {
 
          vector<T> buffer(p->second.begin(), p->second.end());
 
          sort(buffer.begin(), buffer.end());
 
          return buffer[M - n];
        }
      }
 
      THROW(JValueOutOfRange, "Invalid index " << M << ' ' << N);
    }
   
 
    /**
     * Get maximal probability corresponding given minimal value.
     *
     * \param  value      value
     * \return            probability
     */
    double getProbability(const double value) const
    {
      const int M = getKey(value);
      size_t    n = 0;
 
      const_reverse_iterator i = this->rbegin();
 
      for ( ; i != this->rend() && i->first > M; ++i) {
        n += i->second.size();
      }
        
      if (i != this->rend()) {
        for (T x : i->second) {
          if (x >= value) {
            n += 1;
          }
        }
      }
 
      return (double) n / (double) this->getN();
    }
 
 
    const T factor;      //!< scaling factor for rounding value
  };
}
 
#endif