field.py 23.3 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

Martin Reinecke's avatar
Martin Reinecke committed
19
from __future__ import division, print_function
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 .domain_tuple import DomainTuple
Martin Reinecke's avatar
Martin Reinecke committed
26
from functools import reduce
27
from . import dobj
28

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

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

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

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

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

51
52
53
54
    copy: boolean

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

Martin Reinecke's avatar
Martin Reinecke committed
57
    domain : DomainTuple
58
59
60
        See Parameters.
    dtype : type
        Contains the datatype stored in the Field.
Theo Steininger's avatar
Theo Steininger committed
61

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

70
    """
71

Theo Steininger's avatar
Theo Steininger committed
72
    # ---Initialization methods---
73

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

Martin Reinecke's avatar
Martin Reinecke committed
77
        dtype = self._infer_dtype(dtype=dtype, val=val)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
78
        if isinstance(val, Field):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
79
            if self.domain != val.domain:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
80
                raise ValueError("Domain mismatch")
81
            self._val = dobj.from_object(val.val, dtype=dtype, copy=copy)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
82
        elif (np.isscalar(val)):
83
84
            self._val = dobj.full(self.domain.shape, dtype=dtype, fill_value=val)
        elif isinstance(val, dobj.data_object):
Martin Reinecke's avatar
Martin Reinecke committed
85
            if self.domain.shape == val.shape:
86
                self._val = dobj.from_object(val, dtype=dtype, copy=copy)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
87
88
89
            else:
                raise ValueError("Shape mismatch")
        elif val is None:
90
            self._val = dobj.empty(self.domain.shape, dtype=dtype)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
91
92
        else:
            raise TypeError("unknown source type")
csongor's avatar
csongor committed
93

94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
    @staticmethod
    def full(domain, val, dtype=None):
        if not np.isscalar(val):
            raise TypeError("val must be a scalar")
        return Field(DomainTuple.make(domain), val, dtype)

    @staticmethod
    def ones(domain, dtype=None):
        return Field(DomainTuple.make(domain), 1., dtype)

    @staticmethod
    def zeros(domain, dtype=None):
        return Field(DomainTuple.make(domain), 0., dtype)

    @staticmethod
    def empty(domain, dtype=None):
        return Field(DomainTuple.make(domain), None, dtype)

    @staticmethod
    def full_like(field, val, dtype=None):
        if not isinstance(field, Field):
            raise TypeError("field must be of Field type")
        return Field.full(field.domain, val, dtype)

    @staticmethod
    def zeros_like(field, dtype=None):
        if not isinstance(field, Field):
            raise TypeError("field must be of Field type")
        if dtype is None:
            dtype = field.dtype
        return Field.zeros(field.domain, dtype)

    @staticmethod
    def ones_like(field, dtype=None):
        if not isinstance(field, Field):
            raise TypeError("field must be of Field type")
        if dtype is None:
            dtype = field.dtype
        return Field.ones(field.domain, dtype)

    @staticmethod
    def empty_like(field, dtype=None):
        if not isinstance(field, Field):
            raise TypeError("field must be of Field type")
        if dtype is None:
            dtype = field.dtype
        return Field.empty(field.domain, dtype)

Martin Reinecke's avatar
Martin Reinecke committed
142
143
    @staticmethod
    def _parse_domain(domain, val=None):
144
        if domain is None:
145
            if isinstance(val, Field):
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
146
147
                return val.domain
            if np.isscalar(val):
148
                return DomainTuple.make(())  # empty domain tuple
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
149
            raise TypeError("could not infer domain from value")
Martin Reinecke's avatar
Martin Reinecke committed
150
        return DomainTuple.make(domain)
151

Martin Reinecke's avatar
Martin Reinecke committed
152
    # MR: this needs some rethinking ... do we need to have at least float64?
Martin Reinecke's avatar
Martin Reinecke committed
153
154
    @staticmethod
    def _infer_dtype(dtype, val):
Martin Reinecke's avatar
Martin Reinecke committed
155
156
157
158
159
        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)
160

161
    # ---Factory methods---
162

Martin Reinecke's avatar
Martin Reinecke committed
163
164
    @staticmethod
    def from_random(random_type, domain, dtype=np.float64, **kwargs):
165
166
167
168
169
170
171
        """ 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
172

173
174
        domain : DomainObject
            The domain of the output random field
Theo Steininger's avatar
Theo Steininger committed
175

176
177
        dtype : type
            The datatype of the output random field
Theo Steininger's avatar
Theo Steininger committed
178

179
180
181
182
183
184
185
        Returns
        -------
        out : Field
            The output object.

        See Also
        --------
186
        power_synthesize
187
        """
Theo Steininger's avatar
Theo Steininger committed
188

Martin Reinecke's avatar
Martin Reinecke committed
189
        domain = DomainTuple.make(domain)
190
        generator_function = getattr(Random, random_type)
191
        return Field(domain=domain, val=generator_function(dtype=dtype,
Martin Reinecke's avatar
Martin Reinecke committed
192
                     shape=domain.shape, **kwargs))
193

194
195
    # ---Powerspectral methods---

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

Theo Steininger's avatar
Theo Steininger committed
200
201
202
        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
203
        harmonic space. The resulting field has the same units as the initial
Theo Steininger's avatar
Theo Steininger committed
204
        field, corresponding to the square root of the power spectrum.
205
206
207

        Parameters
        ----------
Theo Steininger's avatar
Theo Steininger committed
208
        spaces : int *optional*
Martin Reinecke's avatar
Martin Reinecke committed
209
            The subspace for which the powerspectrum shall be computed.
Theo Steininger's avatar
Theo Steininger committed
210
211
212
            (default : None).
        binbounds : array-like *optional*
            Inner bounds of the bins (default : None).
Martin Reinecke's avatar
Martin Reinecke committed
213
            if binbounds==None : bins are inferred.
214
215
216
217
218
219
220
221
222
223
        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
224

225
226
        Raise
        -----
Martin Reinecke's avatar
Martin Reinecke committed
227
228
        TypeError
            Raised if any of the input field's domains is not harmonic
Theo Steininger's avatar
Theo Steininger committed
229

230
231
        Returns
        -------
Theo Steininger's avatar
Theo Steininger committed
232
        out : Field
Martin Reinecke's avatar
Martin Reinecke committed
233
            The output object. Its domain is a PowerSpace and it contains
234
235
236
237
238
            the power spectrum of 'self's field.

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

240
        """
Theo Steininger's avatar
Theo Steininger committed
241

Theo Steininger's avatar
Theo Steininger committed
242
        # check if all spaces in `self.domain` are either harmonic or
243
244
245
        # power_space instances
        for sp in self.domain:
            if not sp.harmonic and not isinstance(sp, PowerSpace):
Martin Reinecke's avatar
Martin Reinecke committed
246
247
                print("WARNING: Field has a space in `domain` which is "
                      "neither harmonic nor a PowerSpace.")
248
249

        # check if the `spaces` input is valid
250
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
251
            spaces = range(len(self.domain))
Martin Reinecke's avatar
Martin Reinecke committed
252
253
        else:
            spaces = utilities.cast_iseq_to_tuple(spaces)
254
255

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

258
        if keep_phase_information:
259
            parts = [self.real*self.real, self.imag*self.imag]
260
        else:
261
            parts = [self.real*self.real + self.imag*self.imag]
262

Martin Reinecke's avatar
Martin Reinecke committed
263
        parts = [ part.weight(1,spaces) for part in parts ]
264
        for space_index in spaces:
Martin Reinecke's avatar
Martin Reinecke committed
265
266
            parts = [self._single_power_analyze(field=part,
                                                idx=space_index,
267
                                                binbounds=binbounds)
268
                     for part in parts]
269

270
        return parts[0] + 1j*parts[1] if keep_phase_information else parts[0]
271

Martin Reinecke's avatar
Martin Reinecke committed
272
    @staticmethod
Martin Reinecke's avatar
Martin Reinecke committed
273
    def _single_power_analyze(field, idx, binbounds):
274
        from .operators.power_projection_operator import PowerProjectionOperator
Martin Reinecke's avatar
Martin Reinecke committed
275
        power_domain = PowerSpace(field.domain[idx], binbounds)
276
        ppo = PowerProjectionOperator(field.domain, power_domain, idx)
Martin Reinecke's avatar
bug fix    
Martin Reinecke committed
277
        return ppo(field.weight(-1))
278

Martin Reinecke's avatar
Martin Reinecke committed
279
    def _compute_spec(self, spaces):
280
281
        from .operators.power_projection_operator import PowerProjectionOperator
        from .basic_arithmetics import sqrt
282
283
        if spaces is None:
            spaces = range(len(self.domain))
Martin Reinecke's avatar
Martin Reinecke committed
284
285
        else:
            spaces = utilities.cast_iseq_to_tuple(spaces)
286
287
288

        # create the result domain
        result_domain = list(self.domain)
289
290

        spec = sqrt(self)
291
292
        for i in spaces:
            result_domain[i] = self.domain[i].harmonic_partner
293
            ppo = PowerProjectionOperator(result_domain, self.domain[i], i)
294
            spec = ppo.adjoint_times(spec)
295

296
        return spec
297
298

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

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

304
305
306
        Parameters
        ----------
        spaces : {tuple, int, None} *optional*
Theo Steininger's avatar
Theo Steininger committed
307
308
309
            Specifies the subspace containing all the PowerSpaces which
            should be converted (default : None).
            if spaces==None : Tries to convert the whole domain.
310
        real_power : boolean *optional*
Theo Steininger's avatar
Theo Steininger committed
311
312
            Determines whether the power spectrum is treated as intrinsically
            real or complex (default : True).
313
        real_signal : boolean *optional*
Theo Steininger's avatar
Theo Steininger committed
314
315
            True will result in a purely real signal-space field
            (default : True).
Theo Steininger's avatar
Theo Steininger committed
316

317
318
319
320
        Returns
        -------
        out : Field
            The output object. A random field created with the power spectrum
Theo Steininger's avatar
Theo Steininger committed
321
            stored in the `spaces` in `self`.
322

Theo Steininger's avatar
Theo Steininger committed
323
324
325
326
327
328
        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.

329
330
331
        See Also
        --------
        power_analyze
Theo Steininger's avatar
Theo Steininger committed
332
333
334
335
336

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

337
        """
Theo Steininger's avatar
Theo Steininger committed
338

Martin Reinecke's avatar
Martin Reinecke committed
339
        spec = self._compute_spec(spaces)
340
341
342

        # create random samples: one or two, depending on whether the
        # power spectrum is real or complex
343
        result = [self.from_random('normal', mean=0., std=1.,
Martin Reinecke's avatar
Martin Reinecke committed
344
                                   domain=spec.domain,
345
346
347
348
349
350
351
                                   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
352
                             domain=spec.domain, dtype=np.float)
353
354

        # apply the rescaler to the random fields
355
        result[0] *= spec.real
356
        if not real_power:
357
            result[1] *= spec.imag
358

359
        return result[0] if real_power else result[0] + 1j*result[1]
360

361
    def power_synthesize_special(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
362
        spec = self._compute_spec(spaces)
363
364
365

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

        return spec.real
369

Theo Steininger's avatar
Theo Steininger committed
370
    # ---Properties---
371

Theo Steininger's avatar
Theo Steininger committed
372
373
    @property
    def val(self):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
374
        """ Returns the data object associated with this Field.
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
375
        No copy is made.
Theo Steininger's avatar
Theo Steininger committed
376

377
378
        Returns
        -------
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
379
        out : numpy.ndarray
380
        """
Martin Reinecke's avatar
Martin Reinecke committed
381
        return self._val
csongor's avatar
csongor committed
382

Martin Reinecke's avatar
Martin Reinecke committed
383
384
385
386
    @property
    def dtype(self):
        return self._val.dtype

387
388
    @property
    def shape(self):
Theo Steininger's avatar
Theo Steininger committed
389
        """ Returns the total shape of the Field's data array.
Theo Steininger's avatar
Theo Steininger committed
390

391
392
393
        Returns
        -------
        out : tuple
Martin Reinecke's avatar
Martin Reinecke committed
394
            The output object. The tuple contains the dimensions of the spaces
395
            in domain.
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
396
       """
Martin Reinecke's avatar
Martin Reinecke committed
397
        return self.domain.shape
csongor's avatar
csongor committed
398

399
400
    @property
    def dim(self):
Theo Steininger's avatar
Theo Steininger committed
401
        """ Returns the total number of pixel-dimensions the field has.
Theo Steininger's avatar
Theo Steininger committed
402

Theo Steininger's avatar
Theo Steininger committed
403
        Effectively, all values from shape are multiplied.
Theo Steininger's avatar
Theo Steininger committed
404

405
406
407
408
409
        Returns
        -------
        out : int
            The dimension of the Field.
        """
Martin Reinecke's avatar
Martin Reinecke committed
410
        return self.domain.dim
csongor's avatar
csongor committed
411

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

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

Theo Steininger's avatar
Theo Steininger committed
424
    # ---Special unary/binary operations---
425

Martin Reinecke's avatar
Martin Reinecke committed
426
    def copy(self):
427
        """ Returns a full copy of the Field.
Theo Steininger's avatar
Theo Steininger committed
428

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

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

438
439
    def scalar_weight(self, spaces=None):
        if np.isscalar(spaces):
440
            return self.domain[spaces].scalar_dvol()
441
442
443

        if spaces is None:
            spaces = range(len(self.domain))
Martin Reinecke's avatar
Martin Reinecke committed
444
        res = 1.
445
        for i in spaces:
446
            tmp = self.domain[i].scalar_dvol()
447
448
449
450
451
            if tmp is None:
                return None
            res *= tmp
        return res

452
    def weight(self, power=1, spaces=None, out=None):
Theo Steininger's avatar
Theo Steininger committed
453
        """ Weights the pixels of `self` with their invidual pixel-volume.
454
455
456
457

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

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

463
464
465
466
467
        out : Field or None
            if not None, the result is returned in a new Field
            otherwise the contents of "out" are overwritten with the result.
            "out" may be identical to "self"!

468
469
470
        Returns
        -------
        out : Field
Theo Steininger's avatar
Theo Steininger committed
471
            The weighted field.
472
473

        """
474
475
476
477
478
        if out is None:
            out = self.copy()
        else:
            if out is not self:
                out.copy_content_from(self)
csongor's avatar
csongor committed
479

csongor's avatar
csongor committed
480
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
481
            spaces = range(len(self.domain))
Martin Reinecke's avatar
Martin Reinecke committed
482
483
        else:
            spaces = utilities.cast_iseq_to_tuple(spaces)
csongor's avatar
csongor committed
484

485
486
        fct = 1.
        for ind in spaces:
487
            wgt = self.domain[ind].dvol()
488
489
490
            if np.isscalar(wgt):
                fct *= wgt
            else:
491
                new_shape = dobj.ones(len(self.shape), dtype=np.int)
Martin Reinecke's avatar
Martin Reinecke committed
492
493
                new_shape[self.domain.axes[ind][0]:
                          self.domain.axes[ind][-1]+1] = wgt.shape
494
                wgt = wgt.reshape(new_shape)
495
                out *= wgt**power
496
        fct = fct**power
Martin Reinecke's avatar
Martin Reinecke committed
497
        if fct != 1.:
498
            out *= fct
499

500
        return out
csongor's avatar
csongor committed
501

Martin Reinecke's avatar
Martin Reinecke committed
502
    def vdot(self, x=None, spaces=None):
Theo Steininger's avatar
Theo Steininger committed
503
        """ Computes the volume-factor-aware dot product of 'self' with x.
Theo Steininger's avatar
Theo Steininger committed
504

505
506
507
        Parameters
        ----------
        x : Field
Theo Steininger's avatar
Theo Steininger committed
508
            The domain of x must contain `self.domain`
Theo Steininger's avatar
Theo Steininger committed
509

Theo Steininger's avatar
Theo Steininger committed
510
        spaces : tuple of ints
511
512
            If the domain of `self` and `x` are not the same, `spaces` defines
            which domains of `x` are mapped to those of `self`.
Theo Steininger's avatar
Theo Steininger committed
513

514
515
516
        Returns
        -------
        out : float, complex
Theo Steininger's avatar
Theo Steininger committed
517

518
        """
519
520
521
        if not isinstance(x, Field):
            raise ValueError("The dot-partner must be an instance of " +
                             "the NIFTy field class")
Theo Steininger's avatar
Theo Steininger committed
522

Martin Reinecke's avatar
Martin Reinecke committed
523
        # Compute the dot respecting the fact of discrete/continuous spaces
Martin Reinecke's avatar
Martin Reinecke committed
524
525
        tmp = self.scalar_weight(spaces)
        if tmp is None:
526
            fct = 1.
Martin Reinecke's avatar
Martin Reinecke committed
527
            y = self.weight(power=1)
528
        else:
Martin Reinecke's avatar
Martin Reinecke committed
529
530
            y = self
            fct = tmp
Theo Steininger's avatar
Theo Steininger committed
531

532
        if spaces is None:
533
            return fct*dobj.vdot(y.val.ravel(), x.val.ravel())
534
        else:
535
536
537
538
539
540
541
542
            spaces = utilities.cast_iseq_to_tuple(spaces)
            active_axes = []
            for i in spaces:
                active_axes += self.domain.axes[i]
            res = 0.
            for sl in utilities.get_slice_list(self.shape, active_axes):
                res += dobj.vdot(y.val, x.val[sl])
            return res*fct
Theo Steininger's avatar
Theo Steininger committed
543

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

Theo Steininger's avatar
Theo Steininger committed
547
548
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
549
        norm : float
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
550
            The L2-norm of the field values.
csongor's avatar
csongor committed
551
552

        """
553
        return np.sqrt(np.abs(self.vdot(x=self)))
csongor's avatar
csongor committed
554

Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
555
    def conjugate(self):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
556
        """ Returns the complex conjugate of the field.
Theo Steininger's avatar
Theo Steininger committed
557

558
559
560
561
        Returns
        -------
        cc : field
            The complex conjugated field.
csongor's avatar
csongor committed
562
563

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

Theo Steininger's avatar
Theo Steininger committed
566
    # ---General unary/contraction methods---
567

Theo Steininger's avatar
Theo Steininger committed
568
569
    def __pos__(self):
        return self.copy()
570

Theo Steininger's avatar
Theo Steininger committed
571
    def __neg__(self):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
572
        return Field(self.domain, -self.val, self.dtype)
csongor's avatar
csongor committed
573

Theo Steininger's avatar
Theo Steininger committed
574
    def __abs__(self):
575
        return Field(self.domain, dobj.abs(self.val), self.dtype)
csongor's avatar
csongor committed
576

577
    def _contraction_helper(self, op, spaces):
Theo Steininger's avatar
Theo Steininger committed
578
        if spaces is None:
579
            return getattr(self.val, op)()
Martin Reinecke's avatar
Martin Reinecke committed
580
581
        else:
            spaces = utilities.cast_iseq_to_tuple(spaces)
csongor's avatar
csongor committed
582

Martin Reinecke's avatar
Martin Reinecke committed
583
        axes_list = tuple(self.domain.axes[sp_index] for sp_index in spaces)
584

Martin Reinecke's avatar
Martin Reinecke committed
585
        if len(axes_list) > 0:
Theo Steininger's avatar
Theo Steininger committed
586
            axes_list = reduce(lambda x, y: x+y, axes_list)
csongor's avatar
csongor committed
587

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

Theo Steininger's avatar
Theo Steininger committed
591
592
593
        # 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
594
        else:
Theo Steininger's avatar
Theo Steininger committed
595
            return_domain = tuple(self.domain[i]
Martin Reinecke's avatar
Martin Reinecke committed
596
                                  for i in range(len(self.domain))
Theo Steininger's avatar
Theo Steininger committed
597
                                  if i not in spaces)
598

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

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

604
605
606
607
    def integrate(self, spaces=None):
        tmp = self.weight(1, spaces=spaces)
        return tmp.sum(spaces)

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

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

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

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

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

623
624
    def mean(self, spaces=None):
        return self._contraction_helper('mean', spaces)
csongor's avatar
csongor committed
625

626
627
    def var(self, spaces=None):
        return self._contraction_helper('var', spaces)
csongor's avatar
csongor committed
628

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

632
633
634
635
636
637
638
    def copy_content_from(self, other):
        if not isinstance(other, Field):
            raise TypeError("argument must be a Field")
        if other.domain != self.domain:
            raise ValueError("domains are incompatible.")
        self.val[()] = other.val

Theo Steininger's avatar
Theo Steininger committed
639
    # ---General binary methods---
csongor's avatar
csongor committed
640

641
    def _binary_helper(self, other, op):
csongor's avatar
csongor committed
642
        # if other is a field, make sure that the domains match
643
        if isinstance(other, Field):
644
645
            if other.domain != self.domain:
                raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
Martin Reinecke committed
646
647
            tval = getattr(self.val, op)(other.val)
            return self if tval is self.val else Field(self.domain, tval)
csongor's avatar
csongor committed
648

Martin Reinecke's avatar
Martin Reinecke committed
649
650
        tval = getattr(self.val, op)(other)
        return self if tval is self.val else Field(self.domain, tval)
csongor's avatar
csongor committed
651
652

    def __add__(self, other):
Theo Steininger's avatar
Theo Steininger committed
653
        return self._binary_helper(other, op='__add__')
654

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

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

    def __sub__(self, other):
Theo Steininger's avatar
Theo Steininger committed
662
        return self._binary_helper(other, op='__sub__')
csongor's avatar
csongor committed
663
664

    def __rsub__(self, other):
Theo Steininger's avatar
Theo Steininger committed
665
        return self._binary_helper(other, op='__rsub__')
csongor's avatar
csongor committed
666
667

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

    def __mul__(self, other):
Theo Steininger's avatar
Theo Steininger committed
671
        return self._binary_helper(other, op='__mul__')
672

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

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

    def __div__(self, other):
Theo Steininger's avatar
Theo Steininger committed
680
        return self._binary_helper(other, op='__div__')
csongor's avatar
csongor committed
681

Martin Reinecke's avatar
Martin Reinecke committed
682
683
684
    def __truediv__(self, other):
        return self._binary_helper(other, op='__truediv__')

csongor's avatar
csongor committed
685
    def __rdiv__(self, other):
Theo Steininger's avatar
Theo Steininger committed
686
        return self._binary_helper(other, op='__rdiv__')
csongor's avatar
csongor committed
687

Martin Reinecke's avatar
Martin Reinecke committed
688
689
690
    def __rtruediv__(self, other):
        return self._binary_helper(other, op='__rtruediv__')

csongor's avatar
csongor committed
691
    def __idiv__(self, other):
692
        return self._binary_helper(other, op='__idiv__')
693

csongor's avatar
csongor committed
694
    def __pow__(self, other):
Theo Steininger's avatar
Theo Steininger committed
695
        return self._binary_helper(other, op='__pow__')
csongor's avatar
csongor committed
696
697

    def __rpow__(self, other):
Theo Steininger's avatar
Theo Steininger committed
698
        return self._binary_helper(other, op='__rpow__')
csongor's avatar
csongor committed
699
700

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

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

    def __le__(self, other):
Theo Steininger's avatar
Theo Steininger committed
707
        return self._binary_helper(other, op='__le__')
csongor's avatar
csongor committed
708
709
710
711
712

    def __ne__(self, other):
        if other is None:
            return True
        else:
Theo Steininger's avatar
Theo Steininger committed
713
            return self._binary_helper(other, op='__ne__')
csongor's avatar
csongor committed
714
715
716
717
718

    def __eq__(self, other):
        if other is None:
            return False
        else:
Theo Steininger's avatar
Theo Steininger committed
719
            return self._binary_helper(other, op='__eq__')
csongor's avatar
csongor committed
720
721

    def __ge__(self, other):
Theo Steininger's avatar
Theo Steininger committed
722
        return self._binary_helper(other, op='__ge__')
csongor's avatar
csongor committed
723
724

    def __gt__(self, other):
Theo Steininger's avatar
Theo Steininger committed
725
726
727
        return self._binary_helper(other, op='__gt__')

    def __repr__(self):
Martin Reinecke's avatar
Martin Reinecke committed
728
        return "<nifty2go.Field>"
Theo Steininger's avatar
Theo Steininger committed
729
730
731
732

    def __str__(self):
        minmax = [self.min(), self.max()]
        mean = self.mean()
Martin Reinecke's avatar
Martin Reinecke committed
733
        return "nifty2go.Field instance\n- domain      = " + \
Theo Steininger's avatar
Theo Steininger committed
734
               repr(self.domain) + \
735
               "\n- val         = " + repr(self.val) + \
Theo Steininger's avatar
Theo Steininger committed
736
737
               "\n  - min.,max. = " + str(minmax) + \
               "\n  - mean = " + str(mean)