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
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
109
            if isinstance(val, Field):
                domain = val.domain
Martin Reinecke's avatar
Martin Reinecke committed
110
            elif np.isscalar(val):
111
                domain = ()
Martin Reinecke's avatar
Martin Reinecke committed
112
113
            else:
                raise TypeError("could not infer domain from value")
114
        elif isinstance(domain, DomainObject):
115
            domain = (domain,)
116
117
118
        elif not isinstance(domain, tuple):
            domain = tuple(domain)

csongor's avatar
csongor committed
119
        for d in domain:
120
            if not isinstance(d, DomainObject):
121
122
                raise TypeError(
                    "Given domain contains something that is not a "
123
                    "DomainObject instance.")
csongor's avatar
csongor committed
124
125
        return domain

126
127
    def _get_axes_tuple(self, things_with_shape):
        i = 0
Theo Steininger's avatar
Theo Steininger committed
128
129
        axes_list = []
        for thing in things_with_shape:
130
            nax = len(thing.shape)
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
131
            axes_list += [tuple(range(i, i+nax))]
132
            i += nax
Theo Steininger's avatar
Theo Steininger committed
133
        return tuple(axes_list)
134

135
    def _infer_dtype(self, dtype, val):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
136
137
        if val is None:
            return np.float64 if dtype is None else dtype
csongor's avatar
csongor committed
138
        if dtype is None:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
139
            if isinstance(val, Field):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
140
141
                return val.dtype
            return np.result_type(val)
142

143
    # ---Factory methods---
144

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

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

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

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

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

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

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

172
173

        """
Theo Steininger's avatar
Theo Steininger committed
174

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

180
181
    # ---Powerspectral methods---

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

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

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

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

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

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

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

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

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

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

247
248
249
250
251
252
253
254
255
256
        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]
257
258

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

        return result

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

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

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

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

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

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

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

        """
Theo Steininger's avatar
Theo Steininger committed
476

Martin Reinecke's avatar
Martin Reinecke committed
477
478
479
        if new_val is None:
            pass
        elif isinstance(new_val, Field):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
480
            if self.domain != new_val.domain:
Martin Reinecke's avatar
Martin Reinecke committed
481
482
483
484
                raise ValueError("Domain mismatch")
            if copy:
                self._val[()] = new_val.val
            else:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
485
                self._val = np.array(new_val.val, dtype=self.dtype, copy=False)
Martin Reinecke's avatar
Martin Reinecke committed
486
        elif (np.isscalar(new_val)):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
487
            self._val[()] = new_val
Martin Reinecke's avatar
Martin Reinecke committed
488
489
490
491
        elif isinstance(new_val, np.ndarray):
            if copy:
                self._val[()] = new_val
            else:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
492
                if self.shape != new_val.shape:
Martin Reinecke's avatar
Martin Reinecke committed
493
                    raise ValueError("Shape mismatch")
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
494
                self._val = np.array(new_val, dtype=self.dtype, copy=False)
Martin Reinecke's avatar
Martin Reinecke committed
495
496
        else:
            raise TypeError("unknown source type")
497
        return self
csongor's avatar
csongor committed
498

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

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

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

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

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

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

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

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

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

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

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

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

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

566
567
    @property
    def dof(self):
Theo Steininger's avatar
Theo Steininger committed
568
569
570
571
572
573
        """ 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
574
575
576
577
578
579
580
        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
581
582
        """ Returns the total volume of all spaces in the domain.
        """
Theo Steininger's avatar
Theo Steininger committed
583
        volume_tuple = tuple(sp.total_volume for sp in self.domain)
584
        try:
Theo Steininger's avatar
Theo Steininger committed
585
            return reduce(lambda x, y: x * y, volume_tuple)
586
        except TypeError:
Theo Steininger's avatar
Theo Steininger committed
587
            return 0.
588

Theo Steininger's avatar
Theo Steininger committed
589
590
591
592
    @property
    def real(self):
        """ The real part of the field (data is not copied).
        """
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
593
        return Field(self.domain, self.val.real)
Theo Steininger's avatar
Theo Steininger committed
594
595
596
597
598

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

Theo Steininger's avatar
Theo Steininger committed
601
    # ---Special unary/binary operations---
602

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

606
607
        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
608
        able to define the domain and the dtype of the returned Field.
609
610
611
612
613

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

615
616
        dtype : type
            The new dtype the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
617

618
619
620
621
622
623
624
625
626
627
        Returns
        -------
        out : Field
            The output object. An identical copy of 'self'.

        See Also
        --------
        copy_empty

        """
Theo Steininger's avatar
Theo Steininger committed
628

629
630
        if domain is None:
            domain = self.domain
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
631
        return Field(domain=domain, val=self._val, dtype=dtype, copy=True)
csongor's avatar
csongor committed
632

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
633
    def copy_empty(self, domain=None, dtype=None):
634
635
636
        """ 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
637
638
639
        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
640
        to change the domain and the dtype of the returned Field.
Theo Steininger's avatar
Theo Steininger committed
641

642
643
644
645
        Parameters
        ----------
        domain : DomainObject
            The new domain the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
646

647
648
        dtype : type
            The new dtype the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
649

650
651
652
        Returns
        -------
        out : Field
Theo Steininger's avatar
Theo Steininger committed
653
            The output object.
654
655
656
657
658
659

        See Also
        --------
        copy

        """
Theo Steininger's avatar
Theo Steininger committed
660

Theo Steininger's avatar
Theo Steininger committed
661
662
663
664
        if domain is None:
            domain = self.domain
        if dtype is None:
            dtype = self.dtype
665
        return Field(domain=domain, dtype=dtype)
Theo Steininger's avatar
Theo Steininger committed
666
667

    def weight(self, power=1, inplace=False, spaces=None):
Theo Steininger's avatar
Theo Steininger committed
668
        """ Weights the pixels of `self` with their invidual pixel-volume.
669
670
671
672

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

675
        inplace : boolean
Theo Steininger's avatar
Theo Steininger committed
676
677
            If True, `self` will be weighted and returned. Otherwise, a copy
            is made.
Theo Steininger's avatar
Theo Steininger committed
678

Theo Steininger's avatar
Theo Steininger committed
679
680
        spaces : tuple of ints
            Determines on which subspace the operation takes place.
Theo Steininger's avatar
Theo Steininger committed
681

682
683
684
        Returns
        -------
        out : Field
Theo Steininger's avatar
Theo Steininger committed
685
            The weighted field.
686
687

        """
688
        new_val = self.get_val(copy=False)
csongor's avatar
csongor committed
689

690
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
csongor's avatar
csongor committed
691
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
692
            spaces = list(range(len(self.domain)))
csongor's avatar
csongor committed
693

694
        for ind, sp in enumerate(self.domain):
Theo Steininger's avatar
Theo Steininger committed
695
696
697
698
699
            if ind in spaces:
                new_val = sp.weight(new_val,
                                    power=power,
                                    axes=self.domain_axes[ind],
                                    inplace=inplace)
Martin Reinecke's avatar
updates    
Martin Reinecke committed
700
701
                # we need at most one copy, the rest can happen in place
                inplace = True
702

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

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

708
709
710
        Parameters
        ----------
        x : Field
Theo Steininger's avatar
Theo Steininger committed
711
            The domain of x must contain `self.domain`
Theo Steininger's avatar
Theo Steininger committed
712

Theo Steininger's avatar
Theo Steininger committed
713
714
715
        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
716

717
        bare : boolean
Theo Steininger's avatar
Theo Steininger committed
718
            If true, no volume factors will be included in the computation.
Theo Steininger's avatar
Theo Steininger committed
719

720
721
722
        Returns
        -------
        out : float, complex
Theo Steininger's avatar
Theo Steininger committed
723

724
        """
725
726
727
        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
728

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

732
        if spaces is None:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
733
            return np.vdot(y.val.flatten(), x.val.flatten())
734
735
736
        else:
            # create a diagonal operator which is capable of taking care of the
            # axes-matching
737
            from .operators.diagonal_operator import DiagonalOperator
738
739
740
741
742
743
            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
744

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

Theo Steininger's avatar
Theo Steininger committed
748
749
750
        Returns
        -------
        norm : scalar
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
751
            The L2-norm of the field values.
csongor's avatar
csongor committed
752
753

        """
Martin Reinecke's avatar
fixes    
Martin Reinecke committed
754
        return np.sqrt(np.abs(self.vdot(x=self)))
csongor's avatar
csongor committed
755
756

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

759
760
761
        Parameters
        ----------
        inplace : boolean
Theo Steininger's avatar
Theo Steininger committed
762
            Decides whether the conjugation should be performed inplace.
Theo Steininger's avatar
Theo Steininger committed
763

764
765
766
767
        Returns
        -------
        cc : field
            The complex conjugated field.
csongor's avatar
csongor committed
768
769
770

        """
        if inplace:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
771
            self.imag *= -1
772
            return self
csongor's avatar
csongor committed
773
        else:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
774
            return Field(self.domain, np.conj(self.val), self.dtype)
csongor's avatar
csongor committed
775

Theo Steininger's avatar
Theo Steininger committed
776
    # ---General unary/contraction methods---
777

Theo Steininger's avatar
Theo Steininger committed
778
779
    def __pos__(self):
        return self.copy()
780

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

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

787
    def _contraction_helper(self, op, spaces):
Theo Steininger's avatar
Theo Steininger committed
788
        if spaces is None:
789
790
791
            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
792

793
        axes_list = tuple(self.domain_axes[sp_index] for sp_index in spaces)
794
795

        try:
Theo Steininger's avatar
Theo Steininger committed
796
            axes_list = reduce(lambda x, y: x+y, axes_list)
797
        except TypeError:
Theo Steininger's avatar
Theo Steininger committed
798
            axes_list = ()
csongor's avatar
csongor committed
799

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
800
        # perform the contraction on the data
Theo Steininger's avatar
Theo Steininger committed
801
802
        data = self.get_val(copy=False)
        data = getattr(data, op)(axis=axes_list)
csongor's avatar
csongor committed
803

Theo Steininger's avatar
Theo Steininger committed
804
805
806
        # 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
807
        else:
Theo Steininger's avatar
Theo Steininger committed
808
            return_domain = tuple(self.domain[i]
Martin Reinecke's avatar
Martin Reinecke committed
809
                                  for i in range(len(self.domain))
Theo Steininger's avatar
Theo Steininger committed
810
                                  if i not in spaces)
811

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

814
815
    def sum(self, spaces=None):
        return self._contraction_helper('sum', spaces)
csongor's avatar
csongor committed
816

817
818
    def prod(self, spaces=None):
        return self._contraction_helper('prod', spaces)
csongor's avatar
csongor committed
819

820
821
    def all(self, spaces=None):
        return self._contraction_helper('all', spaces)
csongor's avatar
csongor committed
822

823
824
    def any(self, spaces=None):
        return self._contraction_helper('any', spaces)
csongor's avatar
csongor committed
825

826
827
    def min(self, spaces=None):
        return self._contraction_helper('min', spaces)
csongor's avatar
csongor committed
828

829
830
    def nanmin(self, spaces=None):
        return self._contraction_helper('nanmin', spaces)
csongor's avatar
csongor committed
831

832
833
    def max(self, spaces=None):
        return self._contraction_helper('max', spaces)
csongor's avatar
csongor committed
834

835
836
    def nanmax(self, spaces=None):
        return self._contraction_helper('nanmax', spaces)
csongor's avatar
csongor committed
837

838
839
    def mean(self, spaces=None):
        return self._contraction_helper('mean', spaces)