JPlatypus.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 <limits>

#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 "JLang/JFileStream.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 "JTools/JQuantile.hh"

#include "JSupport/JMultipleFileScanner.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/JSuperEvt.hh"
#include "JAcoustics/JSuperEvtToolkit.hh"
#include "JAcoustics/JSupport.hh"

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


namespace {

  using namespace std;
  using namespace JPP;
  
  typedef JKatoomba<JGandalf>            JKatoomba_t;
  //typedef JKatoomba<JAbstractMinimiser>  JKatoomba_t;
  
  /**
   * Thread pool for global fits.
   */
  class JPlatypus {
  public:
    /**
     * Constructor.
     *
     * \param  katoomba       global fit
     * \param  ns             number of threads
     */
    JPlatypus(const JKatoomba_t& katoomba, const size_t ns) :
      stop(false)
    {
      using namespace std;

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

        thread worker([this, katoomba]() {

          vector<JHit> data;

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

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

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

              const JSuperEvt& evt = input.front();

              data.clear();
              
              for (JSuperEvt::rx_type::const_iterator hit = evt.rx.begin(); hit != evt.rx.end(); ++hit) {
                
                data.push_back(JHit(emitters[hit->id],
                                    hit->counter,
                                    JLocation(hit->string, hit->floor),
                                    hit->toa,
                                    parameters.sigma_s,
                                    hit->weight));
              }
            
              platypus.value = getModel(evt);

              input.pop();
            }
            
            const double chi2 = platypus (data.begin(), data.end()) / platypus.estimator->getRho(1.0);
            const double ndf  = getWeight(data.begin(), data.end()) - platypus.value.getN();

            {
              unique_lock<mutex> lock(out);

              Q.put(chi2 / ndf);
            }
          }
        });

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


    /**
     * Destructor.
     */   
    ~JPlatypus()
    {
      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 super event.
     *
     * \param  evt            super event
     */
    void enqueue(const JSuperEvt& evt)
    {
      using namespace std;
    
      {
        unique_lock<mutex> lock(in);

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

        input.emplace(evt);
      }

      cv.notify_one();
    }

    static std::string               oid;             //!< detector identifier
    static JQuantile                 Q;               //!< chi2/NDF
    static JFitParameters            parameters;      //!< fit parameters
    static::JHashMap<int, JEmitter>  emitters;        //!< emitters

    
  private:
    vector<thread>       workers;
    queue <JSuperEvt>    input;
    mutex                in;
    mutex                out;
    condition_variable   cv;
    bool                 stop;
  };

  std::string               JPlatypus::oid     =  "";
  JQuantile                 JPlatypus::Q;
  JFitParameters            JPlatypus::parameters;
  JHashMap<int, JEmitter>   JPlatypus::emitters;
}


/**
 * \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;

  JMultipleFileScanner<JSuperEvt>  inputFile;
  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
  size_t                  jobs;                                 // number of parallel jobs
  int                     debug;

  try {

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

    zap['f'] = make_field(inputFile,       "output of JAcousticEventBuilder[.sh]");
    zap['n'] = make_field(numberOfEvents)                                    = JLimit::max();
    zap['a'] = make_field(detectorFile);
    zap['@'] = make_field(JPlatypus::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['N'] = make_field(jobs,            "number of parallel jobs")        = 1;
    zap['d'] = make_field(debug)                                             = 1;

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

  ROOT::EnableThreadSafety();
  
  const int sleep_us = 100;                                           // sleep time [us]

  JDetector detector;

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

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

  for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) {
    try {
      JPlatypus::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, JPlatypus::parameters.option);

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

  katoomba.estimator.reset(getMEstimator(JPlatypus::parameters.mestimator));

  try {
    
    JPlatypus platypus(katoomba, jobs);

    while (inputFile.hasNext()) {

      STATUS("event: " << setw(8) << inputFile.getCounter() << '\r'); DEBUG(endl);

      const JSuperEvt* evt = inputFile.next();

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

      platypus.enqueue(*evt);
    }
    STATUS(endl);
  }
  catch(const exception& error) {
    FATAL("main " << error.what());
  }

  JFileOutputStream(3) << SCIENTIFIC(1,10) << JPlatypus::Q.getMean(numeric_limits<float>::max()) << endl;
}