Skip to content

GitLab

Commit c56bf1eb authored by theos's avatar theos
Browse files

Fixed lmhptransformation.py and hplmtransformation.py.

parent c1828a0b
Hide whitespace changes
Inline Side-by-side
nifty/operators/fft_operator/transformations/gllmtransformation.py
View file @ c56bf1eb
......@@ -21,8 +21,8 @@ class GLLMTransformation(SlicingTransformation):
# ---Mandatory properties and methods---
@staticmethod
def get_codomain(domain):
@classmethod
def get_codomain(cls, domain):
"""
Generates a compatible codomain to which transformations are
reasonable, i.e.\ an instance of the :py:class:`lm_space` class.
......@@ -37,9 +37,6 @@ class GLLMTransformation(SlicingTransformation):
codomain : LMSpace
A compatible codomain.
"""
if domain is None:
raise ValueError(about._errors.cstring(
"ERROR: cannot generate codomain for None-input"))
if not isinstance(domain, GLSpace):
raise TypeError(about._errors.cstring(
......@@ -53,10 +50,12 @@ class GLLMTransformation(SlicingTransformation):
else:
return_dtype = np.float64
return LMSpace(lmax=lmax, mmax=mmax, dtype=return_dtype)
result = LMSpace(lmax=lmax, mmax=mmax, dtype=return_dtype)
cls.check_codomain(domain, result)
return result
@classmethod
def check_codomain(cls, domain, codomain):
@staticmethod
def check_codomain(domain, codomain):
if not isinstance(domain, GLSpace):
raise TypeError(about._errors.cstring(
"ERROR: domain is not a GLSpace"))
......@@ -71,17 +70,20 @@ class GLLMTransformation(SlicingTransformation):
mmax = codomain.mmax
if lmax != mmax:
raise ValueError('ERROR: codomain has lmax != mmax.')
raise ValueError(about._errors.cstring(
'ERROR: codomain has lmax != mmax.'))
if lmax != nlat - 1:
raise ValueError('ERROR: codomain has lmax != nlat - 1.')
raise ValueError(about._errors.cstring(
'ERROR: codomain has lmax != nlat - 1.'))
if nlon != 2 * nlat - 1:
raise ValueError('ERROR: domain has nlon != 2 * nlat - 1.')
raise ValueError(about._errors.cstring(
'ERROR: domain has nlon != 2 * nlat - 1.'))
return None
def _transformation_of_slice(self, inp):
def _transformation_of_slice(self, inp, **kwargs):
nlat = self.domain.nlat
nlon = self.domain.nlon
lmax = self.codomain.lmax
......@@ -93,19 +95,14 @@ class GLLMTransformation(SlicingTransformation):
nlat=nlat,
nlon=nlon,
lmax=lmax,
mmax=mmax)
mmax=mmax,
**kwargs)
for x in (inp.real, inp.imag)]
[resultReal, resultImag] = [ltf.buildIdx(x, lmax=lmax)
for x in [resultReal, resultImag]]
# construct correct complex dtype
one = resultReal.dtype.type(1)
result_dtype = np.dtype(type(one + 1j))
result = np.empty_like(resultReal, dtype=result_dtype)
result.real = resultReal
result.imag = resultImag
result = self._combine_complex_result(resultReal, resultImag)
else:
result = self.libsharpMap2Alm(inp, nlat=nlat, nlon=nlon, lmax=lmax,
mmax=mmax)
......
nifty/operators/fft_operator/transformations/hplmtransformation.py
View file @ c56bf1eb
import numpy as np
from transformation import Transformation
from d2o import distributed_data_object
from nifty.config import dependency_injector as gdi
import nifty.nifty_utilities as utilities
from nifty.config import dependency_injector as gdi,\
about
from nifty import HPSpace, LMSpace
from slicing_transformation import SlicingTransformation
import lm_transformation_factory as ltf
hp = gdi.get('healpy')
class HPLMTransformation(Transformation):
class HPLMTransformation(SlicingTransformation):
# ---Overwritten properties and methods---
def __init__(self, domain, codomain=None, module=None):
if 'healpy' not in gdi:
raise ImportError("The module healpy is needed but not available")
if codomain is None:
self.domain = domain
self.codomain = self.get_codomain(domain)
elif self.check_codomain(domain, codomain):
self.domain = domain
self.codomain = codomain
else:
raise ValueError("ERROR: Incompatible codomain!")
raise ImportError(about._errors.cstring(
"The module healpy is needed but not available"))
@staticmethod
def get_codomain(domain):
super(HPLMTransformation, self).__init__(domain, codomain, module)
# ---Mandatory properties and methods---
@classmethod
def get_codomain(cls, domain):
"""
Generates a compatible codomain to which transformations are
reasonable, i.e.\ an instance of the :py:class:`lm_space` class.
......@@ -40,100 +38,65 @@ class HPLMTransformation(Transformation):
codomain : LMSpace
A compatible codomain.
"""
if domain is None:
raise ValueError('ERROR: cannot generate codomain for None')
if not isinstance(domain, HPSpace):
raise TypeError('ERROR: domain needs to be a HPSpace')
raise TypeError(about._errors.cstring(
"ERROR: domain needs to be a HPSpace"))
lmax = 3 * domain.nside - 1
mmax = lmax
return LMSpace(lmax=lmax, mmax=mmax, dtype=np.dtype('complex128'))
result = LMSpace(lmax=lmax, mmax=mmax, dtype=np.dtype('float64'))
cls.check_codomain(domain, result)
return result
@staticmethod
def check_codomain(domain, codomain):
if not isinstance(domain, HPSpace):
raise TypeError('ERROR: domain is not a HPSpace')
if codomain is None:
return False
raise TypeError(about._errors.cstring(
'ERROR: domain is not a HPSpace'))
if not isinstance(codomain, LMSpace):
raise TypeError('ERROR: codomain must be a LMSpace.')
raise TypeError(about._errors.cstring(
'ERROR: codomain must be a LMSpace.'))
nside = domain.nside
lmax = codomain.lmax
mmax = codomain.mmax
if (3 * nside - 1 != lmax) or (lmax != mmax):
return False
if 3 * nside - 1 != lmax:
raise ValueError(about._errors.cstring(
'ERROR: codomain has 3*nside-1 != lmax.'))
return True
if lmax != mmax:
raise ValueError(about._errors.cstring(
'ERROR: codomain has lmax != mmax.'))
def transform(self, val, axes=None, **kwargs):
"""
HP -> LM transform method.
Parameters
----------
val : np.ndarray or distributed_data_object
The value array which is to be transformed
return None
axes : None or tuple
The axes along which the transformation should take place
def _transformation_of_slice(self, inp, **kwargs):
lmax = self.codomain.lmax
mmax = self.codomain.mmax
"""
# get by number of iterations from kwargs
niter = kwargs['niter'] if 'niter' in kwargs else 0
if issubclass(inp.dtype.type, np.complexfloating):
[resultReal, resultImag] = [hp.map2alm(x.astype(np.float64,
copy=False),
lmax=lmax,
mmax=mmax,
pol=True,
use_weights=False,
**kwargs)
for x in (inp.real, inp.imag)]
if self.domain.discrete:
val = self.domain.weight(val, power=-0.5, axes=axes)
[resultReal, resultImag] = [ltf.buildIdx(x, lmax=lmax)
for x in [resultReal, resultImag]]
# shorthands for transform parameters
lmax = self.codomain.lmax
mmax = self.codomain.mmax
result = self._combine_complex_result(resultReal, resultImag)
if isinstance(val, distributed_data_object):
temp_val = val.get_full_data()
else:
temp_val = val
return_val = None
for slice_list in utilities.get_slice_list(temp_val.shape, axes):
if slice_list == [slice(None, None)]:
inp = temp_val
else:
if return_val is None:
return_val = np.empty_like(temp_val)
inp = temp_val[slice_list]
if self.domain.dtype == np.dtype('complex128'):
inpReal = hp.map2alm(
np.real(inp).astype(np.float64, copy=False),
lmax=lmax, mmax=mmax, iter=niter, pol=True,
use_weights=False, datapath=None)
inpImg = hp.map2alm(
np.imag(inp).astype(np.float64, copy=False),
lmax=lmax, mmax=mmax, iter=niter, pol=True,
use_weights=False, datapath=None)
inpReal = ltf.buildIdx(inpReal, lmax=lmax)
inpImg = ltf.buildIdx(inpImg, lmax=lmax)
inp = inpReal + inpImg * 1j
else:
inp = hp.map2alm(inp.astype(np.float64, copy=False),
lmax=lmax, mmax=mmax, iter=niter, pol=True,
use_weights=False, datapath=None)
inp = ltf.buildIdx(inp, lmax=lmax)
if slice_list == [slice(None, None)]:
return_val = inp
else:
return_val[slice_list] = inp
if isinstance(val, distributed_data_object):
new_val = val.copy_empty(dtype=self.codomain.dtype)
new_val.set_full_data(return_val, copy=False)
else:
return_val = return_val.astype(self.codomain.dtype, copy=False)
result = hp.map2alm(inp.astype(np.float64, copy=False),
lmax=lmax, mmax=mmax, pol=True,
use_weights=False)
result = ltf.buildIdx(result, lmax=lmax)
return return_val
return result
nifty/operators/fft_operator/transformations/lmgltransformation.py
View file @ c56bf1eb
......@@ -22,8 +22,8 @@ class LMGLTransformation(SlicingTransformation):
# ---Mandatory properties and methods---
@staticmethod
def get_codomain(domain):
@classmethod
def get_codomain(cls, domain):
"""
Generates a compatible codomain to which transformations are
reasonable, i.e.\ a pixelization of the two-sphere.
......@@ -45,10 +45,6 @@ class LMGLTransformation(SlicingTransformation):
harmonic transforms revisited";
`arXiv:1303.4945 <http://www.arxiv.org/abs/1303.4945>`_
"""
if domain is None:
raise ValueError(about._errors.cstring(
'ERROR: cannot generate codomain for None-input'))
if not isinstance(domain, LMSpace):
raise TypeError(about._errors.cstring(
'ERROR: domain needs to be a LMSpace'))
......@@ -61,15 +57,19 @@ class LMGLTransformation(SlicingTransformation):
nlat = domain.lmax + 1
nlon = domain.lmax * 2 + 1
return GLSpace(nlat=nlat, nlon=nlon, dtype=new_dtype)
result = GLSpace(nlat=nlat, nlon=nlon, dtype=new_dtype)
cls.check_codomain(domain, result)
return result
@staticmethod
def check_codomain(domain, codomain):
if not isinstance(domain, LMSpace):
raise TypeError('ERROR: domain is not a LMSpace')
raise TypeError(about._errors.cstring(
'ERROR: domain is not a LMSpace'))
if not isinstance(codomain, GLSpace):
raise TypeError('ERROR: codomain must be a GLSpace.')
raise TypeError(about._errors.cstring(
'ERROR: codomain must be a GLSpace.'))
nlat = codomain.nlat
nlon = codomain.nlon
......@@ -77,17 +77,20 @@ class LMGLTransformation(SlicingTransformation):
mmax = domain.mmax
if lmax != mmax:
raise ValueError('ERROR: domain has lmax != mmax.')
raise ValueError(about._errors.cstring(
'ERROR: domain has lmax != mmax.'))
if nlat != lmax + 1:
raise ValueError('ERROR: codomain has nlat != lmax + 1.')
raise ValueError(about._errors.cstring(
'ERROR: codomain has nlat != lmax + 1.'))
if nlon != 2 * lmax + 1:
raise ValueError('ERROR: domain has nlon != 2 * lmax + 1.')
raise ValueError(about._errors.cstring(
'ERROR: domain has nlon != 2 * lmax + 1.'))
return None
def _transformation_of_slice(self, inp):
def _transformation_of_slice(self, inp, **kwargs):
nlat = self.codomain.nlat
nlon = self.codomain.nlon
lmax = self.domain.lmax
......@@ -102,16 +105,12 @@ class LMGLTransformation(SlicingTransformation):
nlon=nlon,
lmax=lmax,
mmax=mmax,
cl=False)
cl=False,
**kwargs)
for x in [resultReal, resultImag]]
# construct correct complex dtype
one = resultReal.dtype.type(1)
result_dtype = np.dtype(type(one + 1j))
result = self._combine_complex_result(resultReal, resultImag)
result = np.empty_like(resultReal, dtype=result_dtype)
result.real = resultReal
result.imag = resultImag
else:
result = ltf.buildLm(inp, lmax=lmax)
result = self.libsharpAlm2Map(result, nlat=nlat, nlon=nlon,
......
nifty/operators/fft_operator/transformations/lmhptransformation.py
View file @ c56bf1eb
import numpy as np
from transformation import Transformation
from d2o import distributed_data_object
from nifty.config import dependency_injector as gdi
import nifty.nifty_utilities as utilities
from nifty.config import dependency_injector as gdi,\
about
from nifty import HPSpace, LMSpace
from slicing_transformation import SlicingTransformation
import lm_transformation_factory as ltf
hp = gdi.get('healpy')
class LMHPTransformation(Transformation):
class LMHPTransformation(SlicingTransformation):
# ---Overwritten properties and methods---
def __init__(self, domain, codomain=None, module=None):
if gdi.get('healpy') is None:
raise ImportError(
"The module libsharp is needed but not available.")
raise ImportError(about._errors.cstring(
"The module libsharp is needed but not available."))
if self.check_codomain(domain, codomain):
self.domain = domain
self.codomain = codomain
else:
raise ValueError("ERROR: Incompatible codomain!")
super(LMHPTransformation, self).__init__(domain, codomain, module)
@staticmethod
def get_codomain(domain):
# ---Mandatory properties and methods---
@classmethod
def get_codomain(cls, domain):
"""
Generates a compatible codomain to which transformations are
reasonable, i.e.\ a pixelization of the two-sphere.
......@@ -44,100 +43,67 @@ class LMHPTransformation(Transformation):
High-Resolution Discretization and Fast Analysis of Data
Distributed on the Sphere", *ApJ* 622..759G.
"""
if domain is None:
raise ValueError('ERROR: cannot generate codomain for None')
if not isinstance(domain, LMSpace):
raise TypeError('ERROR: domain needs to be a LMSpace')
raise TypeError(about._errors.cstring(
'ERROR: domain needs to be a LMSpace'))
nside = (domain.lmax + 1) // 3
return HPSpace(nside=nside)
result = HPSpace(nside=nside)
cls.check_codomain(domain, result)
return result
@staticmethod
def check_codomain(domain, codomain):
if not isinstance(domain, LMSpace):
raise TypeError('ERROR: domain is not a LMSpace')
if codomain is None:
return False
raise TypeError(about._errors.cstring(
'ERROR: domain is not a LMSpace'))
if not isinstance(codomain, HPSpace):
raise TypeError('ERROR: codomain must be a HPSpace.')
raise TypeError(about._errors.cstring(
'ERROR: codomain must be a HPSpace.'))
nside = codomain.nside
lmax = domain.lmax
mmax = domain.mmax
if (lmax != mmax) or (3 * nside - 1 != lmax):
return False
if lmax != mmax:
raise ValueError(about._errors.cstring(
'ERROR: domain has lmax != mmax.'))
return True
if 3*nside - 1 != lmax:
raise ValueError(about._errors.cstring(
'ERROR: codomain has 3*nside -1 != lmax.'))
def transform(self, val, axes=None, **kwargs):
"""
LM -> HP transform method.
return None
Parameters
----------
val : np.ndarray or distributed_data_object
The value array which is to be transformed
def _transformation_of_slice(self, inp, **kwargs):
nside = self.codomain.nside
lmax = self.domain.lmax
mmax = self.domain.mmax
axes : None or tuple
The axes along which the transformation should take place
if issubclass(inp.dtype.type, np.complexfloating):
[resultReal, resultImag] = [ltf.buildLm(x, lmax=lmax)
for x in (inp.real, inp.imag)]
[resultReal, resultImag] = [hp.map2alm(x.astype(np.complex128,
copy=False),
nside,
lmax=lmax,
mmax=mmax,
pixwin=False,
fwhm=0.0,
sigma=None,
pol=True,
inplace=False,
**kwargs)
for x in [resultReal, resultImag]]
result = self._combine_complex_result(resultReal, resultImag)
"""
if isinstance(val, distributed_data_object):
temp_val = val.get_full_data()
else:
temp_val = val
return_val = None
for slice_list in utilities.get_slice_list(temp_val.shape, axes):
if slice_list == [slice(None, None)]:
inp = temp_val
else:
if return_val is None:
return_val = np.empty_like(temp_val)
inp = temp_val[slice_list]
nside = self.codomain.nside
lmax = self.domain.lmax
mmax = self.domain.mmax
if self.domain.dtype == np.dtype('complex128'):
inpReal = np.real(inp)
inpImag = np.imag(inp)
inpReal = ltf.buildLm(inpReal,lmax=lmax)
inpImag = ltf.buildLm(inpImag,lmax=lmax)
inpReal = hp.alm2map(inpReal.astype(np.complex128, copy=False), nside,
lmax=lmax, mmax=mmax,
pixwin=False, fwhm=0.0, sigma=None,
pol=True, inplace=False)
inpImag = hp.alm2map(inpImag.astype(np.complex128, copy=False), nside,
lmax=lmax, mmax=mmax,
pixwin=False, fwhm=0.0, sigma=None,
pol=True, inplace=False)
inp = inpReal+inpImag*(1j)
else:
inp = ltf.buildLm(inp, lmax=lmax)
inp = hp.alm2map(inp.astype(np.complex128, copy=False), nside,
lmax=lmax, mmax=mmax,
pixwin=False, fwhm=0.0, sigma=None,
pol=True, inplace=False)
if slice_list == [slice(None, None)]:
return_val = inp
else:
return_val[slice_list] = inp
# re-weight if discrete
if self.codomain.discrete:
val = self.codomain.weight(val, power=0.5, axes=axes)
if isinstance(val, distributed_data_object):
new_val = val.copy_empty(dtype=self.codomain.dtype)
new_val.set_full_data(return_val, copy=False)
else:
return_val = return_val.astype(self.codomain.dtype, copy=False)
result = ltf.buildLm(inp, lmax=lmax)
result = hp.alm2map(result.astype(np.complex128, copy=False),
nside, lmax=lmax, mmax=mmax, pixwin=False,
fwhm=0.0, sigma=None, pol=True, inplace=False)
return return_val
return result
nifty/operators/fft_operator/transformations/rgrgtransformation.py
View file @ c56bf1eb
......@@ -33,8 +33,8 @@ class RGRGTransformation(Transformation):
else:
raise ValueError('ERROR: unknow FFT module:' + module)
@staticmethod
def get_codomain(domain, dtype=None, zerocenter=None):
@classmethod
def get_codomain(cls, domain, dtype=None, zerocenter=None):
"""
Generates a compatible codomain to which transformations are
reasonable, i.e.\ either a shifted grid or a Fourier conjugate
......@@ -78,7 +78,7 @@ class RGRGTransformation(Transformation):
distances=distances,
harmonic=(not domain.harmonic),
dtype=dtype)
cls.check_codomain(domain, new_space)
return new_space
@staticmethod
......
nifty/operators/fft_operator/transformations/slicing_transformation.py
View file @ c56bf1eb
......@@ -25,12 +25,23 @@ class SlicingTransformation(Transformation):
data = val.get_data(slice_list, copy=False)
data = data.get_full_data()
data = self._transformation_of_slice(data)
data = self._transformation_of_slice(data, **kwargs)
return_val.set_data(data=data, to_key=slice_list, copy=False)
return return_val
def _combine_complex_result(self, resultReal, resultImag):
# construct correct complex dtype
one = resultReal.dtype.type(1)