From eb5a1957fe9f3cf24bf7bb4d38b9db363cf4d6a1 Mon Sep 17 00:00:00 2001 From: Jouvin Date: Sat, 12 Oct 2019 17:17:13 +0200 Subject: [PATCH] some bugs+docs fix --- gammapy/modeling/models/spectral.py | 167 +++++++++--------- gammapy/modeling/tests/test_serialize_yaml.py | 4 +- gammapy/spectrum/flux_point.py | 4 +- 3 files changed, 92 insertions(+), 83 deletions(-) diff --git a/gammapy/modeling/models/spectral.py b/gammapy/modeling/models/spectral.py index e1200de1877..a7daa4f9c63 100644 --- a/gammapy/modeling/models/spectral.py +++ b/gammapy/modeling/models/spectral.py @@ -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. @@ -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. @@ -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) @@ -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) @@ -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) @@ -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) @@ -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) @@ -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) @@ -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) @@ -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) @@ -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) @@ -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 ---------- @@ -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 @@ -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) @@ -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) @@ -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): @@ -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): @@ -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) ) @@ -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)) ) diff --git a/gammapy/modeling/tests/test_serialize_yaml.py b/gammapy/modeling/tests/test_serialize_yaml.py index b7e8285913c..2091634858c 100644 --- a/gammapy/modeling/tests/test_serialize_yaml.py +++ b/gammapy/modeling/tests/test_serialize_yaml.py @@ -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( diff --git a/gammapy/spectrum/flux_point.py b/gammapy/spectrum/flux_point.py index 4e1a8d0370a..c55103a5bb8 100644 --- a/gammapy/spectrum/flux_point.py +++ b/gammapy/spectrum/flux_point.py @@ -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 @@ -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 -------