field.py 26.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
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
    @staticmethod
    def full(domain, val, dtype=None):
Martin Reinecke's avatar
Martin Reinecke committed
84
85
86
87
88
89
90
91
92
93
94
95
96
97
        """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
        """
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
        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
116
117
118
119
120
121
122
123
124
125
126
127
128
129
        """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
        """
130
131
        if not isinstance(field, Field):
            raise TypeError("field must be of Field type")
Martin Reinecke's avatar
Martin Reinecke committed
132
        return Field.full(field._domain, val, dtype)
133
134
135
136
137
138
139

    @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
140
        return Field.zeros(field._domain, dtype)
141
142
143
144
145
146
147

    @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
148
        return Field.ones(field._domain, dtype)
149
150
151
152
153
154
155

    @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
156
        return Field.empty(field._domain, dtype)
157

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

479
        return out
csongor's avatar
csongor committed
480

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

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

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

493
494
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
495
        float, complex, either scalar (for full dot products)
496
                              or Field (for partial dot products)
497
        """
498
499
500
        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
501

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Theo Steininger's avatar
Theo Steininger committed
763
    def __repr__(self):
Martin Reinecke's avatar
Martin Reinecke committed
764
        return "<nifty4.Field>"
Theo Steininger's avatar
Theo Steininger committed
765
766

    def __str__(self):
Martin Reinecke's avatar
Martin Reinecke committed
767
        return "nifty4.Field instance\n- domain      = " + \
768
               self._domain.__str__() + \
769
               "\n- val         = " + repr(self.val)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
770

771
772
773
774
775
776
777
778
779
780
781
782
783
784
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):
            return self._binary_helper(other, op=op)
        return func2
    setattr(Field, op, func(op))

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
785
786
787
788

# Arithmetic functions working on Fields

def _math_helper(x, function, out):
789
    function = getattr(dobj, function)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
790
791
792
    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
793
        if not isinstance(out, Field) or x._domain != out._domain:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
794
795
796
797
            raise ValueError("Bad 'out' argument")
        function(x.val, out=out.val)
        return out
    else:
Martin Reinecke's avatar
Martin Reinecke committed
798
        return Field(domain=x._domain, val=function(x.val))
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
799

800
_current_module = sys.modules[__name__]
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
801

802
803
804
805
806
807
for f in ["sqrt", "exp", "log", "tanh", "conjugate"]:
    def func(f):
        def func2(x, out=None):
            return _math_helper(x, f, out)
        return func2
    setattr(_current_module, f, func(f))