JKatoomba.hh

Go to the documentation of this file.

#ifndef __JACOUSTICS__JKATOOMBA__
#define __JACOUSTICS__JKATOOMBA__

#include <vector>

#include "JFit/JEstimator.hh"
#include "JFit/JRegressor.hh"
#include "JFit/JMEstimator.hh"
#include "JMath/JMatrixNS.hh"
#include "JLang/JSharedPointer.hh"

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


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

namespace JACOUSTICS {

  using JFIT::JEstimator;
  using JFIT::JAbstractMinimiser;
  using JFIT::JMEstimator;


  /**
   * Auxiliary base class for fit function of acoustic model.
   */
  struct JKatoomba_t : 
    public JGeometry
  {
    /**
     * Constructor
     *
     * \param  detector           detector
     * \param  velocity           sound velocity
     */
    JKatoomba_t(const JDetector&      detector,
                const JSoundVelocity& velocity) :
      detector(detector),
      velocity(velocity)
    {};


    /**
     * Get estimated time-of-arrival for given hit.
     *
     * \param  model              model
     * \param  hit                hit
     * \return                    time-of-arrival
     */
    template<class JPDF_t>
    double getToA(const JModel& model, const JHit<JPDF_t>& hit) const
    {
      const JString&         string     = detector[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.emitter[hit.getEKey()].t1  +  D * Vi;
    }


    const JDetector&       detector;
    const JSoundVelocity&  velocity;

    JLANG::JSharedPointer<JMEstimator> estimator;  //!< M-Estimator function
  };


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


  /**
   * 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_t
  {
    typedef double                                        result_type;


    /**
     * Constructor
     *
     * \param  detector           detector
     * \param  velocity           sound velocity
     */
    JKatoomba(const JDetector&      detector,
              const JSoundVelocity& velocity) :
      JKatoomba_t(detector, velocity)
    {};


    /**
     * 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 and gradient
     */
    template<class JPDF_t>
    result_type operator()(const JModel& model, const JHit<JPDF_t>& hit) const
    {
      const double toa_s = this->getToA(model, hit);
      const double u     = (toa_s - hit.getValue()) / hit.sigma;

      return estimator->getRho(u);
    }


    /**
     * 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 JAbstractMinimiser<JModel>::operator()(*this, __begin, __end);
    }
  };


  /**
   * Template specialisation of fit function of acoustic model based on linear approximation.
   */
  template<>
  struct JKatoomba<JEstimator> : 
    public JKatoomba_t
  {
    /**
     * Constructor
     *
     * \param  detector           detector
     * \param  velocity           sound velocity
     */
    JKatoomba(const JDetector&      detector,
              const JSoundVelocity& velocity) :
      JKatoomba_t(detector, velocity)
    {};


    /**
     * 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>
    const JModel& operator()(T __begin, T __end) const
    {
      using namespace std;
      using namespace JPP;

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

      struct H_t :                                    // H-equation as per hit
            public JModel::JString,
            public JModel::JEmitter
      {} H;

      struct { int t1; int tx; int ty; } i;           // indices of H-equation in global model 

      JMatrixNS      V(value.getN());                 // V = H^T x H
      vector<double> Y(value.getN(), 0.0);            // Y = H^T x y
      

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

        const JString&    string   = detector[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[hit->getFloor()].getHeight();
        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;

        H.t1  =  1.0;
        H.tx  =  Vi * p1.getX() * h1 / ds;
        H.ty  =  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;
        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.t1]  +=  H.t1 * y;
        Y[i.tx]  +=  H.tx * y;
        Y[i.ty]  +=  H.ty * y;
      }

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

      V.invert();

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

  private:
    mutable JModel value;
  };
}

#endif