field.py 23.2 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
        return utilities.parse_domain(domain)

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

Martin Reinecke's avatar
Martin Reinecke committed
121
    # MR: this needs some rethinking ... do we need to have at least float64?
Martin Reinecke's avatar
Martin Reinecke committed
122
123
    @staticmethod
    def _infer_dtype(dtype, val):
Martin Reinecke's avatar
Martin Reinecke committed
124
125
126
127
128
        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)
129

130
    # ---Factory methods---
131

Martin Reinecke's avatar
Martin Reinecke committed
132
133
    @staticmethod
    def from_random(random_type, domain, dtype=np.float64, **kwargs):
134
135
136
137
138
139
140
        """ 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
141

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

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

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

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

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

162
163
    # ---Powerspectral methods---

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

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

        Parameters
        ----------
Theo Steininger's avatar
Theo Steininger committed
176
        spaces : int *optional*
Martin Reinecke's avatar
Martin Reinecke committed
177
            The subspace for which the powerspectrum shall be computed.
Theo Steininger's avatar
Theo Steininger committed
178
179
180
            (default : None).
        binbounds : array-like *optional*
            Inner bounds of the bins (default : None).
Martin Reinecke's avatar
Martin Reinecke committed
181
            if binbounds==None : bins are inferred.
182
183
184
185
186
187
188
189
190
191
        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
192

193
194
        Raise
        -----
Martin Reinecke's avatar
Martin Reinecke committed
195
196
        TypeError
            Raised if any of the input field's domains is not harmonic
Theo Steininger's avatar
Theo Steininger committed
197

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

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

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

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

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

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

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

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

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

Martin Reinecke's avatar
Martin Reinecke committed
239
    @staticmethod
Martin Reinecke's avatar
Martin Reinecke committed
240
241
242
243
244
245
246
247
    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
248

Martin Reinecke's avatar
Martin Reinecke committed
249
        power_spectrum = utilities.bincount_axis(pindex, weights=field.val,
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
250
                                                 axis=axes)
Martin Reinecke's avatar
Martin Reinecke committed
251
252
253
254
255
256
        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)
257

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

        return spec.real
355

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

490
        return new_field
csongor's avatar
csongor committed
491

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Martin Reinecke's avatar
Martin Reinecke committed
632
633
        tval = getattr(self.val, op)(other)
        return self if tval is self.val else Field(self.domain, tval)
csongor's avatar
csongor committed
634
635

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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