diagonal_operator.py 7.39 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

Martin Reinecke's avatar
Martin Reinecke committed
19
20
from __future__ import division
from builtins import range
21
22
import numpy as np

Martin Reinecke's avatar
Martin Reinecke committed
23
24
from ...field import Field
from ..endomorphic_operator import EndomorphicOperator
25
26
27


class DiagonalOperator(EndomorphicOperator):
Theo Steininger's avatar
Theo Steininger committed
28
29
30
31
32
    """ NIFTY class for diagonal operators.

    The NIFTY DiagonalOperator class is a subclass derived from the
    EndomorphicOperator. It multiplies an input field pixel-wise with its
    diagonal.
33

34
35
36

    Parameters
    ----------
Theo Steininger's avatar
Theo Steininger committed
37
38
    domain : tuple of DomainObjects, i.e. Spaces and FieldTypes
        The domain on which the Operator's input Field lives.
Martin Reinecke's avatar
Martin Reinecke committed
39
    diagonal : {scalar, list, array, Field}
40
41
        The diagonal entries of the operator.
    bare : boolean
Theo Steininger's avatar
Theo Steininger committed
42
43
        Indicates whether the input for the diagonal is bare or not
        (default: False).
44
45
    copy : boolean
        Internal copy of the diagonal (default: True)
46
47
48
    default_spaces : tuple of ints *optional*
        Defines on which space(s) of a given field the Operator acts by
        default (default: None)
49
50
51

    Attributes
    ----------
52
53
54
55
56
57
58
59
60
    domain : tuple of DomainObjects, i.e. Spaces and FieldTypes
        The domain on which the Operator's input Field lives.
    target : tuple of DomainObjects, i.e. Spaces and FieldTypes
        The domain in which the outcome of the operator lives. As the Operator
        is endomorphic this is the same as its domain.
    unitary : boolean
        Indicates whether the Operator is unitary or not.
    self_adjoint : boolean
        Indicates whether the operator is self_adjoint or not.
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80

    Raises
    ------

    Notes
    -----
    The ambiguity of bare or non-bare diagonal entries is based on the choice
    of a matrix representation of the operator in question. The naive choice
    of absorbing the volume weights into the matrix leads to a matrix-vector
    calculus with the non-bare entries which seems intuitive, though.
    The choice of keeping matrix entries and volume weights separate
    deals with the bare entries that allow for correct interpretation
    of the matrix entries; e.g., as variance in case of an covariance operator.

    See Also
    --------
    EndomorphicOperator

    """

81
82
    # ---Overwritten properties and methods---

83
    def __init__(self, domain=(), diagonal=None, bare=False, copy=True,
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
84
                 default_spaces=None):
85
86
        super(DiagonalOperator, self).__init__(default_spaces)

87
        self._domain = self._parse_domain(domain)
88

89
90
        self._self_adjoint = None
        self._unitary = None
91
92
        self.set_diagonal(diagonal=diagonal, bare=bare, copy=copy)

93
94
    def _times(self, x, spaces):
        return self._times_helper(x, spaces, operation=lambda z: z.__mul__)
95

96
97
    def _adjoint_times(self, x, spaces):
        return self._times_helper(x, spaces,
98
                                  operation=lambda z: z.adjoint().__mul__)
99

100
    def _inverse_times(self, x, spaces):
Martin Reinecke's avatar
Martin Reinecke committed
101
        return self._times_helper(x, spaces, operation=lambda z: z.__rtruediv__)
102

103
104
    def _adjoint_inverse_times(self, x, spaces):
        return self._times_helper(x, spaces,
Martin Reinecke's avatar
Martin Reinecke committed
105
                                  operation=lambda z: z.adjoint().__rtruediv__)
106

107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
    def diagonal(self, bare=False, copy=True):
        """ Returns the diagonal of the Operator.

        Parameters
        ----------
        bare : boolean
            Whether the returned Field values should be bare or not.
        copy : boolean
            Whether the returned Field should be copied or not.

        Returns
        -------
        out : Field
            The diagonal of the Operator.

        """
        if bare:
Martin Reinecke's avatar
Martin Reinecke committed
124
            return self._diagonal.weight(power=-1)
125
        elif copy:
Martin Reinecke's avatar
Martin Reinecke committed
126
            return self._diagonal.copy()
127
        else:
Martin Reinecke's avatar
Martin Reinecke committed
128
            return self._diagonal
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145

    def inverse_diagonal(self, bare=False):
        """ Returns the inverse-diagonal of the operator.

        Parameters
        ----------
        bare : boolean
            Whether the returned Field values should be bare or not.

        Returns
        -------
        out : Field
            The inverse of the diagonal of the Operator.

        """
        return 1./self.diagonal(bare=bare, copy=False)

146
147
    # ---Mandatory properties and methods---

148
149
150
151
    @property
    def domain(self):
        return self._domain

152
    @property
Martin Reinecke's avatar
Martin Reinecke committed
153
154
155
156
    def self_adjoint(self):
        if self._self_adjoint is None:
            self._self_adjoint = (self._diagonal.val.imag == 0).all()
        return self._self_adjoint
157
158
159

    @property
    def unitary(self):
160
161
162
        if self._unitary is None:
            self._unitary = (self._diagonal.val *
                             self._diagonal.val.conjugate() == 1).all()
163
164
165
166
167
        return self._unitary

    # ---Added properties and methods---

    def set_diagonal(self, diagonal, bare=False, copy=True):
168
169
170
171
        """ Sets the diagonal of the Operator.

        Parameters
        ----------
Martin Reinecke's avatar
Martin Reinecke committed
172
        diagonal : {scalar, list, array, Field}
173
174
            The diagonal entries of the operator.
        bare : boolean
Theo Steininger's avatar
Theo Steininger committed
175
176
            Indicates whether the input for the diagonal is bare or not
            (default: False).
177
        copy : boolean
Theo Steininger's avatar
Theo Steininger committed
178
            Specifies if a copy of the input shall be made (default: True).
179
180
181

        """

182
        # use the casting functionality from Field to process `diagonal`
Martin Reinecke's avatar
Martin Reinecke committed
183
        f = Field(domain=self.domain, val=diagonal, copy=copy)
184

185
        # weight if the given values were `bare` is True
186
        # do inverse weightening if the other way around
187
        if bare:
188
189
190
            # If `copy` is True, we won't change external data by weightening
            # Otherwise, inplace weightening would change the external field
            f.weight(inplace=copy)
191

Martin Reinecke's avatar
Martin Reinecke committed
192
193
        # Reset the self_adjoint property:
        self._self_adjoint = None
194

195
196
        # Reset the unitarity property
        self._unitary = None
197
198
199

        # store the diagonal-field
        self._diagonal = f
200

201
202
    def _times_helper(self, x, spaces, operation):
        # if the domain matches directly
203
        # -> multiply the fields directly
204
        if x.domain == self.domain:
205
206
207
208
209
            # here the actual multiplication takes place
            return operation(self.diagonal(copy=False))(x)

        active_axes = []
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
210
            active_axes = list(range(len(x.shape)))
211
212
213
214
        else:
            for space_index in spaces:
                active_axes += x.domain_axes[space_index]

Martin Reinecke's avatar
Martin Reinecke committed
215
        local_diagonal = self._diagonal.val
216

Martin Reinecke's avatar
Martin Reinecke committed
217
        reshaper = [x.val.shape[i] if i in active_axes else 1
Martin Reinecke's avatar
Martin Reinecke committed
218
                    for i in range(len(x.shape))]
219
220
221
        reshaped_local_diagonal = np.reshape(local_diagonal, reshaper)

        # here the actual multiplication takes place
Martin Reinecke's avatar
Martin Reinecke committed
222
        local_result = operation(reshaped_local_diagonal)(x.val)
223

Martin Reinecke's avatar
Martin Reinecke committed
224
        return Field (x.domain,val=local_result)