Aquila-Consortium/libsharp2
0
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

step 1

Browse source
This commit is contained in:
Martin Reinecke 2018-10-15 17:26:23 +02:00
parent 593d4eba67
commit 31cbd2acc5
50 changed files with 488 additions and 7483 deletions
Download patch file Download diff file Expand all files Collapse all files

1
python/fake_pyrex/Pyrex/Distutils/__init__.py
View file

@ -1 +0,0 @@
# work around broken setuptools monkey patching

1
python/fake_pyrex/Pyrex/Distutils/build_ext.py
View file

@ -1 +0,0 @@
build_ext = "yes, it's there!"

1
python/fake_pyrex/Pyrex/__init__.py
View file

@ -1 +0,0 @@
# work around broken setuptools monkey patching

2
python/fake_pyrex/README
View file

@ -1,2 +0,0 @@
This directory is here to fool setuptools into building .pyx files
even if Pyrex is not installed. See ../setup.py.

1
python/libsharp/__init__.py
View file

@ -1 +0,0 @@
from .libsharp import *

92
python/libsharp/libsharp.pxd
View file

@ -1,92 +0,0 @@
cdef extern from "sharp.h":
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)
void sharp_legendre_roots(int n, double *x, double *w)
# sharp_lowlevel.h
ctypedef struct sharp_alm_info:
# Maximum \a l index of the array
int lmax
# Number of different \a m values in this object
int nm
# Array with \a nm entries containing the individual m values
int *mval
# Combination of flags from sharp_almflags
int flags
# Array with \a nm entries containing the (hypothetical) indices of
# the coefficients with quantum numbers 0,\a mval[i]
long *mvstart
# Stride between a_lm and a_(l+1),m
long stride
ctypedef struct sharp_geom_info:
pass
void sharp_make_alm_info (int lmax, int mmax, int stride,
ptrdiff_t *mvstart, sharp_alm_info **alm_info)
void sharp_make_geom_info (int nrings, int *nph, ptrdiff_t *ofs,
int *stride, double *phi0, double *theta,
double *wgt, sharp_geom_info **geom_info)
void sharp_destroy_alm_info(sharp_alm_info *info)
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)
ctypedef enum sharp_jobtype:
SHARP_YtW
SHARP_Yt
SHARP_WY
SHARP_Y
ctypedef enum:
SHARP_DP
SHARP_ADD
void sharp_execute(sharp_jobtype type_,
int spin,
void *alm,
void *map,
sharp_geom_info *geom_info,
sharp_alm_info *alm_info,
int ntrans,
int flags,
double *time,
unsigned long long *opcnt) nogil
ctypedef enum:
SHARP_ERROR_NO_MPI
int sharp_execute_mpi_maybe (void *pcomm, sharp_jobtype type, int spin,
void *alm, void *map, sharp_geom_info *geom_info,
sharp_alm_info *alm_info, int ntrans, int flags, double *time,
unsigned long long *opcnt) nogil
void sharp_normalized_associated_legendre_table(int m, int spin, int lmax, int ntheta,
double *theta, int theta_stride, int l_stride, int spin_stride, double *out) nogil
cdef extern from "sharp_geomhelpers.h":
void sharp_make_subset_healpix_geom_info(
int nside, int stride, int nrings,
int *rings, double *weight, sharp_geom_info **geom_info)
void sharp_make_gauss_geom_info(
int nrings, int nphi, double phi0,
int stride_lon, int stride_lat, sharp_geom_info **geom_info)
cdef extern from "sharp_almhelpers.h":
void sharp_make_triangular_alm_info (int lmax, int mmax, int stride,
sharp_alm_info **alm_info)
void sharp_make_rectangular_alm_info (int lmax, int mmax, int stride,
sharp_alm_info **alm_info)
void sharp_make_mmajor_real_packed_alm_info (int lmax, int stride,
int nm, const int *ms, sharp_alm_info **alm_info)

324
python/libsharp/libsharp.pyx
View file

@ -1,324 +0,0 @@
import numpy as np
cimport numpy as np
cimport cython
__all__ = ['legendre_transform', 'legendre_roots', 'sht', 'synthesis', 'adjoint_synthesis',
'analysis', 'adjoint_analysis', 'healpix_grid', 'triangular_order', 'rectangular_order',
'packed_real_order', 'normalized_associated_legendre_table']
def legendre_transform(x, bl, out=None):
if out is None:
out = np.empty_like(x)
if out.shape[0] == 0:
return out
elif 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)
def legendre_roots(n):
x = np.empty(n, np.double)
w = np.empty(n, np.double)
cdef double[::1] x_buf = x, w_buf = w
if not (x_buf.shape[0] == w_buf.shape[0] == n):
raise AssertionError()
if n > 0:
sharp_legendre_roots(n, &x_buf[0], &w_buf[0])
return x, w
JOBTYPE_TO_CONST = {
'Y': SHARP_Y,
'Yt': SHARP_Yt,
'WY': SHARP_WY,
'YtW': SHARP_YtW
}
def sht(jobtype, geom_info ginfo, alm_info ainfo, double[:, :, ::1] input,
int spin=0, comm=None, add=False):
cdef void *comm_ptr
cdef int flags = SHARP_DP | (SHARP_ADD if add else 0)
cdef int r
cdef sharp_jobtype jobtype_i
cdef double[:, :, ::1] output_buf
cdef int ntrans = input.shape[0]
cdef int ntotcomp = ntrans * input.shape[1]
cdef int i, j
if spin == 0 and input.shape[1] != 1:
raise ValueError('For spin == 0, we need input.shape[1] == 1')
elif spin != 0 and input.shape[1] != 2:
raise ValueError('For spin != 0, we need input.shape[1] == 2')
cdef size_t[::1] ptrbuf = np.empty(2 * ntotcomp, dtype=np.uintp)
cdef double **alm_ptrs = <double**>&ptrbuf[0]
cdef double **map_ptrs = <double**>&ptrbuf[ntotcomp]
try:
jobtype_i = JOBTYPE_TO_CONST[jobtype]
except KeyError:
raise ValueError('jobtype must be one of: %s' % ', '.join(sorted(JOBTYPE_TO_CONST.keys())))
if jobtype_i == SHARP_Y or jobtype_i == SHARP_WY:
output = np.empty((input.shape[0], input.shape[1], ginfo.local_size()), dtype=np.float64)
output_buf = output
for i in range(input.shape[0]):
for j in range(input.shape[1]):
alm_ptrs[i * input.shape[1] + j] = &input[i, j, 0]
map_ptrs[i * input.shape[1] + j] = &output_buf[i, j, 0]
else:
output = np.empty((input.shape[0], input.shape[1], ainfo.local_size()), dtype=np.float64)
output_buf = output
for i in range(input.shape[0]):
for j in range(input.shape[1]):
alm_ptrs[i * input.shape[1] + j] = &output_buf[i, j, 0]
map_ptrs[i * input.shape[1] + j] = &input[i, j, 0]
if comm is None:
with nogil:
sharp_execute (
jobtype_i,
geom_info=ginfo.ginfo, alm_info=ainfo.ainfo,
spin=spin, alm=alm_ptrs, map=map_ptrs,
ntrans=ntrans, flags=flags, time=NULL, opcnt=NULL)
else:
from mpi4py import MPI
if not isinstance(comm, MPI.Comm):
raise TypeError('comm must be an mpi4py communicator')
from .libsharp_mpi import _addressof
comm_ptr = <void*><size_t>_addressof(comm)
with nogil:
r = sharp_execute_mpi_maybe (
comm_ptr, jobtype_i,
geom_info=ginfo.ginfo, alm_info=ainfo.ainfo,
spin=spin, alm=alm_ptrs, map=map_ptrs,
ntrans=ntrans, flags=flags, time=NULL, opcnt=NULL)
if r == SHARP_ERROR_NO_MPI:
raise Exception('MPI requested, but not available')
return output
def synthesis(*args, **kw):
return sht('Y', *args, **kw)
def adjoint_synthesis(*args, **kw):
return sht('Yt', *args, **kw)
def analysis(*args, **kw):
return sht('YtW', *args, **kw)
def adjoint_analysis(*args, **kw):
return sht('WY', *args, **kw)
#
# geom_info
#
class NotInitializedError(Exception):
pass
cdef class geom_info:
cdef sharp_geom_info *ginfo
def __cinit__(self, *args, **kw):
self.ginfo = NULL
def local_size(self):
if self.ginfo == NULL:
raise NotInitializedError()
return sharp_map_size(self.ginfo)
def __dealloc__(self):
if self.ginfo != NULL:
sharp_destroy_geom_info(self.ginfo)
self.ginfo = NULL
cdef class healpix_grid(geom_info):
_weight_cache = {} # { (nside, 'T'/'Q'/'U') -> numpy array of ring weights cached from file }
def __init__(self, int nside, stride=1, int[::1] rings=None, double[::1] weights=None):
if weights is not None and weights.shape[0] != 2 * nside:
raise ValueError('weights must have length 2 * nside')
sharp_make_subset_healpix_geom_info(nside, stride,
nrings=4 * nside - 1 if rings is None else rings.shape[0],
rings=NULL if rings is None else &rings[0],
weight=NULL if weights is None else &weights[0],
geom_info=&self.ginfo)
@classmethod
def load_ring_weights(cls, nside, fields):
"""
Loads HEALPix ring weights from file. The environment variable
HEALPIX should be set, and this routine will look in the `data`
subdirectory.
Parameters
----------
nside: int
HEALPix nside parameter
fields: tuple of str
Which weights to extract; pass ('T',) to only get scalar
weights back, or ('T', 'Q', 'U') to get all the weights
Returns
-------
List of NumPy arrays, according to fields parameter.
"""
import os
from astropy.io import fits
data_path = os.path.join(os.environ['HEALPIX'], 'data')
fits_field_names = {
'T': 'TEMPERATURE WEIGHTS',
'Q': 'Q-POLARISATION WEIGHTS',
'U': 'U-POLARISATION WEIGHTS'}
must_load = [field for field in fields if (nside, field) not in cls._weight_cache]
if must_load:
hdulist = fits.open(os.path.join(data_path, 'weight_ring_n%05d.fits' % nside))
try:
for field in must_load:
w = hdulist[1].data.field(fits_field_names[field]).ravel().astype(np.double)
w += 1
cls._weight_cache[nside, field] = w
finally:
hdulist.close()
return [cls._weight_cache[(nside, field)].copy() for field in fields]
#
# alm_info
#
cdef class alm_info:
cdef sharp_alm_info *ainfo
def __cinit__(self, *args, **kw):
self.ainfo = NULL
def local_size(self):
if self.ainfo == NULL:
raise NotInitializedError()
return sharp_alm_count(self.ainfo)
def mval(self):
if self.ainfo == NULL:
raise NotInitializedError()
return np.asarray(<int[:self.ainfo.nm]> self.ainfo.mval)
def mvstart(self):
if self.ainfo == NULL:
raise NotInitializedError()
return np.asarray(<long[:self.ainfo.nm]> self.ainfo.mvstart)
def __dealloc__(self):
if self.ainfo != NULL:
sharp_destroy_alm_info(self.ainfo)
self.ainfo = NULL
@cython.boundscheck(False)
def almxfl(self, np.ndarray[double, ndim=3, mode='c'] alm, np.ndarray[double, ndim=2, mode='c'] fl):
"""Multiply Alm by a Ell based array
Parameters
----------
alm : np.ndarray
input alm, 3 dimensions = (different signal x polarizations x lm-ordering)
fl : np.ndarray
either 1 dimension, e.g. gaussian beam, or 2 dimensions e.g. a polarized beam
Returns
-------
None, it modifies alms in-place
"""
cdef int mvstart = 0
cdef bint has_multiple_beams = alm.shape[2] > 1 and fl.shape[1] > 1
cdef int f, i_m, m, num_ells, i_l, i_signal, i_pol, i_mv
for i_m in range(self.ainfo.nm):
m = self.ainfo.mval[i_m]
f = 1 if (m==0) else 2
num_ells = self.ainfo.lmax + 1 - m
if not has_multiple_beams:
for i_signal in range(alm.shape[0]):
for i_pol in range(alm.shape[1]):
for i_l in range(num_ells):
l = m + i_l
for i_mv in range(mvstart + f*i_l, mvstart + f*i_l +f):
alm[i_signal, i_pol, i_mv] *= fl[l, 0]
else:
for i_signal in range(alm.shape[0]):
for i_pol in range(alm.shape[1]):
for i_l in range(num_ells):
l = m + i_l
for i_mv in range(mvstart + f*i_l, mvstart + f*i_l +f):
alm[i_signal, i_pol, i_mv] *= fl[l, i_pol]
mvstart += f * num_ells
cdef class triangular_order(alm_info):
def __init__(self, int lmax, mmax=None, stride=1):
mmax = mmax if mmax is not None else lmax
sharp_make_triangular_alm_info(lmax, mmax, stride, &self.ainfo)
cdef class rectangular_order(alm_info):
def __init__(self, int lmax, mmax=None, stride=1):
mmax = mmax if mmax is not None else lmax
sharp_make_rectangular_alm_info(lmax, mmax, stride, &self.ainfo)
cdef class packed_real_order(alm_info):
def __init__(self, int lmax, stride=1, int[::1] ms=None):
sharp_make_mmajor_real_packed_alm_info(lmax=lmax, stride=stride,
nm=lmax + 1 if ms is None else ms.shape[0],
ms=NULL if ms is None else &ms[0],
alm_info=&self.ainfo)
#
#
#
@cython.boundscheck(False)
def normalized_associated_legendre_table(int lmax, int m, theta):
cdef double[::1] theta_ = np.ascontiguousarray(theta, dtype=np.double)
out = np.zeros((theta_.shape[0], lmax - m + 1), np.double)
cdef double[:, ::1] out_ = out
if lmax < m:
raise ValueError("lmax < m")
with nogil:
sharp_normalized_associated_legendre_table(m, 0, lmax, theta_.shape[0], &theta_[0], lmax - m + 1, 1, 1, &out_[0,0])
return out

17
python/libsharp/libsharp_mpi.pyx
View file

@ -1,17 +0,0 @@
cdef extern from "mpi.h":
ctypedef void *MPI_Comm
cdef extern from "Python.h":
object PyLong_FromVoidPtr(void*)
cdef extern:
ctypedef class mpi4py.MPI.Comm [object PyMPICommObject]:
cdef MPI_Comm ob_mpi
cdef unsigned flags
# For compatibility with mpi4py <= 1.3.1
# Newer versions could use the MPI._addressof function
def _addressof(Comm comm):
cdef void *ptr = NULL
ptr = <void*>&comm.ob_mpi
return PyLong_FromVoidPtr(ptr)

1
python/libsharp/tests/__init__.py
View file

@ -1 +0,0 @@
# empty

58
python/libsharp/tests/test_legendre.py
View file

@ -1,58 +0,0 @@
import numpy as np
from scipy.special import legendre
from scipy.special import p_roots
import libsharp
from numpy.testing import assert_allclose
def check_legendre_transform(lmax, ntheta):
l = np.arange(lmax + 1)
if lmax >= 1:
sigma = -np.log(1e-3) / lmax / (lmax + 1)
bl = np.exp(-sigma*l*(l+1))
bl *= (2 * l + 1)
else:
bl = np.asarray([1], dtype=np.double)
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_allclose(y, y0, rtol=1e-12, atol=1e-12)
# single-precision
y32 = libsharp.legendre_transform(x.astype(np.float32), bl)
assert_allclose(y, y0, rtol=1e-5, atol=1e-5)
def test_legendre_transform():
nthetas_to_try = [0, 9, 17, 19] + list(np.random.randint(500, size=20))
for ntheta in nthetas_to_try:
for lmax in [0, 1, 2, 3, 20] + list(np.random.randint(50, size=4)):
yield check_legendre_transform, lmax, ntheta
def check_legendre_roots(n):
xs, ws = ([], []) if n == 0 else p_roots(n) # from SciPy
xl, wl = libsharp.legendre_roots(n)
assert_allclose(xs, xl, rtol=1e-14, atol=1e-14)
assert_allclose(ws, wl, rtol=1e-14, atol=1e-14)
def test_legendre_roots():
"""
Test the Legendre root-finding algorithm from libsharp by comparing it with
the SciPy version.
"""
yield check_legendre_roots, 0
yield check_legendre_roots, 1
yield check_legendre_roots, 32
yield check_legendre_roots, 33

36
python/libsharp/tests/test_legendre_table.py
View file

@ -1,36 +0,0 @@
from __future__ import print_function
import numpy as np
from numpy.testing import assert_almost_equal
from nose.tools import eq_, ok_
from libsharp import normalized_associated_legendre_table
from scipy.special import sph_harm, p_roots
def test_compare_legendre_table_with_scipy():
def test(theta, m, lmax):
Plm = normalized_associated_legendre_table(lmax, m, theta)
Plm_p = sph_harm(m, np.arange(m, lmax + 1), 0, theta)[None, :]
if not np.allclose(Plm_p, Plm):
print(Plm_p)
print(Plm)
return ok_, np.allclose(Plm_p, Plm)
yield test(np.pi/2, 0, 10)
yield test(np.pi/4, 0, 10)
yield test(3 * np.pi/4, 0, 10)
yield test(np.pi/4, 1, 4)
yield test(np.pi/4, 2, 4)
yield test(np.pi/4, 50, 50)
yield test(np.pi/2, 49, 50)
def test_legendre_table_wrapper_logic():
# tests the SSE 2 logic in the high-level wrapper by using an odd number of thetas
theta = np.asarray([np.pi/2, np.pi/4, 3 * np.pi / 4])
m = 3
lmax = 10
Plm = normalized_associated_legendre_table(lmax, m, theta)
assert np.allclose(Plm[1, :], normalized_associated_legendre_table(lmax, m, np.pi/4)[0, :])
assert np.allclose(Plm[2, :], normalized_associated_legendre_table(lmax, m, 3 * np.pi/4)[0, :])

32
python/libsharp/tests/test_sht.py
View file

@ -1,32 +0,0 @@
import numpy as np
from numpy.testing import assert_allclose
import libsharp
from mpi4py import MPI
def test_basic():
lmax = 10
nside = 8
rank = MPI.COMM_WORLD.Get_rank()
ms = np.arange(rank, lmax + 1, MPI.COMM_WORLD.Get_size(), dtype=np.int32)
order = libsharp.packed_real_order(lmax, ms=ms)
grid = libsharp.healpix_grid(nside)
alm = np.zeros(order.local_size())
if rank == 0:
alm[0] = 1
elif rank == 1:
alm[0] = 1
map = libsharp.synthesis(grid, order, np.repeat(alm[None, None, :], 3, 0), comm=MPI.COMM_WORLD)
assert np.all(map[2, :] == map[1, :]) and np.all(map[1, :] == map[0, :])
map = map[0, 0, :]
print(rank, "shape", map.shape)
print(rank, "mean", map.mean())
if __name__=="__main__":
test_basic()

137
python/libsharp/tests/test_smoothing_noise_pol_mpi.py
View file

@ -1,137 +0,0 @@
# This test needs to be run with:
# mpirun -np X python test_smoothing_noise_pol_mpi.py
from mpi4py import MPI
import numpy as np
import healpy as hp
import libsharp
mpi = True
rank = MPI.COMM_WORLD.Get_rank()
nside = 256
npix = hp.nside2npix(nside)
np.random.seed(100)
input_map = np.random.normal(size=(3, npix))
fwhm_deg = 10
lmax = 512
nrings = 4 * nside - 1 # four missing pixels
if rank == 0:
print("total rings", nrings)
n_mpi_processes = MPI.COMM_WORLD.Get_size()
rings_per_process = nrings // n_mpi_processes + 1
# ring indices are 1-based
ring_indices_emisphere = np.arange(2*nside, dtype=np.int32) + 1
local_ring_indices = ring_indices_emisphere[rank::n_mpi_processes]
# to improve performance, simmetric rings north/south need to be in the same rank
# therefore we use symmetry to create the full ring indexing
if local_ring_indices[-1] == 2 * nside:
# has equator ring
local_ring_indices = np.concatenate(
[local_ring_indices[:-1],
nrings - local_ring_indices[::-1] + 1]
)
else:
# does not have equator ring
local_ring_indices = np.concatenate(
[local_ring_indices,
nrings - local_ring_indices[::-1] + 1]
)
print("rank", rank, "n_rings", len(local_ring_indices))
if not mpi:
local_ring_indices = None
grid = libsharp.healpix_grid(nside, rings=local_ring_indices)
# returns start index of the ring and number of pixels
startpix, ringpix, _, _, _ = hp.ringinfo(nside, local_ring_indices.astype(np.int64))
local_npix = grid.local_size()
def expand_pix(startpix, ringpix, local_npix):
"""Turn first pixel index and number of pixel in full array of pixels
to be optimized with cython or numba
"""
local_pix = np.empty(local_npix, dtype=np.int64)
i = 0
for start, num in zip(startpix, ringpix):
local_pix[i:i+num] = np.arange(start, start+num)
i += num
return local_pix
local_pix = expand_pix(startpix, ringpix, local_npix)
local_map = input_map[:, local_pix]
local_hitmap = np.zeros(npix)
local_hitmap[local_pix] = 1
hp.write_map("hitmap_{}.fits".format(rank), local_hitmap, overwrite=True)
print("rank", rank, "npix", npix, "local_npix", local_npix, "local_map len", len(local_map), "unique pix", len(np.unique(local_pix)))
local_m_indices = np.arange(rank, lmax + 1, MPI.COMM_WORLD.Get_size(), dtype=np.int32)
if not mpi:
local_m_indices = None
order = libsharp.packed_real_order(lmax, ms=local_m_indices)
local_nl = order.local_size()
print("rank", rank, "local_nl", local_nl, "mval", order.mval())
mpi_comm = MPI.COMM_WORLD if mpi else None
# map2alm
# maps in libsharp are 3D, 2nd dimension is IQU, 3rd is pixel
alm_sharp_I = libsharp.analysis(grid, order,
np.ascontiguousarray(local_map[0].reshape((1, 1, -1))),
spin=0, comm=mpi_comm)
alm_sharp_P = libsharp.analysis(grid, order,
np.ascontiguousarray(local_map[1:].reshape((1, 2, -1))),
spin=2, comm=mpi_comm)
beam = hp.gauss_beam(fwhm=np.radians(fwhm_deg), lmax=lmax, pol=True)
print("Smooth")
# smooth in place (zonca implemented this function)
order.almxfl(alm_sharp_I, np.ascontiguousarray(beam[:, 0:1]))
order.almxfl(alm_sharp_P, np.ascontiguousarray(beam[:, (1, 2)]))
# alm2map
new_local_map_I = libsharp.synthesis(grid, order, alm_sharp_I, spin=0, comm=mpi_comm)
new_local_map_P = libsharp.synthesis(grid, order, alm_sharp_P, spin=2, comm=mpi_comm)
# Transfer map to first process for writing
local_full_map = np.zeros(input_map.shape, dtype=np.float64)
local_full_map[0, local_pix] = new_local_map_I
local_full_map[1:, local_pix] = new_local_map_P
output_map = np.zeros(input_map.shape, dtype=np.float64) if rank == 0 else None
mpi_comm.Reduce(local_full_map, output_map, root=0, op=MPI.SUM)
if rank == 0:
# hp.write_map("sharp_smoothed_map.fits", output_map, overwrite=True)
# hp_smoothed = hp.alm2map(hp.map2alm(input_map, lmax=lmax), nside=nside) # transform only
hp_smoothed = hp.smoothing(input_map, fwhm=np.radians(fwhm_deg), lmax=lmax)
std_diff = (hp_smoothed-output_map).std()
print("Std of difference between libsharp and healpy", std_diff)
# hp.write_map(
# "healpy_smoothed_map.fits",
# hp_smoothed,
# overwrite=True
# )
assert std_diff < 1e-5

83
python/setup.py
View file

@ -1,83 +0,0 @@
#! /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.Build import cythonize
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'],
ext_modules = cythonize([
Extension("libsharp.libsharp",
["libsharp/libsharp.pyx"],
libraries=["sharp", "fftpack", "c_utils"],
include_dirs=[libsharp_include, np.get_include()],
library_dirs=[libsharp_lib],
extra_link_args=["-fopenmp"],
),
Extension("libsharp.libsharp_mpi",
["libsharp/libsharp_mpi.pyx"],
libraries=["sharp", "fftpack", "c_utils"],
include_dirs=[libsharp_include, np.get_include()],
library_dirs=[libsharp_lib],
extra_link_args=["-fopenmp"],
),
]),
)