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;
    }
 
 
    /**
     * 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 i : static_cast<const std::vector<double>&>(*this)) {
          if (-i > x1) {
            x1 = -i;
          }
        }
 
        x1 += EPSILON;                           // offset
 
        f2 = f1;
        f1 = getDerivative(x1);
 
      } else {
 
        return { 0.0, 0.0 };
      }
 
      // 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 };
    }
 
 
    /**
     * 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