Add support for Clenshaw-Curtis and Fejér second rule quadrature / sampling schemes

[matt-graham]

Should resolve #251 and also does some work towards #340

Adds implementations of the Clenshaw-Curtis and Fejér's (second) quadrature rules and corresponding sampling schemes.

The quadrature weights are computed using the FFT based approach described in

Waldvogel, J. (2006). Fast construction of the Fejer and Clenshaw-Curtis quadrature rules.
BIT Numerical Mathematics, 46(1), 195–202. https://doi.org/10.1007/s10543-006-0045-4.

Following the definitions in that paper the Clenshaw-Curtis scheme includes samples at both poles ($\theta = \pm \pi$) while the Fejér's second rule scheme excludes both poles. The Clenshaw-Curtis scheme is selected with sampling="cc" and Fejér's second rule with sampling="f2" arguments to relevant functions.

To facilitate adding these new schemes, this PR also does some light refactoring of the existing quadrature and sample related modules to try to make the code more DRY and with a more standardized interfaced across the different schemes. In some cases the rationale for special casing code paths on different sampling schemes was a bit unclear, particularly around handling of pole singularities and setting m_offset, so I've done a best efforts attempt at trying to rationalise this but this probably could do with some careful checking.

Currently for a bandlimit $L$ both schemes have $N_\theta = 2L - 1$ and $N_\phi = 2L$. The former appears to be required to be $\sim 2L$ to give round-trip errors close to machine precision; I believe this related to the comment

An important difference [of Clenshaw-Curtis and Fejér rules] from Gaussian quadrature is that the number of nodes, given the truncation total wavenumber, $N$ required for alias-free exact meridional integration, is
$J \geq 2N + 1$,
which is about twice that in the Gaussian case. This is because Clenshaw–Curtis quadrature expands the integrand into Chebyshev polynomials of up to $(J −1)\text{th}$ degree and thus $J − 1$ needs to be no less than $2N$, the maximum degree of the integrand polynomial.

in the paper Hotte and Ujiie (2018) A nestable, multigrid-friendly grid on a sphere for global spectral models based on Clenshaw–Curtis quadrature with their definition of $N$ corresponding to $L - 1$ and $J$ to $N_\theta$ in our notation and hence we require $N_\theta \geq 2 L - 1$.

Choosing $N_\theta = 2L - 1$ gives an odd number of nodes in the latitude axis and so includes the equator as a node, this matching the co-latitude points used in Hotte and Ujiie (2018). This still leaves $N_\phi$ undetermined - here I have chosen $N_\phi = 2L$ which is sufficient to avoid aliasing in FFT operations (which requires $N_\phi > 2L - 2$), and matches the MW scheme, but gives $N_\theta \approx N_\phi$ while it seems quite common in grids used in practice in climate applications to have $N_\phi \approx 2N_\theta$ which might suggest using $N_\phi = 4L$.

@jasonmcewen I'm tagging you as you mentioned you would be interested in looking through the details of this.

TODO:

• Decide if we want to adjust spherical grid dimensions for given bandlimit $L$ for new schemes
• Add details of new schemes to sampling scheme documentation overview page
• Check if new schemes can be carried through to Wigner transforms + document

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 96.39%. Comparing base (06c2363) to head (1c08016).
⚠️ Report is 1 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #381      +/-   ##
==========================================
+ Coverage   96.10%   96.39%   +0.29%     
==========================================
  Files          34       34              
  Lines        3539     3603      +64     
==========================================
+ Hits         3401     3473      +72     
+ Misses        138      130       -8     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[matt-graham]

@jasonmcewen I arrived at these by trial and error as I could not see the pattern in how the values for the other schemes was derived - if there some underlying relationship here it would be good to document. It might also be worth moving this to one of modules under s2fft.sampling

[matt-graham]

Copying discussion from Slack here so we have a record going forwards:

The Hotte and Ujiie (2018) paper suggests we require n_theta >= 2 * L - 1 to get round trip errors close to machine precision, which underlies the choice in my current PR to use n_theta = 2 * L - 1 (along with n_phi = 2 * L for the Clenshaw-Curtis sampling="cc" and Fejer second rule sampling="f2" sampling schemes.

The analysis_2d function in ducc0 however only requires that n_theta >= nrings_min = lmax + 2 = L + 1 for the Clenshaw-Curtis geometry="CC" scheme (and n_theta >= nrings_min = 2 * lmax + 1 = 2 * L - 1 for the Fejer second rule geometry="F2" scheme), with L = lmax + 1.

Generating a set of random coefficients then running ducc0.sht.synthesis_2d → ducc0.sht.analysis_2d to generate map then recover coefficients, using geometry="CC" with ntheta = L + 1 we do get round-trip errors close to machine precision $\sim 10^{-15}$. Running our corresponding S2FFT implementations s2fft.inverse → s2fft.forward with sampling="cc" and n_theta = L + 1 (enforced by manually editing n_theta function in source code as we don't allow passing this as an argument unlike ducc0) we however get $\sim 1$ round trip errors. The outputs of ducc0.sht.synthesis_2d and s2fft.inverse for equivalent input coefficients do however match to within machine precision with sampling="cc" which confirms as expected this is to do with the implementation of the quadrature specifically.

Looking at the implementation of ducc0.sht.analysis_2d specifically the lines 1763-1791 in commit 6c871af we can see there is a distinct branch for geometry in {"CC", "F1", "MW", "MWflip"} which I think is all the schemes with minimum n_theta ~ L rather than 2 * L. This calls out to a function resample_to_prepared_CC which I suspect is doing something similar to the upsample_by_two_mwss calls within the s2fft.forward implementation. Appendix A in the paper Reinecke, Belkner and Carron (2023) I think explains the double Fourier sphere technique being used here in more detail.

So if we want to support n_theta = L + 1 in S2FFT I think we will need an upsampling / resampling implementation comparable to what we have for MW/MWSS and matching what is done within ducc0 resampling_to_prepared_CC scheme.

[matt-graham]

@jasonmcewen from my perspective this is now ready to merge. I have addressed the points discussed when I did the code walkthrough / review with you and Kevin, most significantly checking the quadrature rule tests for exact integration of polynomial functions tests the rules at the relevant minimum number of points, and also clarified the logic behind m_offset being equal to zero or one. If we want to support the $N_\theta \sim L$ grid dimensions for CC and F2 schemes that ducc0 supports using the double Fourier sphere approach I'd say that would be best to do in a follow up PR.

[mreineck]

Appendix A in the paper Reinecke, Belkner and Carron (2023) I think explains the double Fourier sphere technique being used here in more detail.

I can confirm that this is exactly how it is done. The beautiful aspect about this is that, after upsampling and applying the weights, you can downsample again and still carry out the Legendre transform part on the low number of rings.

[matt-graham]

@claude review

[claude]

⚠️ Code review skipped — your organization has no extra usage available to pay for this review.

If your organization's extra usage balance is empty, an organization admin can add extra usage credits at claude.ai/admin-settings/usage. If its monthly spend limit was reached, an admin can raise it on the same page. If neither applies, contact Anthropic support.

Once extra usage is available, comment @claude review on this pull request to trigger a review.

[mwyau]

Just stumbled upon this. My downstream package spharmgrid uses ducc0 as the CPU backend, and I am recently adding torch-harmonics and s2fft backends.

For common atmospheric science data, it is common to use (L+1) x 2L fixed grids, for example 73x144 for a 2.5-degree fixed lat-lon grid. For Gaussian grids, it is often L x 2L such as 64x128.

I am using the mwss sampling scheme in my s2fft backend for the fixed lat-lon "cc" grid. A proper CC sampling scheme support will be great. I actually implemented the Reinecke, Belkner and Carron (2023) method in PyTorch. I can port it to s2fft in a follow-up PR after this is merged.

As a side note, the current s2fft gl sampling scheme L x (2L-1) does not work for my use case. I am currently using a workaround in my package to preprocess longitude rings. Any interest in adding L x 2L support for gl sampling scheme? Issue: #406