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