field.py 19.1 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.
47
    """
48

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

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

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

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

    @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
134
        return Field.zeros(field._domain, dtype)
135
136
137
138
139
140
141

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

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

152
153
154
155
156
157
158
    @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
159
160
161
162
    @property
    def local_data(self):
        return dobj.local_data(self._val)

Martin Reinecke's avatar
Martin Reinecke committed
163
164
165
    def cast_domain(self, new_domain):
        return Field(new_domain, self._val)

Martin Reinecke's avatar
Martin Reinecke committed
166
    @staticmethod
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
167
    def _infer_domain(domain, val=None):
168
        if domain is None:
169
            if isinstance(val, Field):
Martin Reinecke's avatar
Martin Reinecke committed
170
                return val._domain
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
171
            if np.isscalar(val):
172
                return DomainTuple.make(())  # empty domain tuple
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
173
            raise TypeError("could not infer domain from value")
Martin Reinecke's avatar
Martin Reinecke committed
174
        return DomainTuple.make(domain)
175

Martin Reinecke's avatar
Martin Reinecke committed
176
177
    @staticmethod
    def _infer_dtype(dtype, val):
Martin Reinecke's avatar
Martin Reinecke committed
178
179
180
181
        if dtype is not None:
            return dtype
        if val is None:
            raise ValueError("could not infer dtype")
Martin Reinecke's avatar
Martin Reinecke committed
182
183
        if isinstance(val, Field):
            return val.dtype
Martin Reinecke's avatar
Martin Reinecke committed
184
        return np.result_type(val)
185

Martin Reinecke's avatar
Martin Reinecke committed
186
187
    @staticmethod
    def from_random(random_type, domain, dtype=np.float64, **kwargs):
188
189
190
191
        """ Draws a random field with the given parameters.

        Parameters
        ----------
Martin Reinecke's avatar
Martin Reinecke committed
192
        random_type : str
193
194
            'pm1', 'normal', 'uniform' are the supported arguments for this
            method.
Theo Steininger's avatar
Theo Steininger committed
195

Martin Reinecke's avatar
Martin Reinecke committed
196
        domain : DomainTuple
197
            The domain of the output random field
Theo Steininger's avatar
Theo Steininger committed
198

199
200
        dtype : type
            The datatype of the output random field
Theo Steininger's avatar
Theo Steininger committed
201

202
203
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
204
        Field
205
206
            The output object.
        """
Martin Reinecke's avatar
Martin Reinecke committed
207
        domain = DomainTuple.make(domain)
Martin Reinecke's avatar
Martin Reinecke committed
208
209
210
        return Field(domain=domain,
                     val=dobj.from_random(random_type, dtype=dtype,
                                          shape=domain.shape, **kwargs))
211

Martin Reinecke's avatar
Martin Reinecke committed
212
213
214
    def fill(self, fill_value):
        self._val.fill(fill_value)

Martin Reinecke's avatar
Martin Reinecke committed
215
216
    def lock(self):
        dobj.lock(self._val)
217
        return self
Martin Reinecke's avatar
Martin Reinecke committed
218
219
220
221
222

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

Theo Steininger's avatar
Theo Steininger committed
223
224
    @property
    def val(self):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
225
        """ Returns the data object associated with this Field.
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
226
        No copy is made.
227
        """
Martin Reinecke's avatar
Martin Reinecke committed
228
        return self._val
csongor's avatar
csongor committed
229

Martin Reinecke's avatar
Martin Reinecke committed
230
231
232
233
    @property
    def dtype(self):
        return self._val.dtype

Martin Reinecke's avatar
Martin Reinecke committed
234
235
236
237
    @property
    def domain(self):
        return self._domain

238
239
    @property
    def shape(self):
Theo Steininger's avatar
Theo Steininger committed
240
        """ Returns the total shape of the Field's data array.
Theo Steininger's avatar
Theo Steininger committed
241

242
243
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
244
245
        tuple of int
            the dimensions of the spaces in domain.
Martin Reinecke's avatar
Martin Reinecke committed
246
        """
Martin Reinecke's avatar
Martin Reinecke committed
247
        return self._domain.shape
csongor's avatar
csongor committed
248

249
    @property
Martin Reinecke's avatar
Martin Reinecke committed
250
    def size(self):
Theo Steininger's avatar
Theo Steininger committed
251
        """ Returns the total number of pixel-dimensions the field has.
Theo Steininger's avatar
Theo Steininger committed
252

Theo Steininger's avatar
Theo Steininger committed
253
        Effectively, all values from shape are multiplied.
Theo Steininger's avatar
Theo Steininger committed
254

255
256
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
257
258
        int
            the dimension of the Field.
259
        """
Martin Reinecke's avatar
Martin Reinecke committed
260
        return self._domain.size
csongor's avatar
csongor committed
261

Theo Steininger's avatar
Theo Steininger committed
262
263
    @property
    def real(self):
Martin Reinecke's avatar
Martin Reinecke committed
264
        """ The real part of the field (data is not copied)."""
265
        if not np.issubdtype(self.dtype, np.complexfloating):
266
            return self
Martin Reinecke's avatar
Martin Reinecke committed
267
        return Field(self._domain, self.val.real)
Theo Steininger's avatar
Theo Steininger committed
268
269
270

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

Martin Reinecke's avatar
Martin Reinecke committed
276
    def copy(self):
277
        """ Returns a full copy of the Field.
Theo Steininger's avatar
Theo Steininger committed
278

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

281
282
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
283
        Field
Martin Reinecke's avatar
Martin Reinecke committed
284
            An identical copy of 'self'.
285
        """
Martin Reinecke's avatar
Martin Reinecke committed
286
        return Field(val=self, copy=True)
csongor's avatar
csongor committed
287

Martin Reinecke's avatar
Martin Reinecke committed
288
289
290
291
292
293
294
295
296
297
298
299
300
    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
301
        return Field(val=self, copy=True, locked=True)
Martin Reinecke's avatar
Martin Reinecke committed
302

303
304
    def scalar_weight(self, spaces=None):
        if np.isscalar(spaces):
Martin Reinecke's avatar
Martin Reinecke committed
305
            return self._domain[spaces].scalar_dvol
306
307

        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
308
            spaces = range(len(self._domain))
Martin Reinecke's avatar
Martin Reinecke committed
309
        res = 1.
310
        for i in spaces:
Martin Reinecke's avatar
Martin Reinecke committed
311
            tmp = self._domain[i].scalar_dvol
312
313
314
315
316
            if tmp is None:
                return None
            res *= tmp
        return res

Martin Reinecke's avatar
Martin Reinecke committed
317
318
    def total_volume(self, spaces=None):
        if np.isscalar(spaces):
Martin Reinecke's avatar
Martin Reinecke committed
319
            return self._domain[spaces].total_volume
Martin Reinecke's avatar
Martin Reinecke committed
320
321
322
323
324

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

328
    def weight(self, power=1, spaces=None, out=None):
Theo Steininger's avatar
Theo Steininger committed
329
        """ Weights the pixels of `self` with their invidual pixel-volume.
330
331
332
333

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

Martin Reinecke's avatar
Martin Reinecke committed
336
        spaces : int or tuple of int
Theo Steininger's avatar
Theo Steininger committed
337
            Determines on which subspace the operation takes place.
Theo Steininger's avatar
Theo Steininger committed
338

339
340
341
342
343
        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"!

344
345
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
346
        Field
Theo Steininger's avatar
Theo Steininger committed
347
            The weighted field.
348
        """
349
350
351
352
353
        if out is None:
            out = self.copy()
        else:
            if out is not self:
                out.copy_content_from(self)
csongor's avatar
csongor committed
354

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

357
358
        fct = 1.
        for ind in spaces:
Martin Reinecke's avatar
Martin Reinecke committed
359
            wgt = self._domain[ind].dvol
360
361
362
            if np.isscalar(wgt):
                fct *= wgt
            else:
Martin Reinecke's avatar
Martin Reinecke committed
363
                new_shape = np.ones(len(self.shape), dtype=np.int)
Martin Reinecke's avatar
Martin Reinecke committed
364
365
                new_shape[self._domain.axes[ind][0]:
                          self._domain.axes[ind][-1]+1] = wgt.shape
366
                wgt = wgt.reshape(new_shape)
Martin Reinecke's avatar
Martin Reinecke committed
367
368
                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
369
                    wgt = dobj.local_data(dobj.from_global_data(wgt))
Martin Reinecke's avatar
Martin Reinecke committed
370
                out.local_data[()] *= wgt**power
371
        fct = fct**power
Martin Reinecke's avatar
Martin Reinecke committed
372
        if fct != 1.:
373
            out *= fct
374

375
        return out
csongor's avatar
csongor committed
376

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

380
381
382
        Parameters
        ----------
        x : Field
383
            x must live on the same domain as `self`.
Theo Steininger's avatar
Theo Steininger committed
384

Martin Reinecke's avatar
Martin Reinecke committed
385
        spaces : None, int or tuple of int (default: None)
386
387
            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
388

389
390
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
391
        float, complex, either scalar (for full dot products)
392
                              or Field (for partial dot products)
393
        """
394
395
396
        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
397

Martin Reinecke's avatar
Martin Reinecke committed
398
        if x._domain != self._domain:
399
            raise ValueError("Domain mismatch")
Theo Steininger's avatar
Theo Steininger committed
400

Martin Reinecke's avatar
Martin Reinecke committed
401
        ndom = len(self._domain)
402
403
404
        spaces = utilities.parse_spaces(spaces, ndom)

        if len(spaces) == ndom:
Martin Reinecke's avatar
Martin Reinecke committed
405
            return dobj.vdot(self.val, x.val)
406
407
        # 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
408
        return (self.conjugate()*x).sum(spaces=spaces)
Theo Steininger's avatar
Theo Steininger committed
409

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

Theo Steininger's avatar
Theo Steininger committed
413
414
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
415
        float
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
416
            The L2-norm of the field values.
csongor's avatar
csongor committed
417
        """
418
        return np.sqrt(np.abs(self.vdot(x=self)))
csongor's avatar
csongor committed
419

Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
420
    def conjugate(self):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
421
        """ Returns the complex conjugate of the field.
Theo Steininger's avatar
Theo Steininger committed
422

423
424
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
425
426
        Field
            The complex conjugated field.
csongor's avatar
csongor committed
427
        """
Martin Reinecke's avatar
Martin Reinecke committed
428
        return Field(self._domain, self.val.conjugate())
csongor's avatar
csongor committed
429

Theo Steininger's avatar
Theo Steininger committed
430
    # ---General unary/contraction methods---
431

Theo Steininger's avatar
Theo Steininger committed
432
433
    def __pos__(self):
        return self.copy()
434

Theo Steininger's avatar
Theo Steininger committed
435
    def __neg__(self):
436
        return Field(self._domain, -self.val)
csongor's avatar
csongor committed
437

Theo Steininger's avatar
Theo Steininger committed
438
    def __abs__(self):
439
        return Field(self._domain, dobj.abs(self.val))
csongor's avatar
csongor committed
440

441
    def _contraction_helper(self, op, spaces):
Theo Steininger's avatar
Theo Steininger committed
442
        if spaces is None:
443
            return getattr(self.val, op)()
444

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

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

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

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

Theo Steininger's avatar
Theo Steininger committed
455
456
457
        # 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
458
        else:
Martin Reinecke's avatar
Martin Reinecke committed
459
            return_domain = tuple(dom
Martin Reinecke's avatar
Martin Reinecke committed
460
                                  for i, dom in enumerate(self._domain)
Theo Steininger's avatar
Theo Steininger committed
461
                                  if i not in spaces)
462

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

465
466
    def sum(self, spaces=None):
        return self._contraction_helper('sum', spaces)
csongor's avatar
csongor committed
467

468
    def integrate(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
469
470
471
472
473
        swgt = self.scalar_weight(spaces)
        if swgt is not None:
            res = self.sum(spaces)
            res *= swgt
            return res
474
475
476
        tmp = self.weight(1, spaces=spaces)
        return tmp.sum(spaces)

477
478
    def prod(self, spaces=None):
        return self._contraction_helper('prod', spaces)
csongor's avatar
csongor committed
479

480
481
    def all(self, spaces=None):
        return self._contraction_helper('all', spaces)
csongor's avatar
csongor committed
482

483
484
    def any(self, spaces=None):
        return self._contraction_helper('any', spaces)
csongor's avatar
csongor committed
485

486
487
    def min(self, spaces=None):
        return self._contraction_helper('min', spaces)
csongor's avatar
csongor committed
488

489
490
    def max(self, spaces=None):
        return self._contraction_helper('max', spaces)
csongor's avatar
csongor committed
491

492
    def mean(self, spaces=None):
493
494
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('mean', spaces)
Martin Reinecke's avatar
Martin Reinecke committed
495
496
497
        # MR FIXME: not very efficient
        tmp = self.weight(1)
        return tmp.sum(spaces)*(1./tmp.total_volume(spaces))
csongor's avatar
csongor committed
498

499
    def var(self, spaces=None):
500
501
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('var', spaces)
Martin Reinecke's avatar
Martin Reinecke committed
502
503
504
505
506
507
508
509
510
        # 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
511

512
    def std(self, spaces=None):
513
514
515
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('std', spaces)
        return sqrt(self.var(spaces))
csongor's avatar
csongor committed
516

517
518
519
    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
520
        if other._domain != self._domain:
521
            raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
Martin Reinecke committed
522
        self.local_data[()] = other.local_data[()]
523

524
    def _binary_helper(self, other, op):
csongor's avatar
csongor committed
525
        # if other is a field, make sure that the domains match
526
        if isinstance(other, Field):
Martin Reinecke's avatar
Martin Reinecke committed
527
            if other._domain != self._domain:
528
                raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
Martin Reinecke committed
529
            tval = getattr(self.val, op)(other.val)
Martin Reinecke's avatar
Martin Reinecke committed
530
            return self if tval is self.val else Field(self._domain, tval)
csongor's avatar
csongor committed
531

532
533
        if np.isscalar(other) or isinstance(other, dobj.data_object):
            tval = getattr(self.val, op)(other)
Martin Reinecke's avatar
Martin Reinecke committed
534
            return self if tval is self.val else Field(self._domain, tval)
535

Martin Reinecke's avatar
Martin Reinecke committed
536
        return NotImplemented
csongor's avatar
csongor committed
537
538

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

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

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

    def __sub__(self, other):
Theo Steininger's avatar
Theo Steininger committed
548
        return self._binary_helper(other, op='__sub__')
csongor's avatar
csongor committed
549
550

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

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

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

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

    def __imul__(self, other):
563
        return self._binary_helper(other, op='__imul__')
csongor's avatar
csongor committed
564
565

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

Martin Reinecke's avatar
Martin Reinecke committed
568
569
570
    def __truediv__(self, other):
        return self._binary_helper(other, op='__truediv__')

csongor's avatar
csongor committed
571
    def __rdiv__(self, other):
Theo Steininger's avatar
Theo Steininger committed
572
        return self._binary_helper(other, op='__rdiv__')
csongor's avatar
csongor committed
573

Martin Reinecke's avatar
Martin Reinecke committed
574
575
576
    def __rtruediv__(self, other):
        return self._binary_helper(other, op='__rtruediv__')

csongor's avatar
csongor committed
577
    def __idiv__(self, other):
578
        return self._binary_helper(other, op='__idiv__')
579

csongor's avatar
csongor committed
580
    def __pow__(self, other):
Theo Steininger's avatar
Theo Steininger committed
581
        return self._binary_helper(other, op='__pow__')
csongor's avatar
csongor committed
582
583

    def __rpow__(self, other):
Theo Steininger's avatar
Theo Steininger committed
584
        return self._binary_helper(other, op='__rpow__')
csongor's avatar
csongor committed
585
586

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

Theo Steininger's avatar
Theo Steininger committed
589
    def __repr__(self):
Martin Reinecke's avatar
Martin Reinecke committed
590
        return "<nifty4.Field>"
Theo Steininger's avatar
Theo Steininger committed
591
592

    def __str__(self):
Martin Reinecke's avatar
Martin Reinecke committed
593
        return "nifty4.Field instance\n- domain      = " + \
594
               self._domain.__str__() + \
595
               "\n- val         = " + repr(self.val)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
596
597
598
599
600
601
602
603


# 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
604
        if not isinstance(out, Field) or x._domain != out._domain:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
605
606
607
608
            raise ValueError("Bad 'out' argument")
        function(x.val, out=out.val)
        return out
    else:
Martin Reinecke's avatar
Martin Reinecke committed
609
        return Field(domain=x._domain, val=function(x.val))
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
610
611
612
613
614
615
616
617
618
619
620
621
622
623


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
624
625
626
627
def tanh(x, out=None):
    return _math_helper(x, dobj.tanh, out)


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