Commit da2b4ca8 authored by Gregory Ashton's avatar Gregory Ashton
Browse files

Adds glitch text and all examples to folder

parent b06a2d89
frequency_grid_files = grided_frequency_search_1D.png
fully_coherent_files = fully_coherent_search_using_MCMC_walkers.png
follow_up_files = follow_up_run_setup.tex follow_up_walkers.png
transient_files = transient_search_initial_stage_twoFcumulative.png transient_search_corner.png
glitch_files = single_glitch_F0F1_grid_2D.png single_glitch_glitchSearch_corner.png
all_files = $(frequency_grid_files) $(fully_coherent_files) $(follow_up_files) $(transient_files) $(glitch_files)
cd data; cp $(all_files) ../../
import pyfstat
import numpy as np
# Properties of the GW data
sqrtSX = 1e-23
tstart = 1000000000
duration = 100*86400
tend = tstart + duration
# Properties of the signal
F0 = 30.0
F1 = -1e-10
F2 = 0
Alpha = 5e-3
Delta = 6e-2
tref = .5*(tstart+tend)
depth = 10
h0 = sqrtSX / depth
data_label = 'fully_coherent_search_using_MCMC'
data = pyfstat.Writer(
label=data_label, outdir='data', tref=tref,
tstart=tstart, F0=F0, F1=F1, F2=F2, duration=duration, Alpha=Alpha,
Delta=Delta, h0=h0, sqrtSX=sqrtSX)
# The predicted twoF, given by lalapps_predictFstat can be accessed by
twoF = data.predict_fstat()
print 'Predicted twoF value: {}\n'.format(twoF)
DeltaF0 = 1e-7
DeltaF1 = 1e-13
VF0 = (np.pi * duration * DeltaF0)**2 / 3.0
VF1 = (np.pi * duration**2 * DeltaF1)**2 * 4/45.
print '\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n'.format(VF0*VF1, VF0, VF1)
theta_prior = {'F0': {'type': 'unif',
'lower': F0-DeltaF0/2.,
'upper': F0+DeltaF0/2.},
'F1': {'type': 'unif',
'lower': F1-DeltaF1/2.,
'upper': F1+DeltaF1/2.},
'F2': F2,
'Alpha': Alpha,
'Delta': Delta
ntemps = 1
log10temperature_min = -1
nwalkers = 1000
nsteps = [50, 50]
mcmc = pyfstat.MCMCSearch(
label='fully_coherent_search_using_MCMC', outdir='data',
sftfilepath='data/*'+data_label+'*sft', theta_prior=theta_prior, tref=tref,
minStartTime=tstart, maxStartTime=tend, nsteps=nsteps, nwalkers=nwalkers,
ntemps=ntemps, log10temperature_min=log10temperature_min)'paper', subtractions=[30, -1e-10])
import pyfstat
import numpy as np
import matplotlib.pyplot as plt'paper')
F0 = 30.0
F1 = 0
F2 = 0
Alpha = 1.0
Delta = 1.5
# Properties of the GW data
sqrtSX = 1e-23
tstart = 1000000000
duration = 100*86400
tend = tstart+duration
tref = .5*(tstart+tend)
depth = 70
data_label = 'grided_frequency_depth_{:1.0f}'.format(depth)
h0 = sqrtSX / depth
data = pyfstat.Writer(
label=data_label, outdir='data', tref=tref,
tstart=tstart, F0=F0, F1=F1, F2=F2, duration=duration, Alpha=Alpha,
Delta=Delta, h0=h0, sqrtSX=sqrtSX)
m = 0.001
dF0 = np.sqrt(12*m)/(np.pi*duration)
DeltaF0 = 800*dF0
F0s = [F0-DeltaF0/2., F0+DeltaF0/2., dF0]
F1s = [F1]
F2s = [F2]
Alphas = [Alpha]
Deltas = [Delta]
search = pyfstat.GridSearch(
'grided_frequency_search', 'data', 'data/*'+data_label+'*sft', F0s, F1s,
F2s, Alphas, Deltas, tref, tstart, tend)
fig, ax = plt.subplots()
xidx = search.keys.index('F0')
frequencies = np.unique([:, xidx])
twoF =[:, -1]
#mismatch = np.sign(x-F0)*(duration * np.pi * (x - F0))**2 / 12.0
ax.plot(frequencies, twoF, 'k', lw=1)
DeltaF = frequencies - F0
sinc = np.sin(np.pi*DeltaF*duration)/(np.pi*DeltaF*duration)
A = np.abs((np.max(twoF)-4)*sinc**2 + 4)
ax.plot(frequencies, A, '-r', lw=1)
ax.set_xlim(F0s[0], F0s[1])
dF0 = np.sqrt(12*1)/(np.pi*duration)
xticks = [F0-10*dF0, F0, F0+10*dF0]
xticklabels = ['$f_0 {-} 10\Delta f$', '$f_0$', '$f_0 {+} 10\Delta f$']
fig.savefig('{}/{}_1D.png'.format(search.outdir, search.label), dpi=300)
......@@ -80,8 +80,8 @@ theta_prior = {'F0': {'type': 'unif', 'lower': F0-DeltaF0/2.,
'Delta': Delta,
'tglitch': {'type': 'unif', 'lower': tstart+0.1*duration,
'upper': tend-0.1*duration},
'delta_F0': {'type': 'halfnorm', 'loc': 0, 'scale': 1e-3*F0},
'delta_F1': {'type': 'norm', 'loc': 0, 'scale': 1e-3*abs(F1)},
'delta_F0': {'type': 'halfnorm', 'loc': 0, 'scale': DeltaF0},
'delta_F1': {'type': 'norm', 'loc': 0, 'scale': DeltaF1},
ntemps = 3
log10temperature_min = -0.1
......@@ -94,6 +94,6 @@ glitch_mcmc = pyfstat.MCMCGlitchSearch(
nwalkers=nwalkers, ntemps=ntemps,
glitch_mcmc.plot_corner(figsize=(3.2, 3.2))
glitch_mcmc.plot_corner(figsize=(6, 6))
......@@ -87,5 +87,5 @@ mcmc = pyfstat.MCMCTransientSearch(
nwalkers=nwalkers, ntemps=ntemps,
......@@ -428,3 +428,39 @@ year = {2015}
year = {2012}
archivePrefix = {arXiv},
arxivId = {0710.1021},
author = {Melatos, A. and Peralta, C. and Wyithe, J. S. B.},
doi = {10.1086/523349},
eprint = {0710.1021},
file = {:home/greg/Dropbox/Papers/Mealtos{\_}2007.pdf:pdf},
issn = {0004-637X},
journal = {The Astrophysical Journal},
keywords = {dense matter — pulsars: general — stars: interiors},
number = {Jensen 1998},
pages = {1103},
title = {{Avalanche dynamics of radio pulsar glitches}},
url = {},
volume = {672},
year = {2008}
archivePrefix = {arXiv},
arxivId = {1102.1743},
author = {Espinoza, Crist{\'{o}}bal and Lyne, Andrew and Stappers, Ben and Kramer, Michael},
doi = {10.1063/1.3615093},
eprint = {1102.1743},
file = {:home/greg/Dropbox/Papers/Espinoza{\_}2011.pdf:pdf},
isbn = {9780735409156},
issn = {0094243X},
journal = {AIP Conference Proceedings},
keywords = {general,glitches,pulsars,stars: neutron},
number = {March},
pages = {117--120},
title = {{Glitches in the rotation of pulsars}},
volume = {1357},
year = {2011}
......@@ -3,6 +3,10 @@
\newcommand{\blambdaSignal}{\boldsymbol{\mathbf{\lambda}}^{\rm s}}
\newcommand{\tglitch}{t_{\rm glitch}}
\newcommand{\tstart}{t_{\rm start}}
\newcommand{\tend}{t_{\rm end}}
\newcommand{\Nglitches}{N_{\rm glitches}}
\newcommand{\Tspan}{T_{\rm span}}
\newcommand{\Tcoh}{T_{\rm coh}}
\newcommand{\tref}{t_{\rm ref}}

46.9 KB | W: | H:


42.1 KB | W: | H:

  • 2-up
  • Swipe
  • Onion skin
......@@ -296,6 +296,7 @@ computed analytically as
\mathcal{L}(x ;\A, \blambda)
P(\A| \Hs, I) d\A
& = \frac{C (2\pi)^{2} e^{\F(x| \blambda)}}
{\sqrt{\textrm{det} \mathcal{M}}},
......@@ -1110,14 +1111,84 @@ a simulated transient signal and Gaussian noise.}
Observations of radio pulsars show that occasionally neutron stars undergo
sudden glitch events during which the pulsation frequency suddenly
\emph{increases} quite distinct from the usual spin-down due to EM or GW
torques (see \citet{espinoza2011} for a review). These events are typically
Poisson distributed \citet{melatos2008}, but the exact mechanism is not yet
fully understood. However, it seems plausible that the glitch mechanism will
similarly effect any CW emission, i.e. CW signals will also undergo glitches.
In \citet{ashton2016} we demonstrated that if this is the case, CW sources in
directed and all-sky searches have a reasonable probability of undergoing one
or more glitches during typical observations spans and that these glitches may
be sufficiently large to disrupt the detection; this is particuarly true for
areas of parameter space with large spindowns. However, the effect is mitigated
by the use of incoherent searches. As such, we may expect all-sky and directed
searches to correctly identify glitching signals as candidates, but
subseuqnetly `loose' them during the follow-up process. As such, we will now
discuss how the $\F$-statistic optimisation routine discussed in this work can
be used to search for glitching signals.
We will model a glitch as a sudden discontinuous change in the Doppler
parameters $\delta \blambda$ at some time $\tglitch$: that the change is only
in the frequency and spin-down Doppler parameters, i.e. $\delta\blambda = [0,
0, \delta f, \delta \fdot, \ldots]$. This is to be interpretted, as is done for
radio-pulsar glitches, as the instantaneous change in frequency and spindown at
the time of the glitch. Moreover, multiple glitches can be included (we will
define $\Nglitches$ as the number of glitches) by adding
an index to the glitch magnitude and glitch time.
An ideal glitch search would include these new parameters into the phase model
itself and compute the $\F$-statistic fully-coherently i.e. rewriting
Equation~\eqref{eqn_B_S_N} as
\Bsn(x| \tglitch, \delta\blambda, \Pic, \blambda, I) \equiv \\
\mathcal{L}(x ;\A, \blambda, \tglitch, \delta\blambda)
P(\A| \Hs, I) d\A
where $\tglitch, \delta\blambda$ modify the phase model. However, we propose the
following pragmatic approach: let us instead use an semi-coherent detection
statistic with the epoch between \emph{segments} as $\tglitch$, i.e.
\Bsn(x| \tglitch, \delta\blambda, \Pic, \blambda, I) \equiv \\
\mathcal{L}(x ;\A, \blambda)
P(\A| \Hs, I) d\A \\ +
\mathcal{L}(x ;\A, \blambda{+}\delta\blambda)
P(\A| \Hs, I) d\A
This simplistic method leverages readily available and tested code with a loss
of sensitivity over the fully-coherent search by a factor
$\sqrt{\Nglitches+1}$. Such a loss is acceptable, given that the signals must
be sufficiently strong to have been identified by the initial all-sky or directed
As an illustration of the use of this method, we simulate a signal in Gaussian
noise which undergoes a glitch. Firstly in Figure~\ref{fig_glitch_grid} we show
the fully-coherent $2\F$ value computed for a grid of points: this distinctly
shows two peaks, a feature indicative of a glitching signal, with a maximum
value of $\sim\CHECK{400}$. Using our incoherent glitch statistic, we plot
the resulting posterior distributions in Figure~\ref{fig_glitch_posteriors}.

796 KB | W: | H:


628 KB | W: | H:

  • 2-up
  • Swipe
  • Onion skin
Supports Markdown
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