diff --git a/Paper/paper_cw_mcmc.tex b/Paper/paper_cw_mcmc.tex index 2c91873c1670b28e57ec714147d70b58930987de..c0ed37eac78a5607bb926407875988e518332093 100644 --- a/Paper/paper_cw_mcmc.tex +++ b/Paper/paper_cw_mcmc.tex @@ -792,7 +792,7 @@ We now provide an illustrative example of the follow-up method. We consider a directed search over the sky position and frequency in 100 days of data from a single detector, with $\sqrt{\Sn}=10^{-23}$~Hz$^{-1/2}$ (at the fiducial frequency of the signal). The simulated signal has an amplitude -$h_0=1.4\times10^{25}$ such that the signal has a depth of $\sqrt{\Sn}/h_0=70$ +$h_0=2\times10^{-25}$ such that the signal has a depth of $\sqrt{\Sn}/h_0=50$ in the noise. First, we must define the setup for the run. Using $\mathcal{R}=10$ and @@ -834,14 +834,28 @@ are listed in Table~\ref{tab_weak_signal_follow_up}.} \end{figure} The key advantage to note here is that all walkers succefully convereged to the -signal peak, which occupies $\approx 10^{-6}$ of the initial volume. While it -is possible for this to occur during an ordinary MCMC simulation (with $\Tcoh$ -fixed at $\Tspan$), it would take much longer to converge as the chains explore -the other `noise peaks' in the data. +signal peak, which occupies $\sim 10^{-6}$ of the initial volume. While it is +possible for this to occur during an ordinary MCMC simulation (with $\Tcoh$ +fixed at $\Tspan$), it would take substantially longer to converge as the +chains explore the other `noise peaks' in the data. \section{Alternative waveform models: transients} \label{sec_transients} +\begin{figure}[htb] +\centering +\includegraphics[width=0.5\textwidth]{transient_search_initial_stage_twoFcumulative} +\caption{} +\label{fig:} +\end{figure} + +\begin{figure}[htb] +\centering +\includegraphics[width=0.5\textwidth]{transient_search_corner} +\caption{} +\label{fig:} +\end{figure} + \section{Alternative waveform models: glitches} \label{sec_glitches} diff --git a/Paper/transient_search_corner.png b/Paper/transient_search_corner.png new file mode 100644 index 0000000000000000000000000000000000000000..7507dd26a823a3a5776dcd2db0d23e5458fb719a Binary files /dev/null and b/Paper/transient_search_corner.png differ diff --git a/Paper/transient_search_initial_stage_twoFcumulative.png b/Paper/transient_search_initial_stage_twoFcumulative.png new file mode 100644 index 0000000000000000000000000000000000000000..460e9da5b93382f17fe017898678aed340a8f658 Binary files /dev/null and b/Paper/transient_search_initial_stage_twoFcumulative.png differ diff --git a/Paper/weak_signal_follow_up_run_setup.tex b/Paper/weak_signal_follow_up_run_setup.tex index 6a77ddc95edd2997fc634a628bd434651bca561a..a43fb4a45d09ee85d526acbd63a7e59264aa5d92 100644 --- a/Paper/weak_signal_follow_up_run_setup.tex +++ b/Paper/weak_signal_follow_up_run_setup.tex @@ -1,9 +1,9 @@ \begin{tabular}{c|cccccc} Stage & $\Nseg$ & $\Tcoh^{\rm days}$ &$\Nsteps$ & $\V$ & $\Vsky$ & $\Vpe$ \\ \hline -0 & 93 & 1.1 & 100 & 10.0 & 1.0 & 10.0 \\ -1 & 43 & 2.3 & 100 & $1{\times}10^{2}$ & 6.0 & 20.0 \\ -2 & 20 & 5.0 & 100 & $1{\times}10^{3}$ & 30.0 & 50.0 \\ -3 & 9 & 11.1 & 100 & $1{\times}10^{4}$ & $1{\times}10^{2}$ & $1{\times}10^{2}$ \\ -4 & 4 & 25.0 & 100 & $1{\times}10^{5}$ & $6{\times}10^{2}$ & $2{\times}10^{2}$ \\ -5 & 1 & 100.0 & 100,100 & $1{\times}10^{6}$ & $1{\times}10^{3}$ & $9{\times}10^{2}$ \\ +0 & 93 & 1.1 & 100 & 20.0 & 2.0 & 10.0 \\ +1 & 43 & 2.3 & 100 & $2{\times}10^{2}$ & 10.0 & 20.0 \\ +2 & 20 & 5.0 & 100 & $2{\times}10^{3}$ & 50.0 & 50.0 \\ +3 & 9 & 11.1 & 100 & $2{\times}10^{4}$ & $2{\times}10^{2}$ & $1{\times}10^{2}$ \\ +4 & 4 & 25.0 & 100 & $2{\times}10^{5}$ & $1{\times}10^{3}$ & $2{\times}10^{2}$ \\ +5 & 1 & 100.0 & 100,100 & $2{\times}10^{6}$ & $3{\times}10^{3}$ & $9{\times}10^{2}$ \\ \end{tabular} diff --git a/Paper/weak_signal_follow_up_walkers.png b/Paper/weak_signal_follow_up_walkers.png index f389425ce0d15083b02cf310a0463e0db2a74cfd..d901ea5dce2cc0b572741a0c426c12d435ebe792 100644 Binary files a/Paper/weak_signal_follow_up_walkers.png and b/Paper/weak_signal_follow_up_walkers.png differ diff --git a/examples/make_fake_data.py b/examples/make_fake_data.py index 62c9d95cd54970ac278b148d63290ba36f23f5a3..e645dbd31facb55b5105a371513c372fe26db8ea 100644 --- a/examples/make_fake_data.py +++ b/examples/make_fake_data.py @@ -53,14 +53,3 @@ two_glitch_data = Writer( dtglitch=dtglitch, delta_phi=delta_phi, delta_F0=delta_F0, delta_F1=delta_F1, delta_F2=delta_F2) two_glitch_data.make_data() - - -# Making transient data in the middle third -data_tstart = tstart - duration -data_duration = 3 * duration -transient = 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, - data_tstart=data_tstart, data_duration=data_duration) -transient.make_data() - diff --git a/examples/transient_search_using_MCMC.py b/examples/transient_search_using_MCMC.py index 7fc7202bd39f2b3a0c23f64b784a335f1952fad2..5a0738a80544b69f468a7a125798c7ad308c56f2 100644 --- a/examples/transient_search_using_MCMC.py +++ b/examples/transient_search_using_MCMC.py @@ -1,35 +1,83 @@ -from pyfstat import MCMCTransientSearch +import pyfstat +import numpy as np F0 = 30.0 F1 = -1e-10 F2 = 0 Alpha = 5e-3 Delta = 6e-2 -tref = 362750407.0 tstart = 1000000000 duration = 100*86400 -tstart = 1000000000 - duration -tend = tstart + 3*duration +data_tstart = tstart - duration +data_tend = data_tstart + 3*duration +tref = .5*(data_tstart+data_tend) -theta_prior = {'F0': {'type': 'unif', 'lower': F0*(1-1e-6), 'upper': F0*(1+1e-6)}, - 'F1': {'type': 'unif', 'lower': F1*(1+1e-2), 'upper': F1*(1-1e-2)}, +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() + +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, - 'transient_tstart': {'type': 'unif', 'lower': tstart, 'upper': tend}, - 'transient_duration': {'type': 'halfnorm', 'loc':0, 'scale': duration} + 'Delta': Delta } -ntemps = 4 +ntemps = 3 log10temperature_min = -1 nwalkers = 100 -nsteps = [1000, 1000] +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() + +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 = MCMCTransientSearch( - label='transient_search_using_MCMC', outdir='data', +mcmc = pyfstat.MCMCTransientSearch( + label='transient_search', outdir='data', sftfilepath='data/*transient*sft', theta_prior=theta_prior, tref=tref, - tstart=tstart, tend=tend, nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps, + minStartTime=data_tstart, maxStartTime=data_tend, nsteps=nsteps, + nwalkers=nwalkers, ntemps=ntemps, log10temperature_min=log10temperature_min) mcmc.run() mcmc.plot_corner(add_prior=True) diff --git a/examples/weak_signal_follow_up.py b/examples/weak_signal_follow_up.py index 488fe1a402e2e534071e57a3ed8b2fab8ecc042b..38bb9937ff44e64a8650f2f1a191ce3b1ac94f0d 100644 --- a/examples/weak_signal_follow_up.py +++ b/examples/weak_signal_follow_up.py @@ -13,7 +13,7 @@ duration = 100*86400 tend = tstart+duration tref = .5*(tstart+tend) -depth = 70 +depth = 50 data_label = 'weak_signal_follow_up_depth_{:1.0f}'.format(depth) h0 = sqrtSX / depth @@ -41,8 +41,8 @@ theta_prior = {'F0': {'type': 'unif', 'lower': F0*(1-1e-6), } ntemps = 3 -log10temperature_min = -1 -nwalkers = 200 +log10temperature_min = -0.5 +nwalkers = 100 scatter_val = 1e-10 nsteps = [100, 100] @@ -52,7 +52,6 @@ mcmc = pyfstat.MCMCFollowUpSearch( minStartTime=tstart, maxStartTime=tend, nwalkers=nwalkers, nsteps=nsteps, ntemps=ntemps, log10temperature_min=log10temperature_min, scatter_val=scatter_val) -mcmc.run(R0=10, Vmin=100) +mcmc.run(R0=10, Vmin=100, subtractions=[F0, Alpha, Delta], context='paper') mcmc.plot_corner(add_prior=True) mcmc.print_summary() -#mcmc.generate_loudest() diff --git a/pyfstat.py b/pyfstat.py index 14561909ba43ae4de1962dd2e159b703835af02f..234f22a954db0098e9b35532fbfc3fc68259fc93 100755 --- a/pyfstat.py +++ b/pyfstat.py @@ -655,7 +655,8 @@ class ComputeFstat(object): else: ax.set_ylabel(r'$\widetilde{2\mathcal{F}}_{\rm cumulative}$') ax.set_xlim(0, taus[-1]/86400) - ax.set_title(title) + if title: + ax.set_title(title) if savefig: plt.savefig('{}/{}_twoFcumulative.png'.format(outdir, label)) return taus, twoFs