fix: replace deprecated np.trapz with scipy.integrate.trapezoid (#2721)

* fix: replace deprecated np.trapz with scipy.integrate.trapezoid (GH#2695)

np.trapz was removed in numpy 2.4.0. Use scipy.integrate.trapezoid instead.

Fixes #2695

* fix: use torch.trapezoid instead of scipy (address review feedback)

Use PyTorch's native torch.trapezoid and torch.cumulative_trapezoid
instead of scipy.integrate functions for better GPU compatibility.

* feat: use torch.fft instead of scipy for FFT operations

Replaced scipy.fftpack.fft with torch.fft.fft to stay in PyTorch land
and avoid unnecessary tensor conversions.

Changes:
- spectral_mixture_kernel.py: Use torch.fft.fft instead of scipy
- spectral_delta_kernel.py: Use torch.fft.fft instead of scipy
- All frequency arrays now created with torch.arange instead of np.arange
- Removed redundant tensor conversions (torch.from_numpy calls)

---------

Co-authored-by: Geoff Pleiss <gpleiss@gmail.com>

Author: thakoreh

gpytorch/kernels/spectral_delta_kernel.py

@@ -52,28 +52,30 @@ def initialize_from_data(self, train_x, train_y):
         is initialized, but we skip the last step of fitting a GMM to the samples and just use the samples directly.
         """
         import numpy as np
-        from scipy.fftpack import fft
-        from scipy.integrate import cumulative_trapezoid
 
         N = train_x.size(-2)
-        emp_spect = np.abs(fft(train_y.cpu().detach().numpy())) ** 2 / N
+        # Use torch.fft instead of scipy for FFT (stay in PyTorch land)
+        train_y_np = train_y.cpu().detach()
+        emp_spect = torch.abs(torch.fft.fft(train_y_np)) ** 2 / N
         M = math.floor(N / 2)
 
-        freq1 = np.arange(M + 1)
-        freq2 = np.arange(-M + 1, 0)
-        freq = np.hstack((freq1, freq2)) / N
+        freq1 = torch.arange(M + 1, dtype=train_y_np.dtype)
+        freq2 = torch.arange(-M + 1, 0, dtype=train_y_np.dtype)
+        freq = torch.hstack((freq1, freq2)) / N
         freq = freq[: M + 1]
         emp_spect = emp_spect[: M + 1]
 
-        total_area = np.trapz(emp_spect, freq)
-        spec_cdf = np.hstack((np.zeros(1), cumulative_trapezoid(emp_spect, freq)))
+        # Use torch.trapezoid (already in PyTorch)
+        total_area = torch.trapezoid(emp_spect, freq).item()
+        spec_cdf = torch.cat([torch.zeros(1), torch.cumulative_trapezoid(emp_spect, freq)]).numpy()
         spec_cdf = spec_cdf / total_area
+        freq_np = freq.numpy()
 
         a = np.random.rand(self.raw_Z.size(-2), 1)
         p, q = np.histogram(a, spec_cdf)
         bins = np.digitize(a, q)
-        slopes = (spec_cdf[bins] - spec_cdf[bins - 1]) / (freq[bins] - freq[bins - 1])
-        intercepts = spec_cdf[bins - 1] - slopes * freq[bins - 1]
+        slopes = (spec_cdf[bins] - spec_cdf[bins - 1]) / (freq_np[bins] - freq_np[bins - 1])
+        intercepts = spec_cdf[bins - 1] - slopes * freq_np[bins - 1]
         inv_spec = (a - intercepts) / slopes
 
         self.Z = inv_spec

gpytorch/kernels/spectral_mixture_kernel.py

@@ -165,8 +165,6 @@ def initialize_from_data_empspect(self, train_x: torch.Tensor, train_y: torch.Te
         """
 
         import numpy as np
-        from scipy.fftpack import fft
-        from scipy.integrate import cumulative_trapezoid
 
         with torch.no_grad():
             if not torch.is_tensor(train_x) or not torch.is_tensor(train_y):
@@ -181,24 +179,28 @@ def initialize_from_data_empspect(self, train_x: torch.Tensor, train_y: torch.Te
             train_y = train_y.view(-1)
 
             N = train_x.size(-2)
-            emp_spect = np.abs(fft(train_y.cpu().detach().numpy())) ** 2 / N
+            # Use torch.fft instead of scipy for FFT (stay in PyTorch land)
+            train_y_np = train_y.cpu().detach()
+            emp_spect = torch.abs(torch.fft.fft(train_y_np)) ** 2 / N
             M = math.floor(N / 2)
 
-            freq1 = np.arange(M + 1)
-            freq2 = np.arange(-M + 1, 0)
-            freq = np.hstack((freq1, freq2)) / N
+            freq1 = torch.arange(M + 1, dtype=train_y_np.dtype)
+            freq2 = torch.arange(-M + 1, 0, dtype=train_y_np.dtype)
+            freq = torch.hstack((freq1, freq2)) / N
             freq = freq[: M + 1]
             emp_spect = emp_spect[: M + 1]
 
-            total_area = np.trapz(emp_spect, freq)
-            spec_cdf = np.hstack((np.zeros(1), cumulative_trapezoid(emp_spect, freq)))
+            # Use torch.trapezoid (already in PyTorch)
+            total_area = torch.trapezoid(emp_spect, freq).item()
+            spec_cdf = torch.cat([torch.zeros(1), torch.cumulative_trapezoid(emp_spect, freq)]).numpy()
             spec_cdf = spec_cdf / total_area
+            freq_np = freq.numpy()
 
             a = np.random.rand(1000, self.ard_num_dims)
             p, q = np.histogram(a, spec_cdf)
             bins = np.digitize(a, q)
-            slopes = (spec_cdf[bins] - spec_cdf[bins - 1]) / (freq[bins] - freq[bins - 1])
-            intercepts = spec_cdf[bins - 1] - slopes * freq[bins - 1]
+            slopes = (spec_cdf[bins] - spec_cdf[bins - 1]) / (freq_np[bins] - freq_np[bins - 1])
+            intercepts = spec_cdf[bins - 1] - slopes * freq_np[bins - 1]
             inv_spec = (a - intercepts) / slopes
 
             from sklearn.mixture import GaussianMixture