field.py 23 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
        spaces : int *optional*
Martin Reinecke's avatar
Martin Reinecke committed
178
            The subspace for which the powerspectrum shall be computed.
Theo Steininger's avatar
Theo Steininger committed
179
180
181
            (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
        Raise
        -----
Martin Reinecke's avatar
Martin Reinecke committed
196
197
        TypeError
            Raised if any of the input field's domains is not harmonic
Theo Steininger's avatar
Theo Steininger committed
198

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

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

209
        """
Theo Steininger's avatar
Theo Steininger committed
210

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

        # check if the `spaces` input is valid
219
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
220
            spaces = range(len(self.domain))
Martin Reinecke's avatar
Martin Reinecke committed
221
222
        else:
            spaces = utilities.cast_iseq_to_tuple(spaces)
223
224

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

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

        for space_index in spaces:
Martin Reinecke's avatar
Martin Reinecke committed
233
234
            parts = [self._single_power_analyze(field=part,
                                                idx=space_index,
235
                                                binbounds=binbounds)
236
                     for part in parts]
237

238
        return parts[0] + 1j*parts[1] if keep_phase_information else parts[0]
239

Martin Reinecke's avatar
Martin Reinecke committed
240
    @staticmethod
Martin Reinecke's avatar
Martin Reinecke committed
241
242
243
244
245
246
247
248
    def _single_power_analyze(field, idx, binbounds):
        power_domain = PowerSpace(field.domain[idx], binbounds)
        pindex = power_domain.pindex
        axes = field.domain_axes[idx]
        new_pindex_shape = [1] * len(field.shape)
        for i, ax in enumerate(axes):
            new_pindex_shape[ax] = pindex.shape[i]
        pindex = np.broadcast_to(pindex.reshape(new_pindex_shape), field.shape)
Theo Steininger's avatar
Theo Steininger committed
249

Martin Reinecke's avatar
Martin Reinecke committed
250
        power_spectrum = utilities.bincount_axis(pindex, weights=field.val,
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
251
                                                 axis=axes)
Martin Reinecke's avatar
Martin Reinecke committed
252
253
254
255
256
257
        new_rho_shape = [1] * len(power_spectrum.shape)
        new_rho_shape[axes[0]] = len(power_domain.rho)
        power_spectrum /= power_domain.rho.reshape(new_rho_shape)
        result_domain = list(field.domain)
        result_domain[idx] = power_domain
        return Field(result_domain, power_spectrum)
258

Martin Reinecke's avatar
Martin Reinecke committed
259
    def _compute_spec(self, spaces):
260
261
        if spaces is None:
            spaces = range(len(self.domain))
Martin Reinecke's avatar
Martin Reinecke committed
262
263
        else:
            spaces = utilities.cast_iseq_to_tuple(spaces)
264
265
266
267
268
269
270
271
272
273
274

        # 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
275
276
277
278
279
280
281
282
283
            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)
284
285

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

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

291
292
293
        Parameters
        ----------
        spaces : {tuple, int, None} *optional*
Theo Steininger's avatar
Theo Steininger committed
294
295
296
            Specifies the subspace containing all the PowerSpaces which
            should be converted (default : None).
            if spaces==None : Tries to convert the whole domain.
297
        real_power : boolean *optional*
Theo Steininger's avatar
Theo Steininger committed
298
299
            Determines whether the power spectrum is treated as intrinsically
            real or complex (default : True).
300
        real_signal : boolean *optional*
Theo Steininger's avatar
Theo Steininger committed
301
302
            True will result in a purely real signal-space field
            (default : True).
Theo Steininger's avatar
Theo Steininger committed
303

304
305
306
307
        Returns
        -------
        out : Field
            The output object. A random field created with the power spectrum
Theo Steininger's avatar
Theo Steininger committed
308
            stored in the `spaces` in `self`.
309

Theo Steininger's avatar
Theo Steininger committed
310
311
312
313
314
315
        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.

316
317
318
        See Also
        --------
        power_analyze
Theo Steininger's avatar
Theo Steininger committed
319
320
321
322
323

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

324
        """
Theo Steininger's avatar
Theo Steininger committed
325

Martin Reinecke's avatar
Martin Reinecke committed
326
        spec = self._compute_spec(spaces)
327
328
329

        # create random samples: one or two, depending on whether the
        # power spectrum is real or complex
330
        result = [self.from_random('normal', mean=0., std=1.,
Martin Reinecke's avatar
Martin Reinecke committed
331
                                   domain=spec.domain,
332
333
334
335
336
337
338
                                   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
339
                             domain=spec.domain, dtype=np.float)
340
341

        # apply the rescaler to the random fields
342
        result[0] *= spec.real
343
        if not real_power:
344
            result[1] *= spec.imag
345

346
        return result[0] if real_power else result[0] + 1j*result[1]
347

348
    def power_synthesize_special(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
349
        spec = self._compute_spec(spaces)
350
351
352

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

        return spec.real
356

theos's avatar
theos committed
357
    # ---Properties---
358

theos's avatar
theos committed
359
360
    @property
    def val(self):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
361
        """ Returns the data object associated with this Field.
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
362
        No copy is made.
Theo Steininger's avatar
Theo Steininger committed
363

364
365
        Returns
        -------
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
366
        out : numpy.ndarray
367
        """
Martin Reinecke's avatar
Martin Reinecke committed
368
        return self._val
csongor's avatar
csongor committed
369

Martin Reinecke's avatar
Martin Reinecke committed
370
371
372
373
    @property
    def dtype(self):
        return self._val.dtype

374
375
    @property
    def shape(self):
Theo Steininger's avatar
Theo Steininger committed
376
        """ Returns the total shape of the Field's data array.
Theo Steininger's avatar
Theo Steininger committed
377

378
379
380
        Returns
        -------
        out : tuple
Martin Reinecke's avatar
Martin Reinecke committed
381
            The output object. The tuple contains the dimensions of the spaces
382
            in domain.
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
383
       """
Martin Reinecke's avatar
Martin Reinecke committed
384
        return self._val.shape
csongor's avatar
csongor committed
385

386
387
    @property
    def dim(self):
Theo Steininger's avatar
Theo Steininger committed
388
        """ Returns the total number of pixel-dimensions the field has.
Theo Steininger's avatar
Theo Steininger committed
389

Theo Steininger's avatar
Theo Steininger committed
390
        Effectively, all values from shape are multiplied.
Theo Steininger's avatar
Theo Steininger committed
391

392
393
394
395
396
        Returns
        -------
        out : int
            The dimension of the Field.
        """
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
397
        return self._val.size
csongor's avatar
csongor committed
398

theos's avatar
theos committed
399
400
    @property
    def total_volume(self):
Theo Steininger's avatar
Theo Steininger committed
401
402
        """ Returns the total volume of all spaces in the domain.
        """
Martin Reinecke's avatar
Martin Reinecke committed
403
        if len(self.domain) == 0:
Theo Steininger's avatar
Theo Steininger committed
404
            return 0.
Martin Reinecke's avatar
Martin Reinecke committed
405
406
        volume_tuple = tuple(sp.total_volume for sp in self.domain)
        return reduce(lambda x, y: x * y, volume_tuple)
407

Theo Steininger's avatar
Theo Steininger committed
408
409
410
411
    @property
    def real(self):
        """ The real part of the field (data is not copied).
        """
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
412
        return Field(self.domain, self.val.real)
Theo Steininger's avatar
Theo Steininger committed
413
414
415
416
417

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

theos's avatar
theos committed
420
    # ---Special unary/binary operations---
421

Martin Reinecke's avatar
Martin Reinecke committed
422
    def copy(self):
423
        """ Returns a full copy of the Field.
Theo Steininger's avatar
Theo Steininger committed
424

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

427
428
429
430
431
        Returns
        -------
        out : Field
            The output object. An identical copy of 'self'.
        """
Martin Reinecke's avatar
Martin Reinecke committed
432
        return Field(val=self, copy=True)
csongor's avatar
csongor committed
433

434
435
436
437
438
439
    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
440
        res = 1.
441
442
443
444
445
446
447
        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
448
    def weight(self, power=1, inplace=False, spaces=None):
Theo Steininger's avatar
Theo Steininger committed
449
        """ Weights the pixels of `self` with their invidual pixel-volume.
450
451
452
453

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

456
        inplace : boolean
Theo Steininger's avatar
Theo Steininger committed
457
458
            If True, `self` will be weighted and returned. Otherwise, a copy
            is made.
Theo Steininger's avatar
Theo Steininger committed
459

Theo Steininger's avatar
Theo Steininger committed
460
461
        spaces : tuple of ints
            Determines on which subspace the operation takes place.
Theo Steininger's avatar
Theo Steininger committed
462

463
464
465
        Returns
        -------
        out : Field
Theo Steininger's avatar
Theo Steininger committed
466
            The weighted field.
467
468

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

csongor's avatar
csongor committed
471
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
472
            spaces = range(len(self.domain))
Martin Reinecke's avatar
Martin Reinecke committed
473
474
        else:
            spaces = utilities.cast_iseq_to_tuple(spaces)
csongor's avatar
csongor committed
475

476
477
478
479
480
481
482
        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)
483
484
                new_shape[self.domain_axes[ind][0]:
                          self.domain_axes[ind][-1]+1] = wgt.shape
485
                wgt = wgt.reshape(new_shape)
486
                new_field *= wgt**power
487
        fct = fct**power
Martin Reinecke's avatar
Martin Reinecke committed
488
        if fct != 1.:
489
            new_field *= fct
490

491
        return new_field
csongor's avatar
csongor committed
492

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

496
497
498
        Parameters
        ----------
        x : Field
Theo Steininger's avatar
Theo Steininger committed
499
            The domain of x must contain `self.domain`
Theo Steininger's avatar
Theo Steininger committed
500

Theo Steininger's avatar
Theo Steininger committed
501
502
503
        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
504

505
        bare : boolean
Theo Steininger's avatar
Theo Steininger committed
506
            If true, no volume factors will be included in the computation.
Theo Steininger's avatar
Theo Steininger committed
507

508
509
510
        Returns
        -------
        out : float, complex
Theo Steininger's avatar
Theo Steininger committed
511

512
        """
513
514
515
        if not isinstance(x, Field):
            raise ValueError("The dot-partner must be an instance of " +
                             "the NIFTy field class")
theos's avatar
theos committed
516

Martin Reinecke's avatar
Martin Reinecke committed
517
        # Compute the dot respecting the fact of discrete/continuous spaces
518
519
520
521
522
523
524
525
526
527
        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
528

529
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
530
            return fct*np.vdot(y.val.ravel(), x.val.ravel())
531
532
533
        else:
            # create a diagonal operator which is capable of taking care of the
            # axes-matching
534
            from .operators.diagonal_operator import DiagonalOperator
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
535
536
            diag = DiagonalOperator(y.domain, y.conjugate(), copy=False)
            dotted = diag(x, spaces=spaces)
537
            return fct*dotted.sum(spaces=spaces)
theos's avatar
theos committed
538

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

Theo Steininger's avatar
Theo Steininger committed
542
543
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
544
        norm : float
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
545
            The L2-norm of the field values.
csongor's avatar
csongor committed
546
547

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

Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
550
    def conjugate(self):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
551
        """ Returns the complex conjugate of the field.
Theo Steininger's avatar
Theo Steininger committed
552

553
554
555
556
        Returns
        -------
        cc : field
            The complex conjugated field.
csongor's avatar
csongor committed
557
558

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

theos's avatar
theos committed
561
    # ---General unary/contraction methods---
562

theos's avatar
theos committed
563
564
    def __pos__(self):
        return self.copy()
565

theos's avatar
theos committed
566
    def __neg__(self):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
567
        return Field(self.domain, -self.val, self.dtype)
csongor's avatar
csongor committed
568

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

572
    def _contraction_helper(self, op, spaces):
theos's avatar
theos committed
573
        if spaces is None:
574
            return getattr(self.val, op)()
Martin Reinecke's avatar
Martin Reinecke committed
575
576
        else:
            spaces = utilities.cast_iseq_to_tuple(spaces)
csongor's avatar
csongor committed
577

578
        axes_list = tuple(self.domain_axes[sp_index] for sp_index in spaces)
579

Martin Reinecke's avatar
Martin Reinecke committed
580
        if len(axes_list) > 0:
theos's avatar
theos committed
581
            axes_list = reduce(lambda x, y: x+y, axes_list)
csongor's avatar
csongor committed
582

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

theos's avatar
theos committed
586
587
588
        # 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
589
        else:
theos's avatar
theos committed
590
            return_domain = tuple(self.domain[i]
Martin Reinecke's avatar
Martin Reinecke committed
591
                                  for i in range(len(self.domain))
theos's avatar
theos committed
592
                                  if i not in spaces)
593

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

596
597
    def sum(self, spaces=None):
        return self._contraction_helper('sum', spaces)
csongor's avatar
csongor committed
598

599
600
    def prod(self, spaces=None):
        return self._contraction_helper('prod', spaces)
csongor's avatar
csongor committed
601

602
603
    def all(self, spaces=None):
        return self._contraction_helper('all', spaces)
csongor's avatar
csongor committed
604

605
606
    def any(self, spaces=None):
        return self._contraction_helper('any', spaces)
csongor's avatar
csongor committed
607

608
609
    def min(self, spaces=None):
        return self._contraction_helper('min', spaces)
csongor's avatar
csongor committed
610

611
612
    def max(self, spaces=None):
        return self._contraction_helper('max', spaces)
csongor's avatar
csongor committed
613

614
615
    def mean(self, spaces=None):
        return self._contraction_helper('mean', spaces)
csongor's avatar
csongor committed
616

617
618
    def var(self, spaces=None):
        return self._contraction_helper('var', spaces)
csongor's avatar
csongor committed
619

620
621
    def std(self, spaces=None):
        return self._contraction_helper('std', spaces)
csongor's avatar
csongor committed
622

theos's avatar
theos committed
623
    # ---General binary methods---
csongor's avatar
csongor committed
624

625
    def _binary_helper(self, other, op):
csongor's avatar
csongor committed
626
        # if other is a field, make sure that the domains match
627
        if isinstance(other, Field):
628
629
            if other.domain != self.domain:
                raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
630
            return Field(self.domain, getattr(self.val, op)(other.val))
csongor's avatar
csongor committed
631

Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
632
        return Field(self.domain, getattr(self.val, op)(other))
csongor's avatar
csongor committed
633
634

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

637
    def __radd__(self, other):
theos's avatar
theos committed
638
        return self._binary_helper(other, op='__radd__')
csongor's avatar
csongor committed
639
640

    def __iadd__(self, other):
641
        return self._binary_helper(other, op='__iadd__')
csongor's avatar
csongor committed
642
643

    def __sub__(self, other):
theos's avatar
theos committed
644
        return self._binary_helper(other, op='__sub__')
csongor's avatar
csongor committed
645
646

    def __rsub__(self, other):
theos's avatar
theos committed
647
        return self._binary_helper(other, op='__rsub__')
csongor's avatar
csongor committed
648
649

    def __isub__(self, other):
650
        return self._binary_helper(other, op='__isub__')
csongor's avatar
csongor committed
651
652

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

655
    def __rmul__(self, other):
theos's avatar
theos committed
656
        return self._binary_helper(other, op='__rmul__')
csongor's avatar
csongor committed
657
658

    def __imul__(self, other):
659
        return self._binary_helper(other, op='__imul__')
csongor's avatar
csongor committed
660
661

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

Martin Reinecke's avatar
Martin Reinecke committed
664
665
666
    def __truediv__(self, other):
        return self._binary_helper(other, op='__truediv__')

csongor's avatar
csongor committed
667
    def __rdiv__(self, other):
theos's avatar
theos committed
668
        return self._binary_helper(other, op='__rdiv__')
csongor's avatar
csongor committed
669

Martin Reinecke's avatar
Martin Reinecke committed
670
671
672
    def __rtruediv__(self, other):
        return self._binary_helper(other, op='__rtruediv__')

csongor's avatar
csongor committed
673
    def __idiv__(self, other):
674
        return self._binary_helper(other, op='__idiv__')
675

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

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

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

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

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

    def __ne__(self, other):
        if other is None:
            return True
        else:
theos's avatar
theos committed
695
            return self._binary_helper(other, op='__ne__')
csongor's avatar
csongor committed
696
697
698
699
700

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

    def __ge__(self, other):
theos's avatar
theos committed
704
        return self._binary_helper(other, op='__ge__')
csongor's avatar
csongor committed
705
706

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

    def __repr__(self):
Martin Reinecke's avatar
Martin Reinecke committed
710
        return "<nifty2go.Field>"
theos's avatar
theos committed
711
712
713
714

    def __str__(self):
        minmax = [self.min(), self.max()]
        mean = self.mean()
Martin Reinecke's avatar
Martin Reinecke committed
715
        return "nifty2go.Field instance\n- domain      = " + \
theos's avatar
theos committed
716
               repr(self.domain) + \
717
               "\n- val         = " + repr(self.val) + \
theos's avatar
theos committed
718
719
               "\n  - min.,max. = " + str(minmax) + \
               "\n  - mean = " + str(mean)