Subint.h

//-*-C++-*-
/***************************************************************************
 *
 *   Copyright (C) 2010 by Paul Demorest
 *   Licensed under the Academic Free License version 2.1
 *
 ***************************************************************************/
 
#ifndef __Subint_h
#define __Subint_h
 
#include "dsp/TimeDivide.h"
#include "dsp/PhaseSeries.h"
 
#include "Phase.h"
#include "Callback.h"
 
#define VERBOSE true
// #define VERBOSE Op::verbose
 
namespace dsp {
 
  class PhaseSeriesUnloader;
 
 
  template <class Op>
  class Subint : public Op {
 
  public:
    
    Subint ();
    
    ~Subint ();
    
    Op* clone () const;
 
    void finish ();
 
 
    Callback<PhaseSeries*> complete;
 
    Callback<PhaseSeries*> partial;
 
 
    void set_start_time (const MJD& start_time)
      { divider.set_start_time (start_time); }
 
    MJD get_start_time () const { return divider.get_start_time(); }
 
    void set_subint_seconds (double subint_seconds)
      { divider.set_seconds (subint_seconds); }
 
    double get_subint_seconds () const { return divider.get_seconds(); }
 
    void set_subint_reference_epoch (const MJD& epoch)
    { divider.set_reference_epoch (epoch); }
 
    MJD get_subint_reference_epoch () const
    { return divider.get_reference_epoch (); }
 
    void set_subint_turns (unsigned subint_turns)
      { divider.set_turns (subint_turns); }
 
    unsigned get_subint_turns () const { return unsigned(divider.get_turns());}
 
    void set_fractional_pulses (bool flag) 
      { divider.set_fractional_pulses (flag); }
 
    void set_reference_phase (double phase)
      { divider.set_reference_phase (phase); Op::set_reference_phase(phase); }
 
 
    void set_unloader (PhaseSeriesUnloader* unloader);
 
    PhaseSeriesUnloader* get_unloader () const { return unloader; }
 
    virtual bool keep (PhaseSeries* data) { return true; }
 
 
    const TimeDivide* get_divider () const { return &divider; }
 
    void set_cerr (std::ostream& os) const;
 
  protected:
 
    virtual void unload_partial ();
 
    virtual void transformation ();
 
    virtual void set_limits (const Observation* input);
 
    Reference::To<PhaseSeriesUnloader> unloader;
    
    TimeDivide divider;
 
    void prepare ();
 
    void zero_output ();
 
    bool built;
 
  };
 
}
 
#include "dsp/SignalPath.h"
#include "dsp/PhaseSeriesUnloader.h"
 
#include "Error.h"
 
//using namespace std;
 
template <class Op>
dsp::Subint<Op>::Subint ()
{
  built = false;
  std::string orig_name = Op::get_name();
  Op::set_name("Subint<" + orig_name + ">");
}
 
template <class Op>
dsp::Subint<Op>::~Subint ()
{
  if (Op::verbose)
    this->cerr << "dsp::Subint::~Subint" << std::endl;
}
 
template <class Op>
Op* dsp::Subint<Op>::clone () const
{
  if (Op::verbose)
    this->cerr << "dsp::Subint<Op>::clone" << std::endl;
  return new Subint<Op>(*this);
}
 
template <class Op>
void dsp::Subint<Op>::set_cerr (std::ostream& os) const
{
  Op::set_cerr (os);
  divider.set_cerr (os);
  if (unloader)
    unloader->set_cerr (os);
}
 
template <class Op>
void dsp::Subint<Op>::set_unloader (dsp::PhaseSeriesUnloader* _unloader)
{
  if (Op::verbose)
    this->cerr << "dsp::Subint::set_unloader this=" << this 
         << " unloader=" << _unloader << std::endl;
 
  unloader = _unloader;
  if (unloader)
    unloader->set_cerr (this->cerr);
}
 
template <class Op>
void dsp::Subint<Op>::prepare () try
{
  if (Op::verbose)
    this->cerr << "dsp::Subint::prepare call Op::prepare" << std::endl;
 
  Op::prepare ();
 
  // if unspecified, the first TimeSeries to be folded will define the
  // start time from which to begin cutting up the observation
  if (divider.get_start_time() == MJD::zero)
  {
    if (Op::verbose)
      this->cerr << "dsp::Subint::prepare set divider start time=" 
           << Op::input->get_start_time() << std::endl;
 
    if (Op::input->get_start_time() == MJD::zero)
      throw Error (InvalidState, "dsp::Subint::prepare",
                   "input start time not set");
 
    divider.set_start_time (Op::input->get_start_time());
  }
 
  if (divider.get_turns())
  {
    if (Op::has_folding_predictor())
    {
      divider.set_predictor (Op::get_folding_predictor());
    }
    else
    {
      if (Op::verbose)
        this->cerr << "dsp::Subint::prepare folding period=" << Op::get_folding_period() << std::endl;
      divider.set_period (Op::get_folding_period());
    }
  }
 
  built = true;
}
catch (Error& error)
{
  throw error += "dsp::Subint::prepare";
}
 
 
template <class Op>
void dsp::Subint<Op>::transformation () try
{
  if (Op::verbose)
    this->cerr << "dsp::Subint::transformation" << std::endl;
 
  if (divider.get_turns() == 0 && divider.get_seconds() == 0.0)
    throw Error (InvalidState, "dsp::Subint::tranformation",
                 "sub-integration length not specified");
 
  if (!built)
    prepare ();
 
  // flag that the input TimeSeries contains data for another sub-integration
  bool more_data = true;
  bool first_division = true;
 
  while (more_data)
  {
    divider.set_bounds( Op::get_input() );
 
    if (!divider.get_fractional_pulses())
      Op::get_output()->set_ndat_expected( divider.get_division_ndat() );
 
    more_data = divider.get_in_next ();
 
    if (first_division && divider.get_new_division())
    {
      /* A new division has been started and there is still data in
         the current integration.  This is a sign that the current
         input comes from uncontiguous data, which can arise when
         processing in parallel. */
 
      if (Op::verbose)
        this->cerr << "dsp::Subint::transformation unload partial (first_division and new division)" << std::endl;
 
      unload_partial ();
    }
 
    if (!divider.get_is_valid())
      continue;
 
    Op::transformation ();
 
    if (!divider.get_end_reached())
      continue;
 
    if (first_division)
    {
      /* When the end of the first division is reached, it is not 100%
         certain that a complete sub-integration is available */
 
      if (Op::verbose)
        this->cerr << "dsp::Subint::transformation unload partial (first_division of transformation)" << std::endl;
 
      unload_partial ();
      first_division = false;
      continue;
    }
 
    if (Op::verbose)
      this->cerr << "dsp::Subint::transformation sub-integration completed" << std::endl;
 
    PhaseSeries* result = Op::get_result ();
 
    complete.send (result);
 
    if (unloader && keep(result))
    {
      if (Op::verbose)
        this->cerr << "dsp::Subint::transformation this=" << this << " unloader=" << unloader.get() << std::endl;
 
      unloader->unload (result);
      zero_output ();
    }
  }
}
catch (Error& error)
{
  throw error += "dsp::Subint::transformation";
}
 
 
template <class Op>
void dsp::Subint<Op>::set_limits (const Observation* input)
{
  Op::idat_start = divider.get_idat_start ();
  Op::ndat_fold = divider.get_ndat ();
}
 
template <class Op>
void dsp::Subint<Op>::finish () try
{
  Op::finish ();
 
  if (Op::get_result()->get_integration_length() != 0)
  {
    if (Op::verbose)
      this->cerr << "dsp::Subint::finish unload_partial" << std::endl;
 
    unload_partial ();
  }
 
  if (unloader)
  {
    if (Op::verbose)
      this->cerr << "dsp::Subint::finish call unloader finish" << std::endl;
 
    unloader->finish ();
  }
}
catch (Error& error)
{
  throw error += "dsp::Subint::finish";
}
 
template <class Op>
void dsp::Subint<Op>::unload_partial () try
{
  PhaseSeries* result = Op::get_result ();
 
  if (Op::verbose)
    this->cerr << "dsp::Subint::unload_partial MJD=" << result->get_start_time().printall() << std::endl;
 
  if (Op::verbose)
    this->cerr << "dsp::Subint::unload_partial this=" << (void*) this
               << " result=" << (void*) result << std::endl;
 
  partial.send (result);
 
  if (unloader)
  {
    if (Op::verbose)
      this->cerr << "dsp::Subint::unload_partial this=" << this
           << " unloader=" << unloader.get() << std::endl;
 
    unloader->partial (result);
  }
 
  zero_output ();
}
catch (Error& error)
{
  throw error += "dsp::Subint::unload_partial";
}
 
#define SIGNAL_PATH
 
template <class Op>
void dsp::Subint<Op>::zero_output ()
{
#ifdef SIGNAL_PATH
  SignalPath* path = 0;
  PhaseSeries* out_ptr = Op::output;
 
  if (out_ptr->has_extensions())
  {
    Extensions* ext = out_ptr->get_extensions();
    path = ext->get<SignalPath>();
  }
 
  if (path)
  {
    if (Op::verbose)
      this->cerr << "dsp::Subint::zero_output calling SignalPath:reset" 
                 << std::endl;
    path->reset();
  }
  else
#endif
    Op::reset();
}
 
#endif