field.py 26.7 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

csongor's avatar
csongor committed
19
from __future__ import division
Martin Reinecke's avatar
Martin Reinecke committed
20
from builtins import range
21

csongor's avatar
csongor committed
22
23
import numpy as np

Martin Reinecke's avatar
Martin Reinecke committed
24
from .domain_object import DomainObject
25

Martin Reinecke's avatar
Martin Reinecke committed
26
from .spaces.power_space import PowerSpace
csongor's avatar
csongor committed
27

Martin Reinecke's avatar
Martin Reinecke committed
28
29
from . import nifty_utilities as utilities
from .random import Random
Martin Reinecke's avatar
Martin Reinecke committed
30
from functools import reduce
31

csongor's avatar
csongor committed
32

Martin Reinecke's avatar
Martin Reinecke committed
33
class Field(object):
Theo Steininger's avatar
Theo Steininger committed
34
35
36
    """ The discrete representation of a continuous field over multiple spaces.

    In NIFTY, Fields are used to store data arrays and carry all the needed
37
    metainformation (i.e. the domain) for operators to be able to work on them.
Martin Reinecke's avatar
updates    
Martin Reinecke committed
38
    In addition, Field has methods to work with power spectra.
Theo Steininger's avatar
Theo Steininger committed
39

40
41
42
43
    Parameters
    ----------
    domain : DomainObject
        One of the space types NIFTY supports. RGSpace, GLSpace, HPSpace,
Theo Steininger's avatar
Theo Steininger committed
44
        LMSpace or PowerSpace. It might also be a FieldArray, which is
45
        an unstructured domain.
Theo Steininger's avatar
Theo Steininger committed
46

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
47
    val : scalar, numpy.ndarray, Field
48
49
50
        The values the array should contain after init. A scalar input will
        fill the whole array with this scalar. If an array is provided the
        array's dimensions must match the domain's.
Theo Steininger's avatar
Theo Steininger committed
51

52
    dtype : type
Martin Reinecke's avatar
updates    
Martin Reinecke committed
53
        A numpy.type. Most common are float and complex.
Theo Steininger's avatar
Theo Steininger committed
54

55
56
57
58
    copy: boolean

    Attributes
    ----------
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
59
    val : numpy.ndarray
Theo Steininger's avatar
Theo Steininger committed
60

61
62
63
64
65
66
    domain : DomainObject
        See Parameters.
    domain_axes : tuple of tuples
        Enumerates the axes of the Field
    dtype : type
        Contains the datatype stored in the Field.
Theo Steininger's avatar
Theo Steininger committed
67

68
69
70
71
72
73
74
    Raise
    -----
    TypeError
        Raised if
            *the given domain contains something that is not a DomainObject
             instance
            *val is an array that has a different dimension than the domain
Theo Steininger's avatar
Theo Steininger committed
75

76
    """
77

theos's avatar
theos committed
78
    # ---Initialization methods---
79

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
80
    def __init__(self, domain=None, val=None, dtype=None, copy=False):
81
        self.domain = self._parse_domain(domain=domain, val=val)
82
        self.domain_axes = self._get_axes_tuple(self.domain)
83

Martin Reinecke's avatar
Martin Reinecke committed
84
        shape_tuple = tuple(sp.shape for sp in self.domain)
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
85
        if len(shape_tuple) == 0:
Martin Reinecke's avatar
Martin Reinecke committed
86
            global_shape = ()
87
        else:
Martin Reinecke's avatar
Martin Reinecke committed
88
89
            global_shape = reduce(lambda x, y: x + y, shape_tuple)
        dtype = self._infer_dtype(dtype=dtype, val=val)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
90
        if isinstance(val, Field):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
91
            if self.domain != val.domain:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
92
                raise ValueError("Domain mismatch")
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
93
            self._val = np.array(val.val, dtype=dtype, copy=copy)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
94
        elif (np.isscalar(val)):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
95
            self._val = np.full(global_shape, dtype=dtype, fill_value=val)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
96
        elif isinstance(val, np.ndarray):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
97
98
            if global_shape == val.shape:
                self._val = np.array(val, dtype=dtype, copy=copy)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
99
100
101
            else:
                raise ValueError("Shape mismatch")
        elif val is None:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
102
            self._val = np.empty(global_shape, dtype=dtype)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
103
104
        else:
            raise TypeError("unknown source type")
csongor's avatar
csongor committed
105

106
    def _parse_domain(self, domain, val=None):
107
        if domain is None:
108
            if isinstance(val, Field):
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
109
110
111
112
113
114
115
116
                return val.domain
            if np.isscalar(val):
                return ()  # empty domain tuple
            raise TypeError("could not infer domain from value")
        if isinstance(domain, DomainObject):
            return (domain,)

        if not isinstance(domain, tuple):
117
            domain = tuple(domain)
csongor's avatar
csongor committed
118
        for d in domain:
119
            if not isinstance(d, DomainObject):
120
121
                raise TypeError(
                    "Given domain contains something that is not a "
122
                    "DomainObject instance.")
csongor's avatar
csongor committed
123
124
        return domain

125
126
    def _get_axes_tuple(self, things_with_shape):
        i = 0
theos's avatar
theos committed
127
128
        axes_list = []
        for thing in things_with_shape:
129
            nax = len(thing.shape)
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
130
            axes_list += [tuple(range(i, i+nax))]
131
            i += nax
theos's avatar
theos committed
132
        return tuple(axes_list)
133

134
    def _infer_dtype(self, dtype, val):
Martin Reinecke's avatar
Martin Reinecke committed
135
136
137
138
139
        if val is None or dtype is not None:
            return np.result_type(dtype, np.float64)
        if isinstance(val, Field):
            return val.dtype
        return np.result_type(val, np.float64)
140

141
    # ---Factory methods---
142

143
    @classmethod
Martin Reinecke's avatar
Martin Reinecke committed
144
    def from_random(cls, random_type, domain, dtype=None, **kwargs):
145
146
147
148
149
        """ Draws a random field with the given parameters.

        Parameters
        ----------
        cls : class
Theo Steininger's avatar
Theo Steininger committed
150

151
152
153
        random_type : String
            'pm1', 'normal', 'uniform' are the supported arguments for this
            method.
Theo Steininger's avatar
Theo Steininger committed
154

155
156
        domain : DomainObject
            The domain of the output random field
Theo Steininger's avatar
Theo Steininger committed
157

158
159
        dtype : type
            The datatype of the output random field
Theo Steininger's avatar
Theo Steininger committed
160

161
162
163
164
165
166
167
        Returns
        -------
        out : Field
            The output object.

        See Also
        --------
168
        power_synthesize
169
        """
Theo Steininger's avatar
Theo Steininger committed
170

171
        generator_function = getattr(Random, random_type)
Martin Reinecke's avatar
Martin Reinecke committed
172
173
174
        return Field(domain=domain,
                     val=generator_function(dtype=np.dtype(dtype),
                     shape=utilities.domains2shape(domain), **kwargs))
175

176
177
    # ---Powerspectral methods---

Martin Reinecke's avatar
Martin Reinecke committed
178
179
    def power_analyze(self, spaces=None, binbounds=None,
                      keep_phase_information=False):
Theo Steininger's avatar
Theo Steininger committed
180
        """ Computes the square root power spectrum for a subspace of `self`.
Theo Steininger's avatar
Theo Steininger committed
181

Theo Steininger's avatar
Theo Steininger committed
182
183
184
        Creates a PowerSpace for the space addressed by `spaces` with the given
        binning and computes the power spectrum as a Field over this
        PowerSpace. This can only be done if the subspace to  be analyzed is a
185
        harmonic space. The resulting field has the same units as the initial
Theo Steininger's avatar
Theo Steininger committed
186
        field, corresponding to the square root of the power spectrum.
187
188
189

        Parameters
        ----------
Theo Steininger's avatar
Theo Steininger committed
190
191
192
193
194
        spaces : int *optional*
            The subspace for which the powerspectrum shall be computed
            (default : None).
        binbounds : array-like *optional*
            Inner bounds of the bins (default : None).
Martin Reinecke's avatar
Martin Reinecke committed
195
            if binbounds==None : bins are inferred.
196
197
198
199
200
201
202
203
204
205
        keep_phase_information : boolean, *optional*
            If False, return a real-valued result containing the power spectrum
            of the input Field.
            If True, return a complex-valued result whose real component
            contains the power spectrum computed from the real part of the
            input Field, and whose imaginary component contains the power
            spectrum computed from the imaginary part of the input Field.
            The absolute value of this result should be identical to the output
            of power_analyze with keep_phase_information=False.
            (default : False).
Theo Steininger's avatar
Theo Steininger committed
206

207
208
209
210
        Raise
        -----
        ValueError
            Raised if
Theo Steininger's avatar
Theo Steininger committed
211
212
                *len(domain) is != 1 when spaces==None
                *len(spaces) is != 1 if not None
213
                *the analyzed space is not harmonic
Theo Steininger's avatar
Theo Steininger committed
214

215
216
        Returns
        -------
Theo Steininger's avatar
Theo Steininger committed
217
        out : Field
Martin Reinecke's avatar
Martin Reinecke committed
218
            The output object. Its domain is a PowerSpace and it contains
219
220
221
222
223
            the power spectrum of 'self's field.

        See Also
        --------
        power_synthesize, PowerSpace
Theo Steininger's avatar
Theo Steininger committed
224

225
        """
Theo Steininger's avatar
Theo Steininger committed
226

Theo Steininger's avatar
Theo Steininger committed
227
        # check if all spaces in `self.domain` are either harmonic or
228
229
230
        # power_space instances
        for sp in self.domain:
            if not sp.harmonic and not isinstance(sp, PowerSpace):
Martin Reinecke's avatar
Martin Reinecke committed
231
232
                raise TypeError("Field has a space in `domain` which is "
                                "neither harmonic nor a PowerSpace.")
233
234

        # check if the `spaces` input is valid
235
236
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
237
            spaces = list(range(len(self.domain)))
238
239

        if len(spaces) == 0:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
240
            raise ValueError("No space for analysis specified.")
241

242
243
244
245
        if keep_phase_information:
            parts_val = self._hermitian_decomposition(
                                              val=self.val,
                                              preserve_gaussian_variance=False)
Martin Reinecke's avatar
Martin Reinecke committed
246
            parts = [Field(self.domain, part_val, self.dtype, copy=False)
247
248
249
250
251
                     for part_val in parts_val]
        else:
            parts = [self]

        parts = [abs(part)**2 for part in parts]
252
253

        for space_index in spaces:
254
255
            parts = [self._single_power_analyze(
                                work_field=part,
256
                                space_index=space_index,
257
258
                                binbounds=binbounds)
                     for part in parts]
259

260
        if keep_phase_information:
Martin Reinecke's avatar
updates    
Martin Reinecke committed
261
            return parts[0] + 1j*parts[1]
262
        else:
Martin Reinecke's avatar
updates    
Martin Reinecke committed
263
            return parts[0]
264
265

    @classmethod
Martin Reinecke's avatar
Martin Reinecke committed
266
    def _single_power_analyze(cls, work_field, space_index, binbounds):
267
        if not work_field.domain[space_index].harmonic:
Martin Reinecke's avatar
Martin Reinecke committed
268
            raise ValueError("The analyzed space must be harmonic.")
269

270
271
272
273
274
275
        # Create the target PowerSpace instance:
        # If the associated signal-space field was real, we extract the
        # hermitian and anti-hermitian parts of `self` and put them
        # into the real and imaginary parts of the power spectrum.
        # If it was complex, all the power is put into a real power spectrum.

276
        harmonic_domain = work_field.domain[space_index]
277
        power_domain = PowerSpace(harmonic_partner=harmonic_domain,
Theo Steininger's avatar
Theo Steininger committed
278
                                  binbounds=binbounds)
279
280
        power_spectrum = cls._calculate_power_spectrum(
                                field_val=work_field.val,
Martin Reinecke's avatar
Martin Reinecke committed
281
                                pdomain=power_domain,
282
                                axes=work_field.domain_axes[space_index])
283
284

        # create the result field and put power_spectrum into it
285
        result_domain = list(work_field.domain)
286
287
        result_domain[space_index] = power_domain

Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
288
        return Field(domain=result_domain, val=power_spectrum,
Martin Reinecke's avatar
updates    
Martin Reinecke committed
289
                     dtype=power_spectrum.dtype)
290

291
    @classmethod
Martin Reinecke's avatar
Martin Reinecke committed
292
    def _calculate_power_spectrum(cls, field_val, pdomain, axes=None):
293

Martin Reinecke's avatar
Martin Reinecke committed
294
        pindex = pdomain.pindex
295
        if axes is not None:
296
297
298
299
            pindex = cls._shape_up_pindex(
                            pindex=pindex,
                            target_shape=field_val.shape,
                            axes=axes)
Theo Steininger's avatar
Theo Steininger committed
300

Martin Reinecke's avatar
Martin Reinecke committed
301
        power_spectrum = utilities.bincount_axis(pindex, weights=field_val,
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
302
                                                 axis=axes)
Martin Reinecke's avatar
Martin Reinecke committed
303
        rho = pdomain.rho
304
305
306
307
308
309
310
311
        if axes is not None:
            new_rho_shape = [1, ] * len(power_spectrum.shape)
            new_rho_shape[axes[0]] = len(rho)
            rho = rho.reshape(new_rho_shape)
        power_spectrum /= rho

        return power_spectrum

312
    @staticmethod
Martin Reinecke's avatar
Martin Reinecke committed
313
    def _shape_up_pindex(pindex, target_shape, axes):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
314
        semiscaled_local_shape = [1] * len(target_shape)
Theo Steininger's avatar
Theo Steininger committed
315
        for i in range(len(axes)):
Martin Reinecke's avatar
Martin Reinecke committed
316
317
            semiscaled_local_shape[axes[i]] = pindex.shape[i]
        result_obj = np.empty(target_shape, dtype=pindex.dtype)
Martin Reinecke's avatar
Martin Reinecke committed
318
        result_obj[()] = pindex.reshape(semiscaled_local_shape)
319
320
        return result_obj

321
    def power_synthesize(self, spaces=None, real_power=True, real_signal=True,
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
322
                         mean=None, std=None):
Theo Steininger's avatar
Theo Steininger committed
323
        """ Yields a sampled field with `self`**2 as its power spectrum.
Theo Steininger's avatar
Theo Steininger committed
324

Theo Steininger's avatar
Theo Steininger committed
325
        This method draws a Gaussian random field in the harmonic partner
Martin Reinecke's avatar
typos    
Martin Reinecke committed
326
        domain of this field's domains, using this field as power spectrum.
Theo Steininger's avatar
Theo Steininger committed
327

328
329
330
        Parameters
        ----------
        spaces : {tuple, int, None} *optional*
Theo Steininger's avatar
Theo Steininger committed
331
332
333
            Specifies the subspace containing all the PowerSpaces which
            should be converted (default : None).
            if spaces==None : Tries to convert the whole domain.
334
        real_power : boolean *optional*
Theo Steininger's avatar
Theo Steininger committed
335
336
            Determines whether the power spectrum is treated as intrinsically
            real or complex (default : True).
337
        real_signal : boolean *optional*
Theo Steininger's avatar
Theo Steininger committed
338
339
340
341
342
343
            True will result in a purely real signal-space field
            (default : True).
        mean : float *optional*
            The mean of the Gaussian noise field which is used for the Field
            synthetization (default : None).
            if mean==None : mean will be set to 0
344
        std : float *optional*
Theo Steininger's avatar
Theo Steininger committed
345
346
347
            The standard deviation of the Gaussian noise field which is used
            for the Field synthetization (default : None).
            if std==None : std will be set to 1
Theo Steininger's avatar
Theo Steininger committed
348

349
350
351
352
        Returns
        -------
        out : Field
            The output object. A random field created with the power spectrum
Theo Steininger's avatar
Theo Steininger committed
353
            stored in the `spaces` in `self`.
354

Theo Steininger's avatar
Theo Steininger committed
355
356
357
358
359
360
        Notes
        -----
        For this the spaces specified by `spaces` must be a PowerSpace.
        This expects this field to be the square root of a power spectrum, i.e.
        to have the unit of the field to be sampled.

361
362
363
        See Also
        --------
        power_analyze
Theo Steininger's avatar
Theo Steininger committed
364
365
366
367
368

        Raises
        ------
        ValueError : If domain specified by `spaces` is not a PowerSpace.

369
        """
Theo Steininger's avatar
Theo Steininger committed
370

371
372
        # check if the `spaces` input is valid
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
Theo Steininger's avatar
Theo Steininger committed
373
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
374
            spaces = list(range(len(self.domain)))
Theo Steininger's avatar
Theo Steininger committed
375

Martin Reinecke's avatar
Martin Reinecke committed
376
377
        for i in spaces:
            if not isinstance(self.domain[i], PowerSpace):
378
379
                raise ValueError("A PowerSpace is needed for field "
                                 "synthetization.")
380
381
382

        # create the result domain
        result_domain = list(self.domain)
Martin Reinecke's avatar
Martin Reinecke committed
383
384
        for i in spaces:
            result_domain[i] = self.domain[i].harmonic_partner
385
386
387

        # create random samples: one or two, depending on whether the
        # power spectrum is real or complex
388
389
        result_list = [self.__class__.from_random(
                             'normal',
390
391
392
                             mean=mean,
                             std=std,
                             domain=result_domain,
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
393
                             dtype=np.complex)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
394
                       for x in range(1 if real_power else 2)]
395
396
397
398
399

        # from now on extract the values from the random fields for further
        # processing without killing the fields.
        # if the signal-space field should be real, hermitianize the field
        # components
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
400
        spec = np.sqrt(self.val)
401
        for power_space_index in spaces:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
402
            spec = self._spec_to_rescaler(spec, power_space_index)
403
404

        # apply the rescaler to the random fields
Martin Reinecke's avatar
Martin Reinecke committed
405
        result_list[0] *= spec.real
406
        if not real_power:
Martin Reinecke's avatar
Martin Reinecke committed
407
            result_list[1] *= spec.imag
408

409
        if real_signal:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
410
            result_list = [Field(i.domain, self._hermitian_decomposition(
Martin Reinecke's avatar
Martin Reinecke committed
411
412
                                     i.val,
                                     preserve_gaussian_variance=True)[0])
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
413
                           for i in result_list]
414
415

        if real_power:
Martin Reinecke's avatar
Martin Reinecke committed
416
            return result_list[0]
417
        else:
Martin Reinecke's avatar
Martin Reinecke committed
418
            return result_list[0] + 1j*result_list[1]
419

420
    @staticmethod
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
421
    def _hermitian_decomposition(val, preserve_gaussian_variance=False):
422
        if preserve_gaussian_variance:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
423
424
425
426
427
            if not issubclass(val.dtype.type, np.complexfloating):
                raise TypeError("complex input field is needed here")
            return (val.real*np.sqrt(2.), val.imag*np.sqrt(2.))
        else:
            return (val.real.copy(), val.imag.copy())
428

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
429
    def _spec_to_rescaler(self, spec, power_space_index):
430
        power_space = self.domain[power_space_index]
431

432
433
434
435
        local_blow_up = [slice(None)]*len(spec.shape)
        # it is important to count from behind, since spec potentially grows
        # with every iteration
        index = self.domain_axes[power_space_index][0]-len(self.shape)
436
        local_blow_up[index] = power_space.pindex
437
        # here, the power_spectrum is distributed into the new shape
438
        return spec[local_blow_up]
439

theos's avatar
theos committed
440
    # ---Properties---
441

theos's avatar
theos committed
442
443
    @property
    def val(self):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
444
        """ Returns the data object associated with this Field.
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
445
        No copy is made.
Theo Steininger's avatar
Theo Steininger committed
446

447
448
        Returns
        -------
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
449
        out : numpy.ndarray
450
        """
Martin Reinecke's avatar
Martin Reinecke committed
451
        return self._val
csongor's avatar
csongor committed
452

Martin Reinecke's avatar
Martin Reinecke committed
453
454
455
456
    @property
    def dtype(self):
        return self._val.dtype

457
458
    @property
    def shape(self):
Theo Steininger's avatar
Theo Steininger committed
459
        """ Returns the total shape of the Field's data array.
Theo Steininger's avatar
Theo Steininger committed
460

461
462
463
        Returns
        -------
        out : tuple
Martin Reinecke's avatar
Martin Reinecke committed
464
            The output object. The tuple contains the dimensions of the spaces
465
            in domain.
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
466
       """
Martin Reinecke's avatar
Martin Reinecke committed
467
        return self._val.shape
csongor's avatar
csongor committed
468

469
470
    @property
    def dim(self):
Theo Steininger's avatar
Theo Steininger committed
471
        """ Returns the total number of pixel-dimensions the field has.
Theo Steininger's avatar
Theo Steininger committed
472

Theo Steininger's avatar
Theo Steininger committed
473
        Effectively, all values from shape are multiplied.
Theo Steininger's avatar
Theo Steininger committed
474

475
476
477
478
479
        Returns
        -------
        out : int
            The dimension of the Field.
        """
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
480
        return self._val.size
csongor's avatar
csongor committed
481

theos's avatar
theos committed
482
483
    @property
    def total_volume(self):
Theo Steininger's avatar
Theo Steininger committed
484
485
        """ Returns the total volume of all spaces in the domain.
        """
Martin Reinecke's avatar
Martin Reinecke committed
486
        if len(self.domain) == 0:
Theo Steininger's avatar
Theo Steininger committed
487
            return 0.
Martin Reinecke's avatar
Martin Reinecke committed
488
489
        volume_tuple = tuple(sp.total_volume for sp in self.domain)
        return reduce(lambda x, y: x * y, volume_tuple)
490

Theo Steininger's avatar
Theo Steininger committed
491
492
493
494
    @property
    def real(self):
        """ The real part of the field (data is not copied).
        """
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
495
        return Field(self.domain, self.val.real)
Theo Steininger's avatar
Theo Steininger committed
496
497
498
499
500

    @property
    def imag(self):
        """ The imaginary part of the field (data is not copied).
        """
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
501
        return Field(self.domain, self.val.imag)
Theo Steininger's avatar
Theo Steininger committed
502

theos's avatar
theos committed
503
    # ---Special unary/binary operations---
504

Martin Reinecke's avatar
Martin Reinecke committed
505
    def copy(self):
506
        """ Returns a full copy of the Field.
Theo Steininger's avatar
Theo Steininger committed
507

Martin Reinecke's avatar
Martin Reinecke committed
508
        The returned object will be an identical copy of the original Field.
Theo Steininger's avatar
Theo Steininger committed
509

510
511
512
513
514
        Returns
        -------
        out : Field
            The output object. An identical copy of 'self'.
        """
Martin Reinecke's avatar
Martin Reinecke committed
515
        return Field(val=self, copy=True)
csongor's avatar
csongor committed
516

517
518
519
520
521
522
    def scalar_weight(self, spaces=None):
        if np.isscalar(spaces):
            return self.domain[spaces].scalar_weight()

        if spaces is None:
            spaces = range(len(self.domain))
Martin Reinecke's avatar
Martin Reinecke committed
523
        res = 1.
524
525
526
527
528
529
530
        for i in spaces:
            tmp = self.domain[i].scalar_weight()
            if tmp is None:
                return None
            res *= tmp
        return res

theos's avatar
theos committed
531
    def weight(self, power=1, inplace=False, spaces=None):
Theo Steininger's avatar
Theo Steininger committed
532
        """ Weights the pixels of `self` with their invidual pixel-volume.
533
534
535
536

        Parameters
        ----------
        power : number
Theo Steininger's avatar
Theo Steininger committed
537
            The pixels get weighted with the volume-factor**power.
Theo Steininger's avatar
Theo Steininger committed
538

539
        inplace : boolean
Theo Steininger's avatar
Theo Steininger committed
540
541
            If True, `self` will be weighted and returned. Otherwise, a copy
            is made.
Theo Steininger's avatar
Theo Steininger committed
542

Theo Steininger's avatar
Theo Steininger committed
543
544
        spaces : tuple of ints
            Determines on which subspace the operation takes place.
Theo Steininger's avatar
Theo Steininger committed
545

546
547
548
        Returns
        -------
        out : Field
Theo Steininger's avatar
Theo Steininger committed
549
            The weighted field.
550
551

        """
552
        new_field = Field(val=self, copy=not inplace)
csongor's avatar
csongor committed
553

554
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
csongor's avatar
csongor committed
555
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
556
            spaces = list(range(len(self.domain)))
csongor's avatar
csongor committed
557

558
559
560
561
562
563
564
565
566
        fct = 1.
        for ind in spaces:
            wgt = self.domain[ind].weight()
            if np.isscalar(wgt):
                fct *= wgt
            else:
                new_shape = np.ones(len(self.shape), dtype=np.int)
                new_shape[self.domain_axes[ind][0]:self.domain_axes[ind][-1]+1]=wgt.shape
                wgt = wgt.reshape(new_shape)
567
                new_field *= wgt**power
568
        fct = fct**power
Martin Reinecke's avatar
Martin Reinecke committed
569
        if fct != 1.:
570
            new_field *= fct
571

572
        return new_field
csongor's avatar
csongor committed
573

Martin Reinecke's avatar
Martin Reinecke committed
574
    def vdot(self, x=None, spaces=None, bare=False):
Theo Steininger's avatar
Theo Steininger committed
575
        """ Computes the volume-factor-aware dot product of 'self' with x.
Theo Steininger's avatar
Theo Steininger committed
576

577
578
579
        Parameters
        ----------
        x : Field
Theo Steininger's avatar
Theo Steininger committed
580
            The domain of x must contain `self.domain`
Theo Steininger's avatar
Theo Steininger committed
581

Theo Steininger's avatar
Theo Steininger committed
582
583
584
        spaces : tuple of ints
            If the domain of `self` and `x` are not the same, `spaces` specfies
            the mapping.
Theo Steininger's avatar
Theo Steininger committed
585

586
        bare : boolean
Theo Steininger's avatar
Theo Steininger committed
587
            If true, no volume factors will be included in the computation.
Theo Steininger's avatar
Theo Steininger committed
588

589
590
591
        Returns
        -------
        out : float, complex
Theo Steininger's avatar
Theo Steininger committed
592

593
        """
594
595
596
        if not isinstance(x, Field):
            raise ValueError("The dot-partner must be an instance of " +
                             "the NIFTy field class")
theos's avatar
theos committed
597

Martin Reinecke's avatar
Martin Reinecke committed
598
        # Compute the dot respecting the fact of discrete/continuous spaces
599
600
601
602
603
604
605
606
607
608
        fct = 1.
        if bare:
            y = self
        else:
            tmp = self.scalar_weight(spaces)
            if tmp is None:
                y = self.weight(power=1)
            else:
                y = self
                fct = tmp
theos's avatar
theos committed
609

610
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
611
            return fct*np.vdot(y.val.ravel(), x.val.ravel())
612
613
614
        else:
            # create a diagonal operator which is capable of taking care of the
            # axes-matching
615
            from .operators.diagonal_operator import DiagonalOperator
616
617
618
619
620
            diagonal = y.val.conjugate()
            diagonalOperator = DiagonalOperator(domain=y.domain,
                                                diagonal=diagonal,
                                                copy=False)
            dotted = diagonalOperator(x, spaces=spaces)
621
            return fct*dotted.sum(spaces=spaces)
theos's avatar
theos committed
622

Theo Steininger's avatar
Theo Steininger committed
623
    def norm(self):
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
624
        """ Computes the L2-norm of the field values.
csongor's avatar
csongor committed
625

Theo Steininger's avatar
Theo Steininger committed
626
627
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
628
        norm : float
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
629
            The L2-norm of the field values.
csongor's avatar
csongor committed
630
631

        """
Martin Reinecke's avatar
fixes    
Martin Reinecke committed
632
        return np.sqrt(np.abs(self.vdot(x=self)))
csongor's avatar
csongor committed
633
634

    def conjugate(self, inplace=False):
635
        """ Retruns the complex conjugate of the field.
Theo Steininger's avatar
Theo Steininger committed
636

637
638
639
        Parameters
        ----------
        inplace : boolean
Theo Steininger's avatar
Theo Steininger committed
640
            Decides whether the conjugation should be performed inplace.
Theo Steininger's avatar
Theo Steininger committed
641

642
643
644
645
        Returns
        -------
        cc : field
            The complex conjugated field.
csongor's avatar
csongor committed
646
647
648

        """
        if inplace:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
649
            self.imag *= -1
650
            return self
csongor's avatar
csongor committed
651
        else:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
652
            return Field(self.domain, np.conj(self.val), self.dtype)
csongor's avatar
csongor committed
653

theos's avatar
theos committed
654
    # ---General unary/contraction methods---
655

theos's avatar
theos committed
656
657
    def __pos__(self):
        return self.copy()
658

theos's avatar
theos committed
659
    def __neg__(self):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
660
        return Field(self.domain, -self.val, self.dtype)
csongor's avatar
csongor committed
661

theos's avatar
theos committed
662
    def __abs__(self):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
663
        return Field(self.domain, np.abs(self.val), self.dtype)
csongor's avatar
csongor committed
664

665
    def _contraction_helper(self, op, spaces):
theos's avatar
theos committed
666
        if spaces is None:
667
668
669
            return getattr(self.val, op)()
        # build a list of all axes
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
csongor's avatar
csongor committed
670

671
        axes_list = tuple(self.domain_axes[sp_index] for sp_index in spaces)
672

Martin Reinecke's avatar
Martin Reinecke committed
673
        if len(axes_list) > 0:
theos's avatar
theos committed
674
            axes_list = reduce(lambda x, y: x+y, axes_list)
csongor's avatar
csongor committed
675

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
676
        # perform the contraction on the data
677
        data = getattr(self.val, op)(axis=axes_list)
csongor's avatar
csongor committed
678

theos's avatar
theos committed
679
680
681
        # check if the result is scalar or if a result_field must be constr.
        if np.isscalar(data):
            return data
csongor's avatar
csongor committed
682
        else:
theos's avatar
theos committed
683
            return_domain = tuple(self.domain[i]
Martin Reinecke's avatar
Martin Reinecke committed
684
                                  for i in range(len(self.domain))
theos's avatar
theos committed
685
                                  if i not in spaces)
686

Martin Reinecke's avatar
updates    
Martin Reinecke committed
687
            return Field(domain=return_domain, val=data, copy=False)
csongor's avatar
csongor committed
688

689
690
    def sum(self, spaces=None):
        return self._contraction_helper('sum', spaces)
csongor's avatar
csongor committed
691

692
693
    def prod(self, spaces=None):
        return self._contraction_helper('prod', spaces)
csongor's avatar
csongor committed
694

695
696
    def all(self, spaces=None):
        return self._contraction_helper('all', spaces)
csongor's avatar
csongor committed
697

698
699
    def any(self, spaces=None):
        return self._contraction_helper('any', spaces)
csongor's avatar
csongor committed
700

701
702
    def min(self, spaces=None):
        return self._contraction_helper('min', spaces)
csongor's avatar
csongor committed
703

704
705
    def max(self, spaces=None):
        return self._contraction_helper('max', spaces)
csongor's avatar
csongor committed
706

707
708
    def mean(self, spaces=None):
        return self._contraction_helper('mean', spaces)
csongor's avatar
csongor committed
709

710
711
    def var(self, spaces=None):
        return self._contraction_helper('var', spaces)
csongor's avatar
csongor committed
712

713
714
    def std(self, spaces=None):
        return self._contraction_helper('std', spaces)
csongor's avatar
csongor committed
715

theos's avatar
theos committed
716
    # ---General binary methods---
csongor's avatar
csongor committed
717

718
    def _binary_helper(self, other, op):
csongor's avatar
csongor committed
719
        # if other is a field, make sure that the domains match
720
        if isinstance(other, Field):
721
722
            if other.domain != self.domain:
                raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
723
            return Field(self.domain, getattr(self.val, op)(other.val))
csongor's avatar
csongor committed
724

Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
725
        return Field(self.domain, getattr(self.val, op)(other))
csongor's avatar
csongor committed
726
727

    def __add__(self, other):
theos's avatar
theos committed
728
        return self._binary_helper(other, op='__add__')
729

730
    def __radd__(self, other):
theos's avatar
theos committed
731
        return self._binary_helper(other, op='__radd__')
csongor's avatar
csongor committed
732
733

    def __iadd__(self, other):
734
        return self._binary_helper(other, op='__iadd__')
csongor's avatar
csongor committed
735
736

    def __sub__(self, other):
theos's avatar
theos committed
737
        return self._binary_helper(other, op='__sub__')
csongor's avatar
csongor committed
738
739

    def __rsub__(self, other):
theos's avatar
theos committed
740
        return self._binary_helper(other, op='__rsub__')
csongor's avatar
csongor committed
741
742

    def __isub__(self, other):
743
        return self._binary_helper(other, op='__isub__')
csongor's avatar
csongor committed
744
745

    def __mul__(self, other):
theos's avatar
theos committed
746
        return self._binary_helper(other, op='__mul__')
747

748
    def __rmul__(self, other):
theos's avatar
theos committed
749
        return self._binary_helper(other, op='__rmul__')
csongor's avatar
csongor committed
750
751

    def __imul__(self, other):
752
        return self._binary_helper(other, op='__imul__')
csongor's avatar
csongor committed
753
754

    def __div__(self, other):
theos's avatar
theos committed
755
        return self._binary_helper(other, op='__div__')
csongor's avatar
csongor committed
756

Martin Reinecke's avatar
Martin Reinecke committed
757
758
759
    def __truediv__(self, other):
        return self._binary_helper(other, op='__truediv__')

csongor's avatar
csongor committed
760
    def __rdiv__(self, other):
theos's avatar
theos committed
761
        return self._binary_helper(other, op='__rdiv__')
csongor's avatar
csongor committed
762

Martin Reinecke's avatar
Martin Reinecke committed
763
764
765
    def __rtruediv__(self, other):
        return self._binary_helper(other, op='__rtruediv__')

csongor's avatar
csongor committed
766
    def __idiv__(self, other):
767
        return self._binary_helper(other, op='__idiv__')
768

csongor's avatar
csongor committed
769
    def __pow__(self, other):
theos's avatar
theos committed
770
        return self._binary_helper(other, op='__pow__')
csongor's avatar
csongor committed
771
772

    def __rpow__(self, other):
theos's avatar
theos committed
773
        return self._binary_helper(other, op='__rpow__')
csongor's avatar
csongor committed
774
775

    def __ipow__(self, other):
776
        return self._binary_helper(other, op='__ipow__')
csongor's avatar
csongor committed
777
778

    def __lt__(self, other):
theos's avatar
theos committed
779
        return self._binary_helper(other, op='__lt__')
csongor's avatar
csongor committed
780
781

    def __le__(self, other):
theos's avatar
theos committed
782
        return self._binary_helper(other, op='__le__')
csongor's avatar
csongor committed
783
784
785
786
787

    def __ne__(self, other):
        if other is None:
            return True
        else:
theos's avatar
theos committed
788
            return self._binary_helper(other, op='__ne__')
csongor's avatar
csongor committed
789
790
791
792
793

    def __eq__(self, other):
        if other is None:
            return False
        else:
theos's avatar
theos committed
794
            return self._binary_helper(other, op='__eq__')
csongor's avatar
csongor committed
795
796

    def __ge__(self, other):
theos's avatar
theos committed
797
        return self._binary_helper(other, op='__ge__')
csongor's avatar
csongor committed
798
799

    def __gt__(self, other):
theos's avatar
theos committed
800
801
802
803
804
805
806
807
808
809
        return self._binary_helper(other, op='__gt__')

    def __repr__(self):
        return "<nifty_core.field>"

    def __str__(self):
        minmax = [self.min(), self.max()]
        mean = self.mean()
        return "nifty_core.field instance\n- domain      = " + \
               repr(self.domain) + \
810
               "\n- val         = " + repr(self.val) + \
theos's avatar
theos committed
811
812
               "\n  - min.,max. = " + str(minmax) + \
               "\n  - mean = " + str(mean)