software/JAcoustics/JKatoomba.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <map>
#include <algorithm>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"

#include "JLang/JPredicate.hh"
#include "JLang/JComparator.hh"
#include "JLang/JComparison.hh"

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JTripod.hh"
#include "JDetector/JModule.hh"
#include "JDetector/JHydrophone.hh"

#include "JTools/JHashMap.hh"

#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMeta.hh"

#include "JAcoustics/JSoundVelocity.hh"
#include "JAcoustics/JEmitter.hh"
#include "JAcoustics/JAcousticsToolkit.hh"
#include "JAcoustics/JHit.hh"
#include "JAcoustics/JKatoomba.hh"
#include "JAcoustics/JEvent.hh"
#include "JAcoustics/JEvt.hh"
#include "JAcoustics/JEvtToolkit.hh"
#include "JAcoustics/JSupport.hh"

#include "Jeep/JProperties.hh"
#include "Jeep/JContainer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * \file
 *
 * Application to fit position calibration model to 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<JHydrophone> >  hydrophones_container;

  JMultipleFileScanner<JEvent>  inputFile;
  JFileRecorder<JTYPELIST<JEvt, JEvent, JMeta, TH1D, TH2D>::typelist>  outputFile;
  string                  detectorFile;
  JLimit_t&               numberOfEvents = inputFile.getLimit();
  JSoundVelocity          V              = getSoundVelocity;    // default sound velocity
  tripods_container       tripods;                              // tripods
  hydrophones_container   hydrophones;                          // hydrophones
  map<int, JMechanics>    mechanics;
  int                     mestimator     = EM_LINEAR;           // M-estimator
  double                  sigma_s        =  50.0e-6;            // time-of-arrival resolution [s]
  double                  background     =   1.0e-4;            // probability background transmission
  double                  Tmax_s         = 300.0;               // time window [s]
  size_t                  Nmin           =   3;                 // minimum number of emitters
  double                  stdev          =   5.0;               // standard deviation for outlier removal
  int                     fit;                                  // type of fit
  bool                    unique;                               // one ping per cycle
  int                     debug;

  try {

    JProperties properties;

    properties.insert(gmake_property(mestimator));
    properties.insert(gmake_property(sigma_s));
    properties.insert(gmake_property(background));
    properties.insert(gmake_property(Tmax_s));
    properties.insert(gmake_property(Nmin));
    properties.insert(gmake_property(stdev));

    JParser<> zap("Application to fit position calibration model to acoustic data.");

    zap['f'] = make_field(inputFile,       "output of JAcousticEventBuilder[.sh]");
    zap['o'] = make_field(outputFile);
    zap['n'] = make_field(numberOfEvents)                                    = JLimit::max();
    zap['a'] = make_field(detectorFile);
    zap['@'] = make_field(properties)                                        = JPARSER::initialised();
    zap['V'] = make_field(V,               "sound velocity")                 = JPARSER::initialised();
    zap['T'] = make_field(tripods,         "tripod data");
    zap['H'] = make_field(hydrophones,     "hydrophone data")                = JPARSER::initialised();
    zap['M'] = make_field(mechanics,       "mechanics data")                 = JPARSER::initialised();
    zap['F'] = make_field(fit,             "fit type")                       = linear_t, simplex_t, gandalf_t;
    zap['u'] = make_field(unique,          "one ping per cycle");
    zap['d'] = make_field(debug)                                             = 1;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  cout.tie(&cerr);

  JDetector detector;

  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }


  JHashMap<int, JLocation>   receivers;
  JHashMap<int, JEmitter>    emitters;

  for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
    receivers[i->getID()] = i->getLocation();
  }

  for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
    emitters[i->getID()]  = JEmitter(i->getID(),
                                     i->getUTMPosition() - detector.getUTMPosition());
  }


  JGeometry geometry(detector, hydrophones);

  for (JGeometry::iterator i = geometry.begin(); i != geometry.end(); ++i) {
    i->second.mechanics = mechanics[i->first];
  }
  
  DEBUG(geometry);


  JKatoomba<JEstimator>         estimator(geometry, V[detector.getUTMZ()]);
  JKatoomba<JAbstractMinimiser> evaluator(geometry, V[detector.getUTMZ()]);
  JKatoomba<JSimplex>           simplex  (geometry, V[detector.getUTMZ()]);
  JKatoomba<JGandalf>           gandalf  (geometry, V[detector.getUTMZ()]);

  simplex.estimator.reset(getMEstimator(mestimator));
  gandalf.estimator.reset(getMEstimator(mestimator));

  simplex.debug = debug;
  gandalf.debug = debug;

  JSimplex<JModel>::MAXIMUM_ITERATIONS = 10000;

  typedef JHit<JPDFGauss>  hit_type;

  TH1D h0("chi2/NDF", NULL, 5000,  0.0,  100.0);

  TH2D ha("ha", NULL, 
          getNumberOfStrings(detector) + 0, -0.5, getNumberOfStrings(detector) - 0.5,
          getNumberOfFloors (detector) + 1, -0.5, getNumberOfFloors (detector) + 0.5);
          
  TH2D hb("hb", NULL, 
          getNumberOfStrings(detector) + 0, -0.5, getNumberOfStrings(detector) - 0.5,
          getNumberOfFloors (detector) + 1, -0.5, getNumberOfFloors (detector) + 0.5);

  for (Int_t i = 1; i <= ha.GetXaxis()->GetNbins(); ++i) {
    ha.GetXaxis()->SetBinLabel(i, MAKE_CSTRING(geometry.at(i-1).first));
    hb.GetXaxis()->SetBinLabel(i, MAKE_CSTRING(geometry.at(i-1).first));
  } 

  // input data

  string oid;                                    // detector identifier

  vector<JEvent> zbuf;

  while (inputFile.hasNext()) {
    
    STATUS("input  " << setw(6) << inputFile.getCounter() << '\r'); DEBUG(endl);

    const JEvent* evt = inputFile.next();

    if (oid == "") {
      oid = evt->getOID();
    }

    if (oid != evt->getOID()) {                  // consistency check
      FATAL("Invalid detector identifier " << evt->getOID() << " != " << oid << endl);
    }
 
    zbuf.push_back(*evt);
  }
  STATUS(endl);

  sort(zbuf.begin(), zbuf.end());


  outputFile.open();

  outputFile.put(JMeta(argc, argv));


  for (vector<JEvent>::const_iterator p = zbuf.begin(), q; p != zbuf.end(); p = q) {
  
    STATUS("output " << setw(6) << distance(zbuf.cbegin(), p) << '\r'); DEBUG(endl);

    for (q = p; ++q != zbuf.end() && q->begin()->getToE() <= p->rbegin()->getToE() + Tmax_s; ) {}

    if (getNumberOfEmitters(p,q) >= Nmin) {

      map<int, int> numberOfPings;

      for (vector<JEvent>::const_iterator i = p; i != q; ++i) {
        numberOfPings[i->getID()] += 1;
      }

      int weight = numeric_limits<int>::max();

      for (map<int, int>::const_iterator i = numberOfPings.begin(); i != numberOfPings.end(); ++i) {
        DEBUG("Number of pings " << setw(2) << i->first << ' ' << setw(3) << i->second << endl);
      }

      for (map<int, int>::const_iterator i = numberOfPings.begin(); i != numberOfPings.end(); ++i) {
        if (i->second < weight) {
          weight = i->second;
        }
      }

      vector<hit_type> data;

      for (vector<JEvent>::const_iterator i = p; i != q; ++i) {

        {
          const JEmitter& emitter = emitters[i->getID()];

          const double signal = (unique ? (double) weight / (double) numberOfPings[i->getID()] : 1.0);

          // select first transmission

          map<int, vector<JTransmission> > buffer;

          for (JEvent::const_iterator hit = i->begin(); hit != i->end(); ++hit) {
            buffer[hit->getID()].push_back(*hit);
          }

          for (map<int, vector<JTransmission> >::iterator ps = buffer.begin(); ps != buffer.end(); ++ps) {

            if (receivers.has(ps->first) && geometry.hasLocation(receivers[ps->first])) {

              sort(ps->second.begin(), ps->second.end(), make_comparator(&JTransmission::getToA));

              data.push_back(hit_type(emitter,
                                      i->getCounter(),
                                      receivers[ps->first],
                                      JPDFGauss(ps->second.begin()->getToA(), sigma_s, signal, background)));
            }
          }
        }
      }
      
      for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
        DEBUG("hit: " << *hit << endl);
      }



      JModel result = estimator(data.begin(), data.end());              // pre-fit

      DEBUG("prefit:"                           << endl
            << result                           << endl);


      for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
        ha.Fill(geometry.getIndex(hit->getString()), hit->getFloor(), 1.0);
      }


      // remove outliers

      if (stdev > 0.0) {

        for ( ; ; ) {

          vector<hit_type>::iterator out = data.end();

          double xmax = 0.0;

          for (vector<hit_type>::iterator hit = data.begin(); hit != data.end(); ++hit) {

            const double x = fabs(hit->getValue() - estimator.getToA(result, *hit)) / hit->sigma;

            if (x > xmax) {
              xmax = x;
              out  = hit;
            }
          }

          if (xmax > stdev) {

            DEBUG("Remove outlier " << out->getLocation() << ' ' << xmax << endl);

            data.erase(out);

            result = estimator(data.begin(), data.end());

          } else {

            break;
          }
        }
      }

      for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
        hb.Fill(geometry.getIndex(hit->getString()), hit->getFloor(), 1.0);
      }


      double chi2 = numeric_limits<double>::max();

      switch (fit) {

      case linear_t:

        result = estimator(data.begin(), data.end());
        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;
      }

      double W = 0.0;

      for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
        W += hit->getWeight();
      }

      const int ndf = data.size() - result.getN();

      DEBUG("result:"                           << endl
            << FIXED(9,3) << chi2 << '/' << ndf << endl
            << result                           << endl);

      h0.Fill(chi2 / (W - result.getN()));
      
      // store results

      const JEvt evt = getEvt(JHead(oid, data.begin()->getValue(), data.rbegin()->getValue(), ndf, chi2), result);

      outputFile.put(evt);

      for (vector<JEvent>::const_iterator i = p; i != q; ++i) {

        JEvent out(*i);

        const double toe = result.emitter[JEKey(i->getID(), i->getCounter())].t1;

        for (JEvent::iterator hit = out.begin(); hit != out.end(); ++hit) {
          hit->setToE(toe);
        }

        outputFile.put(out);
      }
    }
  }
  STATUS(endl);

  JMultipleFileScanner<JMeta> io(inputFile);

  io >> outputFile;

  outputFile.put(h0);
  outputFile.put(ha);
  outputFile.put(hb);

  outputFile.close();
}