field.py 30.3 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):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
125
126
        if val is None:
            return np.float64 if dtype is None else dtype
csongor's avatar
csongor committed
127
        if dtype is None:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
128
129
130
            if isinstance(val,Field):
                return val.dtype
            return np.result_type(val)
131

132
    # ---Factory methods---
133

134
    @classmethod
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
135
    def from_random(cls, random_type, domain=None, dtype=None, **kwargs):
136
137
138
139
140
        """ Draws a random field with the given parameters.

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

142
143
144
        random_type : String
            'pm1', 'normal', 'uniform' are the supported arguments for this
            method.
Theo Steininger's avatar
Theo Steininger committed
145

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

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

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

        See Also
        --------
159
        power_synthesize
Theo Steininger's avatar
Theo Steininger committed
160

161
162

        """
Theo Steininger's avatar
Theo Steininger committed
163

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
164
        f = cls(domain=domain, dtype=dtype)
165
        generator_function = getattr(Random, random_type)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
166
        f.val=generator_function(dtype=f.dtype, shape=f.shape, **kwargs)
167
168
        return f

169
170
    # ---Powerspectral methods---

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

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

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

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

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

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

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

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

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

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

236
237
238
239
240
241
242
243
244
245
        if keep_phase_information:
            parts_val = self._hermitian_decomposition(
                                              val=self.val,
                                              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]
246
247

        for space_index in spaces:
248
249
            parts = [self._single_power_analyze(
                                work_field=part,
250
                                space_index=space_index,
251
252
                                binbounds=binbounds)
                     for part in parts]
253

254
255
256
257
258
259
        if keep_phase_information:
            result_field = parts[0] + 1j*parts[1]
        else:
            result_field = parts[0]

        return result_field
260
261

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

264
        if not work_field.domain[space_index].harmonic:
265
266
            raise ValueError(
                "The analyzed space must be harmonic.")
267

268
269
270
271
272
273
        # 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.

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

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

287
        result_field = work_field.copy_empty(
288
                   domain=result_domain,
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
289
                   dtype=result_dtype)
290
291
292
293
        result_field.set_val(new_val=power_spectrum, copy=False)

        return result_field

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

Martin Reinecke's avatar
Martin Reinecke committed
297
298
299
        pindex = pdomain.pindex
        # MR FIXME: how about iterating over slices, instead of replicating
        # pindex? Would save memory and probably isn't slower.
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,
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
cleanup    
Martin Reinecke committed
414
        result_list[0].val *= spec.real
415
416

        if not real_power:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
417
            result_list[1].val *= spec.imag
418

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

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

        return result

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

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

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

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

Theo Steininger's avatar
Theo Steininger committed
457
    def set_val(self, new_val=None, copy=False):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
458
        """ Sets the field's data object.
459
460
461

        Parameters
        ----------
Theo Steininger's avatar
Theo Steininger committed
462
        new_val : scalar, array-like, Field, None *optional*
463
464
            The values to be stored in the field.
            {default : None}
Theo Steininger's avatar
Theo Steininger committed
465

466
        copy : boolean, *optional*
Theo Steininger's avatar
Theo Steininger committed
467
468
            If False, Field tries to not copy the input data but use it
            directly.
469
470
471
472
473
474
            {default : False}
        See Also
        --------
        val

        """
Theo Steininger's avatar
Theo Steininger committed
475

Martin Reinecke's avatar
Martin Reinecke committed
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
        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")
496
        return self
csongor's avatar
csongor committed
497

498
    def get_val(self, copy=False):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
499
        """ Returns the data object associated with this Field.
500
501
502
503

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

507
508
        Returns
        -------
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
509
        out : numpy.ndarray
510
511
512
513
514
515

        See Also
        --------
        val

        """
Theo Steininger's avatar
Theo Steininger committed
516

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

Theo Steininger's avatar
Theo Steininger committed
519
520
    @property
    def val(self):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
521
        """ Returns the data object associated with this Field.
Theo Steininger's avatar
Theo Steininger committed
522

523
524
        Returns
        -------
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
525
        out : numpy.ndarray
526
527
528
529
530
531

        See Also
        --------
        get_val

        """
Theo Steininger's avatar
Theo Steininger committed
532

Martin Reinecke's avatar
Martin Reinecke committed
533
        return self._val
csongor's avatar
csongor committed
534

Theo Steininger's avatar
Theo Steininger committed
535
536
    @val.setter
    def val(self, new_val):
537
        self.set_val(new_val=new_val, copy=False)
csongor's avatar
csongor committed
538

Martin Reinecke's avatar
Martin Reinecke committed
539
540
541
542
    @property
    def dtype(self):
        return self._val.dtype

543
544
    @property
    def shape(self):
Theo Steininger's avatar
Theo Steininger committed
545
        """ Returns the total shape of the Field's data array.
Theo Steininger's avatar
Theo Steininger committed
546

547
548
549
        Returns
        -------
        out : tuple
Martin Reinecke's avatar
Martin Reinecke committed
550
            The output object. The tuple contains the dimensions of the spaces
551
552
553
554
555
556
557
            in domain.

        See Also
        --------
        dim

        """
Martin Reinecke's avatar
Martin Reinecke committed
558
        return self._val.shape
csongor's avatar
csongor committed
559

560
561
    @property
    def dim(self):
Theo Steininger's avatar
Theo Steininger committed
562
        """ Returns the total number of pixel-dimensions the field has.
Theo Steininger's avatar
Theo Steininger committed
563

Theo Steininger's avatar
Theo Steininger committed
564
        Effectively, all values from shape are multiplied.
Theo Steininger's avatar
Theo Steininger committed
565

566
567
568
569
570
571
572
573
574
575
        Returns
        -------
        out : int
            The dimension of the Field.

        See Also
        --------
        shape

        """
Theo Steininger's avatar
Theo Steininger committed
576

577
        dim_tuple = tuple(sp.dim for sp in self.domain)
Theo Steininger's avatar
Theo Steininger committed
578
        try:
Martin Reinecke's avatar
Martin Reinecke committed
579
            return int(reduce(lambda x, y: x * y, dim_tuple))
Theo Steininger's avatar
Theo Steininger committed
580
581
        except TypeError:
            return 0
csongor's avatar
csongor committed
582

583
584
    @property
    def dof(self):
Theo Steininger's avatar
Theo Steininger committed
585
586
587
588
589
590
        """ 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
591
592
593
594
595
596
597
        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
598
599
600
        """ Returns the total volume of all spaces in the domain.
        """

Theo Steininger's avatar
Theo Steininger committed
601
        volume_tuple = tuple(sp.total_volume for sp in self.domain)
602
        try:
Theo Steininger's avatar
Theo Steininger committed
603
            return reduce(lambda x, y: x * y, volume_tuple)
604
        except TypeError:
Theo Steininger's avatar
Theo Steininger committed
605
            return 0.
606

Theo Steininger's avatar
Theo Steininger committed
607
608
609
610
    @property
    def real(self):
        """ The real part of the field (data is not copied).
        """
611
        return Field(self.domain,self.val.real)
Theo Steininger's avatar
Theo Steininger committed
612
613
614
615
616

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

Theo Steininger's avatar
Theo Steininger committed
619
    # ---Special unary/binary operations---
620

csongor's avatar
csongor committed
621

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

625
626
        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
627
        able to define the domain and the dtype of the returned Field.
628
629
630
631
632

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

634
635
        dtype : type
            The new dtype the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
636

637
638
639
640
641
642
643
644
645
646
        Returns
        -------
        out : Field
            The output object. An identical copy of 'self'.

        See Also
        --------
        copy_empty

        """
Theo Steininger's avatar
Theo Steininger committed
647

648
649
650
        if domain is None:
            domain = self.domain
        return Field(domain=domain,val=self._val,dtype=dtype,copy=True)
csongor's avatar
csongor committed
651

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
652
    def copy_empty(self, domain=None, dtype=None):
653
654
655
        """ 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
656
657
658
        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
659
        to change the domain and the dtype of the returned Field.
Theo Steininger's avatar
Theo Steininger committed
660

661
662
663
664
        Parameters
        ----------
        domain : DomainObject
            The new domain the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
665

666
667
        dtype : type
            The new dtype the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
668

669
670
671
        Returns
        -------
        out : Field
Theo Steininger's avatar
Theo Steininger committed
672
            The output object.
673
674
675
676
677
678

        See Also
        --------
        copy

        """
Theo Steininger's avatar
Theo Steininger committed
679

Theo Steininger's avatar
Theo Steininger committed
680
681
682
683
        if domain is None:
            domain = self.domain
        if dtype is None:
            dtype = self.dtype
684
        return Field(domain=domain, dtype=dtype)
Theo Steininger's avatar
Theo Steininger committed
685
686

    def weight(self, power=1, inplace=False, spaces=None):
Theo Steininger's avatar
Theo Steininger committed
687
        """ Weights the pixels of `self` with their invidual pixel-volume.
688
689
690
691

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

694
        inplace : boolean
Theo Steininger's avatar
Theo Steininger committed
695
696
            If True, `self` will be weighted and returned. Otherwise, a copy
            is made.
Theo Steininger's avatar
Theo Steininger committed
697

Theo Steininger's avatar
Theo Steininger committed
698
699
        spaces : tuple of ints
            Determines on which subspace the operation takes place.
Theo Steininger's avatar
Theo Steininger committed
700

701
702
703
        Returns
        -------
        out : Field
Theo Steininger's avatar
Theo Steininger committed
704
            The weighted field.
705
706

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

709
        new_val = self.get_val(copy=False)
csongor's avatar
csongor committed
710

711
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
csongor's avatar
csongor committed
712
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
713
            spaces = list(range(len(self.domain)))
csongor's avatar
csongor committed
714

715
        for ind, sp in enumerate(self.domain):
Theo Steininger's avatar
Theo Steininger committed
716
717
718
719
720
            if ind in spaces:
                new_val = sp.weight(new_val,
                                    power=power,
                                    axes=self.domain_axes[ind],
                                    inplace=inplace)
721
722

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

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

728
729
730
        Parameters
        ----------
        x : Field
Theo Steininger's avatar
Theo Steininger committed
731
            The domain of x must contain `self.domain`
Theo Steininger's avatar
Theo Steininger committed
732

Theo Steininger's avatar
Theo Steininger committed
733
734
735
        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
736

737
        bare : boolean
Theo Steininger's avatar
Theo Steininger committed
738
            If true, no volume factors will be included in the computation.
Theo Steininger's avatar
Theo Steininger committed
739

740
741
742
        Returns
        -------
        out : float, complex
Theo Steininger's avatar
Theo Steininger committed
743

744
        """
745
746
747
        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
748

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

752
        if spaces is None:
753
            return np.vdot(y.val.flatten(),x.val.flatten())
754
755
756
        else:
            # create a diagonal operator which is capable of taking care of the
            # axes-matching
757
            from .operators.diagonal_operator import DiagonalOperator
758
759
760
761
762
763
            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
764

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

Theo Steininger's avatar
Theo Steininger committed
768
769
770
        Returns
        -------
        norm : scalar
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
771
            The L2-norm of the field values.
csongor's avatar
csongor committed
772
773

        """
Martin Reinecke's avatar
fixes    
Martin Reinecke committed
774
        return np.sqrt(np.abs(self.vdot(x=self)))
csongor's avatar
csongor committed
775
776

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

779
780
781
        Parameters
        ----------
        inplace : boolean
Theo Steininger's avatar
Theo Steininger committed
782
            Decides whether the conjugation should be performed inplace.
Theo Steininger's avatar
Theo Steininger committed
783

784
785
786
787
        Returns
        -------
        cc : field
            The complex conjugated field.
csongor's avatar
csongor committed
788
789

        """
Theo Steininger's avatar
Theo Steininger committed
790

csongor's avatar
csongor committed
791
        if inplace:
792
793
            self.imag*=-1
            return self
csongor's avatar
csongor committed
794
        else:
795
            return Field(self.domain,np.conj(self.val),self.dtype)
csongor's avatar
csongor committed
796

Theo Steininger's avatar
Theo Steininger committed
797
    # ---General unary/contraction methods---
798

Theo Steininger's avatar
Theo Steininger committed
799
800
    def __pos__(self):
        return self.copy()
801

Theo Steininger's avatar
Theo Steininger committed
802
    def __neg__(self):
803
        return Field(self.domain,-self.val,self.dtype)
csongor's avatar
csongor committed
804

Theo Steininger's avatar
Theo Steininger committed
805
    def __abs__(self):
806
        return Field(self.domain,np.abs(self.val),self.dtype)
csongor's avatar
csongor committed
807

808
    def _contraction_helper(self, op, spaces):
Theo Steininger's avatar
Theo Steininger committed
809
        if spaces is None:
810
811
812
            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
813

814
        axes_list = tuple(self.domain_axes[sp_index] for sp_index in spaces)
815
816

        try:
Theo Steininger's avatar
Theo Steininger committed
817
            axes_list = reduce(lambda x, y: x+y, axes_list)
818
        except TypeError:
Theo Steininger's avatar
Theo Steininger committed
819
            axes_list = ()
csongor's avatar
csongor committed
820

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
821
        # perform the contraction on the data
Theo Steininger's avatar
Theo Steininger committed
822
823
        data = self.get_val(copy=False)
        data = getattr(data, op)(axis=axes_list)
csongor's avatar
csongor committed
824

Theo Steininger's avatar
Theo Steininger committed
825
826
827
        # 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
828
        else:
Theo Steininger's avatar
Theo Steininger committed
829
            return_domain = tuple(self.domain[i]
Martin Reinecke's avatar
Martin Reinecke committed
830
                                  for i in range(len(self.domain))
Theo Steininger's avatar
Theo Steininger committed
831
                                  if i not in spaces)
832

Theo Steininger's avatar
Theo Steininger committed
833
834
835
836
            return_field = Field(domain=return_domain,
                                 val=data,
                                 copy=False)
            return return_field
csongor's avatar
csongor committed
837

838
839
    def sum(self, spaces=None):
        return self._contraction_helper('sum', spaces)
csongor's avatar
csongor committed
840

841
842
    def prod(self, spaces=None):
        return self._contraction_helper('prod', spaces)
csongor's avatar
csongor committed
843

844
845
    def all(self, spaces=None):
        return self._contraction_helper('all', spaces)
csongor's avatar
csongor committed
846

847
848
    def any(self, spaces=None):
        return self._contraction_helper('any', spaces)
csongor's avatar
csongor committed
849

850
851
    def min(self, spaces=None):
        return self._contraction_helper('min', spaces)
csongor's avatar
csongor committed
852

853
854
    def nanmin(self, spaces=None):
        return self._contraction_helper('nanmin', spaces)
csongor's avatar
csongor committed
855

856
857
    def max(self, spaces=None):
        return self._contraction_helper('max', spaces)
csongor's avatar
csongor committed
858

859
860
    def nanmax(self, spaces=None):
        return self._contraction_helper('nanmax', spaces)
csongor's avatar
csongor committed
861

862
863
    def mean(self, spaces=None):
        return self._contraction_helper('mean', spaces)
csongor's avatar
csongor committed
864

865
866
    def var(self, spaces=None):
        return self._contraction_helper('var', spaces)
csongor's avatar
csongor committed
867

868
869
    def std(self, spaces=None):
        return self._contraction_helper('std', spaces)
csongor's avatar
csongor committed
870

Theo Steininger's avatar
Theo Steininger committed
871
    # ---General binary methods---
csongor's avatar
csongor committed
872

873
    def _binary_helper(self, other, op):
csongor's avatar
csongor committed
874
        # if other is a field, make sure that the domains match
875
        if isinstance(other, Field):
876
877
878
            if other.domain != self.domain:
                raise ValueError("domains are incompatible.")
            return Field(self.domain,getattr(self.val,op)(other.val))
csongor's avatar
csongor committed
879

880
        return Field(self.domain,getattr(self.val,op)(other))
csongor's avatar
csongor committed
881
882

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

885
    def __radd__(self, other):
Theo Steininger's avatar
Theo Steininger committed
886
        return self._binary_helper(other, op='__radd__')
csongor's avatar
csongor committed
887
888

    def __iadd__(self, other):
889
        return self._binary_helper(other, op='__iadd__')
csongor's avatar
csongor committed
890
891

    def __sub__(self, other):
Theo Steininger's avatar
Theo Steininger committed
892
        return self._binary_helper(other, op='__sub__')
csongor's avatar
csongor committed
893
894

    def __rsub__(self, other):
Theo Steininger's avatar
Theo Steininger committed
895
        return self._binary_helper(other, op='__rsub__')
csongor's avatar
csongor committed
896
897

    def __isub__(self, other):
898
        return self._binary_helper(other, op='__isub__')
csongor's avatar
csongor committed
899
900

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

903
    def __rmul__(self, other):
Theo Steininger's avatar
Theo Steininger committed
904
        return self._binary_helper(other, op='__rmul__')
csongor's avatar
csongor committed
905
906

    def __imul__(self, other):