Can't use negative spins with ducc

[moble]

When trying to transform via ducc with negative spins, I get an error.

For example, I modify this line of ducc_demo.py to include -1 as one of the Spin values like so

            for Spin in [0] if Reality else [-1, 0, 1]:

then I get this error:

  File "cpyssht.pyx", line 529, in pyssht.cpyssht.inverse
  File "ducc_interface.pyx", line 129, in pyssht.ducc_interface.inverse
TypeError: synthesis_2d(): incompatible function arguments. The following argument types are supported:
    1. (*, alm: numpy.ndarray, spin: int, lmax: int, geometry: str, ntheta: object = None, nphi: object = None, mmax: object = None, nthreads: int = 1, map: object = None) -> numpy.ndarray

I suspect that Py_ssize_t Spin should actually be changed to int Spin here and here. (I believe Py_ssize_t is a signed typedef in CPython, but it looks like Cython reinterprets them as unsigned when interfacing with C code.)

[jasonmcewen]

Thanks for your comment @moble !

@mreineck I think you added these functions to support ducc. Do the suggested changes look right to you? Is there anywhere else things should be revised?

[mreineck]

Ducc only accepts nonnegative spins for its transforms at the moment. But it should be possible to support negative spins by taking the complex conjugate of both input and output, and use the absolute spin value for the transform, correct? If you agree that this is a viable solution, I can try to implement this.

[mreineck]

BTW, the error message is not due to cython; it's produced by ducc's pybind11 interface, which knows that the accepted spins must be convertible to an unsigned size_t.

[moble]

it should be possible to support negative spins by taking the complex conjugate of both input and output, and use the absolute spin value for the transform, correct?

Unless I'm misunderstanding what you're saying, I don't think that's right. I think you do have to conjugate both input and output, but also multiply each mode by (-1.)**(s+m) and swap m and -m modes. Here's an MWE using s=-2 (the key use case for gravitational-wave astronomy):

import numpy as np
import pyssht

s = -2
lmax = 8
L = lmax+1
backend = "SSHT"
method = "MW"

f_lm = np.random.randn(2 * L**2).view(dtype=complex)
f_lm[:pyssht.elm2ind(abs(s), -abs(s))] = 0

modified_f_lm = f_lm.copy()
for ell in range(abs(s), lmax+1):
    for m in range(-ell, ell+1):
        modified_f_lm[pyssht.elm2ind(ell, m)] = (-1.)**(s+m) * f_lm[pyssht.elm2ind(ell, -m)].conjugate()

f1 = pyssht.inverse(f_lm, L, Spin=s)
f2 = pyssht.inverse(f_lm.conjugate(), L, Spin=-s).conjugate()
f3 = pyssht.inverse(modified_f_lm, L, Spin=-s).conjugate()

np.max(np.abs((f1-f2))), np.max(np.abs((f1-f3)))

Here, f1 is computed the usual way, f2 by just conjugating input and output with spin -s, and f3 by also multiplying by (-1.)**(s+m) and swapping m and -m modes. The difference between f1 and f2 is of order unity; between f1 and f3 is of order 10x machine epsilon.

Similarly, starting with the function values f1 and converting to mode weights:

f_lm1 = pyssht.forward(f1, L, Spin=s)
f_lm2 = pyssht.forward(f1.conjugate(), L, Spin=-s).conjugate()
f_lm3 = pyssht.forward(f1.conjugate(), L, Spin=-s)

modified_f_lm3 = f_lm3.copy()
for ell in range(abs(s), lmax+1):
    for m in range(-ell, ell+1):
        modified_f_lm3[pyssht.elm2ind(ell, m)] = (-1.)**(s+m) * f_lm3[pyssht.elm2ind(ell, -m)].conjugate()

np.max(np.abs((f_lm-f_lm1))), np.max(np.abs((f_lm-f_lm2))), np.max(np.abs((f_lm-modified_f_lm3)))

Again, only f_lm1 and modified_f_lm3 agree with the original data.

[jasonmcewen]

Ok, thanks for the further comments @mreineck about ducc support for negative spin.

@moble I haven't had a chance to review the code you posted but the conjugate symmetry expression for spin functions is as follows (copied from column 1 on page 4 of this paper):

@moble We are also working on a new spin spherical harmonic transform code in JAX, supporting both HEALPix and equiangular sampling, that is differentiable and runs on accelerators, e.g. GPUs/TPUs. We hope to release the new JAX code by the end of the calendar year so please watch this space in case useful for your research...

[moble]

but the conjugate symmetry expression for spin functions is as follows

I agree with that, but since it's only true for functions satisfying ${}_sf^ = {}_{-s}f$ we have to generalize it for arbitrary functions. That generalization — specifically the transformation ${}_sf_{\ell,m} \to (-1)^{s+m}{}_{-s}f_{\ell,-m}^$ — is exactly what the code above implements.

We are also working on a new spin spherical harmonic transform code in JAX, supporting both HEALPix and equiangular sampling, that is differentiable and runs on accelerators, e.g. GPUs/TPUs.

Thanks @jasonmcewen, that's very interesting. Differentiability and running on accelerators are both vital to what I do. I'm actually just extending my Julia package to use something that's grown out of the optimal-dimensionality result you were involved in.

[mreineck]

Thanks @moble for the code! You are right of course, I completely forgot about the (-1)**m component.

I need to think a bit about the best place to put the necessary conversions, but I hope to have a patch ready soon.

We are also working on a new spin spherical harmonic transform code in JAX,

Once there is a public code base, I'd also be very interested to have a look, @jasonmcewen!

[mreineck]

Tentative fix in #67. I don't understand the CI failure of the development build, but it doesn't seem related to this issue.

[mreineck]

Accidental discovery: in pure Python (-1)**(-1) results in -1, in Cython the same expression evaluates to 0(!). This took some time to debug ...