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

Adds new paper specific examples

parent a2c05b4d
import pyfstat
import numpy as np
import matplotlib.pyplot as plt
sqrtSX = 1e-22
tstart = 1000000000
duration = 100*86400
tend = tstart+duration
# Define parameters of the Crab pulsar as an example
tref = .5*(tstart + tend)
F0 = 30.0
F1 = -1e-10
F2 = 0
Alpha = 5e-3
Delta = 6e-2
# Signal strength
depth = 10
h0 = sqrtSX / depth
VF0 = VF1 = 200
dF0 = np.sqrt(3)/(np.pi*duration)
dF1 = np.sqrt(45/4.)/(np.pi*duration**2)
DeltaF0 = VF0 * dF0
DeltaF1 = VF1 * dF1
# Next, taking the same signal parameters, we include a glitch half way through
dtglitch = duration/2.0
delta_F0 = 0.25*DeltaF0
delta_F1 = -0.1*DeltaF1
glitch_data = pyfstat.Writer(
label='single_glitch', outdir='data', 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()
F0s = [F0-DeltaF0/2., F0+DeltaF0/2., 1*dF0]
F1s = [F1-DeltaF1/2., F1+DeltaF1/2., 1*dF1]
F2s = [F2]
Alphas = [Alpha]
Deltas = [Delta]
search = pyfstat.GridSearch(
'single_glitch_F0F1_grid', 'data', 'data/*single_glitch*sft', F0s, F1s,
F2s, Alphas, Deltas, tref, tstart, tend)
search.run()
search.plot_2D('F0', 'F1')
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 = 3
log10temperature_min = -0.05
nwalkers = 100
nsteps = [500, 500]
mcmc = pyfstat.MCMCSearch(
'single_glitch', 'data', sftfilepath='data/*_single_glitch*.sft',
theta_prior=theta_prior, tref=tref, minStartTime=tstart, maxStartTime=tend,
nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps,
log10temperature_min=log10temperature_min)
mcmc.run()
mcmc.plot_corner(figsize=(3.2, 3.2))
mcmc.print_summary()
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,
'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)},
}
ntemps = 3
log10temperature_min = -0.1
nwalkers = 100
nsteps = [1000, 1000]
glitch_mcmc = pyfstat.MCMCGlitchSearch(
'single_glitch_glitchSearch', 'data',
sftfilepath='data/*_single_glitch*.sft', theta_prior=theta_prior,
tref=tref, minStartTime=tstart, maxStartTime=tend, nsteps=nsteps,
nwalkers=nwalkers, ntemps=ntemps,
log10temperature_min=log10temperature_min)
glitch_mcmc.run()
glitch_mcmc.plot_corner(figsize=(3.2, 3.2))
glitch_mcmc.print_summary()
import pyfstat
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('thesis')
F0 = 30.0
F1 = -1e-10
F2 = 0
Alpha = 5e-3
Delta = 6e-2
tstart = 1000000000
duration = 100*86400
data_tstart = tstart - duration
data_tend = data_tstart + 3*duration
tref = .5*(data_tstart+data_tend)
h0 = 1e-23
sqrtSX = 1e-22
transient = pyfstat.Writer(
label='transient', outdir='data', tref=tref, tstart=tstart, F0=F0, F1=F1,
F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX,
minStartTime=data_tstart, maxStartTime=data_tend)
transient.make_data()
print transient.predict_fstat()
DeltaF0 = 6e-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 = 3
log10temperature_min = -1
nwalkers = 100
nsteps = [750, 250]
mcmc = pyfstat.MCMCSearch(
label='transient_search_initial_stage', outdir='data',
sftfilepath='data/*transient*sft', theta_prior=theta_prior, tref=tref,
minStartTime=data_tstart, maxStartTime=data_tend, nsteps=nsteps,
nwalkers=nwalkers, ntemps=ntemps,
log10temperature_min=log10temperature_min)
mcmc.run()
mcmc.plot_cumulative_max()
mcmc.print_summary()
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,
'transient_tstart': {'type': 'unif',
'lower': data_tstart,
'upper': data_tend},
'transient_duration': {'type': 'halfnorm',
'loc': 0,
'scale': 0.5*duration}
}
nwalkers = 500
nsteps = [200, 200]
mcmc = pyfstat.MCMCTransientSearch(
label='transient_search', outdir='data',
sftfilepath='data/*transient*sft', theta_prior=theta_prior, tref=tref,
minStartTime=data_tstart, maxStartTime=data_tend, nsteps=nsteps,
nwalkers=nwalkers, ntemps=ntemps,
log10temperature_min=log10temperature_min)
mcmc.run()
mcmc.plot_corner(add_prior=True)
mcmc.print_summary()
......@@ -1110,8 +1110,14 @@ a simulated transient signal and Gaussian noise.}
\subsection{Glitches}
\label{sec_glitches}
\label{sec_glitches}
\begin{figure}[htb]
\centering
\includegraphics[width=0.5\textwidth]{single_glitch_F0F1_grid_2D}
\caption{}
\label{fig:}
\end{figure}
\section{Conclusion}
\label{sec_conclusion}
......
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