Repository: gammapy/gammapy (python)

Issue:
`TemplatePhaseCurveTemporalModel` does not normalise the phasecurve template
The `TemplatePhaseCurveTemporalModel` does not normalise the template at init. This leads to big error, because the value of the norm per phase bin might be very different from 1, which means that the flux is falsely estimated, in particular when simulating data. 

Pull request description:
Norm `TemplatePhaseCurveTemporalModel`
This PR fixes #5524. 

I create a method to normalise the "NORM" of table. I need to create the interpolator once to compute the integral. 
This way, we can keep the `lazyproperty` on `_interpolator`.  I also added the possibility to use weights with the "WEIGHT" keyword in the table.

I also slightly modify the `integral` method because in the old one, if you integrate over time that correspond to less than a phase (phi_min=0, phi_max=0.4), the integral was wrong (>> than integral over 1 period). 

Patch:
diff --git a/gammapy/modeling/models/temporal.py b/gammapy/modeling/models/temporal.py
index 39988f9da..80ddca86e 100644
--- a/gammapy/modeling/models/temporal.py
+++ b/gammapy/modeling/models/temporal.py
@@ -1,5 +1,6 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 """Time-dependent models."""
+
 import logging
 import numpy as np
 import scipy.interpolate
@@ -7,8 +8,8 @@ from astropy import units as u
 from astropy.io import fits
 from astropy.table import Table
 from astropy.time import Time
-from astropy.utils import lazyproperty
 from gammapy.maps import MapAxis, RegionNDMap, TimeMapAxis
+from astropy.utils import lazyproperty
 from gammapy.modeling import Parameter
 from gammapy.utils.compat import COPY_IF_NEEDED
 from gammapy.utils.random import InverseCDFSampler, get_random_state
@@ -983,12 +984,26 @@ class TemplatePhaseCurveTemporalModel(TemporalModel):
     f2 = Parameter("f2", _f2_default, frozen=True)
 
     def __init__(self, table, filename=None, **kwargs):
-        self.table = table
+        self.table = self._normalise_table(table)
         if filename is not None:
             filename = str(make_path(filename))
         self.filename = filename
         super().__init__(**kwargs)
 
+    @staticmethod
+    def _normalise_table(table):
+        x = table["PHASE"].data
+        y = table["NORM"].data
+
+        interpolator = scipy.interpolate.InterpolatedUnivariateSpline(
+            x, y, k=1, ext=2, bbox=[0.0, 1.0]
+        )
+
+        integral = interpolator.integral(0, 1)
+
+        table["NORM"] *= 1 / integral
+        return table
+
     @classmethod
     def read(
         cls,
@@ -1010,6 +1025,7 @@ class TemplatePhaseCurveTemporalModel(TemporalModel):
             Filename with path.
         """
         filename = str(make_path(path))
+
         return cls(
             Table.read(filename),
             filename=filename,
@@ -1118,11 +1134,15 @@ class TemplatePhaseCurveTemporalModel(TemporalModel):
         ) - self._interpolator.antiderivative()(ph_min)
 
         # Divide by Jacobian (here we neglect variations of frequency during the integration period)
-        total = (phase_integral + start_integral + end_integral) / frequency
+        total = phase_integral + start_integral + end_integral
         # Normalize by total integration time
-        integral_norm = total / self.time_sum(t_min, t_max)
+        n_period = (self.time_sum(t_min, t_max) * frequency).to("")
+        if int(n_period) == 0:
+            n_period = 1
+
+        integral_norm = total / n_period
 
-        return integral_norm.to("")
+        return integral_norm
 
     @classmethod
     def from_dict(cls, data):


Test patch:
diff --git a/gammapy/modeling/models/tests/test_temporal.py b/gammapy/modeling/models/tests/test_temporal.py
index 7ce68ab44..c8e78676a 100644
--- a/gammapy/modeling/models/tests/test_temporal.py
+++ b/gammapy/modeling/models/tests/test_temporal.py
@@ -354,7 +354,6 @@ def test_sine_temporal_model_integral():
 
 @requires_data()
 def test_to_dict(light_curve):
-
     out = light_curve.to_dict()
     assert out["temporal"]["type"] == "LightCurveTemplateTemporalModel"
     assert "lightcrv_PKSB1222+216.fits" in out["temporal"]["filename"]
@@ -362,7 +361,6 @@ def test_to_dict(light_curve):
 
 @requires_data()
 def test_with_skymodel(light_curve):
-
     sky_model = SkyModel(spectral_model=PowerLawSpectralModel())
     out = sky_model.to_dict()
     assert "temporal" not in out
@@ -399,7 +397,7 @@ def test_phase_curve_model(tmp_path):
     )
 
     result = phase_model(t_ref + [0, 0.025, 0.05] * u.s)
-    assert_allclose(result, [0, 0.5, 0], atol=1e-5)
+    assert_allclose(result, [0.000000e00, 1.999989e00, 2.169609e-05], atol=1e-5)
 
     phase_model.filename = str(make_path(tmp_path / "tmp.fits"))
     phase_model.write()
@@ -415,22 +413,22 @@ def test_phase_curve_model(tmp_path):
     )
 
     assert_allclose(new_model.table["PHASE"].data, phase)
-    assert_allclose(new_model.table["NORM"].data, norm)
+    assert_allclose(new_model.table["NORM"].data, norm * 4)
 
     # exact number of phases
     integral = phase_model.integral(t_ref, t_ref + 10 * u.s)
-    assert_allclose(integral, 0.25, rtol=1e-5)
+    assert_allclose(integral, 1.0, rtol=1e-5)
     # long duration. Should be equal to the phase average
     integral = phase_model.integral(t_ref + 1 * u.h, t_ref + 3 * u.h)
-    assert_allclose(integral, 0.25, rtol=1e-5)
+    assert_allclose(integral, 1.0, rtol=1e-5)
     # 1.25 phase
     integral = phase_model.integral(t_ref, t_ref + 62.5 * u.ms)
-    assert_allclose(integral, 0.225, rtol=1e-5)
+    assert_allclose(integral, 0.9, rtol=1e-5)
 
 
 def test_phase_curve_model_sample_time():
     phase = np.linspace(0.0, 1, 51)
-    norm = 1 * (phase < 0.5)
+    norm = 1.0 * (phase < 0.5)
     table = Table(data={"PHASE": phase, "NORM": norm})
 
     t_ref = Time("2020-06-01", scale="utc")
@@ -472,7 +470,7 @@ def test_phasecurve_DC1():
     times = Time(t_ref, format="mjd") + [0.0, 0.5, 0.65, 1.0] * P0
     norm = model(times)
 
-    assert_allclose(norm, [0.05, 0.15, 1.0, 0.05])
+    assert_allclose(norm, [0.294118, 0.882353, 5.882353, 0.294118], atol=1e-5)
 
     with mpl_plot_check():
         model.plot_phasogram(n_points=200)