From cc75c461e3ebddef61c7e110af5d88fdc55494d9 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Wed, 4 Dec 2024 16:53:22 +0000 Subject: [PATCH 01/13] Vectorize generate_flm function --- s2fft/utils/signal_generator.py | 34 +++++++++++++++++++-------------- 1 file changed, 20 insertions(+), 14 deletions(-) diff --git a/s2fft/utils/signal_generator.py b/s2fft/utils/signal_generator.py index 321569c1..6467fb7a 100644 --- a/s2fft/utils/signal_generator.py +++ b/s2fft/utils/signal_generator.py @@ -1,3 +1,5 @@ +from __future__ import annotations + import numpy as np import torch @@ -12,7 +14,7 @@ def generate_flm( spin: int = 0, reality: bool = False, using_torch: bool = False, -) -> np.ndarray: +) -> np.ndarray | torch.Tensor: r""" Generate a 2D set of random harmonic coefficients. @@ -38,18 +40,22 @@ def generate_flm( """ flm = np.zeros(samples.flm_shape(L), dtype=np.complex128) - for el in range(max(L_lower, abs(spin)), L): - if reality: - flm[el, 0 + L - 1] = rng.normal() - else: - flm[el, 0 + L - 1] = rng.normal() + 1j * rng.normal() - - for m in range(1, el + 1): - flm[el, m + L - 1] = rng.normal() + 1j * rng.normal() - if reality: - flm[el, -m + L - 1] = (-1) ** m * np.conj(flm[el, m + L - 1]) - else: - flm[el, -m + L - 1] = rng.normal() + 1j * rng.normal() + min_el = max(L_lower, abs(spin)) + flm[min_el:L, L - 1] = rng.standard_normal(L - min_el) + if not reality: + flm[min_el:L, L - 1] += 1j * rng.standard_normal(L - min_el) + m_indices, el_indices = np.triu_indices(n=L, k=1, m=L) + np.array([[1], [0]]) + if min_el > 0: + in_range_el = el_indices >= min_el + m_indices = m_indices[in_range_el] + el_indices = el_indices[in_range_el] + len_indices = len(m_indices) + flm[el_indices, L - 1 - m_indices] = rng.standard_normal( + len_indices + ) + 1j * rng.standard_normal(len_indices) + flm[el_indices, L - 1 + m_indices] = (-1) ** m_indices * np.conj( + flm[el_indices, L - 1 - m_indices] + ) return torch.from_numpy(flm) if using_torch else flm @@ -61,7 +67,7 @@ def generate_flmn( L_lower: int = 0, reality: bool = False, using_torch: bool = False, -) -> np.ndarray: +) -> np.ndarray | torch.Tensor: r""" Generate a 3D set of random Wigner coefficients. From 07c8c3137cafcb1befb340eea1eecfe887efe8b4 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Thu, 12 Dec 2024 15:25:52 +0000 Subject: [PATCH 02/13] Make coefficients real for m=0 and factor out functions --- s2fft/utils/signal_generator.py | 90 +++++++++++++++++++++++++++------ 1 file changed, 75 insertions(+), 15 deletions(-) diff --git a/s2fft/utils/signal_generator.py b/s2fft/utils/signal_generator.py index 6467fb7a..889ac371 100644 --- a/s2fft/utils/signal_generator.py +++ b/s2fft/utils/signal_generator.py @@ -7,6 +7,66 @@ from s2fft.sampling import so3_samples as wigner_samples +def complex_normal( + rng: np.random.Generator, + size: int | tuple[int], + var: float, +) -> np.ndarray: + """ + Generate array of samples from zero-mean complex normal distribution. + + For `z ~ ComplexNormal(0, var)` we have that `imag(z) ~ Normal(0, var/2)` and + `real(z) ~ Normal(0, var/2)` where `Normal(μ, σ²)` is the (real-valued) normal + distribution with mean parameter `μ` and variance parameter `σ²`. + + Args: + rng: Numpy random generator object to generate samples using. + size: Output shape of array to generate. + var: Variance of complex normal distribution to generate samples from. + + Returns: + Complex-valued array of shape `size` contained generated samples. + + """ + return (rng.standard_normal(size) + 1j * rng.standard_normal(size)) * ( + var / 2 + ) ** 0.5 + + +def complex_el_and_m_indices(L: int, min_el: int) -> tuple[np.ndarray, np.ndarray]: + """ + Generate pairs of el, m indices for accessing complex harmonic coefficients. + + Equivalent to nested list-comprehension based implementation + + ``` + el_indices, m_indices = np.array( + [(el, m) for el in range(min_el, L) for m in range(1, el + 1))] + ) + ``` + + For `L, min_el = 1024, 0`, this implementation is around 80x quicker in + benchmarks compared to list-comprehension implementation. + + Args: + L: Harmonic band-limit. + min_el: Inclusive lower-bound for el indices. + + Returns: + Tuple `(el_indices, m_indices)` with both entries 1D integer-valued NumPy arrays + of same size, with values of corresponding entries corresponding to pairs of + el and m indices. + + """ + el_indices, m_indices = np.tril_indices(m=L, k=-1, n=L) + m_indices += 1 + if min_el > 0: + in_range_el = el_indices >= min_el + el_indices = el_indices[in_range_el] + m_indices = m_indices[in_range_el] + return el_indices, m_indices + + def generate_flm( rng: np.random.Generator, L: int, @@ -39,24 +99,24 @@ def generate_flm( """ flm = np.zeros(samples.flm_shape(L), dtype=np.complex128) - min_el = max(L_lower, abs(spin)) + # m = 0 coefficients are always real flm[min_el:L, L - 1] = rng.standard_normal(L - min_el) - if not reality: - flm[min_el:L, L - 1] += 1j * rng.standard_normal(L - min_el) - m_indices, el_indices = np.triu_indices(n=L, k=1, m=L) + np.array([[1], [0]]) - if min_el > 0: - in_range_el = el_indices >= min_el - m_indices = m_indices[in_range_el] - el_indices = el_indices[in_range_el] + # Construct arrays of m and el indices for entries in flm corresponding to complex- + # valued coefficients (m > 0) + el_indices, m_indices = complex_el_and_m_indices(L, min_el) len_indices = len(m_indices) - flm[el_indices, L - 1 - m_indices] = rng.standard_normal( - len_indices - ) + 1j * rng.standard_normal(len_indices) - flm[el_indices, L - 1 + m_indices] = (-1) ** m_indices * np.conj( - flm[el_indices, L - 1 - m_indices] - ) - + # Generate independent complex coefficients for positive m + flm[el_indices, L - 1 + m_indices] = complex_normal(rng, len_indices, var=2) + if reality: + # Real-valued signal so set complex coefficients for negative m using conjugate + # symmetry such that flm[el, L - 1 - m] = (-1)**m * flm[el, L - 1 + m].conj + flm[el_indices, L - 1 - m_indices] = (-1) ** m_indices * ( + flm[el_indices, L - 1 + m_indices].conj() + ) + else: + # Non-real signal so generate independent complex coefficients for negative m + flm[el_indices, L - 1 - m_indices] = complex_normal(rng, len_indices, var=2) return torch.from_numpy(flm) if using_torch else flm From 6acb1407ca18ccf69943762e74431b71d17597e4 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Thu, 12 Dec 2024 18:04:48 +0000 Subject: [PATCH 03/13] Add initial signal generator tests --- tests/test_signal_generator.py | 48 ++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 tests/test_signal_generator.py diff --git a/tests/test_signal_generator.py b/tests/test_signal_generator.py new file mode 100644 index 00000000..dcdaa921 --- /dev/null +++ b/tests/test_signal_generator.py @@ -0,0 +1,48 @@ +import numpy as np +import pytest + +import s2fft +import s2fft.sampling as smp +import s2fft.utils.signal_generator as gen + +L_values_to_test = [4, 7, 64] +L_lower_to_test = [0] +spin_to_test = [-2, 0, 1] +reality_values_to_test = [False, True] + + +@pytest.mark.parametrize("L", L_values_to_test) +@pytest.mark.parametrize("min_el", [0, 1]) +def test_complex_el_and_m_indices(L, min_el): + expected_el_indices, expected_m_indices = np.array( + [(el, m) for el in range(min_el, L) for m in range(1, el + 1)] + ).T + el_indices, m_indices = gen.complex_el_and_m_indices(L, min_el) + assert (el_indices == expected_el_indices).all() + assert (m_indices == expected_m_indices).all() + + +def check_flm_conjugate_symmetry(flm, L, min_el): + for el in range(min_el, L): + for m in range(el + 1): + assert flm[el, L - 1 - m] == (-1) ** m * flm[el, L - 1 + m].conj() + + +@pytest.mark.parametrize("L", L_values_to_test) +@pytest.mark.parametrize("L_lower", L_lower_to_test) +@pytest.mark.parametrize("spin", spin_to_test) +@pytest.mark.parametrize("reality", reality_values_to_test) +@pytest.mark.filterwarnings("ignore::RuntimeWarning") +def test_generate_flm(rng, L, L_lower, spin, reality): + if reality and spin != 0: + pytest.skip("Reality only valid for scalar fields (spin=0).") + flm = gen.generate_flm(rng, L, L_lower, spin, reality) + assert flm.shape == smp.s2_samples.flm_shape(L) + assert flm.dtype == np.complex128 + assert np.isfinite(flm).all() + if reality: + check_flm_conjugate_symmetry(flm, L, max(L_lower, abs(spin))) + f_complex = s2fft.inverse(flm, L, spin=spin, reality=False, L_lower=L_lower) + assert np.allclose(f_complex.imag, 0) + f_real = s2fft.inverse(flm, L, spin=spin, reality=True, L_lower=L_lower) + assert np.allclose(f_complex.real, f_real) From 7b4c779d37dcb9d301dc1aa9ebef0a8adfcdb371 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Thu, 12 Dec 2024 18:05:26 +0000 Subject: [PATCH 04/13] Create copies of args to avoid in-place updates --- s2fft/transforms/spherical.py | 3 +++ s2fft/transforms/wigner.py | 3 +++ 2 files changed, 6 insertions(+) diff --git a/s2fft/transforms/spherical.py b/s2fft/transforms/spherical.py index 1eb138ac..036a7402 100644 --- a/s2fft/transforms/spherical.py +++ b/s2fft/transforms/spherical.py @@ -155,6 +155,9 @@ def inverse_numpy( m_start_ind = L - 1 if reality else 0 L0 = L_lower + # Copy flm argument to avoid in-place updates being propagated back to caller + flm = flm.copy() + # Apply harmonic normalisation flm[L0:] = np.einsum( "lm,l->lm", flm[L0:], np.sqrt((2 * np.arange(L0, L) + 1) / (4 * np.pi)) diff --git a/s2fft/transforms/wigner.py b/s2fft/transforms/wigner.py index 92090cf8..3e5ad72a 100644 --- a/s2fft/transforms/wigner.py +++ b/s2fft/transforms/wigner.py @@ -156,6 +156,9 @@ def inverse_numpy( ) fban = np.zeros(samples.f_shape(L, N, sampling, nside), dtype=np.complex128) + # Copy flm argument to avoid in-place updates being propagated back to caller + flmn = flmn.copy() + flmn[:, L_lower:] = np.einsum( "...nlm,...l->...nlm", flmn[:, L_lower:], From bd76a11e0795e9594e887de1cdc75c88a383e0d3 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Thu, 12 Dec 2024 18:06:33 +0000 Subject: [PATCH 05/13] Correct code snippet in docstring --- s2fft/utils/signal_generator.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/s2fft/utils/signal_generator.py b/s2fft/utils/signal_generator.py index 889ac371..21a164f5 100644 --- a/s2fft/utils/signal_generator.py +++ b/s2fft/utils/signal_generator.py @@ -41,8 +41,8 @@ def complex_el_and_m_indices(L: int, min_el: int) -> tuple[np.ndarray, np.ndarra ``` el_indices, m_indices = np.array( - [(el, m) for el in range(min_el, L) for m in range(1, el + 1))] - ) + [(el, m) for el in range(min_el, L) for m in range(1, el + 1)] + ).T ``` For `L, min_el = 1024, 0`, this implementation is around 80x quicker in From 6e9a93041fa96109d270bdbba2546e52b0f8dbf6 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Thu, 12 Dec 2024 18:06:56 +0000 Subject: [PATCH 06/13] Add vectorized flmn generator implementation --- s2fft/utils/signal_generator.py | 62 ++++++++++++++++++++++----------- 1 file changed, 41 insertions(+), 21 deletions(-) diff --git a/s2fft/utils/signal_generator.py b/s2fft/utils/signal_generator.py index 21a164f5..c6659dbd 100644 --- a/s2fft/utils/signal_generator.py +++ b/s2fft/utils/signal_generator.py @@ -153,26 +153,46 @@ def generate_flmn( """ flmn = np.zeros(wigner_samples.flmn_shape(L, N), dtype=np.complex128) - for n in range(-N + 1, N): - for el in range(max(L_lower, abs(n)), L): - if reality: - flmn[N - 1 + n, el, 0 + L - 1] = rng.normal() - flmn[N - 1 - n, el, 0 + L - 1] = (-1) ** n * flmn[ - N - 1 + n, - el, - 0 + L - 1, - ] - else: - flmn[N - 1 + n, el, 0 + L - 1] = rng.normal() + 1j * rng.normal() - - for m in range(1, el + 1): - flmn[N - 1 + n, el, m + L - 1] = rng.normal() + 1j * rng.normal() - if reality: - flmn[N - 1 - n, el, -m + L - 1] = (-1) ** (m + n) * np.conj( - flmn[N - 1 + n, el, m + L - 1] - ) - else: - flmn[N - 1 + n, el, -m + L - 1] = rng.normal() + 1j * rng.normal() - + min_el = max(L_lower, abs(n)) + # Separately deal with m = 0 case + if reality: + if n == 0: + # For m = n = 0 + # flmn[N - 1, el, L - 1] = flmn[N - 1, el, L - 1].conj (real-valued) + # Generate independent real coefficients for n = 0 + flmn[N - 1, min_el:L, L - 1] = rng.standard_normal(L - min_el) + elif n > 0: + # Generate independent complex coefficients for positive n + flmn[N - 1 + n, min_el:L, L - 1] = complex_normal( + rng, L - min_el, var=2 + ) + # For m = 0, n > 0 + # flmn[N - 1 - n, el, L - 1] = (-1)**n * flmn[N - 1 + n, el, L - 1].conj + flmn[N - 1 - n, min_el:L, L - 1] = (-1) ** n * ( + flmn[N - 1 + n, min_el:L, L - 1].conj() + ) + else: + flmn[N - 1 + n, min_el:L, L - 1] = complex_normal(rng, L - min_el, var=2) + # Construct arrays of m and el indices for entries in flmn slices for n + # corresponding to complex-valued coefficients (m > 0) + el_indices, m_indices = complex_el_and_m_indices(L, min_el) + len_indices = len(m_indices) + # Generate independent complex coefficients for positive m + flmn[N - 1 + n, el_indices, L - 1 + m_indices] = complex_normal( + rng, len_indices, var=2 + ) + if reality: + # Real-valued signal so set complex coefficients for negative m using + # conjugate symmetry relationship + # flmn[N - 1 - n, el, L - 1 - m] = + # (-1)**(m + n) * flmn[N - 1 + n, el, L - 1 + m].conj + flmn[N - 1 - n, el_indices, L - 1 - m_indices] = (-1) ** (m_indices + n) * ( + flmn[N - 1 + n, el_indices, L - 1 + m_indices].conj() + ) + else: + # Complex signal so generate independent complex coefficients for negative m + flmn[N - 1 + n, el_indices, L - 1 - m_indices] = complex_normal( + rng, len_indices, var=2 + ) return torch.from_numpy(flmn) if using_torch else flmn From 019da14fc0f5d758ef44ed7debeebb1b4524f1eb Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Fri, 13 Dec 2024 11:42:01 +0000 Subject: [PATCH 07/13] Fix incorrect passing through of L_lower to generate_precomputes --- s2fft/transforms/otf_recursions.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/s2fft/transforms/otf_recursions.py b/s2fft/transforms/otf_recursions.py index 50718eb2..bb036e46 100644 --- a/s2fft/transforms/otf_recursions.py +++ b/s2fft/transforms/otf_recursions.py @@ -80,7 +80,14 @@ def inverse_latitudinal_step( half_slices = [el + mm + 1, el - mm + 1] if precomps is None: - precomps = generate_precomputes(L, -mm, sampling, nside, L_lower) + precomps = generate_precomputes( + L=L, + spin=-mm, + sampling=sampling, + nside=nside, + forward=False, + L_lower=L_lower, + ) lrenorm, vsign, cpi, cp2, indices = precomps for i in range(2): From 62633047a8a5aaf4703b98595cc089cc5a85cb67 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Fri, 13 Dec 2024 11:44:07 +0000 Subject: [PATCH 08/13] Add test for complex_normal function --- tests/test_signal_generator.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/tests/test_signal_generator.py b/tests/test_signal_generator.py index dcdaa921..7dced273 100644 --- a/tests/test_signal_generator.py +++ b/tests/test_signal_generator.py @@ -11,6 +11,26 @@ reality_values_to_test = [False, True] +@pytest.mark.parametrize("size", (10, 100, 1000)) +@pytest.mark.parametrize("var", (1, 2)) +def test_complex_normal(rng, size, var): + samples = gen.complex_normal(rng, size, var) + assert samples.dtype == np.complex128 + assert samples.size == size + mean = samples.mean() + # Error in real + imag components of mean estimate ~ Normal(0, (var / 2) / size) + # Therefore difference between mean estimate and true zero value should be + # less than 3 * sqrt(var / (2 * size)) with probability 0.997 + mean_error_tol = 3 * (var / (2 * size)) ** 0.5 + assert abs(mean.imag) < mean_error_tol and abs(mean.real) < mean_error_tol + # If S is (unbiased) sample variance estimate then (size - 1) * S / var is a + # chi-squared distributed random variable with (size - 1) degrees of freedom + # For size >> 1, S ~approx Normal(var, 2 * var**2 / (size - 1)) so error in + # variance estimate should be less than 3 * sqrt(2 * var**2 / (size - 1)) + # with high probability + assert abs(samples.var(ddof=1) - var) < 3 * (2 * var**2 / (size - 1)) ** 0.5 + + @pytest.mark.parametrize("L", L_values_to_test) @pytest.mark.parametrize("min_el", [0, 1]) def test_complex_el_and_m_indices(L, min_el): From d976243dade59c933de0c9f0c165a6b5c4b53976 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Fri, 13 Dec 2024 11:44:29 +0000 Subject: [PATCH 09/13] Check flm zero pattern in test --- tests/test_signal_generator.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/tests/test_signal_generator.py b/tests/test_signal_generator.py index 7dced273..4dbe166a 100644 --- a/tests/test_signal_generator.py +++ b/tests/test_signal_generator.py @@ -42,6 +42,13 @@ def test_complex_el_and_m_indices(L, min_el): assert (m_indices == expected_m_indices).all() +def check_flm_zeros(flm, L, min_el): + for el in range(L): + for m in range(L): + if el < min_el or m > el: + assert flm[el, L - 1 + m] == flm[el, L - 1 - m] == 0 + + def check_flm_conjugate_symmetry(flm, L, min_el): for el in range(min_el, L): for m in range(el + 1): @@ -60,6 +67,7 @@ def test_generate_flm(rng, L, L_lower, spin, reality): assert flm.shape == smp.s2_samples.flm_shape(L) assert flm.dtype == np.complex128 assert np.isfinite(flm).all() + check_flm_zeros(flm, L, max(L_lower, abs(spin))) if reality: check_flm_conjugate_symmetry(flm, L, max(L_lower, abs(spin))) f_complex = s2fft.inverse(flm, L, spin=spin, reality=False, L_lower=L_lower) From 3e4e8ad243ae1416507c6b40bf24c1342a126b92 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Fri, 13 Dec 2024 11:44:40 +0000 Subject: [PATCH 10/13] Add test for generate_flmn --- tests/test_signal_generator.py | 47 ++++++++++++++++++++++++++++++++-- 1 file changed, 45 insertions(+), 2 deletions(-) diff --git a/tests/test_signal_generator.py b/tests/test_signal_generator.py index 4dbe166a..4e872a67 100644 --- a/tests/test_signal_generator.py +++ b/tests/test_signal_generator.py @@ -5,8 +5,8 @@ import s2fft.sampling as smp import s2fft.utils.signal_generator as gen -L_values_to_test = [4, 7, 64] -L_lower_to_test = [0] +L_values_to_test = [6, 7, 16] +L_lower_to_test = [0, 1] spin_to_test = [-2, 0, 1] reality_values_to_test = [False, True] @@ -74,3 +74,46 @@ def test_generate_flm(rng, L, L_lower, spin, reality): assert np.allclose(f_complex.imag, 0) f_real = s2fft.inverse(flm, L, spin=spin, reality=True, L_lower=L_lower) assert np.allclose(f_complex.real, f_real) + + +def check_flmn_zeros(flmn, L, N, L_lower): + for n in range(-N + 1, N): + min_el = max(L_lower, abs(n)) + for el in range(L): + for m in range(L): + if el < min_el or m > el: + assert ( + flmn[N - 1 + n, el, L - 1 + m] + == flmn[N - 1 + n, el, L - 1 - m] + == 0 + ) + + +def check_flmn_conjugate_symmetry(flmn, L, N, L_lower): + for n in range(-N + 1, N): + min_el = max(L_lower, abs(n)) + for el in range(min_el, L): + for m in range(el + 1): + assert ( + flmn[N - 1 - n, el, L - 1 - m] + == (-1) ** (m + n) * flmn[N - 1 + n, el, L - 1 + m].conj() + ) + + +@pytest.mark.parametrize("L", L_values_to_test) +@pytest.mark.parametrize("N", [1, 2]) +@pytest.mark.parametrize("L_lower", L_lower_to_test) +@pytest.mark.parametrize("reality", reality_values_to_test) +@pytest.mark.filterwarnings("ignore::RuntimeWarning") +def test_generate_flmn(rng, L, N, L_lower, reality): + flmn = gen.generate_flmn(rng, L, N, L_lower, reality) + assert flmn.shape == smp.so3_samples.flmn_shape(L, N) + assert flmn.dtype == np.complex128 + assert np.isfinite(flmn).all() + check_flmn_zeros(flmn, L, N, L_lower) + if reality: + check_flmn_conjugate_symmetry(flmn, L, N, L_lower) + f_complex = s2fft.wigner.inverse(flmn, L, N, reality=False, L_lower=L_lower) + assert np.allclose(f_complex.imag, 0) + f_real = s2fft.wigner.inverse(flmn, L, N, reality=True, L_lower=L_lower) + assert np.allclose(f_complex.real, f_real) From b789ac9f19982d441694bf38929b04b2f4908f61 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Fri, 13 Dec 2024 12:03:58 +0000 Subject: [PATCH 11/13] Correct reference to flm rather than flmn in comment [skip ci] --- s2fft/transforms/wigner.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/s2fft/transforms/wigner.py b/s2fft/transforms/wigner.py index 3e5ad72a..7c195737 100644 --- a/s2fft/transforms/wigner.py +++ b/s2fft/transforms/wigner.py @@ -156,7 +156,7 @@ def inverse_numpy( ) fban = np.zeros(samples.f_shape(L, N, sampling, nside), dtype=np.complex128) - # Copy flm argument to avoid in-place updates being propagated back to caller + # Copy flmn argument to avoid in-place updates being propagated back to caller flmn = flmn.copy() flmn[:, L_lower:] = np.einsum( From aefa0f254570dddf3b7a0520adc31e1b245e34f2 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Fri, 13 Dec 2024 13:01:47 +0000 Subject: [PATCH 12/13] Fix issue with raising integer to negative power --- s2fft/utils/signal_generator.py | 9 ++++++--- tests/test_signal_generator.py | 2 +- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/s2fft/utils/signal_generator.py b/s2fft/utils/signal_generator.py index c6659dbd..407e9614 100644 --- a/s2fft/utils/signal_generator.py +++ b/s2fft/utils/signal_generator.py @@ -187,9 +187,12 @@ def generate_flmn( # conjugate symmetry relationship # flmn[N - 1 - n, el, L - 1 - m] = # (-1)**(m + n) * flmn[N - 1 + n, el, L - 1 + m].conj - flmn[N - 1 - n, el_indices, L - 1 - m_indices] = (-1) ** (m_indices + n) * ( - flmn[N - 1 + n, el_indices, L - 1 + m_indices].conj() - ) + # As (m_indices + n) can be negative use floating point value (-1.0) as + # base of exponentation operation to avoid Numpy + # 'ValueError: Integers to negative integer powers are not allowed' error + flmn[N - 1 - n, el_indices, L - 1 - m_indices] = (-1.0) ** ( + m_indices + n + ) * flmn[N - 1 + n, el_indices, L - 1 + m_indices].conj() else: # Complex signal so generate independent complex coefficients for negative m flmn[N - 1 + n, el_indices, L - 1 - m_indices] = complex_normal( diff --git a/tests/test_signal_generator.py b/tests/test_signal_generator.py index 4e872a67..4f9be7d0 100644 --- a/tests/test_signal_generator.py +++ b/tests/test_signal_generator.py @@ -101,7 +101,7 @@ def check_flmn_conjugate_symmetry(flmn, L, N, L_lower): @pytest.mark.parametrize("L", L_values_to_test) -@pytest.mark.parametrize("N", [1, 2]) +@pytest.mark.parametrize("N", [1, 2, 3]) @pytest.mark.parametrize("L_lower", L_lower_to_test) @pytest.mark.parametrize("reality", reality_values_to_test) @pytest.mark.filterwarnings("ignore::RuntimeWarning") From 17055fb12c7dd3828c0e8254008ea9dad2f150ba Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Fri, 13 Dec 2024 14:13:30 +0000 Subject: [PATCH 13/13] Relax tolerance in test to remove spurious fail --- tests/test_spherical_precompute.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_spherical_precompute.py b/tests/test_spherical_precompute.py index 629c4aa9..fee69471 100644 --- a/tests/test_spherical_precompute.py +++ b/tests/test_spherical_precompute.py @@ -298,7 +298,7 @@ def test_transform_inverse_high_spin( kernel = c.spin_spherical_kernel(L, spin, reality, sampling, forward=False) f = inverse(flm, L, spin, kernel, sampling, reality, "numpy") - tol = 1e-8 if sampling.lower() in ["dh", "gl"] else 1e-12 + tol = 1e-8 if sampling.lower() in ["dh", "gl"] else 1e-11 np.testing.assert_allclose(f, f_check, atol=tol, rtol=tol)