semicoherent_glitch_robust_directed_MCMC_search_on_1_glitch.py 2.2 KB
Newer Older
Gregory Ashton's avatar
Gregory Ashton committed
1 2 3
import numpy as np
import matplotlib.pyplot as plt
import pyfstat
4
import gridcorner
Gregory Ashton's avatar
Gregory Ashton committed
5
import time
Gregory Ashton's avatar
Gregory Ashton committed
6
from make_simulated_data import tstart, duration, tref, F0, F1, F2, Alpha, Delta, delta_F0, dtglitch, outdir
Gregory Ashton's avatar
Gregory Ashton committed
7

8
plt.style.use('./paper.mplstyle')
Gregory Ashton's avatar
Gregory Ashton committed
9 10 11

label = 'semicoherent_glitch_robust_directed_MCMC_search_on_1_glitch'

12
Nstar = 1000
Gregory Ashton's avatar
Gregory Ashton committed
13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
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,
    }

35
ntemps = 3
Gregory Ashton's avatar
Gregory Ashton committed
36 37
log10beta_min = -0.5
nwalkers = 100
Gregory Ashton's avatar
Gregory Ashton committed
38
nsteps = [250, 250]
Gregory Ashton's avatar
Gregory Ashton committed
39 40 41 42 43 44 45

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(
46 47
    subtractor=F0, multiplier=1e6, symbol='$f-f_\mathrm{s}$')
mcmc.unit_dictionary['F0'] = '$\mu$Hz'
Gregory Ashton's avatar
Gregory Ashton committed
48
mcmc.transform_dictionary['F1'] = dict(
49 50 51 52 53 54 55 56
    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]'
Gregory Ashton's avatar
Gregory Ashton committed
57

Gregory Ashton's avatar
Gregory Ashton committed
58
t1 = time.time()
Gregory Ashton's avatar
Gregory Ashton committed
59
mcmc.run()
Gregory Ashton's avatar
Gregory Ashton committed
60
dT = time.time() - t1
61
fig_and_axes = gridcorner._get_fig_and_axes(4, 2, 0.05)
62
mcmc.plot_corner(label_offset=0.25, truths=[0, 0, 0, 0],
63
                 fig_and_axes=fig_and_axes)
Gregory Ashton's avatar
Gregory Ashton committed
64
mcmc.print_summary()
Gregory Ashton's avatar
Gregory Ashton committed
65 66 67

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