mcmc_gw

Header file for the Graviational Wave specific MCMC routines

Functions

double maximized_coal_log_likelihood_IMRPhenomD(double *frequencies, int length, std::complex<double> *data, double *noise, double SNR, double chirpmass, double symmetric_mass_ratio, double spin1, double spin2, bool NSflag, fftw_outline *plan)

Function to calculate the log Likelihood as defined by -1/2 (d-h|d-h) maximized over the extrinsic parameters phic and tc.

frequency array must be uniform spacing - this shouldn’t be a problem when working with real data as DFT return uniform spacing

Parameters:: chirpmass – in solar masses

double maximized_coal_log_likelihood_IMRPhenomD(double *frequencies, size_t length, double *real_data, double *imag_data, double *noise, double SNR, double chirpmass, double symmetric_mass_ratio, double spin1, double spin2, bool NSflag)

Parameters:: chirpmass – in solar masses

double maximized_coal_log_likelihood_IMRPhenomD(double *frequencies, size_t length, double *real_data, double *imag_data, double *noise, double SNR, double chirpmass, double symmetric_mass_ratio, double spin1, double spin2, bool NSflag, fftw_outline *plan)

Parameters:: chirpmass – in solar masses

double maximized_coal_log_likelihood_IMRPhenomD_Full_Param(double *frequencies, int length, std::complex<double> *data, double *noise, double chirpmass, double symmetric_mass_ratio, double spin1, double spin2, double Luminosity_Distance, double theta, double phi, double iota, bool NSflag, fftw_outline *plan)

Parameters:: chirpmass – in solar masses

double maximized_coal_log_likelihood_IMRPhenomD_Full_Param(double *frequencies, size_t length, double *real_data, double *imag_data, double *noise, double chirpmass, double symmetric_mass_ratio, double spin1, double spin2, double Luminosity_Distance, double theta, double phi, double iota, bool NSflag)

Parameters:: chirpmass – in solar masses

double maximized_coal_log_likelihood_IMRPhenomD_Full_Param(double *frequencies, size_t length, double *real_data, double *imag_data, double *noise, double chirpmass, double symmetric_mass_ratio, double spin1, double spin2, double Luminosity_Distance, double theta, double phi, double iota, bool NSflag, fftw_outline *plan)

Parameters:: chirpmass – in solar masses

double maximized_Log_Likelihood_aligned_spin_internal(std::complex<double> *data, double *psd, double *frequencies, std::complex<double> *detector_response, size_t length, fftw_outline *plan)

Maximized match over coalescence variables - returns log likelihood NOT NORMALIZED for aligned spins.

Note: this function is not properly normalized for an absolute comparison. This is made for MCMC sampling, so to minimize time, constant terms like (Data|Data), which would cancel in the Metropolis-Hasting ratio, are left out for efficiency

double Log_Likelihood(std::complex<double> *data, double *psd, double *frequencies, size_t length, gen_params_base<double> *params, std::string detector, std::string generation_method): Unmarginalized log of the likelihood.

double maximized_Log_Likelihood_unaligned_spin_internal(std::complex<double> *data, double *psd, double *frequencies, std::complex<double> *hplus, std::complex<double> *hcross, size_t length, fftw_outline *plan)

log likelihood function that maximizes over extrinsic parameters tc, phic, D, and phiRef, the reference frequency - for unaligned spins

Ref: arXiv 1603.02444v2

NOTE: Only works for +/x polarizations

double maximized_Log_Likelihood(std::complex<double> *data, double *psd, double *frequencies, size_t length, gen_params_base<double> *params, std::string detector, std::string generation_method, fftw_outline *plan)

routine to maximize over all extrinsic quantities and return the log likelihood

IMRPhenomD — maximizes over DL, phic, tc, \iota, \phi, \theta IMRPhenomP — maximizes over DL, phic,tc, \psi, \phi , \theta

double maximized_Log_Likelihood(double *data_real, double *data_imag, double *psd, double *frequencies, size_t length, gen_params_base<double> *params, std::string detector, std::string generation_method, fftw_outline *plan)

double maximized_coal_Log_Likelihood(double *data_real, double *data_imag, double *psd, double *frequencies, size_t length, double *template_real, double *template_imag, fftw_outline *plan)

double maximized_coal_Log_Likelihood(std::complex<double> *data, double *psd, double *frequencies, size_t length, gen_params_base<double> *params, std::string detector, std::string generation_method, fftw_outline *plan, double *tc, double *phic): Function to maximize only over coalescence variables tc and phic, returns the maximum values used.

double maximized_coal_Log_Likelihood_internal(std::complex<double> *data, double *psd, double *frequencies, std::complex<double> *detector_response, size_t length, fftw_outline *plan, double *tc, double *phic)

double Log_Likelihood_internal(std::complex<double> *data, double *psd, double *frequencies, double *weights, std::complex<double> *detector_response, int length, bool log10F, std::string integration_method)

Internal function for the unmarginalized log of the likelihood.

.5 * ( ( h | h ) - 2 ( D | h ) )

double Log_Likelihood_internal(std::complex<double> *data, double *psd, std::complex<double> *detector_response, Quadrature *QuadMethod): Compute unmarginalized likelihood with a specified Quadrature method.

double MCMC_likelihood_wrapper_SKYSEARCH(double *param, mcmc_data_interface *interface, void *parameters)

void SkySearch_PTMCMC_MH_dynamic_PT_alloc_uncorrelated_GW(mcmc_sampler_output *sampler_output, double **output, int dimension, int N_steps, int chain_N, int max_chain_N_thermo_ensemble, double *initial_pos, double *seeding_var, double **ensemble_initial_pos, double *chain_temps, int swp_freq, int t0, int nu, int max_chunk_size, std::string chain_distribution_scheme, double (*log_prior)(double *param, mcmc_data_interface *interface, void *parameters), int numThreads, bool pool, bool show_prog, int num_detectors, std::complex<double> **data, double **noise_psd, double **frequencies, int *data_length, double gps_time, std::string *detectors, int Nmod, double *bppe, std::complex<double> *hplus, std::complex<double> *hcross, std::string statistics_filename, std::string chain_filename, std::string likelihood_log_filename, std::string checkpoint_filename)

Takes in intrinsic parameters and conducts a sky search using uncorrelated MCMC.

Searches over RA, sin(DEC), PSI, INCLINATION, phiREF, tc,ln(DL)

Only works with PhenomD

hplus and hcross should be constructed with the value for inclination in initial_pos — its scaled out internally

void PTMCMC_MH_GW(double ***output, int dimension, int N_steps, int chain_N, double *initial_pos, double *seeding_var, double **ensemble_initial_pos, double *chain_temps, int swp_freq, double (*log_prior)(double *param, mcmc_data_interface *interface, void *parameters), int numThreads, bool pool, bool show_prog, int num_detectors, std::complex<double> **data, double **noise_psd, double **frequencies, int *data_length, double gps_time, std::string *detector, MCMC_modification_struct *mod_struct, std::string generation_method, std::string statistics_filename, std::string chain_filename, std::string auto_corr_filename, std::string likelihood_log_filename, std::string checkpoint_filename)

Wrapper for the MCMC_MH function, specifically for GW analysis.

Handles the details of setting up the MCMC sampler and wraps the fisher and log likelihood to conform to the format of the sampler

NOTE — This sampler is NOT thread safe. There is global memory declared for each call to MCMC_MH_GW, so separate samplers should not be run in the same process space

Supported parameter combinations:

IMRPhenomD - 4 dimensions — ln chirpmass, eta, chi1, chi2

IMRPhenomD - 7 dimensions — ln D_L, tc, phic, ln chirpmass, eta, chi1, chi2

IMRPhenomD - 9 dimensions — cos inclination, RA, DEC, ln D_L, ln chirpmass, eta, chi1, chi2, psi

dCS_IMRPhenomD- 8 dimensions — cos inclination, RA, DEC, ln D_L, ln chirpmass, eta, chi1, chi2, \alpha^2 (the coupling parameter)

dCS_IMRPhenomD_root_alpha- 8 dimensions — cos inclination, RA, DEC, ln D_L, ln chirpmass, eta, chi1, chi2, \sqrt \alpha (in km) (the coupling parameter)

IMRPhenomPv2 - 9 dimensions — cos J_N, ln chirpmass, eta, |chi1|, |chi1|, theta_1, theta_2, phi_1, phi_2

Parameters:

statistics_filename – Filename to output sampling statistics, if empty string, not output
chain_filename – Filename to output data (chain 0 only), if empty string, not output
auto_corr_filename – Filename to output auto correlation in some interval, if empty string, not output
checkpoint_filename – Filename to output data for checkpoint, if empty string, not saved

void continue_PTMCMC_MH_dynamic_PT_alloc_uncorrelated_GW(std::string checkpoint_file_start, mcmc_sampler_output *sampler_output, double **output, int N_steps, int max_chain_N_thermo_ensemble, double *chain_temps, int swp_freq, int t0, int nu, int max_chunk_size, std::string chain_distribution_scheme, double (*log_prior)(double *param, mcmc_data_interface *interface, void *parameters), int numThreads, bool pool, bool show_prog, int num_detectors, std::complex<double> **data, double **noise_psd, double **frequencies, int *data_length, double gps_time, std::string *detectors, MCMC_modification_struct *mod_struct, std::string generation_method, std::string statistics_filename, std::string chain_filename, std::string likelihood_log_filename, std::string checkpoint_filename)

Takes in an MCMC checkpoint file and continues the chain.

Obviously, the user must be sure to correctly match the dimension, number of chains, the generation_method, the prior function, the data, psds, freqs, and the detectors (number and name), and the gps_time to the previous run, otherwise the behavior of the sampler is undefined.

numThreads and pool do not necessarily have to be the same

void PTMCMC_MH_dynamic_PT_alloc_uncorrelated_GW(mcmc_sampler_output *sampler_output, double **output, int dimension, int N_steps, int chain_N, int max_chain_N_thermo_ensemble, double *initial_pos, double *seeding_var, double **ensemble_initial_pos, double *chain_temps, int swp_freq, int t0, int nu, int max_chunk_size, std::string chain_distribution_scheme, double (*log_prior)(double *param, mcmc_data_interface *interface, void *parameters), int numThreads, bool pool, bool show_prog, int num_detectors, std::complex<double> **data, double **noise_psd, double **frequencies, int *data_length, double gps_time, std::string *detectors, MCMC_modification_struct *mod_struct, std::string generation_method, std::string statistics_filename, std::string chain_filename, std::string likelihood_log_filename, std::string checkpoint_filename)

Takes in an MCMC checkpoint file and continues the chain.

Obviously, the user must be sure to correctly match the dimension, number of chains, the generation_method, the prior function, the data, psds, freqs, and the detectors (number and name), and the gps_time to the previous run, otherwise the behavior of the sampler is undefined.

numThreads and pool do not necessarily have to be the same

void PTMCMC_MH_dynamic_PT_alloc_GW(double ***output, int dimension, int N_steps, int chain_N, int max_chain_N_thermo_ensemble, double *initial_pos, double *seeding_var, double **ensemble_initial_pos, double *chain_temps, int swp_freq, int t0, int nu, std::string chain_distribution_scheme, double (*log_prior)(double *param, mcmc_data_interface *interface, void *parameters), int numThreads, bool pool, bool show_prog, int num_detectors, std::complex<double> **data, double **noise_psd, double **frequencies, int *data_length, double gps_time, std::string *detectors, MCMC_modification_struct *mod_struct, std::string generation_method, std::string statistics_filename, std::string chain_filename, std::string likelihood_log_filename, std::string checkpoint_filename)

Takes in an MCMC checkpoint file and continues the chain.

Obviously, the user must be sure to correctly match the dimension, number of chains, the generation_method, the prior function, the data, psds, freqs, and the detectors (number and name), and the gps_time to the previous run, otherwise the behavior of the sampler is undefined.

numThreads and pool do not necessarily have to be the same

void continue_PTMCMC_MH_GW(std::string start_checkpoint_file, double ***output, int dimension, int N_steps, int swp_freq, double (*log_prior)(double *param, mcmc_data_interface *interface, void *parameters), int numThreads, bool pool, bool show_prog, int num_detectors, std::complex<double> **data, double **noise_psd, double **frequencies, int *data_length, double gps_time, std::string *detector, MCMC_modification_struct *mod_struct, std::string generation_method, std::string statistics_filename, std::string chain_filename, std::string auto_corr_filename, std::string likelihood_log_filename, std::string final_checkpoint_filename, bool tune)

Takes in an MCMC checkpoint file and continues the chain.

Obviously, the user must be sure to correctly match the dimension, number of chains, the generation_method, the prior function, the data, psds, freqs, and the detectors (number and name), and the gps_time to the previous run, otherwise the behavior of the sampler is undefined.

numThreads and pool do not necessarily have to be the same

void PTMCMC_method_specific_prep(std::string generation_method, int dimension, double **seeding_var_ptr, bool local_seeding)

Unpacks MCMC parameters for method specific initiation.

Populates seeding vector if non supplied, populates mcmc_Nmod, populates mcmc_log_beta, populates mcmc_intrinsic

void RJPTMCMC_method_specific_prep(std::string generation_method, int max_dim, int min_dim, double *seeding_var, bool local_seeding)

Unpacks MCMC parameters for method specific initiation (RJ version)

Populates seeding vector if non supplied, populates mcmc_Nmod, populates mcmc_log_beta, populates mcmc_intrinsic

double MCMC_likelihood_extrinsic(bool save_waveform, gen_params_base<double> *parameters, std::string generation_method, int *data_length, double **frequencies, std::complex<double> **data, double **psd, double **weights, std::string integration_method, bool log10F, std::string *detectors, int num_detectors, Quadrature *QuadMethod = NULL)

void MCMC_fisher_wrapper(double *param, double **output, mcmc_data_interface *interface, void *parameters)

std::string MCMC_prep_params(double *param, double *temp_params, gen_params_base<double> *gen_params, int dimension, std::string generation_method, MCMC_modification_struct *mod_struct)

utility to do MCMC specific transformations on the input param vector before passing to the repacking utillity

Returns the local generation method to be used in the LL functions

double MCMC_likelihood_wrapper(double *param, mcmc_data_interface *interface, void *parameters)

void continue_RJPTMCMC_MH_dynamic_PT_alloc_comprehensive_2WF_GW(std::string checkpoint_file_start, mcmc_sampler_output *sampler_output, double **output, int **status, int *model_status, int nested_model_number, int N_steps, int max_chain_N_thermo_ensemble, double **prior_ranges, double *chain_temps, int swp_freq, int t0, int nu, int max_chunk_size, std::string chain_distribution_scheme, double (*log_prior)(double *param, int *status, int model_status, mcmc_data_interface *interface, void *parameters), int numThreads, bool pool, bool show_prog, int num_detectors, std::complex<double> **data, double **noise_psd, double **frequencies, int *data_length, double gps_time, std::string *detectors, MCMC_modification_struct *mod_struct, std::string generation_method_base, std::string generation_method_extended, std::string statistics_filename, std::string chain_filename, std::string likelihood_log_filename, std::string checkpoint_filename)

Performs RJPTMCMC on two competeting models, waveform 1 and waveform 2. Assumes waveform 2 is a superset including waveform 1 (ie dCS Pv2 and Pv2 or NRT Pv2 and Pv2, but not dCS Pv2 and D).

Parameters:: prior_ranges – Range of priors on MODIFICATIONS — shape : double[N_mods][2] with low in 0 and high in 1*

void RJPTMCMC_MH_dynamic_PT_alloc_comprehensive_2WF_GW(mcmc_sampler_output *sampler_output, double **output, int **status, int *model_status, int nested_model_number, int max_dimension, int min_dimension, int N_steps, int chain_N, int max_chain_N_thermo_ensemble, double *initial_pos, int *initial_status, int initial_model_status, double *seeding_var, double **ensemble_initial_pos, int **ensemble_initial_status, int *ensemble_initial_model_status, double **prior_ranges, double *chain_temps, int swp_freq, int t0, int nu, int max_chunk_size, std::string chain_distribution_scheme, double (*log_prior)(double *param, int *status, int model_status, mcmc_data_interface *interface, void *parameters), int numThreads, bool pool, bool show_prog, int num_detectors, std::complex<double> **data, double **noise_psd, double **frequencies, int *data_length, double gps_time, std::string *detectors, MCMC_modification_struct *mod_struct, std::string generation_method_base, std::string generation_method_extended, std::string statistics_filename, std::string chain_filename, std::string likelihood_log_filename, std::string checkpoint_filename)

Performs RJPTMCMC on two competeting models, waveform 1 and waveform 2. Assumes waveform 2 is a superset including waveform 1 (ie dCS Pv2 and Pv2 or NRT Pv2 and Pv2, but not dCS Pv2 and D).

Parameters:: prior_ranges – Range of priors on MODIFICATIONS — shape : double[N_mods][2] with low in 0 and high in 1*

void RJMCMC_2WF_RJ_proposal_wrapper(double *current_param, double *proposed_param, int *current_status, int *proposed_status, int *current_model_status, int *proposed_model_status, double *MH_corrections, mcmc_data_interface *interface, void *user_param)

void RJMCMC_2WF_fisher_wrapper(double *param, int *status, int model_status, double **fisher, mcmc_data_interface *interface, void *user_param)

double RJMCMC_2WF_likelihood_wrapper(double *param, int *status, int model_status, mcmc_data_interface *interface, void *user_param)

void pack_local_mod_structure(mcmc_data_interface *interface, double *param, int *status, std::string waveform_extended, void *parameters, MCMC_modification_struct *full_struct, MCMC_modification_struct *local_struct)

void MCMC_fisher_transformations(double *param, double **fisher, int dimension, std::string generation_method, bool intrinsic, mcmc_data_interface *interface, MCMC_modification_struct *mod_struct, void *parameters)

Variables

static double **mcmc_noise = NULL

static double *mcmc_init_pos = NULL

static std::complex<double> **mcmc_data = NULL

static double **mcmc_frequencies = NULL

static std::string *mcmc_detectors = NULL

static std::string mcmc_generation_method = ""

static std::string mcmc_generation_method_base = ""

static std::string mcmc_generation_method_extended = ""

static int *mcmc_data_length = NULL

static fftw_outline *mcmc_fftw_plans = NULL

static int mcmc_num_detectors = 2

static double mcmc_gps_time = 0

static double mcmc_gmst = 0

static int mcmc_max_dim

static int mcmc_min_dim

static int mcmc_Nmod

static int mcmc_Nmod_max

static double *mcmc_bppe

static gsl_interp_accel **mcmc_accels = NULL

static gsl_spline **mcmc_splines = NULL

static bool mcmc_log_beta

static bool mcmc_intrinsic

static bool mcmc_save_waveform

static int mcmc_deriv_order = 4

static gsl_rng **mcmc_rvec

static MCMC_modification_struct *mcmc_mod_struct

struct MCMC_modification_struct

#include <mcmc_gw.h>

Public Members

int ppE_Nmod = 0

double *bppe = NULL

int gIMR_Nmod_phi = 0

int *gIMR_phii = NULL

int gIMR_Nmod_sigma = 0

int *gIMR_sigmai = NULL

int gIMR_Nmod_beta = 0

int *gIMR_betai = NULL

int gIMR_Nmod_alpha = 0

int *gIMR_alphai = NULL

bool NSflag1 = false

bool NSflag2 = false

bool tidal_love = true

bool tidal_love_error = false

bool alpha_param = true

bool EA_region1 = false

bool GAUSS_QUAD = false

bool log10F = false

double **weights = NULL

double **fisher_freq = NULL

double **fisher_weights = NULL

double **fisher_PSD = NULL

bool fisher_GAUSS_QUAD = false

int *fisher_length = NULL

bool fisher_log10F = false

Quadrature *QuadMethod = NULL

AdaptiveLikelihood *adaptivell = NULL

double f_ref = 20.

struct MCMC_user_param

#include <mcmc_gw.h>

Public Members

std::complex<double> **burn_data = NULL

double **burn_noise = NULL

double **burn_freqs = NULL

int *burn_lengths = NULL

fftw_outline *burn_plans = NULL

std::mutex *mFish

bool GAUSS_QUAD = false

bool log10F = false

double **weights = NULL

double **fisher_freq = NULL

double **fisher_weights = NULL

double **fisher_PSD = NULL

bool fisher_AD = false

bool fisher_GAUSS_QUAD = false

int *fisher_length = NULL

bool fisher_log10F = false

Quadrature *QuadMethod = NULL

double **mod_prior_ranges = NULL

MCMC_modification_struct *mod_struct

gsl_rng *r