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/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;
  }
 
 
  /**
   * 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.
     * The option corresponds to the use of fit parameters in the model of the detector geometry.\n
     * A negative implies that all strings in the detector use common fit parameters.  
     *
     * \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)
    {
      JModel::hash_evaluator::singularity = (option < 0);
 
      if (option < 0) {
        this->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
     * The option corresponds to the use of fit parameters in the model of the detector geometry.\n
     * A negative implies that all strings in the detector use common fit parameters.  
     *
     * \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
     * The option corresponds to the use of fit parameters in the model of the detector geometry.\n
     * A negative implies that all strings in the detector use common fit parameters.  
     *
     * \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
     * The option corresponds to the use of fit parameters in the model of the detector geometry.\n
     * A negative implies that all strings in the detector use common fit parameters.  
     *
     * \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
     * The option corresponds to the use of fit parameters in the model of the detector geometry.\n
     * A negative implies that all strings in the detector use common fit parameters.  
     *
     * \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