diff --git a/cscripts/csspec.par b/cscripts/csspec.par index f450bed..8a71218 100755 --- a/cscripts/csspec.par +++ b/cscripts/csspec.par @@ -42,6 +42,8 @@ calc_ts, b, h, yes,,, "Compute TS value in each bin" calc_ulim, b, h, yes,,, "Compute upper limit in each bin" fix_srcs, b, h, yes,,, "Fix other skymodel parameters" fix_bkg, b, h, no,,, "Fix background parameters" +dll_sigstep, r, h, 0,0,100, "Step size in standard deviations for log-like profiles" +dll_sigmax, r, h, 5,1,100, "Maximum number of standard deviations for log-like profiles" # # Standard parameters diff --git a/cscripts/csspec.py b/cscripts/csspec.py index af952a6..8b7fc20 100755 --- a/cscripts/csspec.py +++ b/cscripts/csspec.py @@ -19,6 +19,7 @@ # # ========================================================================== import sys +import math import gammalib import ctools from cscripts import mputils @@ -178,6 +179,10 @@ class csspec(ctools.csobservation): self['fix_bkg'].boolean() self['fix_srcs'].boolean() + # Setup dlog-likelihood parameters + self['dll_sigstep'].real() + self['dll_sigmax'].real() + # Read ahead output parameters if self._read_ahead(): self['outfile'].filename() @@ -651,7 +656,10 @@ class csspec(ctools.csobservation): 'flux_err': 0.0, 'TS': 0.0, 'ulimit': 0.0, - 'Npred': 0.0} + 'Npred': 0.0, + 'logL': 0.0, + 'dloglike': [], + 'norm_scan': []} # Write header for fitting self._log_header3(gammalib.EXPLICIT, 'Performing fit in energy bin') @@ -679,6 +687,16 @@ class csspec(ctools.csobservation): fitted_models = like.obs().models() source = fitted_models[self['srcname'].string()] + # Compute delta log-likelihood + if self['dll_sigstep'].real() > 0.0: + + # Extract the normalization scan values + parname = source.spectral()[0].name() + (norm_scan, dlogL, loglike) = self._profile_logL(like.copy(), parname, elogmean) + result['norm_scan'] = norm_scan + result['dloglike'] = dlogL + result['logL'] = loglike + # Extract Test Statistic value if self['calc_ts'].boolean(): result['TS'] = source.ts() @@ -850,7 +868,10 @@ class csspec(ctools.csobservation): 'flux_err': 0.0, 'TS': 0.0, 'ulimit': 0.0, - 'Npred': 0.0} + 'Npred': 0.0, + 'logL': 0.0, + 'dloglike': [], + 'norm_scan': []} # Convert differential flux and flux error to nuFnu norm = elogmean.MeV() * elogmean.MeV() * gammalib.MeV2erg @@ -859,6 +880,17 @@ class csspec(ctools.csobservation): # Compute upper flux limit ulimit_value = -1.0 + parname = 'Intensity%d' % i + + # Compute delta log-likelihood + if self['dll_sigstep'].real() > 0.0: + + # Extract the normalization scan values + (norm_scan, dlogL, loglike) = self._profile_logL(like.copy(), parname, elogmean) + result['norm_scan'] = norm_scan + result['dloglike'] = dlogL + result['logL'] = loglike + if self['calc_ulim'].boolean(): # Logging information @@ -877,7 +909,7 @@ class csspec(ctools.csobservation): # Create upper limit object ulimit = ctools.ctulimit(obs) ulimit['srcname'] = self['srcname'].string() - ulimit['parname'] = 'Intensity%d' % i + ulimit['parname'] = parname ulimit['eref'] = elogmean.TeV() ulimit['tol'] = 1.0e-3 @@ -935,6 +967,85 @@ class csspec(ctools.csobservation): # Return results return results + def _profile_logL(self, like, parname, elogmean): + """ + Computes the delta log-likelihood profile in a single energy bin + + Parameters + ---------- + like : `~ctools.ctlike()` + ctlike fitter containing prefit model + parname : str + Name of the spectral parameter to be fit + elogmean : `~gammalib.GEnergy()` + Energy at which the model is to be evaluated + + Returns + ------- + norm_scan : list + Normalization values + dloglike_scan : list + Computed loglikelihood values + loglike: float + Computed reference log-likelihood for dloglike_scan + """ + # Compute number of steps + dll_sigmax = self['dll_sigmax'].real() + dll_sigstep = self['dll_sigstep'].real() + sigsteps = int(2 * (dll_sigmax/dll_sigstep) + 1) + scale = elogmean.MeV() * elogmean.MeV() * gammalib.MeV2erg + + # Setup the source model for fitting + source = like.obs().models()[self['srcname'].string()] + spectral = source.spectral() + flux_par = spectral[parname] + norm = flux_par.value() + norm_err = flux_par.error() + source.tscalc(False) + + # Fix all parameters in the spectral model + for par in spectral: + par.fix() + + # Re-compute the log-likelihood + like.run() + loglike = like.obs().logL() + + # Store the resulting log-likelihood + norm_scan = [] + dloglike = [] + + # Compute the flux values to evaluate loglike at + log_norm = math.log10(norm) + if (norm_err > 0.0) and (norm > norm_err): + log_step = log_norm - math.log10(norm-norm_err) + log_start = log_norm - (sigsteps/2) * log_step + else: + # For an upper limit bin use a broad range of steps in flux + # from [10^-24, 10^-15] or [norm, 10^-15] whichever is broader + log_start = -24 if (log_norm > -24) else log_norm + log_step = (-15 - log_start) / (sigsteps-1) + + norm_vals = [10 ** (log_start + i*log_step) for i in range(sigsteps)] + + # Loop through normalizations + flux_par.factor_min(0.0) + flux_par.factor_max(1e30) + for new_norm in norm_vals: + flux_par.value(new_norm) + + # Re-run the fit + like.run() + + # Store dlikelihood & norm + dloglike.append(loglike - like.obs().logL()) + if self._method == 'SLICE': + norm_scan.append(spectral.eval(elogmean) * scale) + else: + norm_scan.append(new_norm * scale) + + return (norm_scan, dloglike, -loglike) + def _create_fits(self, results): """ Create FITS file @@ -946,6 +1057,7 @@ class csspec(ctools.csobservation): """ # Create FITS table columns nrows = self._ebounds.size() + ncols = len(results[0]['norm_scan']) energy = gammalib.GFitsTableDoubleCol('Energy', nrows) energy_low = gammalib.GFitsTableDoubleCol('ed_Energy', nrows) energy_high = gammalib.GFitsTableDoubleCol('eu_Energy', nrows) @@ -954,12 +1066,18 @@ class csspec(ctools.csobservation): TSvalues = gammalib.GFitsTableDoubleCol('TS', nrows) ulim_values = gammalib.GFitsTableDoubleCol('UpperLimit', nrows) Npred_values = gammalib.GFitsTableDoubleCol('Npred', nrows) + loglike = gammalib.GFitsTableDoubleCol('loglike', nrows) + norm_scan = gammalib.GFitsTableDoubleCol('norm_scan', nrows, ncols) + dloglike_scan= gammalib.GFitsTableDoubleCol('dloglike_scan', nrows, ncols) energy.unit('TeV') energy_low.unit('TeV') energy_high.unit('TeV') flux.unit('erg/cm2/s') flux_err.unit('erg/cm2/s') ulim_values.unit('erg/cm2/s') + loglike.unit('') + norm_scan.unit('erg/cm2/s') + dloglike_scan.unit('') # File FITS table columns for i, result in enumerate(results): @@ -971,6 +1089,12 @@ class csspec(ctools.csobservation): TSvalues[i] = result['TS'] ulim_values[i] = result['ulimit'] Npred_values[i] = result['Npred'] + loglike[i] = result['logL'] + + # Add likelihood scan values + for fbin in range(ncols): + dloglike_scan[i,fbin] = result['dloglike'][fbin] + norm_scan[i,fbin] = result['norm_scan'][fbin] # Initialise FITS Table with extension "SPECTRUM" table = gammalib.GFitsBinTable(nrows) @@ -990,6 +1114,20 @@ class csspec(ctools.csobservation): table.append(ulim_values) table.append(Npred_values) + # Define the SED type + sed_type = 'e2dnde' + if ncols > 0: + sed_type = 'likelihood' + table.append(loglike) + table.append(norm_scan) + table.append(dloglike_scan) + table.card('SED_TYPE', sed_type, 'SED type') + + # Define the upper limit confidence level + ulimit = ctools.ctulimit() + table.card('UL_CONF', ulimit['confidence'].real(), + 'Confidence level of upper limits') + # Create the FITS file now self._fits = gammalib.GFits() self._fits.append(table) diff --git a/examples/show_spectrum.py b/examples/show_spectrum.py index 0cd9549..9634693 100755 --- a/examples/show_spectrum.py +++ b/examples/show_spectrum.py @@ -19,6 +19,7 @@ # # ========================================================================== import sys +import math try: import matplotlib.pyplot as plt plt.figure() @@ -26,6 +27,11 @@ try: except (ImportError, RuntimeError): print('This script needs the "matplotlib" module') sys.exit() +try: + import numpy as np +except ImportError: + print('This script needs the "numpy" module') + sys.exit() import gammalib import cscripts @@ -92,9 +98,21 @@ def get_spectrum_fits(fits): 'ul_ed_engs' : [], 'ul_eu_engs' : [], 'ul_flux' : [], - 'yerr' : [] + 'yerr' : [], + 'loglike' : [], + 'engs_scan' : [], + 'norm_scan' : [], + 'dll_scan' : [] } + try: + c_loglike = table['loglike'] + c_norm_scan = table['norm_scan'] + c_dll_scan = table['dloglike_scan'] + load_dll = True + except: + load_dll = False + # Loop over rows of the file nrows = table.nrows() for row in range(nrows): @@ -120,6 +138,20 @@ def get_spectrum_fits(fits): spec['ul_ed_engs'].append(c_ed.real(row)) spec['ul_eu_engs'].append(c_eu.real(row)) + # Load the delta log-likelihood values + if load_dll: + spec['norm_scan'] = np.zeros((nrows, c_norm_scan.elements(0))) + spec['dll_scan'] = np.zeros((nrows, c_norm_scan.elements(0))) + + for row in range(nrows): + loglike = c_loglike.real(row) + spec['loglike'].append(loglike) + spec['engs_scan'].append(c_energy.real(row) - c_ed.real(row)) + for col in range(c_norm_scan.elements(row)): + spec['norm_scan'][row,col] = c_norm_scan.real(row,col) + spec['dll_scan'][row,col] = c_dll_scan.real(row,col) + spec['engs_scan'].append(c_energy.real(nrows-1)+c_eu.real(nrows-1)) + # Set upper limit errors spec['yerr'] = [0.6 * x for x in spec['ul_flux']] @@ -127,6 +159,48 @@ def get_spectrum_fits(fits): return spec +def plot_dloglike(spec): + """ + Plot delta log-likelihood + + Parameters + ---------- + spec : dict + Python dictionary defining spectral plot parameters + """ + # Establish the bounds on the x,y plot + ymin,ymax = plt.gca().get_ylim() + ymin = math.log10(ymin) + ymax = math.log10(ymax) + steps = 1000 + y_bounds = np.linspace(ymin, ymax, steps+1) + fluxpnts = [] + for i in range(steps): + mid_pnt = (y_bounds[i] + y_bounds[i+1]) / 2.0 + fluxpnts.append( 10 ** mid_pnt ) + # Scale the points + y_bounds = 10 ** y_bounds + + dll_hist = [] + for ebin in range(len(spec['norm_scan'])): + + # Interpolate the dlogLike values + dll_interp = np.interp(fluxpnts, spec['norm_scan'][ebin], + spec['dll_scan'][ebin]) + dll_hist.append(dll_interp) + + # Transpose the dll_hist array + dll_hist = [[row[i] for row in dll_hist] for i in range(len(dll_hist[0]))] + + plt.pcolormesh(spec['engs_scan'], y_bounds, dll_hist, vmin=-10.0, vmax=0.0, + cmap='Reds', alpha=0.9) + cbar = plt.colorbar() + cbar.set_label(r'$\Delta$ log-likelihood') + plt.grid() + + return + + # ============= # # Plot spectrum # # ============= # @@ -160,6 +234,9 @@ def plot_spectrum(filename, plotfile): yerr=spec['yerr'], xerr=[spec['ul_ed_engs'], spec['ul_eu_engs']], uplims=True, fmt='ro') + if len(spec['dll_scan']) > 0: + plot_dloglike(spec) + # Show figure if len(plotfile) > 0: plt.savefig(plotfile)