JFremantle.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 <type_traits>
#include <functional>
#include <future>
#include <mutex>
#include <thread>
#include <vector>
#include <queue>

#include "TROOT.h"
#include "TFile.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/JTransmitter.hh"
#include "JDetector/JModule.hh"
#include "JDetector/JHydrophone.hh"

#include "JTools/JHashMap.hh"
#include "JTools/JRange.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 {

  using namespace std;
  using namespace JPP;
  
  typedef JHit<JPDFGauss>  hit_type;
  typedef vector<hit_type> data_type;

  /**
   * Thread pool for global fits.
   */
  class JFremantle {
  public:
    /**
     * Constructor.
     *
     * \param  oid            detector identifier
     * \param  katoomba       global fit
     * \param  output         output
     * \param  N              number of threads
     */
    JFremantle(const std::string&   oid,
               const JKatoomba_t&   katoomba,
               JObjectOutput<JEvt>& output,
               const size_t         N) :
      stop(false),
      oid(oid),
      output(output)
    {
      using namespace std;

      for (size_t i = 0; i < N; ++i) {

        thread worker([this, katoomba]() {

          data_type data;

          for (JKatoomba_t fremantle(katoomba); ; ) {
          
            {
              unique_lock<mutex> lock(in);

              cv.wait(lock, [this]() { return stop || !input.empty(); });

              if (stop && input.empty()) {
                return;
              }

              data.swap(input.front());

              input.pop();
            }

            const auto result = fremantle(data.begin(), data.end());

            if (result.chi2 / result.ndf <= fremantle.parameters.chi2perNDF) {

              JRange<double> range(JRange<double>::DEFAULT_RANGE);

              for (data_type::const_iterator i = result.begin; i != result.end; ++i) {
                range.include(i->getValue());
              }

              const JEvt evt = getEvt(JHead(this->oid,
                                            range.getLowerLimit(),
                                            range.getUpperLimit(),
                                            data.size(),
                                            result.value.getN(),
                                            result.ndf,
                                            result.chi2),
                                      result.value);

              {
                unique_lock<mutex> lock(out);

                this->output.put(evt);
              }
            }
          }
        });

        workers.emplace_back(move(worker));
      }
    }


    /**
     * Destructor.
     */   
    ~JFremantle()
    {
      using namespace std;

      {
        unique_lock<mutex> lock(in);

        stop = true;
      }
      
      cv.notify_all();

      for (auto& worker : workers) {
        worker.join();
      }
    }


    /**
     * Get number of pending data.
     *
     * \return                number of pending data
     */
    size_t backlog()
    {
      using namespace std;

      {
        unique_lock<mutex> lock(in);

        return input.size();
      }
    }


    /**
     * Queue data.
     *
     * \param  data           data
     */
    void enqueue(data_type& data)
    {
      using namespace std;
    
      {
        unique_lock<mutex> lock(in);

        if (stop) {
          throw runtime_error("The thread pool has been stopped.");
        }

        input.emplace(move(data));
      }

      cv.notify_one();
    }

  private:
    vector<thread>       workers;
    queue <data_type>    input;
    mutex                in;
    mutex                out;
    condition_variable   cv;
    bool                 stop;
    std::string          oid;
    JObjectOutput<JEvt>& output;
  };
}


/**
 * \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
  bool                    unify;                                // unify weighing of pings
  disable_container       disable;                              // disable tansmissions
  size_t                  jobs;                                 // number of parallel jobs
  int                     sleep_us;                             // sleep time [us]
  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['u'] = make_field(unify,           "unify weighing of pings");
    zap['!'] = make_field(disable,         "disable transmission")           = JPARSER::initialised();
    zap['N'] = make_field(jobs,            "number of parallel jobs")        = 1;
    zap['s'] = make_field(sleep_us,        "sleep time [us]")                = 100;
    zap['d'] = make_field(debug)                                             = 1;

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

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

  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);
  JKatoomba_t katoomba(geometry, V, parameters);

  JKatoomba<JSimplex>::debug = debug;
  JKatoomba<JGandalf>::debug = debug;

  JKatoomba<JSimplex>::MAXIMUM_ITERATIONS = 10000;
  JKatoomba<JGandalf>::MAXIMUM_ITERATIONS = 10000;

  DEBUG(geometry);

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

  try {
    JFremantle fremantle(oid, katoomba, outputFile, jobs);

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

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

          data_type data;

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

            sort(evt->begin(), evt->end(), compare);

            JEvent::iterator __end  = unique(evt->begin(), evt->end(), make_comparator(&JTransmission::getID, JComparison::eq()));

            const JEmitter& emitter = emitters[evt->getID()];
            const double    signal  = (unify ? (double) minimum_number_of_pings / (double) numberOfPings[evt->getID()] : 1.0);

            for (JEvent::const_iterator i = evt->begin(); i != __end; ++i) {

              if (disable.count(JTransmission_t(evt->getID(), i->getID())) == 0 &&
                  disable.count(JTransmission_t(-1,           i->getID())) == 0) {

                if (receivers.has(i->getID()) && geometry.hasLocation(receivers[i->getID()]) && i->getQ() >= parameters.Qmin) {

                  data.push_back(hit_type(emitter,
                                          distance(zbuf.begin(),evt),
                                          receivers[i->getID()],
                                          JPDFGauss(i->getToA(), parameters.sigma_s, signal, parameters.background)));
                }
              }
            }
          }

          while (fremantle.backlog() > jobs) {
            this_thread::sleep_for(chrono::microseconds(sleep_us));
          }

          fremantle.enqueue(data);
        }
      }
    }
    STATUS(endl);
  }
  catch(const exception& error) {
    FATAL("main " << error.what());
  }
  
  JMultipleFileScanner<JMeta> io(inputFile);

  io >> outputFile;

  outputFile.close();
}