SystemCalibrator.h

//-*-C++-*-
/***************************************************************************
 *
 *   Copyright (C) 2008-2025 by Willem van Straten
 *   Licensed under the Academic Free License version 2.1
 *
 ***************************************************************************/
 
// psrchive/More/Polarimetry/Pulsar/SystemCalibrator.h
 
#ifndef __Pulsar_SystemCalibrator_H
#define __Pulsar_SystemCalibrator_H
 
#include "Pulsar/PolnCalibrator.h"
 
#include "Pulsar/ReceptionModelSolver.h"
#include "Pulsar/SignalPath.h"
 
#include "Pulsar/CoherencyMeasurementSet.h"
#include "Pulsar/SourceEstimate.h"
#include "Pulsar/Processor.h"
 
#include "BatchQueue.h"
 
#include <ostream>
 
namespace Pulsar
{
  class ReferenceCalibrator;
  class FluxCalibrator;
  class CalibratorStokes;
  class VariableTransformationManager;
  class VariableFaradayRotation;


 
  class SystemCalibrator : public PolnCalibrator
  {
 
  public:
 
    typedef Calibration::ReceptionModel::Solver Solver;

    class StepFinder;
    
    SystemCalibrator (Archive* archive = 0);

    SystemCalibrator (const SystemCalibrator&);

    virtual ~SystemCalibrator ();


    virtual void share (SystemCalibrator* other);

    Calibrator::Info* get_Info () const;

    void set_name (const std::string& name);

    std::string get_name () const;
    
    MJD get_epoch () const;

    unsigned get_nchan () const;

    unsigned get_ndata (unsigned ichan) const;

    virtual unsigned get_nstate () const;
    
    virtual unsigned get_nstate_pulsar () const;

    virtual unsigned get_state_is_pulsar (unsigned istate) const;

    virtual bool has_cal () const { return calibrator_estimate.size(); }

    virtual Calibrator::Info* new_info_pulsar (unsigned istate) const;

    void set_refcal_through_frontend (bool flag = true);
    
    void set_flux_calibrator (const FluxCalibrator* fluxcal);
    
    void set_calibrators (const std::vector<std::string>& filenames);

    void set_calibrator_preprocessor (Processor*);

    virtual void set_calibrator (const Archive*);

    virtual void set_response( MEAL::Complex2* );

    void set_previous_solution (const PolnCalibrator*);

    void set_previous_cal (const CalibratorStokes*);

    void set_response_fixed (const std::vector<unsigned>&);

    virtual void set_response_variation ( unsigned iparam,
                                          MEAL::Univariate<MEAL::Scalar>* );

    virtual void set_impurity( MEAL::Real4* );

    virtual void set_projection( VariableTransformationManager* );

    virtual void set_faraday_rotation( VariableFaradayRotation* );

    virtual void set_gain( MEAL::Univariate<MEAL::Scalar>* );

    virtual void set_diff_gain( MEAL::Univariate<MEAL::Scalar>* );

    virtual void set_diff_phase( MEAL::Univariate<MEAL::Scalar>* );

    virtual void add_gain_step (const MJD&);

    virtual void add_diff_gain_step (const MJD&);

    virtual void add_diff_phase_step (const MJD&);

    virtual void add_step (const MJD&, Calibration::VariableBackend*);

    void set_step_after_cal (bool val = true) { step_after_cal = val; }
    bool get_step_after_cal () const { return step_after_cal; }

    virtual void set_foreach_calibrator( const MEAL::Complex2* );

    virtual void set_stepeach_calibrator (const Calibration::VariableBackend*);

    virtual void preprocess (Archive* data);

    virtual void add_observation (const Archive* data);

    virtual void set_normalize_by_invariant (bool flag = true);

    virtual void set_normalize_calibrated_by_invariant (bool flag = true);

    virtual void add_pulsar (const Archive* data);

    virtual void add_calibrator (const Archive* data);

    virtual void add_calibrator (const ReferenceCalibrator* polncal);

    virtual bool calibrator_match (const Archive*, std::string& reason_for_not_matching);

    virtual void add_epoch (const MJD& epoch);
   
    MJD get_start_epoch () const { return start_epoch; }

    MJD get_end_epoch () const { return end_epoch; }

    virtual void set_nthread (unsigned nthread);

    virtual void set_equation_configuration (const std::vector<std::string>&);

    virtual bool has_solver () const;

    virtual const Solver* get_solver (unsigned ichan) const;

    virtual void set_solver (Solver*);

    virtual Solver* get_solver ();
    
    virtual void set_retry_reduced_chisq (float);

    virtual void set_invalid_reduced_chisq (float);

    virtual void set_report_projection (bool);

    virtual void set_report_initial_state (bool flag = true);

    virtual void set_report_input_data (bool flag = true);

    virtual void set_report_input_failed (bool flag = true);

    virtual void set_report_data_and_model (bool flag = true);

    virtual void set_report_total_invariant (bool flag = true);

    void set_cal_outlier_threshold (float f) { cal_outlier_threshold = f; }

    float get_cal_outlier_threshold () const { return cal_outlier_threshold; }

    void set_cal_intensity_threshold (float f) { cal_intensity_threshold = f; }

    float get_cal_intensity_threshold () const { return cal_intensity_threshold; }

    void set_cal_polarization_threshold (float f) { cal_polarization_threshold = f; }

    float get_cal_polarization_threshold () const { return cal_polarization_threshold; }

    void set_step_finder (StepFinder*);

    StepFinder* get_step_finder ();

    virtual void solve ();

    virtual void reset ();

    virtual bool get_solved () const;

    virtual bool get_valid (unsigned ichan) const;
    virtual void set_valid (unsigned ichan, bool value, const std::string& reason);

    virtual bool has_valid () const;

    float get_reduced_chisq (unsigned ichan) const;

    virtual void precalibrate (Archive* archive);

    virtual MEAL::Complex2* get_transformation (const Archive* data, unsigned isub, unsigned ichan);

    virtual MEAL::Complex2* get_transformation (const MJD& epoch, Signal::Source source, unsigned ichan);

    virtual Archive* new_solution (const std::string& archive_class) const;

    virtual const CalibratorStokes* get_CalibratorStokes () const;

    virtual const Calibration::SignalPath* get_model (unsigned ichan) const;


    virtual bool match (const Archive*, bool throw_exception = true);
 
    std::string get_mismatch_reason () const { return mismatch_reason; }

    class Unloader;
 
  protected:
 
    friend class SystemCalibratorManager;
    bool match_check_nchan;
    std::string mismatch_reason;
  
    /* and common polarimetric transformations */
    Reference::To<SystemCalibrator> partner;
    void setup_sharing (unsigned ichan);
    
    friend class SystemCalibratorPlotter;
    friend class MatrixTemplateMatching;
 
    std::string instance_name;
    
    virtual void prepare (const Archive* data);

    virtual void calculate_transformation ();

    void create_model ();

    Reference::Vector<Calibration::SignalPath> model;

    std::vector<std::string> equation_configuration;

    Reference::To<Solver> solver;

    Reference::To<const FluxCalibrator> flux_calibrator;
    
    Reference::To<VariableTransformationManager> projection;

    Reference::To<VariableFaradayRotation> faraday_rotation;

    mutable Reference::To<const CalibratorStokes> calibrator_stokes;

    Reference::To< MEAL::Complex2 > response;

    std::vector<unsigned> response_fixed;
    
    Reference::To< MEAL::Real4 > impurity;
 
    typedef MEAL::Univariate<MEAL::Scalar> UniScalar;
    
    std::map< unsigned, Reference::To<UniScalar> > response_variation;

    Reference::To<UniScalar> gain_variation;

    Reference::To<UniScalar> diff_gain_variation;

    Reference::To<UniScalar> diff_phase_variation;
 
    std::vector<MJD> gain_steps;
    std::vector<MJD> diff_gain_steps;
    std::vector<MJD> diff_phase_steps;

    Reference::To< const MEAL::Complex2 > foreach_calibrator;

    Reference::To< const Calibration::VariableBackend > stepeach_calibrator;

    virtual void init_model (unsigned ichan);
 
    typedef Reference::To<Calibration::SourceEstimate> SourceEstimate;
    
    virtual void init_estimates (std::vector<SourceEstimate>&, unsigned ibin = 0);

    virtual void print_input_failed (const std::vector<SourceEstimate>&);
    std::vector<std::ofstream*> input_failed;
    virtual void close_input_failed ();

    virtual void prepare_calibrator_estimate (Signal::Source);

    virtual void create_calibrator_estimate ();

    virtual void copy_calibrator_estimate ();

    bool data_submitted;

    std::vector< std::vector<Calibration::SourceObservation> > calibrator_data;
    
    // submit all calibrator data
    virtual void submit_calibrator_data ();
 
    virtual void submit_calibrator_data (Calibration::CoherencyMeasurementSet&,
                                         const Calibration::SourceObservation&);
 
    virtual void integrate_calibrator_data (const Calibration::SourceObservation&);
 
    virtual void integrate_calibrator_solution (const Calibration::SourceObservation&);

    virtual void load_calibrators ();
    
    Reference::To<StepFinder> step_finder;
    
    std::vector< std::vector<Calibration::CoherencyMeasurementSet> > pulsar_data;

    virtual void add_pulsar (const Archive* data, unsigned isub);


    virtual void add_pulsar (Calibration::CoherencyMeasurementSet&, const Integration*, unsigned ichan) = 0;

    virtual void integrate_pulsar_data (const Calibration::CoherencyMeasurementSet&) { }

    virtual void submit_pulsar_data ();

    virtual void submit_pulsar_data (Calibration::CoherencyMeasurementSet&);

    virtual double get_invariant (Integration* subint, unsigned ichan) = 0;

    std::vector<std::string> calibrator_filenames;

    Reference::To<Processor> calibrator_preprocessor;

    std::vector< Reference::To<Calibration::SourceEstimate> > calibrator_estimate;
    
    MJD start_epoch;

    MJD end_epoch;
    
    bool step_after_cal;

    bool normalize_by_invariant;

    bool normalize_calibrated_by_invariant = false;

    bool solve_in_reverse_channel_order;

    bool correct_interstellar_Faraday_rotation;
    
    bool refcal_through_frontend;

    bool guess_physical_calibrator_stokes;
    
    float retry_chisq;

    float invalid_chisq;

    bool report_projection;

    bool report_initial_state;

    bool report_input_data;

    bool report_input_failed;

    bool report_data_and_model;

    bool report_total_invariant = false;

    double cal_outlier_threshold;

    double cal_intensity_threshold;

    double cal_polarization_threshold;

    virtual void solve_prepare ();

    void configure ( MEAL::Function* equation );

    virtual void export_prepare () const = 0;

    BatchQueue queue;

    bool get_prepared () const;

    Solver* get_solver (unsigned ichan);

    void resolve (unsigned ichan);

    void check_ichan (const char* name, unsigned ichan) const;
 
    unsigned get_data_fail;
    unsigned get_data_call;

    Reference::To<const PolnCalibrator> previous;

    Reference::To<const CalibratorStokes> previous_cal;

    virtual void load_previous ();

    bool previous_loaded = false;

    bool has_pulsar = false;

    Jones<double> invert_basis;

    std::ofstream invint_out;
  
  private:
 
    std::vector<bool> epoch_added;

    bool is_solved;

    bool is_prepared;
 
  };
 
}
 
#endif