Commit ef0fab1b authored by Martin Reinecke's avatar Martin Reinecke
Browse files

consolidation

parent aaaaa209
Pipeline #17996 passed with stage
in 4 minutes and 51 seconds
......@@ -16,5 +16,4 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
from .transformations import *
from .fft_operator import FFTOperator
......@@ -17,17 +17,14 @@
# and financially supported by the Studienstiftung des deutschen Volkes.
from ... import Field, nifty_utilities as utilities
from ...spaces import RGSpace,\
GLSpace,\
HPSpace,\
LMSpace
from ...spaces import RGSpace, GLSpace, HPSpace, LMSpace
from ..linear_operator import LinearOperator
from .transformations import RGRGTransformation,\
LMGLTransformation,\
LMHPTransformation,\
GLLMTransformation,\
HPLMTransformation
from .fft_operator_support import RGRGTransformation,\
LMGLTransformation,\
LMHPTransformation,\
GLLMTransformation,\
HPLMTransformation
class FFTOperator(LinearOperator):
......
......@@ -17,8 +17,25 @@
# and financially supported by the Studienstiftung des deutschen Volkes.
from __future__ import division
import abc
from future.utils import with_metaclass
import numpy as np
from .transformation import Transformation
from ... import nifty_utilities as utilities
class Transformation(with_metaclass(abc.ABCMeta,
type('NewBase', (object,), {}))):
def __init__(self, domain, codomain):
self.domain = domain
self.codomain = codomain
@abc.abstractproperty
def unitary(self):
raise NotImplementedError
def transform(self, val, axes=None):
raise NotImplementedError
class RGRGTransformation(Transformation):
......@@ -60,7 +77,7 @@ class RGRGTransformation(Transformation):
The axes along which the transformation should take place
"""
fct=1.
fct = 1.
if self.codomain.harmonic:
# correct for forward fft.
# naively one would set power to 0.5 here in order to
......@@ -75,10 +92,12 @@ class RGRGTransformation(Transformation):
Tval_imag = self._transform_helper(val.imag, axes)
if self.codomain.harmonic:
Tval_real.real += Tval_real.imag
Tval_real.imag = Tval_imag.real + Tval_imag.imag
Tval_imag.real += Tval_imag.imag
Tval_real.imag = Tval_imag.real
else:
Tval_real.real -= Tval_real.imag
Tval_real.imag = Tval_imag.real - Tval_imag.imag
Tval_imag.real -= Tval_imag.imag
Tval_real.imag = Tval_imag.real
Tval = Tval_real
else:
......@@ -96,3 +115,155 @@ class RGRGTransformation(Transformation):
Tval *= fct
return Tval
class SlicingTransformation(Transformation):
def transform(self, val, axes=None):
return_shape = np.array(val.shape)
return_shape[list(axes)] = self.codomain.shape
return_shape = tuple(return_shape)
return_val = np.empty(return_shape, dtype=val.dtype)
for slice_list in utilities.get_slice_list(val.shape, axes):
return_val[slice_list] = self._transformation_of_slice(
val[slice_list])
return return_val
@abc.abstractmethod
def _transformation_of_slice(self, inp):
raise NotImplementedError
def buildLm(nr, lmax):
new_dtype = np.result_type(nr.dtype, np.complex64)
size = (len(nr)-lmax-1)//2+lmax+1
res = np.empty([size], dtype=new_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):
if nr.dtype == np.dtype('complex64'):
new_dtype = np.float32
elif nr.dtype == np.dtype('complex128'):
new_dtype = np.float64
else:
raise TypeError("dtype of nr not supported.")
size = (lmax+1)*(lmax+1)
final = np.empty(size, dtype=new_dtype)
final[0:lmax+1] = nr[0:lmax+1].real
final[lmax+1::2] = np.sqrt(2)*nr[lmax+1:].real
final[lmax+2::2] = np.sqrt(2)*nr[lmax+1:].imag
return final
class HPLMTransformation(SlicingTransformation):
def __init__(self, domain, codomain=None):
super(HPLMTransformation, self).__init__(domain, codomain)
@property
def unitary(self):
return False
def _transformation_of_slice(self, inp):
from pyHealpix import map2alm
lmax = self.codomain.lmax
mmax = lmax
if issubclass(inp.dtype.type, np.complexfloating):
rr = map2alm(inp.real, lmax, mmax)
rr = buildIdx(rr, lmax=lmax)
ri = map2alm(inp.imag, lmax, mmax)
ri = buildIdx(ri, lmax=lmax)
return rr + 1j*ri
else:
rr = map2alm(inp, lmax, mmax)
return buildIdx(rr, lmax=lmax)
class LMHPTransformation(SlicingTransformation):
def __init__(self, domain, codomain=None):
super(LMHPTransformation, self).__init__(domain, codomain)
@property
def unitary(self):
return False
def _transformation_of_slice(self, inp):
from pyHealpix import alm2map
nside = self.codomain.nside
lmax = self.domain.lmax
mmax = lmax
if issubclass(inp.dtype.type, np.complexfloating):
rr = buildLm(inp.real, lmax=lmax)
ri = buildLm(inp.imag, lmax=lmax)
rr = alm2map(rr, lmax, mmax, nside)
ri = alm2map(ri, lmax, mmax, nside)
return rr + 1j*ri
else:
rr = buildLm(inp, lmax=lmax)
return alm2map(rr, lmax, mmax, nside)
class GLLMTransformation(SlicingTransformation):
def __init__(self, domain, codomain=None):
super(GLLMTransformation, self).__init__(domain, codomain)
@property
def unitary(self):
return False
def _transformation_of_slice(self, inp):
from pyHealpix import sharpjob_d
lmax = self.codomain.lmax
mmax = self.codomain.mmax
sjob = sharpjob_d()
sjob.set_Gauss_geometry(self.domain.nlat, self.domain.nlon)
sjob.set_triangular_alm_info(lmax, mmax)
if issubclass(inp.dtype.type, np.complexfloating):
rr = sjob.map2alm(inp.real)
rr = buildIdx(rr, lmax=lmax)
ri = sjob.map2alm(inp.imag)
ri = buildIdx(ri, lmax=lmax)
return rr + 1j*ri
else:
rr = sjob.map2alm(inp)
return buildIdx(rr, lmax=lmax)
class LMGLTransformation(SlicingTransformation):
def __init__(self, domain, codomain=None):
super(LMGLTransformation, self).__init__(domain, codomain)
@property
def unitary(self):
return False
def _transformation_of_slice(self, inp):
from pyHealpix import sharpjob_d
lmax = self.domain.lmax
mmax = self.domain.mmax
sjob = sharpjob_d()
sjob.set_Gauss_geometry(self.codomain.nlat, self.codomain.nlon)
sjob.set_triangular_alm_info(lmax, mmax)
if issubclass(inp.dtype.type, np.complexfloating):
rr = buildLm(inp.real, lmax=lmax)
ri = buildLm(inp.imag, lmax=lmax)
return sjob.alm2map(rr) + 1j*sjob.alm2map(ri)
else:
result = buildLm(inp, lmax=lmax)
return sjob.alm2map(result)
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2017 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
from .rgrgtransformation import RGRGTransformation
from .gllmtransformation import GLLMTransformation
from .hplmtransformation import HPLMTransformation
from .lmgltransformation import LMGLTransformation
from .lmhptransformation import LMHPTransformation
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2017 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
import numpy as np
from .slicing_transformation import SlicingTransformation
from . import lm_transformation_helper
class GLLMTransformation(SlicingTransformation):
def __init__(self, domain, codomain=None):
super(GLLMTransformation, self).__init__(domain, codomain)
@property
def unitary(self):
return False
def _transformation_of_slice(self, inp):
from pyHealpix import sharpjob_d
lmax = self.codomain.lmax
mmax = self.codomain.mmax
sjob = sharpjob_d()
sjob.set_Gauss_geometry(self.domain.nlat, self.domain.nlon)
sjob.set_triangular_alm_info(lmax, mmax)
if issubclass(inp.dtype.type, np.complexfloating):
rr = sjob.map2alm(inp.real)
rr = lm_transformation_helper.buildIdx(rr, lmax=lmax)
ri = sjob.map2alm(inp.imag)
ri = lm_transformation_helper.buildIdx(ri, lmax=lmax)
return rr + 1j*ri
else:
rr = sjob.map2alm(inp)
return lm_transformation_helper.buildIdx(rr, lmax=lmax)
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2017 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
import numpy as np
from .slicing_transformation import SlicingTransformation
from . import lm_transformation_helper
class HPLMTransformation(SlicingTransformation):
def __init__(self, domain, codomain=None):
super(HPLMTransformation, self).__init__(domain, codomain)
@property
def unitary(self):
return False
def _transformation_of_slice(self, inp):
from pyHealpix import map2alm
lmax = self.codomain.lmax
mmax = lmax
if issubclass(inp.dtype.type, np.complexfloating):
rr = map2alm(inp.real, lmax, mmax)
rr = lm_transformation_helper.buildIdx(rr, lmax=lmax)
ri = map2alm(inp.imag, lmax, mmax)
ri = lm_transformation_helper.buildIdx(ri, lmax=lmax)
return rr + 1j*ri
else:
rr = map2alm(inp, lmax, mmax)
return lm_transformation_helper.buildIdx(rr, lmax=lmax)
from __future__ import division
import numpy as np
def buildLm(nr, lmax):
new_dtype = np.result_type(nr.dtype, np.complex64)
size = (len(nr)-lmax-1)//2+lmax+1
res = np.zeros([size], dtype=new_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):
if nr.dtype == np.dtype('complex64'):
new_dtype = np.float32
elif nr.dtype == np.dtype('complex128'):
new_dtype = np.float64
else:
raise TypeError("dtype of nr not supported.")
size = (lmax+1)*(lmax+1)
final = np.zeros([size], dtype=new_dtype)
final[0:lmax+1] = nr[0:lmax+1].real
final[lmax+1::2] = np.sqrt(2)*nr[lmax+1:].real
final[lmax+2::2] = np.sqrt(2)*nr[lmax+1:].imag
return final
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2017 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
import numpy as np
from .slicing_transformation import SlicingTransformation
from . import lm_transformation_helper
class LMGLTransformation(SlicingTransformation):
def __init__(self, domain, codomain=None):
super(LMGLTransformation, self).__init__(domain, codomain)
@property
def unitary(self):
return False
def _transformation_of_slice(self, inp):
from pyHealpix import sharpjob_d
lmax = self.domain.lmax
mmax = self.domain.mmax
sjob = sharpjob_d()
sjob.set_Gauss_geometry(self.codomain.nlat, self.codomain.nlon)
sjob.set_triangular_alm_info(lmax, mmax)
if issubclass(inp.dtype.type, np.complexfloating):
rr = lm_transformation_helper.buildLm(inp.real, lmax=lmax)
ri = lm_transformation_helper.buildLm(inp.imag, lmax=lmax)
return sjob.alm2map(rr) + 1j*sjob.alm2map(ri)
else:
result = lm_transformation_helper.buildLm(inp, lmax=lmax)
return sjob.alm2map(result)
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2017 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
import numpy as np
from .slicing_transformation import SlicingTransformation
from . import lm_transformation_helper
class LMHPTransformation(SlicingTransformation):
def __init__(self, domain, codomain=None):
super(LMHPTransformation, self).__init__(domain, codomain)
@property
def unitary(self):
return False
def _transformation_of_slice(self, inp):
from pyHealpix import alm2map
nside = self.codomain.nside
lmax = self.domain.lmax
mmax = lmax
if issubclass(inp.dtype.type, np.complexfloating):
rr = lm_transformation_helper.buildLm(inp.real, lmax=lmax)
ri = lm_transformation_helper.buildLm(inp.imag, lmax=lmax)
rr = alm2map(rr, lmax, mmax, nside)
ri = alm2map(ri, lmax, mmax, nside)
return rr + 1j*ri
else:
rr = lm_transformation_helper.buildLm(inp, lmax=lmax)
return alm2map(rr, lmax, mmax, nside)
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2017 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
import abc
import numpy as np
from .... import nifty_utilities as utilities
from .transformation import Transformation
class SlicingTransformation(Transformation):
def transform(self, val, axes=None):
return_shape = np.array(val.shape)
return_shape[list(axes)] = self.codomain.shape
return_shape = tuple(return_shape)
return_val = np.empty(return_shape, dtype=val.dtype)
for slice_list in utilities.get_slice_list(val.shape, axes):
return_val[slice_list] = self._transformation_of_slice(
val[slice_list])
return return_val
@abc.abstractmethod
def _transformation_of_slice(self, inp):
raise NotImplementedError
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2017 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
import abc
from future.utils import with_metaclass
class Transformation(with_metaclass(abc.ABCMeta,
type('NewBase', (object,), {}))):
def __init__(self, domain, codomain):
self.domain = domain
self.codomain = codomain
@abc.abstractproperty
def unitary(self):
raise NotImplementedError
def transform(self, val, axes=None):
raise NotImplementedError
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