diagonal_operator.py 12.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# NIFTy
# Copyright (C) 2017  Theo Steininger
#
# Author: Theo Steininger
#
# 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/>.
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

    def diagonal(self, bare=False, copy=True):
Theo Steininger's avatar
Theo Steininger committed
136
        """ Returns the diagonal of the Operator.
137
138
139
140
141
142
143
144
145
146

        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
        -------
Theo Steininger's avatar
Theo Steininger committed
147
148
        out : Field
            The diagonal of the Operator.
149
150

        """
151
152
153
154
155
156
157
158
159
        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):
160
161
162
163
164
165
166
167
168
        """ Returns the inverse-diagonal of the operator.

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

        Returns
        -------
Theo Steininger's avatar
Theo Steininger committed
169
170
        out : Field
            The inverse of the diagonal of the Operator.
171

172
        """        
173
        return 1./self.diagonal(bare=bare, copy=False)
174
175

    def trace(self, bare=False):
176
177
178
179
180
181
182
183
184
185
        """ Returns the trace the operator.

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

        Returns
        -------
        out : scalar
Theo Steininger's avatar
Theo Steininger committed
186
            The trace of the Operator.
187
188

        """
189
190
191
        return self.diagonal(bare=bare, copy=False).sum()

    def inverse_trace(self, bare=False):
192
193
194
195
196
197
198
199
200
201
        """ Returns the inverse-trace of the operator.

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

        Returns
        -------
        out : scalar
Theo Steininger's avatar
Theo Steininger committed
202
            The inverse of the trace of the Operator.
203
204

        """
205
        return self.inverse_diagonal(bare=bare).sum()
206
207

    def trace_log(self):
208
209
210
211
212
        """ Returns the trave-log of the operator.

        Returns
        -------
        out : scalar
Theo Steininger's avatar
Theo Steininger committed
213
            the trace of the logarithm of the Operator.
214
215

        """
216
        log_diagonal = nifty_log(self.diagonal(copy=False))
217
218
219
        return log_diagonal.sum()

    def determinant(self):
220
221
222
223
        """ Returns the determinant of the operator.

        Returns
        -------
Pumpe, Daniel (dpumpe)'s avatar
update    
Pumpe, Daniel (dpumpe) committed
224
        out : scalar
225
226
227
228
        out : scalar
            the determinant of the Operator

        """
Theo Steininger's avatar
Theo Steininger committed
229

230
231
232
        return self.diagonal(copy=False).val.prod()

    def inverse_determinant(self):
233
234
235
236
237
238
239
240
        """ Returns the inverse-determinant of the operator.

        Returns
        -------
        out : scalar
            the inverse-determinant of the Operator

        """
Theo Steininger's avatar
Theo Steininger committed
241

242
243
244
        return 1/self.determinant()

    def log_determinant(self):
245
246
247
248
249
250
251
252
        """ Returns the log-eterminant of the operator.

        Returns
        -------
        out : scalar
            the log-determinant of the Operator

        """
Theo Steininger's avatar
Theo Steininger committed
253

254
255
256
257
        return np.log(self.determinant())

    # ---Mandatory properties and methods---

258
259
260
261
    @property
    def domain(self):
        return self._domain

262
    @property
Martin Reinecke's avatar
Martin Reinecke committed
263
264
265
266
    def self_adjoint(self):
        if self._self_adjoint is None:
            self._self_adjoint = (self._diagonal.val.imag == 0).all()
        return self._self_adjoint
267
268
269

    @property
    def unitary(self):
270
271
272
        if self._unitary is None:
            self._unitary = (self._diagonal.val *
                             self._diagonal.val.conjugate() == 1).all()
273
274
275
276
277
        return self._unitary

    # ---Added properties and methods---

    @property
278
    def distribution_strategy(self):
279
280
281
        """
        distribution_strategy : string
            Defines the way how the diagonal operator is distributed
Theo Steininger's avatar
Theo Steininger committed
282
283
            among the nodes. Available distribution_strategies are:
            'fftw', 'equal' and 'not'.
284
285
286

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

288
        """
Theo Steininger's avatar
Theo Steininger committed
289

290
        return self._distribution_strategy
291

292
293
    def _parse_distribution_strategy(self, distribution_strategy, val):
        if distribution_strategy is None:
294
            if isinstance(val, distributed_data_object):
295
                distribution_strategy = val.distribution_strategy
296
            elif isinstance(val, Field):
297
                distribution_strategy = val.distribution_strategy
298
            else:
299
                self.logger.info("Datamodel set to default!")
300
301
                distribution_strategy = gc['default_distribution_strategy']
        elif distribution_strategy not in DISTRIBUTION_STRATEGIES['all']:
302
303
            raise ValueError(
                    "Invalid distribution_strategy!")
304
        return distribution_strategy
305
306

    def set_diagonal(self, diagonal, bare=False, copy=True):
307
308
309
310
        """ Sets the diagonal of the Operator.

        Parameters
        ----------
Theo Steininger's avatar
Theo Steininger committed
311
        diagonal : {scalar, list, array, Field, d2o-object}
312
313
            The diagonal entries of the operator.
        bare : boolean
Theo Steininger's avatar
Theo Steininger committed
314
315
            Indicates whether the input for the diagonal is bare or not
            (default: False).
316
        copy : boolean
Theo Steininger's avatar
Theo Steininger committed
317
            Specifies if a copy of the input shall be made (default: True).
318
319
320

        """

321
322
323
        # use the casting functionality from Field to process `diagonal`
        f = Field(domain=self.domain,
                  val=diagonal,
324
                  distribution_strategy=self.distribution_strategy,
325
326
                  copy=copy)

327
        # weight if the given values were `bare` is True
328
        # do inverse weightening if the other way around
329
        if bare:
330
331
332
            # If `copy` is True, we won't change external data by weightening
            # Otherwise, inplace weightening would change the external field
            f.weight(inplace=copy)
333

Martin Reinecke's avatar
Martin Reinecke committed
334
335
        # Reset the self_adjoint property:
        self._self_adjoint = None
336

337
338
        # Reset the unitarity property
        self._unitary = None
339
340
341

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

343
344
    def _times_helper(self, x, spaces, operation):
        # if the domain matches directly
345
        # -> multiply the fields directly
346
        if x.domain == self.domain:
347
348
349
350
351
352
353
            # 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:
354
            active_axes = range(len(x.shape))
355
356
357
358
359
360
361
362
363
364
        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
365
366
            self.logger.warn("The input field is not sub-slice compatible to "
                             "the distribution strategy of the operator.")
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
            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