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

Theo Steininger's avatar
Theo Steininger 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
Theo Steininger's avatar
Theo Steininger 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
Theo Steininger's avatar
Theo Steininger 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
Theo Steininger's avatar
Theo Steininger committed
169

170
171

        """
Theo Steininger's avatar
Theo Steininger committed
172

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

178
179
    # ---Powerspectral methods---

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

Theo Steininger's avatar
Theo Steininger committed
184
185
186
        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
187
        harmonic space. The resulting field has the same units as the initial
Theo Steininger's avatar
Theo Steininger committed
188
        field, corresponding to the square root of the power spectrum.
189
190
191

        Parameters
        ----------
Theo Steininger's avatar
Theo Steininger committed
192
193
194
195
196
        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
197
            if binbounds==None : bins are inferred.
198
199
200
201
202
203
204
205
206
207
        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
208

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

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

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

227
        """
Theo Steininger's avatar
Theo Steininger committed
228

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

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

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

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

        parts = [abs(part)**2 for part in parts]
255
256

        for space_index in spaces:
257
258
            parts = [self._single_power_analyze(
                                work_field=part,
259
                                space_index=space_index,
260
261
                                binbounds=binbounds)
                     for part in parts]
262

263
        if keep_phase_information:
Martin Reinecke's avatar
updates    
Martin Reinecke committed
264
            return parts[0] + 1j*parts[1]
265
        else:
Martin Reinecke's avatar
updates    
Martin Reinecke committed
266
            return parts[0]
267
268

    @classmethod
Martin Reinecke's avatar
Martin Reinecke committed
269
    def _single_power_analyze(cls, work_field, space_index, binbounds):
270

271
        if not work_field.domain[space_index].harmonic:
272
273
            raise ValueError(
                "The analyzed space must be harmonic.")
274

275
276
277
278
279
280
        # 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.

281
        harmonic_domain = work_field.domain[space_index]
282
        power_domain = PowerSpace(harmonic_partner=harmonic_domain,
Theo Steininger's avatar
Theo Steininger committed
283
                                  binbounds=binbounds)
284
285
        power_spectrum = cls._calculate_power_spectrum(
                                field_val=work_field.val,
Martin Reinecke's avatar
Martin Reinecke committed
286
                                pdomain=power_domain,
287
                                axes=work_field.domain_axes[space_index])
288
289

        # create the result field and put power_spectrum into it
290
        result_domain = list(work_field.domain)
291
292
        result_domain[space_index] = power_domain

Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
293
        return Field(domain=result_domain, val=power_spectrum,
Martin Reinecke's avatar
updates    
Martin Reinecke committed
294
                     dtype=power_spectrum.dtype)
295

296
    @classmethod
Martin Reinecke's avatar
Martin Reinecke committed
297
    def _calculate_power_spectrum(cls, field_val, pdomain, axes=None):
298

Martin Reinecke's avatar
Martin Reinecke committed
299
        pindex = pdomain.pindex
300
        if axes is not None:
301
302
303
304
            pindex = cls._shape_up_pindex(
                            pindex=pindex,
                            target_shape=field_val.shape,
                            axes=axes)
Theo Steininger's avatar
Theo Steininger committed
305

Martin Reinecke's avatar
Martin Reinecke committed
306
        power_spectrum = utilities.bincount_axis(pindex, weights=field_val,
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
307
                                                 axis=axes)
Martin Reinecke's avatar
Martin Reinecke committed
308
        rho = pdomain.rho
309
310
311
312
313
314
315
316
        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

317
    @staticmethod
Martin Reinecke's avatar
Martin Reinecke committed
318
    def _shape_up_pindex(pindex, target_shape, axes):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
319
        semiscaled_local_shape = [1] * len(target_shape)
Theo Steininger's avatar
Theo Steininger committed
320
        for i in range(len(axes)):
Martin Reinecke's avatar
Martin Reinecke committed
321
            semiscaled_local_shape[axes[i]] = pindex.shape[i]
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
322
        semiscaled_local_data = pindex.reshape(semiscaled_local_shape)
Martin Reinecke's avatar
Martin Reinecke committed
323
        result_obj = np.empty(target_shape, dtype=pindex.dtype)
324
        result_obj[()] = semiscaled_local_data
325
326
        return result_obj

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

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

334
335
336
        Parameters
        ----------
        spaces : {tuple, int, None} *optional*
Theo Steininger's avatar
Theo Steininger committed
337
338
339
            Specifies the subspace containing all the PowerSpaces which
            should be converted (default : None).
            if spaces==None : Tries to convert the whole domain.
340
        real_power : boolean *optional*
Theo Steininger's avatar
Theo Steininger committed
341
342
            Determines whether the power spectrum is treated as intrinsically
            real or complex (default : True).
343
        real_signal : boolean *optional*
Theo Steininger's avatar
Theo Steininger committed
344
345
346
347
348
349
            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
350
        std : float *optional*
Theo Steininger's avatar
Theo Steininger committed
351
352
353
            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
354

355
356
357
358
        Returns
        -------
        out : Field
            The output object. A random field created with the power spectrum
Theo Steininger's avatar
Theo Steininger committed
359
            stored in the `spaces` in `self`.
360

Theo Steininger's avatar
Theo Steininger committed
361
362
363
364
365
366
        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.

367
368
369
        See Also
        --------
        power_analyze
Theo Steininger's avatar
Theo Steininger committed
370
371
372
373
374

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

375
        """
Theo Steininger's avatar
Theo Steininger committed
376

377
378
        # check if the `spaces` input is valid
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
Theo Steininger's avatar
Theo Steininger committed
379
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
380
            spaces = list(range(len(self.domain)))
Theo Steininger's avatar
Theo Steininger committed
381

382
383
384
385
386
        for power_space_index in spaces:
            power_space = self.domain[power_space_index]
            if not isinstance(power_space, PowerSpace):
                raise ValueError("A PowerSpace is needed for field "
                                 "synthetization.")
387
388
389

        # create the result domain
        result_domain = list(self.domain)
390
391
        for power_space_index in spaces:
            power_space = self.domain[power_space_index]
392
            harmonic_domain = power_space.harmonic_partner
393
            result_domain[power_space_index] = harmonic_domain
394
395
396

        # create random samples: one or two, depending on whether the
        # power spectrum is real or complex
397
398
        result_list = [self.__class__.from_random(
                             'normal',
399
400
401
                             mean=mean,
                             std=std,
                             domain=result_domain,
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
402
                             dtype=np.complex)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
403
                       for x in range(1 if real_power else 2)]
404
405
406
407
408

        # 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
409
        spec = np.sqrt(self.val)
410
        for power_space_index in spaces:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
411
            spec = self._spec_to_rescaler(spec, power_space_index)
412
413

        # apply the rescaler to the random fields
Martin Reinecke's avatar
Martin Reinecke committed
414
        result_list[0] *= spec.real
415
        if not real_power:
Martin Reinecke's avatar
Martin Reinecke committed
416
            result_list[1] *= spec.imag
417

418
        if real_signal:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
419
            result_list = [Field(i.domain, self._hermitian_decomposition(
Martin Reinecke's avatar
Martin Reinecke committed
420
421
                                     i.val,
                                     preserve_gaussian_variance=True)[0])
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
422
                           for i in result_list]
423
424
425

        if real_power:
            result = result_list[0]
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
426
            if not issubclass(result_list[0].dtype.type,
427
428
                              np.complexfloating):
                result = result.real
429
        else:
430
431
432
433
            result = result_list[0] + 1j*result_list[1]

        return result

434
    @staticmethod
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
435
    def _hermitian_decomposition(val, preserve_gaussian_variance=False):
436
        if preserve_gaussian_variance:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
437
438
439
440
441
            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())
442

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
443
    def _spec_to_rescaler(self, spec, power_space_index):
444
        power_space = self.domain[power_space_index]
445

446
447
448
449
        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)
450
        local_blow_up[index] = power_space.pindex
451
        # here, the power_spectrum is distributed into the new shape
452
        return spec[local_blow_up]
453

Theo Steininger's avatar
Theo Steininger committed
454
    # ---Properties---
455

Theo Steininger's avatar
Theo Steininger committed
456
457
    @property
    def val(self):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
458
        """ Returns the data object associated with this Field.
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
459
        No copy is made.
Theo Steininger's avatar
Theo Steininger committed
460

461
462
        Returns
        -------
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
463
        out : numpy.ndarray
464
        """
Martin Reinecke's avatar
Martin Reinecke committed
465
        return self._val
csongor's avatar
csongor committed
466

Martin Reinecke's avatar
Martin Reinecke committed
467
468
469
470
    @property
    def dtype(self):
        return self._val.dtype

471
472
    @property
    def shape(self):
Theo Steininger's avatar
Theo Steininger committed
473
        """ Returns the total shape of the Field's data array.
Theo Steininger's avatar
Theo Steininger committed
474

475
476
477
        Returns
        -------
        out : tuple
Martin Reinecke's avatar
Martin Reinecke committed
478
            The output object. The tuple contains the dimensions of the spaces
479
            in domain.
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
480
       """
Martin Reinecke's avatar
Martin Reinecke committed
481
        return self._val.shape
csongor's avatar
csongor committed
482

483
484
    @property
    def dim(self):
Theo Steininger's avatar
Theo Steininger committed
485
        """ Returns the total number of pixel-dimensions the field has.
Theo Steininger's avatar
Theo Steininger committed
486

Theo Steininger's avatar
Theo Steininger committed
487
        Effectively, all values from shape are multiplied.
Theo Steininger's avatar
Theo Steininger committed
488

489
490
491
492
493
        Returns
        -------
        out : int
            The dimension of the Field.
        """
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
494
        return self._val.size
csongor's avatar
csongor committed
495

Theo Steininger's avatar
Theo Steininger committed
496
497
    @property
    def total_volume(self):
Theo Steininger's avatar
Theo Steininger committed
498
499
        """ Returns the total volume of all spaces in the domain.
        """
Theo Steininger's avatar
Theo Steininger committed
500
        volume_tuple = tuple(sp.total_volume for sp in self.domain)
501
        try:
Theo Steininger's avatar
Theo Steininger committed
502
            return reduce(lambda x, y: x * y, volume_tuple)
503
        except TypeError:
Theo Steininger's avatar
Theo Steininger committed
504
            return 0.
505

Theo Steininger's avatar
Theo Steininger committed
506
507
508
509
    @property
    def real(self):
        """ The real part of the field (data is not copied).
        """
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
510
        return Field(self.domain, self.val.real)
Theo Steininger's avatar
Theo Steininger committed
511
512
513
514
515

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

Theo Steininger's avatar
Theo Steininger committed
518
    # ---Special unary/binary operations---
519

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
520
    def copy(self, domain=None, dtype=None):
521
        """ Returns a full copy of the Field.
Theo Steininger's avatar
Theo Steininger committed
522

523
524
        If no keyword arguments are given, the returned object will be an
        identical copy of the original Field. By explicit specification one is
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
525
        able to define the domain and the dtype of the returned Field.
526
527
528
529
530

        Parameters
        ----------
        domain : DomainObject
            The new domain the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
531

532
533
        dtype : type
            The new dtype the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
534

535
536
537
538
539
540
        Returns
        -------
        out : Field
            The output object. An identical copy of 'self'.

        """
Theo Steininger's avatar
Theo Steininger committed
541

542
543
        if domain is None:
            domain = self.domain
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
544
        return Field(domain=domain, val=self._val, dtype=dtype, copy=True)
csongor's avatar
csongor committed
545

546
547
548
549
550
551
    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
552
        res = 1.
553
554
555
556
557
558
559
        for i in spaces:
            tmp = self.domain[i].scalar_weight()
            if tmp is None:
                return None
            res *= tmp
        return res

Theo Steininger's avatar
Theo Steininger committed
560
    def weight(self, power=1, inplace=False, spaces=None):
Theo Steininger's avatar
Theo Steininger committed
561
        """ Weights the pixels of `self` with their invidual pixel-volume.
562
563
564
565

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

568
        inplace : boolean
Theo Steininger's avatar
Theo Steininger committed
569
570
            If True, `self` will be weighted and returned. Otherwise, a copy
            is made.
Theo Steininger's avatar
Theo Steininger committed
571

Theo Steininger's avatar
Theo Steininger committed
572
573
        spaces : tuple of ints
            Determines on which subspace the operation takes place.
Theo Steininger's avatar
Theo Steininger committed
574

575
576
577
        Returns
        -------
        out : Field
Theo Steininger's avatar
Theo Steininger committed
578
            The weighted field.
579
580

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

583
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
csongor's avatar
csongor committed
584
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
585
            spaces = list(range(len(self.domain)))
csongor's avatar
csongor committed
586

587
588
589
590
591
592
593
594
595
        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)
596
                new_field *= wgt**power
597
        fct = fct**power
Martin Reinecke's avatar
Martin Reinecke committed
598
        if fct != 1.:
599
            new_field *= fct
600

601
        return new_field
csongor's avatar
csongor committed
602

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

606
607
608
        Parameters
        ----------
        x : Field
Theo Steininger's avatar
Theo Steininger committed
609
            The domain of x must contain `self.domain`
Theo Steininger's avatar
Theo Steininger committed
610

Theo Steininger's avatar
Theo Steininger committed
611
612
613
        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
614

615
        bare : boolean
Theo Steininger's avatar
Theo Steininger committed
616
            If true, no volume factors will be included in the computation.
Theo Steininger's avatar
Theo Steininger committed
617

618
619
620
        Returns
        -------
        out : float, complex
Theo Steininger's avatar
Theo Steininger committed
621

622
        """
623
624
625
        if not isinstance(x, Field):
            raise ValueError("The dot-partner must be an instance of " +
                             "the NIFTy field class")
Theo Steininger's avatar
Theo Steininger committed
626

Martin Reinecke's avatar
Martin Reinecke committed
627
        # Compute the dot respecting the fact of discrete/continuous spaces
628
629
630
631
632
633
634
635
636
637
        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
Theo Steininger's avatar
Theo Steininger committed
638

639
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
640
            return fct*np.vdot(y.val.ravel(), x.val.ravel())
641
642
643
        else:
            # create a diagonal operator which is capable of taking care of the
            # axes-matching
644
            from .operators.diagonal_operator import DiagonalOperator
645
646
647
648
649
            diagonal = y.val.conjugate()
            diagonalOperator = DiagonalOperator(domain=y.domain,
                                                diagonal=diagonal,
                                                copy=False)
            dotted = diagonalOperator(x, spaces=spaces)
650
            return fct*dotted.sum(spaces=spaces)
Theo Steininger's avatar
Theo Steininger committed
651

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

Theo Steininger's avatar
Theo Steininger committed
655
656
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
657
        norm : float
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
658
            The L2-norm of the field values.
csongor's avatar
csongor committed
659
660

        """
Martin Reinecke's avatar
fixes    
Martin Reinecke committed
661
        return np.sqrt(np.abs(self.vdot(x=self)))
csongor's avatar
csongor committed
662
663

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

666
667
668
        Parameters
        ----------
        inplace : boolean
Theo Steininger's avatar
Theo Steininger committed
669
            Decides whether the conjugation should be performed inplace.
Theo Steininger's avatar
Theo Steininger committed
670

671
672
673
674
        Returns
        -------
        cc : field
            The complex conjugated field.
csongor's avatar
csongor committed
675
676
677

        """
        if inplace:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
678
            self.imag *= -1
679
            return self
csongor's avatar
csongor committed
680
        else:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
681
            return Field(self.domain, np.conj(self.val), self.dtype)
csongor's avatar
csongor committed
682

Theo Steininger's avatar
Theo Steininger committed
683
    # ---General unary/contraction methods---
684

Theo Steininger's avatar
Theo Steininger committed
685
686
    def __pos__(self):
        return self.copy()
687

Theo Steininger's avatar
Theo Steininger committed
688
    def __neg__(self):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
689
        return Field(self.domain, -self.val, self.dtype)
csongor's avatar
csongor committed
690

Theo Steininger's avatar
Theo Steininger committed
691
    def __abs__(self):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
692
        return Field(self.domain, np.abs(self.val), self.dtype)
csongor's avatar
csongor committed
693

694
    def _contraction_helper(self, op, spaces):
Theo Steininger's avatar
Theo Steininger committed
695
        if spaces is None:
696
697
698
            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
699

700
        axes_list = tuple(self.domain_axes[sp_index] for sp_index in spaces)
701
702

        try:
Theo Steininger's avatar
Theo Steininger committed
703
            axes_list = reduce(lambda x, y: x+y, axes_list)
704
        except TypeError:
Theo Steininger's avatar
Theo Steininger committed
705
            axes_list = ()
csongor's avatar
csongor committed
706

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

Theo Steininger's avatar
Theo Steininger committed
710
711
712
        # 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
713
        else:
Theo Steininger's avatar
Theo Steininger committed
714
            return_domain = tuple(self.domain[i]
Martin Reinecke's avatar
Martin Reinecke committed
715
                                  for i in range(len(self.domain))
Theo Steininger's avatar
Theo Steininger committed
716
                                  if i not in spaces)
717

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

720
721
    def sum(self, spaces=None):
        return self._contraction_helper('sum', spaces)
csongor's avatar
csongor committed
722

723
724
    def prod(self, spaces=None):
        return self._contraction_helper('prod', spaces)
csongor's avatar
csongor committed
725

726
727
    def all(self, spaces=None):
        return self._contraction_helper('all', spaces)
csongor's avatar
csongor committed
728

729
730
    def any(self, spaces=None):
        return self._contraction_helper('any', spaces)
csongor's avatar
csongor committed
731

732
733
    def min(self, spaces=None):
        return self._contraction_helper('min', spaces)
csongor's avatar
csongor committed
734

735
736
    def max(self, spaces=None):
        return self._contraction_helper('max', spaces)
csongor's avatar
csongor committed
737

738
739
    def mean(self, spaces=None):
        return self._contraction_helper('mean', spaces)
csongor's avatar
csongor committed
740

741
742
    def var(self, spaces=None):
        return self._contraction_helper('var', spaces)
csongor's avatar
csongor committed
743

744
745
    def std(self, spaces=None):
        return self._contraction_helper('std', spaces)
csongor's avatar
csongor committed
746

Theo Steininger's avatar
Theo Steininger committed
747
    # ---General binary methods---
csongor's avatar
csongor committed
748

749
    def _binary_helper(self, other, op):
csongor's avatar
csongor committed
750
        # if other is a field, make sure that the domains match
751
        if isinstance(other, Field):
752
753
            if other.domain != self.domain:
                raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
754
            return Field(self.domain, getattr(self.val, op)(other.val))
csongor's avatar
csongor committed
755

Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
756
        return Field(self.domain, getattr(self.val, op)(other))
csongor's avatar
csongor committed
757
758

    def __add__(self, other):
Theo Steininger's avatar
Theo Steininger committed
759
        return self._binary_helper(other, op='__add__')
760

761
    def __radd__(self, other):
Theo Steininger's avatar
Theo Steininger committed
762
        return self._binary_helper(other, op='__radd__')
csongor's avatar
csongor committed
763
764

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

    def __sub__(self, other):
Theo Steininger's avatar
Theo Steininger committed
768
        return self._binary_helper(other, op='__sub__')
csongor's avatar
csongor committed
769
770

    def __rsub__(self, other):
Theo Steininger's avatar
Theo Steininger committed
771
        return self._binary_helper(other, op='__rsub__')
csongor's avatar
csongor committed
772
773

    def __isub__(self, other):
774
        return self._binary_helper(other, op='__isub__')
csongor's avatar
csongor committed
775
776

    def __mul__(self, other):
Theo Steininger's avatar
Theo Steininger committed
777
        return self._binary_helper(other, op='__mul__')
778

779
    def __rmul__(self, other):
Theo Steininger's avatar
Theo Steininger committed
780
        return self._binary_helper(other, op='__rmul__')
csongor's avatar
csongor committed
781
782

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

    def __div__(self, other):
Theo Steininger's avatar
Theo Steininger committed
786
        return self._binary_helper(other, op='__div__')
csongor's avatar
csongor committed
787

Martin Reinecke's avatar
Martin Reinecke committed
788
789
790
    def __truediv__(self, other):
        return self._binary_helper(other, op='__truediv__')

csongor's avatar
csongor committed
791
    def __rdiv__(self, other):
Theo Steininger's avatar
Theo Steininger committed
792
        return self._binary_helper(other, op='__rdiv__')
csongor's avatar
csongor committed
793

Martin Reinecke's avatar
Martin Reinecke committed
794
795
796
    def __rtruediv__(self, other):
        return self._binary_helper(other, op='__rtruediv__')

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

csongor's avatar
csongor committed
800
    def __pow__(self, other):
Theo Steininger's avatar
Theo Steininger committed
801
        return self._binary_helper(other, op='__pow__')
csongor's avatar
csongor committed
802
803

    def __rpow__(self, other):
Theo Steininger's avatar
Theo Steininger committed
804
        return self._binary_helper(other, op='__rpow__')
csongor's avatar
csongor committed
805
806

    def __ipow__(self, other):
807
        return self._binary_helper(other, op='__ipow__')
csongor's avatar
csongor committed
808
809

    def __lt__(self, other):
Theo Steininger's avatar
Theo Steininger committed
810
        return self._binary_helper(other, op='__lt__')
csongor's avatar
csongor committed
811
812

    def __le__(self, other):
Theo Steininger's avatar
Theo Steininger committed
813
        return self._binary_helper(other, op='__le__')
csongor's avatar
csongor committed
814
815
816
817
818

    def __ne__(self, other):
        if other is None:
            return True
        else:
Theo Steininger's avatar
Theo Steininger committed
819
            return self._binary_helper(other