field.py 24.8 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
#
Martin Reinecke's avatar
Martin Reinecke committed
14
# Copyright(C) 2013-2018 Max-Planck-Society
Theo Steininger's avatar
Theo Steininger committed
15
16
17
#
# 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
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
from . import utilities
Martin Reinecke's avatar
Martin Reinecke committed
23
from .domain_tuple import DomainTuple
Martin Reinecke's avatar
Martin Reinecke committed
24
from functools import reduce
25
from . import dobj
26
import sys
27

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
28
29
__all__ = ["Field", "sqrt", "exp", "log", "conjugate"]

30

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

Martin Reinecke's avatar
Martin Reinecke committed
34
    In NIFTy, Fields are used to store data arrays and carry all the needed
35
    metainformation (i.e. the domain) for operators to be able to work on them.
Theo Steininger's avatar
Theo Steininger committed
36

37
38
    Parameters
    ----------
Martin Reinecke's avatar
Martin Reinecke committed
39
    domain : None, DomainTuple, tuple of Domain, or Domain
Theo Steininger's avatar
Theo Steininger committed
40

Martin Reinecke's avatar
Martin Reinecke committed
41
    val : None, Field, data_object, or scalar
42
        The values the array should contain after init. A scalar input will
Martin Reinecke's avatar
Martin Reinecke committed
43
44
        fill the whole array with this scalar. If a data_object is provided,
        its dimensions must match the domain's.
Theo Steininger's avatar
Theo Steininger committed
45

46
    dtype : type
Martin Reinecke's avatar
updates    
Martin Reinecke committed
47
        A numpy.type. Most common are float and complex.
Martin Reinecke's avatar
Martin Reinecke committed
48
49
50
51
52

    Notes
    -----
    If possible, do not invoke the constructor directly, but use one of the
    many convenience functions for Field conatruction!
53
    """
54

Martin Reinecke's avatar
Martin Reinecke committed
55
56
    def __init__(self, domain=None, val=None, dtype=None, copy=False,
                 locked=False):
Martin Reinecke's avatar
Martin Reinecke committed
57
        self._domain = self._infer_domain(domain=domain, val=val)
58

Martin Reinecke's avatar
Martin Reinecke committed
59
        dtype = self._infer_dtype(dtype=dtype, val=val)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
60
        if isinstance(val, Field):
Martin Reinecke's avatar
Martin Reinecke committed
61
            if self._domain != val._domain:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
62
                raise ValueError("Domain mismatch")
Martin Reinecke's avatar
Martin Reinecke committed
63
64
            self._val = dobj.from_object(val.val, dtype=dtype, copy=copy,
                                         set_locked=locked)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
65
        elif (np.isscalar(val)):
Martin Reinecke's avatar
Martin Reinecke committed
66
            self._val = dobj.full(self._domain.shape, dtype=dtype,
67
                                  fill_value=val)
68
        elif isinstance(val, dobj.data_object):
Martin Reinecke's avatar
Martin Reinecke committed
69
            if self._domain.shape == val.shape:
Martin Reinecke's avatar
Martin Reinecke committed
70
71
                self._val = dobj.from_object(val, dtype=dtype, copy=copy,
                                             set_locked=locked)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
72
73
74
            else:
                raise ValueError("Shape mismatch")
        elif val is None:
Martin Reinecke's avatar
Martin Reinecke committed
75
            self._val = dobj.empty(self._domain.shape, dtype=dtype)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
76
77
        else:
            raise TypeError("unknown source type")
csongor's avatar
csongor committed
78

Martin Reinecke's avatar
Martin Reinecke committed
79
80
81
        if locked:
            dobj.lock(self._val)

82
83
84
85
86
87
88
    # prevent implicit conversion to bool
    def __nonzero__(self):
        raise TypeError("Field does not support implicit conversion to bool")

    def __bool__(self):
        raise TypeError("Field does not support implicit conversion to bool")

89
90
    @staticmethod
    def full(domain, val, dtype=None):
Martin Reinecke's avatar
Martin Reinecke committed
91
92
93
94
95
96
97
98
99
100
101
102
103
104
        """Creates a Field with a given domain, filled with a constant value.

        Parameters
        ----------
        domain : Domain, tuple of Domain, or DomainTuple
            domain of the new Field
        val : float/complex/int scalar
            fill value. Data type of the field is inferred from val.

        Returns
        -------
        Field
            the newly created field
        """
105
106
107
108
109
110
111
112
        if not np.isscalar(val):
            raise TypeError("val must be a scalar")
        return Field(DomainTuple.make(domain), val, dtype)

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

113
    @staticmethod
114
    def from_global_data(domain, arr, sum_up=False):
Martin Reinecke's avatar
Martin Reinecke committed
115
116
117
118
119
120
121
122
123
        """Returns a Field constructed from `domain` and `arr`.

        Parameters
        ----------
        domain : DomainTuple, tuple of Domain, or Domain
            the domain of the new Field
        arr : numpy.ndarray
            The data content to be used for the new Field.
            Its shape must match the shape of `domain`.
124
125
126
127
128
        sum_up : bool, optional
            If True, the contents of `arr` are summed up over all MPI tasks
            (if any), and the sum is used as data content.
            If False, the contens of `arr` are used directly, and must be
            identical on all MPI tasks.
Martin Reinecke's avatar
Martin Reinecke committed
129
        """
130
        return Field(domain, dobj.from_global_data(arr, sum_up))
131

Martin Reinecke's avatar
Martin Reinecke committed
132
133
134
135
136
    @staticmethod
    def from_local_data(domain, arr):
        domain = DomainTuple.make(domain)
        return Field(domain, dobj.from_local_data(domain.shape, arr))

137
    def to_global_data(self):
Martin Reinecke's avatar
Martin Reinecke committed
138
139
140
141
142
143
144
145
146
147
148
149
        """Returns an array containing the full data of the field.

        Returns
        -------
        numpy.ndarray : array containing all field entries.
            Its shape is identical to `self.shape`.

        Notes
        -----
        Do not write to the returned array! Depending on whether MPI is
        active or not, this may or may not change the field's data content.
        """
150
151
        return dobj.to_global_data(self._val)

Martin Reinecke's avatar
Martin Reinecke committed
152
153
    @property
    def local_data(self):
Martin Reinecke's avatar
Martin Reinecke committed
154
155
156
157
158
159
160
161
162
163
        """numpy.ndarray : locally residing field data

        Returns a handle to the part of the array data residing on the local
        task (or to the entore array if MPI is not active).

        Notes
        -----
        If the field is not locked, the array data can be modified.
        Use with care!
        """
Martin Reinecke's avatar
Martin Reinecke committed
164
165
        return dobj.local_data(self._val)

Martin Reinecke's avatar
Martin Reinecke committed
166
    def cast_domain(self, new_domain):
Martin Reinecke's avatar
Martin Reinecke committed
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
        """Returns a field with the same data, but a different domain

        Parameters
        ----------
        new_domain : Domain, tuple of Domain, or DomainTuple
            The domain for the returned field. Must be shape-compatible to
            `self`.

        Returns
        -------
        Field
            Field living on `new_domain`, but with the same data as `self`.

        Notes
        -----
        No copy is made. If needed, use an additional copy() invocation.
        """
Martin Reinecke's avatar
Martin Reinecke committed
184
185
        return Field(new_domain, self._val)

Martin Reinecke's avatar
Martin Reinecke committed
186
    @staticmethod
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
187
    def _infer_domain(domain, val=None):
188
        if domain is None:
189
            if isinstance(val, Field):
Martin Reinecke's avatar
Martin Reinecke committed
190
                return val._domain
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
191
            if np.isscalar(val):
192
                return DomainTuple.make(())  # empty domain tuple
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
193
            raise TypeError("could not infer domain from value")
Martin Reinecke's avatar
Martin Reinecke committed
194
        return DomainTuple.make(domain)
195

Martin Reinecke's avatar
Martin Reinecke committed
196
197
    @staticmethod
    def _infer_dtype(dtype, val):
Martin Reinecke's avatar
Martin Reinecke committed
198
199
200
201
        if dtype is not None:
            return dtype
        if val is None:
            raise ValueError("could not infer dtype")
Martin Reinecke's avatar
Martin Reinecke committed
202
203
        if isinstance(val, Field):
            return val.dtype
Martin Reinecke's avatar
Martin Reinecke committed
204
        return np.result_type(val)
205

Martin Reinecke's avatar
Martin Reinecke committed
206
207
    @staticmethod
    def from_random(random_type, domain, dtype=np.float64, **kwargs):
208
209
210
211
        """ Draws a random field with the given parameters.

        Parameters
        ----------
Martin Reinecke's avatar
Martin Reinecke committed
212
213
        random_type : 'pm1', 'normal', or 'uniform'
            The random distribution to use.
Martin Reinecke's avatar
Martin Reinecke committed
214
        domain : DomainTuple
215
216
217
            The domain of the output random field
        dtype : type
            The datatype of the output random field
Theo Steininger's avatar
Theo Steininger committed
218

219
220
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
221
        Field
Martin Reinecke's avatar
Martin Reinecke committed
222
            The newly created Field.
223
        """
Martin Reinecke's avatar
Martin Reinecke committed
224
        domain = DomainTuple.make(domain)
Martin Reinecke's avatar
Martin Reinecke committed
225
226
227
        return Field(domain=domain,
                     val=dobj.from_random(random_type, dtype=dtype,
                                          shape=domain.shape, **kwargs))
228

Martin Reinecke's avatar
Martin Reinecke committed
229
    def fill(self, fill_value):
Martin Reinecke's avatar
Martin Reinecke committed
230
231
232
233
234
235
236
        """Fill `self` uniformly with `fill_value`

        Parameters
        ----------
        fill_value: float or complex or int
            The value to fill the field with.
        """
Martin Reinecke's avatar
Martin Reinecke committed
237
        self._val.fill(fill_value)
238
        return self
Martin Reinecke's avatar
Martin Reinecke committed
239

Martin Reinecke's avatar
Martin Reinecke committed
240
    def lock(self):
Martin Reinecke's avatar
Martin Reinecke committed
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
        """Write-protect the data content of `self`.

        After this call, it will no longer be possible to change the data
        entries of `self`. This is convenient if, for example, a
        DiagonalOperator wants to ensure that its diagonal cannot be modified
        inadvertently, without making copies.

        Notes
        -----
        This will not only prohibit modifications to the entries of `self`, but
        also to the entries of any other Field or numpy array pointing to the
        same data. If an unlocked instance is needed, use copy().

        The fact that there is no `unlock()` method is deliberate.
        """
Martin Reinecke's avatar
Martin Reinecke committed
256
        dobj.lock(self._val)
257
        return self
Martin Reinecke's avatar
Martin Reinecke committed
258
259
260

    @property
    def locked(self):
Martin Reinecke's avatar
Martin Reinecke committed
261
        """bool : True iff the field's data content has been locked"""
Martin Reinecke's avatar
Martin Reinecke committed
262
263
        return dobj.locked(self._val)

theos's avatar
theos committed
264
265
    @property
    def val(self):
Martin Reinecke's avatar
Martin Reinecke committed
266
267
268
269
270
271
        """dobj.data_object : the data object storing the field's entries

        Notes
        -----
        This property is intended for low-level, internal use only. Do not use
        from outside of NIFTy's core; there should be better alternatives.
272
        """
Martin Reinecke's avatar
Martin Reinecke committed
273
        return self._val
csongor's avatar
csongor committed
274

Martin Reinecke's avatar
Martin Reinecke committed
275
276
    @property
    def dtype(self):
Martin Reinecke's avatar
Martin Reinecke committed
277
        """type : the data type of the field's entries"""
Martin Reinecke's avatar
Martin Reinecke committed
278
279
        return self._val.dtype

Martin Reinecke's avatar
Martin Reinecke committed
280
281
    @property
    def domain(self):
Martin Reinecke's avatar
Martin Reinecke committed
282
        """DomainTuple : the field's domain"""
Martin Reinecke's avatar
Martin Reinecke committed
283
284
        return self._domain

285
286
    @property
    def shape(self):
Martin Reinecke's avatar
Martin Reinecke committed
287
        """tuple of int : the concatenated shapes of all sub-domains"""
Martin Reinecke's avatar
Martin Reinecke committed
288
        return self._domain.shape
csongor's avatar
csongor committed
289

290
    @property
Martin Reinecke's avatar
Martin Reinecke committed
291
    def size(self):
Martin Reinecke's avatar
Martin Reinecke committed
292
        """int : total number of pixels in the field"""
Martin Reinecke's avatar
Martin Reinecke committed
293
        return self._domain.size
csongor's avatar
csongor committed
294

Theo Steininger's avatar
Theo Steininger committed
295
296
    @property
    def real(self):
Martin Reinecke's avatar
Martin Reinecke committed
297
        """Field : The real part of the field"""
298
        if not np.issubdtype(self.dtype, np.complexfloating):
299
            return self
Martin Reinecke's avatar
Martin Reinecke committed
300
        return Field(self._domain, self.val.real)
Theo Steininger's avatar
Theo Steininger committed
301
302
303

    @property
    def imag(self):
Martin Reinecke's avatar
Martin Reinecke committed
304
        """Field : The imaginary part of the field"""
305
306
        if not np.issubdtype(self.dtype, np.complexfloating):
            raise ValueError(".imag called on a non-complex Field")
Martin Reinecke's avatar
Martin Reinecke committed
307
        return Field(self._domain, self.val.imag)
Theo Steininger's avatar
Theo Steininger committed
308

Martin Reinecke's avatar
Martin Reinecke committed
309
    def copy(self):
310
        """ Returns a full copy of the Field.
Theo Steininger's avatar
Theo Steininger committed
311

Martin Reinecke's avatar
Martin Reinecke committed
312
        The returned object will be an identical copy of the original Field.
Martin Reinecke's avatar
Martin Reinecke committed
313
        The copy will be writeable, even if `self` was locked.
Theo Steininger's avatar
Theo Steininger committed
314

315
316
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
317
        Field
Martin Reinecke's avatar
Martin Reinecke committed
318
            An identical, but unlocked copy of 'self'.
319
        """
Martin Reinecke's avatar
Martin Reinecke committed
320
        return Field(val=self, copy=True)
csongor's avatar
csongor committed
321

322
323
324
325
326
327
328
329
330
331
332
    def empty_copy(self):
        """ Returns a Field with identical domain and data type, but
        uninitialized data.

        Returns
        -------
        Field
            A copy of 'self', with uninitialized data.
        """
        return Field(self._domain, dtype=self.dtype)

Martin Reinecke's avatar
Martin Reinecke committed
333
334
335
336
337
338
339
340
341
    def locked_copy(self):
        """ Returns a read-only version of the Field.

        If `self` is locked, returns `self`. Otherwise returns a locked copy
        of `self`.

        Returns
        -------
        Field
Martin Reinecke's avatar
Martin Reinecke committed
342
            A read-only version of `self`.
Martin Reinecke's avatar
Martin Reinecke committed
343
        """
Martin Reinecke's avatar
Martin Reinecke committed
344
        return self if self.locked else Field(val=self, copy=True, locked=True)
Martin Reinecke's avatar
Martin Reinecke committed
345

346
    def scalar_weight(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
347
348
349
350
351
352
353
354
355
356
357
358
359
360
        """Returns the uniform volume element for a sub-domain of `self`.

        Parameters
        ----------
        spaces : int, tuple of int or None
            indices of the sub-domains of the field's domain to be considered.
            If `None`, the entire domain is used.

        Returns
        -------
        float or None
            if the requested sub-domain has a uniform volume element, it is
            returned. Otherwise, `None` is returned.
        """
361
        if np.isscalar(spaces):
Martin Reinecke's avatar
Martin Reinecke committed
362
            return self._domain[spaces].scalar_dvol
363
364

        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
365
            spaces = range(len(self._domain))
Martin Reinecke's avatar
Martin Reinecke committed
366
        res = 1.
367
        for i in spaces:
Martin Reinecke's avatar
Martin Reinecke committed
368
            tmp = self._domain[i].scalar_dvol
369
370
371
372
373
            if tmp is None:
                return None
            res *= tmp
        return res

Martin Reinecke's avatar
Martin Reinecke committed
374
    def total_volume(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
375
376
377
378
379
380
381
382
383
384
385
386
387
        """Returns the total volume of a sub-domain of `self`.

        Parameters
        ----------
        spaces : int, tuple of int or None
            indices of the sub-domains of the field's domain to be considered.
            If `None`, the entire domain is used.

        Returns
        -------
        float
            the total volume of the requested sub-domain.
        """
Martin Reinecke's avatar
Martin Reinecke committed
388
        if np.isscalar(spaces):
Martin Reinecke's avatar
Martin Reinecke committed
389
            return self._domain[spaces].total_volume
Martin Reinecke's avatar
Martin Reinecke committed
390
391
392
393
394

        if spaces is None:
            spaces = range(len(self._domain))
        res = 1.
        for i in spaces:
Martin Reinecke's avatar
Martin Reinecke committed
395
            res *= self._domain[i].total_volume
Martin Reinecke's avatar
Martin Reinecke committed
396
397
        return res

398
    def weight(self, power=1, spaces=None, out=None):
Theo Steininger's avatar
Theo Steininger committed
399
        """ Weights the pixels of `self` with their invidual pixel-volume.
400
401
402
403

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

Martin Reinecke's avatar
Martin Reinecke committed
406
407
408
        spaces : None, int or tuple of int
            Determines on which sub-domain the operation takes place.
            If None, the entire domain is used.
Theo Steininger's avatar
Theo Steininger committed
409

410
411
412
413
414
        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"!

415
416
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
417
        Field
Theo Steininger's avatar
Theo Steininger committed
418
            The weighted field.
419
        """
420
421
422
423
424
        if out is None:
            out = self.copy()
        else:
            if out is not self:
                out.copy_content_from(self)
csongor's avatar
csongor committed
425

Martin Reinecke's avatar
Martin Reinecke committed
426
        spaces = utilities.parse_spaces(spaces, len(self._domain))
csongor's avatar
csongor committed
427

428
429
        fct = 1.
        for ind in spaces:
Martin Reinecke's avatar
Martin Reinecke committed
430
            wgt = self._domain[ind].dvol
431
432
433
            if np.isscalar(wgt):
                fct *= wgt
            else:
Martin Reinecke's avatar
Martin Reinecke committed
434
                new_shape = np.ones(len(self.shape), dtype=np.int)
Martin Reinecke's avatar
Martin Reinecke committed
435
436
                new_shape[self._domain.axes[ind][0]:
                          self._domain.axes[ind][-1]+1] = wgt.shape
437
                wgt = wgt.reshape(new_shape)
Martin Reinecke's avatar
Martin Reinecke committed
438
439
                if dobj.distaxis(self._val) >= 0 and ind == 0:
                    # we need to distribute the weights along axis 0
Martin Reinecke's avatar
fixes    
Martin Reinecke committed
440
                    wgt = dobj.local_data(dobj.from_global_data(wgt))
Martin Reinecke's avatar
Martin Reinecke committed
441
                out.local_data[()] *= wgt**power
442
        fct = fct**power
Martin Reinecke's avatar
Martin Reinecke committed
443
        if fct != 1.:
444
            out *= fct
445

446
        return out
csongor's avatar
csongor committed
447

Martin Reinecke's avatar
Martin Reinecke committed
448
    def vdot(self, x=None, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
449
        """ Computes the dot product of 'self' with x.
Theo Steininger's avatar
Theo Steininger committed
450

451
452
453
        Parameters
        ----------
        x : Field
454
            x must live on the same domain as `self`.
Theo Steininger's avatar
Theo Steininger committed
455

Martin Reinecke's avatar
Martin Reinecke committed
456
        spaces : None, int or tuple of int (default: None)
457
458
            The dot product is only carried out over the sub-domains in this
            tuple. If None, it is carried out over all sub-domains.
Theo Steininger's avatar
Theo Steininger committed
459

460
461
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
462
        float, complex, either scalar (for full dot products)
463
                              or Field (for partial dot products)
464
        """
465
        if not isinstance(x, Field):
466
467
            raise TypeError("The dot-partner must be an instance of " +
                            "the NIFTy field class")
theos's avatar
theos committed
468

Martin Reinecke's avatar
Martin Reinecke committed
469
        if x._domain != self._domain:
470
            raise ValueError("Domain mismatch")
theos's avatar
theos committed
471

Martin Reinecke's avatar
Martin Reinecke committed
472
        ndom = len(self._domain)
473
474
475
        spaces = utilities.parse_spaces(spaces, ndom)

        if len(spaces) == ndom:
Martin Reinecke's avatar
Martin Reinecke committed
476
            return dobj.vdot(self.val, x.val)
477
478
        # If we arrive here, we have to do a partial dot product.
        # For the moment, do this the explicit, non-optimized way
Martin Reinecke's avatar
Martin Reinecke committed
479
        return (self.conjugate()*x).sum(spaces=spaces)
theos's avatar
theos committed
480

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

Theo Steininger's avatar
Theo Steininger committed
484
485
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
486
        float
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
487
            The L2-norm of the field values.
csongor's avatar
csongor committed
488
        """
489
        return np.sqrt(abs(self.vdot(x=self)))
csongor's avatar
csongor committed
490

Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
491
    def conjugate(self):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
492
        """ Returns the complex conjugate of the field.
Theo Steininger's avatar
Theo Steininger committed
493

494
495
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
496
497
        Field
            The complex conjugated field.
csongor's avatar
csongor committed
498
        """
Martin Reinecke's avatar
Martin Reinecke committed
499
        return Field(self._domain, self.val.conjugate())
csongor's avatar
csongor committed
500

theos's avatar
theos committed
501
    # ---General unary/contraction methods---
502

theos's avatar
theos committed
503
504
    def __pos__(self):
        return self.copy()
505

theos's avatar
theos committed
506
    def __neg__(self):
507
        return Field(self._domain, -self.val)
csongor's avatar
csongor committed
508

theos's avatar
theos committed
509
    def __abs__(self):
510
        return Field(self._domain, abs(self.val))
csongor's avatar
csongor committed
511

512
    def _contraction_helper(self, op, spaces):
theos's avatar
theos committed
513
        if spaces is None:
514
            return getattr(self.val, op)()
515

Martin Reinecke's avatar
Martin Reinecke committed
516
        spaces = utilities.parse_spaces(spaces, len(self._domain))
csongor's avatar
csongor committed
517

Martin Reinecke's avatar
Martin Reinecke committed
518
        axes_list = tuple(self._domain.axes[sp_index] for sp_index in spaces)
519

Martin Reinecke's avatar
Martin Reinecke committed
520
        if len(axes_list) > 0:
theos's avatar
theos committed
521
            axes_list = reduce(lambda x, y: x+y, axes_list)
csongor's avatar
csongor committed
522

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

theos's avatar
theos committed
526
527
528
        # 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
529
        else:
Martin Reinecke's avatar
Martin Reinecke committed
530
            return_domain = tuple(dom
Martin Reinecke's avatar
Martin Reinecke committed
531
                                  for i, dom in enumerate(self._domain)
theos's avatar
theos committed
532
                                  if i not in spaces)
533

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

536
    def sum(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
537
538
539
540
541
542
543
544
545
546
547
548
549
550
        """Sums up over the sub-domains given by `spaces`.

        Parameters
        ----------
        spaces : None, int or tuple of int (default: None)
            The summation is only carried out over the sub-domains in this
            tuple. If None, it is carried out over all sub-domains.

        Returns
        -------
        Field or scalar
            The result of the summation. If it is carried out over the entire
            domain, this is a scalar, otherwise a Field.
        """
551
        return self._contraction_helper('sum', spaces)
csongor's avatar
csongor committed
552

553
    def integrate(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
        """Integrates over the sub-domains given by `spaces`.

        Integration is performed by summing over `self` multiplied by its
        volume factors.

        Parameters
        ----------
        spaces : None, int or tuple of int (default: None)
            The summation is only carried out over the sub-domains in this
            tuple. If None, it is carried out over all sub-domains.

        Returns
        -------
        Field or scalar
            The result of the integration. If it is carried out over the
            entire domain, this is a scalar, otherwise a Field.
        """
Martin Reinecke's avatar
Martin Reinecke committed
571
572
573
574
575
        swgt = self.scalar_weight(spaces)
        if swgt is not None:
            res = self.sum(spaces)
            res *= swgt
            return res
576
577
578
        tmp = self.weight(1, spaces=spaces)
        return tmp.sum(spaces)

579
    def prod(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
580
581
582
583
584
585
586
587
588
589
590
591
592
593
        """Computes the product over the sub-domains given by `spaces`.

        Parameters
        ----------
        spaces : None, int or tuple of int (default: None)
            The operation is only carried out over the sub-domains in this
            tuple. If None, it is carried out over all sub-domains.

        Returns
        -------
        Field or scalar
            The result of the product. If it is carried out over the entire
            domain, this is a scalar, otherwise a Field.
        """
594
        return self._contraction_helper('prod', spaces)
csongor's avatar
csongor committed
595

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

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

602
    def min(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
603
604
605
606
607
608
609
610
611
612
613
614
615
616
        """Determines the minimum over the sub-domains given by `spaces`.

        Parameters
        ----------
        spaces : None, int or tuple of int (default: None)
            The operation is only carried out over the sub-domains in this
            tuple. If None, it is carried out over all sub-domains.

        Returns
        -------
        Field or scalar
            The result of the operation. If it is carried out over the entire
            domain, this is a scalar, otherwise a Field.
        """
617
        return self._contraction_helper('min', spaces)
csongor's avatar
csongor committed
618

619
    def max(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
620
621
622
623
624
625
626
627
628
629
630
631
632
633
        """Determines the maximum over the sub-domains given by `spaces`.

        Parameters
        ----------
        spaces : None, int or tuple of int (default: None)
            The operation is only carried out over the sub-domains in this
            tuple. If None, it is carried out over all sub-domains.

        Returns
        -------
        Field or scalar
            The result of the operation. If it is carried out over the entire
            domain, this is a scalar, otherwise a Field.
        """
634
        return self._contraction_helper('max', spaces)
csongor's avatar
csongor committed
635

636
    def mean(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
        """Determines the mean over the sub-domains given by `spaces`.

        ``x.mean(spaces)`` is equivalent to
        ``x.integrate(spaces)/x.total_volume(spaces)``.

        Parameters
        ----------
        spaces : None, int or tuple of int (default: None)
            The operation is only carried out over the sub-domains in this
            tuple. If None, it is carried out over all sub-domains.

        Returns
        -------
        Field or scalar
            The result of the operation. If it is carried out over the entire
            domain, this is a scalar, otherwise a Field.
        """
654
655
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('mean', spaces)
Martin Reinecke's avatar
Martin Reinecke committed
656
        # MR FIXME: not very efficient
657
658
        # MR FIXME: do we need "spaces" here?
        tmp = self.weight(1, spaces)
Martin Reinecke's avatar
Martin Reinecke committed
659
        return tmp.sum(spaces)*(1./tmp.total_volume(spaces))
csongor's avatar
csongor committed
660

661
    def var(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
662
663
664
665
666
667
668
669
670
671
672
673
674
675
        """Determines the variance over the sub-domains given by `spaces`.

        Parameters
        ----------
        spaces : None, int or tuple of int (default: None)
            The operation is only carried out over the sub-domains in this
            tuple. If None, it is carried out over all sub-domains.

        Returns
        -------
        Field or scalar
            The result of the operation. If it is carried out over the entire
            domain, this is a scalar, otherwise a Field.
        """
676
677
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('var', spaces)
Martin Reinecke's avatar
Martin Reinecke committed
678
679
680
        # MR FIXME: not very efficient or accurate
        m1 = self.mean(spaces)
        if np.issubdtype(self.dtype, np.complexfloating):
681
            sq = abs(self-m1)**2
Martin Reinecke's avatar
Martin Reinecke committed
682
        else:
683
684
            sq = (self-m1)**2
        return sq.mean(spaces)
csongor's avatar
csongor committed
685

686
    def std(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
        """Determines the standard deviation over the sub-domains given by
        `spaces`.

        ``x.std(spaces)`` is equivalent to ``sqrt(x.var(spaces))``.

        Parameters
        ----------
        spaces : None, int or tuple of int (default: None)
            The operation is only carried out over the sub-domains in this
            tuple. If None, it is carried out over all sub-domains.

        Returns
        -------
        Field or scalar
            The result of the operation. If it is carried out over the entire
            domain, this is a scalar, otherwise a Field.
        """
704
        from .sugar import sqrt
705
706
707
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('std', spaces)
        return sqrt(self.var(spaces))
csongor's avatar
csongor committed
708

709
710
711
    def copy_content_from(self, other):
        if not isinstance(other, Field):
            raise TypeError("argument must be a Field")
Martin Reinecke's avatar
Martin Reinecke committed
712
        if other._domain != self._domain:
713
            raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
Martin Reinecke committed
714
        self.local_data[()] = other.local_data[()]
715

theos's avatar
theos committed
716
    def __repr__(self):
Martin Reinecke's avatar
Martin Reinecke committed
717
        return "<nifty4.Field>"
theos's avatar
theos committed
718
719

    def __str__(self):
Martin Reinecke's avatar
Martin Reinecke committed
720
        return "nifty4.Field instance\n- domain      = " + \
721
               self._domain.__str__() + \
722
               "\n- val         = " + repr(self.val)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
723

724
725
726
727
728
729
730
731
732
    def equivalent(self, other):
        if self is other:
            return True
        if not isinstance(other, Field):
            return False
        if self._domain != other._domain:
            return False
        return (self._val == other._val).all()

733
734
735
736
737
738
739
740
741
742
for op in ["__add__", "__radd__", "__iadd__",
           "__sub__", "__rsub__", "__isub__",
           "__mul__", "__rmul__", "__imul__",
           "__div__", "__rdiv__", "__idiv__",
           "__truediv__", "__rtruediv__", "__itruediv__",
           "__floordiv__", "__rfloordiv__", "__ifloordiv__",
           "__pow__", "__rpow__", "__ipow__",
           "__lt__", "__le__", "__gt__", "__ge__", "__eq__", "__ne__"]:
    def func(op):
        def func2(self, other):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
743
744
745
746
747
748
749
750
751
752
753
754
            # if other is a field, make sure that the domains match
            if isinstance(other, Field):
                if other._domain != self._domain:
                    raise ValueError("domains are incompatible.")
                tval = getattr(self.val, op)(other.val)
                return self if tval is self.val else Field(self._domain, tval)

            if np.isscalar(other) or isinstance(other, dobj.data_object):
                tval = getattr(self.val, op)(other)
                return self if tval is self.val else Field(self._domain, tval)

            return NotImplemented
755
756
        return func2
    setattr(Field, op, func(op))