From daf56ae246be82c6d2b616dc3ce566d6a10a9c66 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Thu, 6 Jun 2024 15:42:33 -0400
Subject: [PATCH 01/21] Add epsilon arg to PowerFilter.

---
 src/aspire/operators/filters.py | 16 ++++++++++++++--
 1 file changed, 14 insertions(+), 2 deletions(-)

diff --git a/src/aspire/operators/filters.py b/src/aspire/operators/filters.py
index 9b910a8fe0..bb7491c780 100644
--- a/src/aspire/operators/filters.py
+++ b/src/aspire/operators/filters.py
@@ -184,9 +184,19 @@ class PowerFilter(Filter):
     A Filter object that is composed of a regular `Filter` object, but evaluates it to a specified power.
     """
 
-    def __init__(self, filter, power=1):
+    def __init__(self, filter, power=1, epsilon=None):
+        """
+        Initialize PowerFilter instance.
+
+        :param filter: A Filter instance.
+        :param power: Exponent to raise filter values.
+        :param epsilon: Threshold on filter values that get raised to a negative power.
+            `filter` values below this threshold will be set to zero during evaluation.
+            Default uses machine epsilon for filter.dtype.
+        """
         self._filter = filter
         self._power = power
+        self._epsilon = epsilon
         super().__init__(dim=filter.dim, radial=filter.radial)
 
     def _evaluate(self, omega):
@@ -204,7 +214,9 @@ def evaluate_grid(self, L, *args, dtype=np.float32, **kwargs):
 
         # Place safeguard on values below machine epsilon for negative powers.
         if self._power < 0:
-            eps = np.finfo(filter_vals.dtype).eps
+            eps = self._epsilon
+            if eps is None:
+                eps = np.finfo(filter_vals.dtype).eps
             condition = abs(filter_vals) < eps
             num_less_eps = np.count_nonzero(condition)
             if num_less_eps > 0:

From 391005d211af49303b62b7ee6fad505e31991c6e Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Thu, 6 Jun 2024 16:02:42 -0400
Subject: [PATCH 02/21] Add threshold to whiten function with matlab default.

---
 src/aspire/source/image.py | 14 ++++++++++++--
 1 file changed, 12 insertions(+), 2 deletions(-)

diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index 473585acb9..f93a40440b 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -798,7 +798,7 @@ def downsample(self, L):
         self.L = L
 
     @_as_copy
-    def whiten(self, noise_estimate=None):
+    def whiten(self, noise_estimate=None, epsilon=None):
         """
         Modify the `ImageSource` in-place by appending a whitening filter to the generation pipeline.
 
@@ -810,6 +810,9 @@ def whiten(self, noise_estimate=None):
             passed a `NoiseEstimator` the `filter` attribute will be
             queried.  Alternatively, the noise PSD may be passed
             directly as a `Filter` object.
+        :param epsilon: Threshold used to determine which frequencies to whiten
+            and which to set to zero. By default all filter values less than
+            100*eps(self.dtype) are zeroed out.
         :return: On return, the `ImageSource` object has been modified in place.
         """
 
@@ -827,8 +830,15 @@ def whiten(self, noise_estimate=None):
                 " instead of `NoiseEstimator` or `Filter`."
             )
 
+        # Set threshold for whiten_filter. All values such that sqrt(noise_filter) < eps
+        # will be set to zero in the whiten_filter.
+        if epsilon is None:
+            epsilon = 100 * np.finfo(self.dtype).eps
+
         logger.info("Whitening source object")
-        whiten_filter = PowerFilter(noise_filter, power=-0.5)
+        # epsilon is squared to account for the PowerFilter applying the threshold
+        # to noise_filter, not sqrt(noise_filter).
+        whiten_filter = PowerFilter(noise_filter, power=-0.5, epsilon=epsilon**2)
 
         logger.info("Transforming all CTF Filters into Multiplicative Filters")
         self.unique_filters = [

From 7990997842a4be07fc7134392499adc99b8b0205 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Mon, 10 Jun 2024 15:27:08 -0400
Subject: [PATCH 03/21] Recast ArrayFilter result after scipy workaround upcast
 occurs.

---
 src/aspire/operators/filters.py |  3 ++-
 tests/test_filters.py           | 16 ++++++++++++++++
 2 files changed, 18 insertions(+), 1 deletion(-)

diff --git a/src/aspire/operators/filters.py b/src/aspire/operators/filters.py
index bb7491c780..df2d33b036 100644
--- a/src/aspire/operators/filters.py
+++ b/src/aspire/operators/filters.py
@@ -349,7 +349,8 @@ def _evaluate(self, omega):
         # Result is 1 x np.prod(self.sz) in shape; convert to a 1-d vector
         result = np.squeeze(result, 0)
 
-        return result
+        # Recast result with correct dtype
+        return result.astype(self.xfer_fn_array.dtype)
 
     def evaluate_grid(self, L, *args, dtype=np.float32, **kwargs):
         """
diff --git a/tests/test_filters.py b/tests/test_filters.py
index 35d7955a9e..33f8dc7349 100644
--- a/tests/test_filters.py
+++ b/tests/test_filters.py
@@ -361,3 +361,19 @@ def test_power_filter_safeguard(dtype, caplog):
     # Check caplog for warning.
     msg = f"setting {num_eps} extremal filter value(s) to zero."
     assert msg in caplog.text
+
+
+@pytest.mark.parametrize("dtype", [np.float32, np.float64])
+def test_array_filter_dtype_passthrough(dtype):
+    """
+    Do to a bug in scipy versions < 1.10.1, scipy's interpolator crashes
+    in singles. We have a workaround that upcasts to doubles. This test
+    ensures that we recast to the correct dtype during calculations.
+    """
+    L = 8
+    arr = np.ones((L, L), dtype=dtype)
+
+    filt = ArrayFilter(arr)
+    filt_vals = filt.evaluate_grid(L, dtype=dtype)
+
+    assert filt_vals.dtype == dtype

From d028e1fff5ea405c4e7784fc2718ef5623ebe170 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Mon, 10 Jun 2024 16:03:53 -0400
Subject: [PATCH 04/21] Revert to original threshold of eps(dtype) on PSD.

---
 src/aspire/source/image.py | 11 +++--------
 1 file changed, 3 insertions(+), 8 deletions(-)

diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index f93a40440b..b23b5c98e4 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -811,8 +811,8 @@ def whiten(self, noise_estimate=None, epsilon=None):
             queried.  Alternatively, the noise PSD may be passed
             directly as a `Filter` object.
         :param epsilon: Threshold used to determine which frequencies to whiten
-            and which to set to zero. By default all filter values less than
-            100*eps(self.dtype) are zeroed out.
+            and which to set to zero. By default all PSD values in the `noise_estimate`
+            less than eps(self.dtype) are zeroed out in the whitening filter.
         :return: On return, the `ImageSource` object has been modified in place.
         """
 
@@ -830,15 +830,10 @@ def whiten(self, noise_estimate=None, epsilon=None):
                 " instead of `NoiseEstimator` or `Filter`."
             )
 
-        # Set threshold for whiten_filter. All values such that sqrt(noise_filter) < eps
-        # will be set to zero in the whiten_filter.
-        if epsilon is None:
-            epsilon = 100 * np.finfo(self.dtype).eps
-
         logger.info("Whitening source object")
         # epsilon is squared to account for the PowerFilter applying the threshold
         # to noise_filter, not sqrt(noise_filter).
-        whiten_filter = PowerFilter(noise_filter, power=-0.5, epsilon=epsilon**2)
+        whiten_filter = PowerFilter(noise_filter, power=-0.5, epsilon=epsilon)
 
         logger.info("Transforming all CTF Filters into Multiplicative Filters")
         self.unique_filters = [

From c6810ffc63daa553a1ba6f1aad82ca1d85f57406 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Tue, 11 Jun 2024 15:33:23 -0400
Subject: [PATCH 05/21] remove comment

---
 src/aspire/source/image.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index b23b5c98e4..c60b64d701 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -831,8 +831,6 @@ def whiten(self, noise_estimate=None, epsilon=None):
             )
 
         logger.info("Whitening source object")
-        # epsilon is squared to account for the PowerFilter applying the threshold
-        # to noise_filter, not sqrt(noise_filter).
         whiten_filter = PowerFilter(noise_filter, power=-0.5, epsilon=epsilon)
 
         logger.info("Transforming all CTF Filters into Multiplicative Filters")

From bf9abbc12151f95974eb0da68df82d4c6252b1fd Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Tue, 11 Jun 2024 15:44:35 -0400
Subject: [PATCH 06/21] set default epsilon inside whiten function.

---
 src/aspire/source/image.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index c60b64d701..fa5be1f7f7 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -830,6 +830,9 @@ def whiten(self, noise_estimate=None, epsilon=None):
                 " instead of `NoiseEstimator` or `Filter`."
             )
 
+        if epsilon is None:
+            epsilon = np.finfo(self.dtype).eps
+
         logger.info("Whitening source object")
         whiten_filter = PowerFilter(noise_filter, power=-0.5, epsilon=epsilon)
 

From 700c0f211502b495759d48eee74ea4dd02f501b4 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Tue, 18 Jun 2024 09:03:49 -0400
Subject: [PATCH 07/21] test PowerFilter argument

---
 tests/test_filters.py | 13 ++++++++++---
 1 file changed, 10 insertions(+), 3 deletions(-)

diff --git a/tests/test_filters.py b/tests/test_filters.py
index 33f8dc7349..8d9cac9e29 100644
--- a/tests/test_filters.py
+++ b/tests/test_filters.py
@@ -1,3 +1,4 @@
+import itertools
 import logging
 import os.path
 from unittest import TestCase
@@ -332,20 +333,26 @@ def testFilterSigns(self):
         self.assertTrue(np.allclose(sign_filter.evaluate(self.omega), signs))
 
 
-@pytest.mark.parametrize("dtype", [np.float32, np.float64])
-def test_power_filter_safeguard(dtype, caplog):
+params = list(itertools.product([np.float32, np.float64], [None, 0.01]))
+
+
+@pytest.mark.parametrize("dtype, epsilon", params)
+def test_power_filter_safeguard(dtype, epsilon, caplog):
     L = 25
     arr = np.ones((L, L), dtype=dtype)
 
     # Set a few values below machine epsilon.
     num_eps = 3
-    eps = np.finfo(dtype).eps
+    eps = epsilon
+    if eps is None:
+        eps = np.finfo(dtype).eps
     arr[L // 2, L // 2 : L // 2 + num_eps] = eps / 2
 
     # For negative powers, values below machine eps will be set to zero.
     filt = PowerFilter(
         filter=ArrayFilter(arr),
         power=-0.5,
+        epsilon=epsilon,
     )
 
     caplog.clear()

From 340130d0161b64e55964f6dea67e38be8ff28711 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Tue, 18 Jun 2024 09:18:42 -0400
Subject: [PATCH 08/21] smoke test for whiten epsilon param.

---
 tests/test_preprocess_pipeline.py | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/tests/test_preprocess_pipeline.py b/tests/test_preprocess_pipeline.py
index fb7d2427ec..a37b19210a 100644
--- a/tests/test_preprocess_pipeline.py
+++ b/tests/test_preprocess_pipeline.py
@@ -126,6 +126,15 @@ def testWhiten2(dtype):
     assert np.allclose(np.eye(2), corr_coef, atol=2e-1)
 
 
+@pytest.mark.parametrize("dtype", [np.float32, np.float64])
+def test_whiten_epsilon(dtype):
+    """Smoke test for epsilon argument"""
+    L = 25
+    sim = get_sim_object(L, dtype)
+    noise_estimator = AnisotropicNoiseEstimator(sim)
+    _ = sim.whiten(noise_estimator.filter, epsilon=0.01)
+
+
 @pytest.mark.parametrize("L, dtype", params)
 def testInvertContrast(L, dtype):
     sim1 = get_sim_object(L, dtype)

From 75dfe87d63deda25154c57b73bed009c171367f2 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Tue, 18 Jun 2024 15:38:45 -0400
Subject: [PATCH 09/21] bump scipy version. remove upcast.

---
 pyproject.toml                  | 2 +-
 src/aspire/operators/filters.py | 7 +++----
 2 files changed, 4 insertions(+), 5 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index c9c25a9976..c674418ec1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -48,7 +48,7 @@ dependencies = [
     "pyshtools<=4.10.4",  # 4.11.7 might have a packaging bug
     "PyWavelets",
     "ray >= 2.9.2",
-    "scipy >= 1.10.0",
+    "scipy >= 1.10.1",
     "scikit-learn",
     "scikit-image",
     "setuptools >= 0.41",
diff --git a/src/aspire/operators/filters.py b/src/aspire/operators/filters.py
index df2d33b036..c93c34f9cb 100644
--- a/src/aspire/operators/filters.py
+++ b/src/aspire/operators/filters.py
@@ -333,8 +333,7 @@ def _evaluate(self, omega):
         #  for values slightly outside the interpolation grid bounds.
         interpolator = RegularGridInterpolator(
             _input_pts,
-            # https://github.com/scipy/scipy/issues/17718
-            self.xfer_fn_array.astype(np.float64),
+            self.xfer_fn_array,
             method="linear",
             bounds_error=False,
             fill_value=None,
@@ -349,8 +348,8 @@ def _evaluate(self, omega):
         # Result is 1 x np.prod(self.sz) in shape; convert to a 1-d vector
         result = np.squeeze(result, 0)
 
-        # Recast result with correct dtype
-        return result.astype(self.xfer_fn_array.dtype)
+        # Scipy's interpolator will upcast singles. Recasting.
+        return result.astype(self.xfer_fn_array.dtype, copy=False)
 
     def evaluate_grid(self, L, *args, dtype=np.float32, **kwargs):
         """

From 8f59c5c24c89e6bf176b8bd0c4bac922d3c117d0 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Tue, 18 Jun 2024 16:18:30 -0400
Subject: [PATCH 10/21] test that whiten safeguard is actually working.

---
 tests/test_filters.py             |  3 +--
 tests/test_preprocess_pipeline.py | 26 +++++++++++++++++++++++---
 2 files changed, 24 insertions(+), 5 deletions(-)

diff --git a/tests/test_filters.py b/tests/test_filters.py
index 8d9cac9e29..8674f13486 100644
--- a/tests/test_filters.py
+++ b/tests/test_filters.py
@@ -373,8 +373,7 @@ def test_power_filter_safeguard(dtype, epsilon, caplog):
 @pytest.mark.parametrize("dtype", [np.float32, np.float64])
 def test_array_filter_dtype_passthrough(dtype):
     """
-    Do to a bug in scipy versions < 1.10.1, scipy's interpolator crashes
-    in singles. We have a workaround that upcasts to doubles. This test
+    scipy's interpolator will upcast singles. This test
     ensures that we recast to the correct dtype during calculations.
     """
     L = 8
diff --git a/tests/test_preprocess_pipeline.py b/tests/test_preprocess_pipeline.py
index a37b19210a..17a4407008 100644
--- a/tests/test_preprocess_pipeline.py
+++ b/tests/test_preprocess_pipeline.py
@@ -127,12 +127,32 @@ def testWhiten2(dtype):
 
 
 @pytest.mark.parametrize("dtype", [np.float32, np.float64])
-def test_whiten_epsilon(dtype):
-    """Smoke test for epsilon argument"""
+def test_whiten_safeguard(dtype):
+    """Test that whitening safeguard works as expected."""
     L = 25
+    epsilon = 0.02
     sim = get_sim_object(L, dtype)
     noise_estimator = AnisotropicNoiseEstimator(sim)
-    _ = sim.whiten(noise_estimator.filter, epsilon=0.01)
+    sim = sim.whiten(noise_estimator.filter, epsilon=epsilon)
+
+    # Get whitening_filter from generation pipeline.
+    whiten_filt = sim.generation_pipeline.xforms[0].filter.evaluate_grid(sim.L)
+
+    # Generate whitening_filter without safeguard directly from noise_estimator.
+    filt_vals = noise_estimator.filter.xfer_fn_array
+    whiten_filt_unsafe = filt_vals**-0.5
+
+    # Get indices where safeguard should be applied
+    # and assert that they are not empty.
+    ind = np.where(filt_vals < epsilon)
+    np.testing.assert_array_less(0, len(ind[0]))
+
+    # Check that whiten_filt and whiten_filt_unsafe agree up to safeguard indices.
+    disagree = np.where(whiten_filt != whiten_filt_unsafe)
+    np.testing.assert_array_equal(ind, disagree)
+
+    # Check that whiten_filt is zero at safeguard indices.
+    np.testing.assert_allclose(whiten_filt[ind], 0.0)
 
 
 @pytest.mark.parametrize("L, dtype", params)

From 2b4b1d9bf6eec6de5ee28ea481cb2a3a7e963fcf Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Thu, 20 Jun 2024 09:27:20 -0400
Subject: [PATCH 11/21] use np.testing in suite that failed on arm.

---
 tests/test_preprocess_pipeline.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/tests/test_preprocess_pipeline.py b/tests/test_preprocess_pipeline.py
index 17a4407008..d7416095b0 100644
--- a/tests/test_preprocess_pipeline.py
+++ b/tests/test_preprocess_pipeline.py
@@ -101,7 +101,7 @@ def testWhiten(dtype):
     corr_coef = np.corrcoef(imgs_wt[:, L - 1, L - 1], imgs_wt[:, L - 2, L - 1])
 
     # correlation matrix should be close to identity
-    assert np.allclose(np.eye(2), corr_coef, atol=1e-1)
+    np.testing.assert_allclose(np.eye(2), corr_coef, atol=1e-1)
     # dtype of returned images should be the same
     assert dtype == imgs_wt.dtype
 
@@ -123,7 +123,7 @@ def testWhiten2(dtype):
     corr_coef = np.corrcoef(imgs_wt[:, L - 1, L - 1], imgs_wt[:, L - 2, L - 1])
 
     # Correlation matrix should be close to identity
-    assert np.allclose(np.eye(2), corr_coef, atol=2e-1)
+    np.testing.assert_allclose(np.eye(2), corr_coef, atol=2e-1)
 
 
 @pytest.mark.parametrize("dtype", [np.float32, np.float64])
@@ -167,7 +167,9 @@ def testInvertContrast(L, dtype):
     imgs2_rc = sim2.images[:num_images]
 
     # all images should be the same after inverting contrast
-    assert np.allclose(imgs1_rc.asnumpy(), imgs2_rc.asnumpy())
+    np.testing.assert_allclose(
+        imgs1_rc.asnumpy(), imgs2_rc.asnumpy(), rtol=1e-05, atol=1e-08
+    )
     # dtype of returned images should be the same
     assert dtype == imgs1_rc.dtype
     assert dtype == imgs2_rc.dtype

From 5a4dec33d8d064327991e8897c2f55f45d209738 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Thu, 20 Jun 2024 09:38:47 -0400
Subject: [PATCH 12/21] use np.testing in test_FLEbasis2D.py

---
 tests/test_FLEbasis2D.py | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/tests/test_FLEbasis2D.py b/tests/test_FLEbasis2D.py
index 7d6b3f5c47..ffb1f8f7d1 100644
--- a/tests/test_FLEbasis2D.py
+++ b/tests/test_FLEbasis2D.py
@@ -142,7 +142,7 @@ def testMatchFBEvaluate(basis):
     fb_images = fb_basis.evaluate(coefs)
     fle_images = basis.evaluate(coefs)
 
-    assert np.allclose(fb_images._data, fle_images._data, atol=1e-4)
+    np.testing.assert_allclose(fb_images._data, fle_images._data, atol=1e-4)
 
 
 @pytest.mark.parametrize("basis", test_bases_match_fb, ids=show_fle_params)
@@ -159,8 +159,8 @@ def testMatchFBDenseEvaluate(basis):
     fle_images = Image(fle_out.T.reshape(-1, basis.nres, basis.nres)).asnumpy()
 
     # Matrix column reording in match_fb mode flips signs of some of the basis functions
-    assert np.allclose(np.abs(fb_images), np.abs(fle_images), atol=1e-3)
-    assert np.allclose(fb_images, fle_images, atol=1e-3)
+    np.testing.assert_allclose(np.abs(fb_images), np.abs(fle_images), atol=1e-3)
+    np.testing.assert_allclose(fb_images, fle_images, atol=1e-3)
 
 
 @pytest.mark.parametrize("basis", test_bases_match_fb, ids=show_fle_params)
@@ -177,7 +177,7 @@ def testMatchFBEvaluate_t(basis):
     fb_coefs = fb_basis.evaluate_t(images)
     fle_coefs = basis.evaluate_t(images)
 
-    assert np.allclose(fb_coefs, fle_coefs, atol=1e-4)
+    np.testing.assert_allclose(fb_coefs, fle_coefs, atol=1e-4)
 
 
 @pytest.mark.parametrize("basis", test_bases_match_fb, ids=show_fle_params)
@@ -197,7 +197,7 @@ def testMatchFBDenseEvaluate_t(basis):
     fle_coefs = basis._create_dense_matrix().T @ vec.T
 
     # Matrix column reording in match_fb mode flips signs of some of the basis coefficients
-    assert np.allclose(np.abs(fb_coefs), np.abs(fle_coefs), atol=1e-4)
+    np.testing.assert_allclose(np.abs(fb_coefs), np.abs(fle_coefs), atol=1e-4)
 
 
 def testLowPass():
@@ -265,4 +265,4 @@ def testRadialConvolution():
             convolution_fft_pad[L // 2 : L // 2 + L, L // 2 : L // 2 + L]
         )
 
-    assert np.allclose(imgs_convolved_fle, imgs_convolved_slow, atol=1e-5)
+    np.testing.assert_allclose(imgs_convolved_fle, imgs_convolved_slow, atol=1e-5)

From 240ba1332d2b7bc22c412b9a89f4e8063cad4899 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Thu, 20 Jun 2024 11:35:10 -0400
Subject: [PATCH 13/21] revert to upcasting.

---
 src/aspire/operators/filters.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/aspire/operators/filters.py b/src/aspire/operators/filters.py
index c93c34f9cb..5720867e49 100644
--- a/src/aspire/operators/filters.py
+++ b/src/aspire/operators/filters.py
@@ -333,7 +333,7 @@ def _evaluate(self, omega):
         #  for values slightly outside the interpolation grid bounds.
         interpolator = RegularGridInterpolator(
             _input_pts,
-            self.xfer_fn_array,
+            self.xfer_fn_array.astype(np.float64),
             method="linear",
             bounds_error=False,
             fill_value=None,

From 9ded5a291a38c22bf0c641b42f36f76b1733cf88 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Thu, 20 Jun 2024 12:08:49 -0400
Subject: [PATCH 14/21] update comments

---
 src/aspire/operators/filters.py | 1 +
 tests/test_filters.py           | 2 +-
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/aspire/operators/filters.py b/src/aspire/operators/filters.py
index 5720867e49..f11202d2f2 100644
--- a/src/aspire/operators/filters.py
+++ b/src/aspire/operators/filters.py
@@ -333,6 +333,7 @@ def _evaluate(self, omega):
         #  for values slightly outside the interpolation grid bounds.
         interpolator = RegularGridInterpolator(
             _input_pts,
+            # scipy requires upcasting to use cython interpolator.
             self.xfer_fn_array.astype(np.float64),
             method="linear",
             bounds_error=False,
diff --git a/tests/test_filters.py b/tests/test_filters.py
index 8674f13486..83b68d86ca 100644
--- a/tests/test_filters.py
+++ b/tests/test_filters.py
@@ -373,7 +373,7 @@ def test_power_filter_safeguard(dtype, epsilon, caplog):
 @pytest.mark.parametrize("dtype", [np.float32, np.float64])
 def test_array_filter_dtype_passthrough(dtype):
     """
-    scipy's interpolator will upcast singles. This test
+    We upcast to use scipy's fast interpolator. This test
     ensures that we recast to the correct dtype during calculations.
     """
     L = 8

From 3a1b0c97dc80f9f34f71aae103c03f1d40546dfe Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Thu, 20 Jun 2024 13:28:50 -0400
Subject: [PATCH 15/21] remove scipy bump

---
 pyproject.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index c674418ec1..c9c25a9976 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -48,7 +48,7 @@ dependencies = [
     "pyshtools<=4.10.4",  # 4.11.7 might have a packaging bug
     "PyWavelets",
     "ray >= 2.9.2",
-    "scipy >= 1.10.1",
+    "scipy >= 1.10.0",
     "scikit-learn",
     "scikit-image",
     "setuptools >= 0.41",

From f15cbfc54203a6fa3179e3e605085a3e2d8c0dce Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Wed, 10 Jul 2024 13:50:36 -0400
Subject: [PATCH 16/21] Revert scipy workaround.

---
 src/aspire/operators/filters.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/src/aspire/operators/filters.py b/src/aspire/operators/filters.py
index f11202d2f2..bb7491c780 100644
--- a/src/aspire/operators/filters.py
+++ b/src/aspire/operators/filters.py
@@ -333,7 +333,7 @@ def _evaluate(self, omega):
         #  for values slightly outside the interpolation grid bounds.
         interpolator = RegularGridInterpolator(
             _input_pts,
-            # scipy requires upcasting to use cython interpolator.
+            # https://github.com/scipy/scipy/issues/17718
             self.xfer_fn_array.astype(np.float64),
             method="linear",
             bounds_error=False,
@@ -349,8 +349,7 @@ def _evaluate(self, omega):
         # Result is 1 x np.prod(self.sz) in shape; convert to a 1-d vector
         result = np.squeeze(result, 0)
 
-        # Scipy's interpolator will upcast singles. Recasting.
-        return result.astype(self.xfer_fn_array.dtype, copy=False)
+        return result
 
     def evaluate_grid(self, L, *args, dtype=np.float32, **kwargs):
         """

From 9c6f90d8d23e60d56bd707465a2f553c25267627 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Wed, 10 Jul 2024 14:01:42 -0400
Subject: [PATCH 17/21] xfail ArrayFilter test for singles.

---
 tests/test_filters.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/tests/test_filters.py b/tests/test_filters.py
index 83b68d86ca..f68e261025 100644
--- a/tests/test_filters.py
+++ b/tests/test_filters.py
@@ -373,9 +373,12 @@ def test_power_filter_safeguard(dtype, epsilon, caplog):
 @pytest.mark.parametrize("dtype", [np.float32, np.float64])
 def test_array_filter_dtype_passthrough(dtype):
     """
-    We upcast to use scipy's fast interpolator. This test
-    ensures that we recast to the correct dtype during calculations.
+    We upcast to use scipy's fast interpolator. We do not recast
+    on exit, so this is an expected fail for singles.
     """
+    if dtype == np.float32:
+        pytest.xfail(reason="ArrayFilter currently upcasts singles.")
+
     L = 8
     arr = np.ones((L, L), dtype=dtype)
 

From c244b42495f13eef0b97e29c0d025b2fe1d10ef8 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Fri, 12 Jul 2024 09:14:23 -0400
Subject: [PATCH 18/21] make xfail strict

---
 tests/test_filters.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_filters.py b/tests/test_filters.py
index f68e261025..e243fc0a18 100644
--- a/tests/test_filters.py
+++ b/tests/test_filters.py
@@ -377,7 +377,7 @@ def test_array_filter_dtype_passthrough(dtype):
     on exit, so this is an expected fail for singles.
     """
     if dtype == np.float32:
-        pytest.xfail(reason="ArrayFilter currently upcasts singles.")
+        pytest.xfail(reason="ArrayFilter currently upcasts singles.", strict=True)
 
     L = 8
     arr = np.ones((L, L), dtype=dtype)

From e397c3155306da790f94862ed5d80c6f79e4153e Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Fri, 12 Jul 2024 10:48:51 -0400
Subject: [PATCH 19/21] remove strict param

---
 tests/test_filters.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_filters.py b/tests/test_filters.py
index e243fc0a18..f68e261025 100644
--- a/tests/test_filters.py
+++ b/tests/test_filters.py
@@ -377,7 +377,7 @@ def test_array_filter_dtype_passthrough(dtype):
     on exit, so this is an expected fail for singles.
     """
     if dtype == np.float32:
-        pytest.xfail(reason="ArrayFilter currently upcasts singles.", strict=True)
+        pytest.xfail(reason="ArrayFilter currently upcasts singles.")
 
     L = 8
     arr = np.ones((L, L), dtype=dtype)

From 622fcdeec47b5be6d5d6ee204296f7062a3bc37b Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Tue, 23 Jul 2024 08:44:42 -0400
Subject: [PATCH 20/21] Use pytest fixtures.

---
 tests/test_filters.py | 15 ++++++++++++---
 1 file changed, 12 insertions(+), 3 deletions(-)

diff --git a/tests/test_filters.py b/tests/test_filters.py
index f68e261025..2abe391071 100644
--- a/tests/test_filters.py
+++ b/tests/test_filters.py
@@ -333,10 +333,20 @@ def testFilterSigns(self):
         self.assertTrue(np.allclose(sign_filter.evaluate(self.omega), signs))
 
 
-params = list(itertools.product([np.float32, np.float64], [None, 0.01]))
+DTYPES = [np.float32, np.float64]
+EPS = [None, 0.01]
+
+
+@pytest.fixture(params=DTYPES, ids=lambda x: f"dtype={x}", scope="module")
+def dtype(request):
+    return request.param
+
+
+@pytest.fixture(params=EPS, ids=lambda x: f"epsilon={x}", scope="module")
+def epsilon(request):
+    return request.param
 
 
-@pytest.mark.parametrize("dtype, epsilon", params)
 def test_power_filter_safeguard(dtype, epsilon, caplog):
     L = 25
     arr = np.ones((L, L), dtype=dtype)
@@ -370,7 +380,6 @@ def test_power_filter_safeguard(dtype, epsilon, caplog):
     assert msg in caplog.text
 
 
-@pytest.mark.parametrize("dtype", [np.float32, np.float64])
 def test_array_filter_dtype_passthrough(dtype):
     """
     We upcast to use scipy's fast interpolator. We do not recast

From 8e02be8c5b96bf3fbec58eeeff0c5366fa97960c Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Tue, 23 Jul 2024 08:48:15 -0400
Subject: [PATCH 21/21] remove unused import

---
 tests/test_filters.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/tests/test_filters.py b/tests/test_filters.py
index 2abe391071..911e3b347b 100644
--- a/tests/test_filters.py
+++ b/tests/test_filters.py
@@ -1,4 +1,3 @@
-import itertools
 import logging
 import os.path
 from unittest import TestCase