field.py 31.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
21
from builtins import zip
from builtins import range
22

23
import ast
csongor's avatar
csongor committed
24
25
import numpy as np

Martin Reinecke's avatar
Martin Reinecke committed
26
from .domain_object import DomainObject
27

Martin Reinecke's avatar
Martin Reinecke committed
28
from .spaces.power_space import PowerSpace
csongor's avatar
csongor committed
29

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

csongor's avatar
csongor committed
34

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

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

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

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
49
    val : scalar, numpy.ndarray, Field
50
51
52
        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
53

54
55
    dtype : type
        A numpy.type. Most common are int, float and complex.
Theo Steininger's avatar
Theo Steininger committed
56

57
58
59
60
    copy: boolean

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

63
64
65
66
67
68
    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
69

70
71
72
73
74
75
76
    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
77

78
    """
79

Theo Steininger's avatar
Theo Steininger committed
80
    # ---Initialization methods---
81

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

Martin Reinecke's avatar
Martin Reinecke committed
86
87
88
        shape_tuple = tuple(sp.shape for sp in self.domain)
        if len(shape_tuple)==0:
            global_shape = ()
89
        else:
Martin Reinecke's avatar
Martin Reinecke committed
90
91
92
93
            global_shape = reduce(lambda x, y: x + y, shape_tuple)
        dtype = self._infer_dtype(dtype=dtype, val=val)
        self._val = np.empty(global_shape,dtype=dtype)
        self.set_val(new_val=val, copy=copy)
csongor's avatar
csongor committed
94

95
    def _parse_domain(self, domain, val=None):
96
        if domain is None:
97
98
            if isinstance(val, Field):
                domain = val.domain
Martin Reinecke's avatar
Martin Reinecke committed
99
            elif np.isscalar(val):
100
                domain = ()
Martin Reinecke's avatar
Martin Reinecke committed
101
102
            else:
                raise TypeError("could not infer domain from value")
103
        elif isinstance(domain, DomainObject):
104
            domain = (domain,)
105
106
107
        elif not isinstance(domain, tuple):
            domain = tuple(domain)

csongor's avatar
csongor committed
108
        for d in domain:
109
            if not isinstance(d, DomainObject):
110
111
                raise TypeError(
                    "Given domain contains something that is not a "
112
                    "DomainObject instance.")
csongor's avatar
csongor committed
113
114
        return domain

115
116
    def _get_axes_tuple(self, things_with_shape):
        i = 0
Theo Steininger's avatar
Theo Steininger committed
117
118
        axes_list = []
        for thing in things_with_shape:
119
120
121
            nax = len(thing.shape)
            axes_list += [tuple(range(i,i+nax))]
            i += nax
Theo Steininger's avatar
Theo Steininger committed
122
        return tuple(axes_list)
123

124
    def _infer_dtype(self, dtype, val):
csongor's avatar
csongor committed
125
        if dtype is None:
126
            try:
127
                dtype = val.dtype
128
            except AttributeError:
Theo Steininger's avatar
Theo Steininger committed
129
130
131
                try:
                    if val is None:
                        raise TypeError
132
                    dtype = np.result_type(val)
Theo Steininger's avatar
Theo Steininger committed
133
                except(TypeError):
Martin Reinecke's avatar
more    
Martin Reinecke committed
134
                    dtype = np.dtype(np.float64)
Theo Steininger's avatar
Theo Steininger committed
135
        else:
136
            dtype = np.dtype(dtype)
137

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
138
        return np.result_type(dtype, np.float)
139

140
    # ---Factory methods---
141

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

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

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

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

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

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

        See Also
        --------
167
        power_synthesize
Theo Steininger's avatar
Theo Steininger committed
168

169
170

        """
Theo Steininger's avatar
Theo Steininger committed
171

172
        # create a initially empty field
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
173
        f = cls(domain=domain, dtype=dtype)
174

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
175
        # extract the data from f and apply the appropriate
176
177
178
        # random number generator to it
        sample = f.get_val(copy=False)
        generator_function = getattr(Random, random_type)
179

180
        sample[()]=generator_function(dtype=f.dtype,
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
181
                                             shape=sample.shape,
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
182
                                             **kwargs)
183
184
        return f

185
186
    # ---Powerspectral methods---

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

Theo Steininger's avatar
Theo Steininger committed
191
192
193
        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
194
        harmonic space. The resulting field has the same units as the initial
Theo Steininger's avatar
Theo Steininger committed
195
        field, corresponding to the square root of the power spectrum.
196
197
198

        Parameters
        ----------
Theo Steininger's avatar
Theo Steininger committed
199
200
201
202
203
        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
204
            if binbounds==None : bins are inferred.
205
206
207
208
209
210
211
212
213
214
        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
215

216
217
218
219
        Raise
        -----
        ValueError
            Raised if
Theo Steininger's avatar
Theo Steininger committed
220
221
                *len(domain) is != 1 when spaces==None
                *len(spaces) is != 1 if not None
222
                *the analyzed space is not harmonic
Theo Steininger's avatar
Theo Steininger committed
223

224
225
        Returns
        -------
Theo Steininger's avatar
Theo Steininger committed
226
        out : Field
Martin Reinecke's avatar
Martin Reinecke committed
227
            The output object. Its domain is a PowerSpace and it contains
228
229
230
231
232
            the power spectrum of 'self's field.

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

234
        """
Theo Steininger's avatar
Theo Steininger committed
235

Theo Steininger's avatar
Theo Steininger committed
236
        # check if all spaces in `self.domain` are either harmonic or
237
238
239
        # power_space instances
        for sp in self.domain:
            if not sp.harmonic and not isinstance(sp, PowerSpace):
Martin Reinecke's avatar
Martin Reinecke committed
240
                raise TypeError(
241
                    "Field has a space in `domain` which is neither "
242
243
244
                    "harmonic nor a PowerSpace.")

        # check if the `spaces` input is valid
245
246
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
247
            spaces = list(range(len(self.domain)))
248
249

        if len(spaces) == 0:
250
251
            raise ValueError(
                "No space for analysis specified.")
252

253
254
255
256
257
258
259
260
261
262
263
264
265
        if keep_phase_information:
            parts_val = self._hermitian_decomposition(
                                              domain=self.domain,
                                              val=self.val,
                                              spaces=spaces,
                                              domain_axes=self.domain_axes,
                                              preserve_gaussian_variance=False)
            parts = [self.copy_empty().set_val(part_val, copy=False)
                     for part_val in parts_val]
        else:
            parts = [self]

        parts = [abs(part)**2 for part in parts]
266
267

        for space_index in spaces:
268
269
            parts = [self._single_power_analyze(
                                work_field=part,
270
                                space_index=space_index,
271
272
                                binbounds=binbounds)
                     for part in parts]
273

274
275
276
277
278
279
        if keep_phase_information:
            result_field = parts[0] + 1j*parts[1]
        else:
            result_field = parts[0]

        return result_field
280
281

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

284
        if not work_field.domain[space_index].harmonic:
285
286
            raise ValueError(
                "The analyzed space must be harmonic.")
287

288
289
290
291
292
293
        # 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.

294
        harmonic_domain = work_field.domain[space_index]
295
        power_domain = PowerSpace(harmonic_partner=harmonic_domain,
Theo Steininger's avatar
Theo Steininger committed
296
                                  binbounds=binbounds)
297
298
        power_spectrum = cls._calculate_power_spectrum(
                                field_val=work_field.val,
Martin Reinecke's avatar
Martin Reinecke committed
299
                                pdomain=power_domain,
300
                                axes=work_field.domain_axes[space_index])
301
302

        # create the result field and put power_spectrum into it
303
        result_domain = list(work_field.domain)
304
        result_domain[space_index] = power_domain
305
        result_dtype = power_spectrum.dtype
306

307
        result_field = work_field.copy_empty(
308
                   domain=result_domain,
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
309
                   dtype=result_dtype)
310
311
312
313
        result_field.set_val(new_val=power_spectrum, copy=False)

        return result_field

314
    @classmethod
Martin Reinecke's avatar
Martin Reinecke committed
315
    def _calculate_power_spectrum(cls, field_val, pdomain, axes=None):
316

Martin Reinecke's avatar
Martin Reinecke committed
317
318
319
        pindex = pdomain.pindex
        # MR FIXME: how about iterating over slices, instead of replicating
        # pindex? Would save memory and probably isn't slower.
320
        if axes is not None:
321
322
323
324
            pindex = cls._shape_up_pindex(
                            pindex=pindex,
                            target_shape=field_val.shape,
                            axes=axes)
Theo Steininger's avatar
Theo Steininger committed
325

Martin Reinecke's avatar
Martin Reinecke committed
326
        power_spectrum = utilities.bincount_axis(pindex, weights=field_val,
327
                                         axis=axes)
Martin Reinecke's avatar
Martin Reinecke committed
328
        rho = pdomain.rho
329
330
331
332
333
334
335
336
        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

337
    @staticmethod
Martin Reinecke's avatar
Martin Reinecke committed
338
    def _shape_up_pindex(pindex, target_shape, axes):
Theo Steininger's avatar
Theo Steininger committed
339
        semiscaled_local_shape = [1, ] * len(target_shape)
Theo Steininger's avatar
Theo Steininger committed
340
        for i in range(len(axes)):
Martin Reinecke's avatar
Martin Reinecke committed
341
342
            semiscaled_local_shape[axes[i]] = pindex.shape[i]
        local_data = pindex
Theo Steininger's avatar
Theo Steininger committed
343
        semiscaled_local_data = local_data.reshape(semiscaled_local_shape)
Martin Reinecke's avatar
Martin Reinecke committed
344
        result_obj = np.empty(target_shape, dtype=pindex.dtype)
345
        result_obj[()] = semiscaled_local_data
346
347
348

        return result_obj

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

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

356
357
358
        Parameters
        ----------
        spaces : {tuple, int, None} *optional*
Theo Steininger's avatar
Theo Steininger committed
359
360
361
            Specifies the subspace containing all the PowerSpaces which
            should be converted (default : None).
            if spaces==None : Tries to convert the whole domain.
362
        real_power : boolean *optional*
Theo Steininger's avatar
Theo Steininger committed
363
364
            Determines whether the power spectrum is treated as intrinsically
            real or complex (default : True).
365
        real_signal : boolean *optional*
Theo Steininger's avatar
Theo Steininger committed
366
367
368
369
370
371
            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
372
        std : float *optional*
Theo Steininger's avatar
Theo Steininger committed
373
374
375
            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
376

377
378
379
380
        Returns
        -------
        out : Field
            The output object. A random field created with the power spectrum
Theo Steininger's avatar
Theo Steininger committed
381
            stored in the `spaces` in `self`.
382

Theo Steininger's avatar
Theo Steininger committed
383
384
385
386
387
388
        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.

389
390
391
        See Also
        --------
        power_analyze
Theo Steininger's avatar
Theo Steininger committed
392
393
394
395
396

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

397
        """
Theo Steininger's avatar
Theo Steininger committed
398

399
400
401
        # check if the `spaces` input is valid
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))

Theo Steininger's avatar
Theo Steininger committed
402
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
403
            spaces = list(range(len(self.domain)))
Theo Steininger's avatar
Theo Steininger committed
404

405
406
407
408
409
        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.")
410
411
412

        # create the result domain
        result_domain = list(self.domain)
413
414
        for power_space_index in spaces:
            power_space = self.domain[power_space_index]
415
            harmonic_domain = power_space.harmonic_partner
416
            result_domain[power_space_index] = harmonic_domain
417
418
419

        # create random samples: one or two, depending on whether the
        # power spectrum is real or complex
420
        if real_power:
421
            result_list = [None]
422
423
        else:
            result_list = [None, None]
424

425
426
        result_list = [self.__class__.from_random(
                             'normal',
427
428
429
                             mean=mean,
                             std=std,
                             domain=result_domain,
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
430
                             dtype=np.complex)
431
432
433
434
435
436
                       for x in result_list]

        # 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
437

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
438
        spec = self.val.copy()
439
440
        spec = np.sqrt(spec)

441
442
443
444
445
446
447
        for power_space_index in spaces:
            spec = self._spec_to_rescaler(spec, result_list, power_space_index)
        local_rescaler = spec

        result_val_list = [x.val for x in result_list]

        # apply the rescaler to the random fields
Martin Reinecke's avatar
Martin Reinecke committed
448
        result_val_list[0] *= local_rescaler.real
449
450

        if not real_power:
Martin Reinecke's avatar
Martin Reinecke committed
451
            result_val_list[1] *= local_rescaler.imag
452

453
        if real_signal:
454
            result_val_list = [self._hermitian_decomposition(
455
456
457
458
459
                                            result_domain,
                                            result_val,
                                            spaces,
                                            result_list[0].domain_axes,
                                            preserve_gaussian_variance=True)[0]
460
                               for result_val in result_val_list]
461
462
463
464
465
466
467

        # store the result into the fields
        [x.set_val(new_val=y, copy=False) for x, y in
            zip(result_list, result_val_list)]

        if real_power:
            result = result_list[0]
468
469
470
            if not issubclass(result_val_list[0].dtype.type,
                              np.complexfloating):
                result = result.real
471
        else:
472
473
474
475
            result = result_list[0] + 1j*result_list[1]

        return result

476
    @staticmethod
477
478
    def _hermitian_decomposition(domain, val, spaces, domain_axes,
                                 preserve_gaussian_variance=False):
479

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
480
481
        h = val.real.copy()
        a = 1j * val.imag.copy()
482
483

        # correct variance
484
        if preserve_gaussian_variance:
Martin Reinecke's avatar
Martin Reinecke committed
485
486
            assert issubclass(val.dtype.type, np.complexfloating),\
                    "complex input field is needed here"
487
488
489
            h *= np.sqrt(2)
            a *= np.sqrt(2)

490
491
        return (h, a)

492
493
    def _spec_to_rescaler(self, spec, result_list, power_space_index):
        power_space = self.domain[power_space_index]
494

495
496
497
498
        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)
499
        local_blow_up[index] = power_space.pindex
500
        # here, the power_spectrum is distributed into the new shape
501
        return spec[local_blow_up]
502

Theo Steininger's avatar
Theo Steininger committed
503
    # ---Properties---
504

Theo Steininger's avatar
Theo Steininger committed
505
    def set_val(self, new_val=None, copy=False):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
506
        """ Sets the field's data object.
507
508
509

        Parameters
        ----------
Theo Steininger's avatar
Theo Steininger committed
510
        new_val : scalar, array-like, Field, None *optional*
511
512
            The values to be stored in the field.
            {default : None}
Theo Steininger's avatar
Theo Steininger committed
513

514
        copy : boolean, *optional*
Theo Steininger's avatar
Theo Steininger committed
515
516
            If False, Field tries to not copy the input data but use it
            directly.
517
518
519
520
521
522
            {default : False}
        See Also
        --------
        val

        """
Theo Steininger's avatar
Theo Steininger committed
523

Martin Reinecke's avatar
Martin Reinecke committed
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
        if new_val is None:
            pass
        elif isinstance(new_val, Field):
            if self.domain!=new_val.domain:
                raise ValueError("Domain mismatch")
            if copy:
                self._val[()] = new_val.val
            else:
                self._val = new_val.val
        elif (np.isscalar(new_val)):
            self._val[()]=new_val
        elif isinstance(new_val, np.ndarray):
            if copy:
                self._val[()] = new_val
            else:
                if self.shape!=new_val.shape:
                    raise ValueError("Shape mismatch")
                self._val = new_val
        else:
            raise TypeError("unknown source type")
544
        return self
csongor's avatar
csongor committed
545

546
    def get_val(self, copy=False):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
547
        """ Returns the data object associated with this Field.
548
549
550
551

        Parameters
        ----------
        copy : boolean
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
552
            If true, a copy of the Field's underlying data object
Theo Steininger's avatar
Theo Steininger committed
553
            is returned.
Theo Steininger's avatar
Theo Steininger committed
554

555
556
        Returns
        -------
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
557
        out : numpy.ndarray
558
559
560
561
562
563

        See Also
        --------
        val

        """
Theo Steininger's avatar
Theo Steininger committed
564

Martin Reinecke's avatar
Martin Reinecke committed
565
        return self._val.copy() if copy else self._val
csongor's avatar
csongor committed
566

Theo Steininger's avatar
Theo Steininger committed
567
568
    @property
    def val(self):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
569
        """ Returns the data object associated with this Field.
Theo Steininger's avatar
Theo Steininger committed
570

571
572
        Returns
        -------
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
573
        out : numpy.ndarray
574
575
576
577
578
579

        See Also
        --------
        get_val

        """
Theo Steininger's avatar
Theo Steininger committed
580

Martin Reinecke's avatar
Martin Reinecke committed
581
        return self._val
csongor's avatar
csongor committed
582

Theo Steininger's avatar
Theo Steininger committed
583
584
    @val.setter
    def val(self, new_val):
585
        self.set_val(new_val=new_val, copy=False)
csongor's avatar
csongor committed
586

Martin Reinecke's avatar
Martin Reinecke committed
587
588
589
590
    @property
    def dtype(self):
        return self._val.dtype

591
592
    @property
    def shape(self):
Theo Steininger's avatar
Theo Steininger committed
593
        """ Returns the total shape of the Field's data array.
Theo Steininger's avatar
Theo Steininger committed
594

595
596
597
        Returns
        -------
        out : tuple
Martin Reinecke's avatar
Martin Reinecke committed
598
            The output object. The tuple contains the dimensions of the spaces
599
600
601
602
603
604
605
            in domain.

        See Also
        --------
        dim

        """
Martin Reinecke's avatar
Martin Reinecke committed
606
        return self._val.shape
csongor's avatar
csongor committed
607

608
609
    @property
    def dim(self):
Theo Steininger's avatar
Theo Steininger committed
610
        """ Returns the total number of pixel-dimensions the field has.
Theo Steininger's avatar
Theo Steininger committed
611

Theo Steininger's avatar
Theo Steininger committed
612
        Effectively, all values from shape are multiplied.
Theo Steininger's avatar
Theo Steininger committed
613

614
615
616
617
618
619
620
621
622
623
        Returns
        -------
        out : int
            The dimension of the Field.

        See Also
        --------
        shape

        """
Theo Steininger's avatar
Theo Steininger committed
624

625
        dim_tuple = tuple(sp.dim for sp in self.domain)
Theo Steininger's avatar
Theo Steininger committed
626
        try:
Martin Reinecke's avatar
Martin Reinecke committed
627
            return int(reduce(lambda x, y: x * y, dim_tuple))
Theo Steininger's avatar
Theo Steininger committed
628
629
        except TypeError:
            return 0
csongor's avatar
csongor committed
630

631
632
    @property
    def dof(self):
Theo Steininger's avatar
Theo Steininger committed
633
634
635
636
637
638
        """ Returns the total number of degrees of freedom the Field has. For
        real Fields this is equal to `self.dim`. For complex Fields it is
        2*`self.dim`.

        """

Theo Steininger's avatar
Theo Steininger committed
639
640
641
642
643
644
645
        dof = self.dim
        if issubclass(self.dtype.type, np.complexfloating):
            dof *= 2
        return dof

    @property
    def total_volume(self):
Theo Steininger's avatar
Theo Steininger committed
646
647
648
        """ Returns the total volume of all spaces in the domain.
        """

Theo Steininger's avatar
Theo Steininger committed
649
        volume_tuple = tuple(sp.total_volume for sp in self.domain)
650
        try:
Theo Steininger's avatar
Theo Steininger committed
651
            return reduce(lambda x, y: x * y, volume_tuple)
652
        except TypeError:
Theo Steininger's avatar
Theo Steininger committed
653
            return 0.
654

Theo Steininger's avatar
Theo Steininger committed
655
656
657
658
    @property
    def real(self):
        """ The real part of the field (data is not copied).
        """
659
        return Field(self.domain,self.val.real)
Theo Steininger's avatar
Theo Steininger committed
660
661
662
663
664

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

Theo Steininger's avatar
Theo Steininger committed
667
    # ---Special unary/binary operations---
668

csongor's avatar
csongor committed
669

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

673
674
        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
675
        able to define the domain and the dtype of the returned Field.
676
677
678
679
680

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

682
683
        dtype : type
            The new dtype the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
684

685
686
687
688
689
690
691
692
693
694
        Returns
        -------
        out : Field
            The output object. An identical copy of 'self'.

        See Also
        --------
        copy_empty

        """
Theo Steininger's avatar
Theo Steininger committed
695

696
697
698
        if domain is None:
            domain = self.domain
        return Field(domain=domain,val=self._val,dtype=dtype,copy=True)
csongor's avatar
csongor committed
699

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
700
    def copy_empty(self, domain=None, dtype=None):
701
702
703
        """ Returns an empty copy of the Field.

        If no keyword arguments are given, the returned object will be an
Theo Steininger's avatar
Theo Steininger committed
704
705
706
        identical copy of the original Field. The memory for the data array
        is only allocated but not actively set to any value
        (c.f. numpy.ndarray.copy_empty). By explicit specification one is able
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
707
        to change the domain and the dtype of the returned Field.
Theo Steininger's avatar
Theo Steininger committed
708

709
710
711
712
        Parameters
        ----------
        domain : DomainObject
            The new domain the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
713

714
715
        dtype : type
            The new dtype the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
716

717
718
719
        Returns
        -------
        out : Field
Theo Steininger's avatar
Theo Steininger committed
720
            The output object.
721
722
723
724
725
726

        See Also
        --------
        copy

        """
Theo Steininger's avatar
Theo Steininger committed
727

Theo Steininger's avatar
Theo Steininger committed
728
729
730
731
        if domain is None:
            domain = self.domain
        if dtype is None:
            dtype = self.dtype
732
        return Field(domain=domain, dtype=dtype)
Theo Steininger's avatar
Theo Steininger committed
733
734

    def weight(self, power=1, inplace=False, spaces=None):
Theo Steininger's avatar
Theo Steininger committed
735
        """ Weights the pixels of `self` with their invidual pixel-volume.
736
737
738
739

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

742
        inplace : boolean
Theo Steininger's avatar
Theo Steininger committed
743
744
            If True, `self` will be weighted and returned. Otherwise, a copy
            is made.
Theo Steininger's avatar
Theo Steininger committed
745

Theo Steininger's avatar
Theo Steininger committed
746
747
        spaces : tuple of ints
            Determines on which subspace the operation takes place.
Theo Steininger's avatar
Theo Steininger committed
748

749
750
751
        Returns
        -------
        out : Field
Theo Steininger's avatar
Theo Steininger committed
752
            The weighted field.
753
754

        """
755
        new_field = self if inplace else self.copy_empty()
csongor's avatar
csongor committed
756

757
        new_val = self.get_val(copy=False)
csongor's avatar
csongor committed
758

759
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
csongor's avatar
csongor committed
760
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
761
            spaces = list(range(len(self.domain)))
csongor's avatar
csongor committed
762

763
        for ind, sp in enumerate(self.domain):
Theo Steininger's avatar
Theo Steininger committed
764
765
766
767
768
            if ind in spaces:
                new_val = sp.weight(new_val,
                                    power=power,
                                    axes=self.domain_axes[ind],
                                    inplace=inplace)
769
770

        new_field.set_val(new_val=new_val, copy=False)
csongor's avatar
csongor committed
771
772
        return new_field

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

776
777
778
        Parameters
        ----------
        x : Field
Theo Steininger's avatar
Theo Steininger committed
779
            The domain of x must contain `self.domain`
Theo Steininger's avatar
Theo Steininger committed
780

Theo Steininger's avatar
Theo Steininger committed
781
782
783
        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
784

785
        bare : boolean
Theo Steininger's avatar
Theo Steininger committed
786
            If true, no volume factors will be included in the computation.
Theo Steininger's avatar
Theo Steininger committed
787

788
789
790
        Returns
        -------
        out : float, complex
Theo Steininger's avatar
Theo Steininger committed
791

792
        """
793
794
795
        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
796

Martin Reinecke's avatar
Martin Reinecke committed
797
        # Compute the dot respecting the fact of discrete/continuous spaces
798
        y = self if bare else self.weight(power=1)
Theo Steininger's avatar
Theo Steininger committed
799

800
        if spaces is None:
801
            return np.vdot(y.val.flatten(),x.val.flatten())
802
803
804
        else:
            # create a diagonal operator which is capable of taking care of the
            # axes-matching
805
            from .operators.diagonal_operator import DiagonalOperator
806
807
808
809
810
811
            diagonal = y.val.conjugate()
            diagonalOperator = DiagonalOperator(domain=y.domain,
                                                diagonal=diagonal,
                                                copy=False)
            dotted = diagonalOperator(x, spaces=spaces)
            return dotted.sum(spaces=spaces)
Theo Steininger's avatar
Theo Steininger committed
812

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

Theo Steininger's avatar
Theo Steininger committed
816
817
818
        Returns
        -------
        norm : scalar
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
819
            The L2-norm of the field values.
csongor's avatar
csongor committed
820
821

        """
Martin Reinecke's avatar
fixes    
Martin Reinecke committed
822
        return np.sqrt(np.abs(self.vdot(x=self)))
csongor's avatar
csongor committed
823
824

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

827
828
829
        Parameters
        ----------
        inplace : boolean
Theo Steininger's avatar
Theo Steininger committed
830
            Decides whether the conjugation should be performed inplace.
Theo Steininger's avatar
Theo Steininger committed
831

832
833
834
835
        Returns
        -------
        cc : field
            The complex conjugated field.
csongor's avatar
csongor committed
836
837

        """
Theo Steininger's avatar
Theo Steininger committed
838

csongor's avatar
csongor committed
839
        if inplace:
840
841
            self.imag*=-1
            return self
csongor's avatar
csongor committed
842
        else:
843
            return Field(self.domain,np.conj(self.val),self.dtype)
csongor's avatar
csongor committed
844

Theo Steininger's avatar
Theo Steininger committed
845
    # ---General unary/contraction methods---
846

Theo Steininger's avatar
Theo Steininger committed
847
848
    def __pos__(self):
        return self.copy()
849

Theo Steininger's avatar
Theo Steininger committed
850
    def __neg__(self):
851
        return Field(self.domain,-self.val,self.dtype)
csongor's avatar
csongor committed
852

Theo Steininger's avatar
Theo Steininger committed
853
    def __abs__(self):
854
        return Field(self.domain,np.abs(self.val),self.dtype)
csongor's avatar
csongor committed
855

856
    def _contraction_helper(self, op, spaces):
Theo Steininger's avatar
Theo Steininger committed
857
        if spaces is None:
858
859
860
            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
861

862
        axes_list = tuple(self.domain_axes[sp_index] for sp_index in spaces)
863
864

        try:
Theo Steininger's avatar
Theo Steininger committed
865
            axes_list = reduce(lambda x, y: x+y, axes_list)
866
        except TypeError:
Theo Steininger's avatar
Theo Steininger committed
867
            axes_list = ()
csongor's avatar
csongor committed
868

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
869
        # perform the contraction on the data
Theo Steininger's avatar
Theo Steininger committed
870
871
        data = self.get_val(copy=False)
        data = getattr(data, op)(axis=axes_list)
csongor's avatar
csongor committed
872

Theo Steininger's avatar
Theo Steininger committed
873
874
875
        # 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
876
        else:
Theo Steininger's avatar
Theo Steininger committed
877
            return_domain = tuple(self.domain[i]
Martin Reinecke's avatar
Martin Reinecke committed
878
                                  for i in range(len(self.domain))
Theo Steininger's avatar
Theo Steininger committed
879
                                  if i not in spaces)
880

Theo Steininger's avatar
Theo Steininger committed
881
882
883
884
            return_field = Field(domain=return_domain,
                                 val=data,
                                 copy=False)
            return return_field
csongor's avatar
csongor committed
885

886
887
    def sum(self, spaces=None):
        return self._contraction_helper('sum', spaces)
csongor's avatar
csongor committed
888

889
890
    def prod(self, spaces=None):
        return self._contraction_helper('prod', spaces)
csongor's avatar
csongor committed
891

892
893
    def all(self, spaces=None):
        return self._contraction_helper('all', spaces)
csongor's avatar
csongor committed
894

895
896
    def any(self, spaces=None):
        return self._contraction_helper('any', spaces)
csongor's avatar
csongor committed
897

898
899
    def min(self, spaces=None):
        return self._contraction_helper('min', spaces)
csongor's avatar
csongor committed
900

901
902
    def nanmin(self, spaces=None):
        return self._contraction_helper('nanmin', spaces)
csongor's avatar
csongor committed
903

904
905
    def max(self, spaces=None):
        return self._contraction_helper('max', spaces)
csongor's avatar
csongor committed
906

907
908
    def nanmax(self, spaces=None):
        return self._contraction_helper('nanmax', spaces)
csongor's avatar
csongor committed
909

910
911
    def mean(self, spaces=None):
        return self._contraction_helper('mean', spaces)
csongor's avatar
csongor committed
912

913
914
    def var(self, spaces=None):
        return self._contraction_helper('var', spaces)
csongor's avatar
csongor committed
915

916
917
    def std(self, spaces=None):
        return self._contraction_helper('std', spaces)
csongor's avatar
csongor committed
918

Theo Steininger's avatar
Theo Steininger committed
919
    # ---General binary methods---
csongor's avatar
csongor committed
920

921
    def _binary_helper(self, other, op):