MHO_MBDelaySearch.hh

Go to the documentation of this file.

#ifndef MHO_MBDelaySearch_HH__
#define MHO_MBDelaySearch_HH__
 
#include <cmath>
#include <complex>
 
#include "MHO_AxisPack.hh"
#include "MHO_ContainerDefinitions.hh"
#include "MHO_CyclicRotator.hh"
#include "MHO_DelayRate.hh"
#include "MHO_ExtremaSearch.hh"
#include "MHO_InspectingOperator.hh"
#include "MHO_TableContainer.hh"
#include "MHO_UniformGridPointsCalculator.hh"
 
#ifdef HOPS_USE_FFTW3
    #include "MHO_MultidimensionalFastFourierTransformFFTW.hh"
#else
    #include "MHO_MultidimensionalFastFourierTransform.hh"
#endif
 
namespace hops
{
 
using mbd_axis_pack = MHO_AxisPack< time_axis_type >;
using mbd_type = MHO_TableContainer< visibility_element_type, mbd_axis_pack >;
using mbd_amp_type = MHO_TableContainer< double, mbd_axis_pack >;
 
class MHO_MBDelaySearch: public MHO_InspectingOperator< visibility_type >
{
    public:
        MHO_MBDelaySearch();
        virtual ~MHO_MBDelaySearch();
 
        void SetWeights(weight_type* wt_data) { fWeights = wt_data; }
 
        void SetReferenceFrequency(double ref_freq) { fRefFreq = ref_freq; }
 
        //configure the search windows (using floating point limits)
        //default is the full range
        void SetSBDWindow(double low, double high);
        void SetMBDWindow(double low, double high);
        void SetDRWindow(double low, double high);
 
        //retrieve the window limits (that were actually used)
        void GetSBDWindow(double& low, double& high) const;
        void GetMBDWindow(double& low, double& high) const;
        void GetDRWindow(double& low, double& high) const;
 
        int GetMBDMaxBin() const { return fMBDMaxBin; }
 
        int GetSBDMaxBin() const { return fSBDMaxBin; }
 
        int GetDRMaxBin() const { return fDRMaxBin; }
 
        double GetCoarseMBD() const { return fCoarseMBD; }
 
        double GetCoarseSBD() const { return fCoarseSBD; }
 
        double GetCoarseDR() const { return fCoarseDR; }
 
        void SetCoherenceTime(double t_cohere)
        {
            fTCohere = t_cohere;
            fTCohereEnabled = (t_cohere > 0.0);
        }
 
        double GetSBDBinSeparation() const { return fSBDBinSep; }
 
        int GetNMBDBins() const { return fNGridPoints; };
 
        int GetNSBDBins() const { return fNSBD; };
 
        int GetNDRBins() const { return fNDR; };
 
        int GetNDRSPBins() const { return fDelayRateCalc.GetDelayRateSearchSpaceSize(); }
 
        double GetNPointsSearched() const;
 
        double GetSearchMaximumAmplitude() const { return fMax; }
 
        double GetFrequencySpacing() const { return fGridSpace; }
 
        double GetAverageFrequency() const { return fAverageFreq; }
 
        time_axis_type* GetMBDAxis() { return &fMBDAxis; };
 
        delay_rate_axis_type* GetDRAxis() { return &fDRAxis; };
 
    protected:
        using XArgType = visibility_type;
 
        virtual bool InitializeImpl(const XArgType* in) override;
        virtual bool ExecuteImpl(const XArgType* in) override;
 
        bool finalize_search();
 
        void apply_dr_boxcar_smooth(mbd_dr_type& buffer, std::vector< double >& scratch);
 
        // std::vector< double > DetermineFrequencyPoints(const XArgType* in);
 
    protected:
        void SetWindow(double* win, double low, double high);
        void GetWindow(const MHO_Axis< double >& axis, bool win_set, const double* win, double bin_width, double& low,
                       double& high) const;
 
        //workspace
        bool fInitialized;
        std::vector< double > fChannelFreqs;
 
        visibility_type fSBDDrWorkspace;
        visibility_type sbd_dr_data;
        weight_type fWeightsWorkspace;
        weight_type* fWeights;
 
        mbd_type fMBDWorkspace;
        mbd_amp_type fMBDAmpWorkspace;
 
        //dims and parameters
        double fGridStart;
        double fGridSpace;
        double fAverageFreq;
        std::size_t fNGridPoints;
        std::size_t fNSBD;
        std::size_t fNDR;
        std::size_t fNDRSP;
        std::map< std::size_t, std::size_t > fMBDBinMap;
        double fRefFreq;
        int fSBDStart;
        int fSBDStop;
 
        //incoherent averaging parameters
        bool fTCohereEnabled;
        double fTCohere;   //coherence time (s)
        double fAccPeriod; //accumulation period (s), derived from time axis spacing
 
        //scratch space for apply_dr_boxcar_smooth (base + CUDA); OpenMP uses fPerThreadSmoothScratch instead
        std::vector< double > fSmoothedAmpBuffer; //one DR amplitude row, size fNDRSP
 
        //the bin width in each dimension
        double fSBDBinSep;
        double fDRBinSep;
        double fMBDBinSep;
 
        //the window limits
        double fSBDWin[2];
        double fMBDWin[2];
        double fDRWin[2];
 
        //indicates if window limits have been set
        bool fSBDWinSet;
        bool fMBDWinSet;
        bool fDRWinSet;
 
        //location and value of the maximum
        double fMax;
        int fMBDMaxBin;
        int fSBDMaxBin;
        int fDRMaxBin;
 
        double fCoarseMBD;
        double fCoarseSBD;
        double fCoarseDR;
 
        //the number of points searched
        double fNPointsSearched;
 
        MHO_Axis< double > fSBDAxis;
        MHO_Axis< double > fMBDAxis;
        MHO_Axis< double > fDRAxis;
 
        MHO_UniformGridPointsCalculator fGridCalc;
        MHO_DelayRate fDelayRateCalc; //delay rate calculator
 
        //associated info for channel index re-mapping (needed for combining double sideband channel pairs)
        std::map< std::size_t, std::size_t > fChannelIndexToFreqPointIndex;
        //precomputed flat lookup: fMBDBinForChannel[ch] = fMBDBinMap[fChannelIndexToFreqPointIndex[ch]]
        //avoids two std::map traversals per (sbd, dr, ch) iteration in the hot loop
        std::vector< std::size_t > fMBDBinForChannel;
 
#ifdef HOPS_USE_FFTW3
        using FFT_ENGINE_TYPE = MHO_MultidimensionalFastFourierTransformFFTW< mbd_type >;
        using FFT_2D_ENGINE_TYPE = MHO_MultidimensionalFastFourierTransformFFTW< mbd_dr_type >;
#else
        using FFT_ENGINE_TYPE = MHO_MultidimensionalFastFourierTransform< mbd_type >;
        using FFT_2D_ENGINE_TYPE = MHO_MultidimensionalFastFourierTransform< mbd_dr_type >;
#endif
 
        FFT_ENGINE_TYPE fFFTEngine;           //1D engine - used only for MBD axis label setup
        FFT_2D_ENGINE_TYPE fBatchedFFTEngine; //2D engine - batched FFT over all DR bins at once
 
        //2D [DR x MBD] host buffer for batched FFT (see MHO_MBDelaySearchCUDA for the GPU analogue)
        mbd_dr_type fSearchBuffer;
 
        MHO_CyclicRotator< mbd_type > fCyclicRotator;
 
        MHO_ExtremaSearch< mbd_amp_type > fMaxSearch;
};
 
} // namespace hops
 
#endif