Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Make jfactor a parameter for DarkMatter models #5097

Open
wants to merge 6 commits into
base: main
Choose a base branch
from
Open

Make jfactor a parameter for DarkMatter models #5097

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
c06570a
make jfactor/dfactor a parameter
StFroese Feb 14, 2024
95fefd4
fix test for jfactor as parameter
StFroese Feb 14, 2024
870e051
fix notebook for jfactor as parameter
StFroese Feb 14, 2024
cbe5aaa
black
StFroese Feb 14, 2024
8495cc5
simplify io
QRemy Jul 12, 2024
ebab4a6
cleanup after rebase
QRemy Jul 12, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 4 additions & 2 deletions examples/tutorials/api/astro_dark_matter.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -219,7 +219,8 @@
massDM = 10 * u.TeV
diff_flux = DarkMatterAnnihilationSpectralModel(mass=massDM, channel=channel)
int_flux = (
jfact * diff_flux.integral(energy_min=0.1 * u.TeV, energy_max=10 * u.TeV)
jfact.to_value("GeV2 cm-5")
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @StFroese, thanks a lot for all your work on this! I'm going to comment on the changes you've made.
The purpose of adding this to_value() here is just to write the Jfactor units explicitly, isn't it? It would not be needed because the jfact units are already GeV2 cm-5, am I right? Anyway, I think it is ok.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it is needed because the change of default unit in jfactor changes the dimensionality of diff_flux.

* diff_flux.integral(energy_min=0.1 * u.TeV, energy_max=10 * u.TeV)
).to("cm-2 s-1")

flux_map = WcsNDMap(geom=geom, data=int_flux.value, unit="cm-2 s-1")
Expand All @@ -243,7 +244,8 @@
massDM = 10 * u.TeV
diff_flux = DarkMatterDecaySpectralModel(mass=massDM, channel=channel)
int_flux = (
jfact_decay * diff_flux.integral(energy_min=0.1 * u.TeV, energy_max=10 * u.TeV)
jfact_decay.to_value("GeV cm-2")
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wouldn't it be better to directly rewrite the definition of jfact_decay to dfact in the tutorial notebook? Anyway, this is just a nomenclature suggestion.

* diff_flux.integral(energy_min=0.1 * u.TeV, energy_max=10 * u.TeV)
).to("cm-2 s-1")

flux_map = WcsNDMap(geom=geom, data=int_flux.value, unit="cm-2 s-1")
Expand Down
92 changes: 29 additions & 63 deletions gammapy/astro/darkmatter/spectra.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -178,8 +178,7 @@ class DarkMatterAnnihilationSpectralModel(SpectralModel):
scale : float
Scale parameter for model fitting.
jfactor : `~astropy.units.Quantity`
Integrated J-Factor needed when `~gammapy.modeling.models.PointSpatialModel`
is used.
Integrated J-Factor for model fitting.
z: float
Redshift value.
k: int
Expand All @@ -194,8 +193,7 @@ class DarkMatterAnnihilationSpectralModel(SpectralModel):

>>> channel = "b"
>>> massDM = 5000*u.Unit("GeV")
>>> jfactor = 3.41e19 * u.Unit("GeV2 cm-5")
>>> modelDM = DarkMatterAnnihilationSpectralModel(mass=massDM, channel=channel, jfactor=jfactor) # noqa: E501
>>> modelDM = DarkMatterAnnihilationSpectralModel(mass=massDM, channel=channel) # noqa: E501
Comment on lines -197 to +196
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why removing the jfactor from the example, this pattern should still work no ?


References
----------
Expand All @@ -212,22 +210,28 @@ class DarkMatterAnnihilationSpectralModel(SpectralModel):
interp="log",
)
scale._is_norm = True

jfactor = Parameter(
"jfactor",
1,
unit="GeV2 cm-5",
interp="log",
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
interp="log",
interp="log",
frozen=True,

I assume this has to be frozne by default.

)
tag = ["DarkMatterAnnihilationSpectralModel", "dm-annihilation"]

def __init__(self, mass, channel, scale=scale.quantity, jfactor=1, z=0, k=2):
def __init__(self, mass, channel, scale=scale.quantity, jfactor=jfactor, z=0, k=2):
self.k = k
self.z = z
self.mass = u.Quantity(mass)
self.channel = channel
self.jfactor = u.Quantity(jfactor)
self.primary_flux = PrimaryFlux(mass, channel=self.channel)
super().__init__(scale=scale)
super().__init__(scale=scale, jfactor=jfactor)
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree on all four changes above, along the philosophy of making the jfactor a parameter


def evaluate(self, energy, scale):
def evaluate(self, energy, scale, jfactor):
"""Evaluate dark matter annihilation model."""
flux = (
scale
* self.jfactor
* jfactor
* self.THERMAL_RELIC_CROSS_SECTION
* self.primary_flux(energy=energy * (1 + self.z))
/ self.k
Expand All @@ -242,31 +246,10 @@ def to_dict(self, full_output=False):
data = super().to_dict(full_output=full_output)
data["spectral"]["channel"] = self.channel
data["spectral"]["mass"] = self.mass.to_string()
data["spectral"]["jfactor"] = self.jfactor.to_string()
data["spectral"]["z"] = self.z
data["spectral"]["k"] = self.k
return data

@classmethod
def from_dict(cls, data):
"""Create spectral model from a dictionary.

Parameters
----------
data : dict
Dictionary with model data.

Returns
-------
model : `DarkMatterAnnihilationSpectralModel`
Dark matter annihilation spectral model.
"""
data = data["spectral"]
data.pop("type")
parameters = data.pop("parameters")
scale = [p["value"] for p in parameters if p["name"] == "scale"][0]
return cls(scale=scale, **data)


class DarkMatterDecaySpectralModel(SpectralModel):
r"""Dark matter decay spectral model.
Expand All @@ -276,7 +259,7 @@ class DarkMatterDecaySpectralModel(SpectralModel):
.. math::
\frac{\mathrm d \phi}{\mathrm d E} =
\frac{\Gamma}{4\pi m_{\mathrm{DM}}}
\frac{\mathrm d N}{\mathrm dE} \times J(\Delta\Omega)
\frac{\mathrm d N}{\mathrm dE} \times D(\Delta\Omega)

Parameters
----------
Expand All @@ -287,23 +270,21 @@ class DarkMatterDecaySpectralModel(SpectralModel):
See `PrimaryFlux.channel_registry` for more.
scale : float
Scale parameter for model fitting
jfactor : `~astropy.units.Quantity`
Integrated J-Factor needed when `~gammapy.modeling.models.PointSpatialModel`
is used.
dfactor : `~astropy.units.Quantity`
Integrated D-Factor for model fitting.
z: float
Redshift value.

Examples
--------
This is how to instantiate a `DarkMatterAnnihilationSpectralModel` model::
This is how to instantiate a `DarkMatterDecaySpectralModel` model::

>>> import astropy.units as u
>>> from gammapy.astro.darkmatter import DarkMatterDecaySpectralModel

>>> channel = "b"
>>> massDM = 5000*u.Unit("GeV")
>>> jfactor = 3.41e19 * u.Unit("GeV cm-2")
>>> modelDM = DarkMatterDecaySpectralModel(mass=massDM, channel=channel, jfactor=jfactor) # noqa: E501
>>> modelDM = DarkMatterDecaySpectralModel(mass=massDM, channel=channel) # noqa: E501
Comment on lines -305 to +287
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same here why removing the jfactor from the example ? this pattern should still work


References
----------
Expand All @@ -321,21 +302,27 @@ class DarkMatterDecaySpectralModel(SpectralModel):
)
scale._is_norm = True

dfactor = Parameter(
"dfactor",
1,
unit="GeV cm-2",
interp="log",
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
interp="log",
interp="log",
frozen=True,

)

tag = ["DarkMatterDecaySpectralModel", "dm-decay"]

def __init__(self, mass, channel, scale=scale.quantity, jfactor=1, z=0):
def __init__(self, mass, channel, scale=scale.quantity, dfactor=dfactor, z=0):
Copy link
Contributor

@QRemy QRemy Jul 12, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The mass could be a parameter frozen by default, but it could be done in a separate PR.

self.z = z
self.mass = u.Quantity(mass)
self.channel = channel
self.jfactor = u.Quantity(jfactor)
self.primary_flux = PrimaryFlux(mass, channel=self.channel)
super().__init__(scale=scale)
super().__init__(scale=scale, dfactor=dfactor)

def evaluate(self, energy, scale):
def evaluate(self, energy, scale, dfactor):
"""Evaluate dark matter decay model."""
flux = (
scale
* self.jfactor
* dfactor
* self.primary_flux(energy=energy * (1 + self.z))
/ self.LIFETIME_AGE_OF_UNIVERSE
/ self.mass
Expand All @@ -347,26 +334,5 @@ def to_dict(self, full_output=False):
data = super().to_dict(full_output=full_output)
data["spectral"]["channel"] = self.channel
data["spectral"]["mass"] = self.mass.to_string()
data["spectral"]["jfactor"] = self.jfactor.to_string()
data["spectral"]["z"] = self.z
return data

@classmethod
def from_dict(cls, data):
"""Create spectral model from dictionary.

Parameters
----------
data : dict
Dictionary with model data.

Returns
-------
model : `DarkMatterDecaySpectralModel`
Dark matter decay spectral model.
"""
data = data["spectral"]
data.pop("type")
parameters = data.pop("parameters")
scale = [p["value"] for p in parameters if p["name"] == "scale"][0]
return cls(scale=scale, **data)
14 changes: 9 additions & 5 deletions gammapy/astro/darkmatter/tests/test_spectra.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,9 @@ def test_dm_annihilation_spectral_model(tmpdir):
integral_flux = model.integral(energy_min=energy_min, energy_max=energy_max).to(
"cm-2 s-1"
)
differential_flux = model.evaluate(energy=1 * u.TeV, scale=1).to("cm-2 s-1 TeV-1")
differential_flux = model.evaluate(energy=1 * u.TeV, scale=1, jfactor=jfactor).to(
"cm-2 s-1 TeV-1"
)

sky_model = SkyModel(
spectral_model=model,
Expand All @@ -71,15 +73,17 @@ def test_dm_annihilation_spectral_model(tmpdir):
def test_dm_decay_spectral_model(tmpdir):
channel = "b"
massDM = 5 * u.TeV
jfactor = 3.41e19 * u.Unit("GeV cm-2")
dfactor = 3.41e19 * u.Unit("GeV cm-2")
energy_min = 0.01 * u.TeV
energy_max = 10 * u.TeV

model = DarkMatterDecaySpectralModel(mass=massDM, channel=channel, jfactor=jfactor)
model = DarkMatterDecaySpectralModel(mass=massDM, channel=channel, dfactor=dfactor)
integral_flux = model.integral(energy_min=energy_min, energy_max=energy_max).to(
"cm-2 s-1"
)
differential_flux = model.evaluate(energy=1 * u.TeV, scale=1).to("cm-2 s-1 TeV-1")
differential_flux = model.evaluate(energy=1 * u.TeV, scale=1, dfactor=dfactor).to(
"cm-2 s-1 TeV-1"
)

sky_model = SkyModel(
spectral_model=model,
Expand All @@ -95,6 +99,6 @@ def test_dm_decay_spectral_model(tmpdir):

assert new_models[0].spectral_model.channel == model.channel
assert new_models[0].spectral_model.z == model.z
assert_allclose(new_models[0].spectral_model.jfactor.value, model.jfactor.value)
assert_allclose(new_models[0].spectral_model.dfactor.value, model.dfactor.value)
assert new_models[0].spectral_model.mass.value == 5
assert new_models[0].spectral_model.mass.unit == u.TeV
5 changes: 3 additions & 2 deletions gammapy/astro/darkmatter/tests/test_utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ def test_dmfluxmap_annihilation(jfact_annihilation):

diff_flux = DarkMatterAnnihilationSpectralModel(mass=massDM, channel=channel)
int_flux = (
jfact_annihilation
jfact_annihilation.to_value("GeV2 cm-5")
* diff_flux.integral(energy_min=energy_min, energy_max=energy_max)
).to("cm-2 s-1")
actual = int_flux[5, 5]
Expand All @@ -63,7 +63,8 @@ def test_dmfluxmap_decay(jfact_decay):

diff_flux = DarkMatterDecaySpectralModel(mass=massDM, channel=channel)
int_flux = (
jfact_decay * diff_flux.integral(energy_min=energy_min, energy_max=energy_max)
jfact_decay.to_value("GeV cm-2")
* diff_flux.integral(energy_min=energy_min, energy_max=energy_max)
).to("cm-2 s-1")
actual = int_flux[5, 5]
desired = 7.01927e-3 / u.cm**2 / u.s
Expand Down
4 changes: 4 additions & 0 deletions gammapy/modeling/models/core.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -279,6 +279,10 @@ def from_dict(cls, data, **kwargs):
data["parameters"], cls.default_parameters
)

for key in data.keys():
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How does this replace other from_dict classmethods that have been implemented on spectral, spatial and temporal models ?
If it simply replace the dark matter from_dict maybe it is not so important to change.
Isn't there a risk that this break version compatibility if additional entries are added/removed.

Copy link
Contributor

@QRemy QRemy Jul 16, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It applies only to DM models.
I have no strong opinion about that so I can drop this commit if you prefer

if key not in ["type", "parameters"] and key not in kwargs:
kwargs[key] = data[key]

return cls.from_parameters(parameters, **kwargs)

def __str__(self):
Expand Down
Loading