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

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

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

        return utilities.parse_domain(domain)

    @staticmethod
    def _get_axes_tuple(things_with_shape):
119
        i = 0
theos's avatar
theos committed
120
121
        axes_list = []
        for thing in things_with_shape:
122
            nax = len(thing.shape)
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
123
            axes_list += [tuple(range(i, i+nax))]
124
            i += nax
theos's avatar
theos committed
125
        return tuple(axes_list)
126

Martin Reinecke's avatar
Martin Reinecke committed
127
128
    @staticmethod
    def _infer_dtype(dtype, val):
Martin Reinecke's avatar
Martin Reinecke committed
129
130
131
132
133
        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)
134

135
    # ---Factory methods---
136

Martin Reinecke's avatar
Martin Reinecke committed
137
138
    @staticmethod
    def from_random(random_type, domain, dtype=np.float64, **kwargs):
139
140
141
142
143
144
145
        """ Draws a random field with the given parameters.

        Parameters
        ----------
        random_type : String
            'pm1', 'normal', 'uniform' are the supported arguments for this
            method.
Theo Steininger's avatar
Theo Steininger committed
146

147
148
        domain : DomainObject
            The domain of the output random field
Theo Steininger's avatar
Theo Steininger committed
149

150
151
        dtype : type
            The datatype of the output random field
Theo Steininger's avatar
Theo Steininger committed
152

153
154
155
156
157
158
159
        Returns
        -------
        out : Field
            The output object.

        See Also
        --------
160
        power_synthesize
161
        """
Theo Steininger's avatar
Theo Steininger committed
162

163
        generator_function = getattr(Random, random_type)
Martin Reinecke's avatar
Martin Reinecke committed
164
        return Field(domain=domain,
Martin Reinecke's avatar
Martin Reinecke committed
165
                     val=generator_function(dtype=dtype,
Martin Reinecke's avatar
Martin Reinecke committed
166
                     shape=utilities.domains2shape(domain), **kwargs))
167

168
169
    # ---Powerspectral methods---

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

Theo Steininger's avatar
Theo Steininger committed
174
175
176
        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
177
        harmonic space. The resulting field has the same units as the initial
Theo Steininger's avatar
Theo Steininger committed
178
        field, corresponding to the square root of the power spectrum.
179
180
181

        Parameters
        ----------
Theo Steininger's avatar
Theo Steininger committed
182
183
184
185
186
        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
187
            if binbounds==None : bins are inferred.
188
189
190
191
192
193
194
195
196
197
        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
198

199
200
201
202
        Raise
        -----
        ValueError
            Raised if
Theo Steininger's avatar
Theo Steininger committed
203
204
                *len(domain) is != 1 when spaces==None
                *len(spaces) is != 1 if not None
205
                *the analyzed space is not harmonic
Theo Steininger's avatar
Theo Steininger committed
206

207
208
        Returns
        -------
Theo Steininger's avatar
Theo Steininger committed
209
        out : Field
Martin Reinecke's avatar
Martin Reinecke committed
210
            The output object. Its domain is a PowerSpace and it contains
211
212
213
214
215
            the power spectrum of 'self's field.

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

217
        """
Theo Steininger's avatar
Theo Steininger committed
218

Theo Steininger's avatar
Theo Steininger committed
219
        # check if all spaces in `self.domain` are either harmonic or
220
221
222
        # power_space instances
        for sp in self.domain:
            if not sp.harmonic and not isinstance(sp, PowerSpace):
Martin Reinecke's avatar
Martin Reinecke committed
223
224
                raise TypeError("Field has a space in `domain` which is "
                                "neither harmonic nor a PowerSpace.")
225
226

        # check if the `spaces` input is valid
227
228
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
229
            spaces = range(len(self.domain))
230
231

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

234
        if keep_phase_information:
Martin Reinecke's avatar
Martin Reinecke committed
235
            parts = self._hermitian_decomposition(self, False)
236
237
238
239
        else:
            parts = [self]

        parts = [abs(part)**2 for part in parts]
240
241

        for space_index in spaces:
242
243
            parts = [self._single_power_analyze(
                                work_field=part,
244
                                space_index=space_index,
245
246
                                binbounds=binbounds)
                     for part in parts]
247

248
        if keep_phase_information:
Martin Reinecke's avatar
updates    
Martin Reinecke committed
249
            return parts[0] + 1j*parts[1]
250
        else:
Martin Reinecke's avatar
updates    
Martin Reinecke committed
251
            return parts[0]
252

Martin Reinecke's avatar
Martin Reinecke committed
253
254
    @staticmethod
    def _single_power_analyze(work_field, space_index, binbounds):
255
        if not work_field.domain[space_index].harmonic:
Martin Reinecke's avatar
Martin Reinecke committed
256
            raise ValueError("The analyzed space must be harmonic.")
257

258
259
260
261
262
263
        # 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.

264
        harmonic_domain = work_field.domain[space_index]
265
        power_domain = PowerSpace(harmonic_partner=harmonic_domain,
Theo Steininger's avatar
Theo Steininger committed
266
                                  binbounds=binbounds)
Martin Reinecke's avatar
Martin Reinecke committed
267
        power_spectrum = Field._calculate_power_spectrum(
268
                                field_val=work_field.val,
Martin Reinecke's avatar
Martin Reinecke committed
269
                                pdomain=power_domain,
270
                                axes=work_field.domain_axes[space_index])
271
272

        # create the result field and put power_spectrum into it
273
        result_domain = list(work_field.domain)
274
275
        result_domain[space_index] = power_domain

Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
276
        return Field(domain=result_domain, val=power_spectrum,
Martin Reinecke's avatar
updates    
Martin Reinecke committed
277
                     dtype=power_spectrum.dtype)
278

Martin Reinecke's avatar
Martin Reinecke committed
279
280
    @staticmethod
    def _calculate_power_spectrum(field_val, pdomain, axes=None):
281

Martin Reinecke's avatar
Martin Reinecke committed
282
        pindex = pdomain.pindex
283
        if axes is not None:
Martin Reinecke's avatar
Martin Reinecke committed
284
            pindex = Field._shape_up_pindex(pindex, field_val.shape, axes)
Theo Steininger's avatar
Theo Steininger committed
285

Martin Reinecke's avatar
Martin Reinecke committed
286
        power_spectrum = utilities.bincount_axis(pindex, weights=field_val,
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
287
                                                 axis=axes)
Martin Reinecke's avatar
Martin Reinecke committed
288
        rho = pdomain.rho
289
290
291
292
293
294
295
296
        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

297
    @staticmethod
Martin Reinecke's avatar
Martin Reinecke committed
298
    def _shape_up_pindex(pindex, target_shape, axes):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
299
        semiscaled_local_shape = [1] * len(target_shape)
Martin Reinecke's avatar
Martin Reinecke committed
300
301
        for i, ax in enumerate(axes):
            semiscaled_local_shape[ax] = pindex.shape[i]
Martin Reinecke's avatar
Martin Reinecke committed
302
        result_obj = np.empty(target_shape, dtype=pindex.dtype)
Martin Reinecke's avatar
Martin Reinecke committed
303
        result_obj[()] = pindex.reshape(semiscaled_local_shape)
304
305
        return result_obj

306
    def power_synthesize(self, spaces=None, real_power=True, real_signal=True,
Martin Reinecke's avatar
Martin Reinecke committed
307
                         mean=0., std=1.):
Theo Steininger's avatar
Theo Steininger committed
308
        """ Yields a sampled field with `self`**2 as its power spectrum.
Theo Steininger's avatar
Theo Steininger committed
309

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

313
314
315
        Parameters
        ----------
        spaces : {tuple, int, None} *optional*
Theo Steininger's avatar
Theo Steininger committed
316
317
318
            Specifies the subspace containing all the PowerSpaces which
            should be converted (default : None).
            if spaces==None : Tries to convert the whole domain.
319
        real_power : boolean *optional*
Theo Steininger's avatar
Theo Steininger committed
320
321
            Determines whether the power spectrum is treated as intrinsically
            real or complex (default : True).
322
        real_signal : boolean *optional*
Theo Steininger's avatar
Theo Steininger committed
323
324
325
326
327
328
            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
329
        std : float *optional*
Theo Steininger's avatar
Theo Steininger committed
330
331
332
            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
333

334
335
336
337
        Returns
        -------
        out : Field
            The output object. A random field created with the power spectrum
Theo Steininger's avatar
Theo Steininger committed
338
            stored in the `spaces` in `self`.
339

Theo Steininger's avatar
Theo Steininger committed
340
341
342
343
344
345
        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.

346
347
348
        See Also
        --------
        power_analyze
Theo Steininger's avatar
Theo Steininger committed
349
350
351
352
353

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

354
        """
Theo Steininger's avatar
Theo Steininger committed
355

356
357
        # check if the `spaces` input is valid
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
Theo Steininger's avatar
Theo Steininger committed
358
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
359
            spaces = range(len(self.domain))
Theo Steininger's avatar
Theo Steininger committed
360

Martin Reinecke's avatar
Martin Reinecke committed
361
362
        for i in spaces:
            if not isinstance(self.domain[i], PowerSpace):
363
364
                raise ValueError("A PowerSpace is needed for field "
                                 "synthetization.")
365
366
367

        # create the result domain
        result_domain = list(self.domain)
Martin Reinecke's avatar
Martin Reinecke committed
368
369
        for i in spaces:
            result_domain[i] = self.domain[i].harmonic_partner
370
371
372

        # create random samples: one or two, depending on whether the
        # power spectrum is real or complex
373
374
        result_list = [self.__class__.from_random(
                             'normal',
375
376
377
                             mean=mean,
                             std=std,
                             domain=result_domain,
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
378
                             dtype=np.complex)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
379
                       for x in range(1 if real_power else 2)]
380
381
382
383
384

        # 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
385
        spec = np.sqrt(self.val)
386
        for power_space_index in spaces:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
387
            spec = self._spec_to_rescaler(spec, power_space_index)
388
389

        # apply the rescaler to the random fields
Martin Reinecke's avatar
Martin Reinecke committed
390
        result_list[0] *= spec.real
391
        if not real_power:
Martin Reinecke's avatar
Martin Reinecke committed
392
            result_list[1] *= spec.imag
393

394
        if real_signal:
Martin Reinecke's avatar
Martin Reinecke committed
395
            result_list = [self._hermitian_decomposition(i, True)[0]
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
396
                           for i in result_list]
397
398

        if real_power:
Martin Reinecke's avatar
Martin Reinecke committed
399
            return result_list[0]
400
        else:
Martin Reinecke's avatar
Martin Reinecke committed
401
            return result_list[0] + 1j*result_list[1]
402

403
    @staticmethod
Martin Reinecke's avatar
Martin Reinecke committed
404
    def _hermitian_decomposition(inp, preserve_gaussian_variance=False):
405
        if preserve_gaussian_variance:
Martin Reinecke's avatar
Martin Reinecke committed
406
            if not issubclass(inp.dtype.type, np.complexfloating):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
407
                raise TypeError("complex input field is needed here")
Martin Reinecke's avatar
Martin Reinecke committed
408
            return (inp.real*np.sqrt(2.), inp.imag*np.sqrt(2.))
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
409
        else:
Martin Reinecke's avatar
Martin Reinecke committed
410
            return (inp.real.copy(), inp.imag.copy())
411

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
412
    def _spec_to_rescaler(self, spec, power_space_index):
413
        power_space = self.domain[power_space_index]
414

415
416
417
418
        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)
419
        local_blow_up[index] = power_space.pindex
420
        # here, the power_spectrum is distributed into the new shape
421
        return spec[local_blow_up]
422

theos's avatar
theos committed
423
    # ---Properties---
424

theos's avatar
theos committed
425
426
    @property
    def val(self):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
427
        """ Returns the data object associated with this Field.
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
428
        No copy is made.
Theo Steininger's avatar
Theo Steininger committed
429

430
431
        Returns
        -------
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
432
        out : numpy.ndarray
433
        """
Martin Reinecke's avatar
Martin Reinecke committed
434
        return self._val
csongor's avatar
csongor committed
435

Martin Reinecke's avatar
Martin Reinecke committed
436
437
438
439
    @property
    def dtype(self):
        return self._val.dtype

440
441
    @property
    def shape(self):
Theo Steininger's avatar
Theo Steininger committed
442
        """ Returns the total shape of the Field's data array.
Theo Steininger's avatar
Theo Steininger committed
443

444
445
446
        Returns
        -------
        out : tuple
Martin Reinecke's avatar
Martin Reinecke committed
447
            The output object. The tuple contains the dimensions of the spaces
448
            in domain.
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
449
       """
Martin Reinecke's avatar
Martin Reinecke committed
450
        return self._val.shape
csongor's avatar
csongor committed
451

452
453
    @property
    def dim(self):
Theo Steininger's avatar
Theo Steininger committed
454
        """ Returns the total number of pixel-dimensions the field has.
Theo Steininger's avatar
Theo Steininger committed
455

Theo Steininger's avatar
Theo Steininger committed
456
        Effectively, all values from shape are multiplied.
Theo Steininger's avatar
Theo Steininger committed
457

458
459
460
461
462
        Returns
        -------
        out : int
            The dimension of the Field.
        """
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
463
        return self._val.size
csongor's avatar
csongor committed
464

theos's avatar
theos committed
465
466
    @property
    def total_volume(self):
Theo Steininger's avatar
Theo Steininger committed
467
468
        """ Returns the total volume of all spaces in the domain.
        """
Martin Reinecke's avatar
Martin Reinecke committed
469
        if len(self.domain) == 0:
Theo Steininger's avatar
Theo Steininger committed
470
            return 0.
Martin Reinecke's avatar
Martin Reinecke committed
471
472
        volume_tuple = tuple(sp.total_volume for sp in self.domain)
        return reduce(lambda x, y: x * y, volume_tuple)
473

Theo Steininger's avatar
Theo Steininger committed
474
475
476
477
    @property
    def real(self):
        """ The real part of the field (data is not copied).
        """
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
478
        return Field(self.domain, self.val.real)
Theo Steininger's avatar
Theo Steininger committed
479
480
481
482
483

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

theos's avatar
theos committed
486
    # ---Special unary/binary operations---
487

Martin Reinecke's avatar
Martin Reinecke committed
488
    def copy(self):
489
        """ Returns a full copy of the Field.
Theo Steininger's avatar
Theo Steininger committed
490

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

493
494
495
496
497
        Returns
        -------
        out : Field
            The output object. An identical copy of 'self'.
        """
Martin Reinecke's avatar
Martin Reinecke committed
498
        return Field(val=self, copy=True)
csongor's avatar
csongor committed
499

500
501
502
503
504
505
    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
506
        res = 1.
507
508
509
510
511
512
513
        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
514
    def weight(self, power=1, inplace=False, spaces=None):
Theo Steininger's avatar
Theo Steininger committed
515
        """ Weights the pixels of `self` with their invidual pixel-volume.
516
517
518
519

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

522
        inplace : boolean
Theo Steininger's avatar
Theo Steininger committed
523
524
            If True, `self` will be weighted and returned. Otherwise, a copy
            is made.
Theo Steininger's avatar
Theo Steininger committed
525

Theo Steininger's avatar
Theo Steininger committed
526
527
        spaces : tuple of ints
            Determines on which subspace the operation takes place.
Theo Steininger's avatar
Theo Steininger committed
528

529
530
531
        Returns
        -------
        out : Field
Theo Steininger's avatar
Theo Steininger committed
532
            The weighted field.
533
534

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

537
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
csongor's avatar
csongor committed
538
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
539
            spaces = range(len(self.domain))
csongor's avatar
csongor committed
540

541
542
543
544
545
546
547
548
549
        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)
550
                new_field *= wgt**power
551
        fct = fct**power
Martin Reinecke's avatar
Martin Reinecke committed
552
        if fct != 1.:
553
            new_field *= fct
554

555
        return new_field
csongor's avatar
csongor committed
556

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

560
561
562
        Parameters
        ----------
        x : Field
Theo Steininger's avatar
Theo Steininger committed
563
            The domain of x must contain `self.domain`
Theo Steininger's avatar
Theo Steininger committed
564

Theo Steininger's avatar
Theo Steininger committed
565
566
567
        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
568

569
        bare : boolean
Theo Steininger's avatar
Theo Steininger committed
570
            If true, no volume factors will be included in the computation.
Theo Steininger's avatar
Theo Steininger committed
571

572
573
574
        Returns
        -------
        out : float, complex
Theo Steininger's avatar
Theo Steininger committed
575

576
        """
577
578
579
        if not isinstance(x, Field):
            raise ValueError("The dot-partner must be an instance of " +
                             "the NIFTy field class")
theos's avatar
theos committed
580

Martin Reinecke's avatar
Martin Reinecke committed
581
        # Compute the dot respecting the fact of discrete/continuous spaces
582
583
584
585
586
587
588
589
590
591
        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
592

593
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
594
            return fct*np.vdot(y.val.ravel(), x.val.ravel())
595
596
597
        else:
            # create a diagonal operator which is capable of taking care of the
            # axes-matching
598
            from .operators.diagonal_operator import DiagonalOperator
599
600
601
602
603
            diagonal = y.val.conjugate()
            diagonalOperator = DiagonalOperator(domain=y.domain,
                                                diagonal=diagonal,
                                                copy=False)
            dotted = diagonalOperator(x, spaces=spaces)
604
            return fct*dotted.sum(spaces=spaces)
theos's avatar
theos committed
605

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

Theo Steininger's avatar
Theo Steininger committed
609
610
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
611
        norm : float
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
612
            The L2-norm of the field values.
csongor's avatar
csongor committed
613
614

        """
Martin Reinecke's avatar
fixes    
Martin Reinecke committed
615
        return np.sqrt(np.abs(self.vdot(x=self)))
csongor's avatar
csongor committed
616
617

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

620
621
622
        Parameters
        ----------
        inplace : boolean
Theo Steininger's avatar
Theo Steininger committed
623
            Decides whether the conjugation should be performed inplace.
Theo Steininger's avatar
Theo Steininger committed
624

625
626
627
628
        Returns
        -------
        cc : field
            The complex conjugated field.
csongor's avatar
csongor committed
629
630
631

        """
        if inplace:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
632
            self.imag *= -1
633
            return self
csongor's avatar
csongor committed
634
        else:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
635
            return Field(self.domain, np.conj(self.val), self.dtype)
csongor's avatar
csongor committed
636

theos's avatar
theos committed
637
    # ---General unary/contraction methods---
638

theos's avatar
theos committed
639
640
    def __pos__(self):
        return self.copy()
641

theos's avatar
theos committed
642
    def __neg__(self):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
643
        return Field(self.domain, -self.val, self.dtype)
csongor's avatar
csongor committed
644

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

648
    def _contraction_helper(self, op, spaces):
theos's avatar
theos committed
649
        if spaces is None:
650
651
652
            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
653

654
        axes_list = tuple(self.domain_axes[sp_index] for sp_index in spaces)
655

Martin Reinecke's avatar
Martin Reinecke committed
656
        if len(axes_list) > 0:
theos's avatar
theos committed
657
            axes_list = reduce(lambda x, y: x+y, axes_list)
csongor's avatar
csongor committed
658

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

theos's avatar
theos committed
662
663
664
        # 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
665
        else:
theos's avatar
theos committed
666
            return_domain = tuple(self.domain[i]
Martin Reinecke's avatar
Martin Reinecke committed
667
                                  for i in range(len(self.domain))
theos's avatar
theos committed
668
                                  if i not in spaces)
669

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

672
673
    def sum(self, spaces=None):
        return self._contraction_helper('sum', spaces)
csongor's avatar
csongor committed
674

675
676
    def prod(self, spaces=None):
        return self._contraction_helper('prod', spaces)
csongor's avatar
csongor committed
677

678
679
    def all(self, spaces=None):
        return self._contraction_helper('all', spaces)
csongor's avatar
csongor committed
680

681
682
    def any(self, spaces=None):
        return self._contraction_helper('any', spaces)
csongor's avatar
csongor committed
683

684
685
    def min(self, spaces=None):
        return self._contraction_helper('min', spaces)
csongor's avatar
csongor committed
686

687
688
    def max(self, spaces=None):
        return self._contraction_helper('max', spaces)
csongor's avatar
csongor committed
689

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

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

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

theos's avatar
theos committed
699
    # ---General binary methods---
csongor's avatar
csongor committed
700

701
    def _binary_helper(self, other, op):
csongor's avatar
csongor committed
702
        # if other is a field, make sure that the domains match
703
        if isinstance(other, Field):
704
705
            if other.domain != self.domain:
                raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
706
            return Field(self.domain, getattr(self.val, op)(other.val))
csongor's avatar
csongor committed
707

Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
708
        return Field(self.domain, getattr(self.val, op)(other))
csongor's avatar
csongor committed
709
710

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

713
    def __radd__(self, other):
theos's avatar
theos committed
714
        return self._binary_helper(other, op='__radd__')
csongor's avatar
csongor committed
715
716

    def __iadd__(self, other):
717
        return self._binary_helper(other, op='__iadd__')
csongor's avatar
csongor committed
718
719

    def __sub__(self, other):
theos's avatar
theos committed
720
        return self._binary_helper(other, op='__sub__')
csongor's avatar
csongor committed
721
722

    def __rsub__(self, other):
theos's avatar
theos committed
723
        return self._binary_helper(other, op='__rsub__')
csongor's avatar
csongor committed
724
725

    def __isub__(self, other):
726
        return self._binary_helper(other, op='__isub__')
csongor's avatar
csongor committed
727
728

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

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

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

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

Martin Reinecke's avatar
Martin Reinecke committed
740
741
742
    def __truediv__(self, other):
        return self._binary_helper(other, op='__truediv__')

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

Martin Reinecke's avatar
Martin Reinecke committed
746
747
748
    def __rtruediv__(self, other):
        return self._binary_helper(other, op='__rtruediv__')

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

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

    def __rpow__(self, other):
theos's avatar
theos committed
756
        return self._binary_helper(other, op='__rpow__')
csongor's avatar
csongor committed
757
758

    def __ipow__(self, other):
759
        return self._binary_helper(other, op='__ipow__')
csongor's avatar
csongor committed
760
761

    def __lt__(self, other):
theos's avatar
theos committed
762
        return self._binary_helper(other, op='__lt__')
csongor's avatar
csongor committed
763
764

    def __le__(self, other):
theos's avatar
theos committed
765
        return self._binary_helper(other, op='__le__')
csongor's avatar
csongor committed
766
767
768
769
770

    def __ne__(self, other):
        if other is None:
            return True
        else:
theos's avatar
theos committed
771
            return self._binary_helper(other, op='__ne__')
csongor's avatar
csongor committed
772
773
774
775
776

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

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

    def __gt__(self, other):
theos's avatar
theos committed
783
784
785
        return self._binary_helper(other, op='__gt__')

    def __repr__(self):
Martin Reinecke's avatar
Martin Reinecke committed
786
        return "<nifty2go.Field>"
theos's avatar
theos committed
787
788
789
790

    def __str__(self):
        minmax = [self.min(), self.max()]
        mean = self.mean()
Martin Reinecke's avatar
Martin Reinecke committed
791
        return "nifty2go.Field instance\n- domain      = " + \
theos's avatar
theos committed
792
               repr(self.domain) + \
793
               "\n- val         = " + repr(self.val) + \
theos's avatar
theos committed
794
795
               "\n  - min.,max. = " + str(minmax) + \
               "\n  - mean = " + str(mean)