Skip to content

GitLab

Commit 0fcd7cdc authored by Theo Steininger's avatar Theo Steininger
Browse files

Merge branch 'lm_space_complex_transformsV2' into 'feature/field_multiple_space' 

Lm space complex transforms v2



See merge request !27
parents 3229c061 0f43abfc
Expand all Hide whitespace changes
Inline Side-by-side
nifty/field.py
View file @ 0fcd7cdc
......@@ -5,7 +5,7 @@ from d2o import distributed_data_object,\
STRATEGIES as DISTRIBUTION_STRATEGIES
from nifty.config import about,\
nifty_configuration as gc,\
nifty_configuration as gc
from nifty.field_types import FieldType
......
nifty/operators/fft_operator/fft_operator.py
View file @ 0fcd7cdc
from nifty.config import about
import nifty.nifty_utilities as utilities
from nifty.spaces import RGSpace,\
GLSpace,\
HPSpace,\
LMSpace
from nifty.operators.linear_operator import LinearOperator
from transformations import TransformationFactory
from transformations import RGRGTransformation,\
LMGLTransformation,\
LMHPTransformation,\
GLLMTransformation,\
HPLMTransformation,\
TransformationCache
class FFTOperator(LinearOperator):
# ---Class attributes---
default_codomain_dictionary = {RGSpace: RGSpace,
HPSpace: LMSpace,
GLSpace: LMSpace,
LMSpace: HPSpace,
}
transformation_dictionary = {(RGSpace, RGSpace): RGRGTransformation,
(HPSpace, LMSpace): HPLMTransformation,
(GLSpace, LMSpace): GLLMTransformation,
(LMSpace, HPSpace): LMHPTransformation,
(LMSpace, GLSpace): LMGLTransformation
}
# ---Overwritten properties and methods---
def __init__(self, domain=(), field_type=(), target=None):
def __init__(self, domain=(), field_type=(), target=None, module=None):
super(FFTOperator, self).__init__(domain=domain,
field_type=field_type)
# Initialize domain and target
if len(self.domain) != 1:
raise ValueError(about._errors.cstring(
'ERROR: TransformationOperator accepts only exactly one '
......@@ -24,17 +51,30 @@ class FFTOperator(LinearOperator):
))
if target is None:
target = utilities.get_default_codomain(self.domain[0])
target = (self.get_default_codomain(self.domain[0]), )
self._target = self._parse_domain(target)
self._forward_transformation = TransformationFactory.create(
self.domain[0], self.target[0]
)
self._inverse_transformation = TransformationFactory.create(
self.target[0], self.domain[0]
)
# Create transformation instances
try:
forward_class = self.transformation_dictionary[
(self.domain[0].__class__, self.target[0].__class__)]
except KeyError:
raise TypeError(about._errors.cstring(
"ERROR: No forward transformation for domain-target pair "
"found."))
try:
backward_class = self.transformation_dictionary[
(self.target[0].__class__, self.domain[0].__class__)]
except KeyError:
raise TypeError(about._errors.cstring(
"ERROR: No backward transformation for domain-target pair "
"found."))
self._forward_transformation = TransformationCache.create(
forward_class, self.domain[0], self.target[0], module=module)
self._backward_transformation = TransformationCache.create(
backward_class, self.target[0], self.domain[0], module=module)
def _times(self, x, spaces, types):
spaces = utilities.cast_axis_to_tuple(spaces, len(x.domain))
......@@ -69,7 +109,7 @@ class FFTOperator(LinearOperator):
else:
axes = x.domain_axes[spaces[0]]
new_val = self._inverse_transformation.transform(x.val, axes=axes)
new_val = self._backward_transformation.transform(x.val, axes=axes)
if spaces is None:
result_domain = self.domain
......@@ -99,3 +139,22 @@ class FFTOperator(LinearOperator):
@property
def unitary(self):
return True
# ---Added properties and methods---
@classmethod
def get_default_codomain(cls, domain):
domain_class = domain.__class__
try:
codomain_class = cls.default_codomain_dictionary[domain_class]
except KeyError:
raise TypeError(about._errors.cstring("ERROR: unknown domain"))
try:
transform_class = cls.transformation_dictionary[(domain_class,
codomain_class)]
except KeyError:
raise TypeError(about._errors.cstring(
"ERROR: No transformation for domain-codomain pair found."))
return transform_class.get_codomain(domain)
nifty/operators/fft_operator/transformations/__init__.py
View file @ 0fcd7cdc
......@@ -4,4 +4,4 @@ from hplmtransformation import HPLMTransformation
from lmgltransformation import LMGLTransformation
from lmhptransformation import LMHPTransformation
from transformation_factory import TransformationFactory
\ No newline at end of file
from transformation_cache import TransformationCache
\ No newline at end of file
nifty/operators/fft_operator/transformations/gllmtransformation.py
View file @ 0fcd7cdc
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 GLSpace, LMSpace
from slicing_transformation import SlicingTransformation
import lm_transformation_factory as ltf
gl = gdi.get('libsharp_wrapper_gl')
libsharp = gdi.get('libsharp_wrapper_gl')
class GLLMTransformation(Transformation):
class GLLMTransformation(SlicingTransformation):
# ---Overwritten properties and methods---
def __init__(self, domain, codomain=None, module=None):
if 'libsharp_wrapper_gl' not in gdi:
raise ImportError("The module libsharp 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 libsharp is needed but not available."))
@staticmethod
def get_codomain(domain):
super(GLLMTransformation, 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.
......@@ -38,96 +37,89 @@ class GLLMTransformation(Transformation):
codomain : LMSpace
A compatible codomain.
"""
if domain is None:
raise ValueError('ERROR: cannot generate codomain for None')
if not isinstance(domain, GLSpace):
raise TypeError('ERROR: domain needs to be a GLSpace')
raise TypeError(about._errors.cstring(
"ERROR: domain needs to be a GLSpace"))
nlat = domain.nlat
lmax = nlat - 1
mmax = nlat - 1
if domain.dtype == np.dtype('float32'):
return LMSpace(lmax=lmax, mmax=mmax, dtype=np.complex64)
return_dtype = np.float32
else:
return LMSpace(lmax=lmax, mmax=mmax, dtype=np.complex128)
return_dtype = np.float64
result = LMSpace(lmax=lmax, mmax=mmax, dtype=return_dtype)
cls.check_codomain(domain, result)
return result
@staticmethod
def check_codomain(domain, codomain):
if not isinstance(domain, GLSpace):
raise TypeError('ERROR: domain is not a GLSpace')
if codomain is None:
return False
raise TypeError(about._errors.cstring(
"ERROR: domain is not a GLSpace"))
if not isinstance(codomain, LMSpace):
raise TypeError('ERROR: codomain must be a LMSpace.')
raise TypeError(about._errors.cstring(
"ERROR: codomain must be a LMSpace."))
nlat = domain.nlat
nlon = domain.nlon
lmax = codomain.lmax
mmax = codomain.mmax
if (nlon != 2 * nlat - 1) or (lmax != nlat - 1) or (lmax != mmax):
return False
if lmax != mmax:
raise ValueError(about._errors.cstring(
'ERROR: codomain has lmax != mmax.'))
return True
def transform(self, val, axes=None, **kwargs):
"""
GL -> LM transform method.
if lmax != nlat - 1:
raise ValueError(about._errors.cstring(
'ERROR: codomain has lmax != nlat - 1.'))
Parameters
----------
val : np.ndarray or distributed_data_object
The value array which is to be transformed
if nlon != 2 * nlat - 1:
raise ValueError(about._errors.cstring(
'ERROR: domain has nlon != 2 * nlat - 1.'))
axes : None or tuple
The axes along which the transformation should take place
return None
"""
if self.domain.discrete:
val = self.domain.weight(val, power=-0.5, axes=axes)
# shorthands for transform parameters
def _transformation_of_slice(self, inp, **kwargs):
nlat = self.domain.nlat
nlon = self.domain.nlon
lmax = self.codomain.lmax
mmax = self.codomain.mmax
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('float32'):
inp = gl.map2alm_f(inp,
nlat=nlat, nlon=nlon,
lmax=lmax, mmax=mmax)
else:
inp = gl.map2alm(inp,
nlat=nlat, nlon=nlon,
lmax=lmax, mmax=mmax)
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)
if issubclass(inp.dtype.type, np.complexfloating):
[resultReal, resultImag] = [self.libsharpMap2Alm(x,
nlat=nlat,
nlon=nlon,
lmax=lmax,
mmax=mmax,
**kwargs)
for x in (inp.real, inp.imag)]
[resultReal, resultImag] = [ltf.buildIdx(x, lmax=lmax)
for x in [resultReal, resultImag]]
result = self._combine_complex_result(resultReal, resultImag)
else:
return_val = return_val.astype(self.codomain.dtype, copy=False)
result = self.libsharpMap2Alm(inp, nlat=nlat, nlon=nlon, lmax=lmax,
mmax=mmax)
result = ltf.buildIdx(result, lmax=lmax)
return result
return return_val
# ---Added properties and methods---
def libsharpMap2Alm(self, inp, **kwargs):
if inp.dtype == np.dtype('float32'):
return libsharp.map2alm_f(inp, **kwargs)
elif inp.dtype == np.dtype('float64'):
return libsharp.map2alm(inp, **kwargs)
else:
about.warnings.cprint("WARNING: performing dtype conversion for "
"libsharp compatibility.")
casted_inp = inp.astype(np.dtype('float64'), copy=False)
result = libsharp.map2alm(casted_inp, **kwargs)
return result
\ No newline at end of file
nifty/operators/fft_operator/transformations/hplmtransformation.py
View file @ 0fcd7cdc
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.
......@@ -38,87 +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]
inp = hp.map2alm(inp.astype(np.float64, copy=False),
lmax=lmax, mmax=mmax, iter=niter, pol=True,
use_weights=False, datapath=None)
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/lm_transformation_factory.pyx 0 → 100644
View file @ 0fcd7cdc
import numpy as np
cimport numpy as np
def buildLm(inp, **kwargs):
if inp.dtype == np.dtype('float32'):
return _buildLm_f(inp, **kwargs)
else:
return _buildLm(inp, **kwargs)
def buildIdx(inp, **kwargs):
if inp.dtype == np.dtype('complex64'):
return _buildIdx_f(inp, **kwargs)
else:
return _buildIdx(inp, **kwargs)
cpdef np.ndarray[np.float32_t, ndim=1] _buildIdx_f(np.ndarray[np
.complex64_t, ndim=1] nr, np.int_t lmax):
cdef np.int size = (lmax+1)*(lmax+1)
cdef np.ndarray res=np.zeros([size], dtype=np.complex64)
cdef np.ndarray final=np.zeros([size], dtype=np.float32)
res[0:lmax+1] = nr[0:lmax+1]
for i in xrange(len(nr)-lmax-1):
res[i*2+lmax+1] = nr[i+lmax+1]
res[i*2+lmax+2] = np.conjugate(nr[i+lmax+1])
final = _realify_f(res, lmax)
return final
cpdef np.ndarray[np.float32_t, ndim=1] _realify_f(np.ndarray[np.complex64_t, ndim=1] nr, np.int_t lmax):
cdef np.ndarray resi=np.zeros([len(nr)], dtype=np.float32)
resi[0:lmax+1] = np.real(nr[0:lmax+1])
for i in xrange(lmax+1,len(nr),2):
mi = int(np.ceil(((2*lmax+1)-np.sqrt((2*lmax+1)*(2*lmax+1)-4*(i-lmax)+1))/2))
resi[i]=np.sqrt(2)*np.real(nr[i])*(-1)**(mi*mi)
resi[i+1]=np.sqrt(2)*np.imag(nr[i])*(-1)**(mi*mi)
return resi
cpdef np.ndarray[np.float64_t, ndim=1] _buildIdx(np.ndarray[np.complex128_t,
ndim=1] nr, np.int_t lmax):
cdef np.int size = (lmax+1)*(lmax+1)
cdef np.ndarray res=np.zeros([size], dtype=np.complex128)
cdef np.ndarray final=np.zeros([size], dtype=np.float64)
res[0:lmax+1] = nr[0:lmax+1]
for i in xrange(len(nr)-lmax-1):
res[i*2+lmax+1] = nr[i+lmax+1]
res[i*2+lmax+2] = np.conjugate(nr[i+lmax+1])
final = _realify(res, lmax)
return final
cpdef np.ndarray[np.float64_t, ndim=1] _realify(np.ndarray[np.complex128_t, ndim=1] nr, np.int_t lmax):
cdef np.ndarray resi=np.zeros([len(nr)], dtype=np.float64)
resi[0:lmax+1] = np.real(nr[0:lmax+1])
for i in xrange(lmax+1,len(nr),2):
mi = int(np.ceil(((2*lmax+1)-np.sqrt((2*lmax+1)*(2*lmax+1)-4*(i-lmax)+1))/2))
resi[i]=np.sqrt(2)*np.real(nr[i])*(-1)**(mi*mi)
resi[i+1]=np.sqrt(2)*np.imag(nr[i])*(-1)**(mi*mi)
return resi
cpdef np.ndarray[np.complex64_t, ndim=1] _buildLm_f(np.ndarray[np.float32_t,
ndim=1] nr, np.int_t lmax):