Repository: gammapy/gammapy
Pull request #5712 (2025-02-26 00:16:44)

Add parameter limits on phi for Spatial models
Having limits on the rotation angle for the asymmetry makes the fit more stable.
Also, adapted the warning message for parameter range to warn when the parameter falls strictly outside the range (rather than equals) because otherwise the fit starts with a warning when unfreezing the `e`.

codecov[bot]: ## [Codecov](https://app.codecov.io/gh/gammapy/gammapy/pull/5712?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 94.62%. Comparing base [(`d75be33`)](https://app.codecov.io/gh/gammapy/gammapy/commit/d75be33c1ab4d5967b2977867320befaf96da3b2?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) to head [(`dbdc91b`)](https://app.codecov.io/gh/gammapy/gammapy/commit/dbdc91bb2f8ad6fc5f323b742934786cd70083ec?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy).
> Report is 54 commits behind head on main.

<details><summary>Additional details and impacted files</summary>


```diff
@@           Coverage Diff           @@
##             main    #5712   +/-   ##
=======================================
  Coverage   94.62%   94.62%           
=======================================
  Files         239      239           
  Lines       37468    37470    +2     
=======================================
+ Hits        35454    35456    +2     
  Misses       2014     2014           
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/gammapy/gammapy/pull/5712?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy).

bkhelifi: Should we then add at least a circular [0,360]deg or [-180,180]deg? This should be better for minimizers...

Patch:
diff --git a/gammapy/modeling/models/spatial.py b/gammapy/modeling/models/spatial.py
index 51b6851a7..965176a33 100644
--- a/gammapy/modeling/models/spatial.py
+++ b/gammapy/modeling/models/spatial.py
@@ -28,6 +28,7 @@ from gammapy.utils.interpolation import interpolation_scale
 from gammapy.utils.regions import region_circle_to_ellipse, region_to_frame
 from gammapy.utils.scripts import make_path
 from .core import ModelBase, _build_parameters_from_dict
+from gammapy.utils.units import wrap_at
 
 __all__ = [
     "ConstantFluxSpatialModel",
@@ -606,10 +607,10 @@ class GaussianSpatialModel(SpatialModel):
         Length of the major semiaxis of the Gaussian, in angular units.
         Default is 1 deg.
     e : `float`
-        Eccentricity of the Gaussian (:math:`0< e< 1`).
+        Eccentricity of the Gaussian (:math:`0<= e<= 1`).
         Default is 0.
     phi : `~astropy.coordinates.Angle`
-        Rotation angle :math:`\phi`: of the major semiaxis.
+        Rotation angle :math:`\phi`: of the major semiaxis (:math:`-180 <= phi <= 180`).
         Increases counter-clockwise from the North direction.
         Default is 0 deg.
     frame : {"icrs", "galactic"}
@@ -622,7 +623,7 @@ class GaussianSpatialModel(SpatialModel):
     lat_0 = Parameter("lat_0", "0 deg", min=-90, max=90)
     sigma = Parameter("sigma", "1 deg", min=0)
     e = Parameter("e", 0, min=0, max=1, frozen=True)
-    phi = Parameter("phi", "0 deg", frozen=True)
+    phi = Parameter("phi", "0 deg", frozen=True, min=-180, max=180)
 
     @property
     def evaluation_bin_size_min(self):
@@ -641,6 +642,7 @@ class GaussianSpatialModel(SpatialModel):
     def evaluate(lon, lat, lon_0, lat_0, sigma, e, phi):
         """Evaluate model."""
         sep = angular_separation(lon, lat, lon_0, lat_0)
+        phi = wrap_at(phi, 0 * u.deg, 180 * u.deg)
 
         if e == 0:
             a = 1.0 - np.cos(sigma)
@@ -736,10 +738,10 @@ class GeneralizedGaussianSpatialModel(SpatialModel):
         Shape parameter within (0, 1). Special cases for disk: ->0, Gaussian: 0.5, Laplace:1
         Default is 0.5.
     e : `float`
-        Eccentricity (:math:`0< e< 1`).
+        Eccentricity (:math:`0<= e< =1`).
         Default is 0.
     phi : `~astropy.coordinates.Angle`
-        Rotation angle :math:`\phi`: of the major semiaxis.
+        Rotation angle :math:`\phi`: of the major semiaxis (:math:`-180<=phi<=180`).
         Increases counter-clockwise from the North direction.
         Default is 0 deg.
     frame : {"icrs", "galactic"}
@@ -752,10 +754,11 @@ class GeneralizedGaussianSpatialModel(SpatialModel):
     r_0 = Parameter("r_0", "1 deg")
     eta = Parameter("eta", 0.5, min=0.01, max=1.0)
     e = Parameter("e", 0.0, min=0.0, max=1.0, frozen=True)
-    phi = Parameter("phi", "0 deg", frozen=True)
+    phi = Parameter("phi", "0 deg", frozen=True, min=-180, max=180)
 
     @staticmethod
     def evaluate(lon, lat, lon_0, lat_0, r_0, eta, e, phi):
+        phi = wrap_at(phi, 0 * u.deg, 180 * u.deg)
         sep = angular_separation(lon, lat, lon_0, lat_0)
         if isinstance(eta, u.Quantity):
             eta = eta.value  # gamma function does not allow quantities
@@ -859,10 +862,10 @@ class DiskSpatialModel(SpatialModel):
         :math:`a`: length of the major semiaxis, in angular units.
         Default is 1 deg.
     e : `float`
-        Eccentricity of the ellipse (:math:`0< e< 1`).
+        Eccentricity of the ellipse (:math:`0<= e<= 1`).
         Default is 0.
     phi : `~astropy.coordinates.Angle`
-        Rotation angle :math:`\phi`: of the major semiaxis.
+        Rotation angle :math:`\phi`: of the major semiaxis (:math:`-180<=phi<=200`).
         Increases counter-clockwise from the North direction.
         Default is 0 deg.
     edge_width : float
@@ -879,7 +882,7 @@ class DiskSpatialModel(SpatialModel):
     lat_0 = Parameter("lat_0", "0 deg", min=-90, max=90)
     r_0 = Parameter("r_0", "1 deg", min=0)
     e = Parameter("e", 0, min=0, max=1, frozen=True)
-    phi = Parameter("phi", "0 deg", frozen=True)
+    phi = Parameter("phi", "0 deg", frozen=True, min=-180, max=180)
     edge_width = Parameter("edge_width", value=0.01, min=0, max=1, frozen=True)
 
     @property
@@ -928,6 +931,7 @@ class DiskSpatialModel(SpatialModel):
     def evaluate(lon, lat, lon_0, lat_0, r_0, e, phi, edge_width):
         """Evaluate model."""
         sep = angular_separation(lon, lat, lon_0, lat_0)
+        phi = wrap_at(phi, 0 * u.deg, 180 * u.deg)
 
         if e == 0:
             sigma_eff = r_0
diff --git a/gammapy/modeling/parameter.py b/gammapy/modeling/parameter.py
index a01078bf9..0516bd1dd 100644
--- a/gammapy/modeling/parameter.py
+++ b/gammapy/modeling/parameter.py
@@ -456,8 +456,8 @@ class Parameter:
     def check_limits(self):
         """Emit a warning or error if value is outside the minimum/maximum range."""
         if not self.frozen:
-            if (~np.isnan(self.min) and (self.value <= self.min)) or (
-                ~np.isnan(self.max) and (self.value >= self.max)
+            if (~np.isnan(self.min) and (self.value < self.min)) or (
+                ~np.isnan(self.max) and (self.value > self.max)
             ):
                 log.warning(
                     f"Value {self.value} is outside bounds [{self.min}, {self.max}]"
diff --git a/gammapy/utils/units.py b/gammapy/utils/units.py
index 285f2d037..34c98e40e 100644
--- a/gammapy/utils/units.py
+++ b/gammapy/utils/units.py
@@ -1,15 +1,41 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 """Units and Quantity related helper functions."""
+
 import logging
 from math import floor
 import numpy as np
 import astropy.units as u
 
-__all__ = ["standardise_unit", "unit_from_fits_image_hdu"]
+__all__ = ["standardise_unit", "unit_from_fits_image_hdu", "wrap_at"]
 
 log = logging.getLogger(__name__)
 
 
+def wrap_at(value, min_value, max_value):
+    """Wraps a value to be within the range [min_value, max_value).
+
+    Parameters
+    ----------
+    value : float or `~astropy.units.Unit`
+        The value to wrap
+    min_value : float or `~astropy.units.Unit`
+        The minimum value of the range (inclusive).
+    max_value : float or `~astropy.units.Unit`
+        The maximum value of the range (exclusive).
+
+    Returns
+    -------
+    wrapped_val : float or `~astropy.units.Unit`
+        Wrapped value
+    """
+
+    # Note: Astropy wrap_at currently works only for wrapping
+    # between [0, 360] and [-180, 180]
+    range_size = max_value - min_value
+    wrapped_val = ((value - min_value) % range_size) + min_value
+    return wrapped_val
+
+
 def standardise_unit(unit):
     """Standardise unit.
 


Test patch:
diff --git a/gammapy/modeling/models/tests/test_spatial.py b/gammapy/modeling/models/tests/test_spatial.py
index 11979268a..675feceac 100644
--- a/gammapy/modeling/models/tests/test_spatial.py
+++ b/gammapy/modeling/models/tests/test_spatial.py
@@ -166,6 +166,15 @@ def test_sky_disk():
     assert_allclose(radius.to_value("deg"), 2.222)
     assert_allclose(model.evaluation_bin_size_min, 0.198 * u.deg)
 
+    assert model.phi.min == -180
+    assert model.phi.max == 180
+
+    model.phi.quantity = 210 * u.deg
+    val1 = model(lon, lat)
+    model.phi.value = 30
+    val2 = model(lon, lat)
+    assert_allclose(val1, val2)
+
     # test the normalization for an elongated ellipse near the Galactic Plane
     m_geom_1 = WcsGeom.create(
         binsz=0.015, width=(20, 20), skydir=(2, 2), frame="galactic", proj="AIT"
diff --git a/gammapy/utils/tests/test_units.py b/gammapy/utils/tests/test_units.py
index 4efd390cf..93e1b6db6 100644
--- a/gammapy/utils/tests/test_units.py
+++ b/gammapy/utils/tests/test_units.py
@@ -3,7 +3,8 @@ import pytest
 import numpy as np
 import astropy.units as u
 from gammapy.maps import MapAxis
-from gammapy.utils.units import energy_unit_format, standardise_unit
+from gammapy.utils.units import energy_unit_format, standardise_unit, wrap_at
+from numpy.testing import assert_allclose
 
 
 def test_standardise_unit():
@@ -39,3 +40,12 @@ values = [
 @pytest.mark.parametrize("q, expect", values)
 def test_energy_unit_format(q, expect):
     assert energy_unit_format(q) == expect
+
+
+def test_wrap_at():
+    computed = wrap_at(6.5 * u.deg, 1.0 * u.deg, 3.5 * u.deg)
+    assert_allclose(1.5 * u.deg, computed, rtol=1e-4)
+
+    computed = wrap_at(np.arange(0, 10, 0.5), 1.0, 3.5)
+    assert_allclose(computed.min(), 1.0, rtol=1e-4)
+    assert_allclose(computed.max(), 3.0, rtol=1e-4)