Enhancements/gw (#42)

* resolve #41

* fixes to levels for contour

* clean up script

fast_response/GWFollowup.py

@@ -193,19 +193,22 @@ def get_best_fit_contour(self, proportions=[0.5,0.9]):
         Makes a zoomed skymap of the ts-space with contours
         '''
         if self.tsd is None: return
-
+
+        from scipy import special
         from . import plotting_utils
-        #import meander
+        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
+        dof = 2
+        delta_llh_levels = special.gammaincinv(dof/2.0, np.array([(1-prop) for prop in proportions]))
+        levels=2*delta_llh_levels
+
+        #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:
@@ -225,7 +228,7 @@ def get_best_fit_contour(self, proportions=[0.5,0.9]):
         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)
+        plotting_utils.plot_color_bar(labels=[0,round(max(self.ts_scan['TS_spatial_prior_0']))], cmap=cmap, col_label=r"TS",offset=-150)
         cont_ls = ['solid', 'dashed']*(int(len(proportions)/2))
         cont_labels=[f'{proportion*100:.0f}/% CL' for proportion in proportions]
 

fast_response/precomputed_background/submit_precomputed_trials.py

@@ -18,8 +18,8 @@
     '--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)')
+    '--n_per_batch', default=50, type=int,
+    help='Number of trials to run in each set (default: 50)')
 args = parser.parse_args()
 
 username = pwd.getpwuid(os.getuid())[0]
@@ -53,41 +53,16 @@
         'when_to_transfer_output = ON_EXIT']
     )
 
-#glob_jobs = pycondor.Job(
-#    'glob_gw_precomp',
-#    './glob_precomputed_trials.py',
-#    error=error,
-#    output=output,
-#    log=log,
-#    submit=submit,
-#    getenv=True,
-#    universe='vanilla',
-#    verbose=2,
-#    request_cpus=5,
-#    request_memory=5000,
-#    extra_lines=[
-#        'should_transfer_files = YES',
-#        'when_to_transfer_output = ON_EXIT']
-#    )
-
 prev_trials = glob.glob('/data/user/jthwaites/FastResponseAnalysis/output/trials/*.npz')
-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)):
+for bg_rate in [6.0, 6.2, 6.4, 6.6, 6.8, 7.0, 7.2]:
+    for seed in range(args.seed_start, int(args.ntrials/100)):#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, args.n_per_batch, seed, bg_rate))
-
-
-#job.add_child(glob_jobs)
+            job.add_arg('--deltaT %s --ntrials %i --seed %i --bkg %s'
+                    %(args.tw, args.n_per_batch, seed+50, bg_rate))
 
 job.build_submit()
 
-#should really use a dag - not working and haven't figured out why
-#dagman = pycondor.Dagman(
-#    'gw_precomp_trials',
-#    submit=submit, verbose=2)
-
-#dagman.add_job(job)
-#dagman.build_submit()

fast_response/scripts/run_gw_followup.py

@@ -34,6 +34,7 @@
 
 gw_time = Time(args.time, format='mjd')
 delta_t = float(args.tw)
+
 if args.tw==1000:
     start_time = gw_time - (delta_t / 86400. / 2.)
     stop_time = gw_time + (delta_t / 86400. / 2.)

requirements.txt

@@ -29,3 +29,4 @@ subprocess32==3.5.4
 zmq==0.0.0
 py27hash==1.0.2
 pygcn==1.1.2
+wheel==0.37.1