Merge pull request #812 from adonath/rewrite_compute_differential_flu…

…x_points

Rewrite compute differential flux points method

docs/tutorials/flux_point/flux_point_demo.py

@@ -3,9 +3,9 @@
 import numpy as np
 import matplotlib.pyplot as plt
 from astropy.table import Table
-from gammapy.spectrum.flux_point import (compute_differential_flux_points,
-                                         _x_lafferty, _integrate)
-from gammapy.spectrum.powerlaw import power_law_integral_flux
+#from gammapy.spectrum.flux_point import (compute_differential_flux_points,
+#                                         _x_lafferty, _integrate)
+#from gammapy.spectrum.powerlaw import power_law_integral_flux
 
 SPECTRAL_INDEX = 4
 

docs/tutorials/flux_point/index.rst

@@ -15,7 +15,7 @@ large bin widths must be chosen to ensure statistical errors remain small. In ca
 instance, exponentially) over wide bins, simply choosing the log bin center does not offer a good representation of the true underlying spectrum.
 
 Instead, Lafferty & Wyatt [Lafferty1994]_ offer an alternative approach for such situations. In a bin of width :math:`\Delta E` between bounds
-:math:`E_1` and :math:`E_2` for energy :math:`E`, the expectation :math:`<g_{meas}>` of the true underlying spectrum, :math:`g(E)` is defined as  
+:math:`E_1` and :math:`E_2` for energy :math:`E`, the expectation :math:`<g_{meas}>` of the true underlying spectrum, :math:`g(E)` is defined as
 
 .. math::
     <g_{meas}> = \frac{1}{\Delta E}\int_{E_1}^{E_2}{g(E) dE}
@@ -25,11 +25,11 @@ which the expectation should be regarded as a measurement of the true spectrum i
 which the expectation value is equal to the mean value of the underlying true spectrum within that bin, noted as :math:`E_{lw}`. Thus knowledge of the true spectrum
 :math:`g(E)` or an estimate for this (determined by fitting) is required.
 
-So it follows that, in setting expectation equal to :math:`g(E)` at :math:`E_{lw}`, the position of :math:`E_{lw}` is given by the following equation: 
+So it follows that, in setting expectation equal to :math:`g(E)` at :math:`E_{lw}`, the position of :math:`E_{lw}` is given by the following equation:
 
 .. math::
     E_{lw} = g^{-1}\left(\frac{1}{\Delta E}\int_{E_1}^{E_2}{g(E) dE}\right)
-    
+
 For instances where a power law of spectral index 2 is taken, it can be analytically shown that the Lafferty & Wyatt method and log center method are
 coincident. In the case of steeper power laws (e.g. spectral index 3), the Lafferty & Wyatt method
 returns a lower coordinate on the energy axis than the log bin center, and the reverse effect is seen for less steep power laws (e.g. spectral index 1).
@@ -46,7 +46,7 @@ The plot demonstrates that in these cases where the true spectrum is not known,
 points offers a closer representation than the log center method. Residuals showing the percentage difference of each data point from the true
 spectrum are also shown.
 
-.. plot:: tutorials/flux_point/flux_point_demo.py plot_power_law
+.. plot:: tutorials/flux_point/flux_point_demo.py
    :include-source:
 
 Method Evaluation
@@ -57,9 +57,9 @@ As before, the underlying spectrum is a simple power law of index 4. The positio
 the bin in Log Energy, for instance the Log Center method always positions the flux point at the center of the bin for a
 Log Energy scale, and so the position is 0.5. The coincidence of the two methods for the power law case of index 2 is evident,
 with the Lafferty method positioning the flux point to the right of the Log Center position for power law spectra of index
-less than 2, and to the left for index greater than 2. 
+less than 2, and to the left for index greater than 2.
 
-.. plot:: tutorials/flux_point/x_position_plot.py make_x_plot
+.. plot:: tutorials/flux_point/x_position_plot.py
    :include-source:
 
 The ability of the two methods to represent an underlying spectrum is demonstrated in the residual images below. Here,
@@ -72,9 +72,9 @@ The residuals log ratio image (right) shows the log ratio of residuals of the La
 indicating the regions of assumed/true parameter space for which the Lafferty method offers an improved representation
 of the underlying frequency over the Log Center method (blue regions), and where it does not (red regions).
 
-.. plot:: tutorials/flux_point/residuals_images.py residuals_image
+.. plot:: tutorials/flux_point/residuals_images.py
    :include-source:
-	
+
 Caveats
 -------
 

docs/tutorials/flux_point/residuals_images.py

@@ -2,7 +2,7 @@
 """
 import numpy as np
 import matplotlib.pyplot as plt
-from gammapy.spectrum.flux_point import compute_differential_flux_points
+#from gammapy.spectrum.flux_point import compute_differential_flux_points
 from gammapy.spectrum.powerlaw import power_law_evaluate, power_law_integral_flux
 
 

docs/tutorials/flux_point/x_position_plot.py

@@ -3,7 +3,7 @@
 import numpy as np
 import matplotlib.pyplot as plt
 from astropy.table import Table
-from gammapy.spectrum.flux_point import _energy_lafferty_power_law
+#from gammapy.spectrum.flux_point import _energy_lafferty_power_law
 from gammapy.spectrum.powerlaw import power_law_evaluate, power_law_integral_flux
 from flux_point_demo import get_flux_tables
 

gammapy/spectrum/flux_point.py

@@ -17,7 +17,7 @@
 __all__ = [
     'DifferentialFluxPoints',
     'IntegralFluxPoints',
-    'compute_differential_flux_points',
+    'compute_flux_points_dnde',
     'FluxPointEstimator',
     'FluxPoints',
     'SEDLikelihoodProfile',
@@ -583,203 +583,106 @@ def to_differential_flux_points(self, x_method='lafferty',
         )
 
 
-def compute_differential_flux_points(x_method='lafferty', y_method='power_law',
-                                     table=None, model=None,
-                                     spectral_index=None, energy_min=None,
-                                     energy_max=None, int_flux=None,
-                                     int_flux_err_hi=None, int_flux_err_lo=None):
-    """Creates differential flux points table from integral flux points table.
+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
     ----------
-    table : `~astropy.table.Table`
-        Integral flux data table in energy bins, including columns
-        'ENERGY_MIN', 'ENERGY_MAX', 'INT_FLUX', 'INT_FLUX_ERR'
-    energy_min : float, array_like
-        If table not defined, minimum energy of bin(s) may be input
-        directly as either a float or array.
-    energy_max : float, array_like
-        If table not defined, maximum energy of bin(s) input directly.
-    int_flux : float, array_like
-        If table not defined, integral flux in bin(s) input directly. If array,
-        energy_min, energy_max must be either arrays of the same shape
-        (for differing energy bins) or floats (for the same energy bin).
-    int_flux_err_hi : float, array_like
-        Type must be the same as for int_flux
-    int_flux_err_lo : float, array_like
-        Type must be the same as for int_flux
-    x_method : {'lafferty', 'log_center', 'table'}
-        Flux point energy computation method; either Lafferty & Wyatt
-        model-based positioning, log bin center positioning
-        or user-defined `~astropy.table.Table` positioning
-        using column heading ['ENERGY']
-    y_method : {'power_law', 'model'}
-        Flux computation method assuming PowerLaw or user defined model function
-    model : callable
-        User-defined model function
-    spectral_index : float, array_like
-        Spectral index if default power law model is used. Either a float
-        or array_like (in which case, energy_min, energy_max and int_flux
-        must be floats to avoid ambiguity)
+    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
+
+    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)
+    >>> result = compute_flux_points_dnde(flux_points, model=model)
+
 
     Returns
     -------
-    differential_flux_table : `~astropy.table.Table`
-        Input table with appended columns 'ENERGY', 'DIFF_FLUX', 'DIFF_FLUX_ERR'
+    flux_points : `FluxPoints`
+        Flux points including differential quantity columns `dnde`
+        and `dnde_err` (optional), `dnde_ul` (optional).
 
-    Notes
-    -----
-    For usage, see this tutorial: :ref:`tutorials-flux_point`.
     """
-    # Use input values if not initially provided with a table
-    # and broadcast quantities to arrays if required
-    if table is None:
-        spectral_index = np.array(spectral_index).reshape(np.array(spectral_index).size, )
-        energy_min = np.array(energy_min).reshape(np.array(energy_min).size, )
-        energy_max = np.array(energy_max).reshape(np.array(energy_max).size, )
-        int_flux = np.array(int_flux).reshape(np.array(int_flux).size, )
-        try:
-            # TODO: unresolved reference `int_flux_err`
-            # and too broad except clause.
-            # This function is spaghetti code ... needs to be refactored!
-            int_flux_err = np.array(int_flux_err).reshape(np.array(int_flux_err).size, )
-        except:
-            pass
-        # TODO: Can a better implementation be found here?
-        lengths = dict(SPECTRAL_INDEX=len(spectral_index),
-                       ENERGY_MIN=len(energy_min),
-                       ENERGY_MAX=len(energy_max),
-                       FLUX=len(int_flux))
-        max_length = np.array(list(lengths.values())).max()
-        int_flux = np.array(int_flux) * np.ones(max_length)
-        spectral_index = np.array(spectral_index) * np.ones(max_length)
-        energy_min = np.array(energy_min) * np.ones(max_length)
-        energy_max = np.array(energy_max) * np.ones(max_length)
-        try:
-            int_flux_err = np.array(int_flux_err) * np.ones(max_length)
-        except:
-            pass
-    # Otherwise use the table provided
-    else:
-        energy_min = np.asanyarray(table['ENERGY_MIN'])
-        energy_max = np.asanyarray(table['ENERGY_MAX'])
-        int_flux = np.asanyarray(table['INT_FLUX'])
-        try:
-            int_flux_err_hi = np.asanyarray(table['INT_FLUX_ERR_HI'])
-            int_flux_err_lo = np.asanyarray(table['INT_FLUX_ERR_LO'])
-        except:
-            pass
-
-    # Compute x point
-    if x_method == 'table':
-        # This is only called if the provided table includes energies
-        energy = np.array(table['ENERGY'])
-    elif x_method == 'log_center':
-        from scipy.stats import gmean
-        energy = np.array(gmean((energy_min, energy_max)))
-    elif x_method == 'lafferty':
-        if y_method == 'power_law':
-            # Uses analytical implementation available for the power law case
-            energy = _energy_lafferty_power_law(energy_min, energy_max,
-                                                spectral_index)
-        else:
-            energy = np.array(_x_lafferty(energy_min,
-                                          energy_max, model))
-    else:
-        raise ValueError('Invalid x_method: {0}'.format(x_method))
-
-    # Compute y point
-    if y_method == 'power_law':
-        g = -1 * np.abs(spectral_index)
-        diff_flux = power_law_flux(int_flux, g, energy, energy_min, energy_max)
-    elif y_method == 'model':
-        diff_flux = _ydiff_excess_equals_expected(int_flux, energy_min,
-                                                  energy_max, energy, model)
-    else:
-        raise ValueError('Invalid y_method: {0}'.format(y_method))
-
-    # Add to table
-    table = Table()
-    table['ENERGY'] = energy
-    table['DIFF_FLUX'] = diff_flux
+    input_table = flux_points.table
+    flux = input_table['flux'].quantity
 
-    # Error processing if required
     try:
-        # TODO: more rigorous implementation of error propagation should be implemented
-        # I.e. based on MC simulation rather than gaussian error assumption
-        err_hi = int_flux_err_hi / int_flux
-        diff_flux_err_hi = err_hi * diff_flux
-        table['DIFF_FLUX_ERR_HI'] = diff_flux_err_hi
-
-        err_lo = int_flux_err_lo / int_flux
-        diff_flux_err_lo = err_lo * diff_flux
-        table['DIFF_FLUX_ERR_LO'] = diff_flux_err_lo
-    except:
-        pass
-
-    table.meta['spectral_index'] = spectral_index
-    table.meta['spectral_index_description'] = "Spectral index assumed in the DIFF_FLUX computation"
-    return table
-
-
-def _x_lafferty(xmin, xmax, function):
-    """The Lafferty & Wyatt method to compute X.
-
-    Pass in a function and bin bounds x_min and x_max i.e. for energy
-    See: Lafferty & Wyatt, Nucl. Instr. and Meth. in Phys. Res. A 355(1995) 541-547
-    See: http://nbviewer.ipython.org/gist/cdeil/bdab5f236640ef52f736
-    """
-    from scipy.optimize import brentq
-    from scipy import integrate
+        flux_err = input_table['flux_err'].quantity
+    except KeyError:
+        flux_err = None
+    try:
+        flux_ul = input_table['flux_ul'].quantity
+    except KeyError:
+        flux_ul = None
+
+    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 x_method: {0}'.format(x_method))
 
-    indices = np.arange(len(xmin))
+    dnde = _dnde_from_flux(flux, model, e_ref, e_min, e_max)
 
-    x_points = []
-    for index in indices:
-        deltax = xmax[index] - xmin[index]
-        I = integrate.quad(function, xmin[index], xmax[index], args=())
-        F = (I[0] / deltax)
+    # Add to result table
+    table = input_table.copy()
+    table['e_ref'] = e_ref
+    table['dnde'] = dnde
 
-        def g(x):
-            return function(x) - F
+    if flux_err:
+        # TODO: implement better error handling, e.g. MC based method
+        table['dnde_err'] = dnde * flux_err / flux
 
-        x_point = brentq(g, xmin[index], xmax[index])
-        x_points.append(x_point)
-    return x_points
+    if flux_ul:
+        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 _ydiff_excess_equals_expected(yint, xmin, xmax, x, model):
-    """The ExcessEqualsExpected method to compute Y (differential).
 
-    y / yint = y_model / yint_model"""
-    yint_model = _integrate(xmin, xmax, model)
-    y_model = model(x)
-    return y_model * (yint / yint_model)
+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 _integrate(xmin, xmax, function, segments=1e3):
-    """Integrates method function using the trapezium rule between xmin and xmax.
-    """
-    indices = np.arange(len(xmin))
-    y_values = []
-    for index in indices:
-        x_vals = np.arange(xmin[index], xmax[index], 1.0 / segments)
-        y_vals = function(x_vals)
-        # Division by number of segments required for correct normalization
-        y_values.append(np.trapz(y_vals) / segments)
-    return y_values
-
-
-def _energy_lafferty_power_law(energy_min, energy_max, spectral_index):
-    """Analytical case for determining lafferty x-position for power law case.
-    """
-    # Cannot call into gammapy.powerlaw as implementation is different
-    # due to different reference energies
-    term0 = 1. - spectral_index
-    term1 = energy_max - energy_min
-    term2 = 1. / term0
-    flux_lw = term2 / term1 * (energy_max ** term0 - energy_min ** term0)
-    return np.exp(-np.log(flux_lw) / np.abs(spectral_index))
+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):