Repository: gammapy/gammapy
Pull request #5686 (2025-01-29 22:16:52)

Correct time sampling issue #5685
<!-- These are hidden commments. You can leave them or delete them. -->

<!-- Don't be shy! Github issues or pull requests are welcome any time. -->
<!-- But if you first want to chat, use the Gammapy Slack. -->
<!-- See https://gammapy.org/contact.html -->

**Description**
<!-- What is the motivation for this pull request? -->
<!-- Briefly: what changes are done here? -->
<!-- If this is to address a Github issue, mention its number and use the specific keywords to create a link -->
<!-- See docs.github.com/en/issues/tracking-your-work-with-issues/linking-a-pull-request-to-an-issue -->
<!-- Example: This PR is to resolve #42 -->

This pull request solves #5685 . The issue was that the time interval was split into n pixels and the model evaluated at the edges of these pixels to build the `InverseCDFSampler`. The latter therefore returned sample pixel coordinates between -0.5  and n+0.5. And all pixel coordinates below 0 and beyond n where put to 0 and n.

This solves this by sampling of pixels centers and shifting by half a pixel the sampled pixel coordinates.

**Dear reviewer**
<!-- Let the reviewer and Gammapy team know what you want: -->
<!-- * Is this ready for review? Or is it work in progress and you want some feedback? -->
<!-- * Do you want to go through review here? Or if someone should just finish this up and merge it? -->
<!-- Do you have any specific questions, e.g. about API or implementation? -->
<!-- Do you include a test executing the new code you're adding (to make sure it runs)? -->
<!-- Do you include some documentation? Is it needed? -->



Patch:
diff --git a/gammapy/modeling/models/temporal.py b/gammapy/modeling/models/temporal.py
index 80ddca86e..19892b517 100644
--- a/gammapy/modeling/models/temporal.py
+++ b/gammapy/modeling/models/temporal.py
@@ -200,18 +200,24 @@ class TemporalModel(ModelBase):
         random_state = get_random_state(random_state)
 
         ontime = (t_max - t_min).to("s")
-        n_step = (ontime / t_delta).to_value("").item()
+        n_step = np.ceil((ontime / t_delta).to_value("")).item()
         t_step = ontime / n_step
 
+        # take n+1 bins
         indices = np.arange(n_step + 1)
-        steps = indices * t_step
-        t = Time(t_min + steps, format="mjd")
 
+        # evaluate model at bin center
+        centers = (indices[:-1] + 0.5) * t_step
+        t = Time(t_min + centers, format="mjd")
         pdf = self(t)
 
+        # build sample list
         sampler = InverseCDFSampler(pdf=pdf, random_state=random_state)
         time_pix = sampler.sample(n_events)[0]
-        time = np.interp(time_pix, indices, steps)
+
+        # transform bin in time after shift by half a pixel
+        steps = indices * t_step
+        time = np.interp(time_pix + 0.5, indices, steps)
         return t_min + time
 
     def integral(self, t_min, t_max, oversampling_factor=100, **kwargs):


Test patch:
diff --git a/gammapy/datasets/tests/test_simulate.py b/gammapy/datasets/tests/test_simulate.py
index d36aeb9ff..d67ba94d1 100644
--- a/gammapy/datasets/tests/test_simulate.py
+++ b/gammapy/datasets/tests/test_simulate.py
@@ -387,7 +387,7 @@ def test_mde_sample_sources(dataset, models):
     assert_allclose(events.table["DEC_TRUE"][0], -28.748145, rtol=1e-5)
     assert events.table["DEC_TRUE"].unit == "deg"
 
-    assert_allclose(events.table["TIME"][0], 120.37471, rtol=1e-5)
+    assert_allclose(events.table["TIME"][0], 120.62471, rtol=1e-5)
     assert events.table["TIME"].unit == "s"
 
     assert_allclose(events.table["MC_ID"][0], 1, rtol=1e-5)
diff --git a/gammapy/modeling/models/tests/test_temporal.py b/gammapy/modeling/models/tests/test_temporal.py
index c8e78676a..ceae2b634 100644
--- a/gammapy/modeling/models/tests/test_temporal.py
+++ b/gammapy/modeling/models/tests/test_temporal.py
@@ -162,6 +162,17 @@ def test_time_sampling_template():
     assert_allclose(std - sigma.to("d").value, 0.0, atol=3e-4)
 
 
+def test_time_sampling_uniform():
+    cst = ConstantTemporalModel()
+    t0 = Time.now()
+    t1 = t0 + u.Quantity("100 s")
+    times = cst.sample_time(10000, t0, t1, t_delta="0.5 s").sort()
+    time_diff = (times[1:] - times[:-1]).to_value("s")
+    n_same_time = np.sum(time_diff < 1e-8)
+
+    assert n_same_time == 0
+
+
 def test_time_sampling_gaussian():
     time_ref = Time(55197.00000000, format="mjd")
     sigma = 0.5 * u.h