field.py 24.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
Theo Steininger's avatar
Theo Steininger committed
13
14
15
16
17
#
# Copyright(C) 2013-2017 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
18

csongor's avatar
csongor committed
19
from __future__ import division
Martin Reinecke's avatar
Martin Reinecke committed
20
from builtins import range
csongor's avatar
csongor committed
21
import numpy as np
Martin Reinecke's avatar
Martin Reinecke committed
22
23
24
from .spaces.power_space import PowerSpace
from . import nifty_utilities as utilities
from .random import Random
Martin Reinecke's avatar
Martin Reinecke committed
25
from functools import reduce
26

csongor's avatar
csongor committed
27

Martin Reinecke's avatar
Martin Reinecke committed
28
class Field(object):
Theo Steininger's avatar
Theo Steininger committed
29
30
31
    """ The discrete representation of a continuous field over multiple spaces.

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

35
36
37
38
    Parameters
    ----------
    domain : DomainObject
        One of the space types NIFTY supports. RGSpace, GLSpace, HPSpace,
Theo Steininger's avatar
Theo Steininger committed
39
        LMSpace or PowerSpace. It might also be a FieldArray, which is
40
        an unstructured domain.
Theo Steininger's avatar
Theo Steininger committed
41

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
42
    val : scalar, numpy.ndarray, Field
43
44
45
        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
46

47
    dtype : type
Martin Reinecke's avatar
updates    
Martin Reinecke committed
48
        A numpy.type. Most common are float and complex.
Theo Steininger's avatar
Theo Steininger committed
49

50
51
52
53
    copy: boolean

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

56
57
58
59
60
61
    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
62

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

71
    """
72

theos's avatar
theos committed
73
    # ---Initialization methods---
74

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
75
    def __init__(self, domain=None, val=None, dtype=None, copy=False):
76
        self.domain = self._parse_domain(domain=domain, val=val)
77
        self.domain_axes = self._get_axes_tuple(self.domain)
78

Martin Reinecke's avatar
Martin Reinecke committed
79
        shape_tuple = tuple(sp.shape for sp in self.domain)
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
80
        if len(shape_tuple) == 0:
Martin Reinecke's avatar
Martin Reinecke committed
81
            global_shape = ()
82
        else:
Martin Reinecke's avatar
Martin Reinecke committed
83
84
            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
85
        if isinstance(val, Field):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
86
            if self.domain != val.domain:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
87
                raise ValueError("Domain mismatch")
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
88
            self._val = np.array(val.val, dtype=dtype, copy=copy)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
89
        elif (np.isscalar(val)):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
90
            self._val = np.full(global_shape, dtype=dtype, fill_value=val)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
91
        elif isinstance(val, np.ndarray):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
92
93
            if global_shape == val.shape:
                self._val = np.array(val, dtype=dtype, copy=copy)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
94
95
96
            else:
                raise ValueError("Shape mismatch")
        elif val is None:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
97
            self._val = np.empty(global_shape, dtype=dtype)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
98
99
        else:
            raise TypeError("unknown source type")
csongor's avatar
csongor committed
100

Martin Reinecke's avatar
Martin Reinecke committed
101
102
    @staticmethod
    def _parse_domain(domain, val=None):
103
        if domain is None:
104
            if isinstance(val, Field):
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
105
106
107
108
                return val.domain
            if np.isscalar(val):
                return ()  # empty domain tuple
            raise TypeError("could not infer domain from value")
Martin Reinecke's avatar
Martin Reinecke committed
109
110
111
112
113

        return utilities.parse_domain(domain)

    @staticmethod
    def _get_axes_tuple(things_with_shape):
114
        i = 0
theos's avatar
theos committed
115
116
        axes_list = []
        for thing in things_with_shape:
117
            nax = len(thing.shape)
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
118
            axes_list += [tuple(range(i, i+nax))]
119
            i += nax
theos's avatar
theos committed
120
        return tuple(axes_list)
121

Martin Reinecke's avatar
Martin Reinecke committed
122
    # MR: this needs some rethinking ... do we need to have at least float64?
Martin Reinecke's avatar
Martin Reinecke committed
123
124
    @staticmethod
    def _infer_dtype(dtype, val):
Martin Reinecke's avatar
Martin Reinecke committed
125
126
127
128
129
        if val is None or dtype is not None:
            return np.result_type(dtype, np.float64)
        if isinstance(val, Field):
            return val.dtype
        return np.result_type(val, np.float64)
130

131
    # ---Factory methods---
132

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

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

143
144
        domain : DomainObject
            The domain of the output random field
Theo Steininger's avatar
Theo Steininger committed
145

146
147
        dtype : type
            The datatype of the output random field
Theo Steininger's avatar
Theo Steininger committed
148

149
150
151
152
153
154
155
        Returns
        -------
        out : Field
            The output object.

        See Also
        --------
156
        power_synthesize
157
        """
Theo Steininger's avatar
Theo Steininger committed
158

159
        generator_function = getattr(Random, random_type)
160
        return Field(domain=domain, val=generator_function(dtype=dtype,
Martin Reinecke's avatar
Martin Reinecke committed
161
                     shape=utilities.domains2shape(domain), **kwargs))
162

163
164
    # ---Powerspectral methods---

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

Theo Steininger's avatar
Theo Steininger committed
169
170
171
        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
172
        harmonic space. The resulting field has the same units as the initial
Theo Steininger's avatar
Theo Steininger committed
173
        field, corresponding to the square root of the power spectrum.
174
175
176

        Parameters
        ----------
Theo Steininger's avatar
Theo Steininger committed
177
178
179
180
181
        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
182
            if binbounds==None : bins are inferred.
183
184
185
186
187
188
189
190
191
192
        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
193

194
195
196
197
        Raise
        -----
        ValueError
            Raised if
Theo Steininger's avatar
Theo Steininger committed
198
199
                *len(domain) is != 1 when spaces==None
                *len(spaces) is != 1 if not None
200
                *the analyzed space is not harmonic
Theo Steininger's avatar
Theo Steininger committed
201

202
203
        Returns
        -------
Theo Steininger's avatar
Theo Steininger committed
204
        out : Field
Martin Reinecke's avatar
Martin Reinecke committed
205
            The output object. Its domain is a PowerSpace and it contains
206
207
208
209
210
            the power spectrum of 'self's field.

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

212
        """
Theo Steininger's avatar
Theo Steininger committed
213

Theo Steininger's avatar
Theo Steininger committed
214
        # check if all spaces in `self.domain` are either harmonic or
215
216
217
        # power_space instances
        for sp in self.domain:
            if not sp.harmonic and not isinstance(sp, PowerSpace):
Martin Reinecke's avatar
Martin Reinecke committed
218
219
                raise TypeError("Field has a space in `domain` which is "
                                "neither harmonic nor a PowerSpace.")
220
221

        # check if the `spaces` input is valid
222
223
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
224
            spaces = range(len(self.domain))
225
226

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

229
        if keep_phase_information:
230
            parts = [self.real*self.real, self.imag*self.imag]
231
        else:
232
            parts = [self.real*self.real + self.imag*self.imag]
233
234

        for space_index in spaces:
235
236
237
            parts = [self._single_power_analyze(work_field=part,
                                                space_index=space_index,
                                                binbounds=binbounds)
238
                     for part in parts]
239

240
        return parts[0] + 1j*parts[1] if keep_phase_information else parts[0]
241

Martin Reinecke's avatar
Martin Reinecke committed
242
243
    @staticmethod
    def _single_power_analyze(work_field, space_index, binbounds):
244
        if not work_field.domain[space_index].harmonic:
Martin Reinecke's avatar
Martin Reinecke committed
245
            raise ValueError("The analyzed space must be harmonic.")
246

247
248
249
250
251
252
        # 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.

253
        harmonic_domain = work_field.domain[space_index]
254
        power_domain = PowerSpace(harmonic_partner=harmonic_domain,
Theo Steininger's avatar
Theo Steininger committed
255
                                  binbounds=binbounds)
Martin Reinecke's avatar
Martin Reinecke committed
256
        power_spectrum = Field._calculate_power_spectrum(
257
                                field_val=work_field.val,
Martin Reinecke's avatar
Martin Reinecke committed
258
                                pdomain=power_domain,
259
                                axes=work_field.domain_axes[space_index])
260
261

        # create the result field and put power_spectrum into it
262
        result_domain = list(work_field.domain)
263
264
        result_domain[space_index] = power_domain

Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
265
        return Field(domain=result_domain, val=power_spectrum,
Martin Reinecke's avatar
updates    
Martin Reinecke committed
266
                     dtype=power_spectrum.dtype)
267

Martin Reinecke's avatar
Martin Reinecke committed
268
269
    @staticmethod
    def _calculate_power_spectrum(field_val, pdomain, axes=None):
Martin Reinecke's avatar
Martin Reinecke committed
270
        pindex = pdomain.pindex
271
        if axes is not None:
Martin Reinecke's avatar
Martin Reinecke committed
272
            pindex = Field._shape_up_pindex(pindex, field_val.shape, axes)
Theo Steininger's avatar
Theo Steininger committed
273

Martin Reinecke's avatar
Martin Reinecke committed
274
        power_spectrum = utilities.bincount_axis(pindex, weights=field_val,
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
275
                                                 axis=axes)
Martin Reinecke's avatar
Martin Reinecke committed
276
        rho = pdomain.rho
277
        if axes is not None:
278
            new_rho_shape = [1] * len(power_spectrum.shape)
279
280
281
282
283
284
            new_rho_shape[axes[0]] = len(rho)
            rho = rho.reshape(new_rho_shape)
        power_spectrum /= rho

        return power_spectrum

285
    @staticmethod
Martin Reinecke's avatar
Martin Reinecke committed
286
    def _shape_up_pindex(pindex, target_shape, axes):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
287
        semiscaled_local_shape = [1] * len(target_shape)
Martin Reinecke's avatar
Martin Reinecke committed
288
289
        for i, ax in enumerate(axes):
            semiscaled_local_shape[ax] = pindex.shape[i]
Martin Reinecke's avatar
Martin Reinecke committed
290
        result_obj = np.empty(target_shape, dtype=pindex.dtype)
Martin Reinecke's avatar
Martin Reinecke committed
291
        result_obj[()] = pindex.reshape(semiscaled_local_shape)
292
293
        return result_obj

Martin Reinecke's avatar
Martin Reinecke committed
294
    def _compute_spec(self, spaces):
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
        # check if the `spaces` input is valid
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
        if spaces is None:
            spaces = range(len(self.domain))

        # create the result domain
        result_domain = list(self.domain)
        for i in spaces:
            if not isinstance(self.domain[i], PowerSpace):
                raise ValueError("A PowerSpace is needed for field "
                                 "synthetization.")
            result_domain[i] = self.domain[i].harmonic_partner

        spec = np.sqrt(self.val)
        for i in spaces:
Martin Reinecke's avatar
Martin Reinecke committed
310
311
312
313
314
315
316
317
318
            power_space = self.domain[i]
            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[i][0]-len(self.shape)
            local_blow_up[index] = power_space.pindex
            # here, the power_spectrum is distributed into the new shape
            spec = spec[local_blow_up]
        return Field(result_domain, val=spec)
319
320

    def power_synthesize(self, spaces=None, real_power=True, real_signal=True):
Theo Steininger's avatar
Theo Steininger committed
321
        """ Yields a sampled field with `self`**2 as its power spectrum.
Theo Steininger's avatar
Theo Steininger committed
322

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

326
327
328
        Parameters
        ----------
        spaces : {tuple, int, None} *optional*
Theo Steininger's avatar
Theo Steininger committed
329
330
331
            Specifies the subspace containing all the PowerSpaces which
            should be converted (default : None).
            if spaces==None : Tries to convert the whole domain.
332
        real_power : boolean *optional*
Theo Steininger's avatar
Theo Steininger committed
333
334
            Determines whether the power spectrum is treated as intrinsically
            real or complex (default : True).
335
        real_signal : boolean *optional*
Theo Steininger's avatar
Theo Steininger committed
336
337
            True will result in a purely real signal-space field
            (default : True).
Theo Steininger's avatar
Theo Steininger committed
338

339
340
341
342
        Returns
        -------
        out : Field
            The output object. A random field created with the power spectrum
Theo Steininger's avatar
Theo Steininger committed
343
            stored in the `spaces` in `self`.
344

Theo Steininger's avatar
Theo Steininger committed
345
346
347
348
349
350
        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.

351
352
353
        See Also
        --------
        power_analyze
Theo Steininger's avatar
Theo Steininger committed
354
355
356
357
358

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

359
        """
Theo Steininger's avatar
Theo Steininger committed
360

Martin Reinecke's avatar
Martin Reinecke committed
361
        spec = self._compute_spec(spaces)
362
363
364

        # create random samples: one or two, depending on whether the
        # power spectrum is real or complex
365
        result = [self.from_random('normal', mean=0., std=1.,
Martin Reinecke's avatar
Martin Reinecke committed
366
                                   domain=spec.domain,
367
368
369
370
371
372
373
                                   dtype=np.float if real_signal
                                   else np.complex)
                  for x in range(1 if real_power else 2)]

        # MR: dummy call - will be removed soon
        if real_signal:
            self.from_random('normal', mean=0., std=1.,
Martin Reinecke's avatar
Martin Reinecke committed
374
                             domain=spec.domain, dtype=np.float)
375
376

        # apply the rescaler to the random fields
377
        result[0] *= spec.real
378
        if not real_power:
379
            result[1] *= spec.imag
380

381
        return result[0] if real_power else result[0] + 1j*result[1]
382

383
    def power_synthesize_special(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
384
        spec = self._compute_spec(spaces)
385
386
387

        # MR: dummy call - will be removed soon
        self.from_random('normal', mean=0., std=1.,
Martin Reinecke's avatar
Martin Reinecke committed
388
                         domain=spec.domain, dtype=np.complex)
389
390

        return spec.real
391

theos's avatar
theos committed
392
    # ---Properties---
393

theos's avatar
theos committed
394
395
    @property
    def val(self):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
396
        """ Returns the data object associated with this Field.
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
397
        No copy is made.
Theo Steininger's avatar
Theo Steininger committed
398

399
400
        Returns
        -------
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
401
        out : numpy.ndarray
402
        """
Martin Reinecke's avatar
Martin Reinecke committed
403
        return self._val
csongor's avatar
csongor committed
404

Martin Reinecke's avatar
Martin Reinecke committed
405
406
407
408
    @property
    def dtype(self):
        return self._val.dtype

409
410
    @property
    def shape(self):
Theo Steininger's avatar
Theo Steininger committed
411
        """ Returns the total shape of the Field's data array.
Theo Steininger's avatar
Theo Steininger committed
412

413
414
415
        Returns
        -------
        out : tuple
Martin Reinecke's avatar
Martin Reinecke committed
416
            The output object. The tuple contains the dimensions of the spaces
417
            in domain.
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
418
       """
Martin Reinecke's avatar
Martin Reinecke committed
419
        return self._val.shape
csongor's avatar
csongor committed
420

421
422
    @property
    def dim(self):
Theo Steininger's avatar
Theo Steininger committed
423
        """ Returns the total number of pixel-dimensions the field has.
Theo Steininger's avatar
Theo Steininger committed
424

Theo Steininger's avatar
Theo Steininger committed
425
        Effectively, all values from shape are multiplied.
Theo Steininger's avatar
Theo Steininger committed
426

427
428
429
430
431
        Returns
        -------
        out : int
            The dimension of the Field.
        """
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
432
        return self._val.size
csongor's avatar
csongor committed
433

theos's avatar
theos committed
434
435
    @property
    def total_volume(self):
Theo Steininger's avatar
Theo Steininger committed
436
437
        """ Returns the total volume of all spaces in the domain.
        """
Martin Reinecke's avatar
Martin Reinecke committed
438
        if len(self.domain) == 0:
Theo Steininger's avatar
Theo Steininger committed
439
            return 0.
Martin Reinecke's avatar
Martin Reinecke committed
440
441
        volume_tuple = tuple(sp.total_volume for sp in self.domain)
        return reduce(lambda x, y: x * y, volume_tuple)
442

Theo Steininger's avatar
Theo Steininger committed
443
444
445
446
    @property
    def real(self):
        """ The real part of the field (data is not copied).
        """
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
447
        return Field(self.domain, self.val.real)
Theo Steininger's avatar
Theo Steininger committed
448
449
450
451
452

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

theos's avatar
theos committed
455
    # ---Special unary/binary operations---
456

Martin Reinecke's avatar
Martin Reinecke committed
457
    def copy(self):
458
        """ Returns a full copy of the Field.
Theo Steininger's avatar
Theo Steininger committed
459

Martin Reinecke's avatar
Martin Reinecke committed
460
        The returned object will be an identical copy of the original Field.
Theo Steininger's avatar
Theo Steininger committed
461

462
463
464
465
466
        Returns
        -------
        out : Field
            The output object. An identical copy of 'self'.
        """
Martin Reinecke's avatar
Martin Reinecke committed
467
        return Field(val=self, copy=True)
csongor's avatar
csongor committed
468

469
470
471
472
473
474
    def scalar_weight(self, spaces=None):
        if np.isscalar(spaces):
            return self.domain[spaces].scalar_weight()

        if spaces is None:
            spaces = range(len(self.domain))
Martin Reinecke's avatar
Martin Reinecke committed
475
        res = 1.
476
477
478
479
480
481
482
        for i in spaces:
            tmp = self.domain[i].scalar_weight()
            if tmp is None:
                return None
            res *= tmp
        return res

theos's avatar
theos committed
483
    def weight(self, power=1, inplace=False, spaces=None):
Theo Steininger's avatar
Theo Steininger committed
484
        """ Weights the pixels of `self` with their invidual pixel-volume.
485
486
487
488

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

491
        inplace : boolean
Theo Steininger's avatar
Theo Steininger committed
492
493
            If True, `self` will be weighted and returned. Otherwise, a copy
            is made.
Theo Steininger's avatar
Theo Steininger committed
494

Theo Steininger's avatar
Theo Steininger committed
495
496
        spaces : tuple of ints
            Determines on which subspace the operation takes place.
Theo Steininger's avatar
Theo Steininger committed
497

498
499
500
        Returns
        -------
        out : Field
Theo Steininger's avatar
Theo Steininger committed
501
            The weighted field.
502
503

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

506
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
csongor's avatar
csongor committed
507
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
508
            spaces = range(len(self.domain))
csongor's avatar
csongor committed
509

510
511
512
513
514
515
516
        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)
517
518
                new_shape[self.domain_axes[ind][0]:
                          self.domain_axes[ind][-1]+1] = wgt.shape
519
                wgt = wgt.reshape(new_shape)
520
                new_field *= wgt**power
521
        fct = fct**power
Martin Reinecke's avatar
Martin Reinecke committed
522
        if fct != 1.:
523
            new_field *= fct
524

525
        return new_field
csongor's avatar
csongor committed
526

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

530
531
532
        Parameters
        ----------
        x : Field
Theo Steininger's avatar
Theo Steininger committed
533
            The domain of x must contain `self.domain`
Theo Steininger's avatar
Theo Steininger committed
534

Theo Steininger's avatar
Theo Steininger committed
535
536
537
        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
538

539
        bare : boolean
Theo Steininger's avatar
Theo Steininger committed
540
            If true, no volume factors will be included in the computation.
Theo Steininger's avatar
Theo Steininger committed
541

542
543
544
        Returns
        -------
        out : float, complex
Theo Steininger's avatar
Theo Steininger committed
545

546
        """
547
548
549
        if not isinstance(x, Field):
            raise ValueError("The dot-partner must be an instance of " +
                             "the NIFTy field class")
theos's avatar
theos committed
550

Martin Reinecke's avatar
Martin Reinecke committed
551
        # Compute the dot respecting the fact of discrete/continuous spaces
552
553
554
555
556
557
558
559
560
561
        fct = 1.
        if bare:
            y = self
        else:
            tmp = self.scalar_weight(spaces)
            if tmp is None:
                y = self.weight(power=1)
            else:
                y = self
                fct = tmp
theos's avatar
theos committed
562

563
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
564
            return fct*np.vdot(y.val.ravel(), x.val.ravel())
565
566
567
        else:
            # create a diagonal operator which is capable of taking care of the
            # axes-matching
568
            from .operators.diagonal_operator import DiagonalOperator
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
569
570
            diag = DiagonalOperator(y.domain, y.conjugate(), copy=False)
            dotted = diag(x, spaces=spaces)
571
            return fct*dotted.sum(spaces=spaces)
theos's avatar
theos committed
572

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

Theo Steininger's avatar
Theo Steininger committed
576
577
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
578
        norm : float
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
579
            The L2-norm of the field values.
csongor's avatar
csongor committed
580
581

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

Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
584
    def conjugate(self):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
585
        """ Returns the complex conjugate of the field.
Theo Steininger's avatar
Theo Steininger committed
586

587
588
589
590
        Returns
        -------
        cc : field
            The complex conjugated field.
csongor's avatar
csongor committed
591
592

        """
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
593
        return Field(self.domain, self.val.conjugate(), self.dtype)
csongor's avatar
csongor committed
594

theos's avatar
theos committed
595
    # ---General unary/contraction methods---
596

theos's avatar
theos committed
597
598
    def __pos__(self):
        return self.copy()
599

theos's avatar
theos committed
600
    def __neg__(self):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
601
        return Field(self.domain, -self.val, self.dtype)
csongor's avatar
csongor committed
602

theos's avatar
theos committed
603
    def __abs__(self):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
604
        return Field(self.domain, np.abs(self.val), self.dtype)
csongor's avatar
csongor committed
605

606
    def _contraction_helper(self, op, spaces):
theos's avatar
theos committed
607
        if spaces is None:
608
609
610
            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
611

612
        axes_list = tuple(self.domain_axes[sp_index] for sp_index in spaces)
613

Martin Reinecke's avatar
Martin Reinecke committed
614
        if len(axes_list) > 0:
theos's avatar
theos committed
615
            axes_list = reduce(lambda x, y: x+y, axes_list)
csongor's avatar
csongor committed
616

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

theos's avatar
theos committed
620
621
622
        # 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
623
        else:
theos's avatar
theos committed
624
            return_domain = tuple(self.domain[i]
Martin Reinecke's avatar
Martin Reinecke committed
625
                                  for i in range(len(self.domain))
theos's avatar
theos committed
626
                                  if i not in spaces)
627

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

630
631
    def sum(self, spaces=None):
        return self._contraction_helper('sum', spaces)
csongor's avatar
csongor committed
632

633
634
    def prod(self, spaces=None):
        return self._contraction_helper('prod', spaces)
csongor's avatar
csongor committed
635

636
637
    def all(self, spaces=None):
        return self._contraction_helper('all', spaces)
csongor's avatar
csongor committed
638

639
640
    def any(self, spaces=None):
        return self._contraction_helper('any', spaces)
csongor's avatar
csongor committed
641

642
643
    def min(self, spaces=None):
        return self._contraction_helper('min', spaces)
csongor's avatar
csongor committed
644

645
646
    def max(self, spaces=None):
        return self._contraction_helper('max', spaces)
csongor's avatar
csongor committed
647

648
649
    def mean(self, spaces=None):
        return self._contraction_helper('mean', spaces)
csongor's avatar
csongor committed
650

651
652
    def var(self, spaces=None):
        return self._contraction_helper('var', spaces)
csongor's avatar
csongor committed
653

654
655
    def std(self, spaces=None):
        return self._contraction_helper('std', spaces)
csongor's avatar
csongor committed
656

theos's avatar
theos committed
657
    # ---General binary methods---
csongor's avatar
csongor committed
658

659
    def _binary_helper(self, other, op):
csongor's avatar
csongor committed
660
        # if other is a field, make sure that the domains match
661
        if isinstance(other, Field):
662
663
            if other.domain != self.domain:
                raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
664
            return Field(self.domain, getattr(self.val, op)(other.val))
csongor's avatar
csongor committed
665

Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
666
        return Field(self.domain, getattr(self.val, op)(other))
csongor's avatar
csongor committed
667
668

    def __add__(self, other):
theos's avatar
theos committed
669
        return self._binary_helper(other, op='__add__')
670

671
    def __radd__(self, other):
theos's avatar
theos committed
672
        return self._binary_helper(other, op='__radd__')
csongor's avatar
csongor committed
673
674

    def __iadd__(self, other):
675
        return self._binary_helper(other, op='__iadd__')
csongor's avatar
csongor committed
676
677

    def __sub__(self, other):
theos's avatar
theos committed
678
        return self._binary_helper(other, op='__sub__')
csongor's avatar
csongor committed
679
680

    def __rsub__(self, other):
theos's avatar
theos committed
681
        return self._binary_helper(other, op='__rsub__')
csongor's avatar
csongor committed
682
683

    def __isub__(self, other):
684
        return self._binary_helper(other, op='__isub__')
csongor's avatar
csongor committed
685
686

    def __mul__(self, other):
theos's avatar
theos committed
687
        return self._binary_helper(other, op='__mul__')
688

689
    def __rmul__(self, other):
theos's avatar
theos committed
690
        return self._binary_helper(other, op='__rmul__')
csongor's avatar
csongor committed
691
692

    def __imul__(self, other):
693
        return self._binary_helper(other, op='__imul__')
csongor's avatar
csongor committed
694
695

    def __div__(self, other):
theos's avatar
theos committed
696
        return self._binary_helper(other, op='__div__')
csongor's avatar
csongor committed
697

Martin Reinecke's avatar
Martin Reinecke committed
698
699
700
    def __truediv__(self, other):
        return self._binary_helper(other, op='__truediv__')

csongor's avatar
csongor committed
701
    def __rdiv__(self, other):
theos's avatar
theos committed
702
        return self._binary_helper(other, op='__rdiv__')
csongor's avatar
csongor committed
703

Martin Reinecke's avatar
Martin Reinecke committed
704
705
706
    def __rtruediv__(self, other):
        return self._binary_helper(other, op='__rtruediv__')

csongor's avatar
csongor committed
707
    def __idiv__(self, other):
708
        return self._binary_helper(other, op='__idiv__')
709

csongor's avatar
csongor committed
710
    def __pow__(self, other):
theos's avatar
theos committed
711
        return self._binary_helper(other, op='__pow__')
csongor's avatar
csongor committed
712
713

    def __rpow__(self, other):
theos's avatar
theos committed
714
        return self._binary_helper(other, op='__rpow__')
csongor's avatar
csongor committed
715
716

    def __ipow__(self, other):
717
        return self._binary_helper(other, op='__ipow__')
csongor's avatar
csongor committed
718
719

    def __lt__(self, other):
theos's avatar
theos committed
720
        return self._binary_helper(other, op='__lt__')
csongor's avatar
csongor committed
721
722

    def __le__(self, other):
theos's avatar
theos committed
723
        return self._binary_helper(other, op='__le__')
csongor's avatar
csongor committed
724
725
726
727
728

    def __ne__(self, other):
        if other is None:
            return True
        else:
theos's avatar
theos committed
729
            return self._binary_helper(other, op='__ne__')
csongor's avatar
csongor committed
730
731
732
733
734

    def __eq__(self, other):
        if other is None:
            return False
        else:
theos's avatar
theos committed
735
            return self._binary_helper(other, op='__eq__')
csongor's avatar
csongor committed
736
737

    def __ge__(self, other):
theos's avatar
theos committed
738
        return self._binary_helper(other, op='__ge__')
csongor's avatar
csongor committed
739
740

    def __gt__(self, other):
theos's avatar
theos committed
741
742
743
        return self._binary_helper(other, op='__gt__')

    def __repr__(self):
Martin Reinecke's avatar
Martin Reinecke committed
744
        return "<nifty2go.Field>"
theos's avatar
theos committed
745
746
747
748

    def __str__(self):
        minmax = [self.min(), self.max()]
        mean = self.mean()
Martin Reinecke's avatar
Martin Reinecke committed
749
        return "nifty2go.Field instance\n- domain      = " + \
theos's avatar
theos committed
750
               repr(self.domain) + \
751
               "\n- val         = " + repr(self.val) + \
theos's avatar
theos committed
752
753
               "\n  - min.,max. = " + str(minmax) + \
               "\n  - mean = " + str(mean)