Skip to content
Snippets Groups Projects
Commit b6403cc8 authored by Gregory Ashton's avatar Gregory Ashton
Browse files

Update glitch examples

parent e67d91d6
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 155 additions and 219 deletions
examples/glitch_examples/fully_coherent_search_using_MCMC_on_glitching_data.py deleted 100644 → 0
+ 0
− 35
View file @ e67d91d6
import numpy as np
from pyfstat import MCMCSearch
F0 = 30.0
F1 = -1e-10
F2 = 0
Alpha = np.radians(83.6292)
Delta = np.radians(22.0144)
tref = 362750407.0
tstart = 1000000000
duration = 100*86400
tend = tstart + duration
theta_prior = {'F0': {'type': 'unif', 'lower': F0-1e-4, 'upper': F0+1e-4},
'F1': {'type': 'unif', 'lower': F1*(1+1e-3), 'upper': F1*(1-1e-3)},
'F2': F2,
'Alpha': Alpha,
'Delta': Delta
}
ntemps = 2
log10beta_min = -0.01
nwalkers = 100
nsteps = [500, 500]
mcmc = MCMCSearch('fully_coherent_search_using_MCMC_on_glitching_data', 'data',
sftfilepattern='data/*_glitch*.sft',
theta_prior=theta_prior, tref=tref, minStartTime=tstart, maxStartTime=tend,
nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps,
log10beta_min=log10beta_min)
mcmc.run()
mcmc.plot_corner(add_prior=True)
mcmc.print_summary()
examples/glitch_examples/glitch_robust_search.py deleted 100644 → 0
+ 0
− 54
View file @ e67d91d6
import numpy as np
import pyfstat
outdir = 'data'
label = 'glitch_robust_search'
# Properties of the GW data
tstart = 1000000000
Tspan = 60 * 86400
# Fixed properties of the signal
F0s = 30
F1s = -1e-8
F2s = 0
Alpha = np.radians(83.6292)
Delta = np.radians(22.0144)
tref = tstart + .5 * Tspan
sftfilepattern = 'data/*glitching_signal*sft'
F0_width = np.sqrt(3)/(np.pi*Tspan)
F1_width = np.sqrt(45/4.)/(np.pi*Tspan**2)
DeltaF0 = 50 * F0_width
DeltaF1 = 50 * F1_width
theta_prior = {'F0': {'type': 'unif',
'lower': F0s-DeltaF0,
'upper': F0s+DeltaF0},
'F1': {'type': 'unif',
'lower': F1s-DeltaF1,
'upper': F1s+DeltaF1},
'F2': F2s,
'delta_F0': {'type': 'unif',
'lower': 0,
'upper': 1e-5},
'delta_F1': {'type': 'unif',
'lower': -1e-11,
'upper': 1e-11},
'tglitch': {'type': 'unif',
'lower': tstart+0.1*Tspan,
'upper': tstart+0.9*Tspan},
'Alpha': Alpha,
'Delta': Delta,
}
search = pyfstat.MCMCGlitchSearch(
label=label, outdir=outdir, sftfilepattern=sftfilepattern,
theta_prior=theta_prior, nglitch=1, tref=tref, nsteps=[500, 500],
ntemps=3, log10beta_min=-0.5, minStartTime=tstart,
maxStartTime=tstart+Tspan)
search.run()
search.plot_corner(label_offset=0.8, add_prior=True)
search.print_summary()
examples/glitch_examples/glitch_robust_search_make_simulated_data.py deleted 100644 → 0
+ 0
− 40
View file @ e67d91d6
import numpy as np
import pyfstat
outdir = 'data'
label = 'simulated_glitching_signal'
# Properties of the GW data
tstart = 1000000000
Tspan = 60 * 86400
tref = tstart + .5 * Tspan
# Fixed properties of the signal
F0s = 30
F1s = -1e-8
F2s = 0
Alpha = np.radians(83.6292)
Delta = np.radians(22.0144)
h0 = 1e-25
sqrtSX = 1e-24
psi = -0.1
phi = 0
cosi = 0.5
# Glitch properties
dtglitch = 0.45 * Tspan # time (in secs) after minStartTime
dF0 = 5e-6
dF1 = 1e-12
detectors = 'H1'
glitch_data = pyfstat.Writer(
label=label, outdir=outdir, tref=tref, tstart=tstart,
F0=F0s, F1=F1s, F2=F2s, duration=Tspan, Alpha=Alpha,
Delta=Delta, sqrtSX=sqrtSX, dtglitch=dtglitch,
h0=h0, cosi=cosi, phi=phi, psi=psi,
delta_F0=dF0, delta_F1=dF1, add_noise=True)
glitch_data.make_data()
examples/glitch_examples/make_fake_data.py → examples/glitch_examples/make_simulated_data.py
+ 7
− 6
View file @ b6403cc8
from pyfstat import Writer, GlitchWriter
import numpy as np
outdir = 'data'
# First, we generate data with a reasonably strong smooth signal
# Define parameters of the Crab pulsar as an example
......@@ -9,19 +10,19 @@ F1 = -1e-10
F2 = 0
Alpha = np.radians(83.6292)
Delta = np.radians(22.0144)
tref = 362750407.0
# Signal strength
h0 = 1e-23
h0 = 5e-24
# Properties of the GW data
sqrtSX = 1e-22
tstart = 1000000000
duration = 100*86400
tend = tstart+duration
tref = tstart + 0.5*duration
data = Writer(
label='basic', outdir='data', tref=tref, tstart=tstart, F0=F0, F1=F1,
label='0_glitch', outdir=outdir, tref=tref, tstart=tstart, F0=F0, F1=F1,
F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX)
data.make_data()
......@@ -31,11 +32,11 @@ 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 = 4e-5
delta_F0 = 5e-6
delta_F1 = 0
glitch_data = GlitchWriter(
label='glitch', outdir='data', tref=tref, tstart=tstart, F0=F0, F1=F1,
label='1_glitch', outdir=outdir, tref=tref, tstart=tstart, F0=F0, F1=F1,
F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX,
dtglitch=dtglitch, delta_F0=delta_F0, delta_F1=delta_F1)
glitch_data.make_data()
......@@ -49,7 +50,7 @@ delta_F1 = [0, 0]
delta_F2 = [0, 0]
two_glitch_data = GlitchWriter(
label='twoglitch', outdir='data', tref=tref, tstart=tstart, F0=F0, F1=F1,
label='2_glitch', outdir=outdir, tref=tref, tstart=tstart, F0=F0, F1=F1,
F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX,
dtglitch=dtglitch, delta_phi=delta_phi, delta_F0=delta_F0,
delta_F1=delta_F1, delta_F2=delta_F2)
......
examples/glitch_examples/semi_coherent_glitch_search_using_MCMC.py deleted 100644 → 0
+ 0
− 41
View file @ e67d91d6
import pyfstat
F0 = 30.0
F1 = -1e-10
F2 = 0
Alpha = 5e-3
Delta = 6e-2
tref = 362750407.0
tstart = 1000000000
duration = 100*86400
tend = tstart + duration
theta_prior = {'F0': {'type': 'norm', 'loc': F0, 'scale': abs(1e-6*F0)},
'F1': {'type': 'norm', 'loc': F1, 'scale': abs(1e-6*F1)},
'F2': F2,
'Alpha': Alpha,
'Delta': Delta,
'delta_F0': {'type': 'halfnorm', 'loc': 0,
'scale': 1e-5*F0},
'delta_F1': 0,
'tglitch': {'type': 'unif',
'lower': tstart+0.1*duration,
'upper': tstart+0.9*duration},
}
ntemps = 4
log10beta_min = -1
nwalkers = 100
nsteps = [5000, 1000, 1000]
mcmc = pyfstat.MCMCGlitchSearch(
'semi_coherent_glitch_search_using_MCMC', 'data',
sftfilepattern='data/*_glitch*sft', theta_prior=theta_prior, tref=tref,
tstart=tstart, tend=tend, nsteps=nsteps, nwalkers=nwalkers,
scatter_val=1e-10, nglitch=1, ntemps=ntemps,
log10beta_min=log10beta_min)
mcmc.run()
mcmc.plot_corner(add_prior=True)
mcmc.print_summary()
examples/glitch_examples/semi_coherent_twoglitch_search_using_MCMC.py deleted 100644 → 0
+ 0
− 43
View file @ e67d91d6
import pyfstat
F0 = 30.0
F1 = -1e-10
F2 = 0
Alpha = 5e-3
Delta = 6e-2
tref = 362750407.0
tstart = 1000000000
duration = 100*86400
tend = tstart + duration
theta_prior = {'F0': {'type': 'norm', 'loc': F0, 'scale': abs(1e-6*F0)},
'F1': {'type': 'norm', 'loc': F1, 'scale': abs(1e-6*F1)},
'F2': F2,
'Alpha': Alpha,
'Delta': Delta,
'delta_F0_0': {'type': 'halfnorm', 'loc': 0,
'scale': 1e-7*F0},
'delta_F0_1': {'type': 'halfnorm', 'loc': 0,
'scale': 1e-7*F0},
'delta_F1_0': 0,
'delta_F1_1': 0,
'tglitch_0': {'type': 'unif',
'lower': tstart+0.01*duration,
'upper': tstart+0.5*duration},
'tglitch_1': {'type': 'unif',
'lower': tstart+0.5*duration,
'upper': tstart+0.99*duration},
}
nwalkers = 100
nsteps = [1000, 1000, 5000]
mcmc = pyfstat.MCMCGlitchSearch(
'semi_coherent_twoglitch_search', 'data', sftfilepattern='data/*twoglitch*sft',
theta_prior=theta_prior, tref=tref, tstart=tstart,
tend=tend, nsteps=nsteps, nwalkers=nwalkers, scatter_val=1e-10, nglitch=2)
mcmc.run()
mcmc.plot_corner(add_prior=True, tglitch_ratio=True, figsize=(10, 10))
mcmc.print_summary()
examples/glitch_examples/semicoherent_glitch_robust_directed_MCMC_search_on_1_glitch.py 0 → 100644
+ 51
− 0
View file @ b6403cc8
import numpy as np
import matplotlib.pyplot as plt
import pyfstat
from make_simulated_data import tstart, duration, tref, F0, F1, F2, Alpha, Delta, outdir
plt.style.use('paper')
label = 'semicoherent_glitch_robust_directed_MCMC_search_on_1_glitch'
Nstar = 100
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 = 2
log10beta_min = -0.5
nwalkers = 100
nsteps = [500, 500]
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, symbol='$f-f^\mathrm{s}$')
mcmc.transform_dictionary['F1'] = dict(
subtractor=F1, symbol='$\dot{f}-\dot{f}^\mathrm{s}$')
mcmc.run()
mcmc.plot_corner()
mcmc.print_summary()
examples/glitch_examples/semicoherent_glitch_robust_directed_grid_search_on_1_glitch.py 0 → 100644
+ 53
− 0
View file @ b6403cc8
import pyfstat
import numpy as np
import matplotlib.pyplot as plt
from make_simulated_data import tstart, duration, tref, F0, F1, F2, Alpha, Delta, outdir
try:
from gridcorner import gridcorner
except ImportError:
raise ImportError(
"Python module 'gridcorner' not found, please install from "
"https://gitlab.aei.uni-hannover.de/GregAshton/gridcorner")
label = 'semicoherent_glitch_robust_directed_grid_search_on_1_glitch'
plt.style.use('paper')
Nstar = 100
F0_width = np.sqrt(Nstar)*np.sqrt(12)/(np.pi*duration)
F1_width = np.sqrt(Nstar)*np.sqrt(180)/(np.pi*duration**2)
N = 41
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]
Alphas = [Alpha]
Deltas = [Delta]
max_delta_F0 = 1e-5
tglitchs = [tstart+0.1*duration, tstart+0.9*duration, float(duration)/N]
delta_F0s = [0, max_delta_F0, max_delta_F0/N]
delta_F1s = [0]
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()
F0_vals = np.unique(search.data[:, 0]) - F0
F1_vals = np.unique(search.data[:, 1]) - F1
delta_F0s_vals = np.unique(search.data[:, 5])
tglitch_vals = np.unique(search.data[:, 7])
tglitch_vals_days = (tglitch_vals-tstart) / 86400.
twoF = search.data[:, -1].reshape((len(F0_vals), len(F1_vals),
len(delta_F0s_vals), len(tglitch_vals)))
xyz = [F0_vals, F1_vals, delta_F0s_vals, tglitch_vals_days]
labels = ['$f - f_0$\n[Hz]', '$\dot{f} - \dot{f}_0$\n[Hz/s]',
'$\delta f$\n[Hz]', '$t^g_0$\n[days]', '$\widehat{2\mathcal{F}}$']
fig, axes = gridcorner(
twoF, xyz, projection='log_mean', whspace=0.1, factor=1.2, labels=labels)
fig.savefig('{}/{}_projection_matrix.png'.format(outdir, label),
bbox_inches='tight')
examples/glitch_examples/standard_directed_MCMC_search_on_1_glitch.py 0 → 100644
+ 44
− 0
View file @ b6403cc8
import numpy as np
import matplotlib.pyplot as plt
import pyfstat
from make_simulated_data import tstart, duration, tref, F0, F1, F2, Alpha, Delta, outdir
plt.style.use('paper')
label = 'standard_directed_MCMC_search_on_1_glitch'
Nstar = 10000
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,
'Alpha': Alpha,
'Delta': Delta,
}
ntemps = 2
log10beta_min = -0.5
nwalkers = 100
nsteps = [500, 2000]
mcmc = pyfstat.MCMCSearch(
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)
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}$')
mcmc.run()
mcmc.plot_corner()
mcmc.print_summary()
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment

AEI Privacy Policy Imprint