Skip to content

Test that PDSignalEnergyPDF splines integrate to 1 and NGC 1068 unblinding result #302

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

Open
tomaskontrimas wants to merge 1 commit into
base: master
Choose a base branch
from
+37 −0
Open

Test that PDSignalEnergyPDF splines integrate to 1 and NGC 1068 unblinding result #302

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
37 changes: 37 additions & 0 deletions tests/publicdata_ps/test_time_integrated_ps.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
import unittest

import numpy as np
from scipy import integrate

from skyllh.analyses.i3.publicdata_ps.time_integrated_ps import create_analysis
from skyllh.core.config import Config
Expand Down Expand Up @@ -54,6 +55,26 @@ def setUpClass(cls):
def setUp(self):
self.tl = TimeLord()

def test_change_source_and_unblind(self):
source = PointLikeSource(name='NGC 1068', ra=np.radians(40.67), dec=np.radians(-0.01), weight=1)
self.ana.change_source(source)
try:
rss = RandomStateService(seed=1)
with self.tl.task_timer('Call unblind for NGC 1068.'):
(TS, params_dict, status) = self.ana.unblind(minimizer_rss=rss, tl=self.tl)

logger.info(f'{self.tl}')
logger.info(f'TS = {TS:.7f}')
logger.info(f'ns_fit = {params_dict["ns"]:.7f}')
logger.info(f'gamma_fit = {params_dict["gamma"]:.7f}')
logger.info(f'minimizer status = {status}')

np.testing.assert_allclose(TS, 19.5117469, rtol=1e-5)
np.testing.assert_allclose(params_dict['ns'], 54.9606574, rtol=1e-5)
np.testing.assert_allclose(params_dict['gamma'], 3.0884301, rtol=1e-5)
Comment on lines +58 to +74
Copy link

Copilot AI Apr 16, 2026

Choose a reason for hiding this comment

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

This test adds a second full unblind() call on the same (expensive-to-construct) analysis. Given CI runs the full test suite across multiple Python versions, this can noticeably increase runtime. Consider consolidating into a single test (e.g., loop over the two sources with subTest) or otherwise marking this as a slow/integration check to avoid unnecessarily extending CI time.

Copilot uses AI. Check for mistakes.
Copy link
Copy Markdown
Collaborator Author

@tomaskontrimas tomaskontrimas Apr 16, 2026

Choose a reason for hiding this comment

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

It is an integration test and I reuse the already constructed analysis instance

finally:
self.ana.change_source(self.__class__.source)

def test_unblind(self):
rss = RandomStateService(seed=1)
with self.tl.task_timer('Call unblind.'):
Expand All @@ -69,6 +90,22 @@ def test_unblind(self):
np.testing.assert_allclose(params_dict['ns'], 12.879558, rtol=1e-5)
np.testing.assert_allclose(params_dict['gamma'], 2.122526, rtol=1e-5)

def test_signal_energy_pdf_integral_normalization(self):
"""Each PDSignalEnergyPDF spline integrates to 1 over its energy range."""
for pdfratio in self.ana._pdfratio_list:
for pdf in pdfratio.pdfratio2.sig_pdf_set.values():
integral = (
integrate.quad(
pdf.f_e_spl.evaluate,
pdf.log10_reco_e_min,
pdf.log10_reco_e_max,
limit=200,
full_output=1,
)[0]
/ pdf.f_e_spl.norm
)
np.testing.assert_allclose(integral, 1.0)
Copy link

Copilot AI Apr 16, 2026

Choose a reason for hiding this comment

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

np.testing.assert_allclose(integral, 1.0) uses the default (very tight) tolerances, which can make this normalization check flaky across platforms/BLAS/SciPy versions—especially if f_e_spl.norm is computed via a different integration scheme than quad. Use an explicit tolerance (rtol/atol) that reflects the expected numerical accuracy, and consider reducing the amount of quad work (e.g., check only a small representative subset of PDFs per dataset) to keep test runtime under control.

Suggested change
np.testing.assert_allclose(integral, 1.0)
np.testing.assert_allclose(integral, 1.0, rtol=1e-6, atol=1e-8)

Copilot uses AI. Check for mistakes.
Copy link
Copy Markdown
Collaborator Author

@tomaskontrimas tomaskontrimas Apr 16, 2026

Choose a reason for hiding this comment

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

By default it's rtol=1e-07, atol=0 and I want it to fail if that's not achieved


def test_do_trial(self):
rss = RandomStateService(seed=1)
with self.tl.task_timer('Call do_trial.'):
Expand Down