JAspera.hh

Go to the documentation of this file.

#ifndef __JASTRONOMY__JASPERA__
#define __JASTRONOMY__JASPERA__
 
#include <vector>
#include <cmath>
 
 
/**
 * \file
 *
 * Per aspera ad astra.
 * \author mdejong
 */
namespace JASTRONOMY {}
namespace JPP { using namespace JASTRONOMY; }
 
namespace JASTRONOMY {
 
  /**
   * Auxiliary data structure to fit signal strength using likelihood ratio.
   */
  struct JAspera :
    public std::vector<double>                      // N/S values
  {
 
    static constexpr double EPSILON = 1.0e-3;       //!< precision determination of signal strength
 
 
    /**
     * Result of fit.
     */
    struct fit_type {
      double likelihood  =  0.0;                    //!< likelihood
      double signal      =  0.0;                    //!< signal strength
    };
 
 
    /**
     * Put signal and background to list of pre-computed N/S values.
     *
     * \param  s                signal
     * \param  b                background
     */
    void put(const double s,
             const double b)
    {
      push_back(b/s);
 
      ws += s;
    }
 
 
    /**
     * Put signal and background to list of pre-computed N/S values.
     *
     * \param  n                data
     * \param  s                signal
     * \param  b                background
     */
    void put(const size_t n,
             const double s,
             const double b)
    {
      for (size_t i = 0; i != n; ++i) {
        push_back(b/s);
      }
 
      ws += s;
    }
 
 
    /**
     * Get likelihood for given signal strength.
     *
     * \param  p                signal strength
     * \return                  likelihood
     */
    double getLikelihood(const double p) const
    {
      double y = -p * ws;
 
      for (const double i : static_cast<const std::vector<double>&>(*this)) {
        y += log1p(p/i);
      }
 
      return y;
    }
 
 
    /**
     * Get derivative of likelihood for given signal strength.
     *
     * \param  p                signal strength
     * \return                  derivative of likelihood
     */
    double getDerivative(const double p) const
    {
      double y = -ws;
 
      for (const double i : static_cast<const std::vector<double>&>(*this)) {
        y += 1.0 / (p + i);
      }
 
      return y;
    }
 
 
#if defined(NEWTON_RAPHSON)
 
    /**
     * Auxiliary data structure for 1st and 2nd derivatives.
     */
    struct derivatives_type {
      double fp  = 0.0;
      double fpp = 0.0;
    };
 
 
    /**
     * Get derivatives of likelihood for given signal strength.
     *
     * \param  p                signal strength
     * \return                  derivatives of likelihood
     */
    derivatives_type getDerivatives(const double p) const
    {
      derivatives_type y;
 
      y.fp = -ws;
 
      for (const double x : static_cast<const std::vector<double>&>(*this)) {
 
        const double z = 1.0 / (p + x);
 
        y.fp  += z;
        y.fpp -= z*z;
      }
 
      return y;
    }
#endif
 
    /**
     * Fit signal strength.
     *
     * \param  ns               allow for negative signal
     * \return                  result
     */
    fit_type operator()(const bool ns = false) const
    {
      using namespace std;
 
      if (this->empty()) {
 
        // nothing to be done
 
        return { 0.0, 0.0 };
 
      } else if (this->size() == 1 ) {
 
        // analytical solution
 
        const double x = 1.0/ws - (*this)[0];
 
        if (x > 0.0 || ns)
          return { getLikelihood(x), x };
        else
          return { 0.0, 0.0 };
      }
 
      double x1 = 0.0;
      double x2 = 0.0;
 
      double f1 = getDerivative(0.0);           // discriminator between positive and negative signal
      double f2 = 0.0;
 
      if        (f1 == 0.0) {
 
        return { 0.0, 0.0 };
 
      } else if (f1 >  0.0) {                   // positive signal
 
        x1 = 0.0;                               // lower limit corresponds to no signal
        x2 = (double) this->size() / ws;        // upper limit corresponds to no background (i.e. all N/S = 0)
 
        f2 = getDerivative(x2);
 
      } else if (ns) {                          // negative signal
 
        x2 = 0.0;                               // upper limit corresponds to no signal
        x1 = -(*this)[0];                       // lower limit corresponds to largest negated N/S ratio
 
        for (const double x : static_cast<const std::vector<double>&>(*this)) {
          if (-x > x1) {
            x1 = -x;
          }
        }
 
        x1 += 1.0 / ws;                         // offset
 
        f2 = f1;
        f1 = getDerivative(x1);
 
      } else {
 
        return { 0.0, 0.0 };
      }
 
#if   defined(BISECTION)
 
      // binary search
 
      while (x2 - x1 > EPSILON) {
 
        const double x = 0.5 * (x1 + x2);
 
        if (getDerivative(x) > 0.0)
          x1 = x;
        else
          x2 = x;
      }
 
      const double x = 0.5 * (x1 + x2);
 
      return { getLikelihood(x), x };
 
#elif defined(NEWTON_RAPHSON)
 
      // Newton-Raphson method
 
      double dx = x2 - x1;
      double x  = 0.5 * (x1 + x2);
 
      derivatives_type y = getDerivatives(x);
 
      struct {
        double dx;
      } old;
 
      old.dx = dx;
 
      for ( ; ; ) {
 
        if ((y.fp + (x2 - x)*y.fpp)*(y.fp + (x1 - x)*y.fpp) > 0.0 || fabs(2.0*y.fp) > fabs(old.dx*y.fpp)) {
 
          old.dx = dx;
 
          dx  = 0.5 * (x2 - x1);
          x   = x1 + dx;
 
        } else {
 
          old.dx = dx;
 
          dx  = y.fp/y.fpp;
          x  -= dx;
        }
 
        if (fabs(dx) <= EPSILON) {
          break;
        }
 
        y = getDerivatives(x);
 
        if (y.fp < 0.0)
          x2 = x;
        else
          x1 = x;
      }
 
      return { getLikelihood(x), x };
 
#else
 
      // Ridder's method
 
      while (x2 - x1 > EPSILON) {
 
        const double xm = 0.5 * (x1 + x2);
        const double fm = getDerivative(xm);
 
        const double s  = sqrt(fm*fm - f1*f2);
 
        if (s == 0.0) {
          break;
        }
 
        const double xn = xm + (xm - x1) * fm/s;
        const double fn = getDerivative(xn);
 
        if (fn == 0.0) {
          return { getLikelihood(xn), xn };
        }
 
        if (signbit(fn) != signbit(fm)) {
 
          x1 = xm;
          f1 = fm;
          x2 = xn;
          f2 = fn;
 
        } else {
 
          if (signbit(fn)) {
 
            x2 = xn;
            f2 = fn;
 
          } else {
 
            x1 = xn;
            f1 = fn;
          }
        }
      }
 
      const double x = 0.5 * (x1 + x2);
 
      return { getLikelihood(x), x };
#endif
    }
 
 
    /**
     * Get test statistic for given signal strength.
     * See formula 16 in this <a href="https://link.springer.com/article/10.1140/epjc/s10052-011-1554-0">reference</a>.
     *
     * \param  ps               signal strength
     * \return                  test statistic
     */
    inline double getTestStatisticForUpperLimit(const double ps) const
    {
      const fit_type result = (*this)(true);
 
      if      (result.signal <= 0.0)
        return 0.0                -  this->getLikelihood(ps);
      else if (result.signal <= ps)
        return result.likelihood  -  this->getLikelihood(ps);
      else
        return 0.0;
    }
 
 
    /**
     * Get total signal strength.
     *
     * \return           signal strength
     */
    double getSignal() const
    {
      return ws;
    }
 
 
    /**
     * Set signal strength.
     *
     * \param  wS        signal strength
     */
    void setSignal(const double wS)
    {
      ws = wS;
    }
 
 
    /**
     * Add signal strength.
     *
     * \param  wS        signal strength
     */
    void addSignal(const double wS)
    {
      ws += wS;
    }
 
  protected:
 
    double ws  =  0.0;                                     //!< total signal strength
  };
}
 
#endif