field.py 27.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
Theo Steininger's avatar
Theo Steininger committed
13
#
Martin Reinecke's avatar
Martin Reinecke committed
14
# Copyright(C) 2013-2018 Max-Planck-Society
Theo Steininger's avatar
Theo Steininger committed
15
16
17
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
18

Martin Reinecke's avatar
Martin Reinecke committed
19
from __future__ import division
Martin Reinecke's avatar
Martin Reinecke committed
20
from builtins import range
csongor's avatar
csongor committed
21
import numpy as np
Martin Reinecke's avatar
Martin Reinecke committed
22
from . import utilities
Martin Reinecke's avatar
Martin Reinecke committed
23
from .domain_tuple import DomainTuple
Martin Reinecke's avatar
Martin Reinecke committed
24
from functools import reduce
25
from . import dobj
26

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
    @staticmethod
    def full(domain, val, dtype=None):
Martin Reinecke's avatar
Martin Reinecke committed
83
84
85
86
87
88
89
90
91
92
93
94
95
96
        """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
        """
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
        if not np.isscalar(val):
            raise TypeError("val must be a scalar")
        return Field(DomainTuple.make(domain), val, dtype)

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

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

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

    @staticmethod
    def full_like(field, val, dtype=None):
Martin Reinecke's avatar
Martin Reinecke committed
115
116
117
118
119
120
121
122
123
124
125
126
127
128
        """Creates a Field from a template, filled with a constant value.

        Parameters
        ----------
        field : Field
            the template field, from which the domain is inferred
        val : float/complex/int scalar
            fill value. Data type of the field is inferred from val.

        Returns
        -------
        Field
            the newly created field
        """
129
130
        if not isinstance(field, Field):
            raise TypeError("field must be of Field type")
Martin Reinecke's avatar
Martin Reinecke committed
131
        return Field.full(field._domain, val, dtype)
132
133
134
135
136
137
138

    @staticmethod
    def zeros_like(field, dtype=None):
        if not isinstance(field, Field):
            raise TypeError("field must be of Field type")
        if dtype is None:
            dtype = field.dtype
Martin Reinecke's avatar
Martin Reinecke committed
139
        return Field.zeros(field._domain, dtype)
140
141
142
143
144
145
146

    @staticmethod
    def ones_like(field, dtype=None):
        if not isinstance(field, Field):
            raise TypeError("field must be of Field type")
        if dtype is None:
            dtype = field.dtype
Martin Reinecke's avatar
Martin Reinecke committed
147
        return Field.ones(field._domain, dtype)
148
149
150
151
152
153
154

    @staticmethod
    def empty_like(field, dtype=None):
        if not isinstance(field, Field):
            raise TypeError("field must be of Field type")
        if dtype is None:
            dtype = field.dtype
Martin Reinecke's avatar
Martin Reinecke committed
155
        return Field.empty(field._domain, dtype)
156

157
    @staticmethod
158
    def from_global_data(domain, arr, sum_up=False):
Martin Reinecke's avatar
Martin Reinecke committed
159
160
161
162
163
164
165
166
167
        """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`.
168
169
170
171
172
        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
173
        """
174
        return Field(domain, dobj.from_global_data(arr, sum_up))
175

Martin Reinecke's avatar
Martin Reinecke committed
176
177
178
179
180
    @staticmethod
    def from_local_data(domain, arr):
        domain = DomainTuple.make(domain)
        return Field(domain, dobj.from_local_data(domain.shape, arr))

181
    def to_global_data(self):
Martin Reinecke's avatar
Martin Reinecke committed
182
183
184
185
186
187
188
189
190
191
192
193
        """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.
        """
194
195
        return dobj.to_global_data(self._val)

Martin Reinecke's avatar
Martin Reinecke committed
196
197
    @property
    def local_data(self):
Martin Reinecke's avatar
Martin Reinecke committed
198
199
200
201
202
203
204
205
206
207
        """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
208
209
        return dobj.local_data(self._val)

Martin Reinecke's avatar
Martin Reinecke committed
210
    def cast_domain(self, new_domain):
Martin Reinecke's avatar
Martin Reinecke committed
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
        """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
228
229
        return Field(new_domain, self._val)

Martin Reinecke's avatar
Martin Reinecke committed
230
    @staticmethod
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
231
    def _infer_domain(domain, val=None):
232
        if domain is None:
233
            if isinstance(val, Field):
Martin Reinecke's avatar
Martin Reinecke committed
234
                return val._domain
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
235
            if np.isscalar(val):
236
                return DomainTuple.make(())  # empty domain tuple
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
237
            raise TypeError("could not infer domain from value")
Martin Reinecke's avatar
Martin Reinecke committed
238
        return DomainTuple.make(domain)
239

Martin Reinecke's avatar
Martin Reinecke committed
240
241
    @staticmethod
    def _infer_dtype(dtype, val):
Martin Reinecke's avatar
Martin Reinecke committed
242
243
244
245
        if dtype is not None:
            return dtype
        if val is None:
            raise ValueError("could not infer dtype")
Martin Reinecke's avatar
Martin Reinecke committed
246
247
        if isinstance(val, Field):
            return val.dtype
Martin Reinecke's avatar
Martin Reinecke committed
248
        return np.result_type(val)
249

Martin Reinecke's avatar
Martin Reinecke committed
250
251
    @staticmethod
    def from_random(random_type, domain, dtype=np.float64, **kwargs):
252
253
254
255
        """ Draws a random field with the given parameters.

        Parameters
        ----------
Martin Reinecke's avatar
Martin Reinecke committed
256
257
        random_type : 'pm1', 'normal', or 'uniform'
            The random distribution to use.
Martin Reinecke's avatar
Martin Reinecke committed
258
        domain : DomainTuple
259
260
261
            The domain of the output random field
        dtype : type
            The datatype of the output random field
Theo Steininger's avatar
Theo Steininger committed
262

263
264
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
265
        Field
Martin Reinecke's avatar
Martin Reinecke committed
266
            The newly created Field.
267
        """
Martin Reinecke's avatar
Martin Reinecke committed
268
        domain = DomainTuple.make(domain)
Martin Reinecke's avatar
Martin Reinecke committed
269
270
271
        return Field(domain=domain,
                     val=dobj.from_random(random_type, dtype=dtype,
                                          shape=domain.shape, **kwargs))
272

Martin Reinecke's avatar
Martin Reinecke committed
273
    def fill(self, fill_value):
Martin Reinecke's avatar
Martin Reinecke committed
274
275
276
277
278
279
280
        """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
281
282
        self._val.fill(fill_value)

Martin Reinecke's avatar
Martin Reinecke committed
283
    def lock(self):
Martin Reinecke's avatar
Martin Reinecke committed
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
        """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
299
        dobj.lock(self._val)
300
        return self
Martin Reinecke's avatar
Martin Reinecke committed
301
302
303

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

Theo Steininger's avatar
Theo Steininger committed
307
308
    @property
    def val(self):
Martin Reinecke's avatar
Martin Reinecke committed
309
310
311
312
313
314
        """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.
315
        """
Martin Reinecke's avatar
Martin Reinecke committed
316
        return self._val
csongor's avatar
csongor committed
317

Martin Reinecke's avatar
Martin Reinecke committed
318
319
    @property
    def dtype(self):
Martin Reinecke's avatar
Martin Reinecke committed
320
        """type : the data type of the field's entries"""
Martin Reinecke's avatar
Martin Reinecke committed
321
322
        return self._val.dtype

Martin Reinecke's avatar
Martin Reinecke committed
323
324
    @property
    def domain(self):
Martin Reinecke's avatar
Martin Reinecke committed
325
        """DomainTuple : the field's domain"""
Martin Reinecke's avatar
Martin Reinecke committed
326
327
        return self._domain

328
329
    @property
    def shape(self):
Martin Reinecke's avatar
Martin Reinecke committed
330
        """tuple of int : the concatenated shapes of all sub-domains"""
Martin Reinecke's avatar
Martin Reinecke committed
331
        return self._domain.shape
csongor's avatar
csongor committed
332

333
    @property
Martin Reinecke's avatar
Martin Reinecke committed
334
    def size(self):
Martin Reinecke's avatar
Martin Reinecke committed
335
        """int : total number of pixels in the field"""
Martin Reinecke's avatar
Martin Reinecke committed
336
        return self._domain.size
csongor's avatar
csongor committed
337

Theo Steininger's avatar
Theo Steininger committed
338
339
    @property
    def real(self):
Martin Reinecke's avatar
Martin Reinecke committed
340
        """Field : The real part of the field"""
341
        if not np.issubdtype(self.dtype, np.complexfloating):
342
            return self
Martin Reinecke's avatar
Martin Reinecke committed
343
        return Field(self._domain, self.val.real)
Theo Steininger's avatar
Theo Steininger committed
344
345
346

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

Martin Reinecke's avatar
Martin Reinecke committed
352
    def copy(self):
353
        """ Returns a full copy of the Field.
Theo Steininger's avatar
Theo Steininger committed
354

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

358
359
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
360
        Field
Martin Reinecke's avatar
Martin Reinecke committed
361
            An identical, but unlocked copy of 'self'.
362
        """
Martin Reinecke's avatar
Martin Reinecke committed
363
        return Field(val=self, copy=True)
csongor's avatar
csongor committed
364

Martin Reinecke's avatar
Martin Reinecke committed
365
366
367
368
369
370
371
372
373
    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
374
            A read-only version of `self`.
Martin Reinecke's avatar
Martin Reinecke committed
375
        """
Martin Reinecke's avatar
Martin Reinecke committed
376
        return self if self.locked else Field(val=self, copy=True, locked=True)
Martin Reinecke's avatar
Martin Reinecke committed
377

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

        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
397
            spaces = range(len(self._domain))
Martin Reinecke's avatar
Martin Reinecke committed
398
        res = 1.
399
        for i in spaces:
Martin Reinecke's avatar
Martin Reinecke committed
400
            tmp = self._domain[i].scalar_dvol
401
402
403
404
405
            if tmp is None:
                return None
            res *= tmp
        return res

Martin Reinecke's avatar
Martin Reinecke committed
406
    def total_volume(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
407
408
409
410
411
412
413
414
415
416
417
418
419
        """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
420
        if np.isscalar(spaces):
Martin Reinecke's avatar
Martin Reinecke committed
421
            return self._domain[spaces].total_volume
Martin Reinecke's avatar
Martin Reinecke committed
422
423
424
425
426

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

430
    def weight(self, power=1, spaces=None, out=None):
Theo Steininger's avatar
Theo Steininger committed
431
        """ Weights the pixels of `self` with their invidual pixel-volume.
432
433
434
435

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

Martin Reinecke's avatar
Martin Reinecke committed
438
439
440
        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
441

442
443
444
445
446
        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"!

447
448
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
449
        Field
Theo Steininger's avatar
Theo Steininger committed
450
            The weighted field.
451
        """
452
453
454
455
456
        if out is None:
            out = self.copy()
        else:
            if out is not self:
                out.copy_content_from(self)
csongor's avatar
csongor committed
457

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

460
461
        fct = 1.
        for ind in spaces:
Martin Reinecke's avatar
Martin Reinecke committed
462
            wgt = self._domain[ind].dvol
463
464
465
            if np.isscalar(wgt):
                fct *= wgt
            else:
Martin Reinecke's avatar
Martin Reinecke committed
466
                new_shape = np.ones(len(self.shape), dtype=np.int)
Martin Reinecke's avatar
Martin Reinecke committed
467
468
                new_shape[self._domain.axes[ind][0]:
                          self._domain.axes[ind][-1]+1] = wgt.shape
469
                wgt = wgt.reshape(new_shape)
Martin Reinecke's avatar
Martin Reinecke committed
470
471
                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
472
                    wgt = dobj.local_data(dobj.from_global_data(wgt))
Martin Reinecke's avatar
Martin Reinecke committed
473
                out.local_data[()] *= wgt**power
474
        fct = fct**power
Martin Reinecke's avatar
Martin Reinecke committed
475
        if fct != 1.:
476
            out *= fct
477

478
        return out
csongor's avatar
csongor committed
479

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

483
484
485
        Parameters
        ----------
        x : Field
486
            x must live on the same domain as `self`.
Theo Steininger's avatar
Theo Steininger committed
487

Martin Reinecke's avatar
Martin Reinecke committed
488
        spaces : None, int or tuple of int (default: None)
489
490
            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
491

492
493
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
494
        float, complex, either scalar (for full dot products)
495
                              or Field (for partial dot products)
496
        """
497
498
499
        if not isinstance(x, Field):
            raise ValueError("The dot-partner must be an instance of " +
                             "the NIFTy field class")
Theo Steininger's avatar
Theo Steininger committed
500

Martin Reinecke's avatar
Martin Reinecke committed
501
        if x._domain != self._domain:
502
            raise ValueError("Domain mismatch")
Theo Steininger's avatar
Theo Steininger committed
503

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

        if len(spaces) == ndom:
Martin Reinecke's avatar
Martin Reinecke committed
508
            return dobj.vdot(self.val, x.val)
509
510
        # 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
511
        return (self.conjugate()*x).sum(spaces=spaces)
Theo Steininger's avatar
Theo Steininger committed
512

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

Theo Steininger's avatar
Theo Steininger committed
516
517
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
518
        float
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
519
            The L2-norm of the field values.
csongor's avatar
csongor committed
520
        """
521
        return np.sqrt(np.abs(self.vdot(x=self)))
csongor's avatar
csongor committed
522

Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
523
    def conjugate(self):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
524
        """ Returns the complex conjugate of the field.
Theo Steininger's avatar
Theo Steininger committed
525

526
527
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
528
529
        Field
            The complex conjugated field.
csongor's avatar
csongor committed
530
        """
Martin Reinecke's avatar
Martin Reinecke committed
531
        return Field(self._domain, self.val.conjugate())
csongor's avatar
csongor committed
532

Theo Steininger's avatar
Theo Steininger committed
533
    # ---General unary/contraction methods---
534

Theo Steininger's avatar
Theo Steininger committed
535
536
    def __pos__(self):
        return self.copy()
537

Theo Steininger's avatar
Theo Steininger committed
538
    def __neg__(self):
539
        return Field(self._domain, -self.val)
csongor's avatar
csongor committed
540

Theo Steininger's avatar
Theo Steininger committed
541
    def __abs__(self):
542
        return Field(self._domain, dobj.abs(self.val))
csongor's avatar
csongor committed
543

544
    def _contraction_helper(self, op, spaces):
Theo Steininger's avatar
Theo Steininger committed
545
        if spaces is None:
546
            return getattr(self.val, op)()
547

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

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

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

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

Theo Steininger's avatar
Theo Steininger committed
558
559
560
        # 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
561
        else:
Martin Reinecke's avatar
Martin Reinecke committed
562
            return_domain = tuple(dom
Martin Reinecke's avatar
Martin Reinecke committed
563
                                  for i, dom in enumerate(self._domain)
Theo Steininger's avatar
Theo Steininger committed
564
                                  if i not in spaces)
565

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

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

585
    def integrate(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
        """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
603
604
605
606
607
        swgt = self.scalar_weight(spaces)
        if swgt is not None:
            res = self.sum(spaces)
            res *= swgt
            return res
608
609
610
        tmp = self.weight(1, spaces=spaces)
        return tmp.sum(spaces)

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

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

631
632
    def any(self, spaces=None):
        return self._contraction_helper('any', spaces)
csongor's avatar
csongor committed
633

634
    def min(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
635
636
637
638
639
640
641
642
643
644
645
646
647
648
        """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.
        """
649
        return self._contraction_helper('min', spaces)
csongor's avatar
csongor committed
650

651
    def max(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
652
653
654
655
656
657
658
659
660
661
662
663
664
665
        """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.
        """
666
        return self._contraction_helper('max', spaces)
csongor's avatar
csongor committed
667

668
    def mean(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
        """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.
        """
686
687
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('mean', spaces)
Martin Reinecke's avatar
Martin Reinecke committed
688
689
690
        # MR FIXME: not very efficient
        tmp = self.weight(1)
        return tmp.sum(spaces)*(1./tmp.total_volume(spaces))
csongor's avatar
csongor committed
691

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

719
    def std(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
        """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.
        """
737
738
739
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('std', spaces)
        return sqrt(self.var(spaces))
csongor's avatar
csongor committed
740

741
742
743
    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
744
        if other._domain != self._domain:
745
            raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
Martin Reinecke committed
746
        self.local_data[()] = other.local_data[()]
747

748
    def _binary_helper(self, other, op):
csongor's avatar
csongor committed
749
        # if other is a field, make sure that the domains match
750
        if isinstance(other, Field):
Martin Reinecke's avatar
Martin Reinecke committed
751
            if other._domain != self._domain:
752
                raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
Martin Reinecke committed
753
            tval = getattr(self.val, op)(other.val)
Martin Reinecke's avatar
Martin Reinecke committed
754
            return self if tval is self.val else Field(self._domain, tval)
csongor's avatar
csongor committed
755

756
757
        if np.isscalar(other) or isinstance(other, dobj.data_object):
            tval = getattr(self.val, op)(other)
Martin Reinecke's avatar
Martin Reinecke committed
758
            return self if tval is self.val else Field(self._domain, tval)
759

Martin Reinecke's avatar
Martin Reinecke committed
760
        return NotImplemented
csongor's avatar
csongor committed
761
762

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

765
    def __radd__(self, other):
Theo Steininger's avatar
Theo Steininger committed
766
        return self._binary_helper(other, op='__radd__')
csongor's avatar
csongor committed
767
768

    def __iadd__(self, other):
769
        return self._binary_helper(other, op='__iadd__')
csongor's avatar
csongor committed
770
771

    def __sub__(self, other):
Theo Steininger's avatar
Theo Steininger committed
772
        return self._binary_helper(other, op='__sub__')
csongor's avatar
csongor committed
773
774

    def __rsub__(self, other):
Theo Steininger's avatar
Theo Steininger committed
775
        return self._binary_helper(other, op='__rsub__')
csongor's avatar
csongor committed
776
777

    def __isub__(self, other):
778
        return self._binary_helper(other, op='__isub__')
csongor's avatar
csongor committed
779
780

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

783
    def __rmul__(self, other):
Theo Steininger's avatar
Theo Steininger committed
784
        return self._binary_helper(other, op='__rmul__')
csongor's avatar
csongor committed
785
786

    def __imul__(self, other):
787
        return self._binary_helper(other, op='__imul__')
csongor's avatar
csongor committed
788
789

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

Martin Reinecke's avatar
Martin Reinecke committed
792
793
794
    def __truediv__(self, other):
        return self._binary_helper(other, op='__truediv__')

csongor's avatar
csongor committed
795
    def __rdiv__(self, other):
Theo Steininger's avatar
Theo Steininger committed
796
        return self._binary_helper(other, op='__rdiv__')
csongor's avatar
csongor committed
797

Martin Reinecke's avatar
Martin Reinecke committed
798
799
800
    def __rtruediv__(self, other):
        return self._binary_helper(other, op='__rtruediv__')

csongor's avatar
csongor committed
801
    def __idiv__(self, other):
802
        return self._binary_helper(other, op='__idiv__')
803

csongor's avatar
csongor committed
804
    def __pow__(self, other):
Theo Steininger's avatar
Theo Steininger committed
805
        return self._binary_helper(other, op='__pow__')
csongor's avatar
csongor committed
806
807

    def __rpow__(self, other):
Theo Steininger's avatar
Theo Steininger committed
808
        return self._binary_helper(other, op='__rpow__')
csongor's avatar
csongor committed
809
810

    def __ipow__(self, other):
811
        return self._binary_helper(other, op='__ipow__')
csongor's avatar
csongor committed
812

Martin Reinecke's avatar
Martin Reinecke committed
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
    def __lt__(self, other):
        return self._binary_helper(other, op='__lt__')

    def __le__(self, other):
        return self._binary_helper(other, op='__le__')

    def __ne__(self, other):
        return self._binary_helper(other, op='__ne__')

    def __eq__(self, other):
        return self._binary_helper(other, op='__eq__')

    def __ge__(self, other):
        return self._binary_helper(other, op='__ge__')

    def __gt__(self, other):
        return self._binary_helper(other, op='__gt__')

Theo Steininger's avatar
Theo Steininger committed
831
    def __repr__(self):
Martin Reinecke's avatar
Martin Reinecke committed
832
        return "<nifty4.Field>"
Theo Steininger's avatar
Theo Steininger committed
833
834

    def __str__(self):
Martin Reinecke's avatar
Martin Reinecke committed
835
        return "nifty4.Field instance\n- domain      = " + \
836
               self._domain.__str__() + \
837
               "\n- val         = " + repr(self.val)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
838
839
840
841
842
843
844
845


# Arithmetic functions working on Fields

def _math_helper(x, function, out):
    if not isinstance(x, Field):
        raise TypeError("This function only accepts Field objects.")
    if out is not None:
Martin Reinecke's avatar
Martin Reinecke committed
846
        if not isinstance(out, Field) or x._domain != out._domain:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
847
848
849
850
            raise ValueError("Bad 'out' argument")
        function(x.val, out=out.val)
        return out
    else:
Martin Reinecke's avatar
Martin Reinecke committed
851
        return Field(domain=x._domain, val=function(x.val))
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
852
853
854
855
856
857
858
859
860
861
862
863
864
865


def sqrt(x, out=None):
    return _math_helper(x, dobj.sqrt, out)


def exp(x, out=None):
    return _math_helper(x, dobj.exp, out)


def log(x, out=None):
    return _math_helper(x, dobj.log, out)


Martin Reinecke's avatar
Martin Reinecke committed
866
867
868
869
def tanh(x, out=None):
    return _math_helper(x, dobj.tanh, out)


Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
870
871
def conjugate(x, out=None):
    return _math_helper(x, dobj.conjugate, out)