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;


  /**
   * Locate maximum or minimun of function.
   * 
   * Golden section search code is taken 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>
  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;
  }


  /**
   * Quantile calculator for a given 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;
  };
}

#endif