Commit c9d1bca8 authored by Gregory Ashton's avatar Gregory Ashton

Update glitch examples

parent 11f162d0
......@@ -17,7 +17,7 @@ h0 = 5e-24
# Properties of the GW data
sqrtSX = 1e-22
tstart = 1000000000
duration = 100*86400
duration = 50*86400
tend = tstart+duration
tref = tstart + 0.5*duration
......@@ -26,10 +26,6 @@ data = Writer(
F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX)
data.make_data()
# The predicted twoF, given by lalapps_predictFstat can be accessed by
twoF = data.predict_fstat()
print 'Predicted twoF value: {}\n'.format(twoF)
# Next, taking the same signal parameters, we include a glitch half way through
dtglitch = duration/2.0
delta_F0 = 5e-6
......
......@@ -2,6 +2,7 @@ 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')
......@@ -34,21 +35,25 @@ theta_prior = {
ntemps = 3
log10beta_min = -0.5
nwalkers = 100
nsteps = [500, 1000]
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)
print delta_F0
mcmc.transform_dictionary['F0'] = dict(
subtractor=F0, symbol='$f-f^\mathrm{s}$')
mcmc.transform_dictionary['F1'] = dict(
subtractor=F1, symbol='$\dot{f}-\dot{f}^\mathrm{s}$')
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.35, truths=[0, 0, delta_F0, 50],
fig_and_axes=fig_and_axes)
mcmc.print_summary()
print('Prior widths =', F0_width, F1_width)
print("Actual run time = {}".format(dT))
import pyfstat
import numpy as np
import matplotlib.pyplot as plt
from make_simulated_data import tstart, duration, tref, F0, F1, F2, Alpha, Delta, delta_F0, outdir
from make_simulated_data import tstart, duration, tref, F0, F1, F2, Alpha, Delta, delta_F0, outdir, dtglitch
import time
try:
from gridcorner import gridcorner
......@@ -17,7 +18,7 @@ plt.style.use('./paper.mplstyle')
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)
N = 30
N = 20
F0s = [F0-F0_width/2., F0+F0_width/2., F0_width/N]
F1s = [F1-F1_width/2., F1+F1_width/2., F1_width/N]
F2s = [F2]
......@@ -29,14 +30,15 @@ tglitchs = [tstart+0.1*duration, tstart+0.9*duration, 0.8*float(duration)/N]
delta_F0s = [0, max_delta_F0, max_delta_F0/N]
delta_F1s = [0]
print 'Prior widths=', F0_width, F1_width
t1 = time.time()
search = pyfstat.GridGlitchSearch(
label, outdir, 'data/*1_glitch*sft', F0s=F0s, F1s=F1s, F2s=F2s,
Alphas=Alphas, Deltas=Deltas, tref=tref, minStartTime=tstart,
maxStartTime=tstart+duration, tglitchs=tglitchs, delta_F0s=delta_F0s,
delta_F1s=delta_F1s)
search.run()
dT = time.time() - t1
F0_vals = np.unique(search.data[:, 0]) - F0
F1_vals = np.unique(search.data[:, 1]) - F1
......@@ -51,6 +53,11 @@ labels = ['$f - f^\mathrm{s}$\n[Hz]', '$\dot{f} - \dot{f}^\mathrm{s}$\n[Hz/s]',
'$\delta f$\n[Hz]', '$t^g_0$\n[days]', '$\widehat{2\mathcal{F}}$']
fig, axes = gridcorner(
twoF, xyz, projection='log_mean', labels=labels,
showDvals=False, lines=[0, 0, delta_F0, 50], label_offset=0.35)
showDvals=False, lines=[0, 0, delta_F0, dtglitch/86400.], label_offset=0.35)
fig.savefig('{}/{}_projection_matrix.png'.format(outdir, label),
bbox_inches='tight')
print('Prior widths =', F0_width, F1_width)
print("Actual run time = {}".format(dT))
print("Actual number of grid points = {}".format(search.data.shape[0]))
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment