From fceed12713ce7aa77ad7bcde9fcb1f82c46e4a63 Mon Sep 17 00:00:00 2001 From: Sam Hori Date: Tue, 30 Jul 2024 15:04:52 -0500 Subject: [PATCH 1/7] added rough posterior map first draft code --- fast_response/FastResponseAnalysis.py | 90 ++++++++++++++++++++++++++- 1 file changed, 89 insertions(+), 1 deletion(-) diff --git a/fast_response/FastResponseAnalysis.py b/fast_response/FastResponseAnalysis.py index 0cb7b3d..fbb6451 100644 --- a/fast_response/FastResponseAnalysis.py +++ b/fast_response/FastResponseAnalysis.py @@ -15,7 +15,7 @@ import matplotlib as mpl import matplotlib.pyplot as plt from astropy.time import Time -from scipy.special import erfinv +from scipy.special import erfinv, logsumexp from matplotlib.lines import Line2D from skylab.datasets import Datasets @@ -992,6 +992,94 @@ def unblind_TS(self, custom_events=None): return ts, ns, gamma_fit else: return ts, ns + + + + def make_posterior_map(self, nsToTest=[1,2,3,4,5], prior=lambda ns,gamma,ra,dec: 1, testingPlotsLocation=None): + # Posterior Approach + # We want a function P(RA, DEC | X) where X is the neutrino data + # Compute the entire likelihood space P(X | RA, DEC, ns, gamma) + # Apply Bayes rule, marginalize over nuisance parameters to get P(RA, DEC | X) + + #ns to test give as list + #Give prior as a function f(ns, gamma, ra, dec) + #TODO: flux/ns priors + val=None + spatial_prior = SpatialPrior( + self.skymap, + containment = self._containment,allow_neg=self._allow_neg) + for nsigT in nsToTest: + temp_val = self.llh.scan(0.0,0.0, scramble = False,spatial_prior=spatial_prior, + time_mask = [self.duration/2., self.centertime], + pixel_scan=[self.nside, self._pixel_scan_nsigma], + fixed=['nsignal'], + nsignal=np.array([nsigT]*12*pow(self.nside,2),dtype=float) #TODO: Check if we can do more than 1 at a time + ) + if val is None: + val=temp_val + else: + val=np.hstack((val,temp_val)) + ras=val["ra"] + decs=val["dec"] + nss=val["nsignal"] + gammas=val["index"] + TS_spatial_prior_0=val["TS_spatial_prior_0"] + + logProbs=TS_spatial_prior_0/2 #not normalized so technically C*ln(P) + finalProbs={} + for i in range(len(ras)): + prior=1 #TODO + if((ras[i],decs[i]) in finalProbs): + finalProbs[(ras[i],decs[i])]=logsumexp([finalProbs[(ras[i],decs[i])],logProbs[i]+np.log(prior)]) + else: + finalProbs[(ras[i],decs[i])]=logProbs[i]+np.log(prior) + #Normalize + totalLog=logsumexp(list(finalProbs.values())) + finalProbs={k:v-totalLog for k,v in finalProbs.items()} + outputProbHPMap=np.zeros(12*self.nside**2) + outputProbHPMap[hp.ang2pix(nside=self.nside,theta=np.array(ras_graph)*180/np.pi, phi=np.array(decs_graph)*180/np.pi,lonlat=True)]=np.exp(np.array(list(finalProbs.values()))) + + #set 0s to nan for plotting + outputProbHPMap[outputProbHPMap == 0] = np.nan + + if (testingPlotsLocation is not None): + fig, ax=plt.subplots(1,1,figsize=(10,10)) + ras_graph,decs_graph=zip(*list(finalProbs.keys())) + fig, ax=plt.subplots(1,1,figsize=(10,10)) + plt.sca(ax) + reso=3 + nPixZoom=200 + sizeZoomInDegs=reso/60*nPixZoom + + """ + hp.visufunc.gnomview(map=outputProbHPMap,hold=True,title="Zoomed Posterior Skymap",rot=(src_ra*180/np.pi,src_dec*180/np.pi,0), + xsize=nPixZoom,ysize=nPixZoom,reso=reso,bgcolor="white",badcolor="white",margins=(.01,.01,0,0),notext=True,) + """ + hp.visufunc.mollview(map=outputProbHPMap,hold=True,title="Posterior Skymap",bgcolor="white",badcolor="white",rot=(180,0,0)) + plt.savefig(testingPlotsLocation+"PosteriorMap.png") + plt.close() + + + for (ra,dec) in finalProbs: + total_val=finalProbs[(ra,dec)] + temp_val=val[(val["ra"]==ra)&(val["dec"]==dec)] + nss_temp=temp_val["nsignal"] + gammas_temp=temp_val["index"] + TS_temp=temp_val["TS_spatial_prior_0"]/2 + + totalLog=logsumexp(TS_temp) + p_temp=np.exp(TS_temp-totalLog) + plt.plot(nss_temp,p_temp,alpha=total_val) + plt.colorbar() + plt.xlabel("ns") + plt.ylabel("TS") + plt.savefig(testingPlotsLocation+"nsProfilePlot.png") + plt.close() + + return outputProbHPMap + + + def upper_limit(self): """ UPPER LIMIT WITH SPATIAL PRIOR NOT YET IMPLEMENTED From 1ddbf1932a460ce0b4640ec8fcd3ac546e49abb6 Mon Sep 17 00:00:00 2001 From: Sam Hori Date: Wed, 8 Jan 2025 12:30:59 -0600 Subject: [PATCH 2/7] Rough draft of posterior skymap code --- fast_response/scripts/test_posterior_bash.sh | 32 +++++++ fast_response/scripts/test_posterior_bash2.sh | 15 +++ .../scripts/test_posterior_cascade.py | 90 ++++++++++++++++++ fast_response/scripts/test_posterior_code.py | 92 +++++++++++++++++++ 4 files changed, 229 insertions(+) create mode 100644 fast_response/scripts/test_posterior_bash.sh create mode 100644 fast_response/scripts/test_posterior_bash2.sh create mode 100755 fast_response/scripts/test_posterior_cascade.py create mode 100644 fast_response/scripts/test_posterior_code.py diff --git a/fast_response/scripts/test_posterior_bash.sh b/fast_response/scripts/test_posterior_bash.sh new file mode 100644 index 0000000..92c5743 --- /dev/null +++ b/fast_response/scripts/test_posterior_bash.sh @@ -0,0 +1,32 @@ +#python test_posterior_code.py --skymap=https://gracedb.ligo.org/apiweb/superevents/S240716b/files/Bilby.fits.gz --time=60000.00 --name="PosteriorTest_RealEventRandomDate_1000s_2" --tw=1000 +#https://gracedb.ligo.org/superevents/S240716b/files/ +#astropy.time.Time(1384415586.09,format='gps').mjd + +#S231119u-4-Update +python test_posterior_code.py --skymap=https://gracedb.ligo.org/api/superevents/S231119u/files/bayestar.fits.gz \ + --time=60267.328334432874 --name="S231119u_Test_Posterior" --tw=1000 +#S230928cb-4-Update +python test_posterior_code.py --skymap=https://gracedb.ligo.org/api/superevents/S230928cb/files/bayestar.fits.gz \ + --time=60215.915591805555 --name="S230928cb_Test_Posterior" --tw=1000 +#S230518h-4-Update +python test_posterior_code.py --skymap=https://gracedb.ligo.org/api/superevents/S230518h/files/bayestar.fits.gz \ + --time=60082.54106674768 --name="S230518h_Test_Posterior" --tw=1000 +#S230522a-2-Preliminary +python test_posterior_code.py --skymap=https://gracedb.ligo.org/api/superevents/S230522a/files/bayestar.fits.gz \ + --time=60086.40144832176 --name="S230522a_Test_Posterior" --tw=1000 + +#S230529ay-4-Update +python test_posterior_code.py --skymap=https://gracedb.ligo.org/api/superevents/S230529ay/files/bayestar.fits.gz \ + --time=60093.76042530093 --name="S230529ay_Test_Posterior" --tw=1000 +#S240615dg-4-Update +python test_posterior_code.py --skymap=https://gracedb.ligo.org/api/superevents/S240615dg/files/bayestar.fits.gz \ + --time=60476.48357319445 --name="S240615dg_Test_Posterior" --tw=1000 +#S231029y-4-Update +python test_posterior_code.py --skymap=https://gracedb.ligo.org/api/superevents/S231029y/files/bayestar.fits.gz \ + --time=60246.46885131944 --name="S231029y_Test_Posterior" --tw=1000 +#S230904n-4-Update +python test_posterior_code.py --skymap=https://gracedb.ligo.org/api/superevents/S230904n/files/bayestar.fits.gz \ + --time=60191.21542972222 --name="S230904n_Test_Posterior" --tw=1000 +#MS240802n-3-Initial_test +python test_posterior_code.py --skymap=https://gracedb.ligo.org/api/superevents/MS240802n/files/bayestar.fits.gz \ + --time=60524.544673460645 --name="MS240802n_Test_Posterior" --tw=1000 diff --git a/fast_response/scripts/test_posterior_bash2.sh b/fast_response/scripts/test_posterior_bash2.sh new file mode 100644 index 0000000..d7a35d6 --- /dev/null +++ b/fast_response/scripts/test_posterior_bash2.sh @@ -0,0 +1,15 @@ +#python test_posterior_code.py --skymap=https://gracedb.ligo.org/apiweb/superevents/S240716b/files/Bilby.fits.gz --time=60000.00 --name="PosteriorTest_RealEventRandomDate_1000s_2" --tw=1000 +#https://gracedb.ligo.org/superevents/S240716b/files/ +#astropy.time.Time(1384415586.09,format='gps').mjd + +python test_posterior_code.py --skymap=https://gracedb.ligo.org/apiweb/superevents/S240716b/files/Bilby.fits.gz --time=60000.00 \ + --name="Posterior_Test_0inj" --tw=1000 --n_inj=0 + +python test_posterior_code.py --skymap=https://gracedb.ligo.org/apiweb/superevents/S240716b/files/Bilby.fits.gz --time=60000.00 \ + --name="Posterior_Test_1inj" --tw=1000 --n_inj=1 + +python test_posterior_code.py --skymap=https://gracedb.ligo.org/apiweb/superevents/S240716b/files/Bilby.fits.gz --time=60000.00 \ + --name="Posterior_Test_2inj" --tw=1000 --n_inj=2 + +python test_posterior_code.py --skymap=https://gracedb.ligo.org/apiweb/superevents/S240716b/files/Bilby.fits.gz --time=60000.00 \ + --name="Posterior_Test_3inj" --tw=1000 --n_inj=3 diff --git a/fast_response/scripts/test_posterior_cascade.py b/fast_response/scripts/test_posterior_cascade.py new file mode 100755 index 0000000..7bdd845 --- /dev/null +++ b/fast_response/scripts/test_posterior_cascade.py @@ -0,0 +1,90 @@ +#!/usr/bin/env python + +r'''Script to run followup to +a realtime cascade alert + +Author: Alex Pizzuto +May 2020''' + +import numpy as np +import os, sys, argparse, subprocess +from astropy.time import Time + +from fast_response.AlertFollowup import CascadeFollowup +import fast_response.web_utils as web_utils + +parser = argparse.ArgumentParser(description='Fast Response Analysis') +parser.add_argument('--skymap', type=str, default=None, + help='path to skymap') +parser.add_argument('--time', type=float, default=None, + help='Time of the alert event (mjd)') +parser.add_argument('--gcn_notice_num', default=0, type=int, + help="Number of GCN circular. If not set, links to automated notice") +parser.add_argument('--alert_id', default=None, + type= lambda z:[ tuple(int(y) for y in x.split(':')) for x in z.split(',')], + help="list of events to exclude from this analysis. " + "such as HESE events that contributed to the trigger." + "Example --alert_id 127853:67093193,128290:6888376") +parser.add_argument('--suffix', type=str, default='A', + help="letter to differentiate multiple alerts on the same day (default = A)." + "Event name given by IceCube-yymmdd + suffix.") +args = parser.parse_args() + +cascade_time = Time(args.time, format='mjd') +year, month, day = cascade_time.iso.split('-') +day = day[:2] +casc_name = 'IceCube-Cascade_{}{}{}{}'.format(year[-2:], month, day, args.suffix) + +if 'https://roc.icecube.wisc.edu' in args.skymap: + print('Downloading skymap from https://roc.icecube.wisc.edu') + saved_skymaps_path = os.environ.get('FAST_RESPONSE_OUTPUT') + '/../cascade_skymaps/' + skymap_filename=args.skymap.split('/')[-1] + if not os.path.isdir(saved_skymaps_path): + subprocess.call(['mkdir', saved_skymaps_path]) + subprocess.call(['wget', args.skymap,'--no-check-certificate']) + subprocess.call(['mv',skymap_filename,saved_skymaps_path]) + print('Done.') + args.skymap=saved_skymaps_path+skymap_filename + +all_results = {} +for delta_t in [2.*86400.]: #1000 + start_time = cascade_time - (delta_t / 86400. / 2.) + stop_time = cascade_time + (delta_t / 86400. / 2.) + start = start_time.iso + stop = stop_time.iso + + name = casc_name + ' {:.1e}_s'.format(delta_t) + name = name.replace('_', ' ') + + run_id = args.alert_id[0][0] + ev_id = args.alert_id[0][1] + + f = CascadeFollowup(name, args.skymap, start, stop, skipped=args.alert_id) + t_mask=(f.llh.exp['time']<=f.stop)&(f.llh.exp['time']>=f.start) + print("__________________________________",f.start,f.stop) + print(len(f.llh.exp),f.llh.exp.dtype.names, + len(f.llh.exp[t_mask])) + print(f.llh.exp['time']) + print(max(f.llh.exp['time'])) + print(max(f.exp['time'])) + + f.unblind_TS() + f.plot_ontime() + f.calc_pvalue() + f.make_dNdE() + f.plot_tsd() + f.upper_limit() + f.find_coincident_events() + results = f.save_results() + f.make_prior_map(plotting_location=f.analysispath + '/',format=True) + + f.generate_report() + + nsToTest=[.25,.5,.75,1,1.25,1.5,1.75,2,2.25,2.50,2.75,3,3.25,3.5,3.75,4,4.25,4.50,4.75,5] + f.make_posterior_map(nsToTest=nsToTest, + prior_func=lambda ns,gamma,ra,dec: 1, + plotting_location=f.analysispath + '/', + ) + + + all_results[delta_t] = results diff --git a/fast_response/scripts/test_posterior_code.py b/fast_response/scripts/test_posterior_code.py new file mode 100644 index 0000000..8bbf67a --- /dev/null +++ b/fast_response/scripts/test_posterior_code.py @@ -0,0 +1,92 @@ +#!/usr/bin/env python + +r''' +Author: Sam Hori +July 2024''' + +import argparse, subprocess +from astropy.time import Time +import pyfiglet + +from fast_response.GWFollowup import GWFollowup +#import fast_response.web_utils as web_utils + +parser = argparse.ArgumentParser(description='GW Followup') +parser.add_argument('--skymap', type=str, default=None, + help='path to skymap (can be a web link for GraceDB (LVK) or a path') +parser.add_argument('--time', type=float, default=None, + help='Time of the GW (mjd)') +parser.add_argument('--name', type=str, + default="name of GW event being followed up") +parser.add_argument('--tw', default = 1000, type=int, #2 week: 1218240 + help = 'Time window for the analysis, in sec (default = 1000)') +parser.add_argument('--n_inj', default = 0, type=int, + help = 'Number to inject') +parser.add_argument('--allow_neg_ts', type=bool, default=False, + help='bool to allow negative TS values in gw analysis.') +args = parser.parse_args() + +#GW message header +message = '*'*80 +message += '\n' + str(pyfiglet.figlet_format("GW Followup")) + '\n' +message += '*'*80 +print(message) + +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.) +else: #2 week followup + 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 + +name = args.name +name = name.replace('_', ' ') + +f = GWFollowup(name, args.skymap, start, stop) +t_mask=(f.llh.exp['time']<=f.stop)&(f.llh.exp['time']>=f.start) + +print("__________________________________",f.start,f.stop) +print(len(f.llh.exp),f.llh.exp.dtype.names, + len(f.llh.exp[t_mask])) +print(f.llh.exp['time']) + +f._allow_neg = args.allow_neg_ts +f.save_items['merger_time'] = Time(gw_time, format='mjd').iso +f.save_items['skymap_link'] = args.skymap +f.initialize_injector() +injected_events=f.inject_events(args.n_inj) +f.unblind_TS(custom_events=injected_events) +f.make_prior_map(plotting_location=f.analysispath + '/',format=True) + +if args.tw>1000: + f.plot_ontime(label_events = False,custom_events=injected_events) +else: + f.plot_ontime(custom_events=injected_events) + +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(custom_events=injected_events) +f.per_event_pvalue() + +results = f.save_results() +nsToTest=[.25,.5,.75,1,1.25,1.5,1.75,2,2.25,2.50,2.75,3,3.25,3.5,3.75,4] +f.make_posterior_map(nsToTest=nsToTest, + prior_func=lambda ns,gamma,ra,dec: 1, + plotting_location=f.analysispath + '/', + custom_events=injected_events) +f.generate_report() + +f.write_circular() + From 6461812a088a84290e95448dbaa5eb41d06f817c Mon Sep 17 00:00:00 2001 From: Sam Hori Date: Wed, 8 Jan 2025 12:37:10 -0600 Subject: [PATCH 3/7] Added changes --- fast_response/AlertFollowup.py | 2 +- fast_response/FastResponseAnalysis.py | 263 ++++++++++++++++++++++---- fast_response/GWFollowup.py | 27 ++- fast_response/plotting_utils.py | 17 +- 4 files changed, 259 insertions(+), 50 deletions(-) diff --git a/fast_response/AlertFollowup.py b/fast_response/AlertFollowup.py index f574778..3b01359 100644 --- a/fast_response/AlertFollowup.py +++ b/fast_response/AlertFollowup.py @@ -370,4 +370,4 @@ class CascadeFollowup(AlertFollowup): _dataset = "GFUOnline_v001p02" _fix_index = True _float_index = not _fix_index - _index = 2.5 \ No newline at end of file + _index = 2.5 diff --git a/fast_response/FastResponseAnalysis.py b/fast_response/FastResponseAnalysis.py index fbb6451..4ca3428 100644 --- a/fast_response/FastResponseAnalysis.py +++ b/fast_response/FastResponseAnalysis.py @@ -7,6 +7,7 @@ from abc import abstractmethod import os, sys, time, subprocess import pickle, dateutil.parser, logging, warnings +from collections import Counter import h5py import healpy as hp @@ -14,8 +15,12 @@ import seaborn as sns import matplotlib as mpl import matplotlib.pyplot as plt +import copy from astropy.time import Time +import astropy.io.fits as fits + from scipy.special import erfinv, logsumexp +import scipy.stats as spstats from matplotlib.lines import Line2D from skylab.datasets import Datasets @@ -186,16 +191,18 @@ def get_data(self, livestream_start=None, livestream_stop=None): grls.append(grl) exp = np.concatenate(exps) grl = np.concatenate(grls) - else: + else: if self._verbose: print("Recent time: querying the i3live database") if livestream_start is None or livestream_stop is None: livestream_start = self.start - 6. livestream_stop = self.stop + print("livestream",livestream_start, livestream_stop) exp, mc, livetime, grl = dset.livestream( livestream_start, livestream_stop, append=self._season_names, floor=self._floor) + print(max(exp["time"])) exp.sort(order='time') grl.sort(order='run') livetime = grl['livetime'].sum() @@ -359,7 +366,7 @@ def calc_pvalue(self, ntrials=1000, run_anyway=False): return p @abstractmethod - def find_coincident_events(self): + def find_coincident_events(self,custom_events=None): pass @abstractmethod @@ -489,7 +496,7 @@ def save_results(self, alt_path=None): print("Results successfully saved") return self.save_items - def plot_ontime(self, with_contour=False, contour_files=None, label_events=False): + def plot_ontime(self, with_contour=False, contour_files=None, label_events=False,custom_events=None): r"""Plots ontime events on the full skymap and a zoomed in version near the scan best-fit @@ -505,16 +512,16 @@ def plot_ontime(self, with_contour=False, contour_files=None, label_events=False """ try: - self.plot_skymap_zoom(with_contour=with_contour, contour_files=contour_files) + self.plot_skymap_zoom(with_contour=with_contour, contour_files=contour_files,custom_events=custom_events) except Exception as e: print('Failed to make skymap zoom plot') try: - self.plot_skymap(with_contour=with_contour, contour_files=contour_files, label_events=label_events) + self.plot_skymap(with_contour=with_contour, contour_files=contour_files, label_events=label_events,custom_events=custom_events) except Exception as e: print('Failed to make FULL skymap plot') - def plot_skymap_zoom(self, with_contour=False, contour_files=None): + def plot_skymap_zoom(self, with_contour=False, contour_files=None,custom_events=None): r"""Make a zoomed in portion of a skymap with all ontime neutrino events within a certain range Outputs a plot (in png and pdf formats) to the analysis path @@ -529,7 +536,9 @@ def plot_skymap_zoom(self, with_contour=False, contour_files=None): """ events = self.llh.exp events = events[(events['time'] < self.stop) & (events['time'] > self.start)] - + print(custom_events) + if(custom_events is not None): + events=np.concatenate([events,custom_events]) col_num = 5000 seq_palette = sns.color_palette("icefire", col_num) lscmap = mpl.colors.ListedColormap(seq_palette) @@ -607,7 +616,7 @@ def plot_skymap_zoom(self, with_contour=False, contour_files=None): plt.savefig(self.analysispath + '/' + self.analysisid + 'unblinded_skymap_zoom.pdf',bbox_inches='tight', dpi=300) plt.close() - def plot_skymap(self, with_contour=False, contour_files=None, label_events=False): + def plot_skymap(self, with_contour=False, contour_files=None, label_events=False,custom_events=None): r""" Make skymap with event localization and all neutrino events on the sky within the given time window Outputs a plot in png format to the analysis path @@ -625,6 +634,8 @@ def plot_skymap(self, with_contour=False, contour_files=None, label_events=False events = self.llh.exp events = events[(events['time'] < self.stop) & (events['time'] > self.start)] + if(custom_events is not None): + events=np.concatenate([events,custom_events]) col_num = 5000 seq_palette = sns.color_palette("icefire", col_num) @@ -872,7 +883,7 @@ def run_background_trials(self, ntrials=1000): self.tsd = tsd self.save_items['tsd'] = tsd - def find_coincident_events(self): + def find_coincident_events(self,custom_events=None): r"""Find coincident events for a skymap based analysis. These are ontime events that are also in the 90% contour of the skymap @@ -897,6 +908,14 @@ def find_coincident_events(self): coincident_events = [] else: coincident_events = [] + if(custom_events is not None): + for ev in custom_events: + self.coincident_events.append({}) + self.coincident_events[-1]['delta_psi'] = -1 + self.coincident_events[-1]['spatial_w'] = -1 + self.coincident_events[-1]['energy_w'] = -1 + for key in ['run', 'event', 'ra', 'dec', 'sigma', 'logE', 'time']: + self.coincident_events[-1][key] = ev[key] for ev in events: coincident_events.append({}) for key in ['run', 'event', 'ra', 'dec', 'sigma', 'logE', 'time']: @@ -938,7 +957,8 @@ def unblind_TS(self, custom_events=None): val = self.llh.scan( 0.0,0.0, scramble = False, spatial_prior=spatial_prior, time_mask = [self.duration/2., self.centertime], - pixel_scan=[self.nside, self._pixel_scan_nsigma] + pixel_scan=[self.nside, self._pixel_scan_nsigma], + inject=custom_events ) self.ts_scan = val t2 = time.time() @@ -992,10 +1012,70 @@ def unblind_TS(self, custom_events=None): return ts, ns, gamma_fit else: return ts, ns - + def inject_events(self,ninj=0): + ni, sample,src_ra, src_dec = self.inj.sample(ninj,poisson=False,return_position =True) + return sample - - def make_posterior_map(self, nsToTest=[1,2,3,4,5], prior=lambda ns,gamma,ra,dec: 1, testingPlotsLocation=None): + + def make_prior_map(self, + plotting_location=None, + format=True + ): + #TODO: Document + """ + Generates a zoomed in skymap around the best fit location from the maximum likelihood analysis + """ + + + #Format the skymap if told to format (usually formatted at construction so shouldn't change anything) + + if(format): + skymap = self.format_skymap(self.skymap) + else: + skymap=self.skymap + + #Plot the prior skymap + fig, ax=plt.subplots(1,1,figsize=(10,10)) + plt.sca(ax) + reso=.1 + sizeZoomInDegs=10 + nPixZoom=sizeZoomInDegs/reso*60 + hp.visufunc.gnomview(map=np.log10(skymap),hold=True,title="Prior Skymap",rot=(np.degrees(self.skymap_fit_ra),np.degrees(self.skymap_fit_dec),0), + xsize=nPixZoom,ysize=nPixZoom,reso=reso,bgcolor="white",badcolor="white",margins=(.01,.01,0,0),notext=True,cbar=False) + #hp.visufunc.cartview(map=np.log10(skymap),hold=True,title="Prior Skymap",rot=(np.degrees(self.skymap_fit_ra),np.degrees(self.skymap_fit_dec),0), + # lonra=[np.degrees(self.skymap_fit_ra)-5,np.degrees(self.skymap_fit_ra)+5],latra=[np.degrees(self.skymap_fit_dec)-5,np.degrees(self.skymap_fit_dec)+5],bgcolor="white",badcolor="white",margins=(.01,.01,0,0),notext=True,) + + hp.graticule(verbose=True) + plotting_utils.plot_labels(self.skymap_fit_dec,self.skymap_fit_ra, reso,nPix=nPixZoom,vertical_height_xlabel=.95) + plotting_utils.plot_color_bar(range=[0,np.nanmax(self.skymap)], cmap="viridis", col_label=r"Skymap Prior PDF", + offset=-50,labels=[0, "{0:.1e}".format(np.nanmax(self.skymap)/2),"{0:.2e}".format(np.nanmax(self.skymap))],loc=[0.90, 0.2, 0.03, 0.6]) + name="PriorSkymap" + if(not format): + name+="Unformatted" + plt.savefig(plotting_location+name+".png") + + def _make_posterior_fits(self,pos_skymap): + #internal function to be called from make_posterior_skymap + extra_header = [('index', self.index), + #('energy_range', ), + #('ns_range', ), + ("Sender","IceCube Collaboration"), + ("TRIGGER_EVENT",self.name), + ("startmjd",self.start), + ("stopmjd",self.stop) + + ] + #startmjd, stopmjd, + hp.write_map(self.analysispath + '/' +'posterior_info.fits',m=pos_skymap, + extra_header=extra_header,partial=True + + ) + + def make_posterior_map(self, nsToTest=[1,2,3,4,5], + prior_func=lambda ns,gamma,ra,dec: 1, + plotting_location=None, + custom_events=None + ): # Posterior Approach # We want a function P(RA, DEC | X) where X is the neutrino data # Compute the entire likelihood space P(X | RA, DEC, ns, gamma) @@ -1003,32 +1083,66 @@ def make_posterior_map(self, nsToTest=[1,2,3,4,5], prior=lambda ns,gamma,ra,dec: #ns to test give as list #Give prior as a function f(ns, gamma, ra, dec) + #As currently implemented,the llh will have gamma must be constant set by the self.index + #Prior flat in ns and spatially by default. May make sense to implement #TODO: flux/ns priors + t1 = time.time() + val=None + skymap = self.skymap + + spatial_prior = SpatialPrior( - self.skymap, + skymap, containment = self._containment,allow_neg=self._allow_neg) + + nside=self.nside + + #The full sky scan with no fixing (maybe better to save this during the unblind TS stage) + full_scan_val = self.llh.scan(0.0,0.0, scramble = False,spatial_prior=spatial_prior, + time_mask = [self.duration/2., self.centertime], + pixel_scan=[nside, self._pixel_scan_nsigma], + fixed=['nsignal'], + inject=custom_events + ) + + + for nsigT in nsToTest: + print('testing nsig=',nsigT) temp_val = self.llh.scan(0.0,0.0, scramble = False,spatial_prior=spatial_prior, time_mask = [self.duration/2., self.centertime], - pixel_scan=[self.nside, self._pixel_scan_nsigma], + pixel_scan=[nside, self._pixel_scan_nsigma], fixed=['nsignal'], - nsignal=np.array([nsigT]*12*pow(self.nside,2),dtype=float) #TODO: Check if we can do more than 1 at a time + nsignal=np.array([nsigT]*12*pow(nside,2),dtype=float), + inject=custom_events ) if val is None: val=temp_val else: val=np.hstack((val,temp_val)) + print(val.dtype.names) ras=val["ra"] decs=val["dec"] nss=val["nsignal"] - gammas=val["index"] + + if "gamma" in val.dtype.names: + #TODO: This is available in some, but not all, analyses -- GW vs Internal alerts + # should always be self.index (unless things get) + # maybe "spectrum" instead + + gammas=val["gamma"] + + #gammas=val["gamma"] TS_spatial_prior_0=val["TS_spatial_prior_0"] + idx_max=np.argmax(TS_spatial_prior_0) + ramax,decmax=ras[idx_max],decs[idx_max] logProbs=TS_spatial_prior_0/2 #not normalized so technically C*ln(P) finalProbs={} for i in range(len(ras)): - prior=1 #TODO + dV=1 #integration factor #TODO: implement if uneven intervals are given in ns + prior=prior_func(nss[i],self.index,ras[i],decs[i])*dV if((ras[i],decs[i]) in finalProbs): finalProbs[(ras[i],decs[i])]=logsumexp([finalProbs[(ras[i],decs[i])],logProbs[i]+np.log(prior)]) else: @@ -1036,45 +1150,108 @@ def make_posterior_map(self, nsToTest=[1,2,3,4,5], prior=lambda ns,gamma,ra,dec: #Normalize totalLog=logsumexp(list(finalProbs.values())) finalProbs={k:v-totalLog for k,v in finalProbs.items()} - outputProbHPMap=np.zeros(12*self.nside**2) - outputProbHPMap[hp.ang2pix(nside=self.nside,theta=np.array(ras_graph)*180/np.pi, phi=np.array(decs_graph)*180/np.pi,lonlat=True)]=np.exp(np.array(list(finalProbs.values()))) + + #Write to healpy + outputProbHPMap=np.zeros(12*nside**2) + ras_graph,decs_graph=zip(*list(finalProbs.keys())) + outputProbHPMap[hp.ang2pix(nside=nside,theta=np.array(ras_graph)*180/np.pi, phi=np.array(decs_graph)*180/np.pi,lonlat=True)]=np.exp(np.array(list(finalProbs.values()))) + t2 = time.time() + print("finished posterior calc, took {} s".format(t2-t1)) + #set 0s to nan for plotting - outputProbHPMap[outputProbHPMap == 0] = np.nan + outputProbHPMapPlotting=copy.deepcopy(outputProbHPMap) + outputProbHPMapPlotting[outputProbHPMapPlotting == 0] = np.nan - if (testingPlotsLocation is not None): + if (plotting_location is not None): fig, ax=plt.subplots(1,1,figsize=(10,10)) ras_graph,decs_graph=zip(*list(finalProbs.keys())) fig, ax=plt.subplots(1,1,figsize=(10,10)) plt.sca(ax) - reso=3 - nPixZoom=200 - sizeZoomInDegs=reso/60*nPixZoom + reso=.1 + sizeZoomInDegs=10 + nPixZoom=sizeZoomInDegs/reso*60 + + + #hp.visufunc.gnomview(map=outputProbHPMap,hold=True,title="Zoomed Posterior Skymap",rot=(self.src_ra*180/np.pi,self.src_dec*180/np.pi,0), + # xsize=nPixZoom,ysize=nPixZoom,reso=reso,bgcolor="white",badcolor="white",margins=(.01,.01,0,0),notext=True,) + + hp.visufunc.mollview(map=outputProbHPMapPlotting,hold=True,title="Posterior Skymap",bgcolor="white",badcolor="white",rot=(180,0,0)) + hp.graticule(verbose=True) + + plt.savefig(plotting_location+"PosteriorMap.png") + plt.close() + + #plot zoomed + fig, ax=plt.subplots(1,1,figsize=(10,10)) + plt.sca(ax) + hp.visufunc.gnomview(map=outputProbHPMapPlotting,hold=True,title="Posterior Skymap",rot=(ramax*180/np.pi,decmax*180/np.pi,0), + xsize=nPixZoom,ysize=nPixZoom,reso=reso,bgcolor="white",badcolor="grey",margins=(.01,.01,0,0),notext=True,cbar=False) + hp.graticule(verbose=True) + plotting_utils.plot_labels(decmax,ramax, reso,nPix=nPixZoom,vertical_height_xlabel=.95) + plotting_utils.plot_color_bar(range=[0,np.nanmax(outputProbHPMapPlotting)], cmap="viridis", col_label=r"Posterior PDF", + offset=-35,labels=[0,"{0:.1e}".format(np.nanmax(outputProbHPMapPlotting)/2),"{0:.2e}".format(np.nanmax(outputProbHPMapPlotting))],loc=[0.90, 0.2, 0.03, 0.6]) + plt.savefig(plotting_location+"PosteriorMapZoomed.png") + plt.close() + + fig, ax=plt.subplots(1,1,figsize=(10,10)) + plt.sca(ax) + hp.visufunc.gnomview(map=np.log10(outputProbHPMapPlotting),hold=True,title="Posterior Skymap (Log p)",rot=(ramax*180/np.pi,decmax*180/np.pi,0), + xsize=nPixZoom,ysize=nPixZoom,reso=reso,bgcolor="white",badcolor="grey",margins=(.01,.01,0,0),notext=True,cbar=False) - """ - hp.visufunc.gnomview(map=outputProbHPMap,hold=True,title="Zoomed Posterior Skymap",rot=(src_ra*180/np.pi,src_dec*180/np.pi,0), - xsize=nPixZoom,ysize=nPixZoom,reso=reso,bgcolor="white",badcolor="white",margins=(.01,.01,0,0),notext=True,) - """ - hp.visufunc.mollview(map=outputProbHPMap,hold=True,title="Posterior Skymap",bgcolor="white",badcolor="white",rot=(180,0,0)) - plt.savefig(testingPlotsLocation+"PosteriorMap.png") + hp.graticule(verbose=True) + plotting_utils.plot_labels(decmax,ramax, reso,nPix=nPixZoom,vertical_height_xlabel=.95) + plotting_utils.plot_color_bar(range=[0,np.nanmax(np.log10(outputProbHPMapPlotting))], cmap="viridis", col_label=r"Posterior PDF", + offset=-35,labels=[0,"{0:.1e}".format(np.nanmax(np.log10(outputProbHPMapPlotting))/2),"{0:.2e}".format(np.nanmax(np.log10(outputProbHPMapPlotting)))],loc=[0.90, 0.2, 0.03, 0.6]) + plt.savefig(plotting_location+"PosteriorMapZoomedLog.png") plt.close() - for (ra,dec) in finalProbs: + #plot ns vs probability for some specific points, could use a cleaner method of choosing points + threshold_prob_for_plots=list(sorted([k for k, v in Counter(list(finalProbs.values())).items() if v == 1],reverse=True)) + #threshold_prob_for_plots=[k for k in threshold_prob_for_plots if k>0] + counterJump=int(len(threshold_prob_for_plots)/30) + print(len(threshold_prob_for_plots),counterJump) + + psToPlot=[] + for i in range(0,len(threshold_prob_for_plots),counterJump): + p=threshold_prob_for_plots[i] + psToPlot.append(p) + if (i>len(threshold_prob_for_plots)): + break + filtered_probs = {k:v for (k,v) in finalProbs.items() + if v in psToPlot} + + cmap= mpl.cm.ScalarMappable( + norm = mpl.colors.Normalize(min(threshold_prob_for_plots),max(threshold_prob_for_plots)), + cmap = plt.get_cmap('viridis')) + + for (ra,dec) in filtered_probs: + #print(ra,dec) total_val=finalProbs[(ra,dec)] temp_val=val[(val["ra"]==ra)&(val["dec"]==dec)] nss_temp=temp_val["nsignal"] - gammas_temp=temp_val["index"] + #gammas_temp=temp_val["gamma"] TS_temp=temp_val["TS_spatial_prior_0"]/2 totalLog=logsumexp(TS_temp) p_temp=np.exp(TS_temp-totalLog) - plt.plot(nss_temp,p_temp,alpha=total_val) - plt.colorbar() - plt.xlabel("ns") - plt.ylabel("TS") - plt.savefig(testingPlotsLocation+"nsProfilePlot.png") + plt.plot(nss_temp,p_temp,color=cmap.to_rgba(total_val)) + + + cmap.set_array([]) + cbar=plt.colorbar(cmap) + cbar.set_label('Pixel Probability') + plt.xlabel("ns",fontsize=15) + plt.ylabel("probability",fontsize=15) + plt.savefig(plotting_location+"nsProfilePlot.png") plt.close() + + print("saving fits") + self._make_posterior_fits(outputProbHPMap) + #return the version with zeros rather than nan + t3 = time.time() + print("finished posterior saving, took {} s".format(t3-t1)) return outputProbHPMap @@ -1306,7 +1483,7 @@ def unblind_TS(self): self.save_items['ns'] = ns return ts, ns - def find_coincident_events(self): + def find_coincident_events(self,custom_events=None): r"""Find "coincident events" for the analysis. These are ontime events that satisfy: @@ -1323,6 +1500,14 @@ def find_coincident_events(self): temporal_weights = self.llh.temporal_model.signal(self.llh._events) msk = spatial_weights * energy_ratio * temporal_weights > 10 self.coincident_events = [] + if(custom_events!=None): + for ev in custom_events: + self.coincident_events.append({}) + self.coincident_events[-1]['delta_psi'] = -1 + self.coincident_events[-1]['spatial_w'] = -1 + self.coincident_events[-1]['energy_w'] = -1 + for key in ['run', 'event', 'ra', 'dec', 'sigma', 'logE', 'time']: + self.coincident_events[-1][key] = ev[key] if len(spatial_weights[msk]) > 0: self.coincident_events = [] for ev, s_w, en_w in zip(self.llh._events[msk], diff --git a/fast_response/GWFollowup.py b/fast_response/GWFollowup.py index 3dbd6f0..1b1620d 100644 --- a/fast_response/GWFollowup.py +++ b/fast_response/GWFollowup.py @@ -351,8 +351,8 @@ def get_best_fit_contour(self, proportions=[0.5,0.9]): 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, label_events=True): - return super().plot_ontime(with_contour=True, contour_files=contour_files, label_events=label_events) + def plot_ontime(self, with_contour=True, contour_files=None, label_events=True,custom_events=None): + return super().plot_ontime(with_contour=True, contour_files=contour_files, label_events=label_events,custom_events=custom_events) def write_circular(self): """ @@ -541,7 +541,7 @@ def inject_scan(self, ra, dec, ns, poisson=True): results = dict([('ts',ts),('ns',ns),('gamma',gamma),('ra',ra),('dec',dec)]) return (results, events) - def find_coincident_events(self): + def find_coincident_events(self,custom_events=None): r""" Find "coincident events" for a skymap based analysis. These are ALL ontime events, @@ -573,6 +573,27 @@ def find_coincident_events(self): events['ts'][i] = self.ts_scan['TS_spatial_prior_0'][idx[0]] events['ns'][i] = self.ts_scan['nsignal'][idx[0]] events['gamma'][i] = self.ts_scan['gamma'][idx[0]] + + if(custom_events is not None): + ''' + custom_events.dtype.names = 'time', 'run', 'event', 'ra', 'dec', 'azimuth', 'zenith', 'sigma', 'logE', 'subevent' + dtype = [('time',np.float32),('run',np.uint32),('event',np.uint32),('ra',np.float32), + ('dec',np.float32),('azimuth',np.float32),('zenith',np.float32),('sigma',np.float32), + ('logE',np.float32),('subevent',np.float64)] + custom_events=custom_events[['time', 'run', 'event', 'ra', 'dec', 'azimuth', 'zenith', 'sigma', 'logE', 'subevent']] + ''' + custom_events = append_fields( + custom_events, names=['in_contour', 'ts', 'ns', 'gamma', 'B'], + data=np.empty((5, custom_events['ra'].size)), + usemask=False) + for i in range(custom_events['ra'].size): + custom_events['in_contour'][i]=True + custom_events['B'][i] = -1 + custom_events['ts'][i] = -1 + custom_events['ns'][i] = -1 + custom_events['gamma'][i] = -1 + events=np.concatenate([events,custom_events]) + self.events_rec_array = events self.coincident_events = [dict(zip(events.dtype.names, x)) for x in events] diff --git a/fast_response/plotting_utils.py b/fast_response/plotting_utils.py index b024025..f383d12 100644 --- a/fast_response/plotting_utils.py +++ b/fast_response/plotting_utils.py @@ -29,11 +29,11 @@ def plot_zoom(scan, ra, dec, title, reso=3, var="pVal", range=[0, 6],cmap=None): cmap = mpl.colors.ListedColormap(pdf_palette) hp.gnomview(scan, rot=(np.degrees(ra), np.degrees(dec), 0), cmap=cmap, + cbar=False, max=max(scan), reso=reso, title=title, notext=True, - cbar=False #unit=r"" ) @@ -41,7 +41,7 @@ def plot_zoom(scan, ra, dec, title, reso=3, var="pVal", range=[0, 6],cmap=None): hp.graticule(verbose=False) plot_labels(dec, ra, reso) -def plot_color_bar(labels=[0.,2.,4.,6.], col_label=r"IceCube Event Time", range=[0,6], cmap=None, offset=-35): +def plot_color_bar(labels=[0.,2.,4.,6.], col_label=r"IceCube Event Time", range=[0,6], cmap=None, offset=-35,loc=[0.95, 0.2, 0.03, 0.6]): """ Adds a color bar to an existing healpy map @@ -55,9 +55,11 @@ def plot_color_bar(labels=[0.,2.,4.,6.], col_label=r"IceCube Event Time", range= colormap to use. if not set default is "Blues" offset: int offset value for colorbar's label. default is -35 + loc: float list of length 4 + location for the colorbar axis in form [x,y,width,height] """ fig = plt.gcf() - ax = fig.add_axes([0.95, 0.2, 0.03, 0.6]) + ax = fig.add_axes(loc) labels = labels cb = mpl.colorbar.ColorbarBase(ax, cmap="Blues" if cmap is None else cmap, #norm=mpl.colors.Normalize(vmin=range[0], vmax=range[1]), @@ -68,7 +70,7 @@ def plot_color_bar(labels=[0.,2.,4.,6.], col_label=r"IceCube Event Time", range= cb.set_ticklabels(labels) cb.update_ticks() -def plot_labels(src_dec, src_ra, reso): +def plot_labels(src_dec, src_ra, reso,nPix=200,vertical_height_xlabel=1): """ Add labels to healpy zoom @@ -82,6 +84,7 @@ def plot_labels(src_dec, src_ra, reso): Resolution (arcmins) """ fontsize = 20 + reso=reso*nPix/200 plt.text(-1*np.radians(1.75*reso),np.radians(0), r"%.1f$^{\circ}$"%(np.degrees(src_dec)), horizontalalignment='right', verticalalignment='center', fontsize=fontsize) @@ -102,7 +105,7 @@ def plot_labels(src_dec, src_ra, reso): verticalalignment='top', fontsize=fontsize) plt.text(-1*np.radians(2.35*reso), np.radians(0), r"declination", ha='center', va='center', rotation=90, fontsize=fontsize) - plt.text(np.radians(0), np.radians(-2.05*reso), r"right ascension", + plt.text(np.radians(0), np.radians(-2.05*reso*vertical_height_xlabel), r"right ascension", ha='center', va='center', fontsize=fontsize) def plot_events(dec, ra, sigmas, src_ra, src_dec, reso, sigma_scale=5., col = 'k', constant_sigma=False, @@ -308,8 +311,8 @@ def make_public_zoom_skymap(skymap, events, ra, dec, with_contour=True, name='te #plt.scatter(0,0, marker='*', c = 'k', s = 130, label = label_str) #plt.legend(loc = 2, ncol=1, mode = 'expand', fontsize = 18.5, framealpha = 0.95) + plot_color_bar(range=[0,6], cmap=cmap, col_label='GW Map Probability', offset = -10, - labels=[f'{min(skymap):.1e}',f'{max(skymap):.1e}']) - + labels=[f'{min(skymap):.1e}',f'{max(skymap):.1e}']) plt.savefig(f'./{name}_skymap_zoom_public.png', bbox_inches='tight', dpi=300) plt.close() \ No newline at end of file From 976a29aa9aca733e7aaaee4605a00fd382199ef0 Mon Sep 17 00:00:00 2001 From: Sam Hori Date: Tue, 18 Mar 2025 00:34:49 -0500 Subject: [PATCH 4/7] Added document+ability to use flux rather than ns --- fast_response/scripts/test_posterior_bash.sh | 3 --- fast_response/scripts/test_posterior_cascade.bash | 2 ++ fast_response/scripts/test_posterior_cascade.py | 6 +++++- fast_response/scripts/test_posterior_code.py | 3 ++- fast_response/scripts/test_posterior_gw.sh | 2 ++ 5 files changed, 11 insertions(+), 5 deletions(-) create mode 100644 fast_response/scripts/test_posterior_cascade.bash create mode 100644 fast_response/scripts/test_posterior_gw.sh diff --git a/fast_response/scripts/test_posterior_bash.sh b/fast_response/scripts/test_posterior_bash.sh index 92c5743..688f8a9 100644 --- a/fast_response/scripts/test_posterior_bash.sh +++ b/fast_response/scripts/test_posterior_bash.sh @@ -27,6 +27,3 @@ python test_posterior_code.py --skymap=https://gracedb.ligo.org/api/superevents/ #S230904n-4-Update python test_posterior_code.py --skymap=https://gracedb.ligo.org/api/superevents/S230904n/files/bayestar.fits.gz \ --time=60191.21542972222 --name="S230904n_Test_Posterior" --tw=1000 -#MS240802n-3-Initial_test -python test_posterior_code.py --skymap=https://gracedb.ligo.org/api/superevents/MS240802n/files/bayestar.fits.gz \ - --time=60524.544673460645 --name="MS240802n_Test_Posterior" --tw=1000 diff --git a/fast_response/scripts/test_posterior_cascade.bash b/fast_response/scripts/test_posterior_cascade.bash new file mode 100644 index 0000000..6120b19 --- /dev/null +++ b/fast_response/scripts/test_posterior_cascade.bash @@ -0,0 +1,2 @@ +python test_posterior_cascade.py --skymap=https://roc.icecube.wisc.edu/public/hese_cascades/hese_60505_run00139644.evt000063900267.fits --time=58000 --alert_id=139644:63900267 +python test_posterior_cascade.py --skymap=https://roc.icecube.wisc.edu/public/hese_cascades/hese_60505_run00139644.evt000063900267.fits --time=60505.623363 --alert_id=139644:63900267 \ No newline at end of file diff --git a/fast_response/scripts/test_posterior_cascade.py b/fast_response/scripts/test_posterior_cascade.py index 7bdd845..74a8358 100755 --- a/fast_response/scripts/test_posterior_cascade.py +++ b/fast_response/scripts/test_posterior_cascade.py @@ -76,15 +76,19 @@ f.upper_limit() f.find_coincident_events() results = f.save_results() + f.make_prior_map(plotting_location=f.analysispath + '/',format=True) f.generate_report() nsToTest=[.25,.5,.75,1,1.25,1.5,1.75,2,2.25,2.50,2.75,3,3.25,3.5,3.75,4,4.25,4.50,4.75,5] - f.make_posterior_map(nsToTest=nsToTest, + #nsToTest=[1e-7,2e-7,3e-7,4e-7,5e-7, 6e-7,7e-7,8e-7] + f.make_posterior_map(fluxToTest=nsToTest, prior_func=lambda ns,gamma,ra,dec: 1, plotting_location=f.analysispath + '/', + useFlux=False ) + f.generate_posterior_report() all_results[delta_t] = results diff --git a/fast_response/scripts/test_posterior_code.py b/fast_response/scripts/test_posterior_code.py index 8bbf67a..f0c76e0 100644 --- a/fast_response/scripts/test_posterior_code.py +++ b/fast_response/scripts/test_posterior_code.py @@ -82,11 +82,12 @@ results = f.save_results() nsToTest=[.25,.5,.75,1,1.25,1.5,1.75,2,2.25,2.50,2.75,3,3.25,3.5,3.75,4] -f.make_posterior_map(nsToTest=nsToTest, +f.make_posterior_map(fluxToTest=nsToTest, prior_func=lambda ns,gamma,ra,dec: 1, plotting_location=f.analysispath + '/', custom_events=injected_events) f.generate_report() +f.generate_posterior_report() f.write_circular() diff --git a/fast_response/scripts/test_posterior_gw.sh b/fast_response/scripts/test_posterior_gw.sh new file mode 100644 index 0000000..427f97a --- /dev/null +++ b/fast_response/scripts/test_posterior_gw.sh @@ -0,0 +1,2 @@ +python test_posterior_code.py --skymap=https://gracedb.ligo.org/apiweb/superevents/S250206dm/files/Bilby.fits.gz,0 \ + --time=60712.892713 --name="S250206dm_Posterior" --tw=1000 From 62b9dd5c283349a419370194ca9a53ec02580489 Mon Sep 17 00:00:00 2001 From: Sam Hori Date: Thu, 29 May 2025 13:15:45 -0500 Subject: [PATCH 5/7] Plotting and flux changes for posterior code --- fast_response/AlertFollowup.py | 2 + fast_response/FastResponseAnalysis.py | 270 ++++++++-- .../latex/posterior_report_general.tex | 70 +++ fast_response/plotting_utils.py | 41 +- .../reports/PosteriorReportGenerator.py | 463 ++++++++++++++++++ fast_response/reports/__init__.py | 3 + fast_response/scripts/test_1.sh | 4 + requirements.txt | 2 +- 8 files changed, 814 insertions(+), 41 deletions(-) create mode 100644 fast_response/latex/posterior_report_general.tex create mode 100644 fast_response/reports/PosteriorReportGenerator.py create mode 100644 fast_response/scripts/test_1.sh diff --git a/fast_response/AlertFollowup.py b/fast_response/AlertFollowup.py index 3b01359..a14d618 100644 --- a/fast_response/AlertFollowup.py +++ b/fast_response/AlertFollowup.py @@ -371,3 +371,5 @@ class CascadeFollowup(AlertFollowup): _fix_index = True _float_index = not _fix_index _index = 2.5 + _nside = 512 + diff --git a/fast_response/FastResponseAnalysis.py b/fast_response/FastResponseAnalysis.py index 4ca3428..c9f6e79 100644 --- a/fast_response/FastResponseAnalysis.py +++ b/fast_response/FastResponseAnalysis.py @@ -22,7 +22,7 @@ from scipy.special import erfinv, logsumexp import scipy.stats as spstats from matplotlib.lines import Line2D - +from matplotlib.text import Text from skylab.datasets import Datasets from skylab.llh_models import EnergyLLH from skylab.priors import SpatialPrior @@ -36,7 +36,7 @@ from . import web_utils from . import sensitivity_utils from . import plotting_utils -from .reports import FastResponseReport +from .reports import FastResponseReport,FastResponsePosteriorReport mpl.use('agg') current_palette = sns.color_palette('colorblind', 10) @@ -753,6 +753,14 @@ def generate_report(self): report.generate_report() report.make_pdf() + def generate_posterior_report(self): + r"""Generates report using class attributes + and the ReportGenerator Class + Makes full report and outputs as pdf + """ + report = FastResponsePosteriorReport(self) + report.generate_posterior_report() + report.make_pdf() class PriorFollowup(FastResponseAnalysis): """ @@ -1028,7 +1036,6 @@ def make_prior_map(self, #Format the skymap if told to format (usually formatted at construction so shouldn't change anything) - if(format): skymap = self.format_skymap(self.skymap) else: @@ -1041,14 +1048,25 @@ def make_prior_map(self, sizeZoomInDegs=10 nPixZoom=sizeZoomInDegs/reso*60 hp.visufunc.gnomview(map=np.log10(skymap),hold=True,title="Prior Skymap",rot=(np.degrees(self.skymap_fit_ra),np.degrees(self.skymap_fit_dec),0), - xsize=nPixZoom,ysize=nPixZoom,reso=reso,bgcolor="white",badcolor="white",margins=(.01,.01,0,0),notext=True,cbar=False) + xsize=nPixZoom,ysize=nPixZoom,reso=reso,bgcolor="white",badcolor="white",margins=(.02,.02,.05,0),notext=True,cbar=False) #hp.visufunc.cartview(map=np.log10(skymap),hold=True,title="Prior Skymap",rot=(np.degrees(self.skymap_fit_ra),np.degrees(self.skymap_fit_dec),0), # lonra=[np.degrees(self.skymap_fit_ra)-5,np.degrees(self.skymap_fit_ra)+5],latra=[np.degrees(self.skymap_fit_dec)-5,np.degrees(self.skymap_fit_dec)+5],bgcolor="white",badcolor="white",margins=(.01,.01,0,0),notext=True,) hp.graticule(verbose=True) - plotting_utils.plot_labels(self.skymap_fit_dec,self.skymap_fit_ra, reso,nPix=nPixZoom,vertical_height_xlabel=.95) + plotting_utils.plot_labels(self.skymap_fit_dec,self.skymap_fit_ra, reso,nPix=nPixZoom,label_scale=.97) + + ax=plt.gca() + text_objs = ax.findobj( + lambda artist: isinstance(artist, Text) + ) + for obj in text_objs: + obj.set_fontsize(30) plotting_utils.plot_color_bar(range=[0,np.nanmax(self.skymap)], cmap="viridis", col_label=r"Skymap Prior PDF", - offset=-50,labels=[0, "{0:.1e}".format(np.nanmax(self.skymap)/2),"{0:.2e}".format(np.nanmax(self.skymap))],loc=[0.90, 0.2, 0.03, 0.6]) + offset=-50,labels=[0, "{0:.1e}".format(np.nanmax(self.skymap)/2),"{0:.2e}".format(np.nanmax(self.skymap))],loc=[0.86, 0.2, 0.03, 0.6]) + theta, phi =plotting_utils.plot_contours([.9],self.skymap) + hp.projplot(theta[0], phi[0], linewidth=2., c='k') + for i in range(1, len(theta)): + hp.projplot(theta[i], phi[i], linewidth=2., c='k', label=None) name="PriorSkymap" if(not format): name+="Unformatted" @@ -1070,24 +1088,43 @@ def _make_posterior_fits(self,pos_skymap): extra_header=extra_header,partial=True ) - - def make_posterior_map(self, nsToTest=[1,2,3,4,5], + + def fluxToAverageNs(self,flux): + if(self.inj==None): + print("Initialize the injector. Can't convert flux to ns without an injector") + return self.inj.flux2mu(flux) + def fluxToNs_decDependent(self,flux,dec): + mean_signal=self.fluxToAverageNs(flux) + acceptance_weights = self.inj.relative_signal_acceptance(np.sin(dec)) + num=[np.array(np.around(w * mean_signal / float(self.prior.nprior)), + dtype=np.int).tolist() + for w in acceptance_weights] + num_sum = np.sum(num, dtype=int) + return num_sum + + def make_posterior_map(self, fluxToTest=[1,2,3,4,5], prior_func=lambda ns,gamma,ra,dec: 1, plotting_location=None, - custom_events=None + custom_events=None, + useFlux=False ): # Posterior Approach # We want a function P(RA, DEC | X) where X is the neutrino data # Compute the entire likelihood space P(X | RA, DEC, ns, gamma) # Apply Bayes rule, marginalize over nuisance parameters to get P(RA, DEC | X) - #ns to test give as list + #flux/ns to test: given as list #Give prior as a function f(ns, gamma, ra, dec) #As currently implemented,the llh will have gamma must be constant set by the self.index #Prior flat in ns and spatially by default. May make sense to implement + + #plotting_location is the directory where plots should be output + #custom events will be injected + #useFlux determines whether fluxToTest is treated as ns or flux + + #TODO: flux/ns priors t1 = time.time() - val=None skymap = self.skymap @@ -1099,6 +1136,7 @@ def make_posterior_map(self, nsToTest=[1,2,3,4,5], nside=self.nside #The full sky scan with no fixing (maybe better to save this during the unblind TS stage) + #Can be removed for speed full_scan_val = self.llh.scan(0.0,0.0, scramble = False,spatial_prior=spatial_prior, time_mask = [self.duration/2., self.centertime], pixel_scan=[nside, self._pixel_scan_nsigma], @@ -1106,22 +1144,35 @@ def make_posterior_map(self, nsToTest=[1,2,3,4,5], inject=custom_events ) + decs_list=full_scan_val["dec"] - for nsigT in nsToTest: - print('testing nsig=',nsigT) + test_flux=1e-17 + if(useFlux): + mean_flux_arr=fluxToTest + nsPerFluxArr=[self.fluxToNs_decDependent(test_flux,d)/test_flux for d in decs_list] + else: + #work in uniform ns + mean_flux_arr=fluxToTest + nsPerFluxArr=[1]*12*pow(nside,2) + + + for fluxT in mean_flux_arr: + #print('testing flux=',fluxT) temp_val = self.llh.scan(0.0,0.0, scramble = False,spatial_prior=spatial_prior, time_mask = [self.duration/2., self.centertime], pixel_scan=[nside, self._pixel_scan_nsigma], fixed=['nsignal'], - nsignal=np.array([nsigT]*12*pow(nside,2),dtype=float), + nsignal=np.array([fluxT*conversion for conversion in nsPerFluxArr],dtype=float), + #nsignal=np.array([fluxToTest*nsPerFlux for nsPerFlux in nsPerFluxArr],dtype=float), + inject=custom_events ) if val is None: val=temp_val else: val=np.hstack((val,temp_val)) - print(val.dtype.names) + #print(val.dtype.names) ras=val["ra"] decs=val["dec"] nss=val["nsignal"] @@ -1140,17 +1191,37 @@ def make_posterior_map(self, nsToTest=[1,2,3,4,5], logProbs=TS_spatial_prior_0/2 #not normalized so technically C*ln(P) finalProbs={} - for i in range(len(ras)): - dV=1 #integration factor #TODO: implement if uneven intervals are given in ns - prior=prior_func(nss[i],self.index,ras[i],decs[i])*dV - if((ras[i],decs[i]) in finalProbs): - finalProbs[(ras[i],decs[i])]=logsumexp([finalProbs[(ras[i],decs[i])],logProbs[i]+np.log(prior)]) - else: - finalProbs[(ras[i],decs[i])]=logProbs[i]+np.log(prior) + + + #for each pixel we have lists [ns_1,ns_2,...ns_i], [logP_1,logP_2,...,logP_3] + #we want to integrate over ns + + pixel_list_full=list(zip(ras,decs)) + for pix in set(pixel_list_full): + idxes=np.where((ras==pix[0])&(decs==pix[1]))[0] + sort_idx=np.argsort(nss[idxes]) + idxes=idxes[sort_idx] + nss_temp=nss[idxes] + probs_temp=logProbs[idxes] + ns_diff=np.diff(nss_temp) + + finalProbs[pix]=-np.inf + if(len(idxes)!=len(fluxToTest)): + print("Not all fluxes have probabilities") + print("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!") + exit() + for i in range(1,len(nss_temp)): + #trapizoid wall + #+(P[i]+P[i-1])/2*prior(center) + ns_val=(nss_temp[i]+nss_temp[i-1])/2 + deltaNs=nss_temp[i]-nss_temp[i-1] + prior=prior_func(ns_val,self.index,pix[0],pix[1]) + finalProbs[pix]=logsumexp([finalProbs[pix],logsumexp([probs_temp[i],probs_temp[i-1]])+np.log(prior)+np.log(deltaNs/2)]) + #Normalize totalLog=logsumexp(list(finalProbs.values())) finalProbs={k:v-totalLog for k,v in finalProbs.items()} - + print("N_pixels",len(finalProbs)) #Write to healpy outputProbHPMap=np.zeros(12*nside**2) ras_graph,decs_graph=zip(*list(finalProbs.keys())) @@ -1163,6 +1234,7 @@ def make_posterior_map(self, nsToTest=[1,2,3,4,5], outputProbHPMapPlotting=copy.deepcopy(outputProbHPMap) outputProbHPMapPlotting[outputProbHPMapPlotting == 0] = np.nan + #Add Plotting if (plotting_location is not None): fig, ax=plt.subplots(1,1,figsize=(10,10)) ras_graph,decs_graph=zip(*list(finalProbs.keys())) @@ -1176,9 +1248,15 @@ def make_posterior_map(self, nsToTest=[1,2,3,4,5], #hp.visufunc.gnomview(map=outputProbHPMap,hold=True,title="Zoomed Posterior Skymap",rot=(self.src_ra*180/np.pi,self.src_dec*180/np.pi,0), # xsize=nPixZoom,ysize=nPixZoom,reso=reso,bgcolor="white",badcolor="white",margins=(.01,.01,0,0),notext=True,) - hp.visufunc.mollview(map=outputProbHPMapPlotting,hold=True,title="Posterior Skymap",bgcolor="white",badcolor="white",rot=(180,0,0)) + hp.visufunc.mollview(map=outputProbHPMapPlotting,hold=True,title="Posterior Skymap", + bgcolor="white",badcolor="white",rot=(180,0,0),) hp.graticule(verbose=True) - + ax=plt.gca() + text_objs = ax.findobj( + lambda artist: isinstance(artist, Text) + ) + for obj in text_objs: + obj.set_fontsize(30) plt.savefig(plotting_location+"PosteriorMap.png") plt.close() @@ -1186,31 +1264,152 @@ def make_posterior_map(self, nsToTest=[1,2,3,4,5], fig, ax=plt.subplots(1,1,figsize=(10,10)) plt.sca(ax) hp.visufunc.gnomview(map=outputProbHPMapPlotting,hold=True,title="Posterior Skymap",rot=(ramax*180/np.pi,decmax*180/np.pi,0), - xsize=nPixZoom,ysize=nPixZoom,reso=reso,bgcolor="white",badcolor="grey",margins=(.01,.01,0,0),notext=True,cbar=False) + xsize=nPixZoom,ysize=nPixZoom,reso=reso,bgcolor="white",badcolor="grey",margins=(.02,.02,.05,0),notext=True,cbar=False) hp.graticule(verbose=True) - plotting_utils.plot_labels(decmax,ramax, reso,nPix=nPixZoom,vertical_height_xlabel=.95) + plotting_utils.plot_labels(decmax,ramax, reso,nPix=nPixZoom,label_scale=.97) + ax=plt.gca() + text_objs = ax.findobj( + lambda artist: isinstance(artist, Text) + ) + for obj in text_objs: + obj.set_fontsize(30) plotting_utils.plot_color_bar(range=[0,np.nanmax(outputProbHPMapPlotting)], cmap="viridis", col_label=r"Posterior PDF", - offset=-35,labels=[0,"{0:.1e}".format(np.nanmax(outputProbHPMapPlotting)/2),"{0:.2e}".format(np.nanmax(outputProbHPMapPlotting))],loc=[0.90, 0.2, 0.03, 0.6]) + offset=-40,labels=[0,"{0:.1e}".format(np.nanmax(outputProbHPMapPlotting)/2),"{0:.2e}".format(np.nanmax(outputProbHPMapPlotting))],loc=[0.86, 0.2, 0.03, 0.6]) + theta, phi =plotting_utils.plot_contours([.9],outputProbHPMap) + hp.projplot(theta[0], phi[0], linewidth=2., c='k') + for i in range(1, len(theta)): + hp.projplot(theta[i], phi[i], linewidth=2., c='k', label=None) plt.savefig(plotting_location+"PosteriorMapZoomed.png") plt.close() + + #plot zoomed with events fig, ax=plt.subplots(1,1,figsize=(10,10)) plt.sca(ax) - hp.visufunc.gnomview(map=np.log10(outputProbHPMapPlotting),hold=True,title="Posterior Skymap (Log p)",rot=(ramax*180/np.pi,decmax*180/np.pi,0), - xsize=nPixZoom,ysize=nPixZoom,reso=reso,bgcolor="white",badcolor="grey",margins=(.01,.01,0,0),notext=True,cbar=False) + hp.visufunc.gnomview(map=outputProbHPMapPlotting,hold=True,title="Posterior Skymap",rot=(ramax*180/np.pi,decmax*180/np.pi,0), + xsize=nPixZoom,ysize=nPixZoom,reso=reso,bgcolor="white",badcolor="grey",margins=(.02,.02,.05,0),notext=True,cbar=False) + hp.graticule(verbose=True) + plotting_utils.plot_labels(decmax,ramax, reso,nPix=nPixZoom,label_scale=.97) + ax=plt.gca() + text_objs = ax.findobj( + lambda artist: isinstance(artist, Text) + ) + for obj in text_objs: + obj.set_fontsize(30) + plotting_utils.plot_color_bar(range=[0,np.nanmax(outputProbHPMapPlotting)], cmap="viridis", col_label=r"Posterior PDF", + offset=-40,labels=[0,"{0:.1e}".format(np.nanmax(outputProbHPMapPlotting)/2),"{0:.2e}".format(np.nanmax(outputProbHPMapPlotting))],loc=[0.86, 0.2, 0.03, 0.6]) + theta, phi =plotting_utils.plot_contours([.9],outputProbHPMap) + hp.projplot(theta[0], phi[0], linewidth=2., c='k') + for i in range(1, len(theta)): + hp.projplot(theta[i], phi[i], linewidth=2., c='k', label=None) + events = self.llh.exp + events = events[(events['time'] < self.stop) & (events['time'] > self.start)] + + plotting_utils.plot_events2(events['dec'], events['ra'], events['sigma']*self._angScale, + ramax,decmax, 2*6, sigma_scale=1.0, + constant_sigma=False, same_marker=True, energy_size=True, col = ["red"]*len(events["dec"]), + resolution=reso + ) + plotting_utils.plot_events2(np.radians(-36.61),np.radians(243.28), np.radians(.3), + ramax,decmax, 2*6, sigma_scale=1.0, + constant_sigma=False, same_marker=True, energy_size=True, col = ["yellow"]) #plot the cascade + plt.savefig(plotting_location+"PosteriorMapZoomedWithEvents.png") + plt.close() + + + + + + fig, ax=plt.subplots(1,1,figsize=(10,10)) + plt.sca(ax) + hp.visufunc.gnomview(map=np.log10(outputProbHPMapPlotting),hold=True,title="Posterior Skymap (Log p)",rot=(ramax*180/np.pi,decmax*180/np.pi,0), + xsize=nPixZoom,ysize=nPixZoom,reso=reso,bgcolor="white",badcolor="grey",margins=(.02,.02,.05,0),notext=True,cbar=False) hp.graticule(verbose=True) - plotting_utils.plot_labels(decmax,ramax, reso,nPix=nPixZoom,vertical_height_xlabel=.95) + plotting_utils.plot_labels(decmax,ramax, reso,nPix=nPixZoom,label_scale=.97) + ax=plt.gca() + text_objs = ax.findobj( + lambda artist: isinstance(artist, Text) + ) + for obj in text_objs: + obj.set_fontsize(30) plotting_utils.plot_color_bar(range=[0,np.nanmax(np.log10(outputProbHPMapPlotting))], cmap="viridis", col_label=r"Posterior PDF", - offset=-35,labels=[0,"{0:.1e}".format(np.nanmax(np.log10(outputProbHPMapPlotting))/2),"{0:.2e}".format(np.nanmax(np.log10(outputProbHPMapPlotting)))],loc=[0.90, 0.2, 0.03, 0.6]) + offset=-40,labels=[0,"{0:.1e}".format(np.nanmax(np.log10(outputProbHPMapPlotting))/2),"{0:.2e}".format(np.nanmax(np.log10(outputProbHPMapPlotting)))],loc=[0.86, 0.2, 0.03, 0.6]) + theta, phi =plotting_utils.plot_contours([.9],outputProbHPMap) + hp.projplot(theta[0], phi[0], linewidth=2., c='k') + for i in range(1, len(theta)): + hp.projplot(theta[i], phi[i], linewidth=2., c='k', label=None) plt.savefig(plotting_location+"PosteriorMapZoomedLog.png") plt.close() + + #Zoom in much further + reso=.1 + sizeZoomInDegs=2 + nPixZoom=sizeZoomInDegs/reso*60 + + fig, ax=plt.subplots(1,1,figsize=(10,10)) + plt.sca(ax) + hp.visufunc.gnomview(map=outputProbHPMapPlotting,hold=True,title="Posterior Skymap",rot=(ramax*180/np.pi,decmax*180/np.pi,0), + xsize=nPixZoom,ysize=nPixZoom,reso=reso,bgcolor="white",badcolor="grey",margins=(.02,.02,.05,0),notext=True,cbar=False) + hp.graticule(verbose=True) + plotting_utils.plot_labels(decmax,ramax, reso,nPix=nPixZoom,label_scale=.97) + ax=plt.gca() + text_objs = ax.findobj( + lambda artist: isinstance(artist, Text) + ) + for obj in text_objs: + obj.set_fontsize(30) + plotting_utils.plot_color_bar(range=[0,np.nanmax(outputProbHPMapPlotting)], cmap="viridis", col_label=r"Posterior PDF", + offset=-40,labels=[0,"{0:.1e}".format(np.nanmax(outputProbHPMapPlotting)/2),"{0:.2e}".format(np.nanmax(outputProbHPMapPlotting))],loc=[0.86, 0.2, 0.03, 0.6]) + theta, phi =plotting_utils.plot_contours([.9],outputProbHPMap) + hp.projplot(theta[0], phi[0], linewidth=2., c='k') + for i in range(1, len(theta)): + hp.projplot(theta[i], phi[i], linewidth=2., c='k', label=None) + events = self.llh.exp + events = events[(events['time'] < self.stop) & (events['time'] > self.start)] + + plotting_utils.plot_events2(events['dec'], events['ra'], events['sigma']*self._angScale, + ramax,decmax, 2*6, sigma_scale=1.0, + constant_sigma=False, same_marker=True, energy_size=True, col = ["red"]*len(events["dec"]), + resolution=reso + ) + plotting_utils.plot_events2(np.radians(-36.61),np.radians(243.28), np.radians(.3), + ramax,decmax, 2*6, sigma_scale=1.0, + constant_sigma=False, same_marker=True, energy_size=True, col = ["yellow"]) #plot the cascade + plt.savefig(plotting_location+"PosteriorMapZoomed1DegWithEvents.png") + plt.close() + + #plot zoomed + fig, ax=plt.subplots(1,1,figsize=(10,10)) + plt.sca(ax) + hp.visufunc.gnomview(map=outputProbHPMapPlotting,hold=True,title="Posterior Skymap",rot=(ramax*180/np.pi,decmax*180/np.pi,0), + xsize=nPixZoom,ysize=nPixZoom,reso=reso,bgcolor="white",badcolor="grey",margins=(.02,.02,.05,0),notext=True,cbar=False) + hp.graticule(verbose=True) + plotting_utils.plot_labels(decmax,ramax, reso,nPix=nPixZoom,label_scale=.97) + ax=plt.gca() + text_objs = ax.findobj( + lambda artist: isinstance(artist, Text) + ) + for obj in text_objs: + obj.set_fontsize(30) + plotting_utils.plot_color_bar(range=[0,np.nanmax(outputProbHPMapPlotting)], cmap="viridis", col_label=r"Posterior PDF", + offset=-40,labels=[0,"{0:.1e}".format(np.nanmax(outputProbHPMapPlotting)/2),"{0:.2e}".format(np.nanmax(outputProbHPMapPlotting))],loc=[0.86, 0.2, 0.03, 0.6]) + theta, phi =plotting_utils.plot_contours([.9],outputProbHPMap) + hp.projplot(theta[0], phi[0], linewidth=2., c='k') + for i in range(1, len(theta)): + hp.projplot(theta[i], phi[i], linewidth=2., c='k', label=None) + plt.savefig(plotting_location+"PosteriorMapZoomed1Deg.png") + plt.close() + + + + #plot ns vs probability for some specific points, could use a cleaner method of choosing points threshold_prob_for_plots=list(sorted([k for k, v in Counter(list(finalProbs.values())).items() if v == 1],reverse=True)) #threshold_prob_for_plots=[k for k in threshold_prob_for_plots if k>0] - counterJump=int(len(threshold_prob_for_plots)/30) + counterJump=max(int(len(threshold_prob_for_plots)/30),5) print(len(threshold_prob_for_plots),counterJump) psToPlot=[] @@ -1231,6 +1430,7 @@ def make_posterior_map(self, nsToTest=[1,2,3,4,5], total_val=finalProbs[(ra,dec)] temp_val=val[(val["ra"]==ra)&(val["dec"]==dec)] nss_temp=temp_val["nsignal"] + #gammas_temp=temp_val["gamma"] TS_temp=temp_val["TS_spatial_prior_0"]/2 @@ -1242,8 +1442,8 @@ def make_posterior_map(self, nsToTest=[1,2,3,4,5], cmap.set_array([]) cbar=plt.colorbar(cmap) cbar.set_label('Pixel Probability') - plt.xlabel("ns",fontsize=15) - plt.ylabel("probability",fontsize=15) + plt.xlabel("ns",fontsize=30) + plt.ylabel("probability",fontsize=30) plt.savefig(plotting_location+"nsProfilePlot.png") plt.close() diff --git a/fast_response/latex/posterior_report_general.tex b/fast_response/latex/posterior_report_general.tex new file mode 100644 index 0000000..f76e224 --- /dev/null +++ b/fast_response/latex/posterior_report_general.tex @@ -0,0 +1,70 @@ +\documentclass[titlepage]{article} +\usepackage{tikz} +\usepackage{xcolor} +\usepackage{graphicx} +\usepackage{grffile} +\usepackage{longtable} +\usepackage[margin=1.0in]{geometry} + +\setlength{\parindent}{0cm} +\include{r_posterior} + +\newcommand{\degree}{$^{\circ}$} +\begin{document} + +\begin{titlepage} + \centering + \vspace{4cm} + {\huge\bfseries IceCube Fast-Response \\ Posterior Supplemental Report\par} + \vspace{1cm} + {\LARGE For Internal Use Only\par} + \vfill + {\Large Source Name: \\ \itshape\sourcename\par} + \vspace{0.5cm} + {\Large Observation Date(s):\\ \obsdate \par} + \vfill + \vspace{1cm} + {\Large Report Generated On:\\ \reportdate \par} +\end{titlepage} + + +\pagebreak +\section{Results} + +{ + \centering + {\Large Prior Zoomed} + + \includegraphics[width=0.9\textwidth]{\PriorSkymap} + + +} +\section{Posterior Maps} + +{ + \centering + {\Large Posterior Skymap} + + \includegraphics[width=0.9\textwidth]{\PosteriorSkymap} + \includegraphics[width=0.9\textwidth]{\PosteriorSkymapZoom} + \includegraphics[width=0.9\textwidth]{\PosteriorSkymapZoomWNeutrinos} + + +} +\pagebreak + +{ + \centering + {\Large $N_s$ Profile Plot} + + \includegraphics[width=0.9\textwidth]{\NsProfilePlot} + + +} + + + + +\vfill + +\end{document} diff --git a/fast_response/plotting_utils.py b/fast_response/plotting_utils.py index f383d12..99c0b57 100644 --- a/fast_response/plotting_utils.py +++ b/fast_response/plotting_utils.py @@ -70,7 +70,7 @@ def plot_color_bar(labels=[0.,2.,4.,6.], col_label=r"IceCube Event Time", range= cb.set_ticklabels(labels) cb.update_ticks() -def plot_labels(src_dec, src_ra, reso,nPix=200,vertical_height_xlabel=1): +def plot_labels(src_dec, src_ra, reso,nPix=200,label_scale=1): """ Add labels to healpy zoom @@ -103,14 +103,45 @@ def plot_labels(src_dec, src_ra, reso,nPix=200,vertical_height_xlabel=1): plt.text(np.radians(-reso),np.radians(-1.75*reso), r"%.1f$^{\circ}$"%(reso+np.degrees(src_ra)), horizontalalignment='center', verticalalignment='top', fontsize=fontsize) - plt.text(-1*np.radians(2.35*reso), np.radians(0), r"declination", + plt.text(-1*np.radians(2.35*reso*label_scale), np.radians(0), r"declination", ha='center', va='center', rotation=90, fontsize=fontsize) - plt.text(np.radians(0), np.radians(-2.05*reso*vertical_height_xlabel), r"right ascension", + plt.text(np.radians(0), np.radians(-2.05*reso*label_scale), r"right ascension", ha='center', va='center', fontsize=fontsize) +def plot_events2(dec,ra,sigmas, src_ra, src_dec, reso, sigma_scale=5., col = 'k', constant_sigma=False, + same_marker=False, energy_size=False, with_mark=True, with_dash=False, + label='',resolution=0): + def compute_ang_err(ra,dec,sigma): + dec = np.pi/2 - dec + sigma = np.rad2deg(sigma) + delta, step, bins = 0, 0, 0 + delta= sigma/180.0*np.pi + step = 1./np.sin(delta)/20. + bins = int(360./step) + Theta = np.zeros(bins+1, dtype=np.double) + Phi = np.zeros(bins+1, dtype=np.double) + # define the contour + for j in range(0,bins): + phi = j*step/180.*np.pi + vx = np.cos(phi)*np.sin(ra)*np.sin(delta) + np.cos(ra)*(np.cos(delta)*np.sin(dec) + np.cos(dec)*np.sin(delta)*np.sin(phi)) + vy = np.cos(delta)*np.sin(dec)*np.sin(ra) + np.sin(delta)*(-np.cos(ra)*np.cos(phi) + np.cos(dec)*np.sin(ra)*np.sin(phi)) + vz = np.cos(dec)*np.cos(delta) - np.sin(dec)*np.sin(delta)*np.sin(phi) + Theta[j], Phi[j] = hp.vec2ang(np.array([vx, vy, vz])) + + Theta[bins] = Theta[0] + Phi[bins] = Phi[0] + + return Theta, Phi + + for i in range(len(ra)): + ev_contour = compute_ang_err(ra[i],dec[i],sigmas[i]) + hp.projplot(ev_contour[0], ev_contour[1], linewidth=1.75, color=col[0], + linestyle="solid",coord='C') + + def plot_events(dec, ra, sigmas, src_ra, src_dec, reso, sigma_scale=5., col = 'k', constant_sigma=False, same_marker=False, energy_size=False, with_mark=True, with_dash=False, - label=''): + label='',resolution=0.025): """ Adds events to a healpy zoom plot. Events are expected to be from self.llh.exp @@ -149,7 +180,7 @@ def plot_events(dec, ra, sigmas, src_ra, src_dec, reso, sigma_scale=5., col = 'k if sigma_scale is not None: sigma = np.degrees(sigmas)/sigma_scale - sizes = 5200*sigma**2 + sizes = 5200*np.power(resolution/.025,.5)*sigma**2 if constant_sigma: sizes = 20*np.ones_like(sizes) if with_dash: diff --git a/fast_response/reports/PosteriorReportGenerator.py b/fast_response/reports/PosteriorReportGenerator.py new file mode 100644 index 0000000..9a3dbc2 --- /dev/null +++ b/fast_response/reports/PosteriorReportGenerator.py @@ -0,0 +1,463 @@ +'''Class for generating reports for Fast Reponse and GW + followup analyses. A lot of code has been taken from Fast Reponse + Analysis. + + Author: Kevin Meagher, Raamis Hussain & Alex Pizzuto + Date: Mar 27, 2019 + Modified: April, 2019 +''' + +import numpy as np +import pandas as pd +import os, json, datetime +import logging as log + +use_urllib2 = True +try: + import urllib2, urllib +except: + use_urllib2 = False + import urllib.request, urllib.parse, urllib.error + +import icecube.realtime_tools.live +import subprocess +from os.path import expanduser + +from icecube import icetray +import skylab +from skylab.datasets import Datasets +from astropy.time import Time, TimeDelta +from fast_response.make_ontime_plots import make_rate_plots +import fast_response + +log.basicConfig(level=log.ERROR) +mpl_logger = log.getLogger('matplotlib') +mpl_logger.setLevel(log.ERROR) + +class PosteriorReportGenerator(object): + r""" + Generate a posterior report from a FastResponseAnalysis object + """ + _figure_list = [] + + def __init__(self, analysis): + self.analysis = analysis + if self.analysis.skymap is None: + self.source_type = 'PS' + else: + self.source_type = 'skymap' + + source = {} + source['trigger_iso'] = Time(analysis.centertime, format='mjd').iso + source['start_iso'] = Time(analysis.start, format='mjd').iso + source['stop_iso'] = Time(analysis.stop, format='mjd').iso + source['realtime'] = (analysis.stop - analysis.start) * 86400. + source['start_mjd'] = analysis.start + source['stop_mjd'] = analysis.stop + source['trigger_mjd'] = analysis.centertime + self.source = source + + if self.source["stop_mjd"] > 58309.74: + self.stream = 'neutrino' + elif self.source["stop_mjd"] > 57891.17: + self.stream = 'neutrino17' + else: + self.stream = 'neutrino16' + + self.analysisid = analysis.analysisid + self.dirname = analysis.analysispath + self.time_window=( + Time(self.source['start_mjd'], format='mjd'), + Time(self.source['stop_mjd'], format='mjd')) + self.trigger = Time(self.source['trigger_mjd'], format='mjd') + + def get_report_source(self): + """ Find base path of the general report TeX document + Looks for the directory location fast_response, returns + fast_response/latex/report_general.tex + """ + base = os.path.dirname(fast_response.__file__) + return os.path.join(base, 'latex/posterior_report_general.tex') + + def write_table(self, file, name, header, table, prefix=""): + """ + Write a table in latex format for the generated reports + + Parameters + ----------- + file: I/O file object + Opened file that is being written to + name: str + Name of the table + header: list + Headers for each column (can be empty if none) + table: array of tuples + set of values to put in the table. tuples give individual rows + prefix: str + prefix to latex command, if needed (default "") + """ + + cols = max([len(r) for r in table if r is not None]+[len(header)]) + + x=''.join(' & '.join(str(x).strip() for x in row)+r'\\'+'\n' for row in table if row is not None) + + file.write( + prefix+ + r"\newcommand{"+"\\"+name+"}{"+"\n"+ + r"\begin{longtable}{" + 'l' * cols + r"}"+'\n'+ + r"\hline"+'\n'+ + ' & '.join(str(x).strip() for x in header)+ + (r"\\ \hline"+'\n' if header else "") + + ''.join(' & '.join(str(x).strip() for x in row)+r'\\'+'\n' if row is not None else r"\hline"+'\n' for row in table)+ + ((r"\hline"+'\n') if table and table[-1] is not None else "") + + r"\end{longtable}"+ + "}\n") + + def query_db_runs(self, time_window): + """ Queries IceCube Live for the ontime window, + plus 8 hours on either side (or as much as possible, if in realtime) + + Parameters + ----------- + time_window: tuple of astropy.time Time objects + (start, stop) for the analysis, in mjd + + Returns + ---------- + dict: + table of runs loaded, and status of the runs + """ + + run_url = 'https://live.icecube.wisc.edu/run_info/' + with open('/home/jthwaites/private/tokens/auth.txt') as f: + user = f.readline().rstrip('\n') + pasw = f.readline().rstrip('\n') + run_query = {'user':user, 'pass':pasw, + 'start':(Time(time_window[0],precision=0) + - TimeDelta(8*3600,format='sec')).iso, + 'stop': (Time(time_window[1],precision=0) + + TimeDelta(8*3600,format='sec')).iso} + now = Time(datetime.datetime.now()+datetime.timedelta(hours=6),scale='utc',precision=0) + if use_urllib2: + run_table = json.loads(urllib2.urlopen( + urllib2.Request(run_url, urllib.urlencode(run_query)), + timeout=500).read()) + else: + req_data = urllib.parse.urlencode(run_query).encode("utf-8") + req = urllib.request.Request(run_url) + with urllib.request.urlopen(req, data=req_data, timeout=500) as fi: + run_table = json.load(fi) + + #add duration to run table + for run in run_table: + runstart = Time(run['start'],format='iso',scale='utc') + if run['stop'] is None: + runstop = now + run['stop']=runstop.iso + else: + runstop = Time(run['stop'],scale='utc',precision=0) + + dt=int((runstop-runstart).sec) + hours = int(dt)//3600 + minutes = int(dt)//60-hours*60 + seconds = dt%60 + + livetime = (min(time_window[1],runstop)- + max(time_window[0],runstart)).sec + run['livetime'] = livetime if livetime > 0 else 0 + if (run['livetime'] > 86400. * 2) or (hours > 100): + run['livetime'] = 0 + run['stop'] = run['start'] + dt=0 + hours = 0 + minutes = 0 + seconds = 0 + + run['duration_sec'] = dt + run['duration'] = "{:01d}:{:02d}:{:02d}".format(hours,minutes,seconds) + run['OK'] = ("OK" if run['status'] in ['SUCCESS'] + and run['lightmode'] in ['dark'] + and run['filter_mode'] in ['PhysicsFiltering'] + and run['run_mode'] in ['PhysicsTrig'] + else "NotOK" + ) + + return run_table + + @staticmethod + def ontime_table(query_dict): + """Makes a DataFrame of events, queried from IceCube Live + + Parameters + ----------- + query_dict: dict + Queried events in a particular time window + + Returns + ---------- + pandas DataFrame: + Event parameters for the given events + """ + newdict=[] + for event in query_dict: + newevent = event['value']['data'] + for key,val in list(newevent['reco']['splinempe'].items()): + newevent['splinempe_'+key]=val + if Time(newevent['eventtime'],scale='utc',format='iso') >= Time("2018-07-10 17:52:03.34", format='iso',scale='utc'): + newevent['muex'] = newevent['reco']['energy']['mpe_muex'] + del newevent['reco'] + newdict.append(newevent) + events = pd.DataFrame(newdict) + + if len(events): + t=Time(list(events['eventtime']),scale='utc',format='iso') + events['t']=t.datetime + events['mjd'] = t.mjd + return events + + def make_coinc_events_table(self, f): + """Make a table of ontime events, in LaTeX format + + Parameters + ----------- + f: I/O file object + Opened file that is being written to + """ + event_table = [] + if self.analysis.coincident_events is not None and self.analysis.coincident_events != []: + if self.analysis.skymap is None: + for event in self.analysis.coincident_events: + event_table+=[ + ("Run:Event",'{}:{}'.format(event['run'], event['event'])), + ("Time","{}".format( + event['time'])), + (r'$\alpha$, $\delta$',"{:3.2f}\degree, {:+3.2f}\degree" + .format(np.rad2deg(event['ra']), np.rad2deg(event['dec']))), + ("Angular Uncertainty (90\%)","{:3.2f}\degree" + .format(np.rad2deg(event["sigma"]*2.145966))), + ("Distance from Source", "{:2.2f}\degree".format(event['delta_psi']*180. / np.pi)), + ("Reconstructed Energy (GeV)","{:2.2e}" + .format(10**event['logE'])), + ("Spatial Weight", "{:.2f}".format(event['spatial_w'])), + ("Energy Weight", "{:.2f}".format(event['energy_w'])), + None, + ] + + self.write_table(f, "event", [], event_table) + else: + if 'pvalue' in self.analysis.coincident_events[0].keys(): + with_p = True + if len(self.analysis.coincident_events)<100: + event_table+=[ + ('Event','$t_{\\nu}$-$t_{trigger}$ (s)','RA','Dec', + '\sigma (90\%)','E_{reco} (GeV)', + 'p-value','In 90\% Contour')] + else: + event_table+=[ + ('$t_{\\nu}$-$t_{trigger}$ (s)','RA','Dec', + '\sigma (90\%)','E_{reco} (GeV)', + 'p-value','In 90\% Contour')] + else: + with_p = False + event_table+=[ + ('$t_{\\nu}$-$t_{trigger}$ (s)','RA','Dec', + '\sigma (90\%)','E_{reco} (GeV)', + 'In 90\% Contour')] + i = 0 + for event in self.analysis.coincident_events: + if with_p: + if len(self.analysis.coincident_events)>100: + if event['in_contour']: + event_table+=[ + ('{:.0f}'.format((event['time']-self.source['trigger_mjd'])*86400.), + "{:3.2f}\degree".format(np.rad2deg(event['ra'])), + '{:3.2f}\degree'.format(np.rad2deg(event['dec'])), + "{:3.2f}\degree".format(np.rad2deg(event["sigma"]*self.analysis._angScale)), + "{:.2E}".format(10**event['logE']), + '{:.4f}'.format(event['pvalue']), + str(bool(event['in_contour'])) + )] + else: + event_table+=[ + ('{}'.format(i+1), + '{:.0f}'.format((event['time']-self.source['trigger_mjd'])*86400.), + "{:3.2f}\degree".format(np.rad2deg(event['ra'])), + '{:3.2f}\degree'.format(np.rad2deg(event['dec'])), + "{:3.2f}\degree".format(np.rad2deg(event["sigma"]*self.analysis._angScale)), + "{:.2E}".format(10**event['logE']), + '{:.4f}'.format(event['pvalue']), + str(bool(event['in_contour'])) + )] + i+=1 + else: + if len(self.analysis.coincident_events)>100: + if event['in_contour']: + event_table+=[ + ('{:.0f}'.format((event['time']-self.source['trigger_mjd'])*86400.), + "{:3.2f}\degree".format(np.rad2deg(event['ra'])), + '{:3.2f}\degree'.format(np.rad2deg(event['dec'])), + "{:3.2f}\degree".format(np.rad2deg(event["sigma"]*self.analysis._angScale)), + "{:.2E}".format(10**event['logE']), + str(bool(event['in_contour'])) + )] + else: + event_table+=[ + ('{:.0f}'.format((event['time']-self.source['trigger_mjd'])*86400.), + "{:3.2f}\degree".format(np.rad2deg(event['ra'])), + '{:3.2f}\degree'.format(np.rad2deg(event['dec'])), + "{:3.2f}\degree".format(np.rad2deg(event["sigma"]*self.analysis._angScale)), + "{:.2E}".format(10**event['logE']), + str(bool(event['in_contour'])) + )] + self.write_table(f,"event",[],event_table) + else: + f.write(r"\newcommand{\event}{[None]}") + + def make_new_command(self, f, name, command): + """Define a new command, in LaTeX, for the report + + Parameters + ----------- + f: I/O file object + Opened file that is being written to + name: str + Name of the command being defined + command: str + Value of the command being defined + """ + f.write(r"\newcommand{"+"\\"+name+"}{"+ + command+"}\n") + + def generate_posterior_report(self): + """Generate the report! This function does most of the heavy lifting. + Starts making the source LaTeX file (named r.tex), queries runs, + makes detector status plots, includes all plots for the analysis, + makes all tables for analysis parameters, runs, events, etc. + """ + + report_fname = os.path.join(self.dirname, "r_posterior.tex") + + self.run_table = self.query_db_runs(self.time_window) + now = datetime.datetime.utcnow() + self.ontime = {} + self.ontime['type'] = 'database' + self.ontime['stream'] = self.stream + self.ontime['runkey'] = 'run_id' + self.ontime['time_start'] = Time( + self.run_table[0]['start'], + format='iso', + scale='utc', + precision=0).iso + self.ontime['time_stop'] =Time( + self.run_table[-1]['stop'], + format='iso', + scale='utc', + precision=0).iso + self.ontime['time_query']=now.strftime('%Y-%m-%d %H:%M:%S') + + self.query_events=icecube.realtime_tools.live.get_events( + self.ontime['stream'], + self.ontime['time_start'], + self.ontime['time_stop']) + + + + s = self.source + + dataset = Datasets[self.analysis._dataset] + # Make rate plots and put them in analysis directory + + with open(report_fname, 'wt') as f: + + d1 = s["start_iso"].split()[0] + d2 = s["stop_iso"].split()[0] + + if d1 == d2: + obsdate = d1 + else: + obsdate = d1+' -- '+d2 + + for name, command in [('sourcename', self.analysis.name), + ('analysisid', self.analysisid), + ('reportdate', datetime.date.today().strftime("%Y-%m-%d")), + ('obsdate', obsdate)]: + self.make_new_command(f, name, command) + print(self.analysisid) + for plot_name, plot_path in [ + ('PriorSkymap', "PriorSkymap.png"), + ('PosteriorSkymap', "PosteriorMap.png"), + ('PosteriorSkymapZoom',"PosteriorMapZoomed.png"), + ('NsProfilePlot',"nsProfilePlot.png"), + ("PosteriorSkymapZoomWNeutrinos","PosteriorMapZoomedWithEvents.png") + ]: + if os.path.isfile(self.dirname + '/' + plot_path): + f.write( + r"\newcommand{"+"\\"+plot_name+"}{"+ self.dirname+"/"+ + plot_path + + "}\n" + ) + + + + def make_pdf(self): + """Makes the PDF compiled version of the LaTeX report + using pdflatex, and saves to the analysis folder + """ + # symlink main report tex file + reportfname = self.analysisid + "_posterior_report.tex" + reportpath = os.path.join(self.dirname, reportfname) + reportsrc = self.get_report_source() + + if os.path.exists(reportpath): + os.unlink(reportpath) + + os.symlink(reportsrc, reportpath) + + # get environment variables + env = dict(os.environ) + subprocess.call( + ['pdflatex', '-interaction=batchmode', + '-output-directory=%s' % self.dirname, + self.dirname + '/' + self.analysisid + "_posterior_report.tex"], + #cwd=self.dirname, + env = env, + ) + +class FastResponsePosteriorReport(PosteriorReportGenerator): + """Generates a general FastResponse Report for a point source. + Includes additional plots of the 1D LLH ns scan + """ + _figure_list = [('nsscan', 'llh_ns_scan.png')] + + def __init__(self, analysis): + super().__init__(analysis) + + def generate_posterior_report(self): + super().generate_posterior_report() + +class GravitationalWavePosteriorReport(PosteriorReportGenerator): + """Generates a FastResponse Report for a GW. + Includes additional plot of skymap PDF versus dec, + with PS sensitivity overlaid + """ + _figure_list = [('decPDF', 'decPDF.png')] + + def get_report_source(self): + """Find base path of the general report TeX document. + Looks for the directory location fast_response, + returns fast_response/latex/posterior_report_gw.tex + """ + base = os.path.dirname(fast_response.__file__) + return os.path.join(base, 'latex/posterior_report_gw.tex') + +class AlertPosteriorReport(PosteriorReportGenerator): + """Generates a FastResponse Report for an Alert event.""" + _figure_list = [] + + def get_report_source(self): + """Find base path of the general report TeX document. + Looks for the directory location fast_response, + returns fast_response/latex/posterior_report_alert.tex + """ + base = os.path.dirname(fast_response.__file__) + return os.path.join(base, 'latex/posterior_report_alert.tex') \ No newline at end of file diff --git a/fast_response/reports/__init__.py b/fast_response/reports/__init__.py index 0cdf2eb..dc5dbb3 100644 --- a/fast_response/reports/__init__.py +++ b/fast_response/reports/__init__.py @@ -2,4 +2,7 @@ from .ReportGenerator import AlertReport from .ReportGenerator import GravitationalWaveReport from .ReportGenerator import FastResponseReport +from .PosteriorReportGenerator import \ + FastResponsePosteriorReport,GravitationalWavePosteriorReport,PosteriorReportGenerator,AlertPosteriorReport + from .report_utils import valid_location, get_times, interptime \ No newline at end of file diff --git a/fast_response/scripts/test_1.sh b/fast_response/scripts/test_1.sh new file mode 100644 index 0000000..64e3b73 --- /dev/null +++ b/fast_response/scripts/test_1.sh @@ -0,0 +1,4 @@ +python test_posterior_cascade.py --skymap=https://roc.icecube.wisc.edu/public/hese_cascades/hese_60505_run00139644.evt000063900267.fits --time=60505.623363 --alert_id=139644:63900267 + +python test_posterior_code.py --skymap=https://gracedb.ligo.org/apiweb/superevents/S240716b/files/Bilby.fits.gz --time=60000.00 \ + --name="Posterior_Test_2inj" --tw=1000 --n_inj=2 diff --git a/requirements.txt b/requirements.txt index 8cedf0d..45555b3 100644 --- a/requirements.txt +++ b/requirements.txt @@ -29,4 +29,4 @@ subprocess32==3.5.4 zmq==0.0.0 py27hash==1.0.2 pygcn==1.1.2 -wheel==0.37.1 +wheel==0.37.1 \ No newline at end of file From 4a6c43a2daf897e50d5f2d5b43c7a4a46b90a945 Mon Sep 17 00:00:00 2001 From: Sam Hori Date: Fri, 30 May 2025 13:17:54 -0500 Subject: [PATCH 6/7] Printing some debug messages temporarily --- fast_response/FastResponseAnalysis.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/fast_response/FastResponseAnalysis.py b/fast_response/FastResponseAnalysis.py index c9f6e79..48711ca 100644 --- a/fast_response/FastResponseAnalysis.py +++ b/fast_response/FastResponseAnalysis.py @@ -1310,7 +1310,7 @@ def make_posterior_map(self, fluxToTest=[1,2,3,4,5], constant_sigma=False, same_marker=True, energy_size=True, col = ["red"]*len(events["dec"]), resolution=reso ) - plotting_utils.plot_events2(np.radians(-36.61),np.radians(243.28), np.radians(.3), + plotting_utils.plot_events2([np.radians(-36.61)],[np.radians(243.28)], [np.radians(.3)], ramax,decmax, 2*6, sigma_scale=1.0, constant_sigma=False, same_marker=True, energy_size=True, col = ["yellow"]) #plot the cascade plt.savefig(plotting_location+"PosteriorMapZoomedWithEvents.png") @@ -1375,7 +1375,7 @@ def make_posterior_map(self, fluxToTest=[1,2,3,4,5], constant_sigma=False, same_marker=True, energy_size=True, col = ["red"]*len(events["dec"]), resolution=reso ) - plotting_utils.plot_events2(np.radians(-36.61),np.radians(243.28), np.radians(.3), + plotting_utils.plot_events2([np.radians(-36.61)],[np.radians(243.28)], [np.radians(.3)], ramax,decmax, 2*6, sigma_scale=1.0, constant_sigma=False, same_marker=True, energy_size=True, col = ["yellow"]) #plot the cascade plt.savefig(plotting_location+"PosteriorMapZoomed1DegWithEvents.png") @@ -1396,10 +1396,12 @@ def make_posterior_map(self, fluxToTest=[1,2,3,4,5], obj.set_fontsize(30) plotting_utils.plot_color_bar(range=[0,np.nanmax(outputProbHPMapPlotting)], cmap="viridis", col_label=r"Posterior PDF", offset=-40,labels=[0,"{0:.1e}".format(np.nanmax(outputProbHPMapPlotting)/2),"{0:.2e}".format(np.nanmax(outputProbHPMapPlotting))],loc=[0.86, 0.2, 0.03, 0.6]) - theta, phi =plotting_utils.plot_contours([.9],outputProbHPMap) + print("90 contour",theta,phi) hp.projplot(theta[0], phi[0], linewidth=2., c='k') for i in range(1, len(theta)): hp.projplot(theta[i], phi[i], linewidth=2., c='k', label=None) + theta, phi =plotting_utils.plot_contours([.9],outputProbHPMap) + plt.savefig(plotting_location+"PosteriorMapZoomed1Deg.png") plt.close() From 1936dbb09294c26438ecd2c02840330e11c9cff7 Mon Sep 17 00:00:00 2001 From: Sam Hori Date: Mon, 16 Jun 2025 20:25:59 -0500 Subject: [PATCH 7/7] Working on code cleanup for posterior work --- fast_response/FastResponseAnalysis.py | 67 +++++++------------ .../latex/posterior_report_general.tex | 4 ++ fast_response/plotting_utils.py | 11 ++- 3 files changed, 32 insertions(+), 50 deletions(-) diff --git a/fast_response/FastResponseAnalysis.py b/fast_response/FastResponseAnalysis.py index 48711ca..65d0f96 100644 --- a/fast_response/FastResponseAnalysis.py +++ b/fast_response/FastResponseAnalysis.py @@ -1029,12 +1029,9 @@ def make_prior_map(self, plotting_location=None, format=True ): - #TODO: Document """ Generates a zoomed in skymap around the best fit location from the maximum likelihood analysis """ - - #Format the skymap if told to format (usually formatted at construction so shouldn't change anything) if(format): skymap = self.format_skymap(self.skymap) @@ -1049,9 +1046,7 @@ def make_prior_map(self, nPixZoom=sizeZoomInDegs/reso*60 hp.visufunc.gnomview(map=np.log10(skymap),hold=True,title="Prior Skymap",rot=(np.degrees(self.skymap_fit_ra),np.degrees(self.skymap_fit_dec),0), xsize=nPixZoom,ysize=nPixZoom,reso=reso,bgcolor="white",badcolor="white",margins=(.02,.02,.05,0),notext=True,cbar=False) - #hp.visufunc.cartview(map=np.log10(skymap),hold=True,title="Prior Skymap",rot=(np.degrees(self.skymap_fit_ra),np.degrees(self.skymap_fit_dec),0), - # lonra=[np.degrees(self.skymap_fit_ra)-5,np.degrees(self.skymap_fit_ra)+5],latra=[np.degrees(self.skymap_fit_dec)-5,np.degrees(self.skymap_fit_dec)+5],bgcolor="white",badcolor="white",margins=(.01,.01,0,0),notext=True,) - + hp.graticule(verbose=True) plotting_utils.plot_labels(self.skymap_fit_dec,self.skymap_fit_ra, reso,nPix=nPixZoom,label_scale=.97) @@ -1072,8 +1067,9 @@ def make_prior_map(self, name+="Unformatted" plt.savefig(plotting_location+name+".png") - def _make_posterior_fits(self,pos_skymap): + def _make_posterior_fits(self): #internal function to be called from make_posterior_skymap + #saves the skymaps as a fits file extra_header = [('index', self.index), #('energy_range', ), #('ns_range', ), @@ -1084,10 +1080,11 @@ def _make_posterior_fits(self,pos_skymap): ] #startmjd, stopmjd, - hp.write_map(self.analysispath + '/' +'posterior_info.fits',m=pos_skymap, + hp.write_map(self.analysispath + '/' +'posterior_info.fits',m=self.PosteriorSkymap, extra_header=extra_header,partial=True ) + def fluxToAverageNs(self,flux): if(self.inj==None): @@ -1116,14 +1113,13 @@ def make_posterior_map(self, fluxToTest=[1,2,3,4,5], #flux/ns to test: given as list #Give prior as a function f(ns, gamma, ra, dec) #As currently implemented,the llh will have gamma must be constant set by the self.index - #Prior flat in ns and spatially by default. May make sense to implement + #Prior flat in ns and spatially by default. #plotting_location is the directory where plots should be output #custom events will be injected #useFlux determines whether fluxToTest is treated as ns or flux - #TODO: flux/ns priors t1 = time.time() val=None skymap = self.skymap @@ -1134,17 +1130,7 @@ def make_posterior_map(self, fluxToTest=[1,2,3,4,5], containment = self._containment,allow_neg=self._allow_neg) nside=self.nside - - #The full sky scan with no fixing (maybe better to save this during the unblind TS stage) - #Can be removed for speed - full_scan_val = self.llh.scan(0.0,0.0, scramble = False,spatial_prior=spatial_prior, - time_mask = [self.duration/2., self.centertime], - pixel_scan=[nside, self._pixel_scan_nsigma], - fixed=['nsignal'], - inject=custom_events - ) - - decs_list=full_scan_val["dec"] + decs_list = hp.pix2ang(nside, np.arange(hp.nside2npix(nside))) test_flux=1e-17 @@ -1158,13 +1144,11 @@ def make_posterior_map(self, fluxToTest=[1,2,3,4,5], for fluxT in mean_flux_arr: - #print('testing flux=',fluxT) temp_val = self.llh.scan(0.0,0.0, scramble = False,spatial_prior=spatial_prior, time_mask = [self.duration/2., self.centertime], pixel_scan=[nside, self._pixel_scan_nsigma], fixed=['nsignal'], nsignal=np.array([fluxT*conversion for conversion in nsPerFluxArr],dtype=float), - #nsignal=np.array([fluxToTest*nsPerFlux for nsPerFlux in nsPerFluxArr],dtype=float), inject=custom_events ) @@ -1172,24 +1156,19 @@ def make_posterior_map(self, fluxToTest=[1,2,3,4,5], val=temp_val else: val=np.hstack((val,temp_val)) - #print(val.dtype.names) ras=val["ra"] decs=val["dec"] nss=val["nsignal"] - + + #Placeholder: not required in fixed index analyses, but future improvements could expand this if "gamma" in val.dtype.names: - #TODO: This is available in some, but not all, analyses -- GW vs Internal alerts - # should always be self.index (unless things get) - # maybe "spectrum" instead - gammas=val["gamma"] - #gammas=val["gamma"] TS_spatial_prior_0=val["TS_spatial_prior_0"] idx_max=np.argmax(TS_spatial_prior_0) ramax,decmax=ras[idx_max],decs[idx_max] - logProbs=TS_spatial_prior_0/2 #not normalized so technically C*ln(P) + logProbs=TS_spatial_prior_0/2 finalProbs={} @@ -1211,17 +1190,16 @@ def make_posterior_map(self, fluxToTest=[1,2,3,4,5], print("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!") exit() for i in range(1,len(nss_temp)): - #trapizoid wall - #+(P[i]+P[i-1])/2*prior(center) - ns_val=(nss_temp[i]+nss_temp[i-1])/2 - deltaNs=nss_temp[i]-nss_temp[i-1] + #add prior*likelihood*dFlux at each value + ns_val=(nss_temp[i]+nss_temp[i-1])/2 # + deltaNs=nss_temp[i]-nss_temp[i-1] prior=prior_func(ns_val,self.index,pix[0],pix[1]) finalProbs[pix]=logsumexp([finalProbs[pix],logsumexp([probs_temp[i],probs_temp[i-1]])+np.log(prior)+np.log(deltaNs/2)]) #Normalize totalLog=logsumexp(list(finalProbs.values())) finalProbs={k:v-totalLog for k,v in finalProbs.items()} - print("N_pixels",len(finalProbs)) + print("N_pixels in posterior",len(finalProbs)) #Write to healpy outputProbHPMap=np.zeros(12*nside**2) ras_graph,decs_graph=zip(*list(finalProbs.keys())) @@ -1245,9 +1223,6 @@ def make_posterior_map(self, fluxToTest=[1,2,3,4,5], nPixZoom=sizeZoomInDegs/reso*60 - #hp.visufunc.gnomview(map=outputProbHPMap,hold=True,title="Zoomed Posterior Skymap",rot=(self.src_ra*180/np.pi,self.src_dec*180/np.pi,0), - # xsize=nPixZoom,ysize=nPixZoom,reso=reso,bgcolor="white",badcolor="white",margins=(.01,.01,0,0),notext=True,) - hp.visufunc.mollview(map=outputProbHPMapPlotting,hold=True,title="Posterior Skymap", bgcolor="white",badcolor="white",rot=(180,0,0),) hp.graticule(verbose=True) @@ -1363,10 +1338,13 @@ def make_posterior_map(self, fluxToTest=[1,2,3,4,5], obj.set_fontsize(30) plotting_utils.plot_color_bar(range=[0,np.nanmax(outputProbHPMapPlotting)], cmap="viridis", col_label=r"Posterior PDF", offset=-40,labels=[0,"{0:.1e}".format(np.nanmax(outputProbHPMapPlotting)/2),"{0:.2e}".format(np.nanmax(outputProbHPMapPlotting))],loc=[0.86, 0.2, 0.03, 0.6]) + theta, phi =plotting_utils.plot_contours([.9],outputProbHPMap) hp.projplot(theta[0], phi[0], linewidth=2., c='k') for i in range(1, len(theta)): hp.projplot(theta[i], phi[i], linewidth=2., c='k', label=None) + for i in range(0, len(theta)): + print(">>>",list(theta[i]),list(phi[i])) events = self.llh.exp events = events[(events['time'] < self.stop) & (events['time'] > self.start)] @@ -1396,12 +1374,12 @@ def make_posterior_map(self, fluxToTest=[1,2,3,4,5], obj.set_fontsize(30) plotting_utils.plot_color_bar(range=[0,np.nanmax(outputProbHPMapPlotting)], cmap="viridis", col_label=r"Posterior PDF", offset=-40,labels=[0,"{0:.1e}".format(np.nanmax(outputProbHPMapPlotting)/2),"{0:.2e}".format(np.nanmax(outputProbHPMapPlotting))],loc=[0.86, 0.2, 0.03, 0.6]) - print("90 contour",theta,phi) + theta, phi =plotting_utils.plot_contours([.9],outputProbHPMap) + hp.projplot(theta[0], phi[0], linewidth=2., c='k') for i in range(1, len(theta)): hp.projplot(theta[i], phi[i], linewidth=2., c='k', label=None) - theta, phi =plotting_utils.plot_contours([.9],outputProbHPMap) - + plt.savefig(plotting_location+"PosteriorMapZoomed1Deg.png") plt.close() @@ -1450,12 +1428,13 @@ def make_posterior_map(self, fluxToTest=[1,2,3,4,5], plt.close() print("saving fits") - self._make_posterior_fits(outputProbHPMap) + #return the version with zeros rather than nan t3 = time.time() print("finished posterior saving, took {} s".format(t3-t1)) - return outputProbHPMap + self.PosteriorSkymap=outputProbHPMap + self._make_posterior_fits() diff --git a/fast_response/latex/posterior_report_general.tex b/fast_response/latex/posterior_report_general.tex index f76e224..0cd33b4 100644 --- a/fast_response/latex/posterior_report_general.tex +++ b/fast_response/latex/posterior_report_general.tex @@ -46,7 +46,11 @@ \section{Posterior Maps} {\Large Posterior Skymap} \includegraphics[width=0.9\textwidth]{\PosteriorSkymap} + {\Large Zoomed Posterior Skymap} + \includegraphics[width=0.9\textwidth]{\PosteriorSkymapZoom} + {\Large Posterior Skymap W/ Nu Overlay} + \includegraphics[width=0.9\textwidth]{\PosteriorSkymapZoomWNeutrinos} diff --git a/fast_response/plotting_utils.py b/fast_response/plotting_utils.py index 99c0b57..fcd6c3b 100644 --- a/fast_response/plotting_utils.py +++ b/fast_response/plotting_utils.py @@ -107,11 +107,7 @@ def plot_labels(src_dec, src_ra, reso,nPix=200,label_scale=1): ha='center', va='center', rotation=90, fontsize=fontsize) plt.text(np.radians(0), np.radians(-2.05*reso*label_scale), r"right ascension", ha='center', va='center', fontsize=fontsize) - -def plot_events2(dec,ra,sigmas, src_ra, src_dec, reso, sigma_scale=5., col = 'k', constant_sigma=False, - same_marker=False, energy_size=False, with_mark=True, with_dash=False, - label='',resolution=0): - def compute_ang_err(ra,dec,sigma): +def compute_ang_err(ra,dec,sigma): dec = np.pi/2 - dec sigma = np.rad2deg(sigma) delta, step, bins = 0, 0, 0 @@ -132,7 +128,10 @@ def compute_ang_err(ra,dec,sigma): Phi[bins] = Phi[0] return Theta, Phi - +def plot_events2(dec,ra,sigmas, src_ra, src_dec, reso, sigma_scale=5., col = 'k', constant_sigma=False, + same_marker=False, energy_size=False, with_mark=True, with_dash=False, + label='',resolution=0): + #plot based on explicitly determined contours rather than using the s parameter. for i in range(len(ra)): ev_contour = compute_ang_err(ra[i],dec[i],sigmas[i]) hp.projplot(ev_contour[0], ev_contour[1], linewidth=1.75, color=col[0],