diff --git a/.gitignore b/.gitignore index cd0bf21..b555caa 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,9 @@ *.o +*.so #* *~ +*.pyc +*.pyo /auto /autom4te.cache @@ -9,3 +12,5 @@ /config/config.auto /configure /sharp_oracle.inc + +/python/libsharp/libsharp.c diff --git a/Makefile b/Makefile index dd9931e..f44ecfd 100644 --- a/Makefile +++ b/Makefile @@ -61,3 +61,10 @@ perftest: compile_all genclean: rm libsharp/sharp_legendre.c || exit 0 + +pytest: + rm python/libsharp/libsharp.so || exit 0 + cd python && LIBSHARP_INCLUDE=$(INCDIR) LIBSHARP_LIB=$(LIBDIR) python setup.py build_ext --inplace + cd python && nosetests libsharp + + diff --git a/python/fake_pyrex/Pyrex/Distutils/__init__.py b/python/fake_pyrex/Pyrex/Distutils/__init__.py new file mode 100644 index 0000000..51c8e16 --- /dev/null +++ b/python/fake_pyrex/Pyrex/Distutils/__init__.py @@ -0,0 +1 @@ +# work around broken setuptools monkey patching diff --git a/python/fake_pyrex/Pyrex/Distutils/build_ext.py b/python/fake_pyrex/Pyrex/Distutils/build_ext.py new file mode 100644 index 0000000..4f846f6 --- /dev/null +++ b/python/fake_pyrex/Pyrex/Distutils/build_ext.py @@ -0,0 +1 @@ +build_ext = "yes, it's there!" diff --git a/python/fake_pyrex/Pyrex/__init__.py b/python/fake_pyrex/Pyrex/__init__.py new file mode 100644 index 0000000..51c8e16 --- /dev/null +++ b/python/fake_pyrex/Pyrex/__init__.py @@ -0,0 +1 @@ +# work around broken setuptools monkey patching diff --git a/python/fake_pyrex/README b/python/fake_pyrex/README new file mode 100644 index 0000000..cf3f3ff --- /dev/null +++ b/python/fake_pyrex/README @@ -0,0 +1,2 @@ +This directory is here to fool setuptools into building .pyx files +even if Pyrex is not installed. See ../setup.py. \ No newline at end of file diff --git a/python/libsharp/__init__.py b/python/libsharp/__init__.py new file mode 100644 index 0000000..dd0fa41 --- /dev/null +++ b/python/libsharp/__init__.py @@ -0,0 +1 @@ +from .libsharp import * diff --git a/python/libsharp/libsharp.pyx b/python/libsharp/libsharp.pyx new file mode 100644 index 0000000..0abfe69 --- /dev/null +++ b/python/libsharp/libsharp.pyx @@ -0,0 +1,41 @@ +import numpy as np + +cdef extern from "sharp.h": + ctypedef long ptrdiff_t + + void sharp_legendre_transform_s(float *bl, float *recfac, ptrdiff_t lmax, float *x, + float *out, ptrdiff_t nx) + void sharp_legendre_transform(double *bl, double *recfac, ptrdiff_t lmax, double *x, + double *out, ptrdiff_t nx) + void sharp_legendre_transform_recfac(double *r, ptrdiff_t lmax) + void sharp_legendre_transform_recfac_s(float *r, ptrdiff_t lmax) + + +def legendre_transform(x, bl, out=None): + if out is None: + out = np.empty_like(x) + if x.dtype == np.float64: + if bl.dtype != np.float64: + bl = bl.astype(np.float64) + return _legendre_transform(x, bl, out=out) + elif x.dtype == np.float32: + if bl.dtype != np.float32: + bl = bl.astype(np.float32) + return _legendre_transform_s(x, bl, out=out) + else: + raise ValueError("unsupported dtype") + + +def _legendre_transform(double[::1] x, double[::1] bl, double[::1] out): + if out.shape[0] != x.shape[0]: + raise ValueError('x and out must have same shape') + sharp_legendre_transform(&bl[0], NULL, bl.shape[0] - 1, &x[0], &out[0], x.shape[0]) + return np.asarray(out) + + +def _legendre_transform_s(float[::1] x, float[::1] bl, float[::1] out): + if out.shape[0] != x.shape[0]: + raise ValueError('x and out must have same shape') + sharp_legendre_transform_s(&bl[0], NULL, bl.shape[0] - 1, &x[0], &out[0], x.shape[0]) + return np.asarray(out) + diff --git a/python/libsharp/tests/__init__.py b/python/libsharp/tests/__init__.py new file mode 100644 index 0000000..1bb8bf6 --- /dev/null +++ b/python/libsharp/tests/__init__.py @@ -0,0 +1 @@ +# empty diff --git a/python/libsharp/tests/test_legendre.py b/python/libsharp/tests/test_legendre.py new file mode 100644 index 0000000..2e96822 --- /dev/null +++ b/python/libsharp/tests/test_legendre.py @@ -0,0 +1,29 @@ +import numpy as np +from scipy.special import legendre +import libsharp + + +def test_legendre_transform(): + lmax = 20 + ntheta = 19 + + l = np.arange(lmax + 1) + bl = np.exp(-l*(l+1)) + bl *= (2 * l + 1) + + theta = np.linspace(0, np.pi, ntheta, endpoint=True) + x = np.cos(theta) + + # Compute truth using scipy.special.legendre + P = np.zeros((ntheta, lmax + 1)) + for l in range(lmax + 1): + P[:, l] = legendre(l)(x) + y0 = np.dot(P, bl) + + # double-precision + y = libsharp.legendre_transform(x, bl) + assert np.max(np.abs(y - y) / np.abs(y)) < 1e-12 + + # single-precision + y32 = libsharp.legendre_transform(x.astype(np.float32), bl) + assert np.max(np.abs(y32 - y) / np.abs(y32)) < 1e-4 diff --git a/python/setup.py b/python/setup.py new file mode 100644 index 0000000..1500f05 --- /dev/null +++ b/python/setup.py @@ -0,0 +1,76 @@ +#! /usr/bin/env python + +descr = """Spherical Harmionic transforms package + +Python API for the libsharp spherical harmonic transforms library +""" + +import os +import sys + +DISTNAME = 'libsharp' +DESCRIPTION = 'libsharp library for fast Spherical Harmonic Transforms' +LONG_DESCRIPTION = descr +MAINTAINER = 'Dag Sverre Seljebotn', +MAINTAINER_EMAIL = 'd.s.seljebotn@astro.uio.no', +URL = 'http://sourceforge.net/projects/libsharp/' +LICENSE = 'GPL' +DOWNLOAD_URL = "http://sourceforge.net/projects/libsharp/" +VERSION = '0.1' + +# Add our fake Pyrex at the end of the Python search path +# in order to fool setuptools into allowing compilation of +# pyx files to C files. Importing Cython.Distutils then +# makes Cython the tool of choice for this rather than +# (the possibly nonexisting) Pyrex. +project_path = os.path.split(__file__)[0] +sys.path.append(os.path.join(project_path, 'fake_pyrex')) + +from setuptools import setup, find_packages, Extension +from Cython.Distutils import build_ext +import numpy as np + +libsharp = os.environ.get('LIBSHARP', None) +libsharp_include = os.environ.get('LIBSHARP_INCLUDE', libsharp and os.path.join(libsharp, 'include')) +libsharp_lib = os.environ.get('LIBSHARP_LIB', libsharp and os.path.join(libsharp, 'lib')) + +if libsharp_include is None or libsharp_lib is None: + sys.stderr.write('Please set LIBSHARP environment variable to the install directly of libsharp, ' + 'this script will refer to the lib and include sub-directories. Alternatively ' + 'set LIBSHARP_INCLUDE and LIBSHARP_LIB\n') + sys.exit(1) + +if __name__ == "__main__": + setup(install_requires = ['numpy'], + packages = find_packages(), + test_suite="nose.collector", + # Well, technically zipping the package will work, but since it's + # all compiled code it'll just get unzipped again at runtime, which + # is pointless: + zip_safe = False, + name = DISTNAME, + version = VERSION, + maintainer = MAINTAINER, + maintainer_email = MAINTAINER_EMAIL, + description = DESCRIPTION, + license = LICENSE, + url = URL, + download_url = DOWNLOAD_URL, + long_description = LONG_DESCRIPTION, + classifiers = + [ 'Development Status :: 3 - Alpha', + 'Environment :: Console', + 'Intended Audience :: Developers', + 'Intended Audience :: Science/Research', + 'License :: OSI Approved :: GNU General Public License (GPL)', + 'Topic :: Scientific/Engineering'], + cmdclass = {"build_ext": build_ext}, + ext_modules = [ + Extension("libsharp.libsharp", + ["libsharp/libsharp.pyx"], + libraries=["sharp", "fftpack", "c_utils"], + include_dirs=[libsharp_include], + library_dirs=[libsharp_lib] + ), + ], + )