Jpp  17.3.0-rc.2
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JKatoomba.hh
Go to the documentation of this file.
1 #ifndef __JACOUSTICS__JKATOOMBA__
2 #define __JACOUSTICS__JKATOOMBA__
3 
4 #include <vector>
5 #include <memory>
6 #ifdef THREAD_SAFE
7 #include <mutex>
8 #endif
9 #include "TMatrixTSym.h"
10 #include "TDecompSVD.h"
11 
12 #include "JFit/JEstimator.hh"
13 #include "JFit/JRegressor.hh"
14 #include "JFit/JSimplex.hh"
15 #include "JFit/JMEstimator.hh"
16 #include "JMath/JMatrixNS.hh"
17 #include "JLang/JType.hh"
18 #include "JLang/JSharedPointer.hh"
19 #include "JMath/JMatrixNS.hh"
20 #include "Jeep/JMessage.hh"
21 
24 #include "JAcoustics/JEKey.hh"
25 #include "JAcoustics/JGeometry.hh"
26 #include "JAcoustics/JModel.hh"
27 #include "JAcoustics/JHit.hh"
29 #include "JAcoustics/JTimeRange.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;
44  using JMATH::JMath;
45  using JFIT::JEstimator;
47  using JFIT::JSimplex;
48  using JFIT::JMEstimator;
49 
50 
51  /**
52  * Get total weight of data points.
53  *
54  * \param __begin begin of data
55  * \param __end end of data
56  * \return total weight
57  */
58  template<class T>
59  inline double getWeight(T __begin, T __end)
60  {
61  double W = 0;
62 
63  for (T i = __begin; i != __end; ++i) {
64  W += i->getWeight();
65  }
66 
67  return W;
68  }
69 
70 
71  /**
72  * Auxiliary data structure to sort transmissions.
73  */
74  static const struct {
75  /**
76  * Sort transmissions in following order.
77  * - ascending identifier;
78  * - descending quality;
79  * - ascending time-of-arrival;
80  *
81  * \param first first transmission
82  * \param second second transmission
83  * \return true if first transmission before second; else false
84  */
85  inline bool operator()(const JTransmission& first, const JTransmission& second) const
86  {
87  if (first.getID() == second.getID()) {
88 
89  if (first.getQ() == second.getQ())
90  return first.getToA() < second.getToA();
91  else
92  return first.getQ() > second.getQ();
93 
94  } else {
95 
96  return first.getID() < second.getID();
97  }
98  }
99  } compare;
100 
101 
102  /**
103  * Enumeration for fit algorithms.
104  */
105  enum JFit_t {
106  linear_t = 0, //!< Linear fit
107  simplex_t, //!< Simplex fit
108  gandalf_t //!< Gandalf fit
109  };
110 
111 
112  /**
113  * Template definition of fit function of acoustic model.
114  */
115  template<template<class T> class JMinimiser_t = JType>
116  struct JKatoomba;
117 
118 
119  /**
120  * Auxiliary base class for fit function of acoustic model.
121  */
122  template<>
123  struct JKatoomba<JType> :
124  public JGeometry
125  {
126  typedef std::shared_ptr<JMEstimator> estimator_type;
127 
128 
129  /**
130  * Constructor
131  *
132  * \param detector detector
133  * \param velocity sound velocity
134  * \param option option
135  */
137  const JSoundVelocity& velocity,
138  const int option) :
139  detector(detector),
140  velocity(velocity),
141  option (option)
142  {};
143 
144 
145  /**
146  * Get estimated time-of-arrival for given hit.
147  *
148  * \param model model
149  * \param hit hit
150  * \return time-of-arrival
151  */
152  template<class JPDF_t>
153  double getToA(const JModel& model, const JHit<JPDF_t>& hit) const
154  {
155  const JGEOMETRY::JString& string = detector[hit.getString()];
156  const JMODEL ::JString& parameters = model.string[hit.getString()];
157  const JPosition3D position = string.getPosition(parameters, hit.getFloor());
158 
159  const double D = hit.getDistance(position);
160  const double Vi = velocity.getInverseVelocity(D, hit.getZ(), position.getZ());
161 
162  return model.emitter[hit.getEKey()].t1 + D * Vi;
163  }
164 
165 
166  /**
167  * Get estimated time-of-emission for given hit.
168  *
169  * \param model model
170  * \param hit hit
171  * \return time-of-arrival
172  */
173  template<class JPDF_t>
174  double getToE(const JModel& model, const JHit<JPDF_t>& hit) const
175  {
176  const JGEOMETRY::JString& string = detector[hit.getString()];
177  const JMODEL ::JString& parameters = model.string[hit.getString()];
178  const JPosition3D position = string.getPosition(parameters, hit.getFloor());
179 
180  const double D = hit.getDistance(position);
181  const double Vi = velocity.getInverseVelocity(D, hit.getZ(), position.getZ());
182 
183  return hit.getValue() - D * Vi;
184  }
185 
186 
187  const JDetector& detector; //!< detector
188  JSoundVelocity velocity; //!< sound velocity
189  int option; //!< fit option
190  estimator_type estimator; //!< M-Estimator function
191 
192  protected:
193  /**
194  * H-equation as per hit.
195  */
196  struct H_t :
197  public JMODEL::JEmitter,
198  public JMODEL::JString,
199  public JMath<H_t>
200  {
201  /**
202  * Default constructor.
203  */
204  H_t() :
205  JMODEL::JEmitter(),
206  JMODEL::JString ()
207  {}
208 
209 
210  /**
211  * Constructor.
212  *
213  * \param emitter emitter
214  * \param string string
215  */
216  H_t(const JMODEL::JEmitter& emitter,
217  const JMODEL::JString& string) :
218  JMODEL::JEmitter(emitter),
219  JMODEL::JString (string)
220  {}
221 
222 
223  /**
224  * Scale H-equation.
225  *
226  * \param factor multiplication factor
227  * \return this H-equation
228  */
229  H_t& mul(const double factor)
230  {
231  static_cast<JMODEL::JEmitter&>(*this).mul(factor);
232  static_cast<JMODEL::JString&> (*this).mul(factor);
233 
234  return *this;
235  }
236  };
237 
238 
239  /**
240  * Indices of H-equation in global model.
241  */
242  struct I_t
243  {
244  /**
245  * Default constructor.
246  */
247  I_t() :
248  t1 (0),
249  tx (0),
250  ty (0),
251  tx2(0),
252  ty2(0),
253  vs (0)
254  {}
255 
256 
257  int t1;
258  int tx;
259  int ty;
260  int tx2;
261  int ty2;
262  int vs;
263  };
264  };
265 
266 
267  /**
268  * Template specialisation of fit function of acoustic model based on JAbstractMinimiser minimiser.\n
269  * This class can be used to evaluate the chi2.
270  */
271  template<>
273  public JAbstractMinimiser<JModel>,
274  public JKatoomba<>
275  {
276  typedef double result_type;
277 
278 
279  /**
280  * Constructor
281  *
282  * \param detector detector
283  * \param velocity sound velocity
284  * \param option option
285  */
287  const JSoundVelocity& velocity,
288  const int option) :
289  JKatoomba<>(detector, velocity, option)
290  {};
291 
292 
293  /**
294  * Fit function.\n
295  * This method is used to determine the chi2 of given hit with respect to actual model.
296  *
297  * \param model model
298  * \param hit hit
299  * \return chi2
300  */
301  template<class JPDF_t>
302  result_type operator()(const JModel& model, const JHit<JPDF_t>& hit) const
303  {
304  const double toa_s = this->getToA(model, hit);
305  const double u = (toa_s - hit.getValue()) / hit.sigma;
306 
307  return estimator->getRho(u);
308  }
309 
310 
311  /**
312  * Fit.
313  *
314  * \param model model
315  * \param __begin begin of hits
316  * \param __end end of hits
317  * \return chi2
318  */
319  template<class T>
320  result_type operator()(const JModel& model, T __begin, T __end)
321  {
322  this->value = model;
323 
324  return JAbstractMinimiser<JModel>::operator()(*this, __begin, __end);
325  }
326  };
327 
328 
329  /**
330  * Template specialisation of fit function of acoustic model based on linear approximation.
331  */
332  template<>
334  public JKatoomba<>
335  {
336  static constexpr double TOLERANCE = 1.0e-8; //!< Tolerance for SVD
337 
338  /**
339  * Constructor
340  *
341  * \param detector detector
342  * \param velocity sound velocity
343  * \param option option
344  */
346  const JSoundVelocity& velocity,
347  const int option) :
348  JKatoomba<>(detector, velocity, option)
349  {};
350 
351 
352  /**
353  * Get start values of string parameters.\n
354  * Note that this method may throw an exception.
355  *
356  * \param __begin begin of hits
357  * \param __end end of hits
358  * \return model
359  */
360  template<class T>
361  const JModel& operator()(T __begin, T __end) const
362  {
363  using namespace std;
364  using namespace JPP;
365  using namespace JGEOMETRY;
366 
367 #ifdef THREAD_SAFE
368  static mutex mtx; // TDecompSVD
369 
370  unique_lock<mutex> lock(mtx);
371 #endif
372  value = JModel(__begin, __end); // set up global model according data
373 
374  if (this->option == JMODEL::fit_emitters_only_t)
375  value.setOption(JMODEL::fit_emitters_only_t);
376  else
378 
379  H_t H; // H-equation as per hit
380  I_t i; // indices of H-equation in global model
381 
382  const size_t N = value.getN();
383 
384  TMatrixD V(N,N); // V = H^T x H
385  vector<double> Y(N, 0.0); // Y = H^T x y
386 
387  for (T hit = __begin; hit != __end; ++hit) {
388 
389  const JString& string = detector[hit->getString()];
390  const JPosition3D position = string.getPosition(hit->getFloor());
391  const double Vi = velocity.getInverseVelocity(hit->getDistance(position), hit->getZ(), position.getZ());
392 
393  const double h1 = string.getHeight(hit->getFloor());
394  const JPosition3D p1 = string.getPosition() - hit->getPosition();
395  const double ds = sqrt(p1.getLengthSquared() + h1*h1 + 2.0*p1.getZ()*h1);
396  const double y = hit->getValue() - Vi*ds;
397  const double W = sqrt(hit->getWeight());
398 
399  H.t1 = W * 1.0;
400  H.tx = W * Vi * p1.getX() * h1 / ds;
401  H.ty = W * Vi * p1.getY() * h1 / ds;
402 
403  i.t1 = value.getIndex(hit->getEKey(), &H_t::t1);
404  i.tx = value.getIndex(hit->getString(), &H_t::tx);
405  i.ty = value.getIndex(hit->getString(), &H_t::ty);
406 
407  V(i.t1, i.t1) += H.t1 * H.t1;
408 
409  Y[i.t1] += W * H.t1 * y;
410 
411  if (hit->getFloor() != 0) {
412 
413  if (this->option != JMODEL::fit_emitters_only_t) {
414 
415  V(i.t1, i.tx) += H.t1 * H.tx; V(i.t1, i.ty) += H.t1 * H.ty;
416  V(i.tx, i.t1) += H.tx * H.t1; V(i.ty, i.t1) += H.ty * H.t1;
417 
418  V(i.tx, i.tx) += H.tx * H.tx; V(i.tx, i.ty) += H.tx * H.ty;
419  V(i.ty, i.tx) += H.ty * H.tx; V(i.ty, i.ty) += H.ty * H.ty;
420 
421  Y[i.tx] += W * H.tx * y;
422  Y[i.ty] += W * H.ty * y;
423  }
424  }
425  }
426 
427  // evaluate (H^T x H)^-1 x H^T x y
428 
429  TDecompSVD svd(V, TOLERANCE);
430 
431  Bool_t status;
432 
433  V = svd.Invert(status);
434 
435  for (size_t row = 0; row != N; ++row) {
436  for (size_t col = 0; col != N; ++col) {
437  value[row] += V(row,col) * Y[col];
438  }
439  }
440 
441  return value;
442  }
443 
444  private:
445  mutable JModel value;
446 #ifdef THREAD_SAFE
447  static std::mutex mtx; // TDecompSVD
448 #endif
449  };
450 
451 #ifdef THREAD_SAFE
452  /**
453  * Declaration of mutex.
454  */
455  std::mutex JKatoomba<JEstimator>::mtx;
456 #endif
457 
458  /**
459  * Template specialisation of fit function of acoustic model based on JSimplex minimiser.
460  */
461  template<>
462  struct JKatoomba<JSimplex> :
463  public JSimplex<JModel>,
464  public JKatoomba<>
465  {
466  /**
467  * Constructor
468  *
469  * \param detector detector
470  * \param velocity sound velocity
471  * \param option option
472  */
474  const JSoundVelocity& velocity,
475  const int option) :
476  JKatoomba<>(detector, velocity, option)
477  {};
478 
479 
480  /**
481  * Fit function.\n
482  * This method is used to determine the chi2 of given hit with respect to actual model.
483  *
484  * \param model model
485  * \param hit hit
486  * \return chi2
487  */
488  template<class JPDF_t>
489  double operator()(const JModel& model, const JHit<JPDF_t>& hit) const
490  {
491  const double toa_s = this->getToA(model, hit);
492  const double u = (toa_s - hit.getValue()) / hit.sigma;
493 
494  return estimator->getRho(u);
495  }
496 
497 
498  /**
499  * Fit.
500  *
501  * \param __begin begin of hits
502  * \param __end end of hits
503  * \return chi2
504  */
505  template<class T>
506  double operator()(T __begin, T __end)
507  {
508  this->step.clear();
509 
510  for (JModel::string_type::const_iterator i = this->value.string.begin(); i != this->value.string.end(); ++i) {
511 
512  JModel model;
513 
514  switch (this->option) {
515 
517  model.string[i->first] = JMODEL::JString(0.0, 0.0, 0.0, 0.0, 1.0e-6); this->step.push_back(model);
518 
520  model.string[i->first] = JMODEL::JString(0.0, 0.0, 1.0e-6, 0.0, 0.0); this->step.push_back(model);
521  model.string[i->first] = JMODEL::JString(0.0, 0.0, 0.0, 1.0e-6, 0.0); this->step.push_back(model);
522 
524  model.string[i->first] = JMODEL::JString(1.0e-3, 0.0, 0.0, 0.0, 0.0); this->step.push_back(model);
525  model.string[i->first] = JMODEL::JString(0.0, 1.0e-3, 0.0, 0.0, 0.0); this->step.push_back(model);
526  break;
527 
528  default:
529  break;
530  }
531  }
532 
533  for (JModel::emitter_type::const_iterator i = this->value.emitter.begin(); i != this->value.emitter.end(); ++i) {
534 
535  JModel model;
536 
537  model.emitter[i->first] = JMODEL::JEmitter(5.0e-5); this->step.push_back(model);
538  }
539 
540  return JSimplex<JModel>::operator()(*this, __begin, __end);
541  }
542  };
543 
544 
545  template<class T>
546  class JGandalf {};
547 
548 
549  /**
550  * Template specialisation of fit function of acoustic model based on JGandalf minimiser.
551  */
552  template<>
553  struct JKatoomba<JGandalf> :
554  public JKatoomba<>
555  {
556  typedef double result_type;
557 
558 
559  /**
560  * Constructor
561  *
562  * \param detector detector
563  * \param velocity sound velocity
564  * \param option option
565  */
567  const JSoundVelocity& velocity,
568  const int option) :
569  JKatoomba<>(detector, velocity, option)
570  {};
571 
572 
573  /**
574  * Fit.
575  *
576  * \param __begin begin of hits
577  * \param __end end of hits
578  * \return chi2 and gradient
579  */
580  template<class T>
581  result_type operator()(T __begin, T __end)
582  {
583  using namespace std;
584  using namespace JPP;
585 
586  value.setOption(this->option);
587 
588  const int N = value.getN();
589 
590  V.resize(N);
591  Y.resize(N);
592  h.resize(N);
593 
594  lambda = LAMBDA_MIN;
595 
596  result_type precessor = numeric_limits<double>::max();
597 
598  for (numberOfIterations = 0; numberOfIterations != MAXIMUM_ITERATIONS; ++numberOfIterations) {
599 
600  DEBUG("step: " << numberOfIterations << endl);
601 
602  evaluate(__begin, __end);
603 
604  DEBUG("lambda: " << FIXED(12,5) << lambda << endl);
605  DEBUG("chi2: " << FIXED(12,5) << successor << endl);
606 
607  if (successor < precessor) {
608 
609  if (numberOfIterations != 0) {
610 
611  if (fabs(precessor - successor) < EPSILON*fabs(precessor)) {
612  return successor;
613  }
614 
615  if (lambda > LAMBDA_MIN) {
616  lambda /= LAMBDA_DOWN;
617  }
618  }
619 
620  precessor = successor;
621  previous = value;
622 
623  } else {
624 
625  value = previous;
626  lambda *= LAMBDA_UP;
627 
628  if (lambda > LAMBDA_MAX) {
629  return precessor; // no improvement found
630  }
631 
632  evaluate(__begin, __end);
633  }
634 
635 
636  // force definite positiveness
637 
638  for (int i = 0; i != N; ++i) {
639 
640  if (V(i,i) < PIVOT) {
641  V(i,i) = PIVOT;
642  }
643 
644  h[i] = 1.0 / sqrt(V(i,i));
645  }
646 
647 
648  // normalisation
649 
650  for (int i = 0; i != N; ++i) {
651  for (int j = 0; j != i; ++j) {
652  V(j,i) *= h[i] * h[j];
653  V(i,j) = V(j,i);
654  }
655  }
656 
657  for (int i = 0; i != N; ++i) {
658  V(i,i) = 1.0 + lambda;
659  }
660 
661 
662  try {
663  V.invert();
664  }
665  catch (const JException& error) {
666  ERROR("JKatoomb<JGandalf>: " << error.what() << endl);
667  return precessor;
668  }
669 
670 
671  for (int i = 0; i != N; ++i) {
672 
673  DEBUG("u[" << noshowpos << setw(3) << i << "] = " << showpos << FIXED(20,5) << value[i]);
674 
675  double y = 0.0;
676 
677  for (int j = 0; j != N; ++j) {
678  y += V(i,j) * Y[j] * h[i] * h[j];
679  }
680 
681  value[i] -= y;
682 
683  DEBUG(' ' << FIXED(20,10) << y << noshowpos << endl);
684  }
685  }
686 
687  return precessor;
688  }
689 
690  static int debug; //!< debug level
691  static int MAXIMUM_ITERATIONS; //!< maximal number of iterations
692  static double EPSILON; //!< maximal distance to minimum
693  static double LAMBDA_MIN; //!< minimal value control parameter
694  static double LAMBDA_MAX; //!< maximal value control parameter
695  static double LAMBDA_UP; //!< multiplication factor control parameter
696  static double LAMBDA_DOWN; //!< multiplication factor control parameter
697  static double PIVOT; //!< minimal value diagonal element of matrix
698 
699  double lambda;
703 
704  private:
705  /**
706  * Evaluation of fit.
707  *
708  * \param __begin begin of data
709  * \param __end end of data
710  */
711  template<class T>
712  inline void evaluate(T __begin, T __end)
713  {
714  using namespace std;
715  using namespace JPP;
716 
717  successor = 0.0;
718 
719  V.reset();
720 
721  for (std::vector<double>::iterator i = Y.begin(); i != Y.end(); ++i) {
722  *i = 0.0;
723  }
724 
725  for (T hit = __begin; hit != __end; ++hit) {
726 
727  const JGEOMETRY::JString& string = detector[hit->getString()];
728  const JMODEL ::JString& parameters = value.string[hit->getString()];
729  const JPosition3D position = string.getPosition(parameters, hit->getFloor());
730 
731  const double D = hit->getDistance(position);
732  const double Vi = velocity.getInverseVelocity(D, hit->getZ(), position.getZ());
733  const double toa_s = value.emitter[hit->getEKey()].t1 + D * Vi;
734 
735  const double u = (toa_s - hit->getValue()) / hit->sigma;
736  const double W = sqrt(hit->getWeight());
737 
738  successor += (W*W) * estimator->getRho(u);
739 
740  H_t H(1.0, string.getGradient(parameters, hit->getPosition(), hit->getFloor()) * Vi);
741 
742  H *= W * estimator->getPsi(u) / hit->sigma;
743 
744  I_t i;
745 
746  i.t1 = value.getIndex(hit->getEKey(), &H_t::t1);
747  i.tx = value.getIndex(hit->getString(), &H_t::tx);
748  i.ty = value.getIndex(hit->getString(), &H_t::ty);
749  i.tx2 = value.getIndex(hit->getString(), &H_t::tx2);
750  i.ty2 = value.getIndex(hit->getString(), &H_t::ty2);
751  i.vs = value.getIndex(hit->getString(), &H_t::vs);
752 
753  V(i.t1, i.t1) += H.t1 * H.t1;
754 
755  Y[i.t1] += W * H.t1;
756 
757  if (hit->getFloor() != 0) {
758 
759  switch (this->option) {
760 
762  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;
763 
764  V(i.vs, i.t1) += H.vs * H.t1;
765  V(i.vs, i.tx) += H.vs * H.tx;
766  V(i.vs, i.ty) += H.vs * H.ty;
767  V(i.vs, i.tx2) += H.vs * H.tx2;
768  V(i.vs, i.ty2) += H.vs * H.ty2;
769 
770  V(i.vs, i.vs) += H.vs * H.vs;
771 
772  Y[i.vs] += W * H.vs;
773 
775  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;
776 
777  V(i.tx2, i.t1) += H.tx2 * H.t1;
778  V(i.tx2, i.tx) += H.tx2 * H.tx;
779  V(i.tx2, i.ty) += H.tx2 * H.ty;
780 
781  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;
782 
783  V(i.ty2, i.t1) += H.ty2 * H.t1;
784  V(i.ty2, i.tx) += H.ty2 * H.tx;
785  V(i.ty2, i.ty) += H.ty2 * H.ty;
786 
787  V(i.tx2, i.tx2) += H.tx2 * H.tx2; V(i.tx2, i.ty2) += H.tx2 * H.ty2;
788  V(i.ty2, i.tx2) += H.ty2 * H.tx2; V(i.ty2, i.ty2) += H.ty2 * H.ty2;
789 
790  Y[i.tx2] += W * H.tx2;
791  Y[i.ty2] += W * H.ty2;
792 
794  V(i.t1, i.tx) += H.t1 * H.tx; V(i.t1, i.ty) += H.t1 * H.ty;
795  V(i.tx, i.t1) += H.tx * H.t1; V(i.ty, i.t1) += H.ty * H.t1;
796 
797  V(i.tx, i.tx) += H.tx * H.tx; V(i.tx, i.ty) += H.tx * H.ty;
798  V(i.ty, i.tx) += H.ty * H.tx; V(i.ty, i.ty) += H.ty * H.ty;
799 
800  Y[i.tx] += W * H.tx;
801  Y[i.ty] += W * H.ty;
802  break;
803 
804  default:
805  break;
806  }
807  }
808  }
809  }
810 
811 
812  std::vector<double> Y; // gradient
815  std::vector<double> h; // normalisation vector
816  };
817 
818 
819  /**
820  * debug level.
821  */
823 
824 
825  /**
826  * maximal number of iterations.
827  */
829 
830 
831  /**
832  * maximal distance to minimum.
833  */
834  double JKatoomba<JGandalf>::EPSILON = 1.0e-3;
835 
836 
837  /**
838  * minimal value control parameter
839  */
840  double JKatoomba<JGandalf>::LAMBDA_MIN = 0.01;
841 
842 
843  /**
844  * maximal value control parameter
845  */
846  double JKatoomba<JGandalf>::LAMBDA_MAX = 100.0;
847 
848 
849  /**
850  * multiplication factor control parameter
851  */
852  double JKatoomba<JGandalf>::LAMBDA_UP = 9.0;
853 
854 
855  /**
856  * multiplication factor control parameter
857  */
858  double JKatoomba<JGandalf>::LAMBDA_DOWN = 11.0;
859 
860 
861  /**
862  * minimal value diagonal element of matrix
863  */
864  double JKatoomba<JGandalf>::PIVOT = 1.0e-3;
865 
866 
867  /**
868  * Worker class for fit function of acoustic model.
869  */
870  struct JKatoomba_t {
871  /**
872  * Constructor.
873  *
874  * \param geometry detector geometry
875  * \param V sound velocity
876  * \param parameters parameters
877  */
878  JKatoomba_t(const JGeometry& geometry, const JSoundVelocity& V, const JFitParameters& parameters) :
879  parameters(parameters),
880  estimator(geometry, V, parameters.option),
881  evaluator(geometry, V, parameters.option),
882  gandalf (geometry, V, parameters.option)
883  {
884  using namespace JPP;
885 
886  evaluator.estimator.reset(getMEstimator(parameters.mestimator));
887  gandalf .estimator.reset(getMEstimator(parameters.mestimator));
888  }
889 
890 
891  /**
892  * Auxiliary data structure for return value of fit.
893  */
894  template<class T>
895  struct result_type {
897  double chi2;
898  double ndf;
901 
902  /**
903  * Get size of data.
904  *
905  * \return size
906  */
907  int size() const
908  {
909  return std::distance(begin, end);
910  }
911 
912  /**
913  * Get time range of data.
914  *
915  * \return time range
916  */
918  {
920 
921  for (T i = begin; i != end; ++i) {
922  range.include(i->getValue());
923  }
924 
925  return range;
926  }
927  };
928 
929 
930  /**
931  * Fit.
932  *
933  * \param __begin begin of data
934  * \param __end end of data
935  * \return fit result
936  */
937  template<class T>
938  result_type<T> operator()(T __begin, T __end)
939  {
940  result = estimator(__begin, __end); // pre-fit
941 
942  if (parameters.stdev > 0.0) { // remove outliers
943 
944  for (double chi2 = evaluator(result, __begin, __end), chi2_old = chi2; ; ) {
945 
946  T out = __end;
947 
948  double xmax = 0.0;
949 
950  for (T hit = __begin; hit != __end; ++hit) {
951 
952  const double x = fabs(hit->getValue() - estimator.getToA(result, *hit)) / hit->sigma;
953 
954  if (x > xmax) {
955  xmax = x;
956  out = hit;
957  }
958  }
959 
960  if (xmax > parameters.stdev) {
961 
962  iter_swap(out, --__end);
963 
964  result = estimator(__begin, __end);
965  chi2 = evaluator(result, __begin, __end);
966 
967  if (chi2_old - chi2 > parameters.stdev * parameters.stdev) {
968 
969  chi2_old = chi2;
970 
971  continue;
972 
973  } else {
974 
975  iter_swap(out, __end++);
976 
977  result = estimator(__begin, __end);
978  chi2 = chi2_old;
979  }
980  }
981 
982  break;
983  }
984  }
985 
986  gandalf.value = result; // start value
987 
988  const double chi2 = gandalf (__begin, __end) / gandalf.estimator->getRho(1.0);
989  const double ndf = getWeight(__begin, __end) - gandalf.value.getN();
990 
991  return { gandalf.value, chi2, ndf, __begin, __end }; // return values
992  }
993 
994 
999 
1000  private:
1002  };
1003 }
1004 
1005 #endif
Worker class for fit function of acoustic model.
Definition: JKatoomba.hh:870
fit times of emission of emitters and tilt angles of strings with second order correction and stretch...
const double xmax
Definition: JQuadrature.cc:24
static const JRange< double, std::less< double > > DEFAULT_RANGE
Default range.
Definition: JRange.hh:555
Linear fit methods.
Acoustic hit.
JFit_t
Enumeration for fit algorithms.
Definition: JKatoomba.hh:105
I_t()
Default constructor.
Definition: JKatoomba.hh:247
General exception.
Definition: JException.hh:23
JKatoomba(const JDetector &detector, const JSoundVelocity &velocity, const int option)
Constructor.
Definition: JKatoomba.hh:286
General purpose data regression method.
JString & mul(const double factor)
Scale string.
Wrapper class around STL string class.
Definition: JString.hh:27
Auxiliary base class for aritmetic operations of derived class types.
Definition: JMath.hh:109
double getQ() const
Get quality.
TPaveText * p1
int size() const
Get size of data.
Definition: JKatoomba.hh:907
JKatoomba< JEstimator > estimator
Definition: JKatoomba.hh:996
double getWeight(T __begin, T __end)
Get total weight of data points.
Definition: JKatoomba.hh:59
Sound velocity.
int getFloor() const
Get floor number.
Definition: JLocation.hh:145
static double LAMBDA_UP
multiplication factor control parameter
Definition: JKatoomba.hh:695
range_type & include(argument_type x)
Include given value to range.
Definition: JRange.hh:397
Interface for maximum likelihood estimator (M-estimator).
Definition: JMEstimator.hh:20
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Detector data structure.
Definition: JDetector.hh:89
JKatoomba< JAbstractMinimiser > evaluator
Definition: JKatoomba.hh:997
result_type operator()(const JFunction_t &fit, T __begin, T __end)
Get chi2.
Definition: JRegressor.hh:46
Template specialisation of fit function of acoustic model based on linear approximation.
Definition: JKatoomba.hh:333
Acoustic geometries.
JKatoomba_t(const JGeometry &geometry, const JSoundVelocity &V, const JFitParameters &parameters)
Constructor.
Definition: JKatoomba.hh:878
then JShowerPostfit f $INPUT_FILE o $OUTPUT_FILE N
static const double H
Planck constant [eV s].
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
Definition: diff-Tuna.sh:38
JEmitter & mul(const double factor)
Scale emitter.
const JDetector & detector
detector
Definition: JKatoomba.hh:187
double operator()(const JModel &model, const JHit< JPDF_t > &hit) const
Fit function.
Definition: JKatoomba.hh:489
H_t & mul(const double factor)
Scale H-equation.
Definition: JKatoomba.hh:229
Auxiliary class for a type holder.
Definition: JType.hh:19
then fatal Wrong number of arguments fi set_variable STRING $argv[1] set_variable DETECTORXY_TXT $WORKDIR $DETECTORXY_TXT tail read X Y CHI2 RMS printf optimum n $X $Y $CHI2 $RMS awk v Y
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
JEKey getEKey() const
Get emitter hash key of this hit.
Abstract minimiser.
Definition: JRegressor.hh:25
Template definition of linear fit.
Definition: JEstimator.hh:25
H_t()
Default constructor.
Definition: JKatoomba.hh:204
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition: JVector3D.hh:270
Auxiliary data structure for return value of fit.
Definition: JKatoomba.hh:895
fit times of emission of emitters and tilt angles of strings
JKatoomba(const JDetector &detector, const JSoundVelocity &velocity, const int option)
Constructor.
Definition: JKatoomba.hh:566
Acoustics hit.
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
Acoustic fit parameters.
Linear fit.
Definition: JKatoomba.hh:106
Model for fit to acoustics data.
JFitParameters parameters
Definition: JKatoomba.hh:995
Template specialisation of fit function of acoustic model based on JGandalf minimiser.
Definition: JKatoomba.hh:553
Acoustic transmission.
static double EPSILON
maximal distance to minimum
Definition: JKatoomba.hh:692
Detector file.
Definition: JHead.hh:226
result_type operator()(T __begin, T __end)
Fit.
Definition: JKatoomba.hh:581
std::vector< double > vs
double getToE(const JModel &model, const JHit< JPDF_t > &hit) const
Get estimated time-of-emission for given hit.
Definition: JKatoomba.hh:174
Acoustic emitter.
Definition: JEmitter.hh:27
static double LAMBDA_MIN
minimal value control parameter
Definition: JKatoomba.hh:693
void evaluate(T __begin, T __end)
Evaluation of fit.
Definition: JKatoomba.hh:712
static int debug
debug level
Definition: JKatoomba.hh:690
return result
Definition: JPolint.hh:764
do set_variable OUTPUT_DIRECTORY $WORKDIR T
double operator()(T __begin, T __end)
Fit.
Definition: JKatoomba.hh:506
JKatoomba(const JDetector &detector, const JSoundVelocity &velocity, const int option)
Constructor.
Definition: JKatoomba.hh:136
#define ERROR(A)
Definition: JMessage.hh:66
double getToA(const JModel &model, const JHit< JPDF_t > &hit) const
Get estimated time-of-arrival for given hit.
Definition: JKatoomba.hh:153
static struct JACOUSTICS::@4 compare
Auxiliary data structure to sort transmissions.
result_type operator()(const JModel &model, const JHit< JPDF_t > &hit) const
Fit function.
Definition: JKatoomba.hh:302
double getY() const
Get y position.
Definition: JVector3D.hh:104
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
N x N symmetric matrix.
Definition: JMatrixNS.hh:27
fit times of emission of emitters and tilt angles of strings with second order correction ...
double getLengthSquared() const
Get length squared.
Definition: JVector3D.hh:235
General purpose messaging.
Implementation for depth dependend velocity of sound.
static double LAMBDA_MAX
maximal value control parameter
Definition: JKatoomba.hh:694
static int MAXIMUM_ITERATIONS
maximal number of iterations
Definition: JKatoomba.hh:691
JKatoomba(const JDetector &detector, const JSoundVelocity &velocity, const int option)
Constructor.
Definition: JKatoomba.hh:473
std::vector< double > h
Definition: JKatoomba.hh:815
double operator()(const JFunction_t &fit, T __begin, T __end)
Multi-dimensional fit.
Definition: JSimplex.hh:71
JTimeRange getTimeRange() const
Get time range of data.
Definition: JKatoomba.hh:917
Emitter hash key.
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable TRIPOD $argv[2] set_variable TX $argv[3] set_variable TY $argv[4] if[[!-f $DETECTOR]]
Definition: JFootprint.sh:28
JACOUSTICS::JModel::string_type string
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
JSoundVelocity velocity
sound velocity
Definition: JKatoomba.hh:188
std::shared_ptr< JMEstimator > estimator_type
Definition: JKatoomba.hh:126
JKatoomba(const JDetector &detector, const JSoundVelocity &velocity, const int option)
Constructor.
Definition: JKatoomba.hh:345
JACOUSTICS::JModel::emitter_type emitter
int getString() const
Get string number.
Definition: JLocation.hh:134
virtual const char * what() const override
Get error message.
Definition: JException.hh:48
const JModel & operator()(T __begin, T __end) const
Get start values of string parameters.
Definition: JKatoomba.hh:361
static double LAMBDA_DOWN
multiplication factor control parameter
Definition: JKatoomba.hh:696
Simple fit method based on Powell&#39;s algorithm, see reference: Numerical Recipes in C++...
Definition: JSimplex.hh:42
Template specialisation of fit function of acoustic model based on JAbstractMinimiser minimiser...
Definition: JKatoomba.hh:272
result_type operator()(const JModel &model, T __begin, T __end)
Fit.
Definition: JKatoomba.hh:320
double getToA() const
Get calibrated time of arrival.
then if[[!-f $DETECTOR]] then JDetector sh $DETECTOR fi cat $WORKDIR trigger_parameters txt<< EOFtrigger3DMuon.enabled=1;trigger3DMuon.numberOfHits=5;trigger3DMuon.gridAngle_deg=1;ctMin=0.0;TMaxLocal_ns=15.0;EOF set_variable TRIGGEREFFICIENCY_TRIGGERED_EVENTS_ONLY INPUT_FILES=() for((i=1;$i<=$NUMBER_OF_RUNS;++i));do JSirene.sh $DETECTOR $JPP_DATA/genhen.km3net_wpd_V2_0.evt.gz $WORKDIR/sirene_ ${i}.root JTriggerEfficiency.sh $DETECTOR $DETECTOR $WORKDIR/sirene_ ${i}.root $WORKDIR/trigger_efficiency_ ${i}.root $WORKDIR/trigger_parameters.txt $JPP_DATA/PMT_parameters.txt INPUT_FILES+=($WORKDIR/trigger_efficiency_ ${i}.root) done for ANGLE_DEG in $ANGLES_DEG[*];do set_variable SIGMA_NS 3.0 set_variable OUTLIERS 3 set_variable OUTPUT_FILE $WORKDIR/matrix\[${ANGLE_DEG}\deg\].root $JPP_DIR/examples/JReconstruction-f"$INPUT_FILES[*]"-o $OUTPUT_FILE-S ${SIGMA_NS}-A ${ANGLE_DEG}-O ${OUTLIERS}-d ${DEBUG}--!fiif[[$OPTION=="plot"]];then if((0));then for H1 in h0 h1;do JPlot1D-f"$WORKDIR/matrix["${^ANGLES_DEG}" deg].root:${H1}"-y"1 2e3"-Y-L TR-T""-\^"number of events [a.u.]"-> o chi2
Definition: JMatrixNZ.sh:106
Acoustic transmission.
static double PIVOT
minimal value diagonal element of matrix
Definition: JKatoomba.hh:697
int getID() const
Get identifier.
double getX() const
Get x position.
Definition: JVector3D.hh:94
int mestimator
M-estimator.
int j
Definition: JPolint.hh:703
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
container_type::const_iterator const_iterator
Definition: JHashMap.hh:86
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
Definition: JMEstimator.hh:203
H_t(const JMODEL::JEmitter &emitter, const JMODEL::JString &string)
Constructor.
Definition: JKatoomba.hh:216
double u[N+1]
Definition: JPolint.hh:776
Model for fit to acoutsics data.
fit only times of emission of emitters
estimator_type estimator
M-Estimator function.
Definition: JKatoomba.hh:190
do echo Generating $dir eval D
Definition: JDrawLED.sh:53
double getZ() const
Get z position.
Definition: JVector3D.hh:115
Maximum likelihood estimator (M-estimators).
Simplex fit.
Definition: JKatoomba.hh:107
Template definition of fit function of acoustic model.
Definition: JKatoomba.hh:116
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
JKatoomba< JGandalf > gandalf
Definition: JKatoomba.hh:998
Gandalf fit.
Definition: JKatoomba.hh:108
result_type< T > operator()(T __begin, T __end)
Fit.
Definition: JKatoomba.hh:938
std::vector< double > Y
Definition: JKatoomba.hh:812