From 98327cb432b4564e8356dfcef82f1ef49243b640 Mon Sep 17 00:00:00 2001
From: "gregory.ashton" <gregory.ashton@ligo.org>
Date: Fri, 23 Dec 2016 10:06:33 +0100
Subject: [PATCH] Add initial analysis files for single glitch MC
---
Paper/Examples/single_glitch.py | 4 +-
.../SingleGlitchMCNoiseOnly.sh | 11 +++
.../SingleGlitchMCNoiseOnly/generate_data.py | 99 +++++++++++++++++++
Paper/SingleGlitchMCNoiseOnly/plot_data.py | 41 ++++++++
Paper/SingleGlitchMCNoiseOnly/submitfile | 12 +++
5 files changed, 165 insertions(+), 2 deletions(-)
create mode 100755 Paper/SingleGlitchMCNoiseOnly/SingleGlitchMCNoiseOnly.sh
create mode 100644 Paper/SingleGlitchMCNoiseOnly/generate_data.py
create mode 100644 Paper/SingleGlitchMCNoiseOnly/plot_data.py
create mode 100644 Paper/SingleGlitchMCNoiseOnly/submitfile
diff --git a/Paper/Examples/single_glitch.py b/Paper/Examples/single_glitch.py
index 61b1212..0abb7c4 100644
--- a/Paper/Examples/single_glitch.py
+++ b/Paper/Examples/single_glitch.py
@@ -12,8 +12,8 @@ tref = .5*(tstart + tend)
F0 = 30.0
F1 = -1e-10
F2 = 0
-Alpha = 5e-3
-Delta = 6e-2
+Alpha = np.radians(83.6292)
+Delta = np.radians(22.0144)
# Signal strength
depth = 10
diff --git a/Paper/SingleGlitchMCNoiseOnly/SingleGlitchMCNoiseOnly.sh b/Paper/SingleGlitchMCNoiseOnly/SingleGlitchMCNoiseOnly.sh
new file mode 100755
index 0000000..ab5264c
--- /dev/null
+++ b/Paper/SingleGlitchMCNoiseOnly/SingleGlitchMCNoiseOnly.sh
@@ -0,0 +1,11 @@
+#!/bin/bash
+
+. /home/gregory.ashton/lalsuite-install/etc/lalapps-user-env.sh
+export PATH="/home/gregory.ashton/anaconda2/bin:$PATH"
+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
+done
+cp /local/user/gregory.ashton/NoiseOnlyMCResults_"$1".txt $(pwd)/CollectedOutput
diff --git a/Paper/SingleGlitchMCNoiseOnly/generate_data.py b/Paper/SingleGlitchMCNoiseOnly/generate_data.py
new file mode 100644
index 0000000..7c2d8b4
--- /dev/null
+++ b/Paper/SingleGlitchMCNoiseOnly/generate_data.py
@@ -0,0 +1,99 @@
+import pyfstat
+import numpy as np
+import os
+import sys
+import time
+
+
+ID = sys.argv[1]
+outdir = sys.argv[2]
+
+label = 'run_{}'.format(ID)
+data_label = '{}_data'.format(label)
+results_file_name = '{}/NoiseOnlyMCResults_{}.txt'.format(outdir, ID)
+
+# Properties of the GW data
+sqrtSX = 1e-23
+tstart = 1000000000
+Tspan = 100*86400
+tend = tstart + Tspan
+
+# Fixed properties of the signal
+F0_center = 30
+F1_center = -1e-10
+F2 = 0
+Alpha = np.radians(83.6292)
+Delta = np.radians(22.0144)
+tref = .5*(tstart+tend)
+
+VF0 = VF1 = 200
+dF0 = np.sqrt(3)/(np.pi*Tspan)
+dF1 = np.sqrt(45/4.)/(np.pi*Tspan**2)
+DeltaF0 = VF0 * dF0
+DeltaF1 = VF1 * dF1
+
+nsteps = 25
+run_setup = [((nsteps, 0), 20, False),
+ ((nsteps, 0), 7, False),
+ ((nsteps, 0), 2, False),
+ ((nsteps, nsteps), 1, False)]
+
+h0 = 0
+F0 = F0_center + np.random.uniform(-0.5, 0.5)*DeltaF0
+F1 = F1_center + np.random.uniform(-0.5, 0.5)*DeltaF1
+
+psi = np.random.uniform(-np.pi/4, np.pi/4)
+phi = np.random.uniform(0, 2*np.pi)
+cosi = np.random.uniform(-1, 1)
+
+# Next, taking the same signal parameters, we include a glitch half way through
+dtglitch = Tspan/2.0
+delta_F0 = 0.25*DeltaF0
+delta_F1 = -0.1*DeltaF1
+
+glitch_data = pyfstat.Writer(
+ label=data_label, outdir=outdir, tref=tref, tstart=tstart, F0=F0,
+ F1=F1, F2=F2, duration=Tspan, Alpha=Alpha, Delta=Delta, h0=h0,
+ sqrtSX=sqrtSX, dtglitch=dtglitch, delta_F0=delta_F0, delta_F1=delta_F1)
+glitch_data.make_data()
+
+
+startTime = time.time()
+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*Tspan,
+ 'upper': tend-0.1*Tspan},
+ 'delta_F0': {'type': 'halfnorm', 'loc': 0, 'scale': DeltaF0},
+ 'delta_F1': {'type': 'norm', 'loc': 0, 'scale': DeltaF1},
+ }
+ntemps = 2
+log10temperature_min = -0.1
+nwalkers = 100
+nsteps = [500, 500]
+glitch_mcmc = pyfstat.MCMCGlitchSearch(
+ label=label, outdir=outdir,
+ sftfilepath='{}/*{}*sft'.format(outdir, data_label),
+ theta_prior=theta_prior,
+ tref=tref, minStartTime=tstart, maxStartTime=tend, nsteps=nsteps,
+ nwalkers=nwalkers, ntemps=ntemps,
+ log10temperature_min=log10temperature_min)
+glitch_mcmc.run(run_setup=run_setup, create_plots=False, log_table=False,
+ gen_tex_table=False)
+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']
+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))
+os.system('rm {}/*{}*'.format(outdir, label))
diff --git a/Paper/SingleGlitchMCNoiseOnly/plot_data.py b/Paper/SingleGlitchMCNoiseOnly/plot_data.py
new file mode 100644
index 0000000..dd2332a
--- /dev/null
+++ b/Paper/SingleGlitchMCNoiseOnly/plot_data.py
@@ -0,0 +1,41 @@
+import pyfstat
+import matplotlib.pyplot as plt
+import pandas as pd
+import numpy as np
+import os
+from tqdm import tqdm
+from oct2py import octave
+import glob
+from scipy.stats import rv_continuous, chi2
+
+filenames = glob.glob("CollectedOutput/*.txt")
+
+plt.style.use('paper')
+
+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'])
+ df['CLUSTER_ID'] = fn.split('_')[1]
+ df_list.append(df)
+df = pd.concat(df_list)
+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')
+
+ax.set_xlabel('$\widetilde{2\mathcal{F}}$')
+ax.set_xlim(0, 90)
+ax.legend(frameon=False, fontsize=6)
+fig.tight_layout()
+fig.savefig('single_glitch_noise_twoF_histogram.png')
+
+#from latex_macro_generator import write_to_macro
+#write_to_macro('DirectedMCNoiseOnlyMaximum', '{:1.1f}'.format(np.max(df.twoF)),
+# '../macros.tex')
+#write_to_macro('DirectedMCNoiseN', len(df), '../macros.tex')
diff --git a/Paper/SingleGlitchMCNoiseOnly/submitfile b/Paper/SingleGlitchMCNoiseOnly/submitfile
new file mode 100644
index 0000000..a11d97e
--- /dev/null
+++ b/Paper/SingleGlitchMCNoiseOnly/submitfile
@@ -0,0 +1,12 @@
+Executable= SingleGlitchMCNoiseOnly.sh
+Arguments= $(Cluster)_$(Process)
+Universe=vanilla
+Input=/dev/null
+accounting_group = ligo.dev.o2.cw.explore.test
+Output=CollectedOutput/out.$(Cluster).$(Process)
+Error=CollectedOutput/err.$(Cluster).$(Process)
+Log=CollectedOutput/log.$(Cluster).$(Process)
+request_cpus = 1
+request_memory = 16 GB
+
+Queue 100
--
GitLab