Skip to content

Commit

Permalink
some bugs+docs fix
Browse files Browse the repository at this point in the history
  • Loading branch information
Jouvin committed Oct 13, 2019
1 parent dc87010 commit eb5a195
Show file tree
Hide file tree
Showing 3 changed files with 92 additions and 83 deletions.
167 changes: 88 additions & 79 deletions gammapy/modeling/models/spectral.py
View file
Expand Up @@ -188,14 +188,14 @@ def f(x):
return self._parse_uarray(uarray) * unit

def plot(
self,
energy_range,
ax=None,
energy_unit="TeV",
flux_unit="cm-2 s-1 TeV-1",
energy_power=0,
n_points=100,
**kwargs,
self,
energy_range,
ax=None,
energy_unit="TeV",
flux_unit="cm-2 s-1 TeV-1",
energy_power=0,
n_points=100,
**kwargs,
):
"""Plot spectral model curve.
Expand Down Expand Up @@ -249,14 +249,14 @@ def plot(
return ax

def plot_error(
self,
energy_range,
ax=None,
energy_unit="TeV",
flux_unit="cm-2 s-1 TeV-1",
energy_power=0,
n_points=100,
**kwargs,
self,
energy_range,
ax=None,
energy_unit="TeV",
flux_unit="cm-2 s-1 TeV-1",
energy_power=0,
n_points=100,
**kwargs,
):
"""Plot spectral model error band.
Expand Down Expand Up @@ -380,7 +380,6 @@ def inverse(self, value, emin=0.1 * u.TeV, emax=100 * u.TeV):

energies = []
for val in np.atleast_1d(value):

def f(x):
# scale by 1e12 to achieve better precision
energy = u.Quantity(x, eunit, copy=False)
Expand Down Expand Up @@ -432,7 +431,7 @@ def __init__(self, model1, model2, operator):
self.model2 = model2
self.operator = operator
parameters = (
self.model1.parameters.parameters + self.model2.parameters.parameters
self.model1.parameters.parameters + self.model2.parameters.parameters
)
super().__init__(parameters)

Expand Down Expand Up @@ -683,7 +682,7 @@ class PowerLaw2SpectralModel(SpectralModel):
tag = "PowerLaw2SpectralModel"

def __init__(
self, amplitude="1e-12 cm-2 s-1", index=2, emin="0.1 TeV", emax="100 TeV"
self, amplitude="1e-12 cm-2 s-1", index=2, emin="0.1 TeV", emax="100 TeV"
):
self.amplitude = Parameter("amplitude", amplitude)
self.index = Parameter("index", index)
Expand Down Expand Up @@ -800,11 +799,11 @@ class ExpCutoffPowerLawSpectralModel(SpectralModel):
tag = "ExpCutoffPowerLawSpectralModel"

def __init__(
self,
index=1.5,
amplitude="1e-12 cm-2 s-1 TeV-1",
reference="1 TeV",
lambda_="0.1 TeV-1",
self,
index=1.5,
amplitude="1e-12 cm-2 s-1 TeV-1",
reference="1 TeV",
lambda_="0.1 TeV-1",
):
self.index = Parameter("index", index)
self.amplitude = Parameter("amplitude", amplitude)
Expand Down Expand Up @@ -868,11 +867,11 @@ class ExpCutoffPowerLaw3FGLSpectralModel(SpectralModel):
tag = "ExpCutoffPowerLaw3FGLSpectralModel"

def __init__(
self,
index=1.5,
amplitude="1e-12 cm-2 s-1 TeV-1",
reference="1 TeV",
ecut="10 TeV",
self,
index=1.5,
amplitude="1e-12 cm-2 s-1 TeV-1",
reference="1 TeV",
ecut="10 TeV",
):
self.index = Parameter("index", index)
self.amplitude = Parameter("amplitude", amplitude)
Expand Down Expand Up @@ -920,12 +919,12 @@ class SuperExpCutoffPowerLaw3FGLSpectralModel(SpectralModel):
tag = "SuperExpCutoffPowerLaw3FGLSpectralModel"

def __init__(
self,
index_1=1.5,
index_2=2,
amplitude="1e-12 cm-2 s-1 TeV-1",
reference="1 TeV",
ecut="10 TeV",
self,
index_1=1.5,
index_2=2,
amplitude="1e-12 cm-2 s-1 TeV-1",
reference="1 TeV",
ecut="10 TeV",
):
self.index_1 = Parameter("index_1", index_1)
self.index_2 = Parameter("index_2", index_2)
Expand Down Expand Up @@ -982,12 +981,12 @@ class SuperExpCutoffPowerLaw4FGLSpectralModel(SpectralModel):
tag = "SuperExpCutoffPowerLaw4FGLSpectralModel"

def __init__(
self,
index_1=1.5,
index_2=2,
amplitude="1e-12 cm-2 s-1 TeV-1",
reference="1 TeV",
expfactor="1e-2",
self,
index_1=1.5,
index_2=2,
amplitude="1e-12 cm-2 s-1 TeV-1",
reference="1 TeV",
expfactor="1e-2",
):
self.index_1 = Parameter("index_1", index_1)
self.index_2 = Parameter("index_2", index_2)
Expand Down Expand Up @@ -1049,7 +1048,7 @@ class LogParabolaSpectralModel(SpectralModel):
tag = "LogParabolaSpectralModel"

def __init__(
self, amplitude="1e-12 cm-2 s-1 TeV-1", reference="10 TeV", alpha=2, beta=1
self, amplitude="1e-12 cm-2 s-1 TeV-1", reference="10 TeV", alpha=2, beta=1
):
self.amplitude = Parameter("amplitude", amplitude)
self.reference = Parameter("reference", reference, frozen=True)
Expand Down Expand Up @@ -1125,14 +1124,14 @@ class TemplateSpectralModel(SpectralModel):
tag = "TemplateSpectralModel"

def __init__(
self,
energy,
values,
norm=1,
tilt=0,
reference="1 TeV",
interp_kwargs=None,
meta=None,
self,
energy,
values,
norm=1,
tilt=0,
reference="1 TeV",
interp_kwargs=None,
meta=None,
):
self.norm = Parameter("norm", norm, unit="")
self.tilt = Parameter("tilt", tilt, unit="", frozen=True)
Expand Down Expand Up @@ -1410,7 +1409,22 @@ def evaluate(self, energy, parameter):


class AbsorbedSpectralModel(SpectralModel):
"""Spectral model with EBL absorption.
r"""Spectral model with EBL absorption.
The absorption of the EBL is described as:
.. math::
\exp{-\alpha \times \tau(E,z)}
where tau(E,z) is the optical depth predicted by the model (`~gammapy.modeling.models.Absorption`), which depends
on the energy of the gamma-rays and the redshift z of the source. The class `~gammapy.modeling.models.Absorption`
computes this correction
.. math:: \exp{- \tau(E,z)}
With the optical depth scaled by a factor alpha, the observed spectrum is formed as
.. math::
\mathrm{spectral \, model} \times \exp{-\alpha \times \tau(E,z)}
Parameters
----------
Expand All @@ -1424,21 +1438,18 @@ class AbsorbedSpectralModel(SpectralModel):
parameter name
alpha_norm: float
Norm of the EBL model, by default 1.
alpha_norm_frozen: bool
False if you want to have it as a free parameter
"""

__slots__ = ["spectral_model", "absorption", "parameter", "parameter_name"]
tag = "AbsorbedSpectralModel"

def __init__(
self,
spectral_model,
absorption,
parameter,
parameter_name="redshift",
alpha_norm=1.0,
alpha_norm_frozen=True,
self,
spectral_model,
absorption,
parameter,
parameter_name="redshift",
alpha_norm=1.0,
):
self.spectral_model = spectral_model
self.absorption = absorption
Expand All @@ -1447,11 +1458,9 @@ def __init__(
min_ = self.absorption.param.min()
max_ = self.absorption.param.max()
par = Parameter(parameter_name, parameter, min=min_, max=max_, frozen=True)

self.alpha_norm = Parameter("alpha_norm", alpha_norm, frozen=True)
parameters = spectral_model.parameters.parameters.copy()
parameters.append(par)
self.alpha_norm = Parameter("alpha_norm", alpha_norm, frozen=alpha_norm_frozen)
parameters.append(self.alpha_norm)
parameters.extend([par, self.alpha_norm])

super().__init__(parameters)

Expand Down Expand Up @@ -1607,7 +1616,7 @@ class GaussianSpectralModel(SpectralModel):
tag = "GaussianSpectralModel"

def __init__(
self, norm=1e-12 * u.Unit("cm-2 s-1"), mean=1 * u.TeV, sigma=2 * u.TeV
self, norm=1e-12 * u.Unit("cm-2 s-1"), mean=1 * u.TeV, sigma=2 * u.TeV
):
self.norm = Parameter("norm", norm)
self.mean = Parameter("mean", mean)
Expand All @@ -1618,9 +1627,9 @@ def __init__(
@staticmethod
def evaluate(energy, norm, mean, sigma):
return (
norm
/ (sigma * np.sqrt(2 * np.pi))
* np.exp(-(energy - mean) ** 2 / (2 * sigma ** 2))
norm
/ (sigma * np.sqrt(2 * np.pi))
* np.exp(-(energy - mean) ** 2 / (2 * sigma ** 2))
)

def integral(self, emin, emax, **kwargs):
Expand All @@ -1638,16 +1647,16 @@ def integral(self, emin, emax, **kwargs):
# this is to get a consistent API with SpectralModel.integral()
pars = self.parameters
u_min = (
(emin - pars["mean"].quantity) / (np.sqrt(2) * pars["sigma"].quantity)
(emin - pars["mean"].quantity) / (np.sqrt(2) * pars["sigma"].quantity)
).to_value("")
u_max = (
(emax - pars["mean"].quantity) / (np.sqrt(2) * pars["sigma"].quantity)
(emax - pars["mean"].quantity) / (np.sqrt(2) * pars["sigma"].quantity)
).to_value("")

return (
pars["norm"].quantity
/ 2
* (scipy.special.erf(u_max) - scipy.special.erf(u_min))
pars["norm"].quantity
/ 2
* (scipy.special.erf(u_max) - scipy.special.erf(u_min))
)

def energy_flux(self, emin, emax):
Expand All @@ -1666,15 +1675,15 @@ def energy_flux(self, emin, emax):
"""
pars = self.parameters
u_min = (
(emin - pars["mean"].quantity) / (np.sqrt(2) * pars["sigma"].quantity)
(emin - pars["mean"].quantity) / (np.sqrt(2) * pars["sigma"].quantity)
).to_value("")
u_max = (
(emax - pars["mean"].quantity) / (np.sqrt(2) * pars["sigma"].quantity)
(emax - pars["mean"].quantity) / (np.sqrt(2) * pars["sigma"].quantity)
).to_value("")
a = pars["norm"].quantity * pars["sigma"].quantity / np.sqrt(2 * np.pi)
b = pars["norm"].quantity * pars["mean"].quantity / 2
return a * (np.exp(-u_min ** 2) - np.exp(-u_max ** 2)) + b * (
scipy.special.erf(u_max) - scipy.special.erf(u_min)
scipy.special.erf(u_max) - scipy.special.erf(u_min)
)


Expand Down Expand Up @@ -1710,7 +1719,7 @@ def __init__(self, norm=1e-12 * u.Unit("cm-2 s-1"), mean=1 * u.TeV, sigma=2):
@staticmethod
def evaluate(energy, norm, mean, sigma):
return (
norm
/ (energy * sigma * np.sqrt(2 * np.pi))
* np.exp(-(np.log(energy / mean)) ** 2 / (2 * sigma ** 2))
norm
/ (energy * sigma * np.sqrt(2 * np.pi))
* np.exp(-(np.log(energy / mean)) ** 2 / (2 * sigma ** 2))
)
4 changes: 2 additions & 2 deletions gammapy/modeling/tests/test_serialize_yaml.py
View file
Expand Up @@ -130,8 +130,8 @@ def test_datasets_to_io(tmpdir):
assert dataset0.model.skymodels[1].name == "gll_iem_v06_cutout"

assert (
dataset0.model.skymodels[0].parameters["reference"]
is dataset1.model.skymodels[1].parameters["reference"]
dataset0.model.skymodels[0].parameters["reference"]
is dataset1.model.skymodels[1].parameters["reference"]
)

assert_allclose(
Expand Down
4 changes: 2 additions & 2 deletions gammapy/spectrum/flux_point.py
View file
Expand Up @@ -1297,7 +1297,7 @@ def peek(self, method="diff/model", plot_error=True, **kwargs):
method : {"diff", "diff/model", "diff/sqrt(model)"}
Method used to compute the residuals, see `MapDataset.residuals()`
plot_error : bool
False to not plot the errors of the model since for some of them it is not handle yet properly in Gammapy
False to not plot the error band of the model
"""
from matplotlib.gridspec import GridSpec
import matplotlib.pyplot as plt
Expand Down Expand Up @@ -1398,7 +1398,7 @@ def plot_spectrum(
model_kwargs : dict
Keywords passed to `SpectralModel.plot` and `SpectralModel.plot_error`
plot_error : bool
False to not plot the errors of the model since for some of them it is not handle yet properly in Gammapy
False to not plot the error band of the model
Returns
-------
Expand Down

0 comments on commit eb5a195

Please sign in to comment.