JQuantiles.hh

Go to the documentation of this file.

#ifndef __JTOOLS__JQUANTILES__
#define __JTOOLS__JQUANTILES__
 
#include <limits>
#include <cmath>
 
#include "JLang/JException.hh"
#include "JTools/JRange.hh"
#include "JTools/JElement.hh"
#include "JTools/JResult.hh"
#include "JTools/JAbstractCollection.hh"
#include "JTools/JToolsToolkit.hh"
 
 
/**
 * \author mdejong
 */
 
namespace JTOOLS {}
namespace JPP { using namespace JTOOLS; }
 
namespace JTOOLS {
 
  using JLANG::JException;
  using JLANG::JEmptyCollection;
 
 
  /**
   * Quantile calculator for a given interpolating function.
   * It is assumed that the function has a single maximum.
   */
  class JQuantiles :
    public JRange<double> 
  {
  public:
 
    typedef JAbstractCollection<double>     JAbscissa_t;
 
    /**
     * Default constructor.
     */
    JQuantiles() :
      Xmax(0.0),
      Ymax(0.0),
      fwhm(0.0),
      sum (0.0)
    {}
   
 
    /**
     * Constructor.
     *
     * \param   f1        functional collection
     * \param   Q         quantile
     * \param   eps       relative precision
     */
    template<class JFunction1D_t>
    JQuantiles(const JFunction1D_t& f1,
               const double         Q   = 1.0,
               const double         eps = 1.0e-6) :
      Xmax(0.0),
      Ymax(0.0),
      fwhm(0.0),
      sum (0.0)
    {
      set(f1, Q, eps);
    }
   
 
    /**
     * Constructor.
     *
     * \param   abscissa  abscissa
     * \param   f1        function
     * \param   Q         quantile
     * \param   eps       relative precision
     */
    template<class JFunction1D_t>
    JQuantiles(const JAbscissa_t&   abscissa,
               const JFunction1D_t& f1,
               const double         Q   = 1.0,
               const double         eps = 1.0e-6) :
      Xmax(0.0),
      Ymax(0.0),
      fwhm(0.0),
      sum (0.0)
    {
      set(abscissa, f1, Q, eps);
    }
 
 
    /**
     * Set quantiles.
     *
     * \param   f1        functional collection
     * \param   Q         quantile
     * \param   eps       relative precision
     */
    template<class JFunction1D_t>
    void set(const JFunction1D_t& f1,
             const double         Q   = 1.0,
             const double         eps = 1.0e-6) 
    {
      typedef typename JFunction1D_t::const_iterator  const_iterator;
 
      if (f1.empty()) {
        throw JEmptyCollection("JQuantiles() no data.");
      }
 
 
      // maximum
 
      const_iterator p = f1.begin();
 
      for (const_iterator i = f1.begin(); i != f1.end(); ++i) {
        if (i->getY() > p->getY()) {
          p = i;
        }
      }
      
 
      // x at maximum
      
      Xmax = p->getX();
        
      if (p != f1.begin()) {
          
        const double xa = (--p)->getX();
        const double xb = (++p)->getX();
          
        if (++p != f1.end()) {
            
          const double xc = p->getX();
            
          Xmax = search(xa, xb, xc, f1, -1, eps);
        }
      }
        
      Ymax = get_value(f1(Xmax));
 
 
      // integral & quantile
 
      if (Q > 0.0 && Q <= 1.0) {
 
        JSplineFunction1D<JSplineElement2D<double, double>, JCollection, double> buffer;
 
        try {
 
          sum = makeCDF(f1, buffer);
 
          setLowerLimit(buffer(0.5 * (1.0 - Q)));
          setUpperLimit(buffer(0.5 * (1.0 + Q)));
        }
        catch(const JException& error) {
          sum = 0.0;
        }
 
      } else {
 
        sum = JTOOLS::getIntegral(f1);
 
        if        (Q >  1.0) {
          setLowerLimit(f1. begin()->getX());
          setUpperLimit(f1.rbegin()->getX());
        } else if (Q <= 0.0) {
          setLowerLimit(Xmax);
          setUpperLimit(Xmax);
        }
      }
 
 
      // FWHM
        
      fwhm = 0.0;
        
      for (double xmin = f1.begin()->getX(), xmax = Xmax, v = 0.5*Ymax; ; ) {
        
        const double x = 0.5 * (xmin + xmax);
        const double y = get_value(f1(x));
        
        if (fabs(y - v) < eps*v || xmax - xmin < eps) {
          fwhm -= x;
          break;
        }
        
        if (y > v)
          xmax = x;
        else
          xmin = x;
      }
      
      for (double xmin = Xmax, xmax = f1.rbegin()->getX(), v = 0.5*Ymax; ; ) {
        
        const double x = 0.5 * (xmin + xmax);
        const double y = get_value(f1(x));
          
        if (fabs(y - v) < eps*v || xmax - xmin < eps) {
          fwhm += x;
          break;
        }
        
        if (y > v)
          xmin = x;
        else
          xmax = x;
      }
    }
 
 
    /**
     * Set quantiles.
     *
     * \param   abscissa  abscissa
     * \param   f1        function
     * \param   Q         quantile
     * \param   eps       relative precision
     */
    template<class JFunction1D_t>
    void set(const JAbscissa_t&   abscissa,
             const JFunction1D_t& f1,
             const double         Q   = 1.0,
             const double         eps = 1.0e-6) 
    {
      JSplineFunction1D<JSplineElement2D<double, double>, JCollection, double> buffer;
 
      buffer.configure(abscissa, f1);
      buffer.compile();
 
      set(buffer, Q, eps);
    }
 
 
    /**
     * Get position of maximum.
     *
     * \return     x value at maximum
     */
    double getX() const 
    { 
      return Xmax;
    }
 
 
    /**
     * Get value of maximum.
     *
     * \return     y value at maximum
     */
    double getY() const 
    {
      return Ymax;
    }
 
 
    /**
     * Get Full Width at Half Maximum.
     *
     * \return     FWHM
     */
    double getFWHM() const
    {
      return fwhm;
    }
    
 
    /**
     * Get integral of function.
     *
     * \return     integral
     */
    double getIntegral() const
    {
      return sum;
    }
 
 
  protected:
 
    double Xmax;
    double Ymax;
    double fwhm;
    double sum;
 
 
    /**
     * Locate maximum or minimun of function.
     * 
     * Golden section search code is adopted from reference:
     * Numerical Recipes in C++, W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery,
     * Cambridge University Press.
     * 
     *             xa < xb < xc
     *
     *             is = +1  ->  There is a minimum, i.e:  f(xb) < min(f(xa),f(xc))
     *             is = -1  ->  There is a maximum, i.e:  f(xb) > max(f(xa),f(xc))
     *
     * \param  xa       
     * \param  xb       
     * \param  xc       
     * \param  f        function
     * \param  is       sign (+1 -> minimim, -1 -> maximum)
     * \param  eps      relative precision
     */
    template<class JFunction1D_t>
    static double search(const double xa,
                         const double xb,
                         const double xc,
                         const JFunction1D_t& f,
                         const int    is,
                         const double eps = 1.0e-6)
    {
      static const double R = 0.61803399;
      static const double C = 1.0 - R;
 
      double x0 = xa;
      double x3 = xc;
      double x1, x2;
 
      if (fabs(xc-xb) > fabs(xb-xa)) {
        x1 = xb;
        x2 = xb + C*(xc-xb);
      } else {
        x2 = xb;
        x1 = xb - C*(xb-xa);
      }
 
      double f1 = is * get_value(f(x1));
      double f2 = is * get_value(f(x2));
 
      while (fabs(x3-x0) > eps*(fabs(x1)+fabs(x2))) {
 
        if (f2 < f1) {
 
          x0 = x1;
          x1 = x2;
          x2 = R*x2 + C*x3;
 
          f1 = f2;
          f2 = is * get_value(f(x2));
 
        } else {
 
          x3 = x2;
          x2 = x1;
          x1 = R*x1 + C*x0;
 
          f2 = f1;
          f1 = is * get_value(f(x1));
        }
      }
 
      if (f1 < f2)
        return x1;
      else
        return x2;
    }
  };
}
 
#endif