Select Git revision
Forked from
Gregory Ashton / PyFstat
Source project has a limited visibility.
semicoherent_glitch_robust_directed_MCMC_search_on_1_glitch.py 2.34 KiB
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.0, "upper": F0 + F0_width / 2.0},
"F1": {"type": "unif", "lower": F1 - F1_width / 2.0, "upper": F1 + F1_width / 2.0},
"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),
hist_kwargs=dict(lw=1.5, zorder=-1),
truth_color="C3",
)
mcmc.print_summary()
print(("Prior widths =", F0_width, F1_width))
print(("Actual run time = {}".format(dT)))