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

Author: jvansanten

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)
 

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())

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))