Merge pull request #58 from mreineck/forward_adjoint

mdavezac · web-flow · commit d4d3d94f31ac · 2021-07-25T16:02:29.000+01:00
Draft: add ducc support for forward adjoint
diff --git a/.github/workflows/python.yml b/.github/workflows/python.yml
@@ -118,14 +118,14 @@ jobs:
           mv source-distribution/*.tar.gz  wheel-*/*.whl dist
 
       - name: Publish distribution 📦 to Test PyPI
-        if: ${{ github.ref != 'refs/tags/v1.4.0' }}
+        if: ${{ github.ref != 'refs/tags/v1.5.0' }}
         uses: pypa/gh-action-pypi-publish@master
         with:
           password: ${{ secrets.TEST_PYPI_TOKEN }}
           repository_url: https://test.pypi.org/legacy/
 
       - name: Publish distribution 📦 to PyPI
-        if: ${{ github.ref == 'refs/tags/v1.4.0' }}
+        if: ${{ github.ref == 'refs/tags/v1.5.0' }}
         uses: pypa/gh-action-pypi-publish@master
         with:
           password: ${{ secrets.PYPI_TOKEN }}
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.12)
 
 project(
     Ssht
-    VERSION "1.4.0"
+    VERSION "1.5.0"
     DESCRIPTION "Fast and exact spin spherical harmonic transforms"
     HOMEPAGE_URL "http://astro-informatics.github.io/ssht/"
     LANGUAGES C)
diff --git a/README.md b/README.md
@@ -23,7 +23,7 @@ spin spherical harmonic transforms based on the sampling theorem on the
 sphere derived in <a href="http://www.jasonmcewen.org/publication/mcewen-fssht/">McEwen & Wiaux (2011)</a>.
 
 **SSHT** can also interface with [ducc0](https://pypi.org/project/ducc0/) and
-*use it as a backend for the forward, inverse and inverse adjoint transforms.
+use it as a backend for spherical harmonic transforms and rotations.
 
 
 ## INSTALLATION
@@ -82,4 +82,4 @@ SSHT is released under the GPL-3 license.  For further details see
 
 ## AUTHORS
 
-SSHT was initially written by [Jason McEwen](http://www.jasonmcewen.org/) but significant contributors have since been made by a number of <a href="https://github.com/astro-informatics/ssht/graphs/contributors">others</a>.
+SSHT was initially written by [Jason McEwen](http://www.jasonmcewen.org/) but significant contributors have since been made by a number of <a href="https://github.com/astro-informatics/ssht/graphs/contributors">others</a>.
diff --git a/docs/index.html b/docs/index.html
@@ -51,6 +51,8 @@ <h2>News</h2>
 
       <ul>
 
+        <li><strong>July 2021</strong><br />
+          Release of SSHT 1.5.0 (<a href="https://pypi.org/project/ducc0/">ducc0</a> forward adjoint)
         <li><strong>July 2021</strong><br />
           Release of SSHT 1.4.0 (<a href="https://pypi.org/project/ducc0/">ducc0</a>interface)
         <li><strong>February 2021</strong><br />
diff --git a/makefile b/makefile
@@ -1,7 +1,7 @@
 # ======== COMPILER ========
 
 CC      = gcc
-OPT		= -std=c99 -pedantic -Wall -O3 -fopenmp -DSSHT_VERSION=\"1.4.0\" -DSSHT_BUILD=\"`git rev-parse HEAD`\"
+OPT		= -std=c99 -pedantic -Wall -O3 -fopenmp -DSSHT_VERSION=\"1.5.0\" -DSSHT_BUILD=\"`git rev-parse HEAD`\"
 #OPT	= -Wall -g -fopenmp -DSSHT_VERSION=\"1.0b1\" -DSSHT_BUILD=\"`git rev-parse HEAD`\"
 
 
diff --git a/setup.cfg b/setup.cfg
@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 1.4.0
+current_version = 1.5.0
 commit = False
 tag = False
 parse = (?P<major>\d+)\.(?P<minor>\d+)\.(?P<patch>\d+)(rc(?P<rc>\d+))?
diff --git a/setup.py b/setup.py
@@ -18,7 +18,7 @@
     "conan",
     "pip!=20.0.0,!=20.0.1",
     "pytest",
-    "ducc0>=0.16",
+    "ducc0>=0.18",
 ]
 
 long_description = (
@@ -27,7 +27,7 @@
 
 setup(
     name="pyssht",
-    version="1.4.0",
+    version="1.5.0",
     author=[
         "J. D. McEwen",
         "C. R. G. Wallis",
@@ -36,7 +36,7 @@
         "Y. Wiaux",
     ],
     install_requires=["numpy", "scipy"],
-    extras_require={"dev": dev_requirements, "ducc0": ["ducc0>=0.16"]},
+    extras_require={"dev": dev_requirements, "ducc0": ["ducc0>=0.18"]},
     description="Fast spin spherical transforms",
     long_description=long_description,
     long_description_content_type="text/markdown",
diff --git a/src/pyssht/SSHT_Python_Documentation.md b/src/pyssht/SSHT_Python_Documentation.md
@@ -98,6 +98,10 @@ Performs the adjoint of the forward spherical harmonic transform.
     1. `'MW'`         [McEwen & Wiaux sampling (default)]
     2. `'MWSS'`       [McEwen & Wiaux symmetric sampling]
 * `Reality`  determines if the signal is real or complex, Boolean (default = False)
+* `backend` the backend that runs the transforms:
+    1. `'SSHT'` this package
+    2. `'ducc'` interface to [ducc0](https://pypi.org/project/ducc0/).
+* `nthreads`: number of threads when calling into the `'ducc'` backend. Ignored otherwise.
 
 #### Output
 
@@ -122,8 +126,7 @@ Performs the adjoint of the inverse spherical harmonic transform.
 * `Reality`  determines if the signal is real or complex, Boolean (default = False)
 * `backend` the backend that runs the transforms:
     1. `'SSHT'` this package
-    2. `'ducc'` interface to [ducc0](https://pypi.org/project/ducc0/). "MW_pole"
-       is not available in this backend.
+    2. `'ducc'` interface to [ducc0](https://pypi.org/project/ducc0/).
 * `nthreads`: number of threads when calling into the `'ducc'` backend. Ignored otherwise.
 
 #### Output
diff --git a/src/pyssht/cpyssht.pyx b/src/pyssht/cpyssht.pyx
@@ -612,11 +612,15 @@ def inverse_adjoint(f, int L, int Spin=0, str Method='MW', bint Reality=False, s
     return flm
 
 
-def forward_adjoint(flm, int L, int Spin=0, str Method='MW', bint Reality=False, str backend="SSHT"):
-    from pyssht.parameters import method
+def forward_adjoint(flm, int L, int Spin=0, str Method='MW', bint Reality=False, str backend="SSHT", **kwargs):
+    from pyssht.ducc_interface import forward_adjoint as _ducc0_forward_adjoint
+    from pyssht.parameters import method, Ducc
 
     params = method(Method, spin=Spin, reality=Reality, backend=backend)
 
+    if isinstance(params, Ducc):
+        return _ducc0_forward_adjoint(flm, L, params.spin, params.method, params.reality, params.nthreads)
+
     if flm.ndim != 1:
       raise ssht_input_error('flm must be 1D numpy array')
 
diff --git a/src/pyssht/ducc_interface.pyx b/src/pyssht/ducc_interface.pyx
@@ -205,3 +205,39 @@ def forward(f, L, Spin, Method, Reality, int nthreads = 1):
             geometry=gdict[Method]), L)
         res *= -1
         return res
+
+
+def forward_adjoint(np.ndarray flm, Py_ssize_t L, Py_ssize_t Spin, str Method, bint Reality, int nthreads = 1):
+    ducc0 = import_ducc0()
+    gdict = {"DH":"F1", "MW":"MW", "MWSS":"CC", "GL":"GL"}
+    theta = _get_theta(L, Method)
+    ntheta = theta.shape[0]
+    nphi = 2*L-1
+    if Method == 'MWSS':
+        nphi += 1
+    if Reality:
+        return ducc0.sht.experimental.adjoint_analysis_2d(
+            alm=_extract_real_alm(flm, L).reshape((1,-1)),
+            ntheta=ntheta,
+            nphi=nphi,
+            lmax=L-1,
+            nthreads=nthreads,
+            spin=0,
+            geometry=gdict[Method])[0]
+    elif Spin == 0:
+        alm = _extract_complex_alm(flm, L)
+        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),
+            ntheta=ntheta,
+            nphi=nphi,
+            lmax=L-1,
+            nthreads=nthreads,
+            spin=Spin,
+            geometry=gdict[Method])
+        res = -1j*tmp[1]
+        res -= tmp[0]
+        return res
diff --git a/tests/test_ducc0.py b/tests/test_ducc0.py
@@ -1,9 +1,9 @@
 import numpy as np
-from pytest import approx, fixture, importorskip
+from pytest import approx, fixture, mark, importorskip
 
 import pyssht as ssht
 
-skip_ducc0 = importorskip("ducc0", minversion="0.16")
+skip_ducc0 = importorskip("ducc0", minversion="0.18")
 
 
 @fixture
@@ -157,6 +157,39 @@ def test_complex_inverse_adjoint_ssht_vs_ducc0(
     assert ssht_adj_coeffs == approx(ducc0_adj_coeffs)
 
 
+@mark.parametrize("method", ["MW", "MWSS"])
+def test_real_forward_adjoint(rng: np.random.Generator, method, order):
+    shape = ssht.sample_shape(order, Method=method)
+    f = rng.standard_normal(shape, dtype="float64")
+    flm = ssht.forward(f, order, Reality=True, Method=method)
+    f = ssht.inverse(flm, order, Reality=True, Method=method)
+
+    f_prime = rng.standard_normal(shape, dtype="float64")
+    flm_prime = ssht.forward(f_prime, order, Reality=True, Method=method)
+    f_prime = ssht.forward_adjoint(
+        flm_prime, order, Reality=True, Method=method, backend="ducc"
+    )
+
+    assert flm_prime.conj() @ flm == approx(f_prime.flatten().conj() @ f.flatten())
+
+
+@mark.parametrize("method", ["MW", "MWSS"])
+def test_forward_adjoint(rng: np.random.Generator, spin, method, order):
+    flm = rng.standard_normal((order * order, 2), dtype="float64") @ np.array([1, 1j])
+    flm[0 : spin * spin] = 0.0
+    f = ssht.inverse(flm, order, Spin=spin, Method=method, backend="ducc")
+
+    flm_prime = rng.standard_normal((order * order, 2), dtype="float64") @ np.array(
+        [1, 1j]
+    )
+    flm_prime[0 : spin * spin] = 0.0
+    f_prime = ssht.forward_adjoint(
+        flm_prime, order, Spin=spin, Method=method, backend="ducc"
+    )
+
+    assert flm_prime.conj() @ flm == approx(f_prime.flatten().conj() @ f.flatten())
+
+
 def test_rot(complex_coeffs, order):
     gamma = np.pi / 5
     beta = np.pi / 7
diff --git a/typings/pyssht/__init__.pyi b/typings/pyssht/__init__.pyi
@@ -12,6 +12,17 @@ def forward(
 ) -> np.ndarray:
     pass
 
+def forward_adjoint(
+    flm: np.ndarray,
+    L: int,
+    Spin: int = 0,
+    Method: str = "MW",
+    Reality: bool = False,
+    backend: str = "SSHT",
+    **kwargs,
+) -> np.ndarray:
+    pass
+
 def inverse(
     flm: np.ndarray,
     L: int,
@@ -57,3 +68,9 @@ def rotate_flms(
     **kwargs,
 ) -> np.ndarray:
     pass
+
+def sample_length(L: int, Method: str = "MW") -> int:
+    pass
+
+def sample_shape(L: int, Method: str = "MW") -> Tuple[int, int]:
+    pass
