Add default parameters for spectral models

gammapy/datasets/fermi.py

@@ -310,7 +310,7 @@ def counts(self):
 
         try:
             cube = SkyCube.read(filename, format='fermi-counts')
-        except ValueError:
+        except (ValueError, KeyError):
             from ..cube.healpix import SkyCubeHealpix
             cube = SkyCubeHealpix.read(filename, format='fermi-counts')
 

gammapy/image/tests/test_catalog.py

@@ -3,6 +3,7 @@
 from numpy.testing import assert_allclose
 from astropy import units as u
 from astropy.wcs import WCS
+from astropy.tests.helper import remote_data
 from ...utils.testing import requires_dependency, requires_data
 from ..catalog import CatalogImageEstimator, catalog_image, catalog_table, _source_image
 from ...image import SkyImage
@@ -18,6 +19,7 @@ def test_extended_image():
     pass
 
 
+@remote_data
 @requires_dependency('scipy')
 @requires_data('gammapy-extra')
 def test_source_image():
@@ -40,6 +42,7 @@ def test_source_image():
     assert_allclose(actual, expected)
 
 
+@remote_data
 @requires_dependency('scipy')
 @requires_data('gammapy-extra')
 def test_catalog_image():
@@ -61,6 +64,7 @@ def test_catalog_image():
     assert_allclose(actual, expected, rtol=0.01)
 
 
+@remote_data
 @requires_data('gammapy-extra')
 def test_catalog_table():
     # Checks catalogs are loaded correctly

gammapy/spectrum/flux_point.py

@@ -718,7 +718,6 @@ def fit_point(self, model, energy_group, energy_ref, sqrt_ts_threshold=1):
 
         # Set reference and remove min amplitude
         model.parameters['reference'].value = energy_ref.to('TeV').value
-        model.parameters['amplitude'].parmin = None
 
         fit = SpectrumFit(self.obs, model)
         erange = energy_group.energy_range

gammapy/spectrum/models.py

@@ -396,13 +396,27 @@ class PowerLaw(SpectralModel):
         :math:`Phi_0`
     reference : `~astropy.units.Quantity`
         :math:`E_0`
+
+
+    Examples
+    --------
+    This is how to plot the default `PowerLaw` model:
+
+        from astropy import units as u
+        from gammapy.spectrum.models import PowerLaw
+
+        pwl = PowerLaw()
+        pwl.plot(energy_range=[0.1, 100] * u.TeV)
+        plt.show()
+
     """
 
-    def __init__(self, index, amplitude, reference):
+    def __init__(self, index=2., amplitude=1E-12 * u.Unit('cm-2 s-1 TeV-1'),
+                 reference=1 * u.TeV):
         self.parameters = ParameterList([
-            Parameter('index', index, parmin=0),
-            Parameter('amplitude', amplitude, parmin=0),
-            Parameter('reference', reference, frozen=True)
+            Parameter('index', index),
+            Parameter('amplitude', amplitude),
+            Parameter('reference', reference, parmin=0, frozen=True)
         ])
 
     @staticmethod
@@ -585,14 +599,26 @@ class PowerLaw2(SpectralModel):
         Lower energy limit :math:`E_{0, min}`.
     emax : `~astropy.units.Quantity`
         Upper energy limit :math:`E_{0, max}`.
+
+    Examples
+    --------
+    This is how to plot the default `PowerLaw2` model:
+
+        from astropy import units as u
+        from gammapy.spectrum.models import PowerLaw2
+
+        pwl2 = PowerLaw2()
+        pwl2.plot(energy_range=[0.1, 100] * u.TeV)
+        plt.show()
     """
 
-    def __init__(self, amplitude, index, emin, emax):
+    def __init__(self, amplitude=1E-12 * u.Unit('cm-2 s-1'), index=2,
+                 emin=0.1 * u.TeV, emax=100 * u.TeV):
         self.parameters = ParameterList([
-            Parameter('amplitude', amplitude, parmin=0),
-            Parameter('index', index, parmin=0),
-            Parameter('emin', emin),
-            Parameter('emax', emax)
+            Parameter('amplitude', amplitude),
+            Parameter('index', index),
+            Parameter('emin', emin, frozen=True),
+            Parameter('emax', emax, frozen=True)
         ])
 
     @staticmethod
@@ -692,14 +718,26 @@ class ExponentialCutoffPowerLaw(SpectralModel):
         :math:`E_0`
     lambda : `~astropy.units.Quantity`
         :math:`\lambda`
+
+    Examples
+    --------
+    This is how to plot the default `ExponentialCutoffPowerLaw` model:
+
+        from astropy import units as u
+        from gammapy.spectrum.models import ExponentialCutoffPowerLaw
+
+        ecpl = ExponentialCutoffPowerLaw()
+        ecpl.plot(energy_range=[0.1, 100] * u.TeV)
+        plt.show()
     """
 
-    def __init__(self, index, amplitude, reference, lambda_):
+    def __init__(self, index=1.5, amplitude=1E-12 * u.Unit('cm-2 s-1 TeV-1'),
+                 reference=1 * u.TeV, lambda_=0.1 / u.TeV):
         self.parameters = ParameterList([
-            Parameter('index', index, parmin=0),
-            Parameter('amplitude', amplitude, parmin=0),
+            Parameter('index', index),
+            Parameter('amplitude', amplitude),
             Parameter('reference', reference, frozen=True),
-            Parameter('lambda_', lambda_, parmin=0)
+            Parameter('lambda_', lambda_)
         ])
 
     @staticmethod
@@ -754,13 +792,25 @@ class ExponentialCutoffPowerLaw3FGL(SpectralModel):
         :math:`E_0`
     ecut : `~astropy.units.Quantity`
         :math:`E_{C}`
+
+    Examples
+    --------
+    This is how to plot the default `ExponentialCutoffPowerLaw3FGL` model:
+
+        from astropy import units as u
+        from gammapy.spectrum.models import ExponentialCutoffPowerLaw3FGL
+
+        ecpl_3fgl = ExponentialCutoffPowerLaw3FGL()
+        ecpl_3fgl.plot(energy_range=[0.1, 100] * u.TeV)
+        plt.show()
     """
 
-    def __init__(self, index, amplitude, reference, ecut):
+    def __init__(self, index=1.5, amplitude=1E-12 * u.Unit('cm-2 s-1 TeV-1'),
+                 reference=1 * u.TeV, ecut=10 * u.TeV):
         self.parameters = ParameterList([
-            Parameter('index', index, parmin=0),
-            Parameter('amplitude', amplitude, parmin=0),
-            Parameter('reference', reference, frozen=0),
+            Parameter('index', index),
+            Parameter('amplitude', amplitude),
+            Parameter('reference', reference, frozen=True),
             Parameter('ecut', ecut)
         ])
 
@@ -798,16 +848,28 @@ class PLSuperExpCutoff3FGL(SpectralModel):
         :math:`E_0`
     ecut : `~astropy.units.Quantity`
         :math:`E_{C}`
+
+    Examples
+    --------
+    This is how to plot the default `PLSuperExpCutoff3FGL` model:
+
+        from astropy import units as u
+        from gammapy.spectrum.models import PLSuperExpCutoff3FGL
+
+        secpl_3fgl = PLSuperExpCutoff3FGL()
+        secpl_3fgl.plot(energy_range=[0.1, 100] * u.TeV)
+        plt.show()
     """
 
-    def __init__(self, index_1, index_2, amplitude, reference, ecut):
+    def __init__(self, index_1=1.5, index_2=2, amplitude=1E-12 * u.Unit('cm-2 s-1 TeV-1'),
+                 reference=1 * u.TeV, ecut=10 * u.TeV):
         # TODO: order or parameters is different from argument list / docstring. Make uniform!
         self.parameters = ParameterList([
-            Parameter('amplitude', amplitude, parmin=0),
-            Parameter('reference', reference, frozen=0),
+            Parameter('amplitude', amplitude),
+            Parameter('reference', reference, frozen=True),
             Parameter('ecut', ecut),
-            Parameter('index_1', index_1, parmin=0),
-            Parameter('index_2', index_2, parmin=0),
+            Parameter('index_1', index_1),
+            Parameter('index_2', index_2),
         ])
 
     @staticmethod
@@ -843,11 +905,23 @@ class LogParabola(SpectralModel):
         :math:`\alpha`
     beta : `~astropy.units.Quantity`
         :math:`\beta`
+
+    Examples
+    --------
+    This is how to plot the default `LogParabola` model:
+
+        from astropy import units as u
+        from gammapy.spectrum.models import LogParabola
+
+        log_parabola = LogParabola()
+        log_parabola.plot(energy_range=[0.1, 100] * u.TeV)
+        plt.show()
     """
 
-    def __init__(self, amplitude, reference, alpha, beta):
+    def __init__(self, amplitude=1E-12 * u.Unit('cm-2 s-1 TeV-1'), reference=10 * u.TeV,
+                 alpha=2, beta=1):
         self.parameters = ParameterList([
-            Parameter('amplitude', amplitude, parmin=0),
+            Parameter('amplitude', amplitude),
             Parameter('reference', reference, frozen=True),
             Parameter('alpha', alpha),
             Parameter('beta', beta)
@@ -872,7 +946,7 @@ def evaluate(energy, amplitude, reference, alpha, beta):
 class TableModel(SpectralModel):
     """A model generated from a table of energy and value arrays.
 
-    The units returned will be the units of the values array provided at
+    the units returned will be the units of the values array provided at
     initialization. The model will return values interpolated in
     log-space, returning 0 for energies outside of the limits of the provided
     energy array.

gammapy/utils/modeling.py

@@ -85,14 +85,14 @@ class Parameter(object):
     def __init__(self, name, value, unit='', parmin=None, parmax=None, frozen=False):
         self.name = name
 
-        if isinstance(value, u.Quantity):
+        if isinstance(value, u.Quantity) or isinstance(value, six.string_types):
             self.quantity = value
         else:
             self.value = value
             self.unit = unit
 
-        self.parmin = parmin
-        self.parmax = parmax
+        self.parmin = parmin or np.nan
+        self.parmax = parmax or np.nan
         self.frozen = frozen
 
     @property
@@ -102,6 +102,7 @@ def quantity(self):
 
     @quantity.setter
     def quantity(self, par):
+        par = u.Quantity(par)
         self.value = par.value
         self.unit = par.unit
 
@@ -153,11 +154,14 @@ def to_xml(self):
     def to_sherpa(self, modelname='Default'):
         """Convert to sherpa parameter"""
         from sherpa.models import Parameter
+
+        parmin = np.finfo(np.float32).min if np.isnan(self.parmin) else self.parmin
+        parmax = np.finfo(np.float32).max if np.isnan(self.parmax) else self.parmax
+
         par = Parameter(modelname=modelname, name=self.name,
                         val=self.value, units=self.unit,
-                        min=self.parmin, max=self.parmax,
+                        min=parmin, max=parmax,
                         frozen=self.frozen)
-
         return par
 
 
@@ -231,7 +235,9 @@ def to_table(self):
         """
         names = ['name', 'value', 'error', 'unit', 'min', 'max', 'frozen']
         formats = {'value': '.3e',
-                   'error': '.3e', }
+                   'error': '.3e',
+                   'min': '.3e',
+                   'max': '.3e'}
         table = Table(self.to_list_of_dict(), names=names)
 
         for name in formats: