field.py 18.2 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
    @staticmethod
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
134
    def _infer_domain(domain, val=None):
135
        if domain is None:
136
            if isinstance(val, Field):
Martin Reinecke's avatar
Martin Reinecke committed
137
                return val._domain
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
138
            if np.isscalar(val):
139
                return DomainTuple.make(())  # empty domain tuple
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
140
            raise TypeError("could not infer domain from value")
Martin Reinecke's avatar
Martin Reinecke committed
141
        return DomainTuple.make(domain)
142

Martin Reinecke's avatar
Martin Reinecke committed
143
144
    @staticmethod
    def _infer_dtype(dtype, val):
Martin Reinecke's avatar
Martin Reinecke committed
145
146
147
148
        if dtype is not None:
            return dtype
        if val is None:
            raise ValueError("could not infer dtype")
Martin Reinecke's avatar
Martin Reinecke committed
149
150
        if isinstance(val, Field):
            return val.dtype
Martin Reinecke's avatar
Martin Reinecke committed
151
        return np.result_type(val)
152

Martin Reinecke's avatar
Martin Reinecke committed
153
154
    @staticmethod
    def from_random(random_type, domain, dtype=np.float64, **kwargs):
155
156
157
158
        """ Draws a random field with the given parameters.

        Parameters
        ----------
Martin Reinecke's avatar
Martin Reinecke committed
159
        random_type : str
160
161
            'pm1', 'normal', 'uniform' are the supported arguments for this
            method.
Theo Steininger's avatar
Theo Steininger committed
162

Martin Reinecke's avatar
Martin Reinecke committed
163
        domain : DomainTuple
164
            The domain of the output random field
Theo Steininger's avatar
Theo Steininger committed
165

166
167
        dtype : type
            The datatype of the output random field
Theo Steininger's avatar
Theo Steininger committed
168

169
170
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
171
        Field
172
173
            The output object.
        """
Martin Reinecke's avatar
Martin Reinecke committed
174
        domain = DomainTuple.make(domain)
Martin Reinecke's avatar
Martin Reinecke committed
175
176
177
        return Field(domain=domain,
                     val=dobj.from_random(random_type, dtype=dtype,
                                          shape=domain.shape, **kwargs))
178

Martin Reinecke's avatar
Martin Reinecke committed
179
180
181
    def fill(self, fill_value):
        self._val.fill(fill_value)

Martin Reinecke's avatar
Martin Reinecke committed
182
183
    def lock(self):
        dobj.lock(self._val)
184
        return self
Martin Reinecke's avatar
Martin Reinecke committed
185
186
187
188
189

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

Theo Steininger's avatar
Theo Steininger committed
190
191
    @property
    def val(self):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
192
        """ Returns the data object associated with this Field.
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
193
        No copy is made.
194
        """
Martin Reinecke's avatar
Martin Reinecke committed
195
        return self._val
csongor's avatar
csongor committed
196

Martin Reinecke's avatar
Martin Reinecke committed
197
198
199
200
    @property
    def dtype(self):
        return self._val.dtype

Martin Reinecke's avatar
Martin Reinecke committed
201
202
203
204
    @property
    def domain(self):
        return self._domain

205
206
    @property
    def shape(self):
Theo Steininger's avatar
Theo Steininger committed
207
        """ Returns the total shape of the Field's data array.
Theo Steininger's avatar
Theo Steininger committed
208

209
210
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
211
212
        tuple of int
            the dimensions of the spaces in domain.
Martin Reinecke's avatar
Martin Reinecke committed
213
        """
Martin Reinecke's avatar
Martin Reinecke committed
214
        return self._domain.shape
csongor's avatar
csongor committed
215

216
    @property
Martin Reinecke's avatar
Martin Reinecke committed
217
    def size(self):
Theo Steininger's avatar
Theo Steininger committed
218
        """ Returns the total number of pixel-dimensions the field has.
Theo Steininger's avatar
Theo Steininger committed
219

Theo Steininger's avatar
Theo Steininger committed
220
        Effectively, all values from shape are multiplied.
Theo Steininger's avatar
Theo Steininger committed
221

222
223
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
224
225
        int
            the dimension of the Field.
226
        """
Martin Reinecke's avatar
Martin Reinecke committed
227
        return self._domain.size
csongor's avatar
csongor committed
228

Theo Steininger's avatar
Theo Steininger committed
229
230
    @property
    def real(self):
Martin Reinecke's avatar
Martin Reinecke committed
231
        """ The real part of the field (data is not copied)."""
232
        if not np.issubdtype(self.dtype, np.complexfloating):
233
            return self
Martin Reinecke's avatar
Martin Reinecke committed
234
        return Field(self._domain, self.val.real)
Theo Steininger's avatar
Theo Steininger committed
235
236
237

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

Martin Reinecke's avatar
Martin Reinecke committed
243
    def copy(self):
244
        """ Returns a full copy of the Field.
Theo Steininger's avatar
Theo Steininger committed
245

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

248
249
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
250
        Field
Martin Reinecke's avatar
Martin Reinecke committed
251
            An identical copy of 'self'.
252
        """
Martin Reinecke's avatar
Martin Reinecke committed
253
        return Field(val=self, copy=True)
csongor's avatar
csongor committed
254

Martin Reinecke's avatar
Martin Reinecke committed
255
256
257
258
259
260
261
262
263
264
265
266
267
    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
268
        return Field(val=self, copy=True, locked=True)
Martin Reinecke's avatar
Martin Reinecke committed
269

270
271
    def scalar_weight(self, spaces=None):
        if np.isscalar(spaces):
Martin Reinecke's avatar
Martin Reinecke committed
272
            return self._domain[spaces].scalar_dvol
273
274

        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
275
            spaces = range(len(self._domain))
Martin Reinecke's avatar
Martin Reinecke committed
276
        res = 1.
277
        for i in spaces:
Martin Reinecke's avatar
Martin Reinecke committed
278
            tmp = self._domain[i].scalar_dvol
279
280
281
282
283
            if tmp is None:
                return None
            res *= tmp
        return res

Martin Reinecke's avatar
Martin Reinecke committed
284
285
    def total_volume(self, spaces=None):
        if np.isscalar(spaces):
Martin Reinecke's avatar
Martin Reinecke committed
286
            return self._domain[spaces].total_volume
Martin Reinecke's avatar
Martin Reinecke committed
287
288
289
290
291

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

295
    def weight(self, power=1, spaces=None, out=None):
Theo Steininger's avatar
Theo Steininger committed
296
        """ Weights the pixels of `self` with their invidual pixel-volume.
297
298
299
300

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

Martin Reinecke's avatar
Martin Reinecke committed
303
        spaces : int or tuple of int
Theo Steininger's avatar
Theo Steininger committed
304
            Determines on which subspace the operation takes place.
Theo Steininger's avatar
Theo Steininger committed
305

306
307
308
309
310
        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"!

311
312
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
313
        Field
Theo Steininger's avatar
Theo Steininger committed
314
            The weighted field.
315
        """
316
317
318
319
320
        if out is None:
            out = self.copy()
        else:
            if out is not self:
                out.copy_content_from(self)
csongor's avatar
csongor committed
321

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

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

343
        return out
csongor's avatar
csongor committed
344

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

348
349
350
        Parameters
        ----------
        x : Field
351
            x must live on the same domain as `self`.
Theo Steininger's avatar
Theo Steininger committed
352

Martin Reinecke's avatar
Martin Reinecke committed
353
        spaces : None, int or tuple of int (default: None)
354
355
            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
356

357
358
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
359
        float, complex, either scalar (for full dot products)
360
                              or Field (for partial dot products)
361
        """
362
363
364
        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
365

Martin Reinecke's avatar
Martin Reinecke committed
366
        if x._domain != self._domain:
367
            raise ValueError("Domain mismatch")
Theo Steininger's avatar
Theo Steininger committed
368

Martin Reinecke's avatar
Martin Reinecke committed
369
        ndom = len(self._domain)
370
371
372
        spaces = utilities.parse_spaces(spaces, ndom)

        if len(spaces) == ndom:
Martin Reinecke's avatar
Martin Reinecke committed
373
            return dobj.vdot(self.val, x.val)
374
375
        # 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
376
        return (self.conjugate()*x).sum(spaces=spaces)
Theo Steininger's avatar
Theo Steininger committed
377

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

Theo Steininger's avatar
Theo Steininger committed
381
382
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
383
        float
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
384
            The L2-norm of the field values.
csongor's avatar
csongor committed
385
        """
386
        return np.sqrt(np.abs(self.vdot(x=self)))
csongor's avatar
csongor committed
387

Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
388
    def conjugate(self):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
389
        """ Returns the complex conjugate of the field.
Theo Steininger's avatar
Theo Steininger committed
390

391
392
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
393
394
        Field
            The complex conjugated field.
csongor's avatar
csongor committed
395
        """
Martin Reinecke's avatar
Martin Reinecke committed
396
        return Field(self._domain, self.val.conjugate())
csongor's avatar
csongor committed
397

Theo Steininger's avatar
Theo Steininger committed
398
    # ---General unary/contraction methods---
399

Theo Steininger's avatar
Theo Steininger committed
400
401
    def __pos__(self):
        return self.copy()
402

Theo Steininger's avatar
Theo Steininger committed
403
    def __neg__(self):
404
        return Field(self._domain, -self.val)
csongor's avatar
csongor committed
405

Theo Steininger's avatar
Theo Steininger committed
406
    def __abs__(self):
407
        return Field(self._domain, dobj.abs(self.val))
csongor's avatar
csongor committed
408

409
    def _contraction_helper(self, op, spaces):
Theo Steininger's avatar
Theo Steininger committed
410
        if spaces is None:
411
            return getattr(self.val, op)()
412

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

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

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

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

Theo Steininger's avatar
Theo Steininger committed
423
424
425
        # 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
426
        else:
Martin Reinecke's avatar
Martin Reinecke committed
427
            return_domain = tuple(dom
Martin Reinecke's avatar
Martin Reinecke committed
428
                                  for i, dom in enumerate(self._domain)
Theo Steininger's avatar
Theo Steininger committed
429
                                  if i not in spaces)
430

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

433
434
    def sum(self, spaces=None):
        return self._contraction_helper('sum', spaces)
csongor's avatar
csongor committed
435

436
    def integrate(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
437
438
439
440
441
        swgt = self.scalar_weight(spaces)
        if swgt is not None:
            res = self.sum(spaces)
            res *= swgt
            return res
442
443
444
        tmp = self.weight(1, spaces=spaces)
        return tmp.sum(spaces)

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

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

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

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

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

460
    def mean(self, spaces=None):
461
462
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('mean', spaces)
Martin Reinecke's avatar
Martin Reinecke committed
463
464
465
        # MR FIXME: not very efficient
        tmp = self.weight(1)
        return tmp.sum(spaces)*(1./tmp.total_volume(spaces))
csongor's avatar
csongor committed
466

467
    def var(self, spaces=None):
468
469
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('var', spaces)
Martin Reinecke's avatar
Martin Reinecke committed
470
471
472
473
474
475
476
477
478
        # 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
479

480
    def std(self, spaces=None):
481
482
483
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('std', spaces)
        return sqrt(self.var(spaces))
csongor's avatar
csongor committed
484

485
486
487
    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
488
        if other._domain != self._domain:
489
            raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
490
        dobj.local_data(self.val)[()] = dobj.local_data(other.val)[()]
491

492
    def _binary_helper(self, other, op):
csongor's avatar
csongor committed
493
        # if other is a field, make sure that the domains match
494
        if isinstance(other, Field):
Martin Reinecke's avatar
Martin Reinecke committed
495
            if other._domain != self._domain:
496
                raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
Martin Reinecke committed
497
            tval = getattr(self.val, op)(other.val)
Martin Reinecke's avatar
Martin Reinecke committed
498
            return self if tval is self.val else Field(self._domain, tval)
csongor's avatar
csongor committed
499

500
501
        if np.isscalar(other) or isinstance(other, dobj.data_object):
            tval = getattr(self.val, op)(other)
Martin Reinecke's avatar
Martin Reinecke committed
502
            return self if tval is self.val else Field(self._domain, tval)
503

Martin Reinecke's avatar
Martin Reinecke committed
504
        return NotImplemented
csongor's avatar
csongor committed
505
506

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

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

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

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

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

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

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

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

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

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

Martin Reinecke's avatar
Martin Reinecke committed
536
537
538
    def __truediv__(self, other):
        return self._binary_helper(other, op='__truediv__')

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

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

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

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

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

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

Theo Steininger's avatar
Theo Steininger committed
557
    def __repr__(self):
Martin Reinecke's avatar
Martin Reinecke committed
558
        return "<nifty4.Field>"
Theo Steininger's avatar
Theo Steininger committed
559
560

    def __str__(self):
Martin Reinecke's avatar
Martin Reinecke committed
561
        return "nifty4.Field instance\n- domain      = " + \
562
               self._domain.__str__() + \
563
               "\n- val         = " + repr(self.val)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
564
565
566
567
568
569
570
571


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


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
592
593
594
595
def tanh(x, out=None):
    return _math_helper(x, dobj.tanh, out)


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