Generalises noise only glitch MCMC to Nglitches

Paper/SingleGlitchMCNoiseOnly/SingleGlitchMCNoiseOnly.sh → Paper/GlitchMCNoiseOnly/SingleGlitchMCNoiseOnly.sh

@@ -6,6 +6,6 @@ export MPLCONFIGDIR=/home/gregory.ashton/.config/matplotlib
 
 for ((n=0;n<90;n++))
 do
-/home/gregory.ashton/anaconda2/bin/python generate_data.py "$1" /local/user/gregory.ashton --no-template-counting --no-interactive
+/home/gregory.ashton/anaconda2/bin/python generate_data.py "$1" /local/user/gregory.ashton --no-interactive
 done
 cp /local/user/gregory.ashton/NoiseOnlyMCResults_"$1".txt $(pwd)/CollectedOutput

Paper/SingleGlitchMCNoiseOnly/generate_data.py → Paper/GlitchMCNoiseOnly/generate_data.py

@@ -4,6 +4,7 @@ import os
 import sys
 import time
 
+nglitch = 2
 
 ID = sys.argv[1]
 outdir = sys.argv[2]
@@ -80,7 +81,7 @@ glitch_mcmc = pyfstat.MCMCGlitchSearch(
     sftfilepath='{}/*{}*sft'.format(outdir, data_label),
     theta_prior=theta_prior,
     tref=tref, minStartTime=tstart, maxStartTime=tend, nsteps=nsteps,
-    nwalkers=nwalkers, ntemps=ntemps,
+    nwalkers=nwalkers, ntemps=ntemps, nglitch=nglitch,
     log10temperature_min=log10temperature_min)
 glitch_mcmc.run(run_setup=run_setup, create_plots=False, log_table=False,
                 gen_tex_table=False)
@@ -88,12 +89,12 @@ glitch_mcmc.print_summary()
 d, maxtwoF = glitch_mcmc.get_max_twoF()
 dF0 = F0 - d['F0']
 dF1 = F1 - d['F1']
-tglitch = d['tglitch']
-R = (tglitch - tstart) / Tspan
-delta_F0 = d['delta_F0']
-delta_F1 = d['delta_F1']
+#tglitch = d['tglitch']
+#R = (tglitch - tstart) / Tspan
+#delta_F0 = d['delta_F0']
+#delta_F1 = d['delta_F1']
 runTime = time.time() - startTime
 with open(results_file_name, 'a') as f:
-    f.write('{:1.8e} {:1.8e} {} {:1.8e} {:1.8e} {:1.8e} {}\n'
-            .format(dF0, dF1, R, delta_F0, delta_F1, maxtwoF, runTime))
+    f.write('{} {:1.8e} {:1.8e} {:1.8e} {:1.1f}\n'
+            .format(nglitch, dF0, dF1, maxtwoF, runTime))
 os.system('rm {}/*{}*'.format(outdir, label))

Paper/SingleGlitchMCNoiseOnly/plot_data.py → Paper/GlitchMCNoiseOnly/plot_data.py

@@ -17,8 +17,7 @@ Tspan = 100 * 86400
 df_list = []
 for fn in filenames:
     df = pd.read_csv(
-        fn, sep=' ', names=['dF0', 'dF1', 'R', 'delta_F0', 'delta_F1',
-                            'twoF', 'runTime'])
+        fn, sep=' ', names=['nglitches', 'dF0', 'dF1', 'twoF', 'runTime'])
     df['CLUSTER_ID'] = fn.split('_')[1]
     df_list.append(df)
 df = pd.concat(df_list)
@@ -26,8 +25,9 @@ print 'Number of samples = ', len(df)
 print 'Max twoF', df.twoF.max()
 
 fig, ax = plt.subplots()
-ax.hist(df.twoF, bins=50, histtype='step', color='k', normed=True, linewidth=1,
-        label='Monte-Carlo histogram')
+for ng in np.unique(df.nglitches.values):
+    ax.hist(df[df.nglitches==ng].twoF, bins=20, histtype='step', normed=True,
+            linewidth=1, label='$N_\mathrm{{glitches}}={}$'.format(ng))
 
 ax.set_xlabel('$\widehat{2\mathcal{F}}$ for 1 glitch')
 ax.set_xlim(0, 90)