Merge pull request #67 from mreineck/ducc_negative_spins

Support negative spins in ducc back-end

Author: mdavezac

src/pyssht/ducc_demo.py

@@ -82,7 +82,7 @@ def test_rot(L, nthreads=1):
 for Method in ["MW", "MWSS", "GL", "DH"]:
     for L in L_list:
         for Reality in [False, True]:
-            for Spin in [0] if Reality else [0, 1]:
+            for Spin in [0] if Reality else [-1, 0, 1]:
                 res = test_SHT(L, Method, Reality, Spin, nthreads=nthreads)
                 print(
                     "{:4}, L={:4}, Reality={:5}, Spin={:1}:  L2 error={:e}, speedup factor={:f}".format(

src/pyssht/ducc_interface.pyx

@@ -57,50 +57,70 @@ def _build_real_flm(alm, Py_ssize_t L):
         ofs += L-m
     return res
 
-cdef _extract_complex_alm(flm, Py_ssize_t L):
+cdef _extract_complex_alm(flm, Py_ssize_t L, Py_ssize_t Spin):
     res = np.empty((2, _nalm(L-1, L-1),), dtype=np.complex128)
     cdef Py_ssize_t ofs=0, m, i
-    cdef double mfac
+    cdef double mfac, sfac=(-1)**abs(Spin)
     cdef complex[:,:] myres=res
     cdef complex[:] myflm=flm
     cdef Py_ssize_t[:] lidx = _get_lidx(L)
     cdef complex fp, fm
-    for m in range(L):
-        mfac = (-1)**m
-        for i in range(m,L):
-            fp = myflm[lidx[i]+m]
-            fm = mfac * (myflm[lidx[i]-m].real - 1j*myflm[lidx[i]-m].imag)
-            myres[0, ofs-m+i] = 0.5*(fp+fm)
-            myres[1, ofs-m+i] = -0.5j*(fp-fm)
-        ofs += L-m
+    if Spin >= 0:
+        for m in range(L):
+            mfac = (-1)**m
+            for i in range(m,L):
+                fp = myflm[lidx[i]+m]
+                fm = mfac * (myflm[lidx[i]-m].real - 1j*myflm[lidx[i]-m].imag)
+                myres[0, ofs-m+i] = 0.5*(fp+fm)
+                myres[1, ofs-m+i] = -0.5j*(fp-fm)
+            ofs += L-m
+    else:
+        for m in range(L):
+            mfac = (-1)**m
+            for i in range(m,L):
+                fp = mfac*sfac*(myflm[lidx[i]-m].real - 1j*myflm[lidx[i]-m].imag)
+                fm = sfac*myflm[lidx[i]+m]
+                myres[0, ofs-m+i] = 0.5*(fp+fm)
+                myres[1, ofs-m+i] = -0.5j*(fp-fm)
+            ofs += L-m
     return res
 
-cdef _build_complex_flm(alm, Py_ssize_t L):
+cdef _build_complex_flm(alm, Py_ssize_t L, Py_ssize_t Spin):
     res = np.empty((L*L), dtype=np.complex128)
     cdef Py_ssize_t ofs=0, m, i
     cdef complex fp, fm
     cdef complex[:] myres=res
     cdef complex[:,:] myalm=alm
     cdef Py_ssize_t[:] lidx = _get_lidx(L)
-    cdef double mfac
-    for m in range(L):
-        mfac = (-1)**m
-        for i in range(m,L):
-            fp = myalm[0, ofs-m+i] + 1j*myalm[1, ofs-m+i]
-            fm = myalm[0, ofs-m+i] - 1j*myalm[1, ofs-m+i]
-            myres[lidx[i]+m] = fp
-            myres[lidx[i]-m] = mfac*(fm.real - 1j*fm.imag)
-        ofs += L-m
+    cdef double mfac, sfac=(-1)**abs(Spin)
+    if Spin >= 0:
+        for m in range(L):
+            mfac = (-1)**m
+            for i in range(m,L):
+                fp = myalm[0, ofs-m+i] + 1j*myalm[1, ofs-m+i]
+                fm = myalm[0, ofs-m+i] - 1j*myalm[1, ofs-m+i]
+                myres[lidx[i]+m] = fp
+                myres[lidx[i]-m] = mfac*(fm.real - 1j*fm.imag)
+            ofs += L-m
+    else:
+        for m in range(L):
+            mfac = (-1)**m
+            for i in range(m,L):
+                fp = myalm[0, ofs-m+i] + 1j*myalm[1, ofs-m+i]
+                fm = myalm[0, ofs-m+i] - 1j*myalm[1, ofs-m+i]
+                myres[lidx[i]+m] = sfac*fm
+                myres[lidx[i]-m] = sfac*mfac*(fp.real -1j *fp.imag)
+            ofs += L-m
     return res
 
 
 def rotate_flms(flm, alpha, beta, gamma, L, int nthreads = 1):
     ducc0 = import_ducc0()
-    alm = _extract_complex_alm(flm, L)
+    alm = _extract_complex_alm(flm, L, 0)
     for i in range(2):
         alm[i] = ducc0.sht.rotate_alm(
             alm[i], L-1, gamma, beta, alpha, nthreads=nthreads)
-    return _build_complex_flm(alm, L)
+    return _build_complex_flm(alm, L, 0)
 
 
 def inverse(np.ndarray flm, Py_ssize_t L, Py_ssize_t Spin, str Method, bint Reality, int nthreads = 1):
@@ -121,20 +141,20 @@ def inverse(np.ndarray flm, Py_ssize_t L, Py_ssize_t Spin, str Method, bint Real
             spin=0,
             geometry=gdict[Method])[0]
     elif Spin == 0:
-        alm = _extract_complex_alm(flm, L)
+        alm = _extract_complex_alm(flm, L,0)
         flmr = _build_real_flm(alm[0], L)
         flmi = _build_real_flm(alm[1], L)
         return inverse(flmr, L, 0, Method, True) + 1j*inverse(flmi, L, 0, Method, True)
     else:
         tmp=ducc0.sht.experimental.synthesis_2d(
-            alm=_extract_complex_alm(flm, L),
+            alm=_extract_complex_alm(flm, L, Spin),
             ntheta=ntheta,
             nphi=nphi,
             lmax=L-1,
             nthreads=nthreads,
-            spin=Spin,
+            spin=abs(Spin),
             geometry=gdict[Method])
-        res = -1j*tmp[1]
+        res = -1j*tmp[1] if Spin >=0 else 1j*tmp[1]
         res -= tmp[0]
         return res
 
@@ -155,22 +175,24 @@ def inverse_adjoint(f, L, Spin, Method, Reality, int nthreads = 1):
             spin=0,
             geometry=gdict[Method])[0], L)
     elif Spin == 0:
-       flmr = inverse_adjoint(f.real, L, Spin, Method, True)
-       flmi = inverse_adjoint(f.imag, L, Spin, Method, True)
-       alm = np.empty((2,_nalm(L-1, L-1)), dtype=np.complex128)
-       alm[0] = _extract_real_alm(flmr, L)
-       alm[1] = _extract_real_alm(flmi, L)
-       return _build_complex_flm(alm, L)
+        flmr = inverse_adjoint(f.real, L, Spin, Method, True)
+        flmi = inverse_adjoint(f.imag, L, Spin, Method, True)
+        alm = np.empty((2,_nalm(L-1, L-1)), dtype=np.complex128)
+        alm[0] = _extract_real_alm(flmr, L)
+        alm[1] = _extract_real_alm(flmi, L)
+        return _build_complex_flm(alm, L, 0)
     else:
-       map = f.astype(np.complex128).view(dtype=np.float64).reshape((f.shape[0],f.shape[1],2)).transpose((2,0,1))
-       res = _build_complex_flm(ducc0.sht.experimental.adjoint_synthesis_2d(
-           map=map,
-           lmax=L-1,
-           nthreads=nthreads,
-           spin=Spin,
-           geometry=gdict[Method]), L)
-       res *= -1
-       return res
+        map = f.astype(np.complex128).view(dtype=np.float64).reshape((f.shape[0],f.shape[1],2)).transpose((2,0,1))
+        if  Spin < 0:
+            map[1]*=-1
+        res = _build_complex_flm(ducc0.sht.experimental.adjoint_synthesis_2d(
+            map=map,
+            lmax=L-1,
+            nthreads=nthreads,
+            spin=abs(Spin),
+            geometry=gdict[Method]), L, Spin)
+        res *= -1
+        return res
 
 
 def forward(f, L, Spin, Method, Reality, int nthreads = 1):
@@ -194,15 +216,17 @@ def forward(f, L, Spin, Method, Reality, int nthreads = 1):
         alm = np.empty((2,_nalm(L-1, L-1)), dtype=np.complex128)
         alm[0] = _extract_real_alm(flmr, L)
         alm[1] = _extract_real_alm(flmi, L)
-        return _build_complex_flm(alm, L)
+        return _build_complex_flm(alm, L, 0)
     else:
         map = f.astype(np.complex128).view(dtype=np.float64).reshape((f.shape[0],f.shape[1],2)).transpose((2,0,1))
+        if Spin < 0:
+            map[1]*=-1
         res = _build_complex_flm(ducc0.sht.experimental.analysis_2d(
             map=map,
             lmax=L-1,
             nthreads=nthreads,
-            spin=Spin,
-            geometry=gdict[Method]), L)
+            spin=abs(Spin),
+            geometry=gdict[Method]), L, Spin)
         res *= -1
         return res
 
@@ -225,19 +249,19 @@ def forward_adjoint(np.ndarray flm, Py_ssize_t L, Py_ssize_t Spin, str Method, b
             spin=0,
             geometry=gdict[Method])[0]
     elif Spin == 0:
-        alm = _extract_complex_alm(flm, L)
+        alm = _extract_complex_alm(flm, L, 0)
         flmr = _build_real_flm(alm[0], L)
         flmi = _build_real_flm(alm[1], L)
         return forward_adjoint(flmr, L, 0, Method, True) + 1j*forward_adjoint(flmi, L, 0, Method, True)
     else:
         tmp=ducc0.sht.experimental.adjoint_analysis_2d(
-            alm=_extract_complex_alm(flm, L),
+            alm=_extract_complex_alm(flm, L, Spin),
             ntheta=ntheta,
             nphi=nphi,
             lmax=L-1,
             nthreads=nthreads,
-            spin=Spin,
+            spin=abs(Spin),
             geometry=gdict[Method])
-        res = -1j*tmp[1]
+        res = -1j*tmp[1] if Spin >=0 else 1j*tmp[1]
         res -= tmp[0]
         return res

tests/test_ducc0.py

@@ -21,7 +21,7 @@ def method(request):
     return request.param
 
 
-@fixture(params=[0, 1, 2])
+@fixture(params=[-2, -1, 0, 1, 2])
 def spin(request):
     return request.param
 
@@ -139,7 +139,7 @@ def test_complex_inverse_adjoint_ssht_vs_ducc0(
 
     try:
         ssht_adj_coeffs = ssht.inverse_adjoint(
-            complex_image, order, Reality=False, Method=method, Spin=0
+            complex_image, order, Reality=False, Method=method, Spin=spin
         )
     except ssht_input_error:
         assert method not in ("MW", "MWSS")
@@ -150,7 +150,7 @@ def test_complex_inverse_adjoint_ssht_vs_ducc0(
         order,
         Reality=False,
         Method=method,
-        Spin=0,
+        Spin=spin,
         backend="ducc",
         nthreads=nthreads,
     )