SpectralKurtosis.h

//-*-C++-*-
/***************************************************************************
 *
 *   Copyright (C) 2016 by Andrew Jameson and Willem van Straten
 *   Licensed under the Academic Free License version 2.1
 *
 ***************************************************************************/
 
#include "dsp/Transformation.h"
#include "dsp/WeightedTimeSeries.h"
#include "dsp/BitSeries.h"
#include "dsp/Memory.h"
#include "EventEmitter.h"
 
#ifndef __SpectralKurtosis_h
#define __SpectralKurtosis_h
 
#define ZAP_ALL  0
#define ZAP_SKFB 1
#define ZAP_FSCR 2
#define ZAP_TSCR 3
 
namespace dsp {
 
  class SKLimits;
 
 
  class SpectralKurtosis: public Transformation<TimeSeries,TimeSeries> {
 
  public:
 
    SpectralKurtosis ();
 
    ~SpectralKurtosis ();
 
    bool get_order_supported (TimeSeries::Order order) const;
 
    void load_configuration (const std::string& filename);
 
    void set_M (unsigned _M) { resolution[0].set_M(_M); }
    void set_M (const std::vector<unsigned>&);
 
    void set_noverlap (unsigned _nover) { resolution[0].set_noverlap(_nover); }
    void set_noverlap (const std::vector<unsigned>&);
 
    void set_thresholds (float _std_devs);
    void set_thresholds (const std::vector<float>&);
 
    void set_detect_time_freq (bool flag=true) { detect_time_freq = flag; }
 
    void set_detect_freq (bool flag=true) { detect_freq = flag; }
 
    void set_detect_time (bool flag=true) { detect_time = flag; }
 
 
    void set_omit_outliers (bool flag=true) { omit_outliers = flag; }
 
    unsigned get_nchan () const { return nchan; }
 
    unsigned get_npol () const { return sums_npol; }
 
    void set_channel_range (unsigned start, unsigned end);
 
    void reserve ();
 
    void prepare ();
 
    void prepare_output ();
 
    double get_delay_time () const override;
    
    unsigned get_M () const
    { return resolution[0].get_M(); }
 
    unsigned get_excision_threshold () const
    { return resolution[0].get_std_devs(); }
 
    void get_filtered_sum (std::vector<float>& sum) const
    { sum = filtered_sum; }
 
    void get_filtered_hits (std::vector<uint64_t>& hits) const
    { hits = filtered_hits; }
 
    void get_unfiltered_sum (std::vector<float>& sum) const
    { sum = unfiltered_sum; }
 
    uint64_t get_unfiltered_hits () const { return unfiltered_hits; }
 
    void reset_count () { unfiltered_hits = 0; }
 
    class Engine;
 
    void set_engine (Engine*);
 
    template<class T>
    class Reporter {
    public:
      virtual void operator() (T*, unsigned, unsigned, unsigned, unsigned) {};
    };
 
    // An event emitter that takes a data array, and the nchan, npol, ndat and ndim
    // associated with the data array
    EventEmitter<Reporter<float> > float_reporter;
 
    // This is for reporting the state of the bit zapmask
    EventEmitter<Reporter<unsigned char> > char_reporter;
 
    bool get_report () const { return report; }
 
    void set_report (bool _report) { report = _report; }
 
    bool has_zero_DM_input () const;
    virtual void set_zero_DM_input (TimeSeries* zero_DM_input);
    virtual const TimeSeries* get_zero_DM_input() const;
    virtual TimeSeries* get_zero_DM_input();
 
    // bool has_zero_DM_input_container () const;
    // virtual void set_zero_DM_input_container (const HasInput<TimeSeries> zero_DM_input_container&);
    // virtual const HasInput<TimeSeries>& get_zero_DM_input_container() const;
    // virtual HasInput<TimeSeries>& get_zero_DM_input_container();
 
    virtual void set_zero_DM_buffering_policy (BufferingPolicy* policy)
    { zero_DM_buffering_policy = policy; }
 
    bool has_zero_DM_buffering_policy() const
    { return zero_DM_buffering_policy; }
 
    BufferingPolicy* get_zero_DM_buffering_policy () const
    { return zero_DM_buffering_policy; }
 
    bool get_zero_DM () const { return zero_DM; }
 
    void set_zero_DM (bool _zero_DM) { zero_DM = _zero_DM; }
 
  protected:
 
    void transformation ();
 
    Reference::To<Engine> engine;
 
  private:
 
    void compute ();
 
    void detect ();
    void detect_tscr ();
    void detect_skfb (unsigned ires);
    void detect_fscr (unsigned ires);
    void count_zapped ();
 
    void mask ();
    void reset_mask ();
 
    void insertsk ();
 
    unsigned debugd = 1;
 
    class Count
    {
      public:
        uint64_t tested = 0;
        uint64_t zapped = 0;
        double fraction_zapped () const { return (tested) ? double(zapped)/tested : 0.0; }
        const Count& operator+= (const Count& that) { tested+=that.tested; zapped+=that.zapped; return *this;}
    };
 
    class Resolution
    {
    private:
      
      mutable std::vector<bool> channels;
 
      mutable std::vector<float> thresholds;
 
      unsigned M = 128;
 
      unsigned Nd = 1;
 
      float std_devs = 3.0;
 
      unsigned noverlap = 1;
 
      unsigned overlap_offset = 128;
 
      uint64_t npart = 0;
 
      uint64_t output_ndat = 0;
 
      Count SK_time_freq;
 
      Count SK_time;
 
      void set_thresholds (bool verbose = false) const;
 
    public:
 
 
      void add_include (const std::pair<unsigned, unsigned>&);
 
 
      void add_exclude (const std::pair<unsigned, unsigned>&);
 
      const std::vector<bool>& get_channels (unsigned nchan) const;
      
      unsigned get_M () const { return M; }
      void set_M (unsigned);
 
      unsigned get_Nd () const { return Nd; }
      void set_Nd (unsigned);
 
      float get_std_devs () const { return std_devs; }
      void set_std_devs (float);
 
      /* blocks used to estimate SK overlap by (noverlap-1)/noverlap */
      unsigned get_noverlap () const { return noverlap; }
      void set_noverlap (unsigned);
 
      void prepare (uint64_t ndat = 0);
 
      unsigned get_overlap_offset () const { return overlap_offset; }
 
      uint64_t get_npart () const { return npart; }
 
      uint64_t get_output_ndat () const { return output_ndat; }
 
      void compatible (Resolution& that);
 
      const std::vector<float>& get_thresholds () const;
      
      std::vector< std::pair<unsigned,unsigned> > include;
 
      std::vector< std::pair<unsigned,unsigned> > exclude;
 
      void increment_time_freq (const Count&);
      const Count& get_count_time_freq () const;
 
      void increment_time (const Count&);
      const Count& get_count_time () const;
    };
 
    std::vector<Resolution> resolution;
 
    void resize_resolution (unsigned);
 
    void tscrunch_sums (Resolution& from, Resolution& to);
 
    // for sorting by M
    static bool by_M (const Resolution& A, const Resolution& B);
 
    unsigned nchan = 0;
    unsigned npol = 0;
    unsigned ndim = 0;
    unsigned integrated_Nd = 1;
 
    unsigned sums_npol = 0;
 
    Reference::To<WeightedTimeSeries> sums;
 
    Reference::To<BitSeries> zapmask;
 
    std::vector<SKLimits*> dynamic_limits;
    const SKLimits* get_limits (unsigned count, float std_devs);
 
    std::vector<float> filtered_sum;
 
    std::vector<uint64_t> filtered_hits;
 
    std::vector<float> unfiltered_sum;
 
    uint64_t unfiltered_hits = 0;
 
    Count SK_freq;
 
    bool detect_time_freq = true;
    bool detect_freq = true; // after tscrunch
    bool detect_time = true; // after fscrunch
 
 
    bool omit_outliers = true;
 
    bool prepared = false;
 
    bool report = false;
 
    // //! Input TimeSeries that has not been dedispersed in some previous operation.
    // Reference::To<dsp::TimeSeries> zero_DM_input;
 
    HasInput<TimeSeries> zero_DM_input_container;
 
    Reference::To<BufferingPolicy> zero_DM_buffering_policy;
 
    bool zero_DM = false;
    double delay_time = 0.0;
  };
 
  class SpectralKurtosis::Engine : public Reference::Able
  {
  public:
 
      virtual void setup () = 0;
 
      virtual void compute (const TimeSeries* input, TimeSeries* output,
                            TimeSeries *output_tscr, unsigned tscrunch) = 0;
 
      virtual void reset_mask (BitSeries* output) = 0;
 
      virtual void detect_ft (const TimeSeries* input, BitSeries* output,
                              float upper_thresh, float lower_thresh) = 0;
 
      virtual void detect_fscr (const TimeSeries* input, BitSeries* output,
                                const float mu2, const float std_devs,
                                unsigned schan, unsigned echan) = 0;
 
      virtual void detect_tscr (const TimeSeries* input,
                                const TimeSeries * input_tscr,
                                BitSeries* output,
                                float upper, float lower) = 0;
 
      virtual int count_mask (const BitSeries* output) = 0;
 
      virtual float * get_estimates (const TimeSeries* input) = 0;
 
      virtual unsigned char * get_zapmask (const BitSeries* input) = 0;
 
      virtual void mask (BitSeries* mask, const TimeSeries * in, TimeSeries* out, unsigned M) = 0;
 
      virtual void insertsk (const TimeSeries* input, TimeSeries* out, unsigned M) = 0;
 
  };
}
 
#endif