Pierrick Martin / ctools

Commit 1729ff4a0643638b1731f00e9ef67335c9b676d8

Authored by Josh Cardenzana
1 parent ad59645c

Add support for generating likelihood profiles in csspec 

Update 'show_spectrum.py' to also plot the spectral profiles
Inline Side-by-side
Showing 3 changed files with 221 additions and 4 deletions
cscripts/csspec.par
  View file @1729ff4
... ... @@ -42,6 +42,8 @@ calc_ts, b, h, yes,,, "Compute TS value in each bin"
42 42 calc_ulim, b, h, yes,,, "Compute upper limit in each bin"
43 43 fix_srcs, b, h, yes,,, "Fix other skymodel parameters"
44 44 fix_bkg, b, h, no,,, "Fix background parameters"
  45 +dll_sigstep, r, h, 0,0,100, "Step size in standard deviations for log-like profiles"
  46 +dll_sigmax, r, h, 5,1,100, "Maximum number of standard deviations for log-like profiles"
45 47  
46 48 #
47 49 # Standard parameters
... ...
cscripts/csspec.py
  View file @1729ff4
... ... @@ -19,6 +19,7 @@
19 19 #
20 20 # ==========================================================================
21 21 import sys
  22 +import math
22 23 import gammalib
23 24 import ctools
24 25 from cscripts import mputils
... ... @@ -178,6 +179,10 @@ class csspec(ctools.csobservation):
178 179 self['fix_bkg'].boolean()
179 180 self['fix_srcs'].boolean()
180 181  
  182 + # Setup dlog-likelihood parameters
  183 + self['dll_sigstep'].real()
  184 + self['dll_sigmax'].real()
  185 +
181 186 # Read ahead output parameters
182 187 if self._read_ahead():
183 188 self['outfile'].filename()
... ... @@ -651,7 +656,10 @@ class csspec(ctools.csobservation):
651 656 'flux_err': 0.0,
652 657 'TS': 0.0,
653 658 'ulimit': 0.0,
654   - 'Npred': 0.0}
  659 + 'Npred': 0.0,
  660 + 'logL': 0.0,
  661 + 'dloglike': [],
  662 + 'norm_scan': []}
655 663  
656 664 # Write header for fitting
657 665 self._log_header3(gammalib.EXPLICIT, 'Performing fit in energy bin')
... ... @@ -679,6 +687,16 @@ class csspec(ctools.csobservation):
679 687 fitted_models = like.obs().models()
680 688 source = fitted_models[self['srcname'].string()]
681 689  
  690 + # Compute delta log-likelihood
  691 + if self['dll_sigstep'].real() > 0.0:
  692 +
  693 + # Extract the normalization scan values
  694 + parname = source.spectral()[0].name()
  695 + (norm_scan, dlogL, loglike) = self._profile_logL(like.copy(), parname, elogmean)
  696 + result['norm_scan'] = norm_scan
  697 + result['dloglike'] = dlogL
  698 + result['logL'] = loglike
  699 +
682 700 # Extract Test Statistic value
683 701 if self['calc_ts'].boolean():
684 702 result['TS'] = source.ts()
... ... @@ -850,7 +868,10 @@ class csspec(ctools.csobservation):
850 868 'flux_err': 0.0,
851 869 'TS': 0.0,
852 870 'ulimit': 0.0,
853   - 'Npred': 0.0}
  871 + 'Npred': 0.0,
  872 + 'logL': 0.0,
  873 + 'dloglike': [],
  874 + 'norm_scan': []}
854 875  
855 876 # Convert differential flux and flux error to nuFnu
856 877 norm = elogmean.MeV() * elogmean.MeV() * gammalib.MeV2erg
... ... @@ -859,6 +880,17 @@ class csspec(ctools.csobservation):
859 880  
860 881 # Compute upper flux limit
861 882 ulimit_value = -1.0
  883 + parname = 'Intensity%d' % i
  884 +
  885 + # Compute delta log-likelihood
  886 + if self['dll_sigstep'].real() > 0.0:
  887 +
  888 + # Extract the normalization scan values
  889 + (norm_scan, dlogL, loglike) = self._profile_logL(like.copy(), parname, elogmean)
  890 + result['norm_scan'] = norm_scan
  891 + result['dloglike'] = dlogL
  892 + result['logL'] = loglike
  893 +
862 894 if self['calc_ulim'].boolean():
863 895  
864 896 # Logging information
... ... @@ -877,7 +909,7 @@ class csspec(ctools.csobservation):
877 909 # Create upper limit object
878 910 ulimit = ctools.ctulimit(obs)
879 911 ulimit['srcname'] = self['srcname'].string()
880   - ulimit['parname'] = 'Intensity%d' % i
  912 + ulimit['parname'] = parname
881 913 ulimit['eref'] = elogmean.TeV()
882 914 ulimit['tol'] = 1.0e-3
883 915  
... ... @@ -935,6 +967,85 @@ class csspec(ctools.csobservation):
935 967 # Return results
936 968 return results
937 969  
  970 + def _profile_logL(self, like, parname, elogmean):
  971 + """
  972 + Computes the delta log-likelihood profile in a single energy bin
  973 +
  974 + Parameters
  975 + ----------
  976 + like : `~ctools.ctlike()`
  977 + ctlike fitter containing prefit model
  978 + parname : str
  979 + Name of the spectral parameter to be fit
  980 + elogmean : `~gammalib.GEnergy()`
  981 + Energy at which the model is to be evaluated
  982 +
  983 + Returns
  984 + -------
  985 + norm_scan : list
  986 + Normalization values
  987 + dloglike_scan : list
  988 + Computed loglikelihood values
  989 + loglike: float
  990 + Computed reference log-likelihood for dloglike_scan
  991 + """
  992 + # Compute number of steps
  993 + dll_sigmax = self['dll_sigmax'].real()
  994 + dll_sigstep = self['dll_sigstep'].real()
  995 + sigsteps = int(2 * (dll_sigmax/dll_sigstep) + 1)
  996 + scale = elogmean.MeV() * elogmean.MeV() * gammalib.MeV2erg
  997 +
  998 + # Setup the source model for fitting
  999 + source = like.obs().models()[self['srcname'].string()]
  1000 + spectral = source.spectral()
  1001 + flux_par = spectral[parname]
  1002 + norm = flux_par.value()
  1003 + norm_err = flux_par.error()
  1004 + source.tscalc(False)
  1005 +
  1006 + # Fix all parameters in the spectral model
  1007 + for par in spectral:
  1008 + par.fix()
  1009 +
  1010 + # Re-compute the log-likelihood
  1011 + like.run()
  1012 + loglike = like.obs().logL()
  1013 +
  1014 + # Store the resulting log-likelihood
  1015 + norm_scan = []
  1016 + dloglike = []
  1017 +
  1018 + # Compute the flux values to evaluate loglike at
  1019 + log_norm = math.log10(norm)
  1020 + if (norm_err > 0.0) and (norm > norm_err):
  1021 + log_step = log_norm - math.log10(norm-norm_err)
  1022 + log_start = log_norm - (sigsteps/2) * log_step
  1023 + else:
  1024 + # For an upper limit bin use a broad range of steps in flux
  1025 + # from [10^-24, 10^-15] or [norm, 10^-15] whichever is broader
  1026 + log_start = -24 if (log_norm > -24) else log_norm
  1027 + log_step = (-15 - log_start) / (sigsteps-1)
  1028 +
  1029 + norm_vals = [10 ** (log_start + i*log_step) for i in range(sigsteps)]
  1030 +
  1031 + # Loop through normalizations
  1032 + flux_par.factor_min(0.0)
  1033 + flux_par.factor_max(1e30)
  1034 + for new_norm in norm_vals:
  1035 + flux_par.value(new_norm)
  1036 +
  1037 + # Re-run the fit
  1038 + like.run()
  1039 +
  1040 + # Store dlikelihood & norm
  1041 + dloglike.append(loglike - like.obs().logL())
  1042 + if self._method == 'SLICE':
  1043 + norm_scan.append(spectral.eval(elogmean) * scale)
  1044 + else:
  1045 + norm_scan.append(new_norm * scale)
  1046 +
  1047 + return (norm_scan, dloglike, -loglike)
  1048 +
938 1049 def _create_fits(self, results):
939 1050 """
940 1051 Create FITS file
... ... @@ -946,6 +1057,7 @@ class csspec(ctools.csobservation):
946 1057 """
947 1058 # Create FITS table columns
948 1059 nrows = self._ebounds.size()
  1060 + ncols = len(results[0]['norm_scan'])
949 1061 energy = gammalib.GFitsTableDoubleCol('Energy', nrows)
950 1062 energy_low = gammalib.GFitsTableDoubleCol('ed_Energy', nrows)
951 1063 energy_high = gammalib.GFitsTableDoubleCol('eu_Energy', nrows)
... ... @@ -954,12 +1066,18 @@ class csspec(ctools.csobservation):
954 1066 TSvalues = gammalib.GFitsTableDoubleCol('TS', nrows)
955 1067 ulim_values = gammalib.GFitsTableDoubleCol('UpperLimit', nrows)
956 1068 Npred_values = gammalib.GFitsTableDoubleCol('Npred', nrows)
  1069 + loglike = gammalib.GFitsTableDoubleCol('loglike', nrows)
  1070 + norm_scan = gammalib.GFitsTableDoubleCol('norm_scan', nrows, ncols)
  1071 + dloglike_scan= gammalib.GFitsTableDoubleCol('dloglike_scan', nrows, ncols)
957 1072 energy.unit('TeV')
958 1073 energy_low.unit('TeV')
959 1074 energy_high.unit('TeV')
960 1075 flux.unit('erg/cm2/s')
961 1076 flux_err.unit('erg/cm2/s')
962 1077 ulim_values.unit('erg/cm2/s')
  1078 + loglike.unit('')
  1079 + norm_scan.unit('erg/cm2/s')
  1080 + dloglike_scan.unit('')
963 1081  
964 1082 # File FITS table columns
965 1083 for i, result in enumerate(results):
... ... @@ -971,6 +1089,12 @@ class csspec(ctools.csobservation):
971 1089 TSvalues[i] = result['TS']
972 1090 ulim_values[i] = result['ulimit']
973 1091 Npred_values[i] = result['Npred']
  1092 + loglike[i] = result['logL']
  1093 +
  1094 + # Add likelihood scan values
  1095 + for fbin in range(ncols):
  1096 + dloglike_scan[i,fbin] = result['dloglike'][fbin]
  1097 + norm_scan[i,fbin] = result['norm_scan'][fbin]
974 1098  
975 1099 # Initialise FITS Table with extension "SPECTRUM"
976 1100 table = gammalib.GFitsBinTable(nrows)
... ... @@ -990,6 +1114,20 @@ class csspec(ctools.csobservation):
990 1114 table.append(ulim_values)
991 1115 table.append(Npred_values)
992 1116  
  1117 + # Define the SED type
  1118 + sed_type = 'e2dnde'
  1119 + if ncols > 0:
  1120 + sed_type = 'likelihood'
  1121 + table.append(loglike)
  1122 + table.append(norm_scan)
  1123 + table.append(dloglike_scan)
  1124 + table.card('SED_TYPE', sed_type, 'SED type')
  1125 +
  1126 + # Define the upper limit confidence level
  1127 + ulimit = ctools.ctulimit()
  1128 + table.card('UL_CONF', ulimit['confidence'].real(),
  1129 + 'Confidence level of upper limits')
  1130 +
993 1131 # Create the FITS file now
994 1132 self._fits = gammalib.GFits()
995 1133 self._fits.append(table)
... ...
examples/show_spectrum.py
  View file @1729ff4
... ... @@ -19,6 +19,7 @@
19 19 #
20 20 # ==========================================================================
21 21 import sys
  22 +import math
22 23 try:
23 24 import matplotlib.pyplot as plt
24 25 plt.figure()
... ... @@ -26,6 +27,11 @@ try:
26 27 except (ImportError, RuntimeError):
27 28 print('This script needs the "matplotlib" module')
28 29 sys.exit()
  30 +try:
  31 + import numpy as np
  32 +except ImportError:
  33 + print('This script needs the "numpy" module')
  34 + sys.exit()
29 35 import gammalib
30 36 import cscripts
31 37  
... ... @@ -92,9 +98,21 @@ def get_spectrum_fits(fits):
92 98 'ul_ed_engs' : [],
93 99 'ul_eu_engs' : [],
94 100 'ul_flux' : [],
95   - 'yerr' : []
  101 + 'yerr' : [],
  102 + 'loglike' : [],
  103 + 'engs_scan' : [],
  104 + 'norm_scan' : [],
  105 + 'dll_scan' : []
96 106 }
97 107  
  108 + try:
  109 + c_loglike = table['loglike']
  110 + c_norm_scan = table['norm_scan']
  111 + c_dll_scan = table['dloglike_scan']
  112 + load_dll = True
  113 + except:
  114 + load_dll = False
  115 +
98 116 # Loop over rows of the file
99 117 nrows = table.nrows()
100 118 for row in range(nrows):
... ... @@ -120,6 +138,20 @@ def get_spectrum_fits(fits):
120 138 spec['ul_ed_engs'].append(c_ed.real(row))
121 139 spec['ul_eu_engs'].append(c_eu.real(row))
122 140  
  141 + # Load the delta log-likelihood values
  142 + if load_dll:
  143 + spec['norm_scan'] = np.zeros((nrows, c_norm_scan.elements(0)))
  144 + spec['dll_scan'] = np.zeros((nrows, c_norm_scan.elements(0)))
  145 +
  146 + for row in range(nrows):
  147 + loglike = c_loglike.real(row)
  148 + spec['loglike'].append(loglike)
  149 + spec['engs_scan'].append(c_energy.real(row) - c_ed.real(row))
  150 + for col in range(c_norm_scan.elements(row)):
  151 + spec['norm_scan'][row,col] = c_norm_scan.real(row,col)
  152 + spec['dll_scan'][row,col] = c_dll_scan.real(row,col)
  153 + spec['engs_scan'].append(c_energy.real(nrows-1)+c_eu.real(nrows-1))
  154 +
123 155 # Set upper limit errors
124 156 spec['yerr'] = [0.6 * x for x in spec['ul_flux']]
125 157  
... ... @@ -127,6 +159,48 @@ def get_spectrum_fits(fits):
127 159 return spec
128 160  
129 161  
  162 +def plot_dloglike(spec):
  163 + """
  164 + Plot delta log-likelihood
  165 +
  166 + Parameters
  167 + ----------
  168 + spec : dict
  169 + Python dictionary defining spectral plot parameters
  170 + """
  171 + # Establish the bounds on the x,y plot
  172 + ymin,ymax = plt.gca().get_ylim()
  173 + ymin = math.log10(ymin)
  174 + ymax = math.log10(ymax)
  175 + steps = 1000
  176 + y_bounds = np.linspace(ymin, ymax, steps+1)
  177 + fluxpnts = []
  178 + for i in range(steps):
  179 + mid_pnt = (y_bounds[i] + y_bounds[i+1]) / 2.0
  180 + fluxpnts.append( 10 ** mid_pnt )
  181 + # Scale the points
  182 + y_bounds = 10 ** y_bounds
  183 +
  184 + dll_hist = []
  185 + for ebin in range(len(spec['norm_scan'])):
  186 +
  187 + # Interpolate the dlogLike values
  188 + dll_interp = np.interp(fluxpnts, spec['norm_scan'][ebin],
  189 + spec['dll_scan'][ebin])
  190 + dll_hist.append(dll_interp)
  191 +
  192 + # Transpose the dll_hist array
  193 + dll_hist = [[row[i] for row in dll_hist] for i in range(len(dll_hist[0]))]
  194 +
  195 + plt.pcolormesh(spec['engs_scan'], y_bounds, dll_hist, vmin=-10.0, vmax=0.0,
  196 + cmap='Reds', alpha=0.9)
  197 + cbar = plt.colorbar()
  198 + cbar.set_label(r'$\Delta$ log-likelihood')
  199 + plt.grid()
  200 +
  201 + return
  202 +
  203 +
130 204 # ============= #
131 205 # Plot spectrum #
132 206 # ============= #
... ... @@ -160,6 +234,9 @@ def plot_spectrum(filename, plotfile):
160 234 yerr=spec['yerr'], xerr=[spec['ul_ed_engs'], spec['ul_eu_engs']],
161 235 uplims=True, fmt='ro')
162 236  
  237 + if len(spec['dll_scan']) > 0:
  238 + plot_dloglike(spec)
  239 +
163 240 # Show figure
164 241 if len(plotfile) > 0:
165 242 plt.savefig(plotfile)
... ...