Jpp  master_rocky-40-g5f0272dcd
the software that should make you happy
JKatoomba_t.hh
Go to the documentation of this file.
1 #ifndef __JACOUSTICS__JKATOOMBA_T__
2 #define __JACOUSTICS__JKATOOMBA_T__
3 
4 #include <vector>
5 #include <memory>
6 #include <mutex>
7 
8 #include "TMatrixTSym.h"
9 #include "TDecompSVD.h"
10 
11 #include "JFit/JEstimator.hh"
12 #include "JFit/JRegressor.hh"
13 #include "JFit/JSimplex.hh"
14 #include "JFit/JMEstimator.hh"
15 #include "JMath/JVectorND.hh"
16 #include "JMath/JMatrixNS.hh"
17 #include "JLang/JType.hh"
18 #include "JLang/JSharedPointer.hh"
19 #include "JLang/JException.hh"
20 #include "JMath/JMatrixNS.hh"
21 #include "Jeep/JMessage.hh"
22 
25 #include "JAcoustics/JEKey.hh"
26 #include "JAcoustics/JGeometry.hh"
27 #include "JAcoustics/JModel.hh"
28 #include "JAcoustics/JHit.hh"
30 
31 
32 /**
33  * \file
34  *
35  * Fit functions of acoustic model.
36  * \author mdejong
37  */
38 namespace JACOUSTICS {}
39 namespace JPP { using namespace JACOUSTICS; }
40 
41 namespace JACOUSTICS {
42 
43  using JLANG::JType;
46  using JMATH::JMath;
47  using JFIT::JEstimator;
49  using JFIT::JSimplex;
50  using JFIT::JMEstimator;
51 
52 
53  /**
54  * Get total weight of data points.
55  *
56  * \param __begin begin of data
57  * \param __end end of data
58  * \return total weight
59  */
60  template<class T>
61  inline double getWeight(T __begin, T __end)
62  {
63  double W = 0;
64 
65  for (T i = __begin; i != __end; ++i) {
66  W += i->getWeight();
67  }
68 
69  return W;
70  }
71 
72 
73  /**
74  * Template definition of fit function of acoustic model.
75  */
76  template<template<class T> class JMinimiser_t = JType>
77  struct JKatoomba;
78 
79 
80  /**
81  * Auxiliary base class for fit function of acoustic model.
82  */
83  template<>
84  struct JKatoomba<JType> {
85 
86  typedef std::shared_ptr<JMEstimator> estimator_type;
87 
88 
89  /**
90  * Auxiliary data structure to sort transmissions.
91  */
92  static const struct compare {
93  /**
94  * Sort transmissions in following order.
95  * - ascending identifier;
96  * - descending quality;
97  * - ascending time-of-arrival;
98  *
99  * \param first first transmission
100  * \param second second transmission
101  * \return true if first transmission before second; else false
102  */
103  inline bool operator()(const JTransmission& first, const JTransmission& second) const
104  {
105  if (first.getID() == second.getID()) {
106 
107  if (first.getQ() == second.getQ())
108  return first.getToA() < second.getToA();
109  else
110  return first.getQ() > second.getQ();
111 
112  } else {
113 
114  return first.getID() < second.getID();
115  }
116  }
117  } compare;
118 
119 
120  /**
121  * Constructor.
122  * The option corresponds to the use of fit parameters in the model of the detector geometry.\n
123  * A negative implies that all strings in the detector use common fit parameters.
124  *
125  * \param geometry detector geometry
126  * \param velocity sound velocity
127  * \param option option
128  */
129  JKatoomba(const JGeometry& geometry,
130  const JSoundVelocity& velocity,
131  const int option) :
132  geometry(geometry),
133  velocity(velocity),
134  option (option)
135  {
137 
138  if (option < 0) {
139  this->option = -option;
140  }
141  };
142 
143 
144  /**
145  * Get estimated time-of-arrival for given hit.
146  *
147  * \param model model
148  * \param hit hit
149  * \return time-of-arrival
150  */
151  double getToA(const JModel& model, const JHit& hit) const
152  {
153  const JGEOMETRY::JString& string = geometry[hit.getString()];
154  const JMODEL ::JString& parameters = model.string[hit.getString()];
155  const JPosition3D position = string.getPosition(parameters, hit.getFloor());
156 
157  const double D = hit.getDistance(position);
158  const double Vi = velocity.getInverseVelocity(D, hit.getZ(), position.getZ());
159 
160  return model.emission[hit.getEKey()].t1 + D * Vi;
161  }
162 
163 
164  /**
165  * Get estimated time-of-emission for given hit.
166  *
167  * \param model model
168  * \param hit hit
169  * \return time-of-arrival
170  */
171  double getToE(const JModel& model, const JHit& hit) const
172  {
173  const JGEOMETRY::JString& string = geometry[hit.getString()];
174  const JMODEL ::JString& parameters = model.string[hit.getString()];
175  const JPosition3D position = string.getPosition(parameters, hit.getFloor());
176 
177  const double D = hit.getDistance(position);
178  const double Vi = velocity.getInverseVelocity(D, hit.getZ(), position.getZ());
179 
180  return hit.getValue() - D * Vi;
181  }
182 
183 
184  const JGeometry& geometry; //!< detector geometry
185  JSoundVelocity velocity; //!< sound velocity
186  int option; //!< fit option
187  estimator_type estimator; //!< M-Estimator function
188 
189  protected:
190  /**
191  * H-equation as per hit.
192  */
193  struct H_t :
194  public JMODEL::JEmission,
195  public JMODEL::JString,
196  public JMath<H_t>
197  {
198  /**
199  * Default constructor.
200  */
201  H_t() :
202  JMODEL::JEmission(),
203  JMODEL::JString ()
204  {}
205 
206 
207  /**
208  * Constructor.
209  *
210  * \param emission emission
211  * \param string string
212  */
213  H_t(const JMODEL::JEmission& emission,
214  const JMODEL::JString& string) :
215  JMODEL::JEmission(emission),
216  JMODEL::JString (string)
217  {}
218 
219 
220  /**
221  * Scale H-equation.
222  *
223  * \param factor multiplication factor
224  * \return this H-equation
225  */
226  H_t& mul(const double factor)
227  {
228  static_cast<JMODEL::JEmission&>(*this).mul(factor);
229  static_cast<JMODEL::JString&> (*this).mul(factor);
230 
231  return *this;
232  }
233  };
234 
235 
236  /**
237  * Indices of H-equation in global model.
238  */
239  struct I_t
240  {
241  /**
242  * Default constructor.
243  */
244  I_t() :
245  t1 (0),
246  tx (0),
247  ty (0),
248  tx2(0),
249  ty2(0),
250  vs (0)
251  {}
252 
253 
254  int t1;
255  int tx;
256  int ty;
257  int tx2;
258  int ty2;
259  int vs;
260  };
261  };
262 
263 
264  /**
265  * Template specialisation of fit function of acoustic model based on JAbstractMinimiser minimiser.\n
266  * This class can be used to evaluate the chi2.
267  */
268  template<>
270  public JAbstractMinimiser<JModel>,
271  public JKatoomba<>
272  {
273  typedef double result_type;
274 
275 
276  /**
277  * Constructor
278  * The option corresponds to the use of fit parameters in the model of the detector geometry.\n
279  * A negative implies that all strings in the detector use common fit parameters.
280  *
281  * \param geometry detector geometry
282  * \param velocity sound velocity
283  * \param option option
284  */
285  JKatoomba(const JGeometry& geometry,
286  const JSoundVelocity& velocity,
287  const int option) :
288  JKatoomba<>(geometry, velocity, option)
289  {};
290 
291 
292  /**
293  * Fit function.\n
294  * This method is used to determine the chi2 of given hit with respect to actual model.
295  *
296  * \param model model
297  * \param hit hit
298  * \return chi2
299  */
300  result_type operator()(const JModel& model, const JHit& hit) const
301  {
302  const double toa_s = this->getToA(model, hit);
303  const double u = (toa_s - hit.getValue()) / hit.getSigma();
304 
305  return estimator->getRho(u) * hit.getWeight();
306  }
307 
308 
309  /**
310  * Fit.
311  *
312  * \param __begin begin of hits
313  * \param __end end of hits
314  * \return chi2
315  */
316  template<class T>
317  result_type operator()(T __begin, T __end)
318  {
319  this->value.setOption(this->option);
320 
321  return JAbstractMinimiser<JModel>::operator()(*this, __begin, __end);
322  }
323 
324 
325  /**
326  * Fit.
327  *
328  * \param model model
329  * \param __begin begin of hits
330  * \param __end end of hits
331  * \return chi2
332  */
333  template<class T>
334  result_type operator()(const JModel& model, T __begin, T __end)
335  {
336  this->value = model;
337 
338  return (*this)(__begin, __end);
339  }
340  };
341 
342 
343  /**
344  * Algorithm for matrix inversion.
345  */
346  enum JMatrix_t {
347  SVD_t = 1, //!< SVD
348  LDU_t = 2 //!< LDU
349  };
350 
351 
352  /**
353  * Auxiliary data structure for compatibility of symmetric matrix.
354  * Inversion is based on SVD algorithm.
355  */
356  struct TMatrixDS :
357  public TMatrixD
358  {
359  static constexpr double TOLERANCE = 1.0e-8; //!< Tolerance for SVD
360 
361  /**
362  * Resize matrix.
363  */
364  void resize(const size_t size)
365  {
366  ResizeTo(size, size);
367  }
368 
369  /**
370  * Set matrix to the null matrix.
371  */
372  void reset()
373  {
374  Zero();
375  }
376 
377  /**
378  * Invert matrix according SVD decomposition.
379  */
380  void invert()
381  {
382 #ifdef THREAD_SAFE
383  using namespace std;
384 
385  unique_lock<mutex> lock(mtx);
386 #endif
387  TDecompSVD svd(*this, TOLERANCE);
388 
389  Bool_t status;
390 
391  static_cast<TMatrixD&>(*this) = svd.Invert(status);
392  }
393 
394  /**
395  * Get solution of equation <tt>A x = b</tt>.
396  *
397  * \param u column vector; b on input, x on output(I/O)
398  */
399  template<class JVectorND_t>
400  void solve(JVectorND_t& u)
401  {
402 #ifdef THREAD_SAFE
403  using namespace std;
404 
405  unique_lock<mutex> lock(mtx);
406 #endif
407  TDecompSVD svd(*this, TOLERANCE);
408 
409  TVectorD b(u.size());
410 
411  for (size_t i = 0; i != u.size(); ++i) {
412  b[i] = u[i];
413  }
414 
415  svd.Solve(b);
416 
417  for (size_t i = 0; i != u.size(); ++i) {
418  u[i] = b[i];
419  }
420  }
421 
422  static std::mutex mtx; //!< TDecompSVD
423  };
424 
425 
426  /**
427  * Using declaration for compatibility of symmetric matrix.
428  * Inversion is based on LDU algorithm.
429  */
430  using JMATH::JMatrixNS;
431 
432 
433  /**
434  * Template specialisation of fit function of acoustic model based on linear approximation.
435  */
436  template<>
438  public JKatoomba<>
439  {
440  /**
441  * Constructor
442  * The option corresponds to the use of fit parameters in the model of the detector geometry.\n
443  * A negative implies that all strings in the detector use common fit parameters.
444  *
445  * \param geometry detector geometry
446  * \param velocity sound velocity
447  * \param option option
448  */
449  JKatoomba(const JGeometry& geometry,
450  const JSoundVelocity& velocity,
451  const int option) :
452  JKatoomba<>(geometry, velocity, option)
453  {};
454 
455 
456  /**
457  * Get start values of string parameters.\n
458  * Note that this method may throw an exception.
459  *
460  * \param __begin begin of hits
461  * \param __end end of hits
462  * \return model
463  */
464  template<class T>
465  JModel& operator()(T __begin, T __end) const
466  {
467  switch (MATRIX_INVERSION) {
468 
469  case SVD_t:
470  return this->evaluate(__begin, __end, svd);
471 
472  case LDU_t:
473  return this->evaluate(__begin, __end, ldu);
474 
475  default:
476  THROW(JValueOutOfRange, "Invalid matrix type " << MATRIX_INVERSION);
477  }
478  }
479 
480 
481  static JMatrix_t MATRIX_INVERSION; // matrix inversion
482 
483  private:
484  mutable JMatrixNS ldu;
485  mutable TMatrixDS svd;
486 
488 
489  /**
490  * Get start values of string parameters.\n
491  * Note that this method may throw an exception.
492  *
493  * \param __begin begin of hits
494  * \param __end end of hits
495  * \param V matrix
496  * \return model
497  */
498  template<class T, class JMatrixNS_t>
499  JModel& evaluate(T __begin, T __end, JMatrixNS_t& V) const
500  {
501  using namespace std;
502  using namespace JPP;
503  using namespace JGEOMETRY;
504 
505  value = JModel(__begin, __end); // set up global model according data
506 
507  if (this->option == JMODEL::FIT_EMITTERS_ONLY_t)
508  value.setOption(JMODEL::FIT_EMITTERS_ONLY_t);
509  else
511 
512  H_t H; // H-equation as per hit
513  I_t i; // indices of H-equation in global model
514 
515  const size_t N = value.getN();
516 
517  V.resize(N);
518  Y.resize(N);
519 
520  V.reset();
521  Y.reset();
522 
523  for (T hit = __begin; hit != __end; ++hit) {
524 
525  const JString& string = geometry[hit->getString()];
526  const JPosition3D position = string.getPosition(hit->getFloor());
527  const double Vi = velocity.getInverseVelocity(hit->getDistance(position), hit->getZ(), position.getZ());
528 
529  const double h1 = string.getHeight(hit->getFloor());
530  const JPosition3D p1 = string.getPosition() - hit->getPosition();
531  const double ds = sqrt(p1.getLengthSquared() + h1*h1 + 2.0*p1.getZ()*h1);
532  const double y = hit->getValue() - Vi*ds;
533  const double W = sqrt(hit->getWeight());
534 
535  H.t1 = W * 1.0;
536  H.tx = W * Vi * p1.getX() * h1 / ds;
537  H.ty = W * Vi * p1.getY() * h1 / ds;
538 
539  i.t1 = value.getIndex(hit->getEKey(), &H_t::t1);
540  i.tx = value.getIndex(hit->getString(), &H_t::tx);
541  i.ty = value.getIndex(hit->getString(), &H_t::ty);
542 
543  V(i.t1, i.t1) += H.t1 * H.t1;
544 
545  Y[i.t1] += W * H.t1 * y;
546 
547  if (hit->getFloor() != 0) {
548 
549  if (this->option != JMODEL::FIT_EMITTERS_ONLY_t) {
550 
551  V(i.t1, i.tx) += H.t1 * H.tx; V(i.t1, i.ty) += H.t1 * H.ty;
552  V(i.tx, i.t1) += H.tx * H.t1; V(i.ty, i.t1) += H.ty * H.t1;
553 
554  V(i.tx, i.tx) += H.tx * H.tx; V(i.tx, i.ty) += H.tx * H.ty;
555  V(i.ty, i.tx) += H.ty * H.tx; V(i.ty, i.ty) += H.ty * H.ty;
556 
557  Y[i.tx] += W * H.tx * y;
558  Y[i.ty] += W * H.ty * y;
559  }
560  }
561  }
562 
563  // solve A x = b
564 
565  V.solve(Y);
566 
567  for (size_t row = 0; row != N; ++row) {
568  value[row] += Y[row];
569  }
570 
571  return value;
572  }
573 
574 
575  mutable JModel value;
576  };
577 
578 
579  /**
580  * Template specialisation of fit function of acoustic model based on JSimplex minimiser.
581  */
582  template<>
583  struct JKatoomba<JSimplex> :
584  public JSimplex<JModel>,
585  public JKatoomba<>
586  {
587  /**
588  * Constructor
589  * The option corresponds to the use of fit parameters in the model of the detector geometry.\n
590  * A negative implies that all strings in the detector use common fit parameters.
591  *
592  * \param geometry detector geometry
593  * \param velocity sound velocity
594  * \param option option
595  */
596  JKatoomba(const JGeometry& geometry,
597  const JSoundVelocity& velocity,
598  const int option) :
599  JKatoomba<>(geometry, velocity, option)
600  {};
601 
602 
603  /**
604  * Fit function.\n
605  * This method is used to determine the chi2 of given hit with respect to actual model.
606  *
607  * \param model model
608  * \param hit hit
609  * \return chi2
610  */
611  double operator()(const JModel& model, const JHit& hit) const
612  {
613  const double toa_s = this->getToA(model, hit);
614  const double u = (toa_s - hit.getValue()) / hit.getSigma();
615 
616  return estimator->getRho(u);
617  }
618 
619 
620  /**
621  * Fit.
622  *
623  * \param __begin begin of hits
624  * \param __end end of hits
625  * \return chi2
626  */
627  template<class T>
628  double operator()(T __begin, T __end)
629  {
630  this->step.clear();
631 
632  for (JModel::string_type::const_iterator i = this->value.string.begin(); i != this->value.string.end(); ++i) {
633 
634  JModel model;
635 
636  switch (this->option) {
637 
639  model.string[i->first] = JMODEL::JString(0.0, 0.0, 0.0, 0.0, 1.0e-6); this->step.push_back(model);
640 
642  model.string[i->first] = JMODEL::JString(0.0, 0.0, 1.0e-6, 0.0, 0.0); this->step.push_back(model);
643  model.string[i->first] = JMODEL::JString(0.0, 0.0, 0.0, 1.0e-6, 0.0); this->step.push_back(model);
644 
646  model.string[i->first] = JMODEL::JString(1.0e-3, 0.0, 0.0, 0.0, 0.0); this->step.push_back(model);
647  model.string[i->first] = JMODEL::JString(0.0, 1.0e-3, 0.0, 0.0, 0.0); this->step.push_back(model);
648  break;
649 
650  default:
651  break;
652  }
653  }
654 
655  for (JModel::emission_type::const_iterator i = this->value.emission.begin(); i != this->value.emission.end(); ++i) {
656 
657  JModel model;
658 
659  model.emission[i->first] = JMODEL::JEmission(5.0e-5); this->step.push_back(model);
660  }
661 
662  return JSimplex<JModel>::operator()(*this, __begin, __end);
663  }
664  };
665 
666 
667  /**
668  * Place holder for custom implementation.
669  */
670  template<class T>
671  class JGandalf {};
672 
673 
674  /**
675  * Template specialisation of fit function of acoustic model based on JGandalf minimiser.
676  */
677  template<>
678  struct JKatoomba<JGandalf> :
679  public JKatoomba<>
680  {
681  typedef double result_type;
682 
683  /**
684  * Type definition internal data structure.
685  */
687 
688 
689  /**
690  * Constructor
691  * The option corresponds to the use of fit parameters in the model of the detector geometry.\n
692  * A negative implies that all strings in the detector use common fit parameters.
693  *
694  * \param geometry detector geometry
695  * \param velocity sound velocity
696  * \param option option
697  */
698  JKatoomba(const JGeometry& geometry,
699  const JSoundVelocity& velocity,
700  const int option) :
701  JKatoomba<>(geometry, velocity, option)
702  {};
703 
704 
705  /**
706  * Fit.
707  *
708  * \param __begin begin of hits
709  * \param __end end of hits
710  * \return chi2 and gradient
711  */
712  template<class T>
713  result_type operator()(T __begin, T __end)
714  {
715  using namespace std;
716  using namespace JPP;
717 
718  value.setOption(this->option);
719 
720  const int N = value.getN();
721 
722  V.resize(N);
723  Y.resize(N);
724  h.resize(N);
725 
726  data_type data;
727 
728  for (T hit = __begin; hit != __end; ++hit) {
729  data[hit->getLocation()][hit->getEmitter()].push_back(*hit);
730  }
731 
732  lambda = LAMBDA_MIN;
733 
734  result_type precessor = numeric_limits<double>::max();
735 
736  for (numberOfIterations = 0; numberOfIterations != MAXIMUM_ITERATIONS; ++numberOfIterations) {
737 
738  DEBUG("step: " << numberOfIterations << endl);
739 
740  evaluate(data);
741 
742  DEBUG("lambda: " << FIXED(12,5) << lambda << endl);
743  DEBUG("chi2: " << FIXED(12,5) << successor << endl);
744 
745  if (successor < precessor) {
746 
747  if (numberOfIterations != 0) {
748 
749  if (fabs(precessor - successor) < EPSILON*fabs(precessor)) {
750  return successor;
751  }
752 
753  if (lambda > LAMBDA_MIN) {
754  lambda /= LAMBDA_DOWN;
755  }
756  }
757 
758  precessor = successor;
759  previous = value;
760 
761  } else {
762 
763  value = previous;
764  lambda *= LAMBDA_UP;
765 
766  if (lambda > LAMBDA_MAX) {
767  return precessor; // no improvement found
768  }
769 
770  evaluate(data);
771  }
772 
773  // force definite positiveness
774 
775  for (int i = 0; i != N; ++i) {
776 
777  if (V(i,i) < PIVOT) {
778  V(i,i) = PIVOT;
779  }
780 
781  h[i] = 1.0 / sqrt(V(i,i));
782  }
783 
784  // normalisation
785 
786  for (int i = 0; i != N; ++i) {
787  for (int j = 0; j != i; ++j) {
788  V(j,i) *= h[i] * h[j];
789  V(i,j) = V(j,i);
790  }
791  }
792 
793  for (int i = 0; i != N; ++i) {
794  V(i,i) = 1.0 + lambda;
795  }
796 
797 
798  // solve A x = b
799 
800  for (int col = 0; col != N; ++col) {
801  Y[col] *= h[col];
802  }
803 
804  try {
805  V.solve(Y);
806  }
807  catch (const exception& error) {
808 
809  ERROR("JGandalf: " << error.what() << endl << V << endl);
810 
811  break;
812  }
813 
814  // update value
815 
816  for (int row = 0; row != N; ++row) {
817 
818  DEBUG("u[" << noshowpos << setw(3) << row << "] = " << showpos << FIXED(20,5) << value[row]);
819 
820  value[row] -= h[row] * Y[row];
821 
822  DEBUG(" -> " << FIXED(20,10) << value[row] << noshowpos << endl);
823  }
824  }
825 
826  return precessor;
827  }
828 
829  static int debug; //!< debug level
830  static int MAXIMUM_ITERATIONS; //!< maximal number of iterations
831  static double EPSILON; //!< maximal distance to minimum
832  static double LAMBDA_MIN; //!< minimal value control parameter
833  static double LAMBDA_MAX; //!< maximal value control parameter
834  static double LAMBDA_UP; //!< multiplication factor control parameter
835  static double LAMBDA_DOWN; //!< multiplication factor control parameter
836  static double PIVOT; //!< minimal value diagonal element of matrix
837 
838  double lambda;
842 
843  private:
844  /**
845  * Evaluation of fit.
846  *
847  * \param data data
848  */
849  inline void evaluate(const data_type& data)
850  {
851  using namespace std;
852  using namespace JPP;
853 
854  successor = 0.0;
855 
856  V.reset();
857  Y.reset();
858 
859  for (data_type::const_iterator p = data.begin(); p != data.end(); ++p) {
860 
861  const JGEOMETRY::JString& string = geometry [p->first.getString()];
862  const JMODEL ::JString& parameters = value.string[p->first.getString()];
863  const JPosition3D position = string.getPosition(parameters, p->first.getFloor());
864 
865  for (data_type::mapped_type::const_iterator emitter = p->second.begin(); emitter != p->second.end(); ++emitter) {
866 
867  const double D = emitter->first.getDistance(position);
868  const double Vi = velocity.getInverseVelocity(D, emitter->first.getZ(), position.getZ());
869 
870  const H_t H0(1.0, string.getGradient(parameters, emitter->first.getPosition(), p->first.getFloor()) * Vi);
871 
872  for (data_type::mapped_type::mapped_type::const_iterator hit = emitter->second.begin(); hit != emitter->second.end(); ++hit) {
873 
874  const double toa_s = value.emission[hit->getEKey()].t1 + D * Vi;
875 
876  const double u = (toa_s - hit->getValue()) / hit->getSigma();
877  const double W = sqrt(hit->getWeight());
878 
879  successor += (W*W) * estimator->getRho(u);
880 
881  const H_t H = H0 * (W * estimator->getPsi(u) / hit->getSigma());
882 
883  I_t i;
884 
885  i.t1 = value.getIndex(hit->getEKey(), &H_t::t1);
886  i.tx = value.getIndex(hit->getString(), &H_t::tx);
887  i.ty = value.getIndex(hit->getString(), &H_t::ty);
888  i.tx2 = value.getIndex(hit->getString(), &H_t::tx2);
889  i.ty2 = value.getIndex(hit->getString(), &H_t::ty2);
890  i.vs = value.getIndex(hit->getString(), &H_t::vs);
891 
892  V(i.t1, i.t1) += H.t1 * H.t1;
893 
894  Y[i.t1] += W * H.t1;
895 
896  if (hit->getFloor() != 0) {
897 
898  switch (this->option) {
899 
901  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;
902 
903  V(i.vs, i.t1) = V(i.t1, i.vs);
904  V(i.vs, i.tx) = V(i.tx, i.vs);
905  V(i.vs, i.ty) = V(i.ty, i.vs);
906  V(i.vs, i.tx2) = V(i.tx2, i.vs);
907  V(i.vs, i.ty2) = V(i.ty2, i.vs);
908 
909  V(i.vs, i.vs) += H.vs * H.vs;
910 
911  Y[i.vs] += W * H.vs;
912 
914  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;
915 
916  V(i.tx2, i.t1) = V(i.t1, i.tx2);
917  V(i.tx2, i.tx) = V(i.tx, i.tx2);
918  V(i.tx2, i.ty) = V(i.ty, i.tx2);
919 
920  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;
921 
922  V(i.ty2, i.t1) = V(i.t1, i.ty2);
923  V(i.ty2, i.tx) = V(i.tx, i.ty2);
924  V(i.ty2, i.ty) = V(i.ty, i.ty2);
925 
926  V(i.tx2, i.tx2) += H.tx2 * H.tx2; V(i.tx2, i.ty2) += H.tx2 * H.ty2;
927  V(i.ty2, i.tx2) = V(i.tx2, i.ty2); V(i.ty2, i.ty2) += H.ty2 * H.ty2;
928 
929  Y[i.tx2] += W * H.tx2;
930  Y[i.ty2] += W * H.ty2;
931 
933  V(i.t1, i.tx) += H.t1 * H.tx; V(i.t1, i.ty) += H.t1 * H.ty;
934  V(i.tx, i.t1) = V(i.t1, i.tx); V(i.ty, i.t1) = V(i.t1, i.ty);
935 
936  V(i.tx, i.tx) += H.tx * H.tx; V(i.tx, i.ty) += H.tx * H.ty;
937  V(i.ty, i.tx) = V(i.tx, i.ty); V(i.ty, i.ty) += H.ty * H.ty;
938 
939  Y[i.tx] += W * H.tx;
940  Y[i.ty] += W * H.ty;
941  break;
942 
943  default:
944  break;
945  }
946  }
947  }
948  }
949  }
950  }
951 
952 
953  JMATH::JVectorND Y; // gradient
956  std::vector<double> h; // normalisation vector
957  };
958 }
959 
960 #endif
Acoustic hit.
Model for fit to acoutsics data.
TPaveText * p1
Emitter hash key.
Linear fit methods.
Exceptions.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:712
Acoustic fit parameters.
Acoustic geometries.
Maximum likelihood estimator (M-estimators).
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
#define ERROR(A)
Definition: JMessage.hh:66
General purpose data regression method.
Sound velocity.
Acoustic transmission.
Place holder for custom implementation.
Definition: JKatoomba_t.hh:671
int getFloor() const
Get floor number.
Definition: JLocation.hh:146
int getString() const
Get string number.
Definition: JLocation.hh:135
Abstract minimiser.
Definition: JRegressor.hh:27
result_type operator()(const JFunction_t &fit, T __begin, T __end)
Get chi2.
Definition: JRegressor.hh:48
Template definition of linear fit.
Definition: JEstimator.hh:25
Simple fit method based on Powell's algorithm, see reference: Numerical Recipes in C++,...
Definition: JSimplex.hh:44
double operator()(const JFunction_t &fit, T __begin, T __end)
Multi-dimensional fit.
Definition: JSimplex.hh:71
Data structure for position in three dimensions.
Definition: JPosition3D.hh:38
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition: JVector3D.hh:270
double getZ() const
Get z position.
Definition: JVector3D.hh:115
Exception for division by zero.
Definition: JException.hh:288
Wrapper class around STL string class.
Definition: JString.hh:29
Exception for accessing a value in a collection that is outside of its range.
Definition: JException.hh:180
@ FIT_EMITTERS_AND_STRINGS_1st_ORDER_t
fit times of emission of emitters and tilt angles of strings
@ FIT_EMITTERS_AND_STRINGS_2nd_ORDER_t
fit times of emission of emitters and tilt angles of strings with second order correction
@ FIT_EMITTERS_ONLY_t
fit only times of emission of emitters
@ FIT_EMITTERS_AND_STRINGS_2nd_ORDER_AND_STRETCHING_t
fit times of emission of emitters and tilt angles of strings with second order correction and stretch...
Auxiliary classes and methods for acoustic position calibration.
double getWeight(T __begin, T __end)
Get total weight of data points.
Definition: JKatoomba_t.hh:61
JMatrix_t
Algorithm for matrix inversion.
Definition: JKatoomba_t.hh:346
static const double H
Planck constant [eV s].
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
std::vector< double > vs
double u[N+1]
Definition: JPolint.hh:865
int j
Definition: JPolint.hh:792
Definition: JSTDTypes.hh:14
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:448
Acoustics hit.
double getSigma() const
Get resolution of time-of-arrival.
double getValue() const
Get expectation value of time-of-arrival.
double getWeight() const
Get weight.
JEKey getEKey() const
Get emitter hash key of this hit.
result_type operator()(T __begin, T __end)
Fit.
Definition: JKatoomba_t.hh:317
JKatoomba(const JGeometry &geometry, const JSoundVelocity &velocity, const int option)
Constructor The option corresponds to the use of fit parameters in the model of the detector geometry...
Definition: JKatoomba_t.hh:285
result_type operator()(const JModel &model, T __begin, T __end)
Fit.
Definition: JKatoomba_t.hh:334
result_type operator()(const JModel &model, const JHit &hit) const
Fit function.
Definition: JKatoomba_t.hh:300
JModel & operator()(T __begin, T __end) const
Get start values of string parameters.
Definition: JKatoomba_t.hh:465
JKatoomba(const JGeometry &geometry, const JSoundVelocity &velocity, const int option)
Constructor The option corresponds to the use of fit parameters in the model of the detector geometry...
Definition: JKatoomba_t.hh:449
static JMatrix_t MATRIX_INVERSION
Matrix inversion.
Definition: JKatoomba_t.hh:481
JModel & evaluate(T __begin, T __end, JMatrixNS_t &V) const
Get start values of string parameters.
Definition: JKatoomba_t.hh:499
static double LAMBDA_DOWN
multiplication factor control parameter
Definition: JKatoomba_t.hh:835
static int debug
debug level
Definition: JKatoomba_t.hh:829
static double PIVOT
minimal value diagonal element of matrix
Definition: JKatoomba_t.hh:836
std::map< JLocation, std::map< JEmitter, std::vector< JHit > > > data_type
Type definition internal data structure.
Definition: JKatoomba_t.hh:686
JKatoomba(const JGeometry &geometry, const JSoundVelocity &velocity, const int option)
Constructor The option corresponds to the use of fit parameters in the model of the detector geometry...
Definition: JKatoomba_t.hh:698
result_type operator()(T __begin, T __end)
Fit.
Definition: JKatoomba_t.hh:713
static double LAMBDA_UP
multiplication factor control parameter
Definition: JKatoomba_t.hh:834
static int MAXIMUM_ITERATIONS
maximal number of iterations
Definition: JKatoomba_t.hh:830
static double LAMBDA_MAX
maximal value control parameter
Definition: JKatoomba_t.hh:833
static double EPSILON
maximal distance to minimum
Definition: JKatoomba_t.hh:831
void evaluate(const data_type &data)
Evaluation of fit.
Definition: JKatoomba_t.hh:849
static double LAMBDA_MIN
minimal value control parameter
Definition: JKatoomba_t.hh:832
JKatoomba(const JGeometry &geometry, const JSoundVelocity &velocity, const int option)
Constructor The option corresponds to the use of fit parameters in the model of the detector geometry...
Definition: JKatoomba_t.hh:596
double operator()(const JModel &model, const JHit &hit) const
Fit function.
Definition: JKatoomba_t.hh:611
double operator()(T __begin, T __end)
Fit.
Definition: JKatoomba_t.hh:628
H_t & mul(const double factor)
Scale H-equation.
Definition: JKatoomba_t.hh:226
H_t(const JMODEL::JEmission &emission, const JMODEL::JString &string)
Constructor.
Definition: JKatoomba_t.hh:213
bool operator()(const JTransmission &first, const JTransmission &second) const
Sort transmissions in following order.
Definition: JKatoomba_t.hh:103
const JGeometry & geometry
detector geometry
Definition: JKatoomba_t.hh:184
double getToA(const JModel &model, const JHit &hit) const
Get estimated time-of-arrival for given hit.
Definition: JKatoomba_t.hh:151
JKatoomba(const JGeometry &geometry, const JSoundVelocity &velocity, const int option)
Constructor.
Definition: JKatoomba_t.hh:129
std::shared_ptr< JMEstimator > estimator_type
Definition: JKatoomba_t.hh:86
estimator_type estimator
M-Estimator function.
Definition: JKatoomba_t.hh:187
double getToE(const JModel &model, const JHit &hit) const
Get estimated time-of-emission for given hit.
Definition: JKatoomba_t.hh:171
JSoundVelocity velocity
sound velocity
Definition: JKatoomba_t.hh:185
Template definition of fit function of acoustic model.
Definition: JKatoomba_t.hh:77
JEmission & mul(const double factor)
Scale emission.
JString & mul(const double factor)
Scale string.
static bool singularity
Option for common fit parameters.
Model for fit to acoustics data.
JACOUSTICS::JModel::emission_type emission
size_t getIndex(int id, double JString::*p) const
Get index of fit parameter for given string.
JACOUSTICS::JModel::string_type string
Implementation for depth dependend velocity of sound.
Acoustic transmission.
double getToA() const
Get calibrated time of arrival.
double getQ() const
Get quality.
int getID() const
Get identifier.
Auxiliary data structure for compatibility of symmetric matrix.
Definition: JKatoomba_t.hh:358
void resize(const size_t size)
Resize matrix.
Definition: JKatoomba_t.hh:364
static std::mutex mtx
TDecompSVD.
Definition: JKatoomba_t.hh:422
void reset()
Set matrix to the null matrix.
Definition: JKatoomba_t.hh:372
static constexpr double TOLERANCE
Tolerance for SVD.
Definition: JKatoomba_t.hh:359
void invert()
Invert matrix according SVD decomposition.
Definition: JKatoomba_t.hh:380
void solve(JVectorND_t &u)
Get solution of equation A x = b.
Definition: JKatoomba_t.hh:400
Interface for maximum likelihood estimator (M-estimator).
Definition: JMEstimator.hh:20
Auxiliary class for a type holder.
Definition: JType.hh:19
Auxiliary base class for aritmetic operations of derived class types.
Definition: JMath.hh:347
JMatrixND & reset()
Set matrix to the null matrix.
Definition: JMatrixND.hh:459
N x N symmetric matrix.
Definition: JMatrixNS.hh:30
Nx1 matrix.
Definition: JVectorND.hh:23
void reset()
Reset.
Definition: JVectorND.hh:45
container_type::const_iterator const_iterator
Definition: JHashMap.hh:86