field.py 29.1 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
cleanup    
Martin Reinecke committed
135
136
        if val is None:
            return np.float64 if dtype is None else dtype
csongor's avatar
csongor committed
137
        if dtype is None:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
138
            if isinstance(val, Field):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
139
140
                return val.dtype
            return np.result_type(val)
141

142
    # ---Factory methods---
143

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

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

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

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

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

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

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

171
172

        """
Theo Steininger's avatar
Theo Steininger committed
173

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
174
        f = cls(domain=domain, dtype=dtype)
175
        generator_function = getattr(Random, random_type)
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
176
        f.val = generator_function(dtype=f.dtype, shape=f.shape, **kwargs)
177
178
        return f

179
180
    # ---Powerspectral methods---

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

383
384
385
386
387
        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.")
388
389
390

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

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

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

        # apply the rescaler to the random fields
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
415
        result_list[0].val *= spec.real
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
PEP8    
Martin Reinecke committed
420
            result_list = [Field(i.domain, self._hermitian_decomposition(
Martin Reinecke's avatar
Martin Reinecke committed
421
422
                                     i.val,
                                     preserve_gaussian_variance=True)[0])
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
423
                           for i in result_list]
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
        if new_val is None:
            pass
        elif isinstance(new_val, Field):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
479
            if self.domain != new_val.domain:
Martin Reinecke's avatar
Martin Reinecke committed
480
481
482
483
                raise ValueError("Domain mismatch")
            if copy:
                self._val[()] = new_val.val
            else:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
484
                self._val = np.array(new_val.val, dtype=self.dtype, copy=False)
Martin Reinecke's avatar
Martin Reinecke committed
485
        elif (np.isscalar(new_val)):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
486
            self._val[()] = new_val
Martin Reinecke's avatar
Martin Reinecke committed
487
488
489
490
        elif isinstance(new_val, np.ndarray):
            if copy:
                self._val[()] = new_val
            else:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
491
                if self.shape != new_val.shape:
Martin Reinecke's avatar
Martin Reinecke committed
492
                    raise ValueError("Shape mismatch")
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
493
                self._val = np.array(new_val, dtype=self.dtype, copy=False)
Martin Reinecke's avatar
Martin Reinecke committed
494
495
        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

        See Also
        --------
        val
        """
Martin Reinecke's avatar
Martin Reinecke committed
515
        return self._val.copy() if copy else self._val
csongor's avatar
csongor committed
516

Theo Steininger's avatar
Theo Steininger committed
517
518
    @property
    def val(self):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
519
        """ Returns the data object associated with this Field.
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
520
        No copy is made.
Theo Steininger's avatar
Theo Steininger committed
521

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

        See Also
        --------
        get_val
        """
Martin Reinecke's avatar
Martin Reinecke committed
530
        return self._val
csongor's avatar
csongor committed
531

Theo Steininger's avatar
Theo Steininger committed
532
533
    @val.setter
    def val(self, new_val):
534
        self.set_val(new_val=new_val, copy=False)
csongor's avatar
csongor committed
535

Martin Reinecke's avatar
Martin Reinecke committed
536
537
538
539
    @property
    def dtype(self):
        return self._val.dtype

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

544
545
546
        Returns
        -------
        out : tuple
Martin Reinecke's avatar
Martin Reinecke committed
547
            The output object. The tuple contains the dimensions of the spaces
548
            in domain.
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
549
       """
Martin Reinecke's avatar
Martin Reinecke committed
550
        return self._val.shape
csongor's avatar
csongor committed
551

552
553
    @property
    def dim(self):
Theo Steininger's avatar
Theo Steininger committed
554
        """ Returns the total number of pixel-dimensions the field has.
Theo Steininger's avatar
Theo Steininger committed
555

Theo Steininger's avatar
Theo Steininger committed
556
        Effectively, all values from shape are multiplied.
Theo Steininger's avatar
Theo Steininger committed
557

558
559
560
561
562
        Returns
        -------
        out : int
            The dimension of the Field.
        """
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
563
        return self._val.size
csongor's avatar
csongor committed
564

Theo Steininger's avatar
Theo Steininger committed
565
566
    @property
    def total_volume(self):
Theo Steininger's avatar
Theo Steininger committed
567
568
        """ Returns the total volume of all spaces in the domain.
        """
Theo Steininger's avatar
Theo Steininger committed
569
        volume_tuple = tuple(sp.total_volume for sp in self.domain)
570
        try:
Theo Steininger's avatar
Theo Steininger committed
571
            return reduce(lambda x, y: x * y, volume_tuple)
572
        except TypeError:
Theo Steininger's avatar
Theo Steininger committed
573
            return 0.
574

Theo Steininger's avatar
Theo Steininger committed
575
576
577
578
    @property
    def real(self):
        """ The real part of the field (data is not copied).
        """
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
579
        return Field(self.domain, self.val.real)
Theo Steininger's avatar
Theo Steininger committed
580
581
582
583
584

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

Theo Steininger's avatar
Theo Steininger committed
587
    # ---Special unary/binary operations---
588

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

592
593
        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
594
        able to define the domain and the dtype of the returned Field.
595
596
597
598
599

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

601
602
        dtype : type
            The new dtype the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
603

604
605
606
607
608
609
        Returns
        -------
        out : Field
            The output object. An identical copy of 'self'.

        """
Theo Steininger's avatar
Theo Steininger committed
610

611
612
        if domain is None:
            domain = self.domain
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
613
        return Field(domain=domain, val=self._val, dtype=dtype, copy=True)
csongor's avatar
csongor committed
614

Theo Steininger's avatar
Theo Steininger committed
615
    def weight(self, power=1, inplace=False, spaces=None):
Theo Steininger's avatar
Theo Steininger committed
616
        """ Weights the pixels of `self` with their invidual pixel-volume.
617
618
619
620

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

623
        inplace : boolean
Theo Steininger's avatar
Theo Steininger committed
624
625
            If True, `self` will be weighted and returned. Otherwise, a copy
            is made.
Theo Steininger's avatar
Theo Steininger committed
626

Theo Steininger's avatar
Theo Steininger committed
627
628
        spaces : tuple of ints
            Determines on which subspace the operation takes place.
Theo Steininger's avatar
Theo Steininger committed
629

630
631
632
        Returns
        -------
        out : Field
Theo Steininger's avatar
Theo Steininger committed
633
            The weighted field.
634
635

        """
636
        new_val = self.get_val(copy=not inplace)
csongor's avatar
csongor committed
637

638
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
csongor's avatar
csongor committed
639
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
640
            spaces = list(range(len(self.domain)))
csongor's avatar
csongor committed
641

642
643
644
645
646
647
648
649
650
651
652
        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)
                new_val *= wgt**power
        fct = fct**power
Martin Reinecke's avatar
Martin Reinecke committed
653
        if fct != 1.:
654
            new_val *= fct
655

Martin Reinecke's avatar
updates    
Martin Reinecke committed
656
        return Field(self.domain, new_val, self.dtype)
csongor's avatar
csongor committed
657

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

661
662
663
        Parameters
        ----------
        x : Field
Theo Steininger's avatar
Theo Steininger committed
664
            The domain of x must contain `self.domain`
Theo Steininger's avatar
Theo Steininger committed
665

Theo Steininger's avatar
Theo Steininger committed
666
667
668
        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
669

670
        bare : boolean
Theo Steininger's avatar
Theo Steininger committed
671
            If true, no volume factors will be included in the computation.
Theo Steininger's avatar
Theo Steininger committed
672

673
674
675
        Returns
        -------
        out : float, complex
Theo Steininger's avatar
Theo Steininger committed
676

677
        """
678
679
680
        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
681

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

685
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
686
            return np.vdot(y.val.reshape(-1), x.val.reshape(-1))
687
688
689
        else:
            # create a diagonal operator which is capable of taking care of the
            # axes-matching
690
            from .operators.diagonal_operator import DiagonalOperator
691
692
693
694
695
696
            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
697

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

Theo Steininger's avatar
Theo Steininger committed
701
702
703
        Returns
        -------
        norm : scalar
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
704
            The L2-norm of the field values.
csongor's avatar
csongor committed
705
706

        """
Martin Reinecke's avatar
fixes    
Martin Reinecke committed
707
        return np.sqrt(np.abs(self.vdot(x=self)))
csongor's avatar
csongor committed
708
709

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

712
713
714
        Parameters
        ----------
        inplace : boolean
Theo Steininger's avatar
Theo Steininger committed
715
            Decides whether the conjugation should be performed inplace.
Theo Steininger's avatar
Theo Steininger committed
716

717
718
719
720
        Returns
        -------
        cc : field
            The complex conjugated field.
csongor's avatar
csongor committed
721
722
723

        """
        if inplace:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
724
            self.imag *= -1
725
            return self
csongor's avatar
csongor committed
726
        else:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
727
            return Field(self.domain, np.conj(self.val), self.dtype)
csongor's avatar
csongor committed
728

Theo Steininger's avatar
Theo Steininger committed
729
    # ---General unary/contraction methods---
730

Theo Steininger's avatar
Theo Steininger committed
731
732
    def __pos__(self):
        return self.copy()
733

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

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

740
    def _contraction_helper(self, op, spaces):
Theo Steininger's avatar
Theo Steininger committed
741
        if spaces is None:
742
743
744
            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
745

746
        axes_list = tuple(self.domain_axes[sp_index] for sp_index in spaces)
747
748

        try:
Theo Steininger's avatar
Theo Steininger committed
749
            axes_list = reduce(lambda x, y: x+y, axes_list)
750
        except TypeError:
Theo Steininger's avatar
Theo Steininger committed
751
            axes_list = ()
csongor's avatar
csongor committed
752

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
753
        # perform the contraction on the data
Theo Steininger's avatar
Theo Steininger committed
754
755
        data = self.get_val(copy=False)
        data = getattr(data, op)(axis=axes_list)
csongor's avatar
csongor committed
756

Theo Steininger's avatar
Theo Steininger committed
757
758
759
        # 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
760
        else:
Theo Steininger's avatar
Theo Steininger committed
761
            return_domain = tuple(self.domain[i]
Martin Reinecke's avatar
Martin Reinecke committed
762
                                  for i in range(len(self.domain))
Theo Steininger's avatar
Theo Steininger committed
763
                                  if i not in spaces)
764

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

767
768
    def sum(self, spaces=None):
        return self._contraction_helper('sum', spaces)
csongor's avatar
csongor committed
769

770
771
    def prod(self, spaces=None):
        return self._contraction_helper('prod', spaces)
csongor's avatar
csongor committed
772

773
774
    def all(self, spaces=None):
        return self._contraction_helper('all', spaces)
csongor's avatar
csongor committed
775

776
777
    def any(self, spaces=None):
        return self._contraction_helper('any', spaces)
csongor's avatar
csongor committed
778

779
780
    def min(self, spaces=None):
        return self._contraction_helper('min', spaces)
csongor's avatar
csongor committed
781

782
783
    def nanmin(self, spaces=None):
        return self._contraction_helper('nanmin', spaces)
csongor's avatar
csongor committed
784

785
786
    def max(self, spaces=None):
        return self._contraction_helper('max', spaces)
csongor's avatar
csongor committed
787

788
789
    def nanmax(self, spaces=None):
        return self._contraction_helper('nanmax', spaces)
csongor's avatar
csongor committed
790

791
792
    def mean(self, spaces=None):
        return self._contraction_helper('mean', spaces)
csongor's avatar
csongor committed
793

794
795
    def var(self, spaces=None):
        return self._contraction_helper('var', spaces)
csongor's avatar
csongor committed
796

797
798
    def std(self, spaces=None):
        return self._contraction_helper('std', spaces)
csongor's avatar
csongor committed
799

Theo Steininger's avatar
Theo Steininger committed
800
    # ---General binary methods---
csongor's avatar
csongor committed
801

802
    def _binary_helper(self, other, op):
csongor's avatar
csongor committed
803
        # if other is a field, make sure that the domains match
804
        if isinstance(other, Field):
805
806
            if other.domain != self.domain:
                raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
807
            return Field(self.domain, getattr(self.val, op)(other.val))
csongor's avatar
csongor committed
808

Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
809
        return Field(self.domain, getattr(self.val, op)(other))
csongor's avatar
csongor committed
810
811

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

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

    def __iadd__(self, other):
818
        return self._binary_helper(other, op='__iadd__')
csongor's avatar
csongor committed
819
820

    def __sub__(self, other):
Theo Steininger's avatar
Theo Steininger committed
821
        return self._binary_helper(other, op='__sub__')
csongor's avatar
csongor committed
822
823

    def __rsub__(self, other):
Theo Steininger's avatar
Theo Steininger committed
824
        return self._binary_helper(other, op='__rsub__')
csongor's avatar
csongor committed
825
826

    def __isub__(self, other):
827
        return self._binary_helper(other, op='__isub__')
csongor's avatar
csongor committed
<