fft_operator.py 8.59 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12
# 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/>.
Theo Steininger's avatar
Theo Steininger committed
13 14 15 16 17
#
# 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.
18

19 20
import numpy as np

Martin Reinecke's avatar
Martin Reinecke committed
21 22 23 24 25
from ... import nifty_utilities as utilities
from ...spaces import RGSpace,\
                      GLSpace,\
                      HPSpace,\
                      LMSpace
26

Martin Reinecke's avatar
Martin Reinecke committed
27
from ..linear_operator import LinearOperator
Martin Reinecke's avatar
Martin Reinecke committed
28
from .transformations import RGRGTransformation,\
29 30 31
                            LMGLTransformation,\
                            LMHPTransformation,\
                            GLLMTransformation,\
32
                            HPLMTransformation
Jait Dixit's avatar
Jait Dixit committed
33 34


Jait Dixit's avatar
Jait Dixit committed
35
class FFTOperator(LinearOperator):
36 37 38
    """Transforms between a pair of position and harmonic domains.

    Built-in domain pairs are
39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
      - a harmonic and a non-harmonic RGSpace (with matching distances)
      - a HPSpace and a LMSpace
      - a GLSpace and a LMSpace
    Within a domain pair, both orderings are possible.

    The operator provides a "times" and an "adjoint_times" operation.
    For a pair of RGSpaces, the "adjoint_times" operation is equivalent to
    "inverse_times"; for the sphere-related domains this is not the case, since
    the operator matrix is not square.

    Parameters
    ----------
    domain: Space or single-element tuple of Spaces
        The domain of the data that is input by "times" and output by
        "adjoint_times".
54
    target: Space or single-element tuple of Spaces (optional)
55 56 57 58
        The domain of the data that is output by "times" and input by
        "adjoint_times".
        If omitted, a co-domain will be chosen automatically.
        Whenever "domain" is an RGSpace, the codomain (and its parameters) are
Martin Reinecke's avatar
Martin Reinecke committed
59
        uniquely determined.
60 61
        For GLSpace, HPSpace, and LMSpace, a sensible (but not unique)
        co-domain is chosen that should work satisfactorily in most situations,
Martin Reinecke's avatar
Martin Reinecke committed
62
        but for full control, the user should explicitly specify a codomain.
63 64 65 66 67 68 69 70 71 72

    Attributes
    ----------
    domain: Tuple of Spaces (with one entry)
        The domain of the data that is input by "times" and output by
        "adjoint_times".
    target: Tuple of Spaces (with one entry)
        The domain of the data that is output by "times" and input by
        "adjoint_times".
    unitary: bool
73 74
        Returns True if the operator is unitary (currently only the case if
        the domain and codomain are RGSpaces), else False.
75 76 77 78 79

    Raises
    ------
    ValueError:
        if "domain" or "target" are not of the proper type.
80

81
    """
82

83
    # ---Class attributes---
84

85 86 87
    default_codomain_dictionary = {RGSpace: RGSpace,
                                   HPSpace: LMSpace,
                                   GLSpace: LMSpace,
88
                                   LMSpace: GLSpace,
89 90 91 92 93 94 95 96 97
                                   }

    transformation_dictionary = {(RGSpace, RGSpace): RGRGTransformation,
                                 (HPSpace, LMSpace): HPLMTransformation,
                                 (GLSpace, LMSpace): GLLMTransformation,
                                 (LMSpace, HPSpace): LMHPTransformation,
                                 (LMSpace, GLSpace): LMGLTransformation
                                 }

Jait Dixit's avatar
Jait Dixit committed
98 99
    # ---Overwritten properties and methods---

Martin Reinecke's avatar
more  
Martin Reinecke committed
100
    def __init__(self, domain, target=None, default_spaces=None):
101
        super(FFTOperator, self).__init__(default_spaces)
102 103

        # Initialize domain and target
104
        self._domain = self._parse_domain(domain)
105
        if len(self.domain) != 1:
106 107
            raise ValueError("TransformationOperator accepts only exactly one "
                             "space as input domain.")
Jait Dixit's avatar
Jait Dixit committed
108

109
        if target is None:
110
            target = (self.get_default_codomain(self.domain[0]), )
Jait Dixit's avatar
Jait Dixit committed
111
        self._target = self._parse_domain(target)
112 113 114
        if len(self.target) != 1:
            raise ValueError("TransformationOperator accepts only exactly one "
                             "space as output target.")
Jait Dixit's avatar
Jait Dixit committed
115

116
        # Create transformation instances
117
        forward_class = self.transformation_dictionary[
118
                (self.domain[0].__class__, self.target[0].__class__)]
119
        backward_class = self.transformation_dictionary[
120 121
                (self.target[0].__class__, self.domain[0].__class__)]

122
        self._forward_transformation = forward_class(
Martin Reinecke's avatar
Martin Reinecke committed
123
            self.domain[0], self.target[0])
124

125
        self._backward_transformation = backward_class(
Martin Reinecke's avatar
Martin Reinecke committed
126
            self.target[0], self.domain[0])
Jait Dixit's avatar
Jait Dixit committed
127

128
    def _times(self, x, spaces):
129
        spaces = utilities.cast_axis_to_tuple(spaces, len(x.domain))
130
        if spaces is None:
131 132 133 134
            # this case means that x lives on only one space, which is
            # identical to the space in the domain of `self`. Otherwise the
            # input check of LinearOperator would have failed.
            axes = x.domain_axes[0]
135 136
        else:
            axes = x.domain_axes[spaces[0]]
137

138
        new_val = self._forward_transformation.transform(x.val, axes=axes)
139

140 141 142 143 144
        if spaces is None:
            result_domain = self.target
        else:
            result_domain = list(x.domain)
            result_domain[spaces[0]] = self.target[0]
145

Martin Reinecke's avatar
more  
Martin Reinecke committed
146
        result_dtype = new_val.dtype
147
        result_field = x.copy_empty(domain=result_domain,
148
                                    dtype=result_dtype)
Theo Steininger's avatar
Theo Steininger committed
149
        result_field.set_val(new_val=new_val, copy=True)
Jait Dixit's avatar
Jait Dixit committed
150

151
        return result_field
Jait Dixit's avatar
Jait Dixit committed
152

153
    def _adjoint_times(self, x, spaces):
Jait Dixit's avatar
Jait Dixit committed
154
        spaces = utilities.cast_axis_to_tuple(spaces, len(x.domain))
155
        if spaces is None:
156 157 158 159
            # this case means that x lives on only one space, which is
            # identical to the space in the domain of `self`. Otherwise the
            # input check of LinearOperator would have failed.
            axes = x.domain_axes[0]
160 161
        else:
            axes = x.domain_axes[spaces[0]]
Jait Dixit's avatar
Jait Dixit committed
162

163
        new_val = self._backward_transformation.transform(x.val, axes=axes)
164 165 166 167 168 169 170

        if spaces is None:
            result_domain = self.domain
        else:
            result_domain = list(x.domain)
            result_domain[spaces[0]] = self.domain[0]

Martin Reinecke's avatar
more  
Martin Reinecke committed
171
        result_dtype = new_val.dtype
172

173
        result_field = x.copy_empty(domain=result_domain,
174
                                    dtype=result_dtype)
Theo Steininger's avatar
Theo Steininger committed
175
        result_field.set_val(new_val=new_val, copy=True)
176 177

        return result_field
Jait Dixit's avatar
Jait Dixit committed
178 179 180

    # ---Mandatory properties and methods---

181 182 183 184
    @property
    def domain(self):
        return self._domain

Jait Dixit's avatar
Jait Dixit committed
185 186 187 188
    @property
    def target(self):
        return self._target

189 190
    @property
    def unitary(self):
191 192
        return (self._forward_transformation.unitary and
                self._backward_transformation.unitary)
193 194 195 196 197

    # ---Added properties and methods---

    @classmethod
    def get_default_codomain(cls, domain):
198
        """Returns a codomain to the given domain.
199 200 201 202 203 204 205 206 207 208 209 210

        Parameters
        ----------
        domain: Space
            An instance of RGSpace, HPSpace, GLSpace or LMSpace.

        Returns
        -------
        target: Space
            A (more or less perfect) counterpart to "domain" with respect
            to a FFT operation.
            Whenever "domain" is an RGSpace, the codomain (and its parameters)
Martin Reinecke's avatar
Martin Reinecke committed
211
            are uniquely determined.
Martin Reinecke's avatar
Martin Reinecke committed
212 213 214 215
            For GLSpace, HPSpace, and LMSpace, a sensible (but not unique)
            co-domain is chosen that should work satisfactorily in most
            situations. For full control however, the user should not rely on
            this method.
216 217 218 219 220

        Raises
        ------
        ValueError:
            if no default codomain is defined for "domain".
221

222
        """
223 224 225 226
        domain_class = domain.__class__
        try:
            codomain_class = cls.default_codomain_dictionary[domain_class]
        except KeyError:
227
            raise ValueError("Unknown domain")
228 229 230 231 232

        try:
            transform_class = cls.transformation_dictionary[(domain_class,
                                                             codomain_class)]
        except KeyError:
233
            raise ValueError(
234
                "No transformation for domain-codomain pair found.")
235 236

        return transform_class.get_codomain(domain)