Commit 900e4084 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

introduce variable mmax to LMSpace

parent f1eeccb4
Pipeline #21144 passed with stage
in 4 minutes and 8 seconds
......@@ -70,48 +70,15 @@ class RGRGTransformation(Transformation):
return Tval
class SlicingTransformation(Transformation):
def transform(self, x):
p2h = x.domain == self.pdom
if p2h:
res = Field(self.hdom, dtype=x.dtype)
for slice in utilities.get_slice_list(x.shape,
x.domain.axes[self.space]):
res.val[slice] = self._slice_p2h(x.val[slice])
else:
res = Field(self.pdom, dtype=x.dtype)
for slice in utilities.get_slice_list(x.shape,
x.domain.axes[self.space]):
res.val[slice] = self._slice_h2p(x.val[slice])
return res
def buildLm(nr, lmax):
res = np.empty((len(nr)+lmax+1)//2, dtype=(nr[0]*1j).dtype)
res[0:lmax+1] = nr[0:lmax+1]
res[lmax+1:] = np.sqrt(0.5)*(nr[lmax+1::2] + 1j*nr[lmax+2::2])
return res
def buildIdx(nr, lmax):
res = np.empty((lmax+1)*(lmax+1), dtype=nr[0].real.dtype)
res[0:lmax+1] = nr[0:lmax+1].real
res[lmax+1::2] = np.sqrt(2)*nr[lmax+1:].real
res[lmax+2::2] = np.sqrt(2)*nr[lmax+1:].imag
return res
class SphericalTransformation(SlicingTransformation):
class SphericalTransformation(Transformation):
def __init__(self, pdom, hdom, space):
super(SphericalTransformation, self).__init__(pdom, hdom, space)
from pyHealpix import sharpjob_d
self.lmax = self.hdom[self.space].lmax
self.mmax = self.hdom[self.space].mmax
self.sjob = sharpjob_d()
self.sjob.set_triangular_alm_info(self.hdom[self.space].lmax,
self.hdom[self.space].mmax)
self.sjob.set_triangular_alm_info(self.lmax, self.mmax)
if isinstance(self.pdom[self.space], GLSpace):
self.sjob.set_Gauss_geometry(self.pdom[self.space].nlat,
self.pdom[self.space].nlon)
......@@ -124,8 +91,32 @@ class SphericalTransformation(SlicingTransformation):
def _slice_p2h(self, inp):
rr = self.sjob.map2alm(inp)
return buildIdx(rr, lmax=self.lmax)
assert len(rr) == ((self.mmax+1)*(self.mmax+2))//2 + \
(self.mmax+1)*(self.lmax-self.mmax)
res = np.empty(2*len(rr)-self.lmax-1, dtype=rr[0].real.dtype)
res[0:self.lmax+1] = rr[0:self.lmax+1].real
res[self.lmax+1::2] = np.sqrt(2)*rr[self.lmax+1:].real
res[self.lmax+2::2] = np.sqrt(2)*rr[self.lmax+1:].imag
return res
def _slice_h2p(self, inp):
result = buildLm(inp, lmax=self.lmax)
return self.sjob.alm2map(result)
res = np.empty((len(inp)+self.lmax+1)//2, dtype=(inp[0]*1j).dtype)
assert len(res) == ((self.mmax+1)*(self.mmax+2))//2 + \
(self.mmax+1)*(self.lmax-self.mmax)
res[0:self.lmax+1] = inp[0:self.lmax+1]
res[self.lmax+1:] = np.sqrt(0.5)*(inp[self.lmax+1::2] +
1j*inp[self.lmax+2::2])
return self.sjob.alm2map(res)
def transform(self, x):
p2h = x.domain == self.pdom
axes = x.domain.axes[self.space]
if p2h:
res = Field(self.hdom if p2h else self.pdom, dtype=x.dtype)
for slice in utilities.get_slice_list(x.shape, axes):
res.val[slice] = self._slice_p2h(x.val[slice])
else:
res = Field(self.pdom, dtype=x.dtype)
for slice in utilities.get_slice_list(x.shape, axes):
res.val[slice] = self._slice_h2p(x.val[slice])
return res
......@@ -112,7 +112,7 @@ class GLSpace(Space):
def get_default_codomain(self):
from .. import LMSpace
return LMSpace(lmax=self.nlat-1)
return LMSpace(lmax=self._nlat-1, mmax=self._nlon//2)
def check_codomain(self, codomain):
from .. import LMSpace
......
......@@ -32,6 +32,13 @@ class LMSpace(Space):
lmax : int
The maximum :math:`l` value of any spherical harmonics
:math:`Y_{lm}` that is represented in this Space.
Must be >=0.
mmax : int *optional*
The maximum :math:`m` value of any spherical harmonics
:math:`Y_{lm}` that is represented in this Space.
If not supplied, it is set to lmax.
Must be >=0 and <=lmax.
See Also
--------
......@@ -39,15 +46,6 @@ class LMSpace(Space):
GLSpace : A class for the Gauss-Legendre discretization of the
sphere [#]_.
Raises
------
ValueError
If given lmax is negative.
Notes
-----
This implementation implicitly sets the mmax parameter to lmax.
References
----------
.. [#] K.M. Gorski et al., 2005, "HEALPix: A Framework for
......@@ -58,13 +56,20 @@ class LMSpace(Space):
`arXiv:1303.4945 <http://www.arxiv.org/abs/1303.4945>`_
"""
def __init__(self, lmax):
def __init__(self, lmax, mmax=None):
super(LMSpace, self).__init__()
self._needed_for_hash += ["_lmax"]
self._lmax = self._parse_lmax(lmax)
self._needed_for_hash += ["_lmax", "_mmax"]
self._lmax = np.int(lmax)
if self._lmax < 0:
raise ValueError("lmax must be >=0.")
if mmax is None:
mmax = self._lmax
self._mmax = np.int(mmax)
if self._mmax < 0 or self._mmax > self._lmax:
raise ValueError("mmax must be >=0 and <=lmax.")
def __repr__(self):
return ("LMSpace(lmax=%r)" % self.lmax)
return ("LMSpace(lmax=%r, mmax=%r)" % (self.lmax, self.mmax))
@property
def harmonic(self):
......@@ -76,28 +81,27 @@ class LMSpace(Space):
@property
def dim(self):
l = self.lmax
l = self._lmax
m = self._mmax
# the LMSpace consists of the full triangle (including -m's!),
# minus two little triangles if mmax < lmax
# dim = (((2*(l+1)-1)+1)**2/4 - 2 * (l-m)(l-m+1)/2
# dim = np.int((l+1)**2 - (l-m)*(l-m+1.))
# We fix l == m
return np.int((l+1)*(l+1))
return (l+1)**2 - (l-m)*(l-m+1)
def scalar_dvol(self):
return 1.
def get_k_length_array(self):
lmax = self.lmax
lmax = self._lmax
mmax = self._mmax
ldist = np.empty((self.dim,), dtype=np.float64)
ldist[0:lmax+1] = np.arange(lmax+1, dtype=np.float64)
tmp = np.empty((2*lmax+2), dtype=np.float64)
tmp[0::2] = np.arange(lmax+1)
tmp[1::2] = np.arange(lmax+1)
idx = lmax+1
for l in range(1, lmax+1):
ldist[idx:idx+2*(lmax+1-l)] = tmp[2*l:]
idx += 2*(lmax+1-l)
for m in range(1, mmax+1):
ldist[idx:idx+2*(lmax+1-m)] = tmp[2*m:]
idx += 2*(lmax+1-m)
return Field((self,), ldist)
def get_unique_k_lengths(self):
......@@ -105,6 +109,9 @@ class LMSpace(Space):
@staticmethod
def _kernel(x, sigma):
# cf. "All-sky convolution for polarimetry experiments"
# by Challinor et al.
# http://arxiv.org/abs/astro-ph/0008228
res = x+1.
res *= x
res *= -0.5*sigma*sigma
......@@ -112,9 +119,6 @@ class LMSpace(Space):
return res
def get_fft_smoothing_kernel_function(self, sigma):
# cf. "All-sky convolution for polarimetry experiments"
# by Challinor et al.
# http://arxiv.org/abs/astro-ph/0008228
return lambda x: self._kernel(x, sigma)
@property
......@@ -128,15 +132,8 @@ class LMSpace(Space):
def mmax(self):
""" Returns the maximum :math:`m` value of any spherical harmonic
:math:`Y_{lm}` that is represented in this Space.
Currently this is identical to lmax.
"""
return self._lmax
def _parse_lmax(self, lmax):
lmax = np.int(lmax)
if lmax < 0:
raise ValueError("lmax must be >=0.")
return lmax
return self._mmax
def get_default_codomain(self):
from .. import GLSpace
......@@ -145,4 +142,4 @@ class LMSpace(Space):
def check_codomain(self, codomain):
from .. import GLSpace, HPSpace
if not isinstance(codomain, (GLSpace, HPSpace)):
raise TypeError("codomain must be a GLSpace.")
raise TypeError("codomain must be a GLSpace or HPSpace.")
......@@ -27,21 +27,21 @@ from test.common import expand
# [lmax, expected]
CONSTRUCTOR_CONFIGS = [
[5, {
[5, None, {
'lmax': 5,
'mmax': 5,
'shape': (36,),
'harmonic': True,
'dim': 36,
}],
[7, {
[7, 4, {
'lmax': 7,
'mmax': 7,
'shape': (64,),
'mmax': 4,
'shape': (52,),
'harmonic': True,
'dim': 64,
'dim': 52,
}],
[-1, {
[-1, 28, {
'error': ValueError
}]
]
......@@ -72,18 +72,18 @@ class LMSpaceInterfaceTests(unittest.TestCase):
['mmax', int],
['dim', int]])
def test_property_ret_type(self, attribute, expected_type):
l = LMSpace(7)
l = LMSpace(7, 5)
assert_(isinstance(getattr(l, attribute), expected_type))
class LMSpaceFunctionalityTests(unittest.TestCase):
@expand(CONSTRUCTOR_CONFIGS)
def test_constructor(self, lmax, expected):
def test_constructor(self, lmax, mmax, expected):
if 'error' in expected:
with assert_raises(expected['error']):
LMSpace(lmax)
LMSpace(lmax, mmax)
else:
l = LMSpace(lmax)
l = LMSpace(lmax, mmax)
for key, value in expected.items():
assert_equal(getattr(l, key), value)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment