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

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

29

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

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

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

Martin Reinecke's avatar
Martin Reinecke committed
40
    val : None, Field, data_object, or scalar
41
        The values the array should contain after init. A scalar input will
Martin Reinecke's avatar
Martin Reinecke committed
42
43
        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
44

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

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

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

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

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

81
82
83
84
85
86
87
    # 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")

88
89
    @staticmethod
    def full(domain, val, dtype=None):
Martin Reinecke's avatar
Martin Reinecke committed
90
91
92
93
94
95
96
97
98
99
100
101
102
103
        """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
        """
104
105
106
107
108
109
110
111
        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)

112
    @staticmethod
113
    def from_global_data(domain, arr, sum_up=False):
Martin Reinecke's avatar
Martin Reinecke committed
114
115
116
117
118
119
120
121
122
        """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`.
123
124
125
126
127
        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
128
        """
129
        return Field(domain, dobj.from_global_data(arr, sum_up))
130

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

136
    def to_global_data(self):
Martin Reinecke's avatar
Martin Reinecke committed
137
138
139
140
141
142
143
144
145
146
147
148
        """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.
        """
149
150
        return dobj.to_global_data(self._val)

Martin Reinecke's avatar
Martin Reinecke committed
151
152
    @property
    def local_data(self):
Martin Reinecke's avatar
Martin Reinecke committed
153
154
155
156
157
158
159
160
161
162
        """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
163
164
        return dobj.local_data(self._val)

Martin Reinecke's avatar
Martin Reinecke committed
165
    def cast_domain(self, new_domain):
Martin Reinecke's avatar
Martin Reinecke committed
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
        """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
183
184
        return Field(new_domain, self._val)

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

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

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

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

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

Martin Reinecke's avatar
Martin Reinecke committed
228
    def fill(self, fill_value):
Martin Reinecke's avatar
Martin Reinecke committed
229
230
231
232
233
234
235
        """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
236
        self._val.fill(fill_value)
237
        return self
Martin Reinecke's avatar
Martin Reinecke committed
238

Martin Reinecke's avatar
Martin Reinecke committed
239
    def lock(self):
Martin Reinecke's avatar
Martin Reinecke committed
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
        """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
255
        dobj.lock(self._val)
256
        return self
Martin Reinecke's avatar
Martin Reinecke committed
257
258
259

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

Theo Steininger's avatar
Theo Steininger committed
263
264
    @property
    def val(self):
Martin Reinecke's avatar
Martin Reinecke committed
265
266
267
268
269
270
        """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.
271
        """
Martin Reinecke's avatar
Martin Reinecke committed
272
        return self._val
csongor's avatar
csongor committed
273

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

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

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

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

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

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

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

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

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

321
322
323
324
325
326
327
328
329
330
331
    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
332
333
334
335
336
337
338
339
340
    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
341
            A read-only version of `self`.
Martin Reinecke's avatar
Martin Reinecke committed
342
        """
Martin Reinecke's avatar
Martin Reinecke committed
343
        return self if self.locked else Field(val=self, copy=True, locked=True)
Martin Reinecke's avatar
Martin Reinecke committed
344

345
    def scalar_weight(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
346
347
348
349
350
351
352
353
354
355
356
357
358
359
        """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.
        """
360
        if np.isscalar(spaces):
Martin Reinecke's avatar
Martin Reinecke committed
361
            return self._domain[spaces].scalar_dvol
362
363

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

Martin Reinecke's avatar
Martin Reinecke committed
373
    def total_volume(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
374
375
376
377
378
379
380
381
382
383
384
385
386
        """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
387
        if np.isscalar(spaces):
Martin Reinecke's avatar
Martin Reinecke committed
388
            return self._domain[spaces].total_volume
Martin Reinecke's avatar
Martin Reinecke committed
389
390
391
392
393

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

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

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

Martin Reinecke's avatar
Martin Reinecke committed
405
406
407
        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
408

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

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

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

427
428
        fct = 1.
        for ind in spaces:
Martin Reinecke's avatar
Martin Reinecke committed
429
            wgt = self._domain[ind].dvol
430
431
432
            if np.isscalar(wgt):
                fct *= wgt
            else:
Martin Reinecke's avatar
Martin Reinecke committed
433
                new_shape = np.ones(len(self.shape), dtype=np.int)
Martin Reinecke's avatar
Martin Reinecke committed
434
435
                new_shape[self._domain.axes[ind][0]:
                          self._domain.axes[ind][-1]+1] = wgt.shape
436
                wgt = wgt.reshape(new_shape)
Martin Reinecke's avatar
Martin Reinecke committed
437
438
                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
439
                    wgt = dobj.local_data(dobj.from_global_data(wgt))
Martin Reinecke's avatar
Martin Reinecke committed
440
                out.local_data[()] *= wgt**power
441
        fct = fct**power
Martin Reinecke's avatar
Martin Reinecke committed
442
        if fct != 1.:
443
            out *= fct
444

445
        return out
csongor's avatar
csongor committed
446

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

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

Martin Reinecke's avatar
Martin Reinecke committed
455
        spaces : None, int or tuple of int (default: None)
456
457
            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
458

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

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

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

        if len(spaces) == ndom:
Martin Reinecke's avatar
Martin Reinecke committed
475
            return dobj.vdot(self.val, x.val)
476
477
        # 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
478
        return (self.conjugate()*x).sum(spaces=spaces)
Theo Steininger's avatar
Theo Steininger committed
479

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

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

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

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

Theo Steininger's avatar
Theo Steininger committed
500
    # ---General unary/contraction methods---
501

Theo Steininger's avatar
Theo Steininger committed
502
503
    def __pos__(self):
        return self.copy()
504

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

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

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

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

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

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

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

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

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

535
    def sum(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
536
537
538
539
540
541
542
543
544
545
546
547
548
549
        """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.
        """
550
        return self._contraction_helper('sum', spaces)
csongor's avatar
csongor committed
551

552
    def integrate(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
        """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
570
571
572
573
574
        swgt = self.scalar_weight(spaces)
        if swgt is not None:
            res = self.sum(spaces)
            res *= swgt
            return res
575
576
577
        tmp = self.weight(1, spaces=spaces)
        return tmp.sum(spaces)

578
    def prod(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
579
580
581
582
583
584
585
586
587
588
589
590
591
592
        """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.
        """
593
        return self._contraction_helper('prod', spaces)
csongor's avatar
csongor committed
594

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

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

601
    def min(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
602
603
604
605
606
607
608
609
610
611
612
613
614
615
        """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.
        """
616
        return self._contraction_helper('min', spaces)
csongor's avatar
csongor committed
617

618
    def max(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
619
620
621
622
623
624
625
626
627
628
629
630
631
632
        """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.
        """
633
        return self._contraction_helper('max', spaces)
csongor's avatar
csongor committed
634

635
    def mean(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
        """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.
        """
653
654
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('mean', spaces)
Martin Reinecke's avatar
Martin Reinecke committed
655
        # MR FIXME: not very efficient
656
657
        # MR FIXME: do we need "spaces" here?
        tmp = self.weight(1, spaces)
Martin Reinecke's avatar
Martin Reinecke committed
658
        return tmp.sum(spaces)*(1./tmp.total_volume(spaces))
csongor's avatar
csongor committed
659

660
    def var(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
661
662
663
664
665
666
667
668
669
670
671
672
673
674
        """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.
        """
675
676
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('var', spaces)
Martin Reinecke's avatar
Martin Reinecke committed
677
678
679
        # MR FIXME: not very efficient or accurate
        m1 = self.mean(spaces)
        if np.issubdtype(self.dtype, np.complexfloating):
680
            sq = abs(self-m1)**2
Martin Reinecke's avatar
Martin Reinecke committed
681
        else:
682
683
            sq = (self-m1)**2
        return sq.mean(spaces)
csongor's avatar
csongor committed
684

685
    def std(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
        """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.
        """
703
        from .sugar import sqrt
704
705
706
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('std', spaces)
        return sqrt(self.var(spaces))
csongor's avatar
csongor committed
707

708
709
710
    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
711
        if other._domain != self._domain:
712
            raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
Martin Reinecke committed
713
        self.local_data[()] = other.local_data[()]
714

Theo Steininger's avatar
Theo Steininger committed
715
    def __repr__(self):
Philipp Arras's avatar
Philipp Arras committed
716
        return "<nifty5.Field>"
Theo Steininger's avatar
Theo Steininger committed
717
718

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

723
724
725
    def isEquivalentTo(self, other):
        """Determines (as quickly as possible) whether `self`'s content is
        identical to `other`'s content."""
726
727
728
729
730
731
732
733
        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()

734
735
736
737
738
    def isSubsetOf(self, other):
        """Identical to `Field.isEquivalentTo()`. This method is provided for
        easier interoperability with `MultiField`."""
        return self.isEquivalentTo(other)

739

740
741
742
743
744
745
746
747
748
749
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
750
751
752
753
754
755
756
757
758
759
760
761
            # 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
762
763
        return func2
    setattr(Field, op, func(op))