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

csongor's avatar
csongor committed
19
from __future__ import division
Martin Reinecke's avatar
Martin Reinecke committed
20
from builtins import range
21

csongor's avatar
csongor committed
22
23
import numpy as np

Martin Reinecke's avatar
Martin Reinecke committed
24
from .domain_object import DomainObject
25

Martin Reinecke's avatar
Martin Reinecke committed
26
from .spaces.power_space import PowerSpace
csongor's avatar
csongor committed
27

Martin Reinecke's avatar
Martin Reinecke committed
28
29
from . import nifty_utilities as utilities
from .random import Random
Martin Reinecke's avatar
Martin Reinecke committed
30
from functools import reduce
31

csongor's avatar
csongor committed
32

Martin Reinecke's avatar
Martin Reinecke committed
33
class Field(object):
Theo Steininger's avatar
Theo Steininger committed
34
35
36
    """ The discrete representation of a continuous field over multiple spaces.

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

40
41
42
43
    Parameters
    ----------
    domain : DomainObject
        One of the space types NIFTY supports. RGSpace, GLSpace, HPSpace,
Theo Steininger's avatar
Theo Steininger committed
44
        LMSpace or PowerSpace. It might also be a FieldArray, which is
45
        an unstructured domain.
Theo Steininger's avatar
Theo Steininger committed
46

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
47
    val : scalar, numpy.ndarray, Field
48
49
50
        The values the array should contain after init. A scalar input will
        fill the whole array with this scalar. If an array is provided the
        array's dimensions must match the domain's.
Theo Steininger's avatar
Theo Steininger committed
51

52
    dtype : type
Martin Reinecke's avatar
updates    
Martin Reinecke committed
53
        A numpy.type. Most common are float and complex.
Theo Steininger's avatar
Theo Steininger committed
54

55
56
57
58
    copy: boolean

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

61
62
63
64
65
66
    domain : DomainObject
        See Parameters.
    domain_axes : tuple of tuples
        Enumerates the axes of the Field
    dtype : type
        Contains the datatype stored in the Field.
Theo Steininger's avatar
Theo Steininger committed
67

68
69
70
71
72
73
74
    Raise
    -----
    TypeError
        Raised if
            *the given domain contains something that is not a DomainObject
             instance
            *val is an array that has a different dimension than the domain
Theo Steininger's avatar
Theo Steininger committed
75

76
    """
77

theos's avatar
theos committed
78
    # ---Initialization methods---
79

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
80
    def __init__(self, domain=None, val=None, dtype=None, copy=False):
81
        self.domain = self._parse_domain(domain=domain, val=val)
82
        self.domain_axes = self._get_axes_tuple(self.domain)
83

Martin Reinecke's avatar
Martin Reinecke committed
84
        shape_tuple = tuple(sp.shape for sp in self.domain)
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
85
        if len(shape_tuple) == 0:
Martin Reinecke's avatar
Martin Reinecke committed
86
            global_shape = ()
87
        else:
Martin Reinecke's avatar
Martin Reinecke committed
88
89
            global_shape = reduce(lambda x, y: x + y, shape_tuple)
        dtype = self._infer_dtype(dtype=dtype, val=val)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
90
        if isinstance(val, Field):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
91
            if self.domain != val.domain:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
92
                raise ValueError("Domain mismatch")
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
93
            self._val = np.array(val.val, dtype=dtype, copy=copy)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
94
        elif (np.isscalar(val)):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
95
            self._val = np.full(global_shape, dtype=dtype, fill_value=val)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
96
        elif isinstance(val, np.ndarray):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
97
98
            if global_shape == val.shape:
                self._val = np.array(val, dtype=dtype, copy=copy)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
99
100
101
            else:
                raise ValueError("Shape mismatch")
        elif val is None:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
102
            self._val = np.empty(global_shape, dtype=dtype)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
103
104
        else:
            raise TypeError("unknown source type")
csongor's avatar
csongor committed
105

Martin Reinecke's avatar
Martin Reinecke committed
106
107
    @staticmethod
    def _parse_domain(domain, val=None):
108
        if domain is None:
109
            if isinstance(val, Field):
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
110
111
112
113
                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
114
115
116
117
118

        return utilities.parse_domain(domain)

    @staticmethod
    def _get_axes_tuple(things_with_shape):
119
        i = 0
theos's avatar
theos committed
120
121
        axes_list = []
        for thing in things_with_shape:
122
            nax = len(thing.shape)
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
123
            axes_list += [tuple(range(i, i+nax))]
124
            i += nax
theos's avatar
theos committed
125
        return tuple(axes_list)
126

Martin Reinecke's avatar
Martin Reinecke committed
127
128
    @staticmethod
    def _infer_dtype(dtype, val):
Martin Reinecke's avatar
Martin Reinecke committed
129
130
131
132
133
        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)
134

135
    # ---Factory methods---
136

Martin Reinecke's avatar
Martin Reinecke committed
137
138
    @staticmethod
    def from_random(random_type, domain, dtype=np.float64, **kwargs):
139
140
141
142
143
144
145
        """ 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
146

147
148
        domain : DomainObject
            The domain of the output random field
Theo Steininger's avatar
Theo Steininger committed
149

150
151
        dtype : type
            The datatype of the output random field
Theo Steininger's avatar
Theo Steininger committed
152

153
154
155
156
157
158
159
        Returns
        -------
        out : Field
            The output object.

        See Also
        --------
160
        power_synthesize
161
        """
Theo Steininger's avatar
Theo Steininger committed
162

163
        generator_function = getattr(Random, random_type)
Martin Reinecke's avatar
Martin Reinecke committed
164
        return Field(domain=domain,
Martin Reinecke's avatar
Martin Reinecke committed
165
                     val=generator_function(dtype=dtype,
Martin Reinecke's avatar
Martin Reinecke committed
166
                     shape=utilities.domains2shape(domain), **kwargs))
167

168
169
    # ---Powerspectral methods---

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

Theo Steininger's avatar
Theo Steininger committed
174
175
176
        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
177
        harmonic space. The resulting field has the same units as the initial
Theo Steininger's avatar
Theo Steininger committed
178
        field, corresponding to the square root of the power spectrum.
179
180
181

        Parameters
        ----------
Theo Steininger's avatar
Theo Steininger committed
182
183
184
185
186
        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
187
            if binbounds==None : bins are inferred.
188
189
190
191
192
193
194
195
196
197
        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
198

199
200
201
202
        Raise
        -----
        ValueError
            Raised if
Theo Steininger's avatar
Theo Steininger committed
203
204
                *len(domain) is != 1 when spaces==None
                *len(spaces) is != 1 if not None
205
                *the analyzed space is not harmonic
Theo Steininger's avatar
Theo Steininger committed
206

207
208
        Returns
        -------
Theo Steininger's avatar
Theo Steininger committed
209
        out : Field
Martin Reinecke's avatar
Martin Reinecke committed
210
            The output object. Its domain is a PowerSpace and it contains
211
212
213
214
215
            the power spectrum of 'self's field.

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

217
        """
Theo Steininger's avatar
Theo Steininger committed
218

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

        # check if the `spaces` input is valid
227
228
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
229
            spaces = range(len(self.domain))
230
231

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

234
        if keep_phase_information:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
235
            parts = [self.real.copy(), self.imag.copy()]
236
237
238
239
        else:
            parts = [self]

        parts = [abs(part)**2 for part in parts]
240
241

        for space_index in spaces:
242
243
            parts = [self._single_power_analyze(
                                work_field=part,
244
                                space_index=space_index,
245
246
                                binbounds=binbounds)
                     for part in parts]
247

248
        if keep_phase_information:
Martin Reinecke's avatar
updates    
Martin Reinecke committed
249
            return parts[0] + 1j*parts[1]
250
        else:
Martin Reinecke's avatar
updates    
Martin Reinecke committed
251
            return parts[0]
252

Martin Reinecke's avatar
Martin Reinecke committed
253
254
    @staticmethod
    def _single_power_analyze(work_field, space_index, binbounds):
255
        if not work_field.domain[space_index].harmonic:
Martin Reinecke's avatar
Martin Reinecke committed
256
            raise ValueError("The analyzed space must be harmonic.")
257

258
259
260
261
262
263
        # 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.

264
        harmonic_domain = work_field.domain[space_index]
265
        power_domain = PowerSpace(harmonic_partner=harmonic_domain,
Theo Steininger's avatar
Theo Steininger committed
266
                                  binbounds=binbounds)
Martin Reinecke's avatar
Martin Reinecke committed
267
        power_spectrum = Field._calculate_power_spectrum(
268
                                field_val=work_field.val,
Martin Reinecke's avatar
Martin Reinecke committed
269
                                pdomain=power_domain,
270
                                axes=work_field.domain_axes[space_index])
271
272

        # create the result field and put power_spectrum into it
273
        result_domain = list(work_field.domain)
274
275
        result_domain[space_index] = power_domain

Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
276
        return Field(domain=result_domain, val=power_spectrum,
Martin Reinecke's avatar
updates    
Martin Reinecke committed
277
                     dtype=power_spectrum.dtype)
278

Martin Reinecke's avatar
Martin Reinecke committed
279
280
    @staticmethod
    def _calculate_power_spectrum(field_val, pdomain, axes=None):
281

Martin Reinecke's avatar
Martin Reinecke committed
282
        pindex = pdomain.pindex
283
        if axes is not None:
Martin Reinecke's avatar
Martin Reinecke committed
284
            pindex = Field._shape_up_pindex(pindex, field_val.shape, axes)
Theo Steininger's avatar
Theo Steininger committed
285

Martin Reinecke's avatar
Martin Reinecke committed
286
        power_spectrum = utilities.bincount_axis(pindex, weights=field_val,
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
287
                                                 axis=axes)
Martin Reinecke's avatar
Martin Reinecke committed
288
        rho = pdomain.rho
289
290
291
292
293
294
295
296
        if axes is not None:
            new_rho_shape = [1, ] * len(power_spectrum.shape)
            new_rho_shape[axes[0]] = len(rho)
            rho = rho.reshape(new_rho_shape)
        power_spectrum /= rho

        return power_spectrum

297
    @staticmethod
Martin Reinecke's avatar
Martin Reinecke committed
298
    def _shape_up_pindex(pindex, target_shape, axes):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
299
        semiscaled_local_shape = [1] * len(target_shape)
Martin Reinecke's avatar
Martin Reinecke committed
300
301
        for i, ax in enumerate(axes):
            semiscaled_local_shape[ax] = pindex.shape[i]
Martin Reinecke's avatar
Martin Reinecke committed
302
        result_obj = np.empty(target_shape, dtype=pindex.dtype)
Martin Reinecke's avatar
Martin Reinecke committed
303
        result_obj[()] = pindex.reshape(semiscaled_local_shape)
304
305
        return result_obj

306
    def power_synthesize(self, spaces=None, real_power=True, real_signal=True,
Martin Reinecke's avatar
Martin Reinecke committed
307
                         mean=0., std=1.):
Theo Steininger's avatar
Theo Steininger committed
308
        """ Yields a sampled field with `self`**2 as its power spectrum.
Theo Steininger's avatar
Theo Steininger committed
309

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

313
314
315
        Parameters
        ----------
        spaces : {tuple, int, None} *optional*
Theo Steininger's avatar
Theo Steininger committed
316
317
318
            Specifies the subspace containing all the PowerSpaces which
            should be converted (default : None).
            if spaces==None : Tries to convert the whole domain.
319
        real_power : boolean *optional*
Theo Steininger's avatar
Theo Steininger committed
320
321
            Determines whether the power spectrum is treated as intrinsically
            real or complex (default : True).
322
        real_signal : boolean *optional*
Theo Steininger's avatar
Theo Steininger committed
323
324
325
326
327
328
            True will result in a purely real signal-space field
            (default : True).
        mean : float *optional*
            The mean of the Gaussian noise field which is used for the Field
            synthetization (default : None).
            if mean==None : mean will be set to 0
329
        std : float *optional*
Theo Steininger's avatar
Theo Steininger committed
330
331
332
            The standard deviation of the Gaussian noise field which is used
            for the Field synthetization (default : None).
            if std==None : std will be set to 1
Theo Steininger's avatar
Theo Steininger committed
333

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

Theo Steininger's avatar
Theo Steininger committed
340
341
342
343
344
345
        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.

346
347
348
        See Also
        --------
        power_analyze
Theo Steininger's avatar
Theo Steininger committed
349
350
351
352
353

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

354
        """
Theo Steininger's avatar
Theo Steininger committed
355

356
357
        # check if the `spaces` input is valid
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
Theo Steininger's avatar
Theo Steininger committed
358
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
359
            spaces = range(len(self.domain))
Theo Steininger's avatar
Theo Steininger committed
360

Martin Reinecke's avatar
Martin Reinecke committed
361
362
        for i in spaces:
            if not isinstance(self.domain[i], PowerSpace):
363
364
                raise ValueError("A PowerSpace is needed for field "
                                 "synthetization.")
365
366
367

        # create the result domain
        result_domain = list(self.domain)
Martin Reinecke's avatar
Martin Reinecke committed
368
369
        for i in spaces:
            result_domain[i] = self.domain[i].harmonic_partner
370
371
372

        # create random samples: one or two, depending on whether the
        # power spectrum is real or complex
373
374
        result_list = [self.__class__.from_random(
                             'normal',
375
376
377
                             mean=mean,
                             std=std,
                             domain=result_domain,
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
378
                             dtype=np.complex)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
379
                       for x in range(1 if real_power else 2)]
380
381
382
383
384

        # from now on extract the values from the random fields for further
        # processing without killing the fields.
        # if the signal-space field should be real, hermitianize the field
        # components
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
385
        spec = np.sqrt(self.val)
386
        for power_space_index in spaces:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
387
            spec = self._spec_to_rescaler(spec, power_space_index)
388
389

        # apply the rescaler to the random fields
Martin Reinecke's avatar
Martin Reinecke committed
390
        result_list[0] *= spec.real
391
        if not real_power:
Martin Reinecke's avatar
Martin Reinecke committed
392
            result_list[1] *= spec.imag
393

394
        if real_signal:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
395
            result_list = [i.real*np.sqrt(2.) for i in result_list]
396
397

        if real_power:
Martin Reinecke's avatar
Martin Reinecke committed
398
            return result_list[0]
399
        else:
Martin Reinecke's avatar
Martin Reinecke committed
400
            return result_list[0] + 1j*result_list[1]
401

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
402
    def _spec_to_rescaler(self, spec, power_space_index):
403
        power_space = self.domain[power_space_index]
404

405
406
407
408
        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)
409
        local_blow_up[index] = power_space.pindex
410
        # here, the power_spectrum is distributed into the new shape
411
        return spec[local_blow_up]
412

theos's avatar
theos committed
413
    # ---Properties---
414

theos's avatar
theos committed
415
416
    @property
    def val(self):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
417
        """ Returns the data object associated with this Field.
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
418
        No copy is made.
Theo Steininger's avatar
Theo Steininger committed
419

420
421
        Returns
        -------
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
422
        out : numpy.ndarray
423
        """
Martin Reinecke's avatar
Martin Reinecke committed
424
        return self._val
csongor's avatar
csongor committed
425

Martin Reinecke's avatar
Martin Reinecke committed
426
427
428
429
    @property
    def dtype(self):
        return self._val.dtype

430
431
    @property
    def shape(self):
Theo Steininger's avatar
Theo Steininger committed
432
        """ Returns the total shape of the Field's data array.
Theo Steininger's avatar
Theo Steininger committed
433

434
435
436
        Returns
        -------
        out : tuple
Martin Reinecke's avatar
Martin Reinecke committed
437
            The output object. The tuple contains the dimensions of the spaces
438
            in domain.
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
439
       """
Martin Reinecke's avatar
Martin Reinecke committed
440
        return self._val.shape
csongor's avatar
csongor committed
441

442
443
    @property
    def dim(self):
Theo Steininger's avatar
Theo Steininger committed
444
        """ Returns the total number of pixel-dimensions the field has.
Theo Steininger's avatar
Theo Steininger committed
445

Theo Steininger's avatar
Theo Steininger committed
446
        Effectively, all values from shape are multiplied.
Theo Steininger's avatar
Theo Steininger committed
447

448
449
450
451
452
        Returns
        -------
        out : int
            The dimension of the Field.
        """
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
453
        return self._val.size
csongor's avatar
csongor committed
454

theos's avatar
theos committed
455
456
    @property
    def total_volume(self):
Theo Steininger's avatar
Theo Steininger committed
457
458
        """ Returns the total volume of all spaces in the domain.
        """
Martin Reinecke's avatar
Martin Reinecke committed
459
        if len(self.domain) == 0:
Theo Steininger's avatar
Theo Steininger committed
460
            return 0.
Martin Reinecke's avatar
Martin Reinecke committed
461
462
        volume_tuple = tuple(sp.total_volume for sp in self.domain)
        return reduce(lambda x, y: x * y, volume_tuple)
463

Theo Steininger's avatar
Theo Steininger committed
464
465
466
467
    @property
    def real(self):
        """ The real part of the field (data is not copied).
        """
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
468
        return Field(self.domain, self.val.real)
Theo Steininger's avatar
Theo Steininger committed
469
470
471
472
473

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

theos's avatar
theos committed
476
    # ---Special unary/binary operations---
477

Martin Reinecke's avatar
Martin Reinecke committed
478
    def copy(self):
479
        """ Returns a full copy of the Field.
Theo Steininger's avatar
Theo Steininger committed
480

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

483
484
485
486
487
        Returns
        -------
        out : Field
            The output object. An identical copy of 'self'.
        """
Martin Reinecke's avatar
Martin Reinecke committed
488
        return Field(val=self, copy=True)
csongor's avatar
csongor committed
489

490
491
492
493
494
495
    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
496
        res = 1.
497
498
499
500
501
502
503
        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
504
    def weight(self, power=1, inplace=False, spaces=None):
Theo Steininger's avatar
Theo Steininger committed
505
        """ Weights the pixels of `self` with their invidual pixel-volume.
506
507
508
509

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

512
        inplace : boolean
Theo Steininger's avatar
Theo Steininger committed
513
514
            If True, `self` will be weighted and returned. Otherwise, a copy
            is made.
Theo Steininger's avatar
Theo Steininger committed
515

Theo Steininger's avatar
Theo Steininger committed
516
517
        spaces : tuple of ints
            Determines on which subspace the operation takes place.
Theo Steininger's avatar
Theo Steininger committed
518

519
520
521
        Returns
        -------
        out : Field
Theo Steininger's avatar
Theo Steininger committed
522
            The weighted field.
523
524

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

527
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
csongor's avatar
csongor committed
528
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
529
            spaces = range(len(self.domain))
csongor's avatar
csongor committed
530

531
532
533
534
535
536
537
538
539
        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)
                new_shape[self.domain_axes[ind][0]:self.domain_axes[ind][-1]+1]=wgt.shape
                wgt = wgt.reshape(new_shape)
540
                new_field *= wgt**power
541
        fct = fct**power
Martin Reinecke's avatar
Martin Reinecke committed
542
        if fct != 1.:
543
            new_field *= fct
544

545
        return new_field
csongor's avatar
csongor committed
546

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

550
551
552
        Parameters
        ----------
        x : Field
Theo Steininger's avatar
Theo Steininger committed
553
            The domain of x must contain `self.domain`
Theo Steininger's avatar
Theo Steininger committed
554

Theo Steininger's avatar
Theo Steininger committed
555
556
557
        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
558

559
        bare : boolean
Theo Steininger's avatar
Theo Steininger committed
560
            If true, no volume factors will be included in the computation.
Theo Steininger's avatar
Theo Steininger committed
561

562
563
564
        Returns
        -------
        out : float, complex
Theo Steininger's avatar
Theo Steininger committed
565

566
        """
567
568
569
        if not isinstance(x, Field):
            raise ValueError("The dot-partner must be an instance of " +
                             "the NIFTy field class")
theos's avatar
theos committed
570

Martin Reinecke's avatar
Martin Reinecke committed
571
        # Compute the dot respecting the fact of discrete/continuous spaces
572
573
574
575
576
577
578
579
580
581
        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
582

583
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
584
            return fct*np.vdot(y.val.ravel(), x.val.ravel())
585
586
587
        else:
            # create a diagonal operator which is capable of taking care of the
            # axes-matching
588
            from .operators.diagonal_operator import DiagonalOperator
589
590
591
592
593
            diagonal = y.val.conjugate()
            diagonalOperator = DiagonalOperator(domain=y.domain,
                                                diagonal=diagonal,
                                                copy=False)
            dotted = diagonalOperator(x, spaces=spaces)
594
            return fct*dotted.sum(spaces=spaces)
theos's avatar
theos committed
595

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

Theo Steininger's avatar
Theo Steininger committed
599
600
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
601
        norm : float
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
602
            The L2-norm of the field values.
csongor's avatar
csongor committed
603
604

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

    def conjugate(self, inplace=False):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
608
        """ Returns the complex conjugate of the field.
Theo Steininger's avatar
Theo Steininger committed
609

610
611
612
        Parameters
        ----------
        inplace : boolean
Theo Steininger's avatar
Theo Steininger committed
613
            Decides whether the conjugation should be performed inplace.
Theo Steininger's avatar
Theo Steininger committed
614

615
616
617
618
        Returns
        -------
        cc : field
            The complex conjugated field.
csongor's avatar
csongor committed
619
620
621

        """
        if inplace:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
622
            self.imag *= -1
623
            return self
csongor's avatar
csongor committed
624
        else:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
625
            return Field(self.domain, np.conj(self.val), self.dtype)
csongor's avatar
csongor committed
626

theos's avatar
theos committed
627
    # ---General unary/contraction methods---
628

theos's avatar
theos committed
629
630
    def __pos__(self):
        return self.copy()
631

theos's avatar
theos committed
632
    def __neg__(self):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
633
        return Field(self.domain, -self.val, self.dtype)
csongor's avatar
csongor committed
634

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

638
    def _contraction_helper(self, op, spaces):
theos's avatar
theos committed
639
        if spaces is None:
640
641
642
            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
643

644
        axes_list = tuple(self.domain_axes[sp_index] for sp_index in spaces)
645

Martin Reinecke's avatar
Martin Reinecke committed
646
        if len(axes_list) > 0:
theos's avatar
theos committed
647
            axes_list = reduce(lambda x, y: x+y, axes_list)
csongor's avatar
csongor committed
648

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

theos's avatar
theos committed
652
653
654
        # 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
655
        else:
theos's avatar
theos committed
656
            return_domain = tuple(self.domain[i]
Martin Reinecke's avatar
Martin Reinecke committed
657
                                  for i in range(len(self.domain))
theos's avatar
theos committed
658
                                  if i not in spaces)
659

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

662
663
    def sum(self, spaces=None):
        return self._contraction_helper('sum', spaces)
csongor's avatar
csongor committed
664

665
666
    def prod(self, spaces=None):
        return self._contraction_helper('prod', spaces)
csongor's avatar
csongor committed
667

668
669
    def all(self, spaces=None):
        return self._contraction_helper('all', spaces)
csongor's avatar
csongor committed
670

671
672
    def any(self, spaces=None):
        return self._contraction_helper('any', spaces)
csongor's avatar
csongor committed
673

674
675
    def min(self, spaces=None):
        return self._contraction_helper('min', spaces)
csongor's avatar
csongor committed
676

677
678
    def max(self, spaces=None):
        return self._contraction_helper('max', spaces)
csongor's avatar
csongor committed
679

680
681
    def mean(self, spaces=None):
        return self._contraction_helper('mean', spaces)
csongor's avatar
csongor committed
682

683
684
    def var(self, spaces=None):
        return self._contraction_helper('var', spaces)
csongor's avatar
csongor committed
685

686
687
    def std(self, spaces=None):
        return self._contraction_helper('std', spaces)
csongor's avatar
csongor committed
688

theos's avatar
theos committed
689
    # ---General binary methods---
csongor's avatar
csongor committed
690

691
    def _binary_helper(self, other, op):
csongor's avatar
csongor committed
692
        # if other is a field, make sure that the domains match
693
        if isinstance(other, Field):
694
695
            if other.domain != self.domain:
                raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
696
            return Field(self.domain, getattr(self.val, op)(other.val))
csongor's avatar
csongor committed
697

Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
698
        return Field(self.domain, getattr(self.val, op)(other))
csongor's avatar
csongor committed
699
700

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

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

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

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

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

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

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

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

    def __imul__(self, other):
725
        return self._binary_helper(other, op='__imul__')
csongor's avatar
csongor committed
726
727

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

Martin Reinecke's avatar
Martin Reinecke committed
730
731
732
    def __truediv__(self, other):
        return self._binary_helper(other, op='__truediv__')

csongor's avatar
csongor committed
733
    def __rdiv__(self, other):
theos's avatar
theos committed
734
        return self._binary_helper(other, op='__rdiv__')
csongor's avatar
csongor committed
735

Martin Reinecke's avatar
Martin Reinecke committed
736
737
738
    def __rtruediv__(self, other):
        return self._binary_helper(other, op='__rtruediv__')

csongor's avatar
csongor committed
739
    def __idiv__(self, other):
740
        return self._binary_helper(other, op='__idiv__')
741

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

    def __rpow__(self, other):
theos's avatar
theos committed
746
        return self._binary_helper(other, op='__rpow__')
csongor's avatar
csongor committed
747
748

    def __ipow__(self, other):
749
        return self._binary_helper(other, op='__ipow__')
csongor's avatar
csongor committed
750
751

    def __lt__(self, other):
theos's avatar
theos committed
752
        return self._binary_helper(other, op='__lt__')
csongor's avatar
csongor committed
753
754

    def __le__(self, other):
theos's avatar
theos committed
755
        return self._binary_helper(other, op='__le__')
csongor's avatar
csongor committed
756
757
758
759
760

    def __ne__(self, other):
        if other is None:
            return True
        else:
theos's avatar
theos committed
761
            return self._binary_helper(other, op='__ne__')
csongor's avatar
csongor committed
762
763
764
765
766

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

    def __ge__(self, other):
theos's avatar
theos committed
770
        return self._binary_helper(other, op='__ge__')
csongor's avatar
csongor committed
771
772

    def __gt__(self, other):
theos's avatar
theos committed
773
774
775
        return self._binary_helper(other, op='__gt__')

    def __repr__(self):
Martin Reinecke's avatar
Martin Reinecke committed
776
        return "<nifty2go.Field>"
theos's avatar
theos committed
777
778
779
780

    def __str__(self):
        minmax = [self.min(), self.max()]
        mean = self.mean()
Martin Reinecke's avatar
Martin Reinecke committed
781
        return "nifty2go.Field instance\n- domain      = " + \
theos's avatar
theos committed
782
               repr(self.domain) + \
783
               "\n- val         = " + repr(self.val) + \
theos's avatar
theos committed
784
785
               "\n  - min.,max. = " + str(minmax) + \
               "\n  - mean = " + str(mean)