#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