Skip to content

Replace integrate.quad with native spline integration method #299

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 4 commits into
base: master
Choose a base branch
from
Open

Replace integrate.quad with native spline integration method #299

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
256cd31
Replace integrate.quad with native spline integration method
tomaskontrimas Apr 13, 2026
24573f5
Remove integrity check after covering it with test
tomaskontrimas Apr 16, 2026
818530b
ruff
tomaskontrimas Apr 16, 2026
7e4cd32
Potential fix for pull request finding
tomaskontrimas Apr 16, 2026
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
14 changes: 3 additions & 11 deletions skyllh/analyses/i3/publicdata_ps/aeff.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,8 +1,5 @@
import numpy as np
from scipy import (
integrate,
interpolate,
)
from scipy import interpolate

from skyllh.analyses.i3.publicdata_ps.utils import (
FctSpline2D,
Expand Down Expand Up @@ -326,18 +323,13 @@ def get_detection_prob_for_decnu(self, decnu, enu_min, enu_max, enu_range_min, e

spl = interpolate.splrep(x, y, k=1, s=0)

def _eval_spl_func(x):
return interpolate.splev(x, spl, der=0, ext=1)

norm = integrate.quad(_eval_spl_func, enu_range_min, enu_range_max, limit=200, full_output=1)[0]
norm = interpolate.splint(enu_range_min, enu_range_max, spl)

enu_min = np.atleast_1d(enu_min)
enu_max = np.atleast_1d(enu_max)

det_prob = np.empty((len(enu_min),), dtype=np.double)
for i in range(len(enu_min)):
integral = integrate.quad(_eval_spl_func, enu_min[i], enu_max[i], limit=200, full_output=1)[0]

det_prob[i] = integral / norm
det_prob[i] = interpolate.splint(enu_min[i], enu_max[i], spl) / norm
Comment on lines 324 to +333
Copy link
Copy Markdown
Collaborator Author

@tomaskontrimas tomaskontrimas Apr 16, 2026 •
edited
Loading

Choose a reason for hiding this comment

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

From docs: splint silently assumes that the spline function is zero outside the data interval (a, b).

which I think is better than the suggested implementation and it still pass existing analysis tests


return det_prob
11 changes: 0 additions & 11 deletions skyllh/analyses/i3/publicdata_ps/signalpdf.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,4 @@
import numpy as np
from scipy import integrate

from skyllh.analyses.i3.publicdata_ps.aeff import (
PDAeff,
Expand Down Expand Up @@ -83,16 +82,6 @@ def __init__(
# Add the PDF axes.
self.add_axis(PDFAxis(name='log_energy', vmin=self.log10_reco_e_min, vmax=self.log10_reco_e_max))

Comment on lines 82 to 84
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.

Covered it in the integration test

# Check integrity.
integral = (
integrate.quad(
self.f_e_spl.evaluate, self.log10_reco_e_min, self.log10_reco_e_max, limit=200, full_output=1
)[0]
/ self.f_e_spl.norm
)
if not np.isclose(integral, 1):
raise ValueError(f'The integral over log10_reco_e of the energy term must be unity! But it is {integral}!')

def assert_is_valid_for_trial_data(self, tdm, tl=None):
pass

Expand Down
9 changes: 4 additions & 5 deletions skyllh/analyses/i3/publicdata_ps/utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,8 +1,5 @@
import numpy as np
from scipy import (
integrate,
interpolate,
)
from scipy import interpolate

from skyllh.core.binning import (
get_bincenters_from_binedges,
Expand Down Expand Up @@ -42,7 +39,9 @@ class from scipy.

self.norm = None
if norm:
self.norm = integrate.quad(self.__call__, self.x_min, self.x_max, limit=200, full_output=1)[0]
# The spline is defined only in the `x` (bincenters) interval by construction.
# We choose not to extrapolate out-of-range values.
self.norm = float(self.spl_f.integrate(x[0], x[-1]))

def __call__(self, x, oor_value=0):
"""Evaluates the spline at the given x values. For x-values
Expand Down