field.py 26.9 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
Martin Reinecke's avatar
Martin Reinecke committed
158
159
160
161
162
163
164
165
166
167
168
169
170
171
    def from_global_data(domain, arr):
        """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`.
            If MPI is active, the contents of `arr` must be the same on all
            MPI tasks.
        """
        return Field(domain, dobj.from_global_data(arr))
172
173

    def to_global_data(self):
Martin Reinecke's avatar
Martin Reinecke committed
174
175
176
177
178
179
180
181
182
183
184
185
        """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.
        """
186
187
        return dobj.to_global_data(self._val)

Martin Reinecke's avatar
Martin Reinecke committed
188
189
    @property
    def local_data(self):
Martin Reinecke's avatar
Martin Reinecke committed
190
191
192
193
194
195
196
197
198
199
        """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
200
201
        return dobj.local_data(self._val)

Martin Reinecke's avatar
Martin Reinecke committed
202
    def cast_domain(self, new_domain):
Martin Reinecke's avatar
Martin Reinecke committed
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
        """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
220
221
        return Field(new_domain, self._val)

Martin Reinecke's avatar
Martin Reinecke committed
222
    @staticmethod
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
223
    def _infer_domain(domain, val=None):
224
        if domain is None:
225
            if isinstance(val, Field):
Martin Reinecke's avatar
Martin Reinecke committed
226
                return val._domain
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
227
            if np.isscalar(val):
228
                return DomainTuple.make(())  # empty domain tuple
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
229
            raise TypeError("could not infer domain from value")
Martin Reinecke's avatar
Martin Reinecke committed
230
        return DomainTuple.make(domain)
231

Martin Reinecke's avatar
Martin Reinecke committed
232
233
    @staticmethod
    def _infer_dtype(dtype, val):
Martin Reinecke's avatar
Martin Reinecke committed
234
235
236
237
        if dtype is not None:
            return dtype
        if val is None:
            raise ValueError("could not infer dtype")
Martin Reinecke's avatar
Martin Reinecke committed
238
239
        if isinstance(val, Field):
            return val.dtype
Martin Reinecke's avatar
Martin Reinecke committed
240
        return np.result_type(val)
241

Martin Reinecke's avatar
Martin Reinecke committed
242
243
    @staticmethod
    def from_random(random_type, domain, dtype=np.float64, **kwargs):
244
245
246
247
        """ Draws a random field with the given parameters.

        Parameters
        ----------
Martin Reinecke's avatar
Martin Reinecke committed
248
249
        random_type : 'pm1', 'normal', or 'uniform'
            The random distribution to use.
Martin Reinecke's avatar
Martin Reinecke committed
250
        domain : DomainTuple
251
252
253
            The domain of the output random field
        dtype : type
            The datatype of the output random field
Theo Steininger's avatar
Theo Steininger committed
254

255
256
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
257
        Field
Martin Reinecke's avatar
Martin Reinecke committed
258
            The newly created Field.
259
        """
Martin Reinecke's avatar
Martin Reinecke committed
260
        domain = DomainTuple.make(domain)
Martin Reinecke's avatar
Martin Reinecke committed
261
262
263
        return Field(domain=domain,
                     val=dobj.from_random(random_type, dtype=dtype,
                                          shape=domain.shape, **kwargs))
264

Martin Reinecke's avatar
Martin Reinecke committed
265
    def fill(self, fill_value):
Martin Reinecke's avatar
Martin Reinecke committed
266
267
268
269
270
271
272
        """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
273
274
        self._val.fill(fill_value)

Martin Reinecke's avatar
Martin Reinecke committed
275
    def lock(self):
Martin Reinecke's avatar
Martin Reinecke committed
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
        """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
291
        dobj.lock(self._val)
292
        return self
Martin Reinecke's avatar
Martin Reinecke committed
293
294
295

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

Theo Steininger's avatar
Theo Steininger committed
299
300
    @property
    def val(self):
Martin Reinecke's avatar
Martin Reinecke committed
301
302
303
304
305
306
        """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.
307
        """
Martin Reinecke's avatar
Martin Reinecke committed
308
        return self._val
csongor's avatar
csongor committed
309

Martin Reinecke's avatar
Martin Reinecke committed
310
311
    @property
    def dtype(self):
Martin Reinecke's avatar
Martin Reinecke committed
312
        """type : the data type of the field's entries"""
Martin Reinecke's avatar
Martin Reinecke committed
313
314
        return self._val.dtype

Martin Reinecke's avatar
Martin Reinecke committed
315
316
    @property
    def domain(self):
Martin Reinecke's avatar
Martin Reinecke committed
317
        """DomainTuple : the field's domain"""
Martin Reinecke's avatar
Martin Reinecke committed
318
319
        return self._domain

320
321
    @property
    def shape(self):
Martin Reinecke's avatar
Martin Reinecke committed
322
        """tuple of int : the concatenated shapes of all sub-domains"""
Martin Reinecke's avatar
Martin Reinecke committed
323
        return self._domain.shape
csongor's avatar
csongor committed
324

325
    @property
Martin Reinecke's avatar
Martin Reinecke committed
326
    def size(self):
Martin Reinecke's avatar
Martin Reinecke committed
327
        """int : total number of pixels in the field"""
Martin Reinecke's avatar
Martin Reinecke committed
328
        return self._domain.size
csongor's avatar
csongor committed
329

Theo Steininger's avatar
Theo Steininger committed
330
331
    @property
    def real(self):
Martin Reinecke's avatar
Martin Reinecke committed
332
        """Field : The real part of the field"""
333
        if not np.issubdtype(self.dtype, np.complexfloating):
334
            return self
Martin Reinecke's avatar
Martin Reinecke committed
335
        return Field(self._domain, self.val.real)
Theo Steininger's avatar
Theo Steininger committed
336
337
338

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

Martin Reinecke's avatar
Martin Reinecke committed
344
    def copy(self):
345
        """ Returns a full copy of the Field.
Theo Steininger's avatar
Theo Steininger committed
346

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

350
351
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
352
        Field
Martin Reinecke's avatar
Martin Reinecke committed
353
            An identical, but unlocked copy of 'self'.
354
        """
Martin Reinecke's avatar
Martin Reinecke committed
355
        return Field(val=self, copy=True)
csongor's avatar
csongor committed
356

Martin Reinecke's avatar
Martin Reinecke committed
357
358
359
360
361
362
363
364
365
    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
366
            A read-only version of `self`.
Martin Reinecke's avatar
Martin Reinecke committed
367
        """
Martin Reinecke's avatar
Martin Reinecke committed
368
        return self if self.locked else Field(val=self, copy=True, locked=True)
Martin Reinecke's avatar
Martin Reinecke committed
369

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

        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
389
            spaces = range(len(self._domain))
Martin Reinecke's avatar
Martin Reinecke committed
390
        res = 1.
391
        for i in spaces:
Martin Reinecke's avatar
Martin Reinecke committed
392
            tmp = self._domain[i].scalar_dvol
393
394
395
396
397
            if tmp is None:
                return None
            res *= tmp
        return res

Martin Reinecke's avatar
Martin Reinecke committed
398
    def total_volume(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
399
400
401
402
403
404
405
406
407
408
409
410
411
        """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
412
        if np.isscalar(spaces):
Martin Reinecke's avatar
Martin Reinecke committed
413
            return self._domain[spaces].total_volume
Martin Reinecke's avatar
Martin Reinecke committed
414
415
416
417
418

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

422
    def weight(self, power=1, spaces=None, out=None):
Theo Steininger's avatar
Theo Steininger committed
423
        """ Weights the pixels of `self` with their invidual pixel-volume.
424
425
426
427

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

Martin Reinecke's avatar
Martin Reinecke committed
430
431
432
        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
433

434
435
436
437
438
        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"!

439
440
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
441
        Field
Theo Steininger's avatar
Theo Steininger committed
442
            The weighted field.
443
        """
444
445
446
447
448
        if out is None:
            out = self.copy()
        else:
            if out is not self:
                out.copy_content_from(self)
csongor's avatar
csongor committed
449

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

452
453
        fct = 1.
        for ind in spaces:
Martin Reinecke's avatar
Martin Reinecke committed
454
            wgt = self._domain[ind].dvol
455
456
457
            if np.isscalar(wgt):
                fct *= wgt
            else:
Martin Reinecke's avatar
Martin Reinecke committed
458
                new_shape = np.ones(len(self.shape), dtype=np.int)
Martin Reinecke's avatar
Martin Reinecke committed
459
460
                new_shape[self._domain.axes[ind][0]:
                          self._domain.axes[ind][-1]+1] = wgt.shape
461
                wgt = wgt.reshape(new_shape)
Martin Reinecke's avatar
Martin Reinecke committed
462
463
                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
464
                    wgt = dobj.local_data(dobj.from_global_data(wgt))
Martin Reinecke's avatar
Martin Reinecke committed
465
                out.local_data[()] *= wgt**power
466
        fct = fct**power
Martin Reinecke's avatar
Martin Reinecke committed
467
        if fct != 1.:
468
            out *= fct
469

470
        return out
csongor's avatar
csongor committed
471

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

475
476
477
        Parameters
        ----------
        x : Field
478
            x must live on the same domain as `self`.
Theo Steininger's avatar
Theo Steininger committed
479

Martin Reinecke's avatar
Martin Reinecke committed
480
        spaces : None, int or tuple of int (default: None)
481
482
            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
483

484
485
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
486
        float, complex, either scalar (for full dot products)
487
                              or Field (for partial dot products)
488
        """
489
490
491
        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
492

Martin Reinecke's avatar
Martin Reinecke committed
493
        if x._domain != self._domain:
494
            raise ValueError("Domain mismatch")
Theo Steininger's avatar
Theo Steininger committed
495

Martin Reinecke's avatar
Martin Reinecke committed
496
        ndom = len(self._domain)
497
498
499
        spaces = utilities.parse_spaces(spaces, ndom)

        if len(spaces) == ndom:
Martin Reinecke's avatar
Martin Reinecke committed
500
            return dobj.vdot(self.val, x.val)
501
502
        # 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
503
        return (self.conjugate()*x).sum(spaces=spaces)
Theo Steininger's avatar
Theo Steininger committed
504

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

Theo Steininger's avatar
Theo Steininger committed
508
509
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
510
        float
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
511
            The L2-norm of the field values.
csongor's avatar
csongor committed
512
        """
513
        return np.sqrt(np.abs(self.vdot(x=self)))
csongor's avatar
csongor committed
514

Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
515
    def conjugate(self):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
516
        """ Returns the complex conjugate of the field.
Theo Steininger's avatar
Theo Steininger committed
517

518
519
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
520
521
        Field
            The complex conjugated field.
csongor's avatar
csongor committed
522
        """
Martin Reinecke's avatar
Martin Reinecke committed
523
        return Field(self._domain, self.val.conjugate())
csongor's avatar
csongor committed
524

Theo Steininger's avatar
Theo Steininger committed
525
    # ---General unary/contraction methods---
526

Theo Steininger's avatar
Theo Steininger committed
527
528
    def __pos__(self):
        return self.copy()
529

Theo Steininger's avatar
Theo Steininger committed
530
    def __neg__(self):
531
        return Field(self._domain, -self.val)
csongor's avatar
csongor committed
532

Theo Steininger's avatar
Theo Steininger committed
533
    def __abs__(self):
534
        return Field(self._domain, dobj.abs(self.val))
csongor's avatar
csongor committed
535

536
    def _contraction_helper(self, op, spaces):
Theo Steininger's avatar
Theo Steininger committed
537
        if spaces is None:
538
            return getattr(self.val, op)()
539

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

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

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

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

Theo Steininger's avatar
Theo Steininger committed
550
551
552
        # 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
553
        else:
Martin Reinecke's avatar
Martin Reinecke committed
554
            return_domain = tuple(dom
Martin Reinecke's avatar
Martin Reinecke committed
555
                                  for i, dom in enumerate(self._domain)
Theo Steininger's avatar
Theo Steininger committed
556
                                  if i not in spaces)
557

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

560
    def sum(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
561
562
563
564
565
566
567
568
569
570
571
572
573
574
        """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.
        """
575
        return self._contraction_helper('sum', spaces)
csongor's avatar
csongor committed
576

577
    def integrate(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
        """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
595
596
597
598
599
        swgt = self.scalar_weight(spaces)
        if swgt is not None:
            res = self.sum(spaces)
            res *= swgt
            return res
600
601
602
        tmp = self.weight(1, spaces=spaces)
        return tmp.sum(spaces)

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

620
621
    def all(self, spaces=None):
        return self._contraction_helper('all', spaces)
csongor's avatar
csongor committed
622

623
624
    def any(self, spaces=None):
        return self._contraction_helper('any', spaces)
csongor's avatar
csongor committed
625

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

643
    def max(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
644
645
646
647
648
649
650
651
652
653
654
655
656
657
        """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.
        """
658
        return self._contraction_helper('max', spaces)
csongor's avatar
csongor committed
659

660
    def mean(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
        """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.
        """
678
679
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('mean', spaces)
Martin Reinecke's avatar
Martin Reinecke committed
680
681
682
        # MR FIXME: not very efficient
        tmp = self.weight(1)
        return tmp.sum(spaces)*(1./tmp.total_volume(spaces))
csongor's avatar
csongor committed
683

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

711
    def std(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
        """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.
        """
729
730
731
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('std', spaces)
        return sqrt(self.var(spaces))
csongor's avatar
csongor committed
732

733
734
735
    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
736
        if other._domain != self._domain:
737
            raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
Martin Reinecke committed
738
        self.local_data[()] = other.local_data[()]
739

740
    def _binary_helper(self, other, op):
csongor's avatar
csongor committed
741
        # if other is a field, make sure that the domains match
742
        if isinstance(other, Field):
Martin Reinecke's avatar
Martin Reinecke committed
743
            if other._domain != self._domain:
744
                raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
Martin Reinecke committed
745
            tval = getattr(self.val, op)(other.val)
Martin Reinecke's avatar
Martin Reinecke committed
746
            return self if tval is self.val else Field(self._domain, tval)
csongor's avatar
csongor committed
747

748
749
        if np.isscalar(other) or isinstance(other, dobj.data_object):
            tval = getattr(self.val, op)(other)
Martin Reinecke's avatar
Martin Reinecke committed
750
            return self if tval is self.val else Field(self._domain, tval)
751

Martin Reinecke's avatar
Martin Reinecke committed
752
        return NotImplemented
csongor's avatar
csongor committed
753
754

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

757
    def __radd__(self, other):
Theo Steininger's avatar
Theo Steininger committed
758
        return self._binary_helper(other, op='__radd__')
csongor's avatar
csongor committed
759
760

    def __iadd__(self, other):
761
        return self._binary_helper(other, op='__iadd__')
csongor's avatar
csongor committed
762
763

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

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

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

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

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

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

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

Martin Reinecke's avatar
Martin Reinecke committed
784
785
786
    def __truediv__(self, other):
        return self._binary_helper(other, op='__truediv__')

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

Martin Reinecke's avatar
Martin Reinecke committed
790
791
792
    def __rtruediv__(self, other):
        return self._binary_helper(other, op='__rtruediv__')

csongor's avatar
csongor committed
793
    def __idiv__(self, other):
794
        return self._binary_helper(other, op='__idiv__')
795

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

    def __rpow__(self, other):
Theo Steininger's avatar
Theo Steininger committed
800
        return self._binary_helper(other, op='__rpow__')
csongor's avatar
csongor committed
801
802

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

Theo Steininger's avatar
Theo Steininger committed
805
    def __repr__(self):
Martin Reinecke's avatar
Martin Reinecke committed
806
        return "<nifty4.Field>"
Theo Steininger's avatar
Theo Steininger committed
807
808

    def __str__(self):
Martin Reinecke's avatar
Martin Reinecke committed
809
        return "nifty4.Field instance\n- domain      = " + \
810
               self._domain.__str__() + \
811
               "\n- val         = " + repr(self.val)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
812
813
814
815
816
817
818
819


# 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
820
        if not isinstance(out, Field) or x._domain != out._domain:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
821
822
823
824
            raise ValueError("Bad 'out' argument")
        function(x.val, out=out.val)
        return out
    else:
Martin Reinecke's avatar
Martin Reinecke committed
825
        return Field(domain=x._domain, val=function(x.val))
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
826
827
828
829
830
831
832
833
834
835
836
837
838
839


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
840
841
842
843
def tanh(x, out=None):
    return _math_helper(x, dobj.tanh, out)


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