field.py 18.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

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.
Theo Steininger's avatar
Theo Steininger committed
47

Martin Reinecke's avatar
Martin Reinecke committed
48
    copy : bool
49
    """
50

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

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

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

78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
    @staticmethod
    def full(domain, val, dtype=None):
        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):
        if not isinstance(field, Field):
            raise TypeError("field must be of Field type")
Martin Reinecke's avatar
Martin Reinecke committed
100
        return Field.full(field._domain, val, dtype)
101
102
103
104
105
106
107

    @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
108
        return Field.zeros(field._domain, dtype)
109
110
111
112
113
114
115

    @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
116
        return Field.ones(field._domain, dtype)
117
118
119
120
121
122
123

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

126
127
128
129
130
131
132
    @staticmethod
    def from_global_data(domain, dobject):
        return Field(domain, dobj.from_global_data(dobject))

    def to_global_data(self):
        return dobj.to_global_data(self._val)

Martin Reinecke's avatar
Martin Reinecke committed
133
134
135
136
    @property
    def local_data(self):
        return dobj.local_data(self._val)

Martin Reinecke's avatar
Martin Reinecke committed
137
138
139
    def cast_domain(self, new_domain):
        return Field(new_domain, self._val)

Martin Reinecke's avatar
Martin Reinecke committed
140
    @staticmethod
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
141
    def _infer_domain(domain, val=None):
142
        if domain is None:
143
            if isinstance(val, Field):
Martin Reinecke's avatar
Martin Reinecke committed
144
                return val._domain
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
145
            if np.isscalar(val):
146
                return DomainTuple.make(())  # empty domain tuple
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
147
            raise TypeError("could not infer domain from value")
Martin Reinecke's avatar
Martin Reinecke committed
148
        return DomainTuple.make(domain)
149

Martin Reinecke's avatar
Martin Reinecke committed
150
151
    @staticmethod
    def _infer_dtype(dtype, val):
Martin Reinecke's avatar
Martin Reinecke committed
152
153
154
155
        if dtype is not None:
            return dtype
        if val is None:
            raise ValueError("could not infer dtype")
Martin Reinecke's avatar
Martin Reinecke committed
156
157
        if isinstance(val, Field):
            return val.dtype
Martin Reinecke's avatar
Martin Reinecke committed
158
        return np.result_type(val)
159

Martin Reinecke's avatar
Martin Reinecke committed
160
161
    @staticmethod
    def from_random(random_type, domain, dtype=np.float64, **kwargs):
162
163
164
165
        """ Draws a random field with the given parameters.

        Parameters
        ----------
Martin Reinecke's avatar
Martin Reinecke committed
166
        random_type : str
167
168
            'pm1', 'normal', 'uniform' are the supported arguments for this
            method.
Theo Steininger's avatar
Theo Steininger committed
169

Martin Reinecke's avatar
Martin Reinecke committed
170
        domain : DomainTuple
171
            The domain of the output random field
Theo Steininger's avatar
Theo Steininger committed
172

173
174
        dtype : type
            The datatype of the output random field
Theo Steininger's avatar
Theo Steininger committed
175

176
177
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
178
        Field
179
180
            The output object.
        """
Martin Reinecke's avatar
Martin Reinecke committed
181
        domain = DomainTuple.make(domain)
Martin Reinecke's avatar
Martin Reinecke committed
182
183
184
        return Field(domain=domain,
                     val=dobj.from_random(random_type, dtype=dtype,
                                          shape=domain.shape, **kwargs))
185

Martin Reinecke's avatar
Martin Reinecke committed
186
187
188
    def fill(self, fill_value):
        self._val.fill(fill_value)

Martin Reinecke's avatar
Martin Reinecke committed
189
190
    def lock(self):
        dobj.lock(self._val)
191
        return self
Martin Reinecke's avatar
Martin Reinecke committed
192
193
194
195
196

    @property
    def locked(self):
        return dobj.locked(self._val)

Theo Steininger's avatar
Theo Steininger committed
197
198
    @property
    def val(self):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
199
        """ Returns the data object associated with this Field.
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
200
        No copy is made.
201
        """
Martin Reinecke's avatar
Martin Reinecke committed
202
        return self._val
csongor's avatar
csongor committed
203

Martin Reinecke's avatar
Martin Reinecke committed
204
205
206
207
    @property
    def dtype(self):
        return self._val.dtype

Martin Reinecke's avatar
Martin Reinecke committed
208
209
210
211
    @property
    def domain(self):
        return self._domain

212
213
    @property
    def shape(self):
Theo Steininger's avatar
Theo Steininger committed
214
        """ Returns the total shape of the Field's data array.
Theo Steininger's avatar
Theo Steininger committed
215

216
217
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
218
219
        tuple of int
            the dimensions of the spaces in domain.
Martin Reinecke's avatar
Martin Reinecke committed
220
        """
Martin Reinecke's avatar
Martin Reinecke committed
221
        return self._domain.shape
csongor's avatar
csongor committed
222

223
    @property
Martin Reinecke's avatar
Martin Reinecke committed
224
    def size(self):
Theo Steininger's avatar
Theo Steininger committed
225
        """ Returns the total number of pixel-dimensions the field has.
Theo Steininger's avatar
Theo Steininger committed
226

Theo Steininger's avatar
Theo Steininger committed
227
        Effectively, all values from shape are multiplied.
Theo Steininger's avatar
Theo Steininger committed
228

229
230
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
231
232
        int
            the dimension of the Field.
233
        """
Martin Reinecke's avatar
Martin Reinecke committed
234
        return self._domain.size
csongor's avatar
csongor committed
235

Theo Steininger's avatar
Theo Steininger committed
236
237
    @property
    def real(self):
Martin Reinecke's avatar
Martin Reinecke committed
238
        """ The real part of the field (data is not copied)."""
239
        if not np.issubdtype(self.dtype, np.complexfloating):
240
            return self
Martin Reinecke's avatar
Martin Reinecke committed
241
        return Field(self._domain, self.val.real)
Theo Steininger's avatar
Theo Steininger committed
242
243
244

    @property
    def imag(self):
Martin Reinecke's avatar
Martin Reinecke committed
245
        """ The imaginary part of the field (data is not copied)."""
246
247
        if not np.issubdtype(self.dtype, np.complexfloating):
            raise ValueError(".imag called on a non-complex Field")
Martin Reinecke's avatar
Martin Reinecke committed
248
        return Field(self._domain, self.val.imag)
Theo Steininger's avatar
Theo Steininger committed
249

Martin Reinecke's avatar
Martin Reinecke committed
250
    def copy(self):
251
        """ Returns a full copy of the Field.
Theo Steininger's avatar
Theo Steininger committed
252

Martin Reinecke's avatar
Martin Reinecke committed
253
        The returned object will be an identical copy of the original 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
            An identical copy of 'self'.
259
        """
Martin Reinecke's avatar
Martin Reinecke committed
260
        return Field(val=self, copy=True)
csongor's avatar
csongor committed
261

Martin Reinecke's avatar
Martin Reinecke committed
262
263
264
265
266
267
268
269
270
271
272
273
274
    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
            A read-only version of 'self'.
        """
        if self.locked:
            return self
Martin Reinecke's avatar
Martin Reinecke committed
275
        return Field(val=self, copy=True, locked=True)
Martin Reinecke's avatar
Martin Reinecke committed
276

277
278
    def scalar_weight(self, spaces=None):
        if np.isscalar(spaces):
Martin Reinecke's avatar
Martin Reinecke committed
279
            return self._domain[spaces].scalar_dvol
280
281

        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
282
            spaces = range(len(self._domain))
Martin Reinecke's avatar
Martin Reinecke committed
283
        res = 1.
284
        for i in spaces:
Martin Reinecke's avatar
Martin Reinecke committed
285
            tmp = self._domain[i].scalar_dvol
286
287
288
289
290
            if tmp is None:
                return None
            res *= tmp
        return res

Martin Reinecke's avatar
Martin Reinecke committed
291
292
    def total_volume(self, spaces=None):
        if np.isscalar(spaces):
Martin Reinecke's avatar
Martin Reinecke committed
293
            return self._domain[spaces].total_volume
Martin Reinecke's avatar
Martin Reinecke committed
294
295
296
297
298

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

302
    def weight(self, power=1, spaces=None, out=None):
Theo Steininger's avatar
Theo Steininger committed
303
        """ Weights the pixels of `self` with their invidual pixel-volume.
304
305
306
307

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

Martin Reinecke's avatar
Martin Reinecke committed
310
        spaces : int or tuple of int
Theo Steininger's avatar
Theo Steininger committed
311
            Determines on which subspace the operation takes place.
Theo Steininger's avatar
Theo Steininger committed
312

313
314
315
316
317
        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"!

318
319
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
320
        Field
Theo Steininger's avatar
Theo Steininger committed
321
            The weighted field.
322
        """
323
324
325
326
327
        if out is None:
            out = self.copy()
        else:
            if out is not self:
                out.copy_content_from(self)
csongor's avatar
csongor committed
328

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

331
332
        fct = 1.
        for ind in spaces:
Martin Reinecke's avatar
Martin Reinecke committed
333
            wgt = self._domain[ind].dvol
334
335
336
            if np.isscalar(wgt):
                fct *= wgt
            else:
Martin Reinecke's avatar
Martin Reinecke committed
337
                new_shape = np.ones(len(self.shape), dtype=np.int)
Martin Reinecke's avatar
Martin Reinecke committed
338
339
                new_shape[self._domain.axes[ind][0]:
                          self._domain.axes[ind][-1]+1] = wgt.shape
340
                wgt = wgt.reshape(new_shape)
Martin Reinecke's avatar
Martin Reinecke committed
341
342
                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
343
                    wgt = dobj.local_data(dobj.from_global_data(wgt))
Martin Reinecke's avatar
Martin Reinecke committed
344
                out.local_data[()] *= wgt**power
345
        fct = fct**power
Martin Reinecke's avatar
Martin Reinecke committed
346
        if fct != 1.:
347
            out *= fct
348

349
        return out
csongor's avatar
csongor committed
350

Martin Reinecke's avatar
Martin Reinecke committed
351
    def vdot(self, x=None, spaces=None):
Theo Steininger's avatar
Theo Steininger committed
352
        """ Computes the volume-factor-aware dot product of 'self' with x.
Theo Steininger's avatar
Theo Steininger committed
353

354
355
356
        Parameters
        ----------
        x : Field
357
            x must live on the same domain as `self`.
Theo Steininger's avatar
Theo Steininger committed
358

Martin Reinecke's avatar
Martin Reinecke committed
359
        spaces : None, int or tuple of int (default: None)
360
361
            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
362

363
364
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
365
        float, complex, either scalar (for full dot products)
366
                              or Field (for partial dot products)
367
        """
368
369
370
        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
371

Martin Reinecke's avatar
Martin Reinecke committed
372
        if x._domain != self._domain:
373
            raise ValueError("Domain mismatch")
Theo Steininger's avatar
Theo Steininger committed
374

Martin Reinecke's avatar
Martin Reinecke committed
375
        ndom = len(self._domain)
376
377
378
        spaces = utilities.parse_spaces(spaces, ndom)

        if len(spaces) == ndom:
Martin Reinecke's avatar
Martin Reinecke committed
379
            return dobj.vdot(self.val, x.val)
380
381
        # 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
382
        return (self.conjugate()*x).sum(spaces=spaces)
Theo Steininger's avatar
Theo Steininger committed
383

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

Theo Steininger's avatar
Theo Steininger committed
387
388
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
389
        float
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
390
            The L2-norm of the field values.
csongor's avatar
csongor committed
391
        """
392
        return np.sqrt(np.abs(self.vdot(x=self)))
csongor's avatar
csongor committed
393

Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
394
    def conjugate(self):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
395
        """ Returns the complex conjugate of the field.
Theo Steininger's avatar
Theo Steininger committed
396

397
398
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
399
400
        Field
            The complex conjugated field.
csongor's avatar
csongor committed
401
        """
Martin Reinecke's avatar
Martin Reinecke committed
402
        return Field(self._domain, self.val.conjugate())
csongor's avatar
csongor committed
403

Theo Steininger's avatar
Theo Steininger committed
404
    # ---General unary/contraction methods---
405

Theo Steininger's avatar
Theo Steininger committed
406
407
    def __pos__(self):
        return self.copy()
408

Theo Steininger's avatar
Theo Steininger committed
409
    def __neg__(self):
410
        return Field(self._domain, -self.val)
csongor's avatar
csongor committed
411

Theo Steininger's avatar
Theo Steininger committed
412
    def __abs__(self):
413
        return Field(self._domain, dobj.abs(self.val))
csongor's avatar
csongor committed
414

415
    def _contraction_helper(self, op, spaces):
Theo Steininger's avatar
Theo Steininger committed
416
        if spaces is None:
417
            return getattr(self.val, op)()
418

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

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

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

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

Theo Steininger's avatar
Theo Steininger committed
429
430
431
        # 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
432
        else:
Martin Reinecke's avatar
Martin Reinecke committed
433
            return_domain = tuple(dom
Martin Reinecke's avatar
Martin Reinecke committed
434
                                  for i, dom in enumerate(self._domain)
Theo Steininger's avatar
Theo Steininger committed
435
                                  if i not in spaces)
436

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

439
440
    def sum(self, spaces=None):
        return self._contraction_helper('sum', spaces)
csongor's avatar
csongor committed
441

442
    def integrate(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
443
444
445
446
447
        swgt = self.scalar_weight(spaces)
        if swgt is not None:
            res = self.sum(spaces)
            res *= swgt
            return res
448
449
450
        tmp = self.weight(1, spaces=spaces)
        return tmp.sum(spaces)

451
452
    def prod(self, spaces=None):
        return self._contraction_helper('prod', spaces)
csongor's avatar
csongor committed
453

454
455
    def all(self, spaces=None):
        return self._contraction_helper('all', spaces)
csongor's avatar
csongor committed
456

457
458
    def any(self, spaces=None):
        return self._contraction_helper('any', spaces)
csongor's avatar
csongor committed
459

460
461
    def min(self, spaces=None):
        return self._contraction_helper('min', spaces)
csongor's avatar
csongor committed
462

463
464
    def max(self, spaces=None):
        return self._contraction_helper('max', spaces)
csongor's avatar
csongor committed
465

466
    def mean(self, spaces=None):
467
468
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('mean', spaces)
Martin Reinecke's avatar
Martin Reinecke committed
469
470
471
        # MR FIXME: not very efficient
        tmp = self.weight(1)
        return tmp.sum(spaces)*(1./tmp.total_volume(spaces))
csongor's avatar
csongor committed
472

473
    def var(self, spaces=None):
474
475
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('var', spaces)
Martin Reinecke's avatar
Martin Reinecke committed
476
477
478
479
480
481
482
483
484
        # 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
485

486
    def std(self, spaces=None):
487
488
489
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('std', spaces)
        return sqrt(self.var(spaces))
csongor's avatar
csongor committed
490

491
492
493
    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
494
        if other._domain != self._domain:
495
            raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
Martin Reinecke committed
496
        self.local_data[()] = other.local_data[()]
497

498
    def _binary_helper(self, other, op):
csongor's avatar
csongor committed
499
        # if other is a field, make sure that the domains match
500
        if isinstance(other, Field):
Martin Reinecke's avatar
Martin Reinecke committed
501
            if other._domain != self._domain:
502
                raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
Martin Reinecke committed
503
            tval = getattr(self.val, op)(other.val)
Martin Reinecke's avatar
Martin Reinecke committed
504
            return self if tval is self.val else Field(self._domain, tval)
csongor's avatar
csongor committed
505

506
507
        if np.isscalar(other) or isinstance(other, dobj.data_object):
            tval = getattr(self.val, op)(other)
Martin Reinecke's avatar
Martin Reinecke committed
508
            return self if tval is self.val else Field(self._domain, tval)
509

Martin Reinecke's avatar
Martin Reinecke committed
510
        return NotImplemented
csongor's avatar
csongor committed
511
512

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

515
    def __radd__(self, other):
Theo Steininger's avatar
Theo Steininger committed
516
        return self._binary_helper(other, op='__radd__')
csongor's avatar
csongor committed
517
518

    def __iadd__(self, other):
519
        return self._binary_helper(other, op='__iadd__')
csongor's avatar
csongor committed
520
521

    def __sub__(self, other):
Theo Steininger's avatar
Theo Steininger committed
522
        return self._binary_helper(other, op='__sub__')
csongor's avatar
csongor committed
523
524

    def __rsub__(self, other):
Theo Steininger's avatar
Theo Steininger committed
525
        return self._binary_helper(other, op='__rsub__')
csongor's avatar
csongor committed
526
527

    def __isub__(self, other):
528
        return self._binary_helper(other, op='__isub__')
csongor's avatar
csongor committed
529
530

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

533
    def __rmul__(self, other):
Theo Steininger's avatar
Theo Steininger committed
534
        return self._binary_helper(other, op='__rmul__')
csongor's avatar
csongor committed
535
536

    def __imul__(self, other):
537
        return self._binary_helper(other, op='__imul__')
csongor's avatar
csongor committed
538
539

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

Martin Reinecke's avatar
Martin Reinecke committed
542
543
544
    def __truediv__(self, other):
        return self._binary_helper(other, op='__truediv__')

csongor's avatar
csongor committed
545
    def __rdiv__(self, other):
Theo Steininger's avatar
Theo Steininger committed
546
        return self._binary_helper(other, op='__rdiv__')
csongor's avatar
csongor committed
547

Martin Reinecke's avatar
Martin Reinecke committed
548
549
550
    def __rtruediv__(self, other):
        return self._binary_helper(other, op='__rtruediv__')

csongor's avatar
csongor committed
551
    def __idiv__(self, other):
552
        return self._binary_helper(other, op='__idiv__')
553

csongor's avatar
csongor committed
554
    def __pow__(self, other):
Theo Steininger's avatar
Theo Steininger committed
555
        return self._binary_helper(other, op='__pow__')
csongor's avatar
csongor committed
556
557

    def __rpow__(self, other):
Theo Steininger's avatar
Theo Steininger committed
558
        return self._binary_helper(other, op='__rpow__')
csongor's avatar
csongor committed
559
560

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

Theo Steininger's avatar
Theo Steininger committed
563
    def __repr__(self):
Martin Reinecke's avatar
Martin Reinecke committed
564
        return "<nifty4.Field>"
Theo Steininger's avatar
Theo Steininger committed
565
566

    def __str__(self):
Martin Reinecke's avatar
Martin Reinecke committed
567
        return "nifty4.Field instance\n- domain      = " + \
568
               self._domain.__str__() + \
569
               "\n- val         = " + repr(self.val)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
570
571
572
573
574
575
576
577


# 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
578
        if not isinstance(out, Field) or x._domain != out._domain:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
579
580
581
582
            raise ValueError("Bad 'out' argument")
        function(x.val, out=out.val)
        return out
    else:
Martin Reinecke's avatar
Martin Reinecke committed
583
        return Field(domain=x._domain, val=function(x.val))
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
584
585
586
587
588
589
590
591
592
593
594
595
596
597


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
598
599
600
601
def tanh(x, out=None):
    return _math_helper(x, dobj.tanh, out)


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