First draft of SHTs in Python wrapper

[dagss]

This is a first draft around a Python wrapper that supports MPI. There is a lot of lacking features.

The current testcase splits the alms, but does synthesis to all the nodes; this should be changed so that every node has one ring each, then there is MPI_gather in the end to gather the map for checking at the root node. Also obviously the testcase should compare to healpy or similar, rather than do a visual plot that needs human intervention.

So far only the layout that as built-in (i.e., takes an ms parameter where you can tell it how to split) is the realpacked format described here:

http://commander.readthedocs.org/en/latest/library_guide/sphere.html

I couldn't figure out how to (quickly...) link to MPI when using Python using Python distutils, so the wrapper extension .so is for now built like this:

SHARP_TARGET=your_mpi_config make python/libsharp/libsharp.so

after which you can do

cd python
mpiexec -np 4 nosetests libsharp/tests/test_sht.py

to run the incomplete, manual, test. Since libsharp needs to be built anyway, and thus the make system involved, I think it is probably smarter to keep building the Python extension this way -- we can emit a Wheel package in the build directory which can then be installed by setuptools/pip.

[dagss]

The HEALPix ring weight loading is not tested at all...

[pchanial]

I could compile the extension successfully after some hacking in the Makefile and fixing the typo:

diff --git a/python/libsharp/libsharp.pxd b/python/libsharp/libsharp.pxd
index fc27f19..9209b2b 100644
--- a/python/libsharp/libsharp.pxd
+++ b/python/libsharp/libsharp.pxd
@@ -27,7 +27,7 @@ cdef extern from "sharp.h":
     void sharp_destroy_geom_info(sharp_geom_info *info)

     ptrdiff_t sharp_map_size(sharp_geom_info *info)
-    ptrdiff_t sharp_alm_count(sharp_alm_info *self);
+    ptrdiff_t sharp_alm_count(sharp_alm_info *self)

But then, it looks like you've modified mpi4py to expose the pointer to the C communicator, by adding an function _addressof. Haven't you? I've tried to use the Fortran handle from py2f but so far I've not been successful...

[dagss]

On 05/18/2015 03:54 PM, Pierre Chanial wrote:

I could compile the extension successfully after some hacking in the Makefile and fixing the typo:

|diff --git a/python/libsharp/libsharp.pxd b/python/libsharp/libsharp.pxd index fc27f19..9209b2b 100644 --- a/python/libsharp/libsharp.pxd +++ b/python/libsharp/libsharp.pxd @@ -27,7 +27,7 @@ cdef extern from "sharp.h": void sharp_destroy_geom_info(sharp_geom_info *info)

  ptrdiff_t sharp_map_size(sharp_geom_info *info)

• ptrdiff_t sharp_alm_count(sharp_alm_info *self);
• ptrdiff_t sharp_alm_count(sharp_alm_info *self)

But then, it looks like you've modified mpi4py to expose the pointer to the C communicator, by adding an function _addressof. I've tried to use the Fortran handle from py2f but so far I've not been successful...

I think you just have to upgrade your mpi4py, it's definitely official, I was tipsed about the _addressof my lead mpi4py author on mpi4py mailing list.

I guess different Cython versions might account for the typo too, I'm not sure..

Dag Sverre

[pchanial]

I got it running, great! Good job! I'll spend more time testing it.

[dagss]

@pchanial please check the commit above where I added support for multiple SHTs and spin-weighted transforms. Input/output now have 3 dimensions; first dimension is number of independent transforms, second dimension is 1 for spin=0 and 2 for spin=2, and third is pixels/alm.

Note that I didn't pull in the stuff from your branch yet, I made a comment.

[dagss]

PS not sure if you know already, so just in case, to do polarized transforms you need to do one on one map spin=0 and one on two maps at spin=2. Adding a wrapper sht_pol for this would be a welcome contribution. Look in sharp_healpix_f_inc.c file in the HEALPix sources for example, I don't really know the details myself so I don't trust myself to get it perfectly right.

[dagss]

Note to self: The NumPy include path should be added to the Makefile.

[dagss]

Thanks a lot @pchanial !