Enhancements/gw (#40)

* new gw-specific version of initialize_llh

* add help for gcn num flag

* allow call to get_data with start and stop times for gw

* write in tests, clean up

* globbing trials for gw

* corrections to requirements

* avoiding getting absurdly long reports in long tws

* adding options for 2 week gw followup

* script to run 2 week followup given xml

* changed label on p-val hist

* print exception if error importing FRA

* optons to change number of trials

* load ps sens, bg trials from /data/ana

* adding code to find contours in 2 week case

* keep np version fixed to avoid install errors

* starting code for CL calc

Author: jessiethw

fast_response/GWFollowup.py

@@ -8,6 +8,7 @@
 import matplotlib as mpl
 mpl.use('agg')
 import matplotlib.pyplot as plt
+import pickle
 
 from .FastResponseAnalysis import PriorFollowup
 from .reports import GravitationalWaveReport
@@ -47,17 +48,13 @@ def run_background_trials(self, month=None, ntrials=1000):
             from scipy import sparse
             current_rate = self.llh.nbackground / (self.duration * 86400.) * 1000.
 
-            #TODO: replace here once done
-            #closest_rate = sensitivity_utils.find_nearest(np.linspace(6.0, 7.2, 7), current_rate)
-            rates = [6.0,6.2,6.4,6.8,7.0]
-            closest_rate = sensitivity_utils.find_nearest(rates, current_rate)
+            closest_rate = sensitivity_utils.find_nearest(np.linspace(6.0, 7.2, 7), current_rate)
             print(f'Loading 2 week bg trials, rate: {closest_rate}')
 
-            #bg_trial_dir = '/data/ana/analyses/NuSources/' \
-            #    + '2023_realtime_gw_analysis/fast_response/' \
-            #    + 'precomputed_background/'
-            #TODO: change to permanent storage once assigned
-            bg_trial_dir = '/data/user/jthwaites/FastResponseAnalysis/output/trials/glob_trials/'
+            bg_trial_dir = '/data/ana/analyses/NuSources/' \
+                + '2023_realtime_gw_analysis/fast_response/' \
+                + 'precomputed_background/'
+
             pre_ts_array = sparse.load_npz(
                 bg_trial_dir
                 + 'gw_precomputed_trials_delta_t_'
@@ -188,7 +185,60 @@ def initialize_llh(self, skipped=None, scramble=False):
         llh.set_temporal_model(box)
 
         return llh
+    
+    def get_best_fit_contour(self, proportions=[0.5,0.9]):
+        '''Get a contour for the error around the best-fit point
+        proportions = confidence levels to calculate
+        -------
+        Makes a zoomed skymap of the ts-space with contours
+        '''
+        if self.tsd is None: return
+
+        from . import plotting_utils
+        #import meander
+        import seaborn as sns
+        
+        print('Calculating contour around best-fit TS')
+
+        #get threshold TS value for that level in the bg distribution
+        levels = [np.percentile(self.tsd, 100*(1-proportion)) for proportion in proportions]
+        sample_points = np.array(hp.pix2ang(self.nside, np.arange(len(self.skymap)))).T
+        loc=np.array((np.pi/2 - self.ts_scan['dec'], self.ts_scan['ra'])).T
+        contours_by_level = meander.spherical_contours(loc, self.ts_scan['TS_spatial_prior_0'], levels)
+        #print(contours_by_level)
+
+        thetas = []; phis=[]
+        for contours in contours_by_level:
+            for contour in contours:
+                theta, phi = contour.T
+                phi[phi<0] += 2.0*np.pi
+                thetas.append(theta)
+                phis.append(phi)
+
+        #norm_ts = self.ts_scan['TS_spatial_prior_0'] / sum(self.ts_scan['TS_spatial_prior_0'])
+        #thetas, phis = plotting_utils.plot_contours(proportions, norm_ts)
+        
+        #make the plot
+        pdf_palette = sns.color_palette("Blues", 500)
+        cmap = mpl.colors.ListedColormap(pdf_palette)
+
+        plotting_utils.plot_zoom(self.ts_scan['TS_spatial_prior_0'], self.skymap_fit_ra, self.skymap_fit_dec,
+                                 "", range = [0,10], reso=3., cmap = cmap)
+        
+        plotting_utils.plot_color_bar(range=[0,6], cmap=cmap, col_label=r"TS",offset=-50)
+        cont_ls = ['solid', 'dashed']*(int(len(proportions)/2))
+        cont_labels=[f'{proportion*100:.0f}/% CL' for proportion in proportions]
+
+        for i in range(len(thetas)):
+            hp.projplot(thetas[i], phis[i], linewidth=2., c='k')#, ls=cont_ls[i], label=cont_labels[i])
 
+        plt.scatter(0,0, marker='*', c = 'k', s = 130, label = "Scan Hot Spot") 
+        plt.legend(loc = 2, ncol=1, fontsize = 16, framealpha = 0.95)
+
+        plt.savefig(self.analysispath + '/' + self.analysisid + 'ts_contours.png',bbox_inches='tight')
+        plt.savefig(self.analysispath + '/' + self.analysisid + 'ts_contours.pdf',bbox_inches='tight', dpi=300)
+        plt.close()
+        
     def plot_ontime(self, with_contour=True, contour_files=None):
         return super().plot_ontime(with_contour=True, contour_files=contour_files)
 
@@ -410,15 +460,22 @@ def ps_sens_range(self):
         --------
         low: float
             lowest sensitivity within dec range
-        high: floaot
+        high: float
             highest sensitivity wihtin dec range
         '''
-        dec_range = np.linspace(-85,85,35)
-        sens = [1.15, 1.06, .997, .917, .867, .802, .745, .662,
-                .629, .573, .481, .403, .332, .250, .183, .101,
-                .035, .0286, .0311, .0341, .0361, .0394, .0418,
-                .0439, .0459, .0499, .0520, .0553, .0567, .0632,
-                .0679, .0732, .0788, .083, .0866]
+        sens_dir = '/data/ana/analyses/NuSources/2023_realtime_gw_analysis/' \
+                +  'fast_response/ps_sensitivities'
+
+        with open(f'{sens_dir}/ps_sensitivities_deltaT_{self.duration*86400.:.2e}s.pkl','rb') as f:
+            saved_sens=pickle.load(f)
+            dec_range=saved_sens['dec']
+            sens=saved_sens['sens_flux']
+        #dec_range = np.linspace(-85,85,35)
+        #sens = [1.15, 1.06, .997, .917, .867, .802, .745, .662,
+        #        .629, .573, .481, .403, .332, .250, .183, .101,
+        #        .035, .0286, .0311, .0341, .0361, .0394, .0418,
+        #        .0439, .0459, .0499, .0520, .0553, .0567, .0632,
+        #        .0679, .0732, .0788, .083, .0866]
 
         src_theta, src_phi = hp.pix2ang(self.nside, self.ipix_90)
         src_dec = np.pi/2. - src_theta
@@ -477,13 +534,19 @@ def make_dec_pdf(self):
 
         sinDec_bins = np.linspace(-1,1,30)
         bin_centers = (sinDec_bins[:-1] + sinDec_bins[1:]) / 2
-
-        dec_range = np.linspace(-1,1,35)
-        sens = [1.15, 1.06, .997, .917, .867, .802, .745, .662,
-                .629, .573, .481, .403, .332, .250, .183, .101,
-                .035, .0286, .0311, .0341, .0361, .0394, .0418,
-                .0439, .0459, .0499, .0520, .0553, .0567, .0632,
-                .0679, .0732, .0788, .083, .0866]
+        sens_dir = '/data/ana/analyses/NuSources/2023_realtime_gw_analysis/' \
+                +  'fast_response/ps_sensitivities'
+        with open(f'{sens_dir}/ps_sensitivities_deltaT_{self.duration*86400.:.2e}s.pkl','rb') as f:
+            saved_sens=pickle.load(f)
+            dec_range=np.sin(saved_sens['dec']*np.pi/180)
+            sens=saved_sens['sens_flux']
+
+        #dec_range = np.linspace(-1,1,35)
+        #sens = [1.15, 1.06, .997, .917, .867, .802, .745, .662,
+        #        .629, .573, .481, .403, .332, .250, .183, .101,
+        #        .035, .0286, .0311, .0341, .0361, .0394, .0418,
+        #        .0439, .0459, .0499, .0520, .0553, .0567, .0632,
+        #        .0679, .0732, .0788, .083, .0866]
         sens = np.array(sens)
 
         pixels = np.arange(len(self.skymap))
@@ -530,5 +593,8 @@ def generate_report(self):
         r'''Generates report using class attributes
         and the ReportGenerator Class'''
         report = GravitationalWaveReport(self)
+        if self.duration > 1.:
+            report._figure_list = [('decPDF', 'decPDF.png'),('ts_contours', 'ts_contours.png')]
+
         report.generate_report()
         report.make_pdf()

fast_response/listeners/gw_gcn_listener.py

@@ -53,14 +53,17 @@ def process_gcn(payload, root):
               for elem in root.iterfind('.//Param')}
 
     # Read trigger time of event
-    if root.attrib['role'] == 'observation':
-        eventtime = root.find('.//ISOTime').text
-        event_mjd = Time(eventtime, format='isot').mjd
-    else:
+    eventtime = root.find('.//ISOTime').text
+    event_mjd = Time(eventtime, format='isot').mjd
+    
+    if root.attrib['role'] == 'test':
         #if testing, want to query livestream rather than load archival, so use recent time
         eventtime = '2023-01-13T21:51:25.506'
         event_mjd = Time(eventtime, format='isot').mjd - 400./86400.
         eventtime = Time(event_mjd, format = 'mjd').isot
+    else:
+        eventtime = root.find('.//ISOTime').text
+        event_mjd = Time(eventtime, format='isot').mjd
 
     print('GW merger time: %s \n' % Time(eventtime, format='isot').iso)
     log_file.flush()
@@ -179,7 +182,8 @@ def process_gcn(payload, root):
                         help='bool to decide if we should run an already unblinded skymap with unblinded data')
     args = parser.parse_args()
 
-    logfile=args.log_path +'log.log'
+    logfile=args.log_path +'/log.log'
+    print(f'Logging to file: {logfile}')
     original_stdout=sys.stdout
     log_file = open(logfile, "a+")
     sys.stdout=log_file
@@ -198,6 +202,7 @@ def process_gcn(payload, root):
             #sample_skymap_path='/data/user/jthwaites/o3-gw-skymaps/'
             sample_skymap_path=os.path.dirname(fast_response.__file__) +'/sample_skymaps/'
         except Exception as e:
+            print(e)
             sample_skymap_path='/data/user/jthwaites/o3-gw-skymaps/'
 
         payload = open(sample_skymap_path+args.test_path, 'rb').read()

fast_response/precomputed_background/glob_precomputed_trials.py

@@ -50,7 +50,7 @@ def glob_allsky_scans(delta_t, rate, dir, low_stats=False):
 
     return maps
 
-for rate in [6.0, 6.2, 6.4, 6.8, 7.0]: # 6.6, 7.2
+for rate in [6.0, 6.2, 6.4, 6.6, 6.8, 7.0, 7.2]:
     for low_stats in [True]:#, False]:
         print("Rate: {} mHz, low stats: {}".format(rate, low_stats))
         maps = glob_allsky_scans(args.deltaT, rate, args.outdir, low_stats=low_stats)

fast_response/precomputed_background/submit_precomputed_trials.py

@@ -14,6 +14,12 @@
 parser.add_argument(
     '--ntrials', type=int, default=30000,
     help='Number of precomputed trials to run, default 30,000')
+parser.add_argument(
+    '--seed_start', type=int,default=0,
+    help='Seed to start with when running trials')
+parser.add_argument(
+    '--n_per_batch', default=100, type=int,
+    help='Number of trials to run in each set (default: 100)')
 args = parser.parse_args()
 
 username = pwd.getpwuid(os.getuid())[0]
@@ -65,13 +71,14 @@
 #    )
 
 prev_trials = glob.glob('/data/user/jthwaites/FastResponseAnalysis/output/trials/*.npz')
-for bg_rate in [7.2]:#[6.0, 6.2, 6.4, 6.6, 6.8, 7.0, 7.2]:
-    for seed in range(int(args.ntrials/100)):
+for bg_rate in [6.6]:#[6.0, 6.2, 6.4, 6.6, 6.8, 7.0, 7.2]:
+    for seed in range(args.seed_start, int(args.ntrials/args.n_per_batch)):
         seed = seed*100
         if f'/data/user/jthwaites/FastResponseAnalysis/output/trials/gw_{bg_rate}_mHz_seed_{seed}_delta_t_1.2e+06.npz' not in prev_trials:
             #deltaT {args.tw} --ntrials 100 --seed {seed} --bkg {bg}
             job.add_arg('--deltaT %s --ntrials %i --seed %i --bkg %s'
-                    %(args.tw, 100, seed, bg_rate))
+                    %(args.tw, args.n_per_batch, seed, bg_rate))
+
 
 #job.add_child(glob_jobs)
 

fast_response/scripts/run_external_followup.py

@@ -44,7 +44,7 @@ def run_analysis(args):
     f.calc_pvalue(ntrials = args.ntrials)
     f.make_dNdE()
     f.plot_tsd()
-    f.upper_limit()
+    f.upper_limit(n_per_sig = args.n_per_sig)
     results = f.save_results()
 
     print("Beginning to make the report")
@@ -80,6 +80,8 @@ def run_analysis(args):
                         "Example --skip-events=127853:67093193")
     parser.add_argument('--ntrials', default=100, type=int,
                         help="Number of background trials to perform")
+    parser.add_argument('--n_per_sig', default=100, type=int,
+                        help="Number of signal trials per injected ns to perform")
     parser.add_argument('--document', default=False, action='store_true')
     args = parser.parse_args()
 

fast_response/scripts/run_gw_followup.py

@@ -41,6 +41,7 @@
     print('Beginning 2 week NS followup')
     start_time = gw_time - 0.1
     stop_time = gw_time + 14.
+
 start = start_time.iso
 stop = stop_time.iso
 
@@ -56,6 +57,9 @@
 f.calc_pvalue()
 f.make_dNdE()
 f.plot_tsd(allow_neg=f._allow_neg)
+#if delta_t/86400. > 1.:
+#    f.get_best_fit_contour()
+
 f.upper_limit()
 f.find_coincident_events()
 f.per_event_pvalue()

fast_response/web_utils.py

@@ -253,7 +253,7 @@ def updateFastResponsePlots(gw=False):
     p_x_vals = np.logspace(-3,0.,15)
     plt.figure(figsize = (10,6), dpi=300)
     plt.hist(df['Pre-trial p_val'], weights = np.ones(len(df)) / len(df), bins = p_x_vals)
-    plt.step(p_x_vals[1:], np.diff(p_x_vals), label = 'Uniform p-value expectation', lw = 3.)
+    plt.step(p_x_vals[1:], np.diff(p_x_vals), label = 'Uniform p-value distribution', lw = 3.)
     plt.xscale('log')
     plt.yscale('log')
     plt.gca().invert_xaxis()

requirements.txt

@@ -11,7 +11,7 @@ kiwisolver==1.1.0
 matplotlib==2.2.5
 meander==0.0.3
 more-itertools==5.0.0
-numpy>=1.21.5
+numpy==1.21.5
 numexpr==2.8.1
 pandas>=1.2.2
 pluggy==0.7.1
@@ -28,4 +28,4 @@ six==1.14.0
 subprocess32==3.5.4
 zmq==0.0.0
 py27hash==1.0.2
-pygcn>=1.1.2
+pygcn==1.1.2