add gl sampling (f00 quad test failing)

Author: CosmoMatt

s2fft/sampling/s2_samples.py

@@ -10,7 +10,7 @@ def ntheta(L: int = None, sampling: str = "mw", nside: int = None) -> int:
             Defaults to None.
 
         sampling (str, optional): Sampling scheme.  Supported sampling schemes include
-            {"mw", "mwss", "dh", "healpix"}.  Defaults to "mw".
+            {"mw", "mwss", "dh", "gl", "healpix"}.  Defaults to "mw".
 
         nside (int, optional): HEALPix Nside resolution parameter.  Only required
             if sampling="healpix".  Defaults to None.
@@ -31,7 +31,7 @@ def ntheta(L: int = None, sampling: str = "mw", nside: int = None) -> int:
             f"Sampling scheme sampling={sampling} with L={L} not supported"
         )
 
-    if sampling.lower() == "mw":
+    if sampling.lower() in ["mw", "gl"]:
         return L
 
     elif sampling.lower() == "mwss":
@@ -93,7 +93,7 @@ def nphi_equiang(L: int, sampling: str = "mw") -> int:
         L (int): Harmonic band-limit.
 
         sampling (str, optional): Sampling scheme.  Supported sampling schemes include
-            {"mw", "mwss", "dh"}.  Defaults to "mw".
+            {"mw", "mwss", "dh", "gl"}.  Defaults to "mw".
 
     Raises:
         ValueError: HEALPix sampling scheme.
@@ -104,7 +104,7 @@ def nphi_equiang(L: int, sampling: str = "mw") -> int:
         int: Number of :math:`\phi` samples.
     """
 
-    if sampling.lower() == "mw":
+    if sampling.lower() in ["mw", "gl"]:
         return 2 * L - 1
 
     elif sampling.lower() == "mwss":
@@ -129,7 +129,7 @@ def ftm_shape(L: int, sampling: str = "mw", nside: int = None) -> Tuple[int, int
         L (int): Harmonic band-limit.
 
         sampling (str, optional): Sampling scheme.  Supported sampling schemes include
-            {"mw", "mwss", "dh", "healpix"}.  Defaults to "mw".
+            {"mw", "mwss", "dh", "gl", "healpix"}.  Defaults to "mw".
 
         nside (int, optional): HEALPix Nside resolution parameter.
 
@@ -144,7 +144,7 @@ def ftm_shape(L: int, sampling: str = "mw", nside: int = None) -> Tuple[int, int
     if sampling.lower() in ["mwss", "healpix"]:
         return ntheta(L, sampling, nside), 2 * L
 
-    elif sampling.lower() in ["mw", "dh"]:
+    elif sampling.lower() in ["mw", "dh", "gl"]:
         return ntheta(L, sampling, nside), 2 * L - 1
 
     else:
@@ -203,14 +203,17 @@ def thetas(L: int = None, sampling: str = "mw", nside: int = None) -> np.ndarray
             Defaults to None.
 
         sampling (str, optional): Sampling scheme.  Supported sampling schemes include
-            {"mw", "mwss", "dh", "healpix"}.  Defaults to "mw".
+            {"mw", "mwss", "dh", "gl", "healpix"}.  Defaults to "mw".
 
         nside (int, optional): HEALPix Nside resolution parameter.  Only required
             if sampling="healpix".  Defaults to None.
 
     Returns:
         np.ndarray: Array of :math:`\theta` samples for given sampling scheme.
     """
+    if sampling.lower() == "gl":
+        return np.flip(np.arccos(np.polynomial.legendre.leggauss(L)[0]))
+
     t = np.arange(0, ntheta(L=L, sampling=sampling, nside=nside)).astype(np.float64)
 
     return t2theta(t, L, sampling, nside)
@@ -228,7 +231,7 @@ def t2theta(
             Defaults to None.
 
         sampling (str, optional): Sampling scheme.  Supported sampling schemes include
-            {"mw", "mwss", "dh", "healpix"}.  Defaults to "mw".
+            {"mw", "mwss", "dh", "gl", "healpix"}.  Defaults to "mw".
 
         nside (int, optional): HEALPix Nside resolution parameter.  Only required
             if sampling="healpix".  Defaults to None.
@@ -346,7 +349,7 @@ def phis_equiang(L: int, sampling: str = "mw") -> np.ndarray:
         L (int, optional): Harmonic band-limit.
 
         sampling (str, optional): Sampling scheme.  Supported equiangular sampling
-            schemes include {"mw", "mwss", "dh"}.  Defaults to "mw".
+            schemes include {"mw", "mwss", "dh", "gl"}.  Defaults to "mw".
 
     Returns:
         np.ndarray: Array of :math:`\phi` samples for given sampling scheme.
@@ -365,7 +368,7 @@ def p2phi_equiang(L: int, p: int, sampling: str = "mw") -> np.ndarray:
         p (int): :math:`\phi` index.
 
         sampling (str, optional): Sampling scheme.  Supported equiangular sampling
-            schemes include {"mw", "mwss", "dh"}.  Defaults to "mw".
+            schemes include {"mw", "mwss", "dh", "gl"}.  Defaults to "mw".
 
     Raises:
         ValueError: HEALPix sampling not support (only equiangular schemes supported).
@@ -376,7 +379,7 @@ def p2phi_equiang(L: int, p: int, sampling: str = "mw") -> np.ndarray:
         np.ndarray: :math:`\phi` sample(s) for given sampling scheme.
     """
 
-    if sampling.lower() == "mw":
+    if sampling.lower() in ["mw", "gl"]:
         return 2 * p * np.pi / (2 * L - 1)
 
     elif sampling.lower() == "mwss":
@@ -431,7 +434,7 @@ def f_shape(L: int = None, sampling: str = "mw", nside: int = None) -> Tuple[int
         L (int, optional): Harmonic band-limit.
 
         sampling (str, optional): Sampling scheme.  Supported sampling schemes include
-            {"mw", "mwss", "dh", "healpix"}.  Defaults to "mw".
+            {"mw", "mwss", "dh", "gl", "healpix"}.  Defaults to "mw".
 
         nside (int, optional): HEALPix Nside resolution parameter.  Only required
             if sampling="healpix".  Defaults to None.

s2fft/utils/quadrature.py

@@ -16,7 +16,7 @@ def quad_weights_transform(
         L (int): Harmonic band-limit.
 
         sampling (str, optional): Sampling scheme.  Supported sampling schemes include
-            {"mwss", "dh", "healpix}.  Defaults to "mwss".
+            {"mwss", "dh", "gl", "healpix}.  Defaults to "mwss".
 
         spin (int, optional): Harmonic spin. Defaults to 0.
 
@@ -38,6 +38,9 @@ def quad_weights_transform(
     elif sampling.lower() == "dh":
         return quad_weights_dh(L)
 
+    elif sampling.lower() == "gl":
+        return quad_weights_gl(L)
+
     elif sampling.lower() == "healpix":
         return quad_weights_hp(nside)
 
@@ -56,7 +59,7 @@ def quad_weights(
             Defaults to None.
 
         sampling (str, optional): Sampling scheme.  Supported sampling schemes include
-            {"mw", "mwss", "dh", "healpix"}.  Defaults to "mw".
+            {"mw", "mwss", "dh", "gl", "healpix"}.  Defaults to "mw".
 
         spin (int, optional): Harmonic spin. Defaults to 0.
 
@@ -80,6 +83,9 @@ def quad_weights(
     elif sampling.lower() == "dh":
         return quad_weights_dh(L)
 
+    elif sampling.lower() == "gl":
+        return quad_weights_gl(L)
+
     elif sampling.lower() == "healpix":
         return quad_weights_hp(nside)
 
@@ -112,6 +118,43 @@ def quad_weights_hp(nside: int) -> np.ndarray:
     return hp_weights
 
 
+def quad_weights_gl(L: int) -> np.ndarray:
+    r"""Compute GL quadrature weights for :math:`\theta` and :math:`\phi` integration.
+
+    Args:
+        L (int): Harmonic band-limit.
+
+    Returns:
+        np.ndarray: Weights computed for each :math:`\theta` (weights are identical
+        as :math:`\phi` varies for given :math:`\theta`).
+    """
+    x1, x2 = -1.0, 1.0
+    ntheta = samples.ntheta(L, "gl")
+    weights = np.zeros(ntheta, dtype=np.float64)
+
+    m = int((L + 1) / 2)
+    x1 = 0.5 * (x2 - x1)
+
+    for i in range(1, m + 1):
+        z = np.cos(np.pi * (i - 0.25) / (L + 0.5))
+        z1 = 2.0
+        while np.abs(z - z1) > 1e-14:
+            p1 = 1.0
+            p2 = 0.0
+            for j in range(1, L + 1):
+                p3 = p2
+                p2 = p1
+                p1 = ((2.0 * j - 1.0) * z * p2 - (j - 1.0) * p3) / j
+            pp = L * (z * p1 - p2) / (z * z - 1.0)
+            z1 = z
+            z = z1 - p1 / pp
+
+        weights[i - 1] = 2.0 * x1 / ((1.0 - z**2) * pp * pp)
+        weights[L + 1 - i - 1] = weights[i - 1]
+
+    return weights
+
+
 def quad_weights_dh(L: int) -> np.ndarray:
     r"""Compute DH quadrature weights for :math:`\theta` and :math:`\phi` integration.
 

tests/test_quadrature.py

@@ -6,12 +6,12 @@
 
 
 @pytest.mark.parametrize("L", [5, 6])
-@pytest.mark.parametrize("sampling", ["mw", "mwss"])
+@pytest.mark.parametrize("sampling", ["gl"])
+# @pytest.mark.parametrize("sampling", ["mw", "mwss", "dh", "gl"])
 def test_quadrature_mw_weights(flm_generator, L: int, sampling: str):
     spin = 0
 
     q = quadrature.quad_weights(L, sampling, spin)
-
     flm = flm_generator(L, spin, reality=False)
 
     f = spherical.inverse(flm, L, spin, sampling)

tests/test_samples.py

@@ -9,7 +9,7 @@
 
 
 @pytest.mark.parametrize("L", [15, 16])
-@pytest.mark.parametrize("sampling", ["mw", "mwss", "dh"])
+@pytest.mark.parametrize("sampling", ["mw", "mwss", "dh", "gl"])
 def test_samples_n_and_angles(L: int, sampling: str):
     # Test ntheta and nphi
     ntheta = samples.ntheta(L, sampling)
@@ -18,8 +18,11 @@ def test_samples_n_and_angles(L: int, sampling: str):
     assert (ntheta, nphi) == pytest.approx((ntheta_ssht, nphi_ssht))
 
     # Test thetas and phis
-    t = np.arange(0, ntheta)
-    thetas = samples.t2theta(t, L, sampling)
+    if sampling.lower() == "gl":
+        thetas = samples.thetas(L, sampling)
+    else:
+        t = np.arange(0, ntheta)
+        thetas = samples.t2theta(t, L, sampling)
     p = np.arange(0, nphi)
     phis = samples.p2phi_equiang(L, p, sampling)
     thetas_ssht, phis_ssht = ssht.sample_positions(L, sampling.upper())