[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 (d332189) to head (098e5cc).
⚠️ 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        3544     3611      +67     
==========================================
+ Hits         3406     3481      +75     
+ 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]

Hi @mwyau thanks for your comment and apologies for the delay getting to responding to this.

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

In general we would like to improve the range of sampling schemes we support and in particular support more schemes commonly used in climate / atmospheric science settings. Ideally to make this easier we need to do some refactoring of the current sampling scheme / grid support to make it more extensible (#340)- this PR makes a few incidental changes along this line e.g. adding in functions to indicate which properties different schemes have rather than listing separately each time we branch on these properties, but hopefully we'll do a broader refactor as a follow up. In particular to avoid issues like you mention around around the specific choice of grid size we use for sampling="gl" we'd like to allow freedom in choosing the grid size for a given quadrature rule / set of sample locations, under some constraints on grid sizes allow exact round-tripping through transforms were relevant.

[matt-graham]

Following previous discussion about reviewing this and the inability to get the Claude Review integration working due to lack of usage credits I got a local Claude instance to review the changes in the PR and flag any inconsistencies and potential readability improvements. The main things flagged and hopefully now fixed:

• For new sampling = "cc" scheme, while the quadrature weight calculation special cased L = 1 which corresponds to ntheta = 1 the theta sample positions would compute as NaN in this case. Given we define CC scheme as including both poles as so requiring n_theta >= 2 and correspondingly L > 1 if this is violated an error is now raised early.
• The sets defining which schemes include poles / are equiangular in theta did not make it clear why these properties were important and required repetitively converting the sampling argument to lowercase. These have now been refactored in to boolean valued set membership functions with explanatory docstrings.
• Various minor inconsistencies in docstrings and in particular lists of sampling schemes supported have been fixed.

[mscroggs]

Question probably for an issue unless it's been discussed before: would an enum be better than string comparisons here? Additionally, an enum would mean a user could use tab completion to see what values sampling can take

[matt-graham]

Agreed that string comparisons here and the pattern of always converting to lowercase should be improved! An enum would be nicer in terms of being explicit what the allowable values are and making this tab-completable so could be a first step.

In the longer run I would ideally like to move to a more modular / object-oriented sampling scheme interface (#340) with each scheme implementing a common set of methods which we dispatch to based on type rather than the current approach of having lots of functions which branch on a string argument, as this makes it less clear which functions need to be extended to add a new scheme and means such additions touch lots of codebase rather than being contained to e.g. a single new class. Some of the changes in this PR were with this in mind, and try to move more of the sample scheme specific logic inside functions within one module, though there are still some remaining branches on the sampling argument in the transform code paths.