field.py 24.2 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
238
        self._val.fill(fill_value)

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)

theos's avatar
theos 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

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

334
    def scalar_weight(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
335
336
337
338
339
340
341
342
343
344
345
346
347
348
        """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.
        """
349
        if np.isscalar(spaces):
Martin Reinecke's avatar
Martin Reinecke committed
350
            return self._domain[spaces].scalar_dvol
351
352

        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
353
            spaces = range(len(self._domain))
Martin Reinecke's avatar
Martin Reinecke committed
354
        res = 1.
355
        for i in spaces:
Martin Reinecke's avatar
Martin Reinecke committed
356
            tmp = self._domain[i].scalar_dvol
357
358
359
360
361
            if tmp is None:
                return None
            res *= tmp
        return res

Martin Reinecke's avatar
Martin Reinecke committed
362
    def total_volume(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
363
364
365
366
367
368
369
370
371
372
373
374
375
        """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
376
        if np.isscalar(spaces):
Martin Reinecke's avatar
Martin Reinecke committed
377
            return self._domain[spaces].total_volume
Martin Reinecke's avatar
Martin Reinecke committed
378
379
380
381
382

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

386
    def weight(self, power=1, spaces=None, out=None):
Theo Steininger's avatar
Theo Steininger committed
387
        """ Weights the pixels of `self` with their invidual pixel-volume.
388
389
390
391

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

Martin Reinecke's avatar
Martin Reinecke committed
394
395
396
        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
397

398
399
400
401
402
        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"!

403
404
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
405
        Field
Theo Steininger's avatar
Theo Steininger committed
406
            The weighted field.
407
        """
408
409
410
411
412
        if out is None:
            out = self.copy()
        else:
            if out is not self:
                out.copy_content_from(self)
csongor's avatar
csongor committed
413

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

416
417
        fct = 1.
        for ind in spaces:
Martin Reinecke's avatar
Martin Reinecke committed
418
            wgt = self._domain[ind].dvol
419
420
421
            if np.isscalar(wgt):
                fct *= wgt
            else:
Martin Reinecke's avatar
Martin Reinecke committed
422
                new_shape = np.ones(len(self.shape), dtype=np.int)
Martin Reinecke's avatar
Martin Reinecke committed
423
424
                new_shape[self._domain.axes[ind][0]:
                          self._domain.axes[ind][-1]+1] = wgt.shape
425
                wgt = wgt.reshape(new_shape)
Martin Reinecke's avatar
Martin Reinecke committed
426
427
                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
428
                    wgt = dobj.local_data(dobj.from_global_data(wgt))
Martin Reinecke's avatar
Martin Reinecke committed
429
                out.local_data[()] *= wgt**power
430
        fct = fct**power
Martin Reinecke's avatar
Martin Reinecke committed
431
        if fct != 1.:
432
            out *= fct
433

434
        return out
csongor's avatar
csongor committed
435

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

439
440
441
        Parameters
        ----------
        x : Field
442
            x must live on the same domain as `self`.
Theo Steininger's avatar
Theo Steininger committed
443

Martin Reinecke's avatar
Martin Reinecke committed
444
        spaces : None, int or tuple of int (default: None)
445
446
            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
447

448
449
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
450
        float, complex, either scalar (for full dot products)
451
                              or Field (for partial dot products)
452
        """
453
454
455
        if not isinstance(x, Field):
            raise ValueError("The dot-partner must be an instance of " +
                             "the NIFTy field class")
theos's avatar
theos committed
456

Martin Reinecke's avatar
Martin Reinecke committed
457
        if x._domain != self._domain:
458
            raise ValueError("Domain mismatch")
theos's avatar
theos committed
459

Martin Reinecke's avatar
Martin Reinecke committed
460
        ndom = len(self._domain)
461
462
463
        spaces = utilities.parse_spaces(spaces, ndom)

        if len(spaces) == ndom:
Martin Reinecke's avatar
Martin Reinecke committed
464
            return dobj.vdot(self.val, x.val)
465
466
        # 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
467
        return (self.conjugate()*x).sum(spaces=spaces)
theos's avatar
theos committed
468

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

Theo Steininger's avatar
Theo Steininger committed
472
473
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
474
        float
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
475
            The L2-norm of the field values.
csongor's avatar
csongor committed
476
        """
477
        return np.sqrt(abs(self.vdot(x=self)))
csongor's avatar
csongor committed
478

Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
479
    def conjugate(self):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
480
        """ Returns the complex conjugate of the field.
Theo Steininger's avatar
Theo Steininger committed
481

482
483
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
484
485
        Field
            The complex conjugated field.
csongor's avatar
csongor committed
486
        """
Martin Reinecke's avatar
Martin Reinecke committed
487
        return Field(self._domain, self.val.conjugate())
csongor's avatar
csongor committed
488

theos's avatar
theos committed
489
    # ---General unary/contraction methods---
490

theos's avatar
theos committed
491
492
    def __pos__(self):
        return self.copy()
493

theos's avatar
theos committed
494
    def __neg__(self):
495
        return Field(self._domain, -self.val)
csongor's avatar
csongor committed
496

theos's avatar
theos committed
497
    def __abs__(self):
498
        return Field(self._domain, abs(self.val))
csongor's avatar
csongor committed
499

500
    def _contraction_helper(self, op, spaces):
theos's avatar
theos committed
501
        if spaces is None:
502
            return getattr(self.val, op)()
503

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

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

Martin Reinecke's avatar
Martin Reinecke committed
508
        if len(axes_list) > 0:
theos's avatar
theos committed
509
            axes_list = reduce(lambda x, y: x+y, axes_list)
csongor's avatar
csongor committed
510

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

theos's avatar
theos committed
514
515
516
        # 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
517
        else:
Martin Reinecke's avatar
Martin Reinecke committed
518
            return_domain = tuple(dom
Martin Reinecke's avatar
Martin Reinecke committed
519
                                  for i, dom in enumerate(self._domain)
theos's avatar
theos committed
520
                                  if i not in spaces)
521

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

524
    def sum(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
525
526
527
528
529
530
531
532
533
534
535
536
537
538
        """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.
        """
539
        return self._contraction_helper('sum', spaces)
csongor's avatar
csongor committed
540

541
    def integrate(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
        """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
559
560
561
562
563
        swgt = self.scalar_weight(spaces)
        if swgt is not None:
            res = self.sum(spaces)
            res *= swgt
            return res
564
565
566
        tmp = self.weight(1, spaces=spaces)
        return tmp.sum(spaces)

567
    def prod(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
568
569
570
571
572
573
574
575
576
577
578
579
580
581
        """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.
        """
582
        return self._contraction_helper('prod', spaces)
csongor's avatar
csongor committed
583

584
585
    def all(self, spaces=None):
        return self._contraction_helper('all', spaces)
csongor's avatar
csongor committed
586

587
588
    def any(self, spaces=None):
        return self._contraction_helper('any', spaces)
csongor's avatar
csongor committed
589

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

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

624
    def mean(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
        """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.
        """
642
643
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('mean', spaces)
Martin Reinecke's avatar
Martin Reinecke committed
644
645
646
        # MR FIXME: not very efficient
        tmp = self.weight(1)
        return tmp.sum(spaces)*(1./tmp.total_volume(spaces))
csongor's avatar
csongor committed
647

648
    def var(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
649
650
651
652
653
654
655
656
657
658
659
660
661
662
        """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.
        """
663
664
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('var', spaces)
Martin Reinecke's avatar
Martin Reinecke committed
665
666
667
668
669
670
671
672
673
        # MR FIXME: not very efficient or accurate
        m1 = self.mean(spaces)
        if np.issubdtype(self.dtype, np.complexfloating):
            sq = abs(self)**2
            m1 = abs(m1)**2
        else:
            sq = self**2
            m1 **= 2
        return sq.mean(spaces) - m1
csongor's avatar
csongor committed
674

675
    def std(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
        """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.
        """
693
694
695
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('std', spaces)
        return sqrt(self.var(spaces))
csongor's avatar
csongor committed
696

697
698
699
    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
700
        if other._domain != self._domain:
701
            raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
Martin Reinecke committed
702
        self.local_data[()] = other.local_data[()]
703

theos's avatar
theos committed
704
    def __repr__(self):
Martin Reinecke's avatar
Martin Reinecke committed
705
        return "<nifty4.Field>"
theos's avatar
theos committed
706
707

    def __str__(self):
Martin Reinecke's avatar
Martin Reinecke committed
708
        return "nifty4.Field instance\n- domain      = " + \
709
               self._domain.__str__() + \
710
               "\n- val         = " + repr(self.val)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
711

712
713
714
715
716
717
718
719
720
721
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
722
723
724
725
726
727
728
729
730
731
732
733
            # 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
734
735
        return func2
    setattr(Field, op, func(op))