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

Martin Reinecke's avatar
Martin Reinecke committed
126
    @staticmethod
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
127
    def _infer_domain(domain, val=None):
128
        if domain is None:
129
            if isinstance(val, Field):
Martin Reinecke's avatar
Martin Reinecke committed
130
                return val._domain
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
131
            if np.isscalar(val):
132
                return DomainTuple.make(())  # empty domain tuple
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
133
            raise TypeError("could not infer domain from value")
Martin Reinecke's avatar
Martin Reinecke committed
134
        return DomainTuple.make(domain)
135

Martin Reinecke's avatar
Martin Reinecke committed
136
137
    @staticmethod
    def _infer_dtype(dtype, val):
Martin Reinecke's avatar
Martin Reinecke committed
138
139
140
141
        if dtype is not None:
            return dtype
        if val is None:
            raise ValueError("could not infer dtype")
Martin Reinecke's avatar
Martin Reinecke committed
142
143
        if isinstance(val, Field):
            return val.dtype
Martin Reinecke's avatar
Martin Reinecke committed
144
        return np.result_type(val)
145

Martin Reinecke's avatar
Martin Reinecke committed
146
147
    @staticmethod
    def from_random(random_type, domain, dtype=np.float64, **kwargs):
148
149
150
151
        """ Draws a random field with the given parameters.

        Parameters
        ----------
Martin Reinecke's avatar
Martin Reinecke committed
152
        random_type : str
153
154
            'pm1', 'normal', 'uniform' are the supported arguments for this
            method.
Theo Steininger's avatar
Theo Steininger committed
155

Martin Reinecke's avatar
Martin Reinecke committed
156
        domain : DomainTuple
157
            The domain of the output random field
Theo Steininger's avatar
Theo Steininger committed
158

159
160
        dtype : type
            The datatype of the output random field
Theo Steininger's avatar
Theo Steininger committed
161

162
163
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
164
        Field
165
166
            The output object.
        """
Martin Reinecke's avatar
Martin Reinecke committed
167
        domain = DomainTuple.make(domain)
Martin Reinecke's avatar
Martin Reinecke committed
168
169
170
        return Field(domain=domain,
                     val=dobj.from_random(random_type, dtype=dtype,
                                          shape=domain.shape, **kwargs))
171

Martin Reinecke's avatar
Martin Reinecke committed
172
173
174
    def fill(self, fill_value):
        self._val.fill(fill_value)

Martin Reinecke's avatar
Martin Reinecke committed
175
176
177
178
179
180
181
    def lock(self):
        dobj.lock(self._val)

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

Theo Steininger's avatar
Theo Steininger committed
182
183
    @property
    def val(self):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
184
        """ Returns the data object associated with this Field.
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
185
        No copy is made.
186
        """
Martin Reinecke's avatar
Martin Reinecke committed
187
        return self._val
csongor's avatar
csongor committed
188

Martin Reinecke's avatar
Martin Reinecke committed
189
190
191
192
    @property
    def dtype(self):
        return self._val.dtype

Martin Reinecke's avatar
Martin Reinecke committed
193
194
195
196
    @property
    def domain(self):
        return self._domain

197
198
    @property
    def shape(self):
Theo Steininger's avatar
Theo Steininger committed
199
        """ Returns the total shape of the Field's data array.
Theo Steininger's avatar
Theo Steininger committed
200

201
202
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
203
204
        tuple of int
            the dimensions of the spaces in domain.
Martin Reinecke's avatar
Martin Reinecke committed
205
        """
Martin Reinecke's avatar
Martin Reinecke committed
206
        return self._domain.shape
csongor's avatar
csongor committed
207

208
    @property
Martin Reinecke's avatar
Martin Reinecke committed
209
    def size(self):
Theo Steininger's avatar
Theo Steininger committed
210
        """ Returns the total number of pixel-dimensions the field has.
Theo Steininger's avatar
Theo Steininger committed
211

Theo Steininger's avatar
Theo Steininger committed
212
        Effectively, all values from shape are multiplied.
Theo Steininger's avatar
Theo Steininger committed
213

214
215
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
216
217
        int
            the dimension of the Field.
218
        """
Martin Reinecke's avatar
Martin Reinecke committed
219
        return self._domain.size
csongor's avatar
csongor committed
220

Theo Steininger's avatar
Theo Steininger committed
221
222
    @property
    def real(self):
Martin Reinecke's avatar
Martin Reinecke committed
223
        """ The real part of the field (data is not copied)."""
224
        if not np.issubdtype(self.dtype, np.complexfloating):
225
            return self
Martin Reinecke's avatar
Martin Reinecke committed
226
        return Field(self._domain, self.val.real)
Theo Steininger's avatar
Theo Steininger committed
227
228
229

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

Martin Reinecke's avatar
Martin Reinecke committed
235
    def copy(self):
236
        """ Returns a full copy of the Field.
Theo Steininger's avatar
Theo Steininger committed
237

Martin Reinecke's avatar
Martin Reinecke committed
238
        The returned object will be an identical copy of the original Field.
Theo Steininger's avatar
Theo Steininger committed
239

240
241
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
242
        Field
Martin Reinecke's avatar
Martin Reinecke committed
243
            An identical copy of 'self'.
244
        """
Martin Reinecke's avatar
Martin Reinecke committed
245
        return Field(val=self, copy=True)
csongor's avatar
csongor committed
246

Martin Reinecke's avatar
Martin Reinecke committed
247
248
249
250
251
252
253
254
255
256
257
258
259
    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
260
        return Field(val=self, copy=True, locked=True)
Martin Reinecke's avatar
Martin Reinecke committed
261

262
263
    def scalar_weight(self, spaces=None):
        if np.isscalar(spaces):
Martin Reinecke's avatar
Martin Reinecke committed
264
            return self._domain[spaces].scalar_dvol
265
266

        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
267
            spaces = range(len(self._domain))
Martin Reinecke's avatar
Martin Reinecke committed
268
        res = 1.
269
        for i in spaces:
Martin Reinecke's avatar
Martin Reinecke committed
270
            tmp = self._domain[i].scalar_dvol
271
272
273
274
275
            if tmp is None:
                return None
            res *= tmp
        return res

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

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

287
    def weight(self, power=1, spaces=None, out=None):
Theo Steininger's avatar
Theo Steininger committed
288
        """ Weights the pixels of `self` with their invidual pixel-volume.
289
290
291
292

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

Martin Reinecke's avatar
Martin Reinecke committed
295
        spaces : int or tuple of int
Theo Steininger's avatar
Theo Steininger committed
296
            Determines on which subspace the operation takes place.
Theo Steininger's avatar
Theo Steininger committed
297

298
299
300
301
302
        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"!

303
304
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
305
        Field
Theo Steininger's avatar
Theo Steininger committed
306
            The weighted field.
307
        """
308
309
310
311
312
        if out is None:
            out = self.copy()
        else:
            if out is not self:
                out.copy_content_from(self)
csongor's avatar
csongor committed
313

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

316
317
        fct = 1.
        for ind in spaces:
Martin Reinecke's avatar
Martin Reinecke committed
318
            wgt = self._domain[ind].dvol
319
320
321
            if np.isscalar(wgt):
                fct *= wgt
            else:
Martin Reinecke's avatar
Martin Reinecke committed
322
                new_shape = np.ones(len(self.shape), dtype=np.int)
Martin Reinecke's avatar
Martin Reinecke committed
323
324
                new_shape[self._domain.axes[ind][0]:
                          self._domain.axes[ind][-1]+1] = wgt.shape
325
                wgt = wgt.reshape(new_shape)
Martin Reinecke's avatar
Martin Reinecke committed
326
327
                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
328
                    wgt = dobj.local_data(dobj.from_global_data(wgt))
329
330
                lout = dobj.local_data(out.val)
                lout *= wgt**power
331
        fct = fct**power
Martin Reinecke's avatar
Martin Reinecke committed
332
        if fct != 1.:
333
            out *= fct
334

335
        return out
csongor's avatar
csongor committed
336

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

340
341
342
        Parameters
        ----------
        x : Field
343
            x must live on the same domain as `self`.
Theo Steininger's avatar
Theo Steininger committed
344

Martin Reinecke's avatar
Martin Reinecke committed
345
        spaces : None, int or tuple of int (default: None)
346
347
            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
348

349
350
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
351
        float, complex, either scalar (for full dot products)
352
                              or Field (for partial dot products)
353
        """
354
355
356
        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
357

Martin Reinecke's avatar
Martin Reinecke committed
358
        if x._domain != self._domain:
359
            raise ValueError("Domain mismatch")
Theo Steininger's avatar
Theo Steininger committed
360

Martin Reinecke's avatar
Martin Reinecke committed
361
        ndom = len(self._domain)
362
363
364
        spaces = utilities.parse_spaces(spaces, ndom)

        if len(spaces) == ndom:
Martin Reinecke's avatar
Martin Reinecke committed
365
            return dobj.vdot(self.val, x.val)
366
367
        # 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
368
        return (self.conjugate()*x).sum(spaces=spaces)
Theo Steininger's avatar
Theo Steininger committed
369

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

Theo Steininger's avatar
Theo Steininger committed
373
374
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
375
        float
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
376
            The L2-norm of the field values.
csongor's avatar
csongor committed
377
        """
378
        return np.sqrt(np.abs(self.vdot(x=self)))
csongor's avatar
csongor committed
379

Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
380
    def conjugate(self):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
381
        """ Returns the complex conjugate of the field.
Theo Steininger's avatar
Theo Steininger committed
382

383
384
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
385
386
        Field
            The complex conjugated field.
csongor's avatar
csongor committed
387
        """
Martin Reinecke's avatar
Martin Reinecke committed
388
        return Field(self._domain, self.val.conjugate())
csongor's avatar
csongor committed
389

Theo Steininger's avatar
Theo Steininger committed
390
    # ---General unary/contraction methods---
391

Theo Steininger's avatar
Theo Steininger committed
392
393
    def __pos__(self):
        return self.copy()
394

Theo Steininger's avatar
Theo Steininger committed
395
    def __neg__(self):
396
        return Field(self._domain, -self.val)
csongor's avatar
csongor committed
397

Theo Steininger's avatar
Theo Steininger committed
398
    def __abs__(self):
399
        return Field(self._domain, dobj.abs(self.val))
csongor's avatar
csongor committed
400

401
    def _contraction_helper(self, op, spaces):
Theo Steininger's avatar
Theo Steininger committed
402
        if spaces is None:
403
            return getattr(self.val, op)()
404

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

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

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

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

Theo Steininger's avatar
Theo Steininger committed
415
416
417
        # 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
418
        else:
Martin Reinecke's avatar
Martin Reinecke committed
419
            return_domain = tuple(dom
Martin Reinecke's avatar
Martin Reinecke committed
420
                                  for i, dom in enumerate(self._domain)
Theo Steininger's avatar
Theo Steininger committed
421
                                  if i not in spaces)
422

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

425
426
    def sum(self, spaces=None):
        return self._contraction_helper('sum', spaces)
csongor's avatar
csongor committed
427

428
    def integrate(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
429
430
431
432
433
        swgt = self.scalar_weight(spaces)
        if swgt is not None:
            res = self.sum(spaces)
            res *= swgt
            return res
434
435
436
        tmp = self.weight(1, spaces=spaces)
        return tmp.sum(spaces)

437
438
    def prod(self, spaces=None):
        return self._contraction_helper('prod', spaces)
csongor's avatar
csongor committed
439

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

443
444
    def any(self, spaces=None):
        return self._contraction_helper('any', spaces)
csongor's avatar
csongor committed
445

446
447
    def min(self, spaces=None):
        return self._contraction_helper('min', spaces)
csongor's avatar
csongor committed
448

449
450
    def max(self, spaces=None):
        return self._contraction_helper('max', spaces)
csongor's avatar
csongor committed
451

452
    def mean(self, spaces=None):
453
454
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('mean', spaces)
Martin Reinecke's avatar
Martin Reinecke committed
455
456
457
        # MR FIXME: not very efficient
        tmp = self.weight(1)
        return tmp.sum(spaces)*(1./tmp.total_volume(spaces))
csongor's avatar
csongor committed
458

459
    def var(self, spaces=None):
460
461
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('var', spaces)
Martin Reinecke's avatar
Martin Reinecke committed
462
463
464
465
466
467
468
469
470
        # 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
471

472
    def std(self, spaces=None):
473
474
475
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('std', spaces)
        return sqrt(self.var(spaces))
csongor's avatar
csongor committed
476

477
478
479
    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
480
        if other._domain != self._domain:
481
            raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
482
        dobj.local_data(self.val)[()] = dobj.local_data(other.val)[()]
483

484
    def _binary_helper(self, other, op):
csongor's avatar
csongor committed
485
        # if other is a field, make sure that the domains match
486
        if isinstance(other, Field):
Martin Reinecke's avatar
Martin Reinecke committed
487
            if other._domain != self._domain:
488
                raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
Martin Reinecke committed
489
            tval = getattr(self.val, op)(other.val)
Martin Reinecke's avatar
Martin Reinecke committed
490
            return self if tval is self.val else Field(self._domain, tval)
csongor's avatar
csongor committed
491

492
493
        if np.isscalar(other) or isinstance(other, dobj.data_object):
            tval = getattr(self.val, op)(other)
Martin Reinecke's avatar
Martin Reinecke committed
494
            return self if tval is self.val else Field(self._domain, tval)
495

Martin Reinecke's avatar
Martin Reinecke committed
496
        return NotImplemented
csongor's avatar
csongor committed
497
498

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

501
    def __radd__(self, other):
Theo Steininger's avatar
Theo Steininger committed
502
        return self._binary_helper(other, op='__radd__')
csongor's avatar
csongor committed
503
504

    def __iadd__(self, other):
505
        return self._binary_helper(other, op='__iadd__')
csongor's avatar
csongor committed
506
507

    def __sub__(self, other):
Theo Steininger's avatar
Theo Steininger committed
508
        return self._binary_helper(other, op='__sub__')
csongor's avatar
csongor committed
509
510

    def __rsub__(self, other):
Theo Steininger's avatar
Theo Steininger committed
511
        return self._binary_helper(other, op='__rsub__')
csongor's avatar
csongor committed
512
513

    def __isub__(self, other):
514
        return self._binary_helper(other, op='__isub__')
csongor's avatar
csongor committed
515
516

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

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

    def __imul__(self, other):
523
        return self._binary_helper(other, op='__imul__')
csongor's avatar
csongor committed
524
525

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

Martin Reinecke's avatar
Martin Reinecke committed
528
529
530
    def __truediv__(self, other):
        return self._binary_helper(other, op='__truediv__')

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

Martin Reinecke's avatar
Martin Reinecke committed
534
535
536
    def __rtruediv__(self, other):
        return self._binary_helper(other, op='__rtruediv__')

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

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

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

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

Theo Steininger's avatar
Theo Steininger committed
549
    def __repr__(self):
Martin Reinecke's avatar
Martin Reinecke committed
550
        return "<nifty4.Field>"
Theo Steininger's avatar
Theo Steininger committed
551
552

    def __str__(self):
Martin Reinecke's avatar
Martin Reinecke committed
553
        return "nifty4.Field instance\n- domain      = " + \
554
               self._domain.__str__() + \
555
               "\n- val         = " + repr(self.val)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
556
557
558
559
560
561
562
563


# 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
564
        if not isinstance(out, Field) or x._domain != out._domain:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
565
566
567
568
            raise ValueError("Bad 'out' argument")
        function(x.val, out=out.val)
        return out
    else:
Martin Reinecke's avatar
Martin Reinecke committed
569
        return Field(domain=x._domain, val=function(x.val))
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
570
571
572
573
574
575
576
577
578
579
580
581
582
583


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
584
585
586
587
def tanh(x, out=None):
    return _math_helper(x, dobj.tanh, out)


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