feat: add Gauss-Legendre quadrature backend for exchange kernels and update documentation and dispatcher.

skilledwolf · skilledwolf · commit cd05cbdd445d · 2025-11-18T23:16:14.000-06:00
diff --git a/README.md b/README.md
@@ -7,10 +7,33 @@ related calculations into a small, reusable package. It provides:
 
 - Analytic Landau-level plane-wave form factors $F_{n',n}(\mathbf{q})$.
 - Exchange kernels $X_{n_1 m_1 n_2 m_2}(\mathbf{G})$ computed via:
+  - Gauss-Legendre quadrature with rational mapping (`gausslegendre` backend, default).
   - Generalized Gauss–Laguerre quadrature (`gausslag` backend).
   - Hankel-transform based integration (`hankel` backend).
 - Symmetry diagnostics for verifying kernel implementations on a given G-grid.
 
+## Backends and Reliability
+
+The package provides three backends for computing exchange kernels, each with different performance and stability characteristics:
+
+1.  **`gausslegendre` (Default)**:
+    -   **Method**: Gauss-Legendre quadrature mapped from $[-1, 1]$ to $[0, \infty)$ via a rational mapping.
+    -   **Pros**: Fast and numerically stable for all Landau level indices ($n$).
+    -   **Cons**: May require tuning `nquad` for extremely large momenta or indices ($n > 100$).
+    -   **Recommended for**: General usage, especially for large $n$ ($n \ge 10$).
+
+2.  **`gausslag`**:
+    -   **Method**: Generalized Gauss-Laguerre quadrature.
+    -   **Pros**: Very fast for small $n$.
+    -   **Cons**: Numerically unstable for large $n$ ($n \ge 12$) due to high-order Laguerre polynomial roots.
+    -   **Recommended for**: Small systems ($n < 10$) where speed is critical.
+
+3.  **`hankel`**:
+    -   **Method**: Discrete Hankel transform.
+    -   **Pros**: High precision and stability.
+    -   **Cons**: Significantly slower than quadrature methods.
+    -   **Recommended for**: Reference calculations and verifying other backends.
+
 ## Mathematical Definitions
 
 ### Plane-Wave Form Factors
@@ -66,7 +89,7 @@ thetas = np.array([0.0, 0.0, np.pi])
 nmax = 2
 
 F = get_form_factors(Gs_dimless, thetas, nmax)          # shape (nG, nmax, nmax)
-X = get_exchange_kernels(Gs_dimless, thetas, nmax)      # default 'gausslag' backend
+X = get_exchange_kernels(Gs_dimless, thetas, nmax)      # default 'gausslegendre' backend
 
 print("F shape:", F.shape)
 print("X shape:", X.shape)
diff --git a/src/quantumhall_matrixelements/__init__.py b/src/quantumhall_matrixelements/__init__.py
@@ -17,6 +17,7 @@
 from .planewave import get_form_factors
 from .exchange_gausslag import get_exchange_kernels_GaussLag
 from .exchange_hankel import get_exchange_kernels_hankel
+from .exchange_legendre import get_exchange_kernels_GaussLegendre
 
 if TYPE_CHECKING:
     from numpy.typing import NDArray
@@ -31,6 +32,7 @@ def get_exchange_kernels(
     nmax: int,
     *,
     method: str | None = None,
+    **kwargs,
 ) -> "ComplexArray":
     """Dispatcher for exchange kernels.
 
@@ -44,27 +46,35 @@ def get_exchange_kernels(
     method :
         Backend selector:
 
-        - ``'gausslag'`` (default): generalized Gauss–Laguerre quadrature
-          over the radial integral, using the analytic angular dependence.
-        - ``'hankel'``          : Hankel-transform based implementation.
+        - ``'gausslegendre'`` (default): Gauss-Legendre quadrature with rational mapping.
+          Recommended for all nmax.
+        - ``'gausslag'``: generalized Gauss–Laguerre quadrature.
+          Fast for small nmax (< 10), but unstable for large nmax.
+        - ``'hankel'``: Hankel-transform based implementation.
+
+    **kwargs :
+        Additional arguments passed to the backend (e.g. ``nquad``, ``scale``).
 
     Notes
     -----
     Both backends return kernels normalized for :math:`\\kappa = 1`. Any
     physical interaction strength should be applied by the caller.
     """
-    chosen = (method or "gausslag").strip().lower()
+    chosen = (method or "gausslegendre").strip().lower()
     if chosen in {"gausslag", "gauss-lag", "gausslaguerre", "gauss-laguerre", "gl"}:
-        return get_exchange_kernels_GaussLag(G_magnitudes, G_angles, nmax)
+        return get_exchange_kernels_GaussLag(G_magnitudes, G_angles, nmax, **kwargs)
     if chosen in {"hankel", "hk"}:
-        return get_exchange_kernels_hankel(G_magnitudes, G_angles, nmax)
-    raise ValueError(f"Unknown exchange-kernel method: {method!r}. Use 'gausslag' or 'hankel'.")
+        return get_exchange_kernels_hankel(G_magnitudes, G_angles, nmax, **kwargs)
+    if chosen in {"gausslegendre", "gauss-legendre", "legendre", "leg"}:
+        return get_exchange_kernels_GaussLegendre(G_magnitudes, G_angles, nmax, **kwargs)
+    raise ValueError(f"Unknown exchange-kernel method: {method!r}. Use 'gausslegendre', 'gausslag', or 'hankel'.")
 
 
 __all__ = [
     "get_form_factors",
     "get_exchange_kernels",
     "get_exchange_kernels_GaussLag",
     "get_exchange_kernels_hankel",
+    "get_exchange_kernels_GaussLegendre",
 ]
 
diff --git a/src/quantumhall_matrixelements/exchange_hankel.py b/src/quantumhall_matrixelements/exchange_hankel.py
@@ -75,6 +75,7 @@ def get_exchange_kernels_hankel(
     G_magnitudes: "RealArray",
     G_angles: "RealArray",
     nmax: int,
+    **kwargs,
 ) -> "ComplexArray":
     """Compute X_{n1,m1,n2,m2}(G) via Hankel transforms (κ=1 convention).
 
diff --git a/src/quantumhall_matrixelements/exchange_legendre.py b/src/quantumhall_matrixelements/exchange_legendre.py
@@ -0,0 +1,186 @@
+"""Exchange kernels via Gauss-Legendre quadrature with rational mapping."""
+from __future__ import annotations
+
+from functools import lru_cache
+from typing import TYPE_CHECKING
+
+import numpy as np
+import scipy.special as sps
+from scipy.special import roots_legendre
+
+if TYPE_CHECKING:
+    from numpy.typing import NDArray
+
+    ComplexArray = NDArray[np.complex128]
+    RealArray = NDArray[np.float64]
+
+
+def _N_order(n1: int, m1: int, n2: int, m2: int) -> int:
+    return (n1 - m1) - (m2 - n2)
+
+
+def _parity_factor(N: int) -> int:
+    """(-1)^((N+|N|)/2) → (-1)^N for N>=0, and 1 for N<0."""
+    return (-1) ** ((N + abs(N)) // 2)
+
+
+@lru_cache(maxsize=None)
+def _logfact(n: int) -> float:
+    return float(sps.gammaln(n + 1))
+
+
+def _C_and_indices(n1: int, m1: int, n2: int, m2: int):
+    """Constants and Laguerre parameters for f_{n1,m1} * f_{m2,n2}."""
+    p, d1 = min(n1, m1), abs(n1 - m1)
+    q, d2 = min(m2, n2), abs(m2 - n2)
+    logC = 0.5 * ((_logfact(p) - _logfact(p + d1)) + (_logfact(q) - _logfact(q + d2)))
+    C = np.exp(logC)
+    return C, p, d1, q, d2
+
+
+@lru_cache(maxsize=None)
+def _legendre_nodes_weights_mapped(nquad: int, scale: float):
+    """
+    Gauss-Legendre nodes/weights mapped from [-1, 1] to [0, inf).
+    Mapping: z = scale * (1+x)/(1-x)
+    Jacobian: dz/dx = scale * 2/(1-x)^2
+    """
+    x, w_leg = roots_legendre(nquad)
+    denom = 1.0 - x
+    z = scale * (1.0 + x) / denom
+    w = w_leg * (scale * 2.0 / (denom * denom))
+    return z, w
+
+
+def get_exchange_kernels_GaussLegendre(
+    G_magnitudes,
+    G_angles,
+    nmax: int,
+    *,
+    potential: str = "coulomb",
+    kappa: float = 1.0,
+    V_of_q=None,
+    nquad: int = 1000,
+    scale: float = 0.5,
+    ell: float = 1.0,
+) -> "ComplexArray":
+    """Compute X_{n1,m1,n2,m2}(G) using Gauss-Legendre quadrature with rational mapping.
+
+    This backend maps the semi-infinite radial integral to the finite interval [-1, 1]
+    using the rational mapping z = scale * (1+x)/(1-x). It avoids the numerical instability
+    of Gauss-Laguerre quadrature for large quantum numbers while remaining faster than
+    Hankel transforms.
+
+    Parameters
+    ----------
+    G_magnitudes, G_angles :
+        Arrays of the same shape describing |G| and polar angle θ_G.
+    nmax :
+        Number of Landau levels.
+    potential :
+        Either ``'coulomb'`` (default) or ``'general'``.
+    kappa :
+        Interaction strength prefactor.
+    V_of_q :
+        Callable ``V_of_q(q) -> V(q)`` used when ``potential='general'``.
+    nquad :
+        Number of quadrature points (default 1000).
+    scale :
+        Mapping scale factor (default 0.5). Controls the distribution of points.
+        Smaller values cluster points near the peak of the integrand for large n.
+    ell :
+        Magnetic length ℓ_B (default 1.0).
+
+    Returns
+    -------
+    Xs : (nG, nmax, nmax, nmax, nmax) complex array
+    """
+    G_magnitudes = np.asarray(G_magnitudes, dtype=float)
+    G_angles = np.asarray(G_angles, dtype=float)
+    if G_magnitudes.shape != G_angles.shape:
+        raise ValueError("G_magnitudes and G_angles must have the same shape.")
+    nG = G_magnitudes.size
+
+    Gscaled = G_magnitudes * float(ell)
+    Xs = np.zeros((nG, nmax, nmax, nmax, nmax), dtype=np.complex128)
+
+    # Get mapped grid
+    z, w = _legendre_nodes_weights_mapped(nquad, scale)
+
+    # Precompute Bessel functions J_N(sqrt(2z)*G)
+    # We cache by absN to avoid recomputing
+    J_cache: dict[int, np.ndarray] = {}
+
+    # Cache for Laguerre evaluations
+    # We evaluate L_n^d(z) for many n, d.
+    # Since n, d are small integers, we can just compute on the fly or use sps.eval_genlaguerre
+    # sps.eval_genlaguerre is efficient enough.
+
+    for n1 in range(nmax):
+        for m1 in range(nmax):
+            for n2 in range(nmax):
+                for m2 in range(nmax):
+                    N = _N_order(n1, m1, n2, m2)
+                    absN = abs(N)
+                    C, p, d1, q, d2 = _C_and_indices(n1, m1, n2, m2)
+
+                    # Compute radial integral
+                    if potential == "coulomb":
+                        # Integrand factor: exp(-z) * z^alpha * L * L * J
+                        # alpha = (d1 + d2 - 1) / 2
+                        alpha = 0.5 * (d1 + d2 - 1)
+                        
+                        L1 = sps.eval_genlaguerre(p, d1, z)
+                        L2 = sps.eval_genlaguerre(q, d2, z)
+                        
+                        # Bessel part
+                        if absN not in J_cache:
+                            arg = np.sqrt(2.0 * z)[None, :] * Gscaled[:, None]
+                            J_cache[absN] = sps.jv(absN, arg)
+                        J_abs = J_cache[absN]
+                        
+                        # Full integrand term (excluding J and weights)
+                        # exp(-z) handles x->1 decay
+                        # z^alpha handles x->-1 behavior
+                        term = np.exp(-z) * (z**alpha) * L1 * L2
+                        
+                        # Sum over quadrature points
+                        # J_abs is (nG, nquad), term is (nquad,), w is (nquad,)
+                        # Result is (nG,)
+                        radial = (J_abs * term) @ w
+                        
+                        signN = _parity_factor(N)
+                        phase_factor = (1j) ** (d1 - d2)
+                        pref = (kappa * C / np.sqrt(2.0)) * phase_factor
+                        
+                    else:
+                        # General potential
+                        if not callable(V_of_q):
+                            raise ValueError("Provide V_of_q for potential='general'")
+                        
+                        alpha = 0.5 * (d1 + d2)
+                        L1 = sps.eval_genlaguerre(p, d1, z)
+                        L2 = sps.eval_genlaguerre(q, d2, z)
+                        
+                        qvals = np.sqrt(2.0 * z) / float(ell)
+                        Veff = V_of_q(qvals) / (2.0 * np.pi * float(ell) ** 2)
+                        
+                        if absN not in J_cache:
+                            arg = np.sqrt(2.0 * z)[None, :] * Gscaled[:, None]
+                            J_cache[absN] = sps.jv(absN, arg)
+                        J_abs = J_cache[absN]
+                        
+                        term = np.exp(-z) * (z**alpha) * L1 * L2 * Veff
+                        radial = (J_abs * term) @ w
+                        
+                        signN = _parity_factor(N)
+                        phase_factor = (1j) ** (d1 - d2)
+                        pref = C * phase_factor
+
+                    phase = np.exp(-1j * N * G_angles)
+                    Xs[:, n1, m1, n2, m2] = (pref * phase) * (signN * radial)
+
+    return Xs
+
+
+__all__ = ["get_exchange_kernels_GaussLegendre"]
diff --git a/tests/test_exchange_legendre.py b/tests/test_exchange_legendre.py
@@ -0,0 +1,32 @@
+import numpy as np
+import pytest
+from quantumhall_matrixelements import get_exchange_kernels_GaussLegendre, get_exchange_kernels
+
+def test_legendre_basic_shape():
+    nmax = 2
+    Gs_dimless = np.array([0.0, 1.0])
+    thetas = np.array([0.0, np.pi])
+    X = get_exchange_kernels_GaussLegendre(Gs_dimless, thetas, nmax, nquad=100)
+    assert X.shape == (2, nmax, nmax, nmax, nmax)
+    assert np.isfinite(X).all()
+
+def test_legendre_vs_hankel_small_n():
+    """Verify agreement with Hankel for small n."""
+    nmax = 2
+    Gs_dimless = np.array([0.5, 1.5])
+    thetas = np.array([0.0, 0.2])
+    
+    X_leg = get_exchange_kernels(Gs_dimless, thetas, nmax, method="gausslegendre", nquad=500)
+    X_hk = get_exchange_kernels(Gs_dimless, thetas, nmax, method="hankel")
+    
+    assert np.allclose(X_leg, X_hk, rtol=1e-3, atol=1e-3)
+
+def test_legendre_large_n_stability():
+    """Verify that it runs without error for large n (where gausslag fails)."""
+    nmax = 15
+    Gs_dimless = np.array([1.0])
+    thetas = np.array([0.0])
+    
+    # This should not raise an error
+    X = get_exchange_kernels_GaussLegendre(Gs_dimless, thetas, nmax, nquad=500)
+    assert np.isfinite(X).all()
diff --git a/tests/test_validation.py b/tests/test_validation.py
