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
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
177
178
179
180
        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
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
195
196
        Raise
        -----
        ValueError
            Raised if
Theo Steininger's avatar
Theo Steininger committed
197
198
                *len(domain) is != 1 when spaces==None
                *len(spaces) is != 1 if not None
199
                *the analyzed space is not harmonic
Theo Steininger's avatar
Theo Steininger committed
200

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

        return power_spectrum

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

293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
    def _prep_powersynth(self, spaces):
        # 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:
            spec = self._spec_to_rescaler(spec, i)
        return (result_domain, spec)

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

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

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

331
332
333
334
        Returns
        -------
        out : Field
            The output object. A random field created with the power spectrum
Theo Steininger's avatar
Theo Steininger committed
335
            stored in the `spaces` in `self`.
336

Theo Steininger's avatar
Theo Steininger committed
337
338
339
340
341
342
        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.

343
344
345
        See Also
        --------
        power_analyze
Theo Steininger's avatar
Theo Steininger committed
346
347
348
349
350

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

351
        """
Theo Steininger's avatar
Theo Steininger committed
352

353
        result_domain, spec = self._prep_powersynth(spaces)
354
355
356

        # create random samples: one or two, depending on whether the
        # power spectrum is real or complex
357
358
359
360
361
362
363
364
365
366
        result = [self.from_random('normal', mean=0., std=1.,
                                   domain=result_domain,
                                   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.,
                             domain=result_domain, dtype=np.float)
367
368

        # apply the rescaler to the random fields
369
        result[0] *= spec.real
370
        if not real_power:
371
            result[1] *= spec.imag
372

373
        return result[0] if real_power else result[0] + 1j*result[1]
374

375
376
377
378
379
380
381
382
    def power_synthesize_special(self, spaces=None):
        result_domain, spec = self._prep_powersynth(spaces)

        # MR: dummy call - will be removed soon
        self.from_random('normal', mean=0., std=1.,
                         domain=result_domain, dtype=np.complex)

        return spec.real
383

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
384
    def _spec_to_rescaler(self, spec, power_space_index):
385
        power_space = self.domain[power_space_index]
386

387
388
389
390
        local_blow_up = [slice(None)]*len(spec.shape)
        # it is important to count from behind, since spec potentially grows
        # with every iteration
        index = self.domain_axes[power_space_index][0]-len(self.shape)
391
        local_blow_up[index] = power_space.pindex
392
        # here, the power_spectrum is distributed into the new shape
393
        return spec[local_blow_up]
394

theos's avatar
theos committed
395
    # ---Properties---
396

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

402
403
        Returns
        -------
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
404
        out : numpy.ndarray
405
        """
Martin Reinecke's avatar
Martin Reinecke committed
406
        return self._val
csongor's avatar
csongor committed
407

Martin Reinecke's avatar
Martin Reinecke committed
408
409
410
411
    @property
    def dtype(self):
        return self._val.dtype

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

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

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

Theo Steininger's avatar
Theo Steininger committed
428
        Effectively, all values from shape are multiplied.
Theo Steininger's avatar
Theo Steininger committed
429

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

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

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

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

theos's avatar
theos committed
458
    # ---Special unary/binary operations---
459

Martin Reinecke's avatar
Martin Reinecke committed
460
    def copy(self):
461
        """ Returns a full copy of the Field.
Theo Steininger's avatar
Theo Steininger committed
462

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

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

472
473
474
475
476
477
    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
478
        res = 1.
479
480
481
482
483
484
485
        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
486
    def weight(self, power=1, inplace=False, spaces=None):
Theo Steininger's avatar
Theo Steininger committed
487
        """ Weights the pixels of `self` with their invidual pixel-volume.
488
489
490
491

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

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

Theo Steininger's avatar
Theo Steininger committed
498
499
        spaces : tuple of ints
            Determines on which subspace the operation takes place.
Theo Steininger's avatar
Theo Steininger committed
500

501
502
503
        Returns
        -------
        out : Field
Theo Steininger's avatar
Theo Steininger committed
504
            The weighted field.
505
506

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

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

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

528
        return new_field
csongor's avatar
csongor committed
529

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

533
534
535
        Parameters
        ----------
        x : Field
Theo Steininger's avatar
Theo Steininger committed
536
            The domain of x must contain `self.domain`
Theo Steininger's avatar
Theo Steininger committed
537

Theo Steininger's avatar
Theo Steininger committed
538
539
540
        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
541

542
        bare : boolean
Theo Steininger's avatar
Theo Steininger committed
543
            If true, no volume factors will be included in the computation.
Theo Steininger's avatar
Theo Steininger committed
544

545
546
547
        Returns
        -------
        out : float, complex
Theo Steininger's avatar
Theo Steininger committed
548

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

Martin Reinecke's avatar
Martin Reinecke committed
554
        # Compute the dot respecting the fact of discrete/continuous spaces
555
556
557
558
559
560
561
562
563
564
        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
565

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

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

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

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

Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
587
    def conjugate(self):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
588
        """ Returns the complex conjugate of the field.
Theo Steininger's avatar
Theo Steininger committed
589

590
591
592
593
        Returns
        -------
        cc : field
            The complex conjugated field.
csongor's avatar
csongor committed
594
595

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

theos's avatar
theos committed
598
    # ---General unary/contraction methods---
599

theos's avatar
theos committed
600
601
    def __pos__(self):
        return self.copy()
602

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

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

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

615
        axes_list = tuple(self.domain_axes[sp_index] for sp_index in spaces)
616

Martin Reinecke's avatar
Martin Reinecke committed
617
        if len(axes_list) > 0:
theos's avatar
theos committed
618
            axes_list = reduce(lambda x, y: x+y, axes_list)
csongor's avatar
csongor committed
619

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

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

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

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

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

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

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

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

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

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

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

657
658
    def std(self, spaces=None):
        return self._contraction_helper('std', spaces)
csongor's avatar
csongor committed
659

theos's avatar
theos committed
660
    # ---General binary methods---
csongor's avatar
csongor committed
661

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

Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
669
        return Field(self.domain, getattr(self.val, op)(other))
csongor's avatar
csongor committed
670
671

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

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

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

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

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

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

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

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

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

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

Martin Reinecke's avatar
Martin Reinecke committed
701
702
703
    def __truediv__(self, other):
        return self._binary_helper(other, op='__truediv__')

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

Martin Reinecke's avatar
Martin Reinecke committed
707
708
709
    def __rtruediv__(self, other):
        return self._binary_helper(other, op='__rtruediv__')

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

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

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

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

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

    def __le__(self, other):
theos's avatar
theos committed
726
        return self._binary_helper(other, op='__le__')
csongor's avatar
csongor committed
727
728
729
730
731

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

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

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

    def __gt__(self, other):
theos's avatar
theos committed
744
745
746
        return self._binary_helper(other, op='__gt__')

    def __repr__(self):
Martin Reinecke's avatar
Martin Reinecke committed
747
        return "<nifty2go.Field>"
theos's avatar
theos committed
748
749
750
751

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