diff --git a/Paper/Examples/single_glitch.py b/Paper/Examples/single_glitch.py new file mode 100644 index 0000000000000000000000000000000000000000..3d76a824acacdeee6c36bcb6537ce32dedcf10c5 --- /dev/null +++ b/Paper/Examples/single_glitch.py @@ -0,0 +1,99 @@ +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() + diff --git a/Paper/Examples/transient_search_using_MCMC.py b/Paper/Examples/transient_search_using_MCMC.py new file mode 100644 index 0000000000000000000000000000000000000000..c5c54c3c2fc318318e7b11c7b1c14025b6ded97f --- /dev/null +++ b/Paper/Examples/transient_search_using_MCMC.py @@ -0,0 +1,91 @@ +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() diff --git a/Paper/paper_cw_mcmc.tex b/Paper/paper_cw_mcmc.tex index 158f6cb27e89a3ac3358dec303fc76a7152fe3fa..83a18b854cf46a18e1978c62492b74da814672d3 100644 --- a/Paper/paper_cw_mcmc.tex +++ b/Paper/paper_cw_mcmc.tex @@ -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} diff --git a/Paper/single_glitch_F0F1_grid_2D.png b/Paper/single_glitch_F0F1_grid_2D.png new file mode 100644 index 0000000000000000000000000000000000000000..e13b851e3e86640cf2eef82b7ae93939f233211b Binary files /dev/null and b/Paper/single_glitch_F0F1_grid_2D.png differ