ProfileShiftFit.h

//-*-C++-*-
/***************************************************************************
 *
 *   Copyright (C) 2008 by Paul Demorest
 *   Licensed under the Academic Free License version 2.1
 *
 ***************************************************************************/
 
#ifndef __Pulsar_ProfileShiftFit_h
#define __Pulsar_ProfileShiftFit_h
 
#include "Pulsar/Algorithm.h"
#include "FTransform.h"
#include "Estimate.h"
#include "toa.h"
#include "BoxMuller.h"
 
namespace Pulsar
{
  class Profile;
  class Integration;


  class ProfileShiftFit : public Algorithm
  {   
 
  public:

    enum Error_Method {
      Traditional_Chi2,     // Chi-squared 2nd derivatives
      MCMC_Variance,        // Shift variance from MCMC
      Numerical,            // Shift variance from numerical integration
    };

    ProfileShiftFit ();

    ~ProfileShiftFit();

    void reset();

    void set_nharm(unsigned nh);

    unsigned get_nharm();

    void set_error_method(Error_Method e) { err_meth = e; }

    void set_mcmc_iterations(int nit) { mcmc_it = nit; }

    int get_mcmc_iterations() const { return mcmc_it; }

    void set_standard (const Profile* p);

    void set_Profile (const Profile* p);

    void compute();

    Tempo::toa toa (const Integration*);

    Estimate<double> get_shift();

    Estimate<double> get_scale();

    double get_mse();

    double get_sigma2();

    double get_snr();

    double get_reduced_chisq () const;


    double get_effective_duty_cycle () const;

    void apply_scale_and_shift(Profile *p);

    int mcmc_trials;
    int mcmc_accept;

    bool choose_maximum_harmonic;
 
  protected:

    void init();

    void choose_nharm ();

    unsigned nharm;

    unsigned effective_nharm;

    Reference::To<const Profile> std;

    unsigned nbins_std;

    float *fstd;

    double std_pow;

    Reference::To<const Profile> prof;

    unsigned nbins_prof;

    float *fprof;

    float *fccf;

    unsigned nbins_ccf;
 
    
    double ccf(double phi);

    double dccf(double phi);

    double d2ccf(double phi);

    double log_shift_pdf(double phi);

    double log_shift_pdf_pos(double phi);

    Error_Method err_meth;

    void error_traditional();

    void error_numerical();

    void error_mcmc_pdf_var();

    int mcmc_it;

    bool computed;

    double shift;

    double eshift;

    double correction;

    double scale;

    double escale;

    double sigma2;

    double mse;

    double chi2;

    unsigned dof;

    double snr;
 
  private:

    BoxMuller nrand;

    double mcmc_state;

    double mcmc_log_pdf;

    void mcmc_init();

    double mcmc_sample();

    double max_log_pdf;

    static double f_pdf(double phi, void *_psf);
    static double f_pdf_x2(double phi, void *_psf);
 
  };
 
}
 
 
#endif