diff --git a/docs/spectrum/flux_points.rst b/docs/spectrum/flux_points.rst deleted file mode 100644 index 54535d1e12..0000000000 --- a/docs/spectrum/flux_points.rst +++ /dev/null @@ -1,71 +0,0 @@ -.. include:: ../references.txt - -.. _flux-point-computation: - -********************** -Flux point computation -********************** - -.. currentmodule:: gammapy.spectrum - -In the following you will see how to compute -`~gammapy.spectrum.DifferentialFluxPoints` given a global model and a -`~gammapy.spectrum.SpectrumObservation`. We will use the example dataset in -`gammapy-extra -`_. -The flux points binning is chosen as 5 equally log-spaced bins between the -observation thresholds. In order to obtain the global model we first perform -the global fit again, for more info see :ref:`spectral_fitting`. - -.. code-block:: python - - import astropy.units as u - from gammapy.spectrum import SpectrumObservation, SpectrumFit, DifferentialFluxPoints - from gammapy.spectrum.models import PowerLaw - from gammapy.utils.energy import EnergyBounds - - pha = "$GAMMAPY_EXTRA/datasets/hess-crab4_pha/pha_obs23523.fits" - obs = SpectrumObservation.read(pha) - - model = PowerLaw(index = 2 * u.Unit(''), - amplitude = 10 ** -12 * u.Unit('cm-2 s-1 TeV-1'), - reference = 1 * u.TeV) - - global_fit = SpectrumFit(obs_list=[obs], model=model) - global_fit.fit() - - global_model = global_fit.result[0].fit.model - binning = EnergyBounds.equal_log_spacing(obs.lo_threshold, obs.hi_threshold, 5) - - points = DifferentialFluxPoints.compute(model=global_model, - binning=binning, - obs_list = [obs]) - - - -Note, that in this case (where we just performed the global fit) we can get the -flux points more conveniently as - -.. code-block:: python - - import astropy.units as u - from gammapy.spectrum import SpectrumObservation, SpectrumFit - from gammapy.spectrum.models import PowerLaw - from gammapy.utils.energy import EnergyBounds - import matplotlib.pyplot as plt - - pha = "$GAMMAPY_EXTRA/datasets/hess-crab4_pha/pha_obs23523.fits" - obs = SpectrumObservation.read(pha) - - model = PowerLaw(index = 2 * u.Unit(''), - amplitude = 10 ** -12 * u.Unit('cm-2 s-1 TeV-1'), - reference = 1 * u.TeV) - - fit = SpectrumFit(obs_list=[obs], model=model) - fit.fit() - - binning = EnergyBounds.equal_log_spacing(obs.lo_threshold, obs.hi_threshold, 5) - fit.compute_fluxpoints(binning=binning) - fit.result[0].plot_spectrum() - plt.show() - diff --git a/docs/spectrum/index.rst b/docs/spectrum/index.rst index 4188052aa7..0d1e6297c0 100644 --- a/docs/spectrum/index.rst +++ b/docs/spectrum/index.rst @@ -78,7 +78,6 @@ Content :maxdepth: 1 fitting - flux_points energy_group plotting_fermi_spectra diff --git a/gammapy/catalog/fermi.py b/gammapy/catalog/fermi.py index b38866c57b..c2584fcec5 100644 --- a/gammapy/catalog/fermi.py +++ b/gammapy/catalog/fermi.py @@ -17,7 +17,7 @@ from ..utils.table import table_standardise_units_inplace from ..image import SkyImage from ..image.models import Delta2D, Template2D -from ..spectrum import FluxPoints, compute_flux_points_dnde +from ..spectrum import FluxPoints from ..spectrum.models import ( PowerLaw, PowerLaw2, @@ -453,9 +453,9 @@ def flux_points(self): flux_points = FluxPoints(table) - flux_points_dnde = compute_flux_points_dnde( - flux_points, model=self.spectral_model) - return flux_points_dnde + # TODO: change this and leave it up to the caller to convert to dnde + # See https://github.com/gammapy/gammapy/issues/1034 + return flux_points.to_sed_type('dnde', model=self.spectral_model) @property def spectral_model(self): @@ -526,9 +526,9 @@ def flux_points(self): flux_points = FluxPoints(table) - flux_points_dnde = compute_flux_points_dnde( - flux_points, model=self.spectral_model) - return flux_points_dnde + # TODO: change this and leave it up to the caller to convert to dnde + # See https://github.com/gammapy/gammapy/issues/1034 + return flux_points.to_sed_type('dnde', model=self.spectral_model) @property def spectral_model(self): diff --git a/gammapy/scripts/tests/test_spectrum_pipe.py b/gammapy/scripts/tests/test_spectrum_pipe.py index bdc1aefcc2..3adbe8ec51 100644 --- a/gammapy/scripts/tests/test_spectrum_pipe.py +++ b/gammapy/scripts/tests/test_spectrum_pipe.py @@ -30,16 +30,16 @@ def get_config(): @requires_dependency('scipy') @requires_dependency('sherpa') def test_spectrum_analysis_iact(tmpdir): - conf = get_config() - conf['outdir'] = tmpdir + config = get_config() + config['outdir'] = tmpdir - ana = SpectrumAnalysisIACT( - observations=obs_list(), - config=conf - ) + analysis = SpectrumAnalysisIACT(observations=obs_list(), config=config) + analysis.run() + flux_points = analysis.flux_point_estimator.flux_points + + print(flux_points) + print(flux_points.table) - assert 'IACT' in str(ana) - ana.run() - actual = ana.flux_point_estimator.flux_points.table[0]['dnde'] - desired = 7.208219787928114e-08 + actual = flux_points.table['dnde'].quantity[0] + desired = 7.208219780210792e-12 * u.Unit('cm-2 s-1 TeV-1') assert_quantity_allclose(actual, desired) diff --git a/gammapy/spectrum/__init__.py b/gammapy/spectrum/__init__.py index 9aee46974f..baf2cd10ed 100644 --- a/gammapy/spectrum/__init__.py +++ b/gammapy/spectrum/__init__.py @@ -9,7 +9,6 @@ from .diffuse import * from .energy_group import * from .flux_point import * -from .sherpa_chi2asym import * from .utils import * from .extract import * from .simulation import * diff --git a/gammapy/spectrum/crab.py b/gammapy/spectrum/crab.py index 92fd864c3e..0e24d5ac0f 100644 --- a/gammapy/spectrum/crab.py +++ b/gammapy/spectrum/crab.py @@ -1,17 +1,4 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst -"""Published Crab nebula reference spectra. - -The Crab is often used as a standard candle in gamma-ray astronomy. -Statements like "this source has a flux of 10 % Crab" or -"our sensitivity is 2 % Crab" are common. - -Here we provide a reference of what the Crab flux actually is according -to several publications: - -* 'HEGRA' : http://adsabs.harvard.edu/abs/2000ApJ...539..317A -* 'HESS PL' and 'HESS ECPL' : http://adsabs.harvard.edu/abs/2006A%26A...457..899A -* 'Meyer' : http://adsabs.harvard.edu/abs/2010A%26A...523A...2M, Appendix D -""" from __future__ import absolute_import, division, print_function, unicode_literals import numpy as np from astropy import units as u @@ -23,39 +10,27 @@ ] # HESS publication: 2006A&A...457..899A -#'int_flux' = 2.26e-11 hess_pl = {'amplitude': 3.45e-11 * u.Unit('1 / (cm2 s TeV)'), 'index': 2.63, - 'reference' : 1 * u.TeV} + 'reference': 1 * u.TeV} -# Note that for ecpl, the diff_flux is not -# the differential flux at 1 TeV, that you -# get by multiplying with exp(-e / cutoff) -# int_flux: 2.27e-11 hess_ecpl = {'amplitude': 3.76e-11 * u.Unit('1 / (cm2 s TeV)'), 'index': 2.39, 'lambda_': 1 / (14.3 * u.TeV), 'reference': 1 * u.TeV} # HEGRA publication : 2004ApJ...614..897A -# int_flux': 1.75e-11 hegra = {'amplitude': 2.83e-11 * u.Unit('1 / (cm2 s TeV)'), 'index': 2.62, 'reference': 1 * u.TeV} -# Meyer et al. publication: 2010arXiv1008.4524M -# diff_flux and index were calculated numerically -# by hand at 1 TeV as a finite differene -meyer = {'diff_flux': 3.3457e-11 * u.Unit('1 / (cm2 s TeV)'), - 'index': 2.5362, - 'int_flux': 2.0744e-11} - -#TODO: make this a general LogPolynomial spectral model class MeyerCrabModel(SpectralModel): + """Meyer 2010 log polynomial Crab spectral model. + + See 2010A%26A...523A...2M, Appendix D. """ - Log polynomial model as used by 2010arXiv1008.4524M. - """ + def __init__(self): coefficients = np.array([-0.00449161, 0, 0.0473174, -0.179475, -0.53616, -10.2708]) @@ -73,15 +48,16 @@ def evaluate(energy, coefficients): class CrabSpectrum(object): - """ - Crab spectral model. + """Crab nebula spectral model. + + The Crab nebula is often used as a standard candle in gamma-ray astronomy. + Fluxes and sensitivities are often quoted relative to the Crab spectrum. The following references are available: - * 'meyer', 2010arXiv1008.4524M - * 'hegra', 2004ApJ...614..897A - * 'hess_pl', 2006A&A...457..899A - * 'hess_ecpl', 2006A&A...457..899A + * 'meyer', http://adsabs.harvard.edu/abs/2010A%26A...523A...2M, Appendix D + * 'hegra', http://adsabs.harvard.edu/abs/2000ApJ...539..317A + * 'hess_pl' and 'hess_ecpl': http://adsabs.harvard.edu/abs/2006A%26A...457..899A Parameters ---------- @@ -90,44 +66,42 @@ class CrabSpectrum(object): Examples -------- - Plot the 'hess_ecpl' reference crab spectrum between 1 TeV and 100 TeV: - - >>> from gammapy.spectrum import CrabSpectrum - >>> import astropy.units as u - >>> crab_hess_ecpl = CrabSpectrum('hess_ecpl') - >>> crab_hess_ecpl.model.plot([1, 100] * u.TeV) - - Use a reference crab spectrum as unit to measure flux: - - >>> from astropy import units as u - >>> from gammapy.spectrum import CrabSpectrum - >>> from gammapy.spectrum.models import PowerLaw - >>> - >>> # define power law model - >>> amplitude = 1e-12 * u.Unit('1 / (cm2 s TeV)') - >>> index = 2.3 - >>> reference = 1 * u.TeV - >>> pwl = PowerLaw(index, amplitude, reference) - >>> - >>> # define crab reference + Let's first import what we need:: + + import astropy.units as u + from gammapy.spectrum import CrabSpectrum + from gammapy.spectrum.models import PowerLaw + + Plot the 'hess_ecpl' reference Crab spectrum between 1 TeV and 100 TeV:: + + crab_hess_ecpl = CrabSpectrum('hess_ecpl') + crab_hess_ecpl.model.plot([1, 100] * u.TeV) + + Use a reference crab spectrum as unit to measure a differential flux (at 10 TeV):: + + >>> pwl = PowerLaw(index=2.3, amplitude=1e-12 * u.Unit('1 / (cm2 s TeV)'), reference=1 * u.TeV) >>> crab = CrabSpectrum('hess_pl').model - >>> - >>> # compute flux at 10 TeV in crab units >>> energy = 10 * u.TeV - >>> flux_cu = (pwl(energy) / crab(energy)).to('%') - >>> print(flux_cu) + >>> dnde_cu = (pwl(energy) / crab(energy)).to('%') + >>> print(dnde_cu) 6.19699156377 % - And the same for integral fluxes: + And the same for integral fluxes (between 1 and 10 TeV):: >>> # compute integral flux in crab units >>> emin, emax = [1, 10] * u.TeV >>> flux_int_cu = (pwl.integral(emin, emax) / crab.integral(emin, emax)).to('%') >>> print(flux_int_cu) 3.5350582166 % - """ + + references = [ + 'meyer', 'hegra', 'hess_pl', 'hess_ecpl', + ] + """Available references (see class docstring).""" + def __init__(self, reference='meyer'): + if reference == 'meyer': model = MeyerCrabModel() elif reference == 'hegra': @@ -137,7 +111,8 @@ def __init__(self, reference='meyer'): elif reference == 'hess_ecpl': model = ExponentialCutoffPowerLaw(**hess_ecpl) else: - raise ValueError("Unknown reference, choose one of the following:" - "'meyer', 'hegra', 'hess_pl' or 'hess_ecpl'") + fmt = 'Invalid reference: {!r}. Choices: {!r}' + raise ValueError(fmt.format(reference, self.references)) + self.model = model self.reference = reference diff --git a/gammapy/spectrum/extract.py b/gammapy/spectrum/extract.py index e1ce9f4f70..7351a427d9 100755 --- a/gammapy/spectrum/extract.py +++ b/gammapy/spectrum/extract.py @@ -18,7 +18,7 @@ class SpectrumExtraction(object): - """Creating input data to 1D spectrum fitting + """Creating input data to 1D spectrum fitting. This class is responsible for extracting a `~gammapy.spectrum.SpectrumObservation` from a diff --git a/gammapy/spectrum/flux_point.py b/gammapy/spectrum/flux_point.py index 85d35aac76..cf31a07a70 100644 --- a/gammapy/spectrum/flux_point.py +++ b/gammapy/spectrum/flux_point.py @@ -1,5 +1,4 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst -"""Differential and integral flux point computations.""" from __future__ import absolute_import, division, print_function, unicode_literals import logging from collections import OrderedDict @@ -10,13 +9,13 @@ from ..utils.scripts import make_path from ..utils.fits import table_from_row_data from ..utils.table import table_standardise_units_copy +from .models import PowerLaw __all__ = [ - 'compute_flux_points_dnde', 'FluxPoints', + # 'FluxPointProfiles', 'FluxPointEstimator', - 'FluxPointsFitter', - 'SEDLikelihoodProfile', + 'FluxPointFitter', ] log = logging.getLogger(__name__) @@ -67,10 +66,111 @@ class FluxPoints(object): """ def __init__(self, table): - table = table_standardise_units_copy(table) + self.table = table_standardise_units_copy(table) # validate that the table is a valid representation # of the given flux point sed type - self.table = self._validate_table(table) + self._validate_table(self.table) + + def __repr__(self): + fmt = '{}(sed_type="{}", n_points={})' + return fmt.format(self.__class__.__name__, self.sed_type, len(self.table)) + + @classmethod + def read(cls, filename, **kwargs): + """Read flux points. + + Parameters + ---------- + filename : str + Filename + kwargs : dict + Keyword arguments passed to `~astropy.table.Table.read`. + """ + filename = make_path(filename) + try: + table = Table.read(str(filename), **kwargs) + except IORegistryError: + kwargs.setdefault('format', 'ascii.ecsv') + table = Table.read(str(filename), **kwargs) + + if 'SED_TYPE' not in table.meta.keys(): + sed_type = cls._guess_sed_type(table) + table.meta['SED_TYPE'] = sed_type + + return cls(table=table) + + def write(self, filename, **kwargs): + """Write flux points. + + Parameters + ---------- + filename : str + Filename + kwargs : dict + Keyword arguments passed to `~astropy.table.Table.write`. + """ + filename = make_path(filename) + try: + self.table.write(str(filename), **kwargs) + except IORegistryError: + kwargs.setdefault('format', 'ascii.ecsv') + self.table.write(str(filename), **kwargs) + + @classmethod + def stack(cls, flux_points): + """Create flux points by stacking list of flux points. + + The first `FluxPoints` object in the list is taken as a reference to infer + column names and units for the stacked object. + + Parameters + ---------- + flux_points : list of `FluxPoints` + List of flux points to stack. + + Returns + ------- + flux_points : `FluxPoints` + Flux points without upper limit points. + """ + reference = flux_points[0].table + + tables = [] + for _ in flux_points: + table = _.table + for colname in reference.colnames: + column = reference[colname] + if column.unit: + table[colname] = table[colname].quantity.to(column.unit) + tables.append(table[reference.colnames]) + + table_stacked = vstack(tables) + table_stacked.meta['SED_TYPE'] = reference.meta['SED_TYPE'] + + return cls(table_stacked) + + def drop_ul(self): + """Drop upper limit flux points. + + Returns + ------- + flux_points : `FluxPoints` + Flux points with upper limit points removed. + + Examples + -------- + + :: + + from gammapy.spectrum import FluxPoints + filename = '$GAMMAPY_EXTRA/test_datasets/spectrum/flux_points/flux_points.fits' + flux_points = FluxPoints.read(filename) + print(flux_points) + print(flux_points.drop_ul()) + """ + table = self.table.copy() + table_drop_ul = table[~self._is_ul] + return self.__class__(table_drop_ul) @property def sed_type(self): @@ -96,9 +196,12 @@ def _guess_sed_type_from_unit(unit): if unit.is_equivalent(default_unit): return sed_type - def _validate_table(self, table): + @staticmethod + def _validate_table(table): """Validate input table.""" - table = Table(table) + # TODO: do we really want to error out on tables that don't have `SED_TYPE` in meta? + # If yes, then this needs to be documented in the docstring, + # and the workaround pointed out (to add the meta key before creating FluxPoints). sed_type = table.meta['SED_TYPE'] required = set(REQUIRED_COLUMNS[sed_type]) @@ -106,7 +209,6 @@ def _validate_table(self, table): missing = required.difference(table.colnames) raise ValueError("Missing columns for sed type '{}':" " {}".format(sed_type, missing)) - return table def _get_y_energy_unit(self, y_unit): """Get energy part of the given y unit.""" @@ -155,62 +257,6 @@ def _is_ul(self): except KeyError: return np.isnan(self.table[self.sed_type]) - def __str__(self): - """ - String representation of the flux points class. - """ - info = '' - info += "Flux points of type '{}'".format(self.sed_type) - return info - - def info(self): - """ - Print flux points info. - """ - print(self) - - @classmethod - def read(cls, filename, **kwargs): - """Read flux points. - - Parameters - ---------- - filename : str - Filename - kwargs : dict - Keyword arguments passed to `~astropy.table.Table.read`. - - """ - filename = make_path(filename) - try: - table = Table.read(str(filename), **kwargs) - except IORegistryError: - kwargs.setdefault('format', 'ascii.ecsv') - table = Table.read(str(filename), **kwargs) - - if 'SED_TYPE' not in table.meta.keys(): - sed_type = cls._guess_sed_type(table) - table.meta['SED_TYPE'] = sed_type - - return cls(table=table) - - def write(self, filename, **kwargs): - """Write flux points. - - Parameters - ---------- - filename : str - Filename - kwargs : dict - Keyword arguments passed to `~astropy.table.Table.write`. - """ - filename = make_path(filename) - try: - self.table.write(str(filename), **kwargs) - except IORegistryError: - kwargs.setdefault('format', 'ascii.ecsv') - self.table.write(str(filename), **kwargs) - @property def e_ref(self): """Reference energy. @@ -254,60 +300,114 @@ def e_max(self): """ return self.table['e_max'].quantity - def drop_ul(self): - """Drop upper limit flux points. - + def to_sed_type(self, sed_type, method='log_center', model=None): + """Convert to a different SED type (return a new `FluxPoints` object). + + See: http://adsabs.harvard.edu/abs/1995NIMPA.355..541L for details + on the `'lafferty'` method. + + Parameters + ---------- + sed_type : {'dnde'} + SED type to convert to. + model : `~gammapy.spectrum.SpectralModel` + Spectral model assumption. Note that the value of the amplitude parameter + does not matter. Still it is recommended to use something with the right + scale and units. E.g. `amplitude = 1e-12 * u.Unit('cm-2 s-1 TeV-1')` + method : {'lafferty', 'log_center', 'table'} + Flux points `e_ref` estimation method: + + * `'laferty'` Lafferty & Wyatt model-based e_ref + * `'log_center'` log bin center e_ref + * `'table'` using column 'e_ref' from input flux_points + Returns ------- flux_points : `FluxPoints` - Flux points with upper limit points removed. - + Flux points including differential quantity columns `dnde` + and `dnde_err` (optional), `dnde_ul` (optional). + Examples -------- + >>> from astropy import units as u + >>> from gammapy.spectrum import FluxPoints + >>> from gammapy.spectrum.models import PowerLaw + >>> filename = '$GAMMAPY_EXTRA/test_datasets/spectrum/flux_points/flux_points.fits' + >>> flux_points = FluxPoints.read(filename) + >>> model = PowerLaw(2.2 * u.Unit(''), 1e-12 * u.Unit('cm-2 s-1 TeV-1'), 1 * u.TeV) + >>> flux_point_dnde = flux_points.to_sed_type('dnde', model=model) + """ + # TODO: implement other directions. Refactor! + if sed_type != 'dnde': + raise NotImplementedError + + if model is None: + model = PowerLaw( + index=2 * u.Unit(''), + amplitude=1 * u.Unit('cm-2 s-1 TeV-1'), + reference=1 * u.TeV, + ) + + input_table = self.table.copy() + + e_min, e_max = self.e_min, self.e_max + + # Compute e_ref + if method == 'table': + e_ref = input_table['e_ref'].quantity + elif method == 'log_center': + e_ref = np.sqrt(e_min * e_max) + elif method == 'lafferty': + # set e_ref that it represents the mean dnde in the given energy bin + e_ref = self._e_ref_lafferty(model, e_min, e_max) + else: + raise ValueError('Invalid method: {}'.format(method)) - :: + flux = input_table['flux'].quantity + dnde = self._dnde_from_flux(flux, model, e_ref, e_min, e_max) - from gammapy.spectrum import FluxPoints - filename = '$GAMMAPY_EXTRA/test_datasets/spectrum/flux_points/flux_points.fits' - flux_points = FluxPoints.read(filename) - print(flux_points) - print(flux_points.drop_ul()) - """ - table = self.table.copy() - table_drop_ul = table[~self._is_ul] - return self.__class__(table_drop_ul) + # Add to result table + table = input_table.copy() + table['e_ref'] = e_ref + table['dnde'] = dnde - @classmethod - def stack(cls, flux_points): - """Create flux points by stacking list of flux points. + if 'flux_err' in table.colnames: + # TODO: implement better error handling, e.g. MC based method + table['dnde_err'] = dnde * table['flux_err'].quantity / flux - The first `FluxPoints` object in the list is taken as a reference to infer - column names and units for the stacked object. + if 'flux_errn' in table.colnames: + table['dnde_errn'] = dnde * table['flux_errn'].quantity / flux + table['dnde_errp'] = dnde * table['flux_errp'].quantity / flux - Parameters - ---------- - flux_points : list of `FluxPoints` - List of flux points to stack. + if 'flux_ul' in table.colnames: + flux_ul = table['flux_ul'].quantity + dnde_ul = self._dnde_from_flux(flux_ul, model, e_ref, e_min, e_max) + table['dnde_ul'] = dnde_ul - Returns - ------- - flux_points : `FluxPoints` - Flux points without upper limit points. + table.meta['SED_TYPE'] = 'dnde' + return FluxPoints(table) + + @staticmethod + def _e_ref_lafferty(model, e_min, e_max): + """Helper for `to_sed_type`. + + Compute e_ref that the value at e_ref corresponds + to the mean value between e_min and e_max. """ - tables = [] - reference = flux_points[0].table + flux = model.integral(e_min, e_max) + dnde_mean = flux / (e_max - e_min) + return model.inverse(dnde_mean) - for _ in flux_points: - table = _.table - for colname in reference.colnames: - column = reference[colname] - if column.unit: - table[colname] = table[colname].quantity.to(column.unit) - tables.append(table[reference.colnames]) - table_stacked = vstack(tables) - table_stacked.meta['SED_TYPE'] = reference.meta['SED_TYPE'] + @staticmethod + def _dnde_from_flux(flux, model, e_ref, e_min, e_max): + """Helper for `to_sed_type`. - return cls(table_stacked) + Compute dnde under the assumption that flux equals expected + flux from model. + """ + flux_model = model.integral(e_min, e_max) + dnde_model = model(e_ref) + return dnde_model * (flux / flux_model) def plot(self, ax=None, sed_type=None, energy_unit='TeV', flux_unit=None, energy_power=0, **kwargs): @@ -396,103 +496,12 @@ def plot(self, ax=None, sed_type=None, energy_unit='TeV', flux_unit=None, return ax -def compute_flux_points_dnde(flux_points, model, method='lafferty'): - """ - Compute differential flux points quantities. - - See: http://adsabs.harvard.edu/abs/1995NIMPA.355..541L for details - on the `'lafferty'` method. - - Parameters - ---------- - flux_points : `FluxPoints` - Input integral flux points. - model : `~gammapy.spectrum.SpectralModel` - Spectral model assumption. Note that the value of the amplitude parameter - does not matter. Still it is recommended to use something with the right - scale and units. E.g. `amplitude = 1e-12 * u.Unit('cm-2 s-1 TeV-1')` - method : {'lafferty', 'log_center', 'table'} - Flux points `e_ref` estimation method: - - * `'laferty'` Lafferty & Wyatt model-based e_ref - * `'log_center'` log bin center e_ref - * `'table'` using column 'e_ref' from input flux_points - - Returns - ------- - flux_points : `FluxPoints` - Flux points including differential quantity columns `dnde` - and `dnde_err` (optional), `dnde_ul` (optional). - - Examples - -------- - >>> from astropy import units as u - >>> from gammapy.spectrum import FluxPoints, compute_flux_points_dnde - >>> from gammapy.spectrum.models import PowerLaw - >>> filename = '$GAMMAPY_EXTRA/test_datasets/spectrum/flux_points/flux_points.fits' - >>> flux_points = FluxPoints.read(filename) - >>> model = PowerLaw(2.2 * u.Unit(''), 1e-12 * u.Unit('cm-2 s-1 TeV-1'), 1 * u.TeV) - >>> flux_point_dnde = compute_flux_points_dnde(flux_points, model=model) - """ - input_table = flux_points.table - flux = input_table['flux'].quantity - e_min = flux_points.e_min - e_max = flux_points.e_max - - # Compute e_ref - if method == 'table': - e_ref = input_table['e_ref'].quantity - elif method == 'log_center': - e_ref = np.sqrt(e_min * e_max) - elif method == 'lafferty': - # set e_ref that it represents the mean dnde in the given energy bin - e_ref = _e_ref_lafferty(model, e_min, e_max) - else: - raise ValueError('Invalid method: {}'.format(method)) - - dnde = _dnde_from_flux(flux, model, e_ref, e_min, e_max) - - # Add to result table - table = input_table.copy() - table['e_ref'] = e_ref - table['dnde'] = dnde - - if 'flux_err' in table.colnames: - # TODO: implement better error handling, e.g. MC based method - table['dnde_err'] = dnde * table['flux_err'].quantity / flux - - if 'flux_errn' in table.colnames: - table['dnde_errn'] = dnde * table['flux_errn'].quantity / flux - table['dnde_errp'] = dnde * table['flux_errp'].quantity / flux - - if 'flux_ul' in table.colnames: - flux_ul = table['flux_ul'].quantity - dnde_ul = _dnde_from_flux(flux_ul, model, e_ref, e_min, e_max) - table['dnde_ul'] = dnde_ul - - table.meta['SED_TYPE'] = 'dnde' - return FluxPoints(table) - - -def _e_ref_lafferty(model, e_min, e_max): - # compute e_ref that the value at e_ref corresponds to the mean value - # between e_min and e_max - flux = model.integral(e_min, e_max) - dnde_mean = flux / (e_max - e_min) - return model.inverse(dnde_mean) - - -def _dnde_from_flux(flux, model, e_ref, e_min, e_max): - # Compute dnde under the assumption that flux equals expected - # flux from model - flux_model = model.integral(e_min, e_max) - dnde_model = model(e_ref) - return dnde_model * (flux / flux_model) - class FluxPointEstimator(object): - """ - Flux point estimator. + """Flux point estimator. + + Computes flux points for a given spectrum observation dataset + (a 1-dim on/off observation), energy binning and spectral model. Parameters ---------- @@ -511,17 +520,13 @@ def __init__(self, obs, groups, model): self.flux_points = None def __str__(self): - s = 'FluxPointEstimator:\n' + s = '{}:\n'.format(self.__class__.__name__) s += str(self.obs) + '\n' s += str(self.groups) + '\n' s += str(self.model) + '\n' return s def compute_points(self): - meta = OrderedDict( - METHOD='TODO', - SED_TYPE='dnde' - ) rows = [] for group in self.groups: if group.bin_type != 'normal': @@ -531,8 +536,12 @@ def compute_points(self): row = self.compute_flux_point(group) rows.append(row) - self.flux_points = FluxPoints(table_from_row_data(rows=rows, - meta=meta)) + meta = OrderedDict([ + ('method', 'TODO'), + ('SED_TYPE', 'dnde'), + ]) + table = table_from_row_data(rows=rows, meta=meta) + self.flux_points = FluxPoints(table) def compute_flux_point(self, energy_group): log.debug('Computing flux point for energy group:\n{}'.format(energy_group)) @@ -544,7 +553,9 @@ def compute_flux_point(self, energy_group): energy_ref = self.compute_energy_ref(energy_group) return self.fit_point( - model=model, energy_group=energy_group, energy_ref=energy_ref, + model=model, + energy_group=energy_group, + energy_ref=energy_ref, ) def compute_energy_ref(self, energy_group): @@ -620,24 +631,30 @@ def fit_point(self, model, energy_group, energy_ref): e_max = energy_group.energy_range.max e_min = energy_group.energy_range.min - diff_flux, diff_flux_err = res.model.evaluate_error(energy_ref) - return OrderedDict( - e_ref=energy_ref, - e_min=e_min, - e_max=e_max, - dnde=diff_flux.to('m-2 s-1 TeV-1'), - dnde_err=diff_flux_err.to('m-2 s-1 TeV-1'), - ) + dnde, dnde_err = res.model.evaluate_error(energy_ref) + return OrderedDict([ + ('e_ref', energy_ref), + ('e_min', e_min), + ('e_max', e_max), + ('dnde', dnde.to(DEFAULT_UNIT['dnde'])), + ('dnde_err', dnde_err.to(DEFAULT_UNIT['dnde'])), + ]) -class SEDLikelihoodProfile(object): - """SED likelihood profile. + +class FluxPointProfiles(object): + """Flux point likelihood profiles. See :ref:`gadf:likelihood_sed`. TODO: merge this class with the classes in ``fermipy/castro.py``, which are much more advanced / feature complete. This is just a temp solution because we don't have time for that. + + Parameters + ---------- + table : `~astropy.table.Table` + Table holding the data """ def __init__(self, table): @@ -645,20 +662,14 @@ def __init__(self, table): @classmethod def read(cls, filename, **kwargs): + """Read from file.""" filename = make_path(filename) table = Table.read(str(filename), **kwargs) return cls(table=table) - def __str__(self): - s = self.__class__.__name__ + '\n' - s += str(self.table) - return s - - def plot(self, ax=None): - import matplotlib.pyplot as plt - if ax is None: - ax = plt.gca() - # TODO + def plot(self, ax): + # TODO: copy code from fermipy, don't start from scratch! + raise NotImplementedError def chi2_flux_points(flux_points, gp_model): @@ -688,7 +699,7 @@ def chi2_flux_points_assym(flux_points, gp_model): return np.nansum(stat_per_bin), stat_per_bin -class FluxPointsFitter(object): +class FluxPointFitter(object): """ Fit a set of flux points with a parametric model. @@ -734,12 +745,14 @@ def __init__(self, stat='chi2', optimizer='simplex', error_estimator='covar', raise ValueError("'{stat}' is not a valid fit statistic, please choose" " either 'chi2' or 'chi2assym'") - if not ul_handling == 'ignore': + if ul_handling != 'ignore': raise NotImplementedError('No handling of upper limits implemented.') - self.parameters = OrderedDict(optimizer=optimizer, - error_estimator=error_estimator, - ul_handling=ul_handling) + self.parameters = OrderedDict([ + ('optimizer', optimizer), + ('error_estimator', error_estimator), + ('ul_handling', ul_handling), + ]) def _setup_sherpa_fit(self, data, model): """Fit flux point using sherpa""" @@ -838,12 +851,14 @@ def run(self, data, model): result : `~collections.OrderedDict` Dictionary with fit results and debug output. """ - result = OrderedDict() - best_fit_model = self.fit(data, model) best_fit_model = self.estimate_errors(data, best_fit_model) - result['best_fit_model'] = best_fit_model - result['dof'] = self.dof(data, model) - result['statval'] = self.statval(data, model)[0] - result['statval/dof'] = result['statval'] / result['dof'] - return result + dof = self.dof(data, model) + statval = self.statval(data, model)[0] + + return OrderedDict([ + ('best_fit_model', best_fit_model), + ('dof', dof), + ('statval', statval), + ('statval/dof', statval / dof), + ]) diff --git a/gammapy/spectrum/sherpa_chi2asym.py b/gammapy/spectrum/sherpa_chi2asym.py deleted file mode 100644 index 819babddad..0000000000 --- a/gammapy/spectrum/sherpa_chi2asym.py +++ /dev/null @@ -1,64 +0,0 @@ -# Licensed under a 3-clause BSD style license - see LICENSE.rst -"""Implement asymmetric chi-square fit statistic in Sherpa. - -To load the ``chi2asym`` fit statistic in your sherpa session:: - - import sherpa_chi2asym - sherpa_chi2asym.load_chi2asym_stat() -""" -from __future__ import absolute_import, division, print_function, unicode_literals -import numpy as np - -__all__ = [ - 'check_chi2', - 'chi2asym_err_func', - 'chi2asym_stat_func', - 'load_chi2asym_stat', -] - - -def chi2asym_stat_func(data, model, staterror=None, - syserror=None, weight=None): - """Define asymmetric chi-square errors. - - TODO: reference ROOT TGraphAsymErrors and add test against ROOT result. - - To make it fit into the Sherpa scheme we do this hack: - * staterror = statistical down error - * syserror = statistical up error - """ - # The error is attached to the data point, so if model > data, - # we have to use the up error, represented by syserror - error = np.where(model > data, syserror, staterror) - - chi = ((data - model) / error) # Chi per bin - chi2 = chi ** 2 # Chi^2 per bin - return chi2.sum(), chi2 - - -def chi2asym_err_func(data): - """Compute statistical error per bin from the data.""" - error = np.ones_like(data) - return error - - -def check_chi2(): - """Execute this function after fitting to see if the - best-fit chi2 reported matches the formula coded here""" - import sherpa.astro.ui as sau - chi2 = sau.get_fit_results().statval - print('chi2 from fit: {}'.format(chi2)) - data = sau.get_dep() - model = sau.get_model_plot().y - error = np.where(model > data, sau.get_syserror(), sau.get_staterror()) - - chi = ((data - model) / error) # Chi per bin - chi2 = chi ** 2 # Chi^2 per bin - print('chi2 re-computed: {}'.format(chi2.sum())) - - -def load_chi2asym_stat(): - """"Load and set the chi2asym statistic""" - import sherpa.astro.ui as sau - sau.load_user_stat("chi2asym", chi2asym_stat_func, chi2asym_err_func) - sau.set_stat(chi2asym) diff --git a/gammapy/spectrum/tests/test_flux_point.py b/gammapy/spectrum/tests/test_flux_point.py index 2b428d8948..193e4859f7 100644 --- a/gammapy/spectrum/tests/test_flux_point.py +++ b/gammapy/spectrum/tests/test_flux_point.py @@ -10,18 +10,19 @@ from ...utils.modeling import ParameterList from ...spectrum import SpectrumResult, SpectrumFit from ...spectrum.models import PowerLaw, SpectralModel -from ..flux_point import (_e_ref_lafferty, _dnde_from_flux, - compute_flux_points_dnde, - FluxPointEstimator, FluxPoints, FluxPointsFitter) -from ..flux_point import SEDLikelihoodProfile +from ..flux_point import FluxPoints +from ..flux_point import FluxPointProfiles +from ..flux_point import FluxPointFitter +from ..flux_point import FluxPointEstimator +from .test_energy_group import seg, obs -E_REF_METHODS = ['table', 'lafferty', 'log_center'] -indices = [0, 1, 2, 3] -FLUX_POINTS_FILES = ['diff_flux_points.ecsv', - 'diff_flux_points.fits', - 'flux_points.ecsv', - 'flux_points.fits'] +FLUX_POINTS_FILES = [ + 'diff_flux_points.ecsv', + 'diff_flux_points.fits', + 'flux_points.ecsv', + 'flux_points.fits', +] class LWTestModel(SpectralModel): @@ -80,7 +81,7 @@ def test_e_ref_lafferty(): model = LWTestModel() e_min = np.array([0.0, 0.1, 0.3, 0.6]) e_max = np.array([0.1, 0.3, 0.6, 1.0]) - actual = _e_ref_lafferty(model, e_min, e_max) + actual = FluxPoints._e_ref_lafferty(model, e_min, e_max) assert_allclose(actual, desired, atol=1e-3) @@ -94,8 +95,8 @@ def test_dnde_from_flux(): # Get values model = XSqrTestModel() - e_ref = _e_ref_lafferty(model, e_min, e_max) - dnde = _dnde_from_flux(flux, model, e_ref, e_min, e_max) + e_ref = FluxPoints._e_ref_lafferty(model, e_min, e_max) + dnde = FluxPoints._dnde_from_flux(flux, model, e_ref, e_min, e_max) # Set up test case comparison dnde_model = model(e_ref) @@ -109,7 +110,7 @@ def test_dnde_from_flux(): @requires_dependency('scipy') -@pytest.mark.parametrize('method', E_REF_METHODS) +@pytest.mark.parametrize('method', ['table', 'lafferty', 'log_center']) def test_compute_flux_points_dnde_exp(method): """ Tests against analytical result or result from gammapy.spectrum.powerlaw. @@ -118,7 +119,6 @@ def test_compute_flux_points_dnde_exp(method): e_min = [1.0, 10.0] * u.TeV e_max = [10.0, 100.0] * u.TeV - spectral_index = 2.0 table = Table() table.meta['SED_TYPE'] = 'flux' @@ -134,9 +134,9 @@ def test_compute_flux_points_dnde_exp(method): e_ref = [2.0, 20.0] * u.TeV table['e_ref'] = e_ref elif method == 'lafferty': - e_ref = _e_ref_lafferty(model, e_min, e_max) + e_ref = FluxPoints._e_ref_lafferty(model, e_min, e_max) - result = compute_flux_points_dnde(FluxPoints(table), model, method) + result = FluxPoints(table).to_sed_type('dnde', model=model, method=method) # Test energy actual = result.e_ref @@ -148,39 +148,35 @@ def test_compute_flux_points_dnde_exp(method): assert_quantity_allclose(actual, desired, rtol=1e-8) -def get_test_cases(): - # TODO: add ECPL model - try: - from .test_energy_group import seg, obs - test_cases = [] - test_cases.append( - dict(model=PowerLaw(index=Quantity(2, ''), - amplitude=Quantity(1e-11, 'm-2 s-1 TeV-1'), - reference=Quantity(1, 'TeV')), - obs=obs(), - seg=seg(obs()), - dnde=2.7465439050126e-11 * u.Unit('cm-2 s-1 TeV-1'), - dnde_err=4.755502901867284e-12 * u.Unit('cm-2 s-1 TeV-1'), - res=-0.11262182922477647, - res_err=0.1536450758523701, - ) - ) - return test_cases - except IOError: - return [] - - @requires_data('gammapy-extra') @requires_dependency('sherpa') @requires_dependency('matplotlib') @requires_dependency('scipy') -@pytest.mark.parametrize('config', get_test_cases()) +@pytest.mark.parametrize('config', ['pl']) def test_flux_points(config): + if config == 'pl': + config = dict( + model=PowerLaw( + index=Quantity(2, ''), + amplitude=Quantity(1e-11, 'm-2 s-1 TeV-1'), + reference=Quantity(1, 'TeV') + ), + obs=obs(), + seg=seg(obs()), + dnde=2.7465439050126e-11 * u.Unit('cm-2 s-1 TeV-1'), + dnde_err=4.755502901867284e-12 * u.Unit('cm-2 s-1 TeV-1'), + res=-0.11262182922477647, + res_err=0.1536450758523701, + ) + tester = FluxPointTester(config) tester.test_all() -@pytest.mark.parametrize('case', get_test_cases()) +@requires_data('gammapy-extra') +@requires_dependency('sherpa') +@requires_dependency('matplotlib') +@requires_dependency('scipy') class FluxPointTester: def __init__(self, config): self.config = config @@ -194,7 +190,8 @@ def setup(self): self.fpe = FluxPointEstimator( obs=self.config['obs'], groups=self.config['seg'].groups, - model=self.best_fit_model) + model=self.best_fit_model, + ) self.fpe.compute_points() def test_all(self): @@ -208,10 +205,11 @@ def test_basic(self): def test_approx_model(self): approx_model = self.fpe.compute_approx_model( - self.config['model'], self.fpe.groups[3]) - assert approx_model.parameters['index'].frozen == True - assert approx_model.parameters['amplitude'].frozen == False - assert approx_model.parameters['reference'].frozen == True + self.config['model'], self.fpe.groups[3] + ) + assert approx_model.parameters['index'].frozen is True + assert approx_model.parameters['amplitude'].frozen is False + assert approx_model.parameters['reference'].frozen is True def test_values(self): flux_points = self.fpe.flux_points @@ -225,8 +223,10 @@ def test_values(self): assert_quantity_allclose(actual, desired) def test_spectrum_result(self): - result = SpectrumResult(model=self.best_fit_model, - points=self.fpe.flux_points) + result = SpectrumResult( + model=self.best_fit_model, + points=self.fpe.flux_points, + ) actual = result.flux_point_residuals[0][0] desired = self.config['res'] @@ -240,14 +240,12 @@ def test_spectrum_result(self): @requires_data('gammapy-extra') -class TestSEDLikelihoodProfile: +class TestFluxPointProfiles: def setup(self): filename = '$GAMMAPY_EXTRA/datasets/spectrum/llsed_hights.fits' - self.sed = SEDLikelihoodProfile.read(filename) - - def test_basics(self): - assert 'SEDLikelihoodProfile' in str(self.sed) + self.sed = FluxPointProfiles.read(filename) + @pytest.mark.xfail @requires_dependency('matplotlib') def test_plot(self): self.sed.plot() @@ -262,11 +260,6 @@ def flux_points(request): @requires_dependency('yaml') @requires_data('gammapy-extra') class TestFluxPoints: - - @requires_dependency('matplotlib') - def test_plot(self, flux_points): - flux_points.plot() - def test_info(self, flux_points): info = str(flux_points) assert flux_points.sed_type in info @@ -316,6 +309,10 @@ def test_stack(self, flux_points): assert len(stacked.table) == 2 * len(flux_points.table) assert stacked.sed_type == flux_points.sed_type + @requires_dependency('matplotlib') + def test_plot(self, flux_points): + flux_points.plot() + @requires_data('gammapy-extra') def test_compute_flux_points_dnde(): @@ -328,7 +325,7 @@ def test_compute_flux_points_dnde(): # TODO: verify index=2.2, but it seems to give reasonable values model = PowerLaw(2.2 * u.Unit(''), 1e-12 * u.Unit('cm-2 s-1 TeV-1'), 1 * u.TeV) - actual_fp = compute_flux_points_dnde(flux_points, model=model, method='log_center') + actual_fp = flux_points.to_sed_type('dnde', model=model, method='log_center') for column in ['dnde', 'dnde_err', 'dnde_ul']: actual = actual_fp.table[column].quantity @@ -338,19 +335,19 @@ def test_compute_flux_points_dnde(): @requires_data('gammapy-extra') @requires_dependency('sherpa') -class TestFluxPointsFitter: +class TestFluxPointFitter: def setup(self): path = '$GAMMAPY_EXTRA/test_datasets/spectrum/flux_points/diff_flux_points.fits' self.flux_points = FluxPoints.read(path) def test_fit_pwl(self): - fitter = FluxPointsFitter() + fitter = FluxPointFitter() model = PowerLaw(2.3 * u.Unit(''), 1e-12 * u.Unit('cm-2 s-1 TeV-1'), 1 * u.TeV) result = fitter.run(self.flux_points, model) index = result['best_fit_model'].parameters['index'] - amplitude = result['best_fit_model'].parameters['amplitude'] assert_quantity_allclose(index.quantity, 2.216 * u.Unit(''), rtol=1e-3) + amplitude = result['best_fit_model'].parameters['amplitude'] assert_quantity_allclose(amplitude.quantity, 2.149E-13 * u.Unit('cm-2 s-1 TeV-1'), rtol=1e-3) assert_allclose(result['statval'], 27.183618, rtol=1e-3) assert_allclose(result['dof'], 22)