From a1306279bd6184e509a670c71fcb1a3f2275a593 Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Sat, 22 Apr 2017 10:35:55 +0200
Subject: [PATCH] use iterative map analysis for Healpix
---
.../transformations/hplmtransformation.py | 7 ++-----
.../transformations/lmhptransformation.py | 9 +++------
test/test_misc.py | 13 +++++++++++++
3 files changed, 18 insertions(+), 11 deletions(-)
diff --git a/nifty/operators/fft_operator/transformations/hplmtransformation.py b/nifty/operators/fft_operator/transformations/hplmtransformation.py
index 8fe7031be..f32e0db86 100644
--- a/nifty/operators/fft_operator/transformations/hplmtransformation.py
+++ b/nifty/operators/fft_operator/transformations/hplmtransformation.py
@@ -88,11 +88,8 @@ class HPLMTransformation(SlicingTransformation):
mmax = lmax
nside= self.domain.nside
- sjob = pyHealpix.sharpjob_d()
- sjob.set_Healpix_geometry(nside)
- sjob.set_triangular_alm_info(lmax, mmax)
if issubclass(inp.dtype.type, np.complexfloating):
- [resultReal, resultImag] = [sjob.map2alm(x)
+ [resultReal, resultImag] = [pyHealpix.map2alm_iter(x,lmax,mmax,3)
for x in (inp.real, inp.imag)]
[resultReal,
@@ -102,7 +99,7 @@ class HPLMTransformation(SlicingTransformation):
result = self._combine_complex_result(resultReal, resultImag)
else:
- result = sjob.map2alm(inp)
+ result = pyHealpix.map2alm_iter(inp,lmax,mmax,3)
result = lm_transformation_factory.buildIdx(result, lmax=lmax)
return result
diff --git a/nifty/operators/fft_operator/transformations/lmhptransformation.py b/nifty/operators/fft_operator/transformations/lmhptransformation.py
index d17d52293..0fc9e425b 100644
--- a/nifty/operators/fft_operator/transformations/lmhptransformation.py
+++ b/nifty/operators/fft_operator/transformations/lmhptransformation.py
@@ -64,7 +64,7 @@ class LMHPTransformation(SlicingTransformation):
if not isinstance(domain, LMSpace):
raise TypeError("domain needs to be a LMSpace.")
- nside = np.max(domain.lmax//2, 1)
+ nside = max((domain.lmax+1) // 2, 1)
result = HPSpace(nside=nside, dtype=domain.dtype)
return result
@@ -89,21 +89,18 @@ class LMHPTransformation(SlicingTransformation):
lmax = self.domain.lmax
mmax = lmax
- sjob = pyHealpix.sharpjob_d()
- sjob.set_Healpix_geometry(nside)
- sjob.set_triangular_alm_info(lmax, mmax)
if issubclass(inp.dtype.type, np.complexfloating):
[resultReal,
resultImag] = [lm_transformation_factory.buildLm(x, lmax=lmax)
for x in (inp.real, inp.imag)]
- [resultReal, resultImag] = [sjob.alm2map(x)
+ [resultReal, resultImag] = [pyHealpix.alm2map(x,lmax,mmax,nside)
for x in [resultReal, resultImag]]
result = self._combine_complex_result(resultReal, resultImag)
else:
result = lm_transformation_factory.buildLm(inp, lmax=lmax)
- result = sjob.alm2map(result)
+ result = pyHealpix.alm2map(result,lmax,mmax,nside)
return result
diff --git a/test/test_misc.py b/test/test_misc.py
index 930be0aba..b8a8f8592 100644
--- a/test/test_misc.py
+++ b/test/test_misc.py
@@ -31,6 +31,7 @@ from nifty import Field,\
FieldArray, \
RGRGTransformation, \
LMGLTransformation, \
+ LMHPTransformation, \
FFTOperator
from nose.plugins.skip import SkipTest
@@ -103,3 +104,15 @@ class Misc_Tests(unittest.TestCase):
inp = Field.from_random(domain=a,random_type='normal',std=7,mean=3,dtype=tp)
out = fft.inverse_times(fft.times(inp))
assert_allclose(inp.val, out.val)
+
+ def test_sht2(self):
+ if 'pyHealpix' not in di:
+ raise SkipTest
+ for lm in [128,256]:
+ for tp in [np.float64,np.complex128,np.float32,np.complex64]:
+ a = LMSpace(lmax=lm)
+ b = LMHPTransformation.get_codomain(a)
+ fft = FFTOperator(domain=a, target=b, domain_dtype=tp, target_dtype=tp)
+ inp = Field.from_random(domain=a,random_type='normal',std=1,mean=0,dtype=tp)
+ out = fft.inverse_times(fft.times(inp))
+ assert_allclose(inp.val, out.val,rtol=1e-3,atol=1e-3)
--
GitLab