From 11b8d2e6f3b0706ee64bc5e25f00f8a3990b743f Mon Sep 17 00:00:00 2001 From: Jakob van Santen Date: Wed, 7 May 2025 10:03:55 +0200 Subject: [PATCH] Replace interp2d with RectBivariateSpline interp2d is gone in scipy 1.14, and RecBivariateSpline is 5x faster to evaluate. While the numerical results are identical, RectBivariateSpline returns a 2D array for scalar inputs, whereas interp2d returned a 1D array. Work around this by removing all 1-entry dimensions from the result. A notable exception is the acceptance in StandardLLH, which appears to have relied on the result being 1D (i.e. tests fail if the interpolation returns a scalar). Fixes #270 --- flarestack/core/llh.py | 9 +++-- .../icecube_utils/reference_sensitivity.py | 33 +++++++++++-------- flarestack/utils/percentile_SoB.py | 4 ++- 3 files changed, 27 insertions(+), 19 deletions(-) diff --git a/flarestack/core/llh.py b/flarestack/core/llh.py index ecfb98ab..7173637f 100644 --- a/flarestack/core/llh.py +++ b/flarestack/core/llh.py @@ -826,7 +826,7 @@ def create_acceptance_function(self): with open(acc_path, "rb") as f: [dec_bins, gamma_bins, acc] = pickle.load(f) - f = scipy.interpolate.interp2d(dec_bins, gamma_bins, acc.T, kind="linear") + f = scipy.interpolate.RectBivariateSpline(dec_bins, gamma_bins, acc, kx=1, ky=1) return f def new_acceptance(self, source, params=None): @@ -845,7 +845,7 @@ def new_acceptance(self, source, params=None): dec = source["dec_rad"] gamma = params[-1] - return self.acceptance_f(dec, gamma) + return self.acceptance_f(dec, gamma)[0, :] def create_kwargs(self, data, pull_corrector, weight_f=None): kwargs = dict() @@ -1480,9 +1480,8 @@ def create_kwargs(self, data, pull_corrector, weight_f=None): coincident_data = data[coincident_nu_mask] coincident_sources = self.sources[coincident_source_mask] - season_weight = lambda x: weight_f([1.0, x], self.season)[ - coincident_source_mask - ] + def season_weight(x): + return weight_f([1.0, x], self.season)[coincident_source_mask] SoB_energy_cache = self.create_SoB_energy_cache(coincident_data) diff --git a/flarestack/icecube_utils/reference_sensitivity.py b/flarestack/icecube_utils/reference_sensitivity.py index fc14d251..a3836f18 100644 --- a/flarestack/icecube_utils/reference_sensitivity.py +++ b/flarestack/icecube_utils/reference_sensitivity.py @@ -1,8 +1,7 @@ import logging -import os import numpy as np -from scipy.interpolate import interp1d, interp2d +from scipy.interpolate import RectBivariateSpline from flarestack.data.icecube.ic_season import get_published_sens_ref_dir @@ -62,10 +61,12 @@ def reference_7year_sensitivity(sindec=np.array(0.0), gamma=2.0): sens = np.vstack((sens[0], sens)) sens = np.vstack((sens, sens[-1])) - sens_ref = interp2d(np.array(sindecs), np.array(gammas), np.log(sens.T)) + sens_ref = RectBivariateSpline( + np.array(sindecs), np.array(gammas), np.log(sens), kx=1, ky=1 + ) if np.array(sindec).ndim > 0: - return np.array([np.exp(sens_ref(x, gamma))[0] for x in sindec]) + return np.array([np.exp(sens_ref(x, gamma).squeeze()) for x in sindec]) else: return np.exp(sens_ref(sindec, gamma)) @@ -98,12 +99,14 @@ def reference_7year_discovery_potential(sindec=0.0, gamma=2.0): disc = np.vstack((disc[0], disc)) disc = np.vstack((disc, disc[-1])) - disc_ref = interp2d(np.array(sindecs), np.array(gammas), np.log(disc.T)) + disc_ref = RectBivariateSpline( + np.array(sindecs), np.array(gammas), np.log(disc), kx=1, ky=1 + ) if np.array(sindec).ndim > 0: - return np.array([np.exp(disc_ref(x, gamma))[0] for x in sindec]) + return np.array([np.exp(disc_ref(x, gamma).squeeze()) for x in sindec]) else: - return np.exp(disc_ref(sindec, gamma)) + return np.exp(disc_ref(sindec, gamma).squeeze()) def reference_10year_sensitivity(sindec=np.array(0.0), gamma=2.0): @@ -130,12 +133,14 @@ def reference_10year_sensitivity(sindec=np.array(0.0), gamma=2.0): scaling = np.array([10 ** (3 * (i)) for i in range(2)]) sens *= scaling - sens_ref = interp2d(np.array(sindecs), np.array(gammas), np.log(sens.T)) + sens_ref = RectBivariateSpline( + np.array(sindecs), np.array(gammas), np.log(sens), kx=1, ky=1 + ) if np.array(sindec).ndim > 0: - return np.array([np.exp(sens_ref(x, gamma))[0] for x in sindec]) + return np.array([np.exp(sens_ref(x, gamma).squeeze()) for x in sindec]) else: - return np.exp(sens_ref(sindec, gamma)) + return np.exp(sens_ref(sindec, gamma).squeeze()) def reference_10year_discovery_potential(sindec=np.array(0.0), gamma=2.0): @@ -162,9 +167,11 @@ def reference_10year_discovery_potential(sindec=np.array(0.0), gamma=2.0): scaling = np.array([10 ** (3 * i) for i in range(2)]) sens *= scaling - sens_ref = interp2d(np.array(sindecs), np.array(gammas), np.log(sens.T)) + sens_ref = RectBivariateSpline( + np.array(sindecs), np.array(gammas), np.log(sens), kx=1, ky=1 + ) if np.array(sindec).ndim > 0: - return np.array([np.exp(sens_ref(x, gamma))[0] for x in sindec]) + return np.array([np.exp(sens_ref(x, gamma).squeeze()) for x in sindec]) else: - return np.exp(sens_ref(sindec, gamma)) + return np.exp(sens_ref(sindec, gamma).squeeze()) diff --git a/flarestack/utils/percentile_SoB.py b/flarestack/utils/percentile_SoB.py index 34f3ecc6..6caca1a4 100644 --- a/flarestack/utils/percentile_SoB.py +++ b/flarestack/utils/percentile_SoB.py @@ -244,7 +244,9 @@ def make_plot(hist, savepath, normed=True): order = 1 - spline = scipy.interpolate.interp2d(x, y, np.log(ratio)) + spline = scipy.interpolate.RectBivariateSpline( + x, y, np.log(ratio.T), kx=order, ky=order + ) for x_val in [2.0, 3.0, 7.0]: print(x_val, spline(0.0, x_val), spline(0.5, x_val))