diagonal_operator.py 10.8 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
21
22
23

import numpy as np

from d2o import distributed_data_object,\
                STRATEGIES as DISTRIBUTION_STRATEGIES

24
from nifty.basic_arithmetics import log as nifty_log
25
from nifty.config import nifty_configuration as gc
26
27
28
29
30
from nifty.field import Field
from nifty.operators.endomorphic_operator import EndomorphicOperator


class DiagonalOperator(EndomorphicOperator):
Theo Steininger's avatar
Theo Steininger committed
31
32
33
34
35
    """ 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.
36

37
38
39

    Parameters
    ----------
Theo Steininger's avatar
Theo Steininger committed
40
41
42
    domain : tuple of DomainObjects, i.e. Spaces and FieldTypes
        The domain on which the Operator's input Field lives.
    diagonal : {scalar, list, array, Field, d2o-object}
43
44
        The diagonal entries of the operator.
    bare : boolean
Theo Steininger's avatar
Theo Steininger committed
45
46
        Indicates whether the input for the diagonal is bare or not
        (default: False).
47
48
49
50
51
    copy : boolean
        Internal copy of the diagonal (default: True)
    distribution_strategy : string
        setting the prober distribution_strategy of the
        diagonal (default : None). In case diagonal is d2o-object or Field,
Theo Steininger's avatar
Theo Steininger committed
52
        their distribution_strategy is used as a fallback.
53
54
55
    default_spaces : tuple of ints *optional*
        Defines on which space(s) of a given field the Operator acts by
        default (default: None)
56
57
58

    Attributes
    ----------
59
60
61
62
63
64
65
66
67
    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.
68
    distribution_strategy : string
Theo Steininger's avatar
Theo Steininger committed
69
70
        Defines the distribution_strategy of the distributed_data_object
        in which the diagonal entries are stored in.
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87

    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.

    Examples
    --------
    >>> x_space = RGSpace(5)
Theo Steininger's avatar
Theo Steininger committed
88
89
    >>> D = DiagonalOperator(x_space, diagonal=[1., 3., 2., 4., 6.])
    >>> f = Field(x_space, val=2.)
90
91
92
    >>> res = D.times(f)
    >>> res.val
    <distributed_data_object>
Theo Steininger's avatar
Theo Steininger committed
93
    array([ 2.,  6.,  4.,  8.,  12.])
94
95
96
97
98
99
100

    See Also
    --------
    EndomorphicOperator

    """

101
102
    # ---Overwritten properties and methods---

103
104
105
106
    def __init__(self, domain=(), diagonal=None, bare=False, copy=True,
                 distribution_strategy=None, default_spaces=None):
        super(DiagonalOperator, self).__init__(default_spaces)

107
        self._domain = self._parse_domain(domain)
108

109
        if distribution_strategy is None:
110
            if isinstance(diagonal, distributed_data_object):
111
                distribution_strategy = diagonal.distribution_strategy
112
            elif isinstance(diagonal, Field):
113
                distribution_strategy = diagonal.distribution_strategy
114

115
        self._distribution_strategy = self._parse_distribution_strategy(
116
117
                               distribution_strategy=distribution_strategy,
                               val=diagonal)
118
119
120

        self.set_diagonal(diagonal=diagonal, bare=bare, copy=copy)

121
122
    def _times(self, x, spaces):
        return self._times_helper(x, spaces, operation=lambda z: z.__mul__)
123

124
125
    def _adjoint_times(self, x, spaces):
        return self._times_helper(x, spaces,
126
                                  operation=lambda z: z.adjoint().__mul__)
127

128
129
    def _inverse_times(self, x, spaces):
        return self._times_helper(x, spaces, operation=lambda z: z.__rdiv__)
130

131
132
    def _adjoint_inverse_times(self, x, spaces):
        return self._times_helper(x, spaces,
133
                                  operation=lambda z: z.adjoint().__rdiv__)
134

135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
    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:
            diagonal = self._diagonal.weight(power=-1)
        elif copy:
            diagonal = self._diagonal.copy()
        else:
            diagonal = self._diagonal
        return diagonal

    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)

175
176
    # ---Mandatory properties and methods---

177
178
179
180
    @property
    def domain(self):
        return self._domain

181
    @property
Martin Reinecke's avatar
Martin Reinecke committed
182
183
184
185
    def self_adjoint(self):
        if self._self_adjoint is None:
            self._self_adjoint = (self._diagonal.val.imag == 0).all()
        return self._self_adjoint
186
187
188

    @property
    def unitary(self):
189
190
191
        if self._unitary is None:
            self._unitary = (self._diagonal.val *
                             self._diagonal.val.conjugate() == 1).all()
192
193
194
195
196
        return self._unitary

    # ---Added properties and methods---

    @property
197
    def distribution_strategy(self):
198
199
200
        """
        distribution_strategy : string
            Defines the way how the diagonal operator is distributed
Theo Steininger's avatar
Theo Steininger committed
201
202
            among the nodes. Available distribution_strategies are:
            'fftw', 'equal' and 'not'.
203
204
205

        Notes :
            https://arxiv.org/abs/1606.05385
Theo Steininger's avatar
Theo Steininger committed
206

207
        """
Theo Steininger's avatar
Theo Steininger committed
208

209
        return self._distribution_strategy
210

211
212
    def _parse_distribution_strategy(self, distribution_strategy, val):
        if distribution_strategy is None:
213
            if isinstance(val, distributed_data_object):
214
                distribution_strategy = val.distribution_strategy
215
            elif isinstance(val, Field):
216
                distribution_strategy = val.distribution_strategy
217
            else:
218
                self.logger.info("Datamodel set to default!")
219
220
                distribution_strategy = gc['default_distribution_strategy']
        elif distribution_strategy not in DISTRIBUTION_STRATEGIES['all']:
221
222
            raise ValueError(
                    "Invalid distribution_strategy!")
223
        return distribution_strategy
224
225

    def set_diagonal(self, diagonal, bare=False, copy=True):
226
227
228
229
        """ Sets the diagonal of the Operator.

        Parameters
        ----------
Theo Steininger's avatar
Theo Steininger committed
230
        diagonal : {scalar, list, array, Field, d2o-object}
231
232
            The diagonal entries of the operator.
        bare : boolean
Theo Steininger's avatar
Theo Steininger committed
233
234
            Indicates whether the input for the diagonal is bare or not
            (default: False).
235
        copy : boolean
Theo Steininger's avatar
Theo Steininger committed
236
            Specifies if a copy of the input shall be made (default: True).
237
238
239

        """

240
241
242
        # use the casting functionality from Field to process `diagonal`
        f = Field(domain=self.domain,
                  val=diagonal,
243
                  distribution_strategy=self.distribution_strategy,
244
245
                  copy=copy)

246
        # weight if the given values were `bare` is True
247
        # do inverse weightening if the other way around
248
        if bare:
249
250
251
            # If `copy` is True, we won't change external data by weightening
            # Otherwise, inplace weightening would change the external field
            f.weight(inplace=copy)
252

Martin Reinecke's avatar
Martin Reinecke committed
253
254
        # Reset the self_adjoint property:
        self._self_adjoint = None
255

256
257
        # Reset the unitarity property
        self._unitary = None
258
259
260

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

262
263
    def _times_helper(self, x, spaces, operation):
        # if the domain matches directly
264
        # -> multiply the fields directly
265
        if x.domain == self.domain:
266
267
268
269
270
271
272
            # here the actual multiplication takes place
            return operation(self.diagonal(copy=False))(x)

        # if the distribution_strategy of self is sub-slice compatible to
        # the one of x, reshape the local data of self and apply it directly
        active_axes = []
        if spaces is None:
273
            active_axes = range(len(x.shape))
274
275
276
277
278
279
280
281
282
283
        else:
            for space_index in spaces:
                active_axes += x.domain_axes[space_index]

        axes_local_distribution_strategy = \
            x.val.get_axes_local_distribution_strategy(active_axes)
        if axes_local_distribution_strategy == self.distribution_strategy:
            local_diagonal = self._diagonal.val.get_local_data(copy=False)
        else:
            # create an array that is sub-slice compatible
284
285
            self.logger.warn("The input field is not sub-slice compatible to "
                             "the distribution strategy of the operator.")
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
            redistr_diagonal_val = self._diagonal.val.copy(
                distribution_strategy=axes_local_distribution_strategy)
            local_diagonal = redistr_diagonal_val.get_local_data(copy=False)

        reshaper = [x.shape[i] if i in active_axes else 1
                    for i in xrange(len(x.shape))]
        reshaped_local_diagonal = np.reshape(local_diagonal, reshaper)

        # here the actual multiplication takes place
        local_result = operation(reshaped_local_diagonal)(
                           x.val.get_local_data(copy=False))

        result_field = x.copy_empty(dtype=local_result.dtype)
        result_field.val.set_local_data(local_result, copy=False)
        return result_field