software/JAcoustics/JKatoomba.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include <map>
#include <deque>
#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 "JROOT/JManager.hh"

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JTripod.hh"
#include "JDetector/JTransmitter.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/JFitParameters.hh"
#include "JAcoustics/JKatoomba.hh"
#include "JAcoustics/JEvent.hh"
#include "JAcoustics/JEvt.hh"
#include "JAcoustics/JEvtToolkit.hh"
#include "JAcoustics/JSupport.hh"
#include "JAcoustics/JTransmission_t.hh"

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

namespace {

  /**
   * Get total weight of data points.
   *
   * \param  __begin        begin of data
   * \param  __end          end   of data
   * \return                total weight
   */
  template<class T>
  inline double getWeight(T __begin, T __end)
  {
    double W = 0;

    for (T i = __begin; i != __end; ++i) {
      W += i->getWeight();
    }

    return W;
  }
}


/**
 * \file
 *
 * Application to make a global fit of the detector geometry to acoustic data.\n
 * \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< set<JTransmission_t> > disable_container;

  JMultipleFileScanner<JEvent>  inputFile;
  JFileRecorder<JTYPELIST<JEvt, JEvent, JFitParameters, JMeta, TH1D, TH2D>::typelist>  outputFile;
  string                  detectorFile;
  JLimit_t&               numberOfEvents = inputFile.getLimit();
  JSoundVelocity          V              = getSoundVelocity;    // default sound velocity
  tripods_container       tripods;                              // tripods
  transmitters_container  transmitters;                         // transmitters
  hydrophones_container   hydrophones;                          // hydrophones
  JFitParameters          parameters;                           // fit parameters
  int                     fit;                                  // type of fit
  bool                    unify;                                // unify weighing of pings
  disable_container       disable;                              // disable tansmissions
  int                     debug;

  try {

    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(parameters)                                        = 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['M'] = make_field(getMechanics,    "mechanics 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['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);
  }

  const floor_range range = getRangeOfFloors(detector);

  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());
  }

  for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) {
    try {
      emitters[i->getID()]  = JEmitter(i->getID(),
                                       i->getPosition() + detector.getModule(i->getLocation()).getPosition());
    }
    catch(const exception&) {
      continue;                                                       // if no module available, discard transmitter
    }
  }

  V.set(detector.getUTMZ());                                          // sound velocity at detector depth

  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);

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

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

  JSimplex <JModel>  ::MAXIMUM_ITERATIONS = 10000;
  JKatoomba<JGandalf>::MAXIMUM_ITERATIONS = 10000;
  
  typedef JHit<JPDFGauss>  hit_type;

  TH1D h0("chi2/NDF", NULL, 50000,  0.0, 1000.0);
  TH1D h1("h1",       NULL,    51, -0.5,   50.5);
  TH1D hn("hn",       NULL,   100,  0.0,    6.0);

  JManager<int, TH2D> HA(new TH2D("ha[%]", NULL, 
                                  getNumberOfStrings(detector) + 0, -0.5, getNumberOfStrings(detector) - 0.5,
                                  range.getLength() + 1, range.getLowerLimit() - 0.5, range.getUpperLimit() + 0.5));
          
  JManager<int, TH2D> HB(new TH2D("hb[%]", NULL, 
                                  getNumberOfStrings(detector) + 0, -0.5, getNumberOfStrings(detector) - 0.5,
                                  range.getLength() + 1, range.getLowerLimit() - 0.5, range.getUpperLimit() + 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));
  } 

  string oid;                                      // detector identifier

  {                                                // sort input files
    map<double, string> zmap;

    for (const string& file_name : inputFile) {

      STATUS(file_name << '\r'); DEBUG(endl);

      for (JMultipleFileScanner<JEvent> in(file_name, 1); in.hasNext(); ) {

        const JEvent* evt = in.next();

        if      (oid == "")
          oid = evt->getOID();
        else if (oid != evt->getOID())             // consistency check
          FATAL("Invalid detector identifier " << evt->getOID() << " != " << oid << endl);
 
        if (!evt->empty()) {
          zmap[evt->begin()->getToE()] = file_name;
        }
      }
    }
    STATUS(endl);

    inputFile.clear();

    for (map<double, string>::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
      inputFile.push_back(i->second);
    }
  }


  outputFile.open();

  outputFile.put(JMeta(argc, argv));
  outputFile.put(parameters);
  
  int counter[] = { 0, 0 };

  typedef deque<JEvent> buffer_type;

  for (buffer_type zbuf; inputFile.hasNext(); ) {

    STATUS(inputFile.getFilename() << '\r'); DEBUG(endl);

    // read one file at a time
    
    for (const string file_name = inputFile.getFilename(); inputFile.hasNext() && file_name == inputFile.getFilename(); ) {

      const JEvent* evt = inputFile.next();

      if (emitters.has(evt->getID())) {
        zbuf.push_back(*evt);
      }
    }

    sort(zbuf.begin(), zbuf.end());                // sort according first time-of-emission

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

      if (q == zbuf.end()) {
        
        if (inputFile.hasNext()) {
          
          zbuf.erase(zbuf.begin(), p);             // remove processed data and continue reading

          break;
        }
      }
      
      if (getNumberOfEmitters(p,q) >= parameters.Nmin) {

        h1.Fill((double) getNumberOfEmitters(p,q));

        map<int, int> numberOfPings;

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

        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);
        }

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

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

        vector<hit_type> data;

        for (buffer_type::iterator i = p; i != q; ++i) {

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

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

            // select transmission with highest quality

            map<int, vector<JTransmission> > buffer;

            for (JEvent::const_iterator hit = i->begin(); hit != i->end(); ++hit) {
              if (disable.count(JTransmission_t(i->getID(), hit->getID())) == 0 &&
                  disable.count(JTransmission_t(-1,         hit->getID())) == 0) {
                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::getQ, JComparison::gt()));

                if (ps->second.begin()->getQ() >= parameters.Qmin) {

                  data.push_back(hit_type(emitter,
                                          distance(p,i),
                                          receivers[ps->first],
                                          JPDFGauss(ps->second.begin()->getToA(), parameters.sigma_s, signal, parameters.background)));
                }
              }
            }
          }
        }


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

        double chi2   = evaluator(result, data.begin(), data.end());
        double ndf    = getWeight(data.begin(), data.end())  -  result.getN();

        DEBUG("prefit:"                                          << ' '
              << FIXED(12,3) << chi2 << '/' << FIXED(6,1) << ndf << endl
              << result                                          << endl);

        for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
          DEBUG("hit: " << *hit << ' ' << FIXED(9,3) << evaluator(result, *hit) << endl);
        }

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


        // remove outliers

        if (parameters.stdev > 0.0) {

          double chi2_old =  evaluator(result, data.begin(), data.end());

          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 > parameters.stdev) {

              DEBUG("Remove outlier " << out->getEKey() << ' ' << out->getLocation() << ' ' << FIXED(12,3) << xmax << "?" << flush);

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

              iter_swap(out, --__end);

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

              double chi2_new =  evaluator(result, data.begin(), __end);

              if (chi2_old - chi2_new  >  parameters.stdev * parameters.stdev) {

                DEBUG(" yes " << FIXED(12,3) << (chi2_old - chi2_new) << endl);

                data.pop_back();

                chi2_old = chi2_new;

                continue;

              } else {

                DEBUG(" no " << endl);

                result = estimator(data.begin(), ++__end);
              }
            }
           
            break;
          }

          chi2   = evaluator(result, data.begin(), data.end());
          ndf    = getWeight(data.begin(), data.end())  -  result.getN();

          DEBUG("prefit:"                                          << ' '
                << FIXED(12,3) << chi2 << '/' << FIXED(6,1) << ndf << endl
                << result                                          << endl);

          for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
            DEBUG("hit: " << *hit << ' ' << FIXED(9,3) << evaluator(result, *hit) << endl);
          }
        }

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


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

        switch (fit) {

        case linear_t:

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

          hn.Fill(log10(1));
          break;

        case simplex_t:

          simplex.value = result;                                         // start value

          chi2   = simplex(data.begin(), data.end()) / simplex.estimator->getRho(1.0);
          result = simplex.value;

          hn.Fill(log10(simplex.numberOfIterations));
          break;

        case gandalf_t:

          gandalf.value = result;                                         // start value

          chi2   = gandalf(data.begin(), data.end()) / gandalf.estimator->getRho(1.0);
          result = gandalf.value;

          hn.Fill(log10(gandalf.numberOfIterations));
          break;

        default:
          break;
        }

        ndf = getWeight(data.begin(), data.end())  -  result.getN();

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

        for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
          DEBUG("hit: " << *hit << ' ' << FIXED(9,3) << evaluator(result, *hit) << endl);
        }

        h0.Fill(chi2 / ndf);
      
        // store results

        if (chi2 / ndf <= parameters.chi2perNDF) {
      
          const JEvt evt = getEvt(JHead(oid,
                                        data. begin()->getValue(),
                                        data.rbegin()->getValue(),
                                        data.size(),
                                        result.getN(),
                                        ndf,
                                        chi2),
                                  result);

          outputFile.put(evt);

          for (buffer_type::iterator i = p; i != q; ++i) {

            JEvent out(*i);                                                 // event with fitted times of emission

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

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

            outputFile.put(out);
          }

          counter[0] += 1;
        
        } else {

          WARNING(endl << "Event not written: " << FIXED(12,3) << chi2 << '/' << FIXED(6,1) << ndf << endl);

          counter[1] += 1;
        }
      }
    }
  }
  STATUS(endl);

  STATUS("Number of events written / rejected: " << counter[0] << " / " << counter[1] << endl);

  JMultipleFileScanner<JMeta> io(inputFile);

  io >> outputFile;

  outputFile.put(h0);
  outputFile.put(h1);
  outputFile.put(hn);

  for (JManager<int, TH2D>::const_iterator i = HA.begin(); i != HA.end(); ++i) { outputFile.put(*(i->second)); }
  for (JManager<int, TH2D>::const_iterator i = HB.begin(); i != HB.end(); ++i) { outputFile.put(*(i->second)); }

  outputFile.close();
}