semicoherent_glitch_robust_directed_MCMC_search_on_1_glitch.py

import numpy as np
import matplotlib.pyplot as plt
import pyfstat
import gridcorner
import time
from make_simulated_data import tstart, duration, tref, F0, F1, F2, Alpha, Delta, delta_F0, dtglitch, outdir

plt.style.use('./paper.mplstyle')

label = 'semicoherent_glitch_robust_directed_MCMC_search_on_1_glitch'

Nstar = 1000
F0_width = np.sqrt(Nstar)*np.sqrt(12)/(np.pi*duration)
F1_width = np.sqrt(Nstar)*np.sqrt(180)/(np.pi*duration**2)

theta_prior = {
    'F0': {'type': 'unif',
           'lower': F0-F0_width/2.,
           'upper': F0+F0_width/2.},
    'F1': {'type': 'unif',
           'lower': F1-F1_width/2.,
           'upper': F1+F1_width/2.},
    'F2': F2,
    'delta_F0': {'type': 'unif',
                 'lower': 0,
                 'upper': 1e-5},
    'delta_F1': 0,
    'tglitch': {'type': 'unif',
                'lower': tstart+0.1*duration,
                'upper': tstart+0.9*duration},
    'Alpha': Alpha,
    'Delta': Delta,
    }

ntemps = 3
log10beta_min = -0.5
nwalkers = 100
nsteps = [250, 250]

mcmc = pyfstat.MCMCGlitchSearch(
    label=label, sftfilepattern='data/*1_glitch*sft', theta_prior=theta_prior,
    tref=tref, minStartTime=tstart, maxStartTime=tstart+duration,
    nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps,
    log10beta_min=log10beta_min, nglitch=1)
mcmc.transform_dictionary['F0'] = dict(
    subtractor=F0, multiplier=1e6, symbol='$f-f_\mathrm{s}$')
mcmc.unit_dictionary['F0'] = '$\mu$Hz'
mcmc.transform_dictionary['F1'] = dict(
    subtractor=F1, multiplier=1e12, symbol='$\dot{f}-\dot{f}_\mathrm{s}$')
mcmc.unit_dictionary['F1'] = '$p$Hz/s'
mcmc.transform_dictionary['delta_F0'] = dict(
    multiplier=1e6, subtractor=delta_F0,
    symbol='$\delta f-\delta f_\mathrm{s}$')
mcmc.unit_dictionary['delta_F0'] = '$\mu$Hz/s'
mcmc.transform_dictionary['tglitch']['subtractor'] = tstart + dtglitch
mcmc.transform_dictionary['tglitch']['label'] ='$t^\mathrm{g}-t^\mathrm{g}_\mathrm{s}$\n[d]'

t1 = time.time()
mcmc.run()
dT = time.time() - t1
fig_and_axes = gridcorner._get_fig_and_axes(4, 2, 0.05)
mcmc.plot_corner(label_offset=0.25, truths=[0, 0, 0, 0],
                 fig_and_axes=fig_and_axes, quantiles=(0.16, 0.84))

mcmc.print_summary()

print('Prior widths =', F0_width, F1_width)
print("Actual run time = {}".format(dT))