JKatoomba_t.hh

Go to the documentation of this file.

#ifndef __JACOUSTICS__JKATOOMBA_T__
#define __JACOUSTICS__JKATOOMBA_T__

#include <vector>
#include <memory>
#include <mutex>

#include "TMatrixTSym.h"
#include "TDecompSVD.h"

#include "JFit/JEstimator.hh"
#include "JFit/JRegressor.hh"
#include "JFit/JSimplex.hh"
#include "JFit/JMEstimator.hh"
#include "JMath/JMatrixNS.hh"
#include "JLang/JType.hh"
#include "JLang/JSharedPointer.hh"
#include "JLang/JException.hh"
#include "JMath/JMatrixNS.hh"
#include "Jeep/JMessage.hh"

#include "JAcoustics/JTransmission.hh"
#include "JAcoustics/JSoundVelocity.hh"
#include "JAcoustics/JEKey.hh"
#include "JAcoustics/JGeometry.hh"
#include "JAcoustics/JModel.hh"
#include "JAcoustics/JHit.hh"
#include "JAcoustics/JFitParameters.hh"


/**
 * \file
 *
 * Fit functions of acoustic model.
 * \author mdejong
 */
namespace JACOUSTICS {}
namespace JPP { using namespace JACOUSTICS; }

namespace JACOUSTICS {

  using JLANG::JType;
  using JLANG::JValueOutOfRange;
  using JMATH::JMath;
  using JFIT::JEstimator;
  using JFIT::JAbstractMinimiser;
  using JFIT::JSimplex;
  using JFIT::JMEstimator;


  /**
   * Get total weight of data points.
   *
   * \param  __begin        begin of data
   * \param  __end          end   of data
   * \return                total weight
   */
  template<class T>
  inline double getWeight(T __begin, T __end)
  {
    double W = 0;

    for (T i = __begin; i != __end; ++i) {
      W += i->getWeight();
    }

    return W;
  }


  /**
   * Template definition of fit function of acoustic model.
   */
  template<template<class T> class JMinimiser_t = JType>
  struct JKatoomba;


  /**
   * Auxiliary base class for fit function of acoustic model.
   */
  template<>
  struct JKatoomba<JType> {

    typedef std::shared_ptr<JMEstimator>  estimator_type;


    /**
     * Auxiliary data structure to sort transmissions.
     */
    static const struct compare {
      /**
       * Sort transmissions in following order.
       *  - ascending  identifier;
       *  - descending quality;
       *  - ascending  time-of-arrival;
       *
       * \param  first              first  transmission
       * \param  second             second transmission
       * \return                    true if first transmission before second; else false
       */
      inline bool operator()(const JTransmission& first, const JTransmission& second) const
      {
        if (first.getID() == second.getID()) {

          if (first.getQ() == second.getQ())
            return first.getToA() < second.getToA();
          else
            return first.getQ()   > second.getQ();

        } else {

          return first.getID() < second.getID();
        }
      }
    } compare;


    /**
     * Constructor
     *
     * \param  geometry           detector geometry
     * \param  velocity           sound velocity
     * \param  option             option
     */
    JKatoomba(const JGeometry&      geometry,
              const JSoundVelocity& velocity,
              const int             option) :
      geometry(geometry),
      velocity(velocity),
      option  (option)
    {};
 

    /**
     * Get estimated time-of-arrival for given hit.
     *
     * \param  model              model
     * \param  hit                hit
     * \return                    time-of-arrival
     */
    double getToA(const JModel& model, const JHit& hit) const
    {
      const JGEOMETRY::JString& string     = geometry[hit.getString()];
      const JMODEL   ::JString& parameters = model.string[hit.getString()];
      const JPosition3D         position   = string.getPosition(parameters, hit.getFloor());

      const double D     = hit.getDistance(position);
      const double Vi    = velocity.getInverseVelocity(D, hit.getZ(), position.getZ());

      return model.emission[hit.getEKey()].t1  +  D * Vi;
    }
 

    /**
     * Get estimated time-of-emission for given hit.
     *
     * \param  model              model
     * \param  hit                hit
     * \return                    time-of-arrival
     */
    double getToE(const JModel& model, const JHit& hit) const
    {
      const JGEOMETRY::JString& string     = geometry[hit.getString()];
      const JMODEL   ::JString& parameters = model.string[hit.getString()];
      const JPosition3D         position   = string.getPosition(parameters, hit.getFloor());

      const double D     = hit.getDistance(position);
      const double Vi    = velocity.getInverseVelocity(D, hit.getZ(), position.getZ());

      return hit.getValue()  -  D * Vi;
    }


    const JGeometry& geometry;                     //!< detector geometry
    JSoundVelocity   velocity;                     //!< sound velocity
    int              option;                       //!< fit option
    estimator_type   estimator;                    //!< M-Estimator function

  protected:
    /**
     * H-equation as per hit.
     */
    struct H_t :
      public JMODEL::JEmission,
      public JMODEL::JString,
      public JMath<H_t>
    {
      /**
       * Default constructor.
       */
      H_t() :
        JMODEL::JEmission(),
        JMODEL::JString ()
      {}


      /**
       * Constructor.
       *
       * \param  emission      emission
       * \param  string        string
       */
      H_t(const JMODEL::JEmission& emission,
          const JMODEL::JString&   string) :
        JMODEL::JEmission(emission),
        JMODEL::JString  (string)
      {}


      /**
       * Scale H-equation.
       *
       * \param  factor        multiplication factor
       * \return               this H-equation
       */
      H_t& mul(const double factor)
      {
        static_cast<JMODEL::JEmission&>(*this).mul(factor);
        static_cast<JMODEL::JString&>  (*this).mul(factor);

        return *this;
      }
    };


    /**
     * Indices of H-equation in global model.
     */
    struct I_t 
    {
      /**
       * Default constructor.
       */
      I_t() :
        t1 (0),
        tx (0),
        ty (0),
        tx2(0),
        ty2(0),
        vs (0)
      {}


      int t1;
      int tx;
      int ty;
      int tx2;
      int ty2;
      int vs;
    };
  };


  /**
   * Template specialisation of fit function of acoustic model based on JAbstractMinimiser minimiser.\n
   * This class can be used to evaluate the chi2.
   */
  template<>
  struct JKatoomba<JAbstractMinimiser> : 
    public JAbstractMinimiser<JModel>,
    public JKatoomba<>
  {
    typedef double                                        result_type;


    /**
     * Constructor
     *
     * \param  geometry           detector geometry
     * \param  velocity           sound velocity
     * \param  option             option
     */
    JKatoomba(const JGeometry&      geometry,
              const JSoundVelocity& velocity,
              const int             option) :
      JKatoomba<>(geometry, velocity, option)
    {};
 

    /**
     * Fit function.\n
     * This method is used to determine the chi2 of given hit with respect to actual model.
     *
     * \param  model              model
     * \param  hit                hit
     * \return                    chi2
     */
    result_type operator()(const JModel& model, const JHit& hit) const
    {
      const double toa_s = this->getToA(model, hit);
      const double u     = (toa_s - hit.getValue()) / hit.getSigma();

      return estimator->getRho(u) * hit.getWeight();
    }


    /**
     * Fit.
     *
     * \param  __begin            begin of hits
     * \param  __end              end   of hits
     * \return                    chi2
     */
    template<class T>
    result_type operator()(T __begin, T __end)
    {
      this->value.setOption(this->option);

      return JAbstractMinimiser<JModel>::operator()(*this, __begin, __end);
    }


    /**
     * Fit.
     *
     * \param  model              model
     * \param  __begin            begin of hits
     * \param  __end              end   of hits
     * \return                    chi2
     */
    template<class T>
    result_type operator()(const JModel& model, T __begin, T __end)
    {
      this->value = model;

      return (*this)(__begin, __end);
    }
  };


  /**
   * Algorithm for matrix inversion.
   */
  enum JMatrix_t {
    SVD_t   =   1,      //!< SVD
    LDU_t   =   2       //!< LDU
  };


  /**
   * Auxiliary data structure for compatibility of symmetric matrix.
   * Inversion is based on SVD algorithm.
   */
  struct TMatrixDS :
    public TMatrixD
  {
    TMatrixDS(int N) :
      TMatrixD(N, N)
    {}
  };


  /**
   * Using declaration for compatibility of symmetric matrix.
   * Inversion is based on LDU algorithm.
   */
  using JMATH::JMatrixNS;


  /**
   * Template specialisation of fit function of acoustic model based on linear approximation.
   */
  template<>
  struct JKatoomba<JEstimator> : 
    public JKatoomba<>
  {
    static constexpr double TOLERANCE = 1.0e-8;       //!< Tolerance for SVD

    /**
     * Constructor
     *
     * \param  geometry           detector geometry
     * \param  velocity           sound velocity
     * \param  option             option
     */
    JKatoomba(const JGeometry&      geometry,
              const JSoundVelocity& velocity,
              const int             option) :
      JKatoomba<>(geometry, velocity, option)
    {};


    /**
     * Get start values of string parameters.\n
     * Note that this method may throw an exception.
     *
     * \param  __begin            begin of hits
     * \param  __end              end   of hits
     * \return                    model
     */
    template<class T>
    JModel& operator()(T __begin, T __end) const
    {
      switch (MATRIX_INVERSION) {

      case SVD_t:
        return this->evaluate(__begin, __end, JType<TMatrixDS>());

      case LDU_t:
        return this->evaluate(__begin, __end, JType<JMatrixNS>());

      default:
        THROW(JValueOutOfRange, "Invalid matrix type " << MATRIX_INVERSION); 
      }
    }


    static JMatrix_t  MATRIX_INVERSION;                 // matrix inversion

  private:
    /**
     * Get start values of string parameters.\n
     * Note that this method may throw an exception.
     *
     * \param  __begin            begin of hits
     * \param  __end              end   of hits
     * \param  type               matrix type
     * \return                    model
     */
    template<class T, class JMatrixNS_t>
    JModel& evaluate(T __begin, T __end, const JType<JMatrixNS_t>& type) const
    {
      using namespace std;
      using namespace JPP;
      using namespace JGEOMETRY;

      value = move(JModel(__begin, __end));           // set up global model according data

      if (this->option == JMODEL::fit_emitters_only_t)
        value.setOption(JMODEL::fit_emitters_only_t);
      else
        value.setOption(JMODEL::fit_emitters_and_strings_1st_order_t);

      H_t H;                                          // H-equation as per hit
      I_t i;                                          // indices of H-equation in global model 

      const size_t N = value.getN();

      JMatrixNS_t    V(N);                            // V = H^T x H
      vector<double> Y(N, 0.0);                       // Y = H^T x y

      for (T hit = __begin; hit != __end; ++hit) {

        const JString&    string   = geometry[hit->getString()];
        const JPosition3D position = string.getPosition(hit->getFloor());
        const double      Vi       = velocity.getInverseVelocity(hit->getDistance(position), hit->getZ(), position.getZ());

        const double      h1 = string.getHeight(hit->getFloor());
        const JPosition3D p1 = string.getPosition()  -  hit->getPosition();
        const double      ds = sqrt(p1.getLengthSquared()  +  h1*h1  +  2.0*p1.getZ()*h1);
        const double      y  = hit->getValue()  -  Vi*ds;
        const double      W  = sqrt(hit->getWeight());

        H.t1  =  W * 1.0;
        H.tx  =  W * Vi * p1.getX() * h1 / ds;
        H.ty  =  W * Vi * p1.getY() * h1 / ds;
        
        i.t1  =  value.getIndex(hit->getEKey(),   &H_t::t1);
        i.tx  =  value.getIndex(hit->getString(), &H_t::tx);
        i.ty  =  value.getIndex(hit->getString(), &H_t::ty);

        V(i.t1, i.t1) += H.t1 * H.t1;

        Y[i.t1] += W * H.t1 * y;

        if (hit->getFloor() != 0) {

          if (this->option != JMODEL::fit_emitters_only_t) {
          
            V(i.t1, i.tx) += H.t1 * H.tx;  V(i.t1, i.ty) += H.t1 * H.ty;
            V(i.tx, i.t1) += H.tx * H.t1;  V(i.ty, i.t1) += H.ty * H.t1;

            V(i.tx, i.tx) += H.tx * H.tx;  V(i.tx, i.ty) += H.tx * H.ty;
            V(i.ty, i.tx) += H.ty * H.tx;  V(i.ty, i.ty) += H.ty * H.ty;

            Y[i.tx] += W * H.tx * y;
            Y[i.ty] += W * H.ty * y;
          }
        }
      }

      // evaluate (H^T x H)^-1 x H^T x y 

      invert(V);
      
      for (size_t row = 0; row != N; ++row) {
        for (size_t col = 0; col != N; ++col) {
          value[row] += V(row,col) * Y[col];
        }
      }
      
      return value;
    }


    /**
     * Invert matrix using SVD algorithm.
     *
     * \param  V                  matrix
     */
    static void invert(TMatrixD& V)
    {
#ifdef THREAD_SAFE
      using namespace std;

      unique_lock<mutex> lock(mtx);
#endif
      TDecompSVD svd(V, TOLERANCE);

      Bool_t status;

      V = svd.Invert(status);
    }
    

    /**
     * Invert matrix using LDU algorithm.
     *
     * \param  V                  matrix
     */
    static void invert(JMatrixNS& V)
    {
      V.invert();
    }


    mutable JModel value;
    static std::mutex mtx;                              // TDecompSVD
  };


  /**
   * Template specialisation of fit function of acoustic model based on JSimplex minimiser.
   */
  template<>
  struct JKatoomba<JSimplex> : 
    public JSimplex<JModel>,
    public JKatoomba<>
  {
    /**
     * Constructor
     *
     * \param  geometry           detector geometry
     * \param  velocity           sound velocity
     * \param  option             option
     */
    JKatoomba(const JGeometry&      geometry,
              const JSoundVelocity& velocity,
              const int             option) :
      JKatoomba<>(geometry, velocity, option)
    {};


    /**
     * Fit function.\n
     * This method is used to determine the chi2 of given hit with respect to actual model.
     *
     * \param  model              model
     * \param  hit                hit
     * \return                    chi2
     */
    double operator()(const JModel& model, const JHit& hit) const
    {
      const double toa_s = this->getToA(model, hit);
      const double u     = (toa_s - hit.getValue()) / hit.getSigma();

      return estimator->getRho(u);
    }


    /**
     * Fit.
     *
     * \param  __begin            begin of hits
     * \param  __end              end   of hits
     * \return                    chi2
     */
    template<class T>
    double operator()(T __begin, T __end)
    {
      this->step.clear();

      for (JModel::string_type::const_iterator i = this->value.string.begin(); i != this->value.string.end(); ++i) {

        JModel model;

        switch (this->option) {

        case JMODEL::fit_emitters_and_strings_2nd_order_and_stretching_t:
          model.string[i->first] = JMODEL::JString(0.0, 0.0, 0.0, 0.0, 1.0e-6);  this->step.push_back(model);

        case JMODEL::fit_emitters_and_strings_2nd_order_t:
          model.string[i->first] = JMODEL::JString(0.0, 0.0, 1.0e-6, 0.0, 0.0);  this->step.push_back(model);
          model.string[i->first] = JMODEL::JString(0.0, 0.0, 0.0, 1.0e-6, 0.0);  this->step.push_back(model);

        case JMODEL::fit_emitters_and_strings_1st_order_t:
          model.string[i->first] = JMODEL::JString(1.0e-3, 0.0, 0.0, 0.0, 0.0);  this->step.push_back(model);
          model.string[i->first] = JMODEL::JString(0.0, 1.0e-3, 0.0, 0.0, 0.0);  this->step.push_back(model);
          break;
          
        default:
          break;
        }
      }
      
      for (JModel::emission_type::const_iterator i = this->value.emission.begin(); i != this->value.emission.end(); ++i) {

        JModel model;

        model.emission[i->first] = JMODEL::JEmission(5.0e-5);                    this->step.push_back(model);
      }

      return JSimplex<JModel>::operator()(*this, __begin, __end);
    }
  };


  /**
   * Place holder for custom implementation.
   */
  template<class T>
  class JGandalf {};


  /**
   * Template specialisation of fit function of acoustic model based on JGandalf minimiser.
   */
  template<>
  struct JKatoomba<JGandalf> : 
    public JKatoomba<>
  {
    typedef double result_type;


    /**
     * Constructor
     *
     * \param  geometry           detector geometry
     * \param  velocity           sound velocity
     * \param  option             option
     */
    JKatoomba(const JGeometry&      geometry,
              const JSoundVelocity& velocity,
              const int             option) :
      JKatoomba<>(geometry, velocity, option)
    {};


    /**
     * Fit.
     *
     * \param  __begin            begin of hits
     * \param  __end              end   of hits
     * \return                    chi2 and gradient
     */
    template<class T>
    result_type operator()(T __begin, T __end)
    {
      using namespace std;
      using namespace JPP;

      value.setOption(this->option);

      const int N = value.getN();

      V.resize(N);
      Y.resize(N);
      h.resize(N);

      lambda    = LAMBDA_MIN;

      result_type precessor = numeric_limits<double>::max();

      for (numberOfIterations = 0; numberOfIterations != MAXIMUM_ITERATIONS; ++numberOfIterations) {

        DEBUG("step:     " << numberOfIterations << endl);

        evaluate(__begin, __end);

        DEBUG("lambda:   " << FIXED(12,5) << lambda    << endl);
        DEBUG("chi2:     " << FIXED(12,5) << successor << endl);

        if (successor < precessor) {

          if (numberOfIterations != 0) {

            if (fabs(precessor - successor) < EPSILON*fabs(precessor)) {
              return successor;
            }

            if (lambda > LAMBDA_MIN) {
              lambda /= LAMBDA_DOWN;
            }
          }

          precessor = successor;
          previous  = value;

        } else {

          value   = previous;
          lambda *= LAMBDA_UP;

          if (lambda > LAMBDA_MAX) {
            return precessor;                                   // no improvement found
          }

          evaluate(__begin, __end);
        }


        // force definite positiveness

        for (int i = 0; i != N; ++i) {

          if (V(i,i) < PIVOT) {
            V(i,i) = PIVOT;
          }

          h[i] = 1.0 / sqrt(V(i,i));
        }


        // normalisation

        for (int i = 0; i != N; ++i) {
          for (int j = 0; j != i; ++j) {
            V(j,i) *= h[i] * h[j];
            V(i,j)  = V(j,i);
          }
        }

        for (int i = 0; i != N; ++i) {
          V(i,i) = 1.0 + lambda;
        }


        try {
          V.invert();
        }
        catch (const JException& error) {
          ERROR("JKatoomb<JGandalf>: " << error.what() << endl);
          return precessor;
        }


        for (int i = 0; i != N; ++i) {

          DEBUG("u[" << noshowpos << setw(3) << i << "] = " << showpos << FIXED(20,5) << value[i]);

          double y = 0.0;

          for (int j = 0; j != N; ++j) {
            y += V(i,j) * Y[j] * h[i] * h[j];
          }

          value[i] -= y;

          DEBUG(' ' << FIXED(20,10) << y << noshowpos << endl);
        }
      }

      return precessor;
    }

    static int     debug;                                      //!< debug level
    static int     MAXIMUM_ITERATIONS;                         //!< maximal number of iterations
    static double  EPSILON;                                    //!< maximal distance to minimum
    static double  LAMBDA_MIN;                                 //!< minimal value control parameter
    static double  LAMBDA_MAX;                                 //!< maximal value control parameter
    static double  LAMBDA_UP;                                  //!< multiplication factor control parameter
    static double  LAMBDA_DOWN;                                //!< multiplication factor control parameter
    static double  PIVOT;                                      //!< minimal value diagonal element of matrix

    double                      lambda;
    JModel                      value;
    int                         numberOfIterations;
    JMATH::JMatrixNS            V;

  private:
    /**
     * Evaluation of fit.
     *
     * \param  __begin        begin of data
     * \param  __end          end   of data
     */
    template<class T>
    inline void evaluate(T __begin, T __end)
    {
      using namespace std;
      using namespace JPP;

      successor = 0.0;

      V.reset();

      for (std::vector<double>::iterator i = Y.begin(); i != Y.end(); ++i) {
        *i = 0.0;
      }

      for (T hit = __begin; hit != __end; ++hit) {

        const JGEOMETRY::JString& string     = geometry[hit->getString()];
        const JMODEL   ::JString& parameters = value.string[hit->getString()];
        const JPosition3D         position   = string.getPosition(parameters, hit->getFloor());

        const double D     = hit->getDistance(position);
        const double Vi    = velocity.getInverseVelocity(D, hit->getZ(), position.getZ());
        const double toa_s = value.emission[hit->getEKey()].t1  +  D * Vi;
          
        const double u     = (toa_s - hit->getValue()) / hit->getSigma();
        const double W     = sqrt(hit->getWeight());

        successor += (W*W) * estimator->getRho(u);
        
        H_t H(1.0, string.getGradient(parameters, hit->getPosition(), hit->getFloor()) * Vi);

        H    *= W * estimator->getPsi(u) / hit->getSigma();

        I_t i;

        i.t1  =  value.getIndex(hit->getEKey(),   &H_t::t1);
        i.tx  =  value.getIndex(hit->getString(), &H_t::tx);
        i.ty  =  value.getIndex(hit->getString(), &H_t::ty);
        i.tx2 =  value.getIndex(hit->getString(), &H_t::tx2);
        i.ty2 =  value.getIndex(hit->getString(), &H_t::ty2);
        i.vs  =  value.getIndex(hit->getString(), &H_t::vs);

        V(i.t1, i.t1) += H.t1 * H.t1;

        Y[i.t1] += W * H.t1;

        if (hit->getFloor() != 0) {

          switch (this->option) {

          case JMODEL::fit_emitters_and_strings_2nd_order_and_stretching_t:
            V(i.t1,  i.vs)  += H.t1  * H.vs;  V(i.tx,  i.vs)  += H.tx  * H.vs;  V(i.ty,  i.vs)  += H.ty  * H.vs;  V(i.tx2, i.vs)  += H.tx2 * H.vs;  V(i.ty2, i.vs)  += H.ty2 * H.vs;

            V(i.vs,  i.t1)  += H.vs  * H.t1;
            V(i.vs,  i.tx)  += H.vs  * H.tx;
            V(i.vs,  i.ty)  += H.vs  * H.ty;
            V(i.vs,  i.tx2) += H.vs  * H.tx2;
            V(i.vs,  i.ty2) += H.vs  * H.ty2;
            
            V(i.vs,  i.vs)  += H.vs  * H.vs;

            Y[i.vs]  += W * H.vs;
            
          case JMODEL::fit_emitters_and_strings_2nd_order_t:
            V(i.t1,  i.tx2) += H.t1  * H.tx2;  V(i.tx,  i.tx2) += H.tx  * H.tx2;  V(i.ty,  i.tx2) += H.ty  * H.tx2;

            V(i.tx2, i.t1)  += H.tx2 * H.t1;
            V(i.tx2, i.tx)  += H.tx2 * H.tx;
            V(i.tx2, i.ty)  += H.tx2 * H.ty;

            V(i.t1,  i.ty2) += H.t1  * H.ty2;  V(i.tx,  i.ty2) += H.tx  * H.ty2;  V(i.ty,  i.ty2) += H.ty  * H.ty2;

            V(i.ty2, i.t1)  += H.ty2 * H.t1;
            V(i.ty2, i.tx)  += H.ty2 * H.tx;
            V(i.ty2, i.ty)  += H.ty2 * H.ty;

            V(i.tx2, i.tx2) += H.tx2 * H.tx2;  V(i.tx2, i.ty2) += H.tx2 * H.ty2;
            V(i.ty2, i.tx2) += H.ty2 * H.tx2;  V(i.ty2, i.ty2) += H.ty2 * H.ty2;

            Y[i.tx2] += W * H.tx2;
            Y[i.ty2] += W * H.ty2;

          case JMODEL::fit_emitters_and_strings_1st_order_t:
            V(i.t1,  i.tx)  += H.t1  * H.tx;   V(i.t1,  i.ty)  += H.t1  * H.ty;
            V(i.tx,  i.t1)  += H.tx  * H.t1;   V(i.ty,  i.t1)  += H.ty  * H.t1;

            V(i.tx,  i.tx)  += H.tx  * H.tx;   V(i.tx,  i.ty)  += H.tx  * H.ty;
            V(i.ty,  i.tx)  += H.ty  * H.tx;   V(i.ty,  i.ty)  += H.ty  * H.ty;

            Y[i.tx]  += W * H.tx;
            Y[i.ty]  += W * H.ty;
            break;
          
          default:
            break;
          }
        }
      }
    }


    std::vector<double>         Y;                             // gradient
    result_type                 successor;
    JModel                      previous;
    std::vector<double>         h;                             // normalisation vector
  };
}

#endif