#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include <map>
#include <algorithm>
#include "TRandom3.h"
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "JLang/JPredicate.hh"
#include "JLang/JSTDToolkit.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JTripod.hh"
#include "JDetector/JTransmitter.hh"
#include "JDetector/JHydrophone.hh"
#include "JDetector/JModule.hh"
#include "JROOT/JManager.hh"
#include "JROOT/JRootToolkit.hh"
#include "JAcoustics/JSoundVelocity.hh"
#include "JAcoustics/JEmitter.hh"
#include "JAcoustics/JAcousticsToolkit.hh"
#include "JAcoustics/JAcousticsSupportkit.hh"
#include "JAcoustics/JHit.hh"
#include "JAcoustics/JFitParameters.hh"
#include "JAcoustics/JKatoomba.hh"
#include "JAcoustics/JTransmission_t.hh"
#include "Jeep/JContainer.hh"
#include "Jeep/JTimer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
namespace {
using namespace JPP;
static const int WILDCARD = -1; // wild card for emitter identifier
static const int NUMBER_OF_PINGS = 11; // default number of pings
static const JWaveform EMITTER_WAVEFORM(2.5e6, // default emitter power [quality]
30.0); // default emitter frequency [kHz]
/**
* Enumeration for fit algorithms.
*/
enum JFit_t {
linear_t = 0, //!< linear fit
simplex_t, //!< simplex fit
gandalf_t //!< gandalf fit
};
}
/**
* \file
*
* Example application to test fit of model to simulated acoustic data.
* \author mdejong
*/
int main(int argc, char **argv)
{
using namespace std;
using namespace JPP;
typedef JContainer< vector<JTripod> > tripods_container;
typedef JContainer< vector<JTransmitter> > transmitters_container;
typedef JContainer< vector<JHydrophone> > hydrophones_container;
typedef JContainer< map<int, JWaveform> > waveform_container;
typedef JContainer< set<JTransmission_t> > disable_container;
string detectorFile;
string outputFile;
int numberOfEvents;
map<int, int> numberOfPings;
JSoundVelocity V = getSoundVelocity; // default sound velocity
tripods_container tripods; // tripods
transmitters_container transmitters; // transmitters
hydrophones_container hydrophones; // hydrophones
JTimeRange T_s(-1.0e-4, +1.0e-4); // time range of background [s]
JFitParameters parameters; // fit parameters
int fit; // type of fit
bool unify; // unify weighing of pings
disable_container disable; // disable tansmissions
waveform_container waveforms; // emitter power and frequency
double background;
UInt_t seed;
int debug;
numberOfPings[WILDCARD] = NUMBER_OF_PINGS;
waveforms [WILDCARD] = EMITTER_WAVEFORM;
try {
JParser<> zap("Example application to test fit of model to simulated acoustic data.");
zap['a'] = make_field(detectorFile);
zap['o'] = make_field(outputFile);
zap['n'] = make_field(numberOfEvents);
zap['@'] = make_field(parameters) = JPARSER::initialised();
zap['N'] = make_field(numberOfPings, "number of pings per emitter") = JPARSER::initialised();
zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
zap['T'] = make_field(tripods, "tripod data");
zap['Y'] = make_field(transmitters, "transmitter data") = JPARSER::initialised();
zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
zap['F'] = make_field(fit, "fit type") = linear_t, simplex_t, gandalf_t;
zap['u'] = make_field(unify, "unify weighing of pings");
zap['!'] = make_field(disable, "disable transmission") = JPARSER::initialised();
zap['W'] = make_field(waveforms, "emitter power and frequency") = JPARSER::initialised();
zap['B'] = make_field(background, "background probability") = 0.0;
zap['S'] = make_field(seed) = 0;
zap['d'] = make_field(debug) = 1;
zap(argc, argv);
}
catch(const exception &error) {
FATAL(error.what() << endl);
}
gRandom->SetSeed(seed);
if (numberOfEvents <= 0) { FATAL("Invalid number of events " << numberOfEvents << endl); }
JDetector detector;
try {
load(detectorFile, detector);
}
catch(const JException& error) {
FATAL(error);
}
vector<JEmitter> emitters;
for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
emitters.push_back(JEmitter(i->getID(),
i->getUTMPosition() - detector.getUTMPosition()));
}
for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) {
try {
emitters.push_back(JEmitter(i->getID(),
i->getPosition() + detector.getModule(i->getLocation()).getPosition()));
}
catch(const exception&) {
continue; // if no string available, discard transmitter
}
}
if (detector.empty()) { FATAL("No modules in detector." << endl); }
if (emitters.empty()) { FATAL("No emitters in system." << endl); }
V.set(detector.getUTMZ()); // sound velocity at detector depth
const double D_m = -detector.getUTMZ(); // depth [m]
const JGeometry geometry(detector, hydrophones);
DEBUG(geometry);
JKatoomba<JEstimator> estimator(geometry, V, parameters.option);
JKatoomba<JAbstractMinimiser> evaluator(geometry, V, parameters.option);
JKatoomba<JSimplex> simplex (geometry, V, parameters.option);
JKatoomba<JGandalf> gandalf (geometry, V, parameters.option);
simplex.estimator.reset(getMEstimator(parameters.mestimator));
gandalf.estimator.reset(getMEstimator(parameters.mestimator));
simplex.debug = debug;
gandalf.debug = debug;
JSimplex<JModel>::MAXIMUM_ITERATIONS = 10000;
TH1D h0("cpu", NULL, 100, 1.0, 7.0);
TH1D h1("chi2/NDF", NULL, 100, 0.0, 5.0);
JManager<int, TH1D> H1(new TH1D("emitter[%]", NULL, 100, -5.0e-5, +5.0e-5));
JManager<int, TH2D> H2(new TH2D("string[%]", NULL, 100, -1.0, +1.0, 100, -1.0, +1.0));
JTimer timer;
for (int number_of_events = 0, count = 0; number_of_events != numberOfEvents; ++number_of_events) {
STATUS("event: " << setw(10) << number_of_events << '\r'); DEBUG(endl);
JModel model; // randomised model data
for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
if (model.string[module->getString()] == JMODEL::JString()) {
model.string[module->getString()] = JMODEL::JString(gRandom->Uniform(-2.0e-2, +2.0e-2),
gRandom->Uniform(-2.0e-2, +2.0e-2));
}
}
DEBUG("Model" << endl << model << endl);
// set module positions according actual model data
for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
if (geometry.hasLocation(module->getLocation())) {
module->set(estimator.geometry[module->getString()].getPosition(model.string[module->getString()], module->getFloor()));
}
}
// generate hits
int minimum_number_of_pings = numeric_limits<int>::max();
for (vector<JEmitter>::const_iterator emitter = emitters.begin(); emitter != emitters.end(); ++emitter) {
minimum_number_of_pings = min(minimum_number_of_pings, getValue(numberOfPings, emitter->getID(), numberOfPings[WILDCARD]));
}
vector<JHit> data;
for (vector<JEmitter>::const_iterator emitter = emitters.begin(); emitter != emitters.end(); ++emitter) {
const int number_of_pings = getValue(numberOfPings, emitter->getID(), numberOfPings[WILDCARD]);
const double weight = (unify ? (double) minimum_number_of_pings / (double) number_of_pings : 1.0);
for (int ping_counter = 0; ping_counter != number_of_pings; ++ping_counter) {
++count;
const double toe_s = ping_counter * 10.0 + gRandom->Uniform(-1.0, +1.0);
model.emission[emitter->getID()][count] = toe_s;
for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
if (disable.count(JTransmission_t(emitter->getID(), module->getID())) == 0) {
if (geometry.hasLocation(module->getLocation())) {
const JWaveform& waveform = getValue(waveforms, emitter->getID(), waveforms[WILDCARD]);
const double d_m = getDistance(module->getPosition(), emitter->getPosition());
const double toa_s = toe_s + V.getTime(d_m, emitter->getZ(), module->getZ());
const double Q = waveform.getQ(D_m, d_m);
if (Q >= parameters.Qmin) {
double t1_s = toa_s;
if (gRandom->Rndm() >= background)
t1_s = gRandom->Gaus(toa_s, parameters.sigma_s);
else
t1_s = gRandom->Uniform(toa_s + T_s.getLowerLimit(),
toa_s + T_s.getUpperLimit());
const JHit hit(*emitter,
count,
module->getLocation(),
t1_s,
parameters.sigma_s,
weight);
DEBUG("hit: " << hit << ' ' << FIXED(7,1) << Q << endl);
data.push_back(hit);
}
}
}
}
}
}
timer.reset();
timer.start();
JModel result = estimator(data.begin(), data.end()); // pre-fit
double chi2 = numeric_limits<double>::max();
switch (fit) {
case linear_t:
chi2 = evaluator(result, data.begin(), data.end());
break;
case simplex_t:
simplex.value = result; // start value
chi2 = simplex(data.begin(), data.end()) / simplex.estimator->getRho(1.0);
result = simplex.value;
break;
case gandalf_t:
gandalf.value = result; // start value
chi2 = gandalf(data.begin(), data.end()) / gandalf.estimator->getRho(1.0);
result = gandalf.value;
break;
default:
break;
}
timer.stop();
double W = 0.0;
for (vector<JHit>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
W += hit->getWeight();
}
const int ndf = data.size() - result.getN();
DEBUG("Final values" << endl
<< FIXED(9,3) << chi2 << '/' << ndf << endl
<< result << endl);
h0.Fill(log10((double) timer.usec_wall));
h1.Fill(chi2 / (double) (W - result.getN()));
for (JHashMap<int, JMODEL::JString>::const_iterator i = model.string.begin(); i != model.string.end(); ++i) {
H2[i->first]->Fill((i->second.tx - result.string [i->first].tx) * 1.0e3, // mrad
(i->second.ty - result.string [i->first].ty) * 1.0e3); // mrad
}
for (JHashMap<JEKey, JMODEL::JEmission>::const_iterator i = model.emission.begin(); i != model.emission.end(); ++i) {
H1[i->first.getID()]->Fill(i->second.t1 - result.emission[i->first].t1);
}
}
STATUS(endl);
TFile out(outputFile.c_str(), "recreate");
out << h0
<< h1
<< H1
<< H2;
out.Write();
out.Close();
}