field.py 25.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
#
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

27

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

407
408
409
410
411
        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"!

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

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

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

443
        return out
csongor's avatar
csongor committed
444

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

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

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

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

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

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

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

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

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

488
489
490
491
492
493
494
495
496
497
    def squared_norm(self):
        """ Computes the square of the L2-norm of the field values.

        Returns
        -------
        float
            The square of the L2-norm of the field values.
        """
        return abs(self.vdot(x=self))

Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
498
    def conjugate(self):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
499
        """ Returns the complex conjugate of the field.
Theo Steininger's avatar
Theo Steininger committed
500

501
502
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
503
504
        Field
            The complex conjugated field.
csongor's avatar
csongor committed
505
        """
Martin Reinecke's avatar
Martin Reinecke committed
506
        return Field(self._domain, self.val.conjugate())
csongor's avatar
csongor committed
507

Theo Steininger's avatar
Theo Steininger committed
508
    # ---General unary/contraction methods---
509

Theo Steininger's avatar
Theo Steininger committed
510
511
    def __pos__(self):
        return self.copy()
512

Theo Steininger's avatar
Theo Steininger committed
513
    def __neg__(self):
514
        return Field(self._domain, -self.val)
csongor's avatar
csongor committed
515

Theo Steininger's avatar
Theo Steininger committed
516
    def __abs__(self):
517
        return Field(self._domain, abs(self.val))
csongor's avatar
csongor committed
518

519
    def _contraction_helper(self, op, spaces):
Theo Steininger's avatar
Theo Steininger committed
520
        if spaces is None:
521
            return getattr(self.val, op)()
522

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

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

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

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

Theo Steininger's avatar
Theo Steininger committed
533
534
535
        # 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
536
        else:
Martin Reinecke's avatar
Martin Reinecke committed
537
            return_domain = tuple(dom
Martin Reinecke's avatar
Martin Reinecke committed
538
                                  for i, dom in enumerate(self._domain)
Theo Steininger's avatar
Theo Steininger committed
539
                                  if i not in spaces)
540

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

543
    def sum(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
544
545
546
547
548
549
550
551
552
553
554
555
556
557
        """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.
        """
558
        return self._contraction_helper('sum', spaces)
csongor's avatar
csongor committed
559

560
    def integrate(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
        """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
578
579
580
581
582
        swgt = self.scalar_weight(spaces)
        if swgt is not None:
            res = self.sum(spaces)
            res *= swgt
            return res
583
584
585
        tmp = self.weight(1, spaces=spaces)
        return tmp.sum(spaces)

586
    def prod(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
587
588
589
590
591
592
593
594
595
596
597
598
599
600
        """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.
        """
601
        return self._contraction_helper('prod', spaces)
csongor's avatar
csongor committed
602

603
604
    def all(self, spaces=None):
        return self._contraction_helper('all', spaces)
csongor's avatar
csongor committed
605

606
607
    def any(self, spaces=None):
        return self._contraction_helper('any', spaces)
csongor's avatar
csongor committed
608

609
    def min(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
610
611
612
613
614
615
616
617
618
619
620
621
622
623
        """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.
        """
624
        return self._contraction_helper('min', spaces)
csongor's avatar
csongor committed
625

626
    def max(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
627
628
629
630
631
632
633
634
635
636
637
638
639
640
        """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.
        """
641
        return self._contraction_helper('max', spaces)
csongor's avatar
csongor committed
642

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

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

693
    def std(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
        """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.
        """
711
        from .sugar import sqrt
712
713
714
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('std', spaces)
        return sqrt(self.var(spaces))
csongor's avatar
csongor committed
715

716
717
718
    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
719
        if other._domain != self._domain:
720
            raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
Martin Reinecke committed
721
        self.local_data[()] = other.local_data[()]
722

Theo Steininger's avatar
Theo Steininger committed
723
    def __repr__(self):
Philipp Arras's avatar
Philipp Arras committed
724
        return "<nifty5.Field>"
Theo Steininger's avatar
Theo Steininger committed
725
726

    def __str__(self):
Philipp Arras's avatar
Philipp Arras committed
727
        return "nifty5.Field instance\n- domain      = " + \
728
               self._domain.__str__() + \
729
               "\n- val         = " + repr(self.val)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
730

731
732
733
    def isEquivalentTo(self, other):
        """Determines (as quickly as possible) whether `self`'s content is
        identical to `other`'s content."""
734
735
736
737
738
739
740
741
        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()

742
743
744
745
746
    def isSubsetOf(self, other):
        """Identical to `Field.isEquivalentTo()`. This method is provided for
        easier interoperability with `MultiField`."""
        return self.isEquivalentTo(other)

747

748
749
750
751
752
753
754
755
756
757
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
Martin Reinecke committed
758
            global COUNTER
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
759
760
761
762
763
764
765
766
767
768
769
770
            # 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
771
772
        return func2
    setattr(Field, op, func(op))