Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
examples/JAcoustics/JKatoomba.cc File Reference

Example application to test fit of model to simulated acoustic data. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <map>
#include <algorithm>
#include "TRandom3.h"
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "JLang/JPredicate.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JTripod.hh"
#include "JDetector/JModule.hh"
#include "JDetector/JHydrophone.hh"
#include "JROOT/JManager.hh"
#include "JROOT/JRootToolkit.hh"
#include "JAcoustics/JSoundVelocity.hh"
#include "JAcoustics/JEmitter.hh"
#include "JAcoustics/JAcousticsToolkit.hh"
#include "JAcoustics/JHit.hh"
#include "JAcoustics/JKatoomba.hh"
#include "Jeep/JProperties.hh"
#include "Jeep/JContainer.hh"
#include "Jeep/JTimer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Example application to test fit of model to simulated acoustic data.

Author
mdejong

Definition in file examples/JAcoustics/JKatoomba.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 44 of file examples/JAcoustics/JKatoomba.cc.

45 {
46  using namespace std;
47  using namespace JPP;
48 
49  typedef JContainer< vector<JTripod> > tripods_container;
50  typedef JContainer< vector<JHydrophone> > hydrophones_container;
51 
52  string detectorFile;
53  string outputFile;
54  int numberOfEvents;
55  map<int, int> numberOfPings;
56  JSoundVelocity V = getSoundVelocity; // default sound velocity
57  tripods_container tripods; // tripods
58  hydrophones_container hydrophones; // hydrophones
59  int mestimator = EM_LINEAR; // M-estimator
60  double sigma_s = 50.0e-6; // time-of-arrival resolution [s]
61  double background = 1.0e-4; // probability background
62  JTimeRange T_s(-1.0e-4, +1.0e-4); // time range of background [s]
63  int fit; // type of fit
64  bool unique; // one ping per cycle
65  UInt_t seed;
66  int debug;
67 
68  numberOfPings[-1] = 11; // default
69 
70  try {
71 
72  JProperties properties;
73 
74  properties.insert(gmake_property(mestimator));
75  properties.insert(gmake_property(sigma_s));
76  properties.insert(gmake_property(background));
77  properties.insert(gmake_property(T_s));
78 
79  JParser<> zap("Example application to test fit of model to simukated acoustic data.");
80 
81  zap['a'] = make_field(detectorFile);
82  zap['o'] = make_field(outputFile);
83  zap['n'] = make_field(numberOfEvents);
84  zap['@'] = make_field(properties) = JPARSER::initialised();
85  zap['N'] = make_field(numberOfPings, "number of pings per emitter") = JPARSER::initialised();
86  zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
87  zap['T'] = make_field(tripods, "tripod data");
88  zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
89  zap['F'] = make_field(fit) = linear_t, simplex_t, gandalf_t;
90  zap['u'] = make_field(unique, "one ping per cycle");
91  zap['S'] = make_field(seed) = 0;
92  zap['d'] = make_field(debug) = 1;
93 
94  zap(argc, argv);
95  }
96  catch(const exception &error) {
97  FATAL(error.what() << endl);
98  }
99 
100 
101  cout.tie(&cerr);
102 
103  gRandom->SetSeed(seed);
104 
105  if (numberOfEvents <= 0) { FATAL("Invalid number of events " << numberOfEvents << endl); }
106 
108 
109  try {
110  load(detectorFile, detector);
111  }
112  catch(const JException& error) {
113  FATAL(error);
114  }
115 
116  vector<JEmitter> emitters;
117 
118  for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
119  emitters.push_back(JEmitter(i->getID(), i->getUTMPosition() - detector.getUTMPosition()));
120  }
121 
122  if (detector.empty()) { FATAL("No modules in detector." << endl); }
123  if (emitters.empty()) { FATAL("No emitters in system." << endl); }
124 
125 
126  const JGeometry geometry(detector, hydrophones);
127 
128  DEBUG(geometry);
129 
130 
131  JKatoomba<JEstimator> estimator(geometry, V);
132  JKatoomba<JAbstractMinimiser> evaluator(geometry, V);
133  JKatoomba<JSimplex> simplex (geometry, V);
134  JKatoomba<JGandalf> gandalf (geometry, V);
135 
136  simplex.estimator.reset(getMEstimator(mestimator));
137  gandalf.estimator.reset(getMEstimator(mestimator));
138 
139  simplex.debug = debug;
140  gandalf.debug = debug;
141 
143 
144 
145  typedef JHit<JPDFGauss> hit_type;
146 
147 
148  TH1D h0("cpu", NULL, 100, 1.0, 7.0);
149  TH1D h1("chi2/NDF", NULL, 100, 0.0, 5.0);
150 
151  JManager<int, TH1D> H1(new TH1D("emitter[%]", NULL, 100, -5.0e-5, +5.0e-5));
152  JManager<int, TH2D> H2(new TH2D("string[%]", NULL, 100, -1.0, +1.0, 100, -1.0, +1.0));
153 
154 
155  JTimer timer;
156 
157  for (int number_of_events = 0, count = 0; number_of_events != numberOfEvents; ++number_of_events) {
158 
159  STATUS("event: " << setw(10) << number_of_events << '\r'); DEBUG(endl);
160 
161  JModel model; // randomised model data
162 
163  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
164 
165  if (model.string[module->getString()] == JMODEL::JString()) {
166 
167  model.string[module->getString()] = JMODEL::JString(gRandom->Uniform(-2.0e-2, +2.0e-2),
168  gRandom->Uniform(-2.0e-2, +2.0e-2));
169  }
170  }
171 
172  DEBUG(model);
173 
174  // set module positions according actual model data
175 
176  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
177  if (geometry.hasLocation(module->getLocation())) {
178  module->set(estimator.detector[module->getString()].getPosition(model.string[module->getString()], module->getFloor()));
179  }
180  }
181 
182 
183  // generate hits
184 
185  vector<hit_type> data;
186 
187  for (vector<JEmitter>::const_iterator emitter = emitters.begin(); emitter != emitters.end(); ++emitter) {
188 
189  int number_of_pings = numberOfPings[-1];
190 
191  if (numberOfPings.find(emitter->getID()) != numberOfPings.end()) {
192  number_of_pings = numberOfPings[emitter->getID()];
193  }
194 
195  int weight = numeric_limits<int>::max();
196 
197  for (map<int, int>::const_iterator i = numberOfPings.begin(); i != numberOfPings.end(); ++i) {
198  if (i->second < weight) {
199  weight = i->second;
200  }
201  }
202 
203  const double signal = (unique ? (double) weight / (double) number_of_pings : 1.0);
204 
205  for (int ping_counter = 0; ping_counter != number_of_pings; ++ping_counter) {
206 
207  ++count;
208 
209  const double toe_s = ping_counter * 10.0 + gRandom->Uniform(-1.0, +1.0);
210 
211  model.emitter[emitter->getID()][count] = toe_s;
212 
213  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
214 
215  if (geometry.hasLocation(module->getLocation())) {
216 
217  const double toa_s = toe_s + V.getTime(getDistance(module->getPosition(), emitter->getPosition()), emitter->getZ(), module->getZ());
218 
219  double t1_s = toa_s;
220 
221  if (gRandom->Rndm() >= background)
222  t1_s = gRandom->Gaus(toa_s, sigma_s);
223  else
224  t1_s = gRandom->Uniform(toa_s + T_s.getLowerLimit(),
225  toa_s + T_s.getUpperLimit());
226 
227  data.push_back(hit_type(*emitter,
228  count,
229  module->getLocation(),
230  JPDFGauss(t1_s, sigma_s, signal, background)));
231  }
232  }
233  }
234  }
235 
236  DEBUG("Model" << endl << model << endl);
237 
238  for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
239  DEBUG("hit: " << *hit << endl);
240  }
241 
242 
243  timer.reset();
244  timer.start();
245 
246  JModel result = estimator(data.begin(), data.end()); // pre-fit
247  double chi2 = numeric_limits<double>::max();
248 
249  switch (fit) {
250 
251  case linear_t:
252 
253  chi2 = evaluator(result, data.begin(), data.end());
254  break;
255 
256  case simplex_t:
257 
258  simplex.value = result; // start value
259 
260  chi2 = simplex(data.begin(), data.end()) / simplex.estimator->getRho(1.0);
261  result = simplex.value;
262  break;
263 
264  case gandalf_t:
265 
266  gandalf.value = result; // start value
267 
268  chi2 = gandalf(data.begin(), data.end()) / gandalf.estimator->getRho(1.0);
269  result = gandalf.value;
270  break;
271 
272  default:
273  break;
274  }
275 
276  timer.stop();
277 
278  double W = 0.0;
279 
280  for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
281  W += hit->getWeight();
282  }
283 
284  const int ndf = data.size() - model.getN();
285 
286  DEBUG("Final values" << endl
287  << FIXED(9,3) << chi2 << '/' << ndf << endl
288  << result << endl);
289 
290 
291  h0.Fill(log10((double) timer.usec_wall));
292  h1.Fill(chi2 / (double) (W - model.getN()));
293 
294  for (JHashMap<int, JMODEL::JString>::const_iterator i = model.string.begin(); i != model.string.end(); ++i) {
295  H2[i->first]->Fill((i->second.tx - result.string [i->first].tx) * 1.0e3, // mrad
296  (i->second.ty - result.string [i->first].ty) * 1.0e3); // mrad
297  }
298 
299  for (JHashMap<JEKey, JMODEL::JEmitter>::const_iterator i = model.emitter.begin(); i != model.emitter.end(); ++i) {
300  H1[i->first.getID()]->Fill(i->second.t1 - result.emitter[i->first].t1);
301  }
302  }
303  STATUS(endl);
304 
305 
306  TFile out(outputFile.c_str(), "recreate");
307 
308  out << h0
309  << h1
310  << H1
311  << H2;
312 
313  out.Write();
314  out.Close();
315 }
Utility class to parse command line options.
Definition: JParser.hh:1500
size_t getN() const
Get number of fit parameters.
General exception.
Definition: JException.hh:23
#define gmake_property(A)
macro to convert (template) parameter to JPropertiesElement object
do rm f tmp H1
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:80
Template specialisation of fit function of acoustic model based on linear approximation.
Definition: JKatoomba.hh:250
Template specialisation of fit function of acoustic model based on JSimplex minimiser.
Definition: JKatoomba.hh:345
Utility class to parse parameter values.
Definition: JProperties.hh:496
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:69
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:66
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:445
string outputFile
Acoustics hit.
double getDistance(const JFirst_t &first, const JSecond_t &second)
Get distance between objects.
Model for fit to acoustics data.
static const JSoundVelocity getSoundVelocity(1541.0,-17.0e-3,-2000.0)
Function object for velocity of sound.
Template specialisation of fit function of acoustic model based on JGandalf minimiser.
Definition: JKatoomba.hh:426
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
Definition: JManager.hh:43
Detector file.
Definition: JHead.hh:196
Acoustic emitter.
Definition: JEmitter.hh:27
Auxiliary wrapper for I/O of container with optional comment (see JComment).
Definition: JContainer.hh:39
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
return result
Definition: JPolint.hh:727
Auxiliary class for CPU timing and usage.
Definition: JTimer.hh:32
int debug
debug level
Definition: JSirene.cc:63
Implementation for velocity of sound.
#define FATAL(A)
Definition: JMessage.hh:67
JACOUSTICS::JModel::string_type string
JACOUSTICS::JModel::emitter_type emitter
then $DIR JKatoomba a $DETECTOR o $WORKDIR katoomba root T $TRIPOD n sigma_s
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
std::vector< int > count
Definition: JAlgorithm.hh:184
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:191
Custom probability density function of time-of-arrival.
container_type::const_iterator const_iterator
Definition: JHashMap.hh:85
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
Definition: JMEstimator.hh:166
do set_variable DETECTOR_TXT $WORKDIR detector
do for((RUN=${RANGE%%-*};$RUN<=${RANGE##*-};RUN+=1))
std::vector< double > weight
Definition: JAlgorithm.hh:428
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62