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/JVectorND.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 JLANG::JDivisionByZero;
  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;
  }


  /**
   * Auxiliary data structure for average.
   */
  struct JQuantile {
    /**
     * Default constructor.
     */
    JQuantile()
    {
      reset();
    }


    /**
     * Reset.
     */
    void reset()
    {
      total = 0.0;
      count = 0;
    }


    /**
     * Put value.
     *
     * \param  x           value
     */
    void put(const double x)
    {
      total += x;
      count += 1;
    }


    /**
     * Get mean value.
     *
     * \return             mean value
     */
    long double getMean() const
    {
      if (count != 0.0)
        return total / count;
      else
        THROW(JDivisionByZero, "JQuantile::getMean()");
    }


    /**
     * Get mean value.
     *
     * \param  value       default value
     * \return             mean value
     */
    long double getMean(const double value) const
    {
      if (count != 0.0)
        return getMean();
      else
        return value;
    }

  private:
    long double   total;
    long long int count;
  };


  /**
   * 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
  {
    static constexpr double TOLERANCE = 1.0e-8;         //!< Tolerance for SVD

    /**
     * Resize matrix.
     */
    void resize(const size_t size)
    {
      ResizeTo(size, size);
    }

    /**
     * Set matrix to the null matrix.
     */
    void reset()
    {
      Zero();
    }

    /**
     * Invert matrix according SVD decomposition.
     */
    void invert()
    {
#ifdef THREAD_SAFE
      using namespace std;

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

      Bool_t status;

      static_cast<TMatrixD&>(*this) = svd.Invert(status);
    }

    /**
     * Get solution of equation <tt>A x = b</tt>.
     *
     * \param  u        column vector; b on input, x on output(I/O)
     */
    template<class JVectorND_t>
    void solve(JVectorND_t& u)
    {
#ifdef THREAD_SAFE
      using namespace std;

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

      TVectorD b(u.size());

      for (size_t i = 0; i != u.size(); ++i) {
        b[i] = u[i];
      }

      svd.Solve(b);

      for (size_t i = 0; i != u.size(); ++i) {
        u[i] = b[i];
      }
    }

    static std::mutex mtx;                              //!< TDecompSVD
  };


  /**
   * 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<>
  {
    /**
     * 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, svd);

      case LDU_t:
        return this->evaluate(__begin, __end, ldu);

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


    static JMatrix_t  MATRIX_INVERSION;                 // matrix inversion

  private:
    mutable JMatrixNS ldu;
    mutable TMatrixDS svd;

    mutable JMATH::JVectorND Y;

    /**
     * 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  V                  matrix
     * \return                    model
     */
    template<class T, class JMatrixNS_t>
    JModel& evaluate(T __begin, T __end, JMatrixNS_t& V) const
    {
      using namespace std;
      using namespace JPP;
      using namespace JGEOMETRY;

      value = 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();

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

      V.reset();
      Y.reset();

      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;
          }
        }
      }

      // solve A x = b

      V.solve(Y);
      
      for (size_t row = 0; row != N; ++row) {
        value[row] += Y[row];
      }

      return value;
    }


    mutable JModel value;
  };


  /**
   * 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;

    /**
     * Type definition internal data structure.
     */
    typedef std::map<JLocation, std::map<JEmitter, std::vector<JHit> > >  data_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);

      data_type data;

      for (T hit = __begin; hit != __end; ++hit) {
        data[hit->getLocation()][hit->getEmitter()].push_back(*hit);
      }

      lambda    = LAMBDA_MIN;

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

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

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

        evaluate(data);

        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(data);
        }

        // 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;
        }


        // solve A x = b

        for (int col = 0; col != N; ++col) {
          Y[col] *= h[col];
        }

        try {
          V.solve(Y);
        }
        catch (const exception& error) {

          ERROR("JGandalf: " << error.what() << endl << V << endl);

          break;
        }

        // update value

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

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

          value[row] -= h[row] * Y[row];

          DEBUG(" -> " << FIXED(20,10) << value[row] << 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  data           data
     */
    inline void evaluate(const data_type& data)
    {
      using namespace std;
      using namespace JPP;

      successor = 0.0;

      V.reset();
      Y.reset();

      for (data_type::const_iterator p = data.begin(); p != data.end(); ++p) {

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

        for (data_type::mapped_type::const_iterator emitter = p->second.begin(); emitter != p->second.end(); ++emitter) {

          const double D  = emitter->first.getDistance(position);
          const double Vi = velocity.getInverseVelocity(D, emitter->first.getZ(), position.getZ());

          const H_t    H0(1.0, string.getGradient(parameters, emitter->first.getPosition(), p->first.getFloor()) * Vi);

          for (data_type::mapped_type::mapped_type::const_iterator hit = emitter->second.begin(); hit != emitter->second.end(); ++hit) {
        
            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);
        
            const H_t    H  = H0 * (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)   = V(i.t1,  i.vs);
                V(i.vs,  i.tx)   = V(i.tx,  i.vs);
                V(i.vs,  i.ty)   = V(i.ty,  i.vs);
                V(i.vs,  i.tx2)  = V(i.tx2, i.vs);
                V(i.vs,  i.ty2)  = V(i.ty2, i.vs);

                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)   = V(i.t1,  i.tx2);
                V(i.tx2, i.tx)   = V(i.tx,  i.tx2);
                V(i.tx2, i.ty)   = V(i.ty,  i.tx2);
            
                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)   = V(i.t1,  i.ty2);
                V(i.ty2, i.tx)   = V(i.tx,  i.ty2);
                V(i.ty2, i.ty)   = V(i.ty,  i.ty2);

                V(i.tx2, i.tx2) += H.tx2 * H.tx2;    V(i.tx2, i.ty2) += H.tx2 * H.ty2;
                V(i.ty2, i.tx2)  = V(i.tx2, i.ty2);  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)   = V(i.t1, i.tx);  V(i.ty,  i.t1)   = V(i.t1, i.ty);

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

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


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

#endif