field.py 24.8 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
import sys
27

Martin Reinecke's avatar
cleanup  
Martin Reinecke committed
28 29
__all__ = ["Field", "sqrt", "exp", "log", "conjugate"]

30

Martin Reinecke's avatar
Martin Reinecke committed
31
class Field(object):
Theo Steininger's avatar
Theo Steininger committed
32 33
    """ The discrete representation of a continuous field over multiple spaces.

Martin Reinecke's avatar
Martin Reinecke committed
34
    In NIFTy, Fields are used to store data arrays and carry all the needed
35
    metainformation (i.e. the domain) for operators to be able to work on them.
Theo Steininger's avatar
Theo Steininger committed
36

37 38
    Parameters
    ----------
Martin Reinecke's avatar
Martin Reinecke committed
39
    domain : None, DomainTuple, tuple of Domain, or Domain
Theo Steininger's avatar
Theo Steininger committed
40

Martin Reinecke's avatar
Martin Reinecke committed
41
    val : None, Field, data_object, or scalar
42
        The values the array should contain after init. A scalar input will
Martin Reinecke's avatar
Martin Reinecke committed
43 44
        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
45

46
    dtype : type
Martin Reinecke's avatar
updates  
Martin Reinecke committed
47
        A numpy.type. Most common are float and complex.
Martin Reinecke's avatar
Martin Reinecke committed
48 49 50 51 52

    Notes
    -----
    If possible, do not invoke the constructor directly, but use one of the
    many convenience functions for Field conatruction!
53
    """
54

Martin Reinecke's avatar
Martin Reinecke committed
55 56
    def __init__(self, domain=None, val=None, dtype=None, copy=False,
                 locked=False):
Martin Reinecke's avatar
Martin Reinecke committed
57
        self._domain = self._infer_domain(domain=domain, val=val)
58

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

Martin Reinecke's avatar
Martin Reinecke committed
79 80 81
        if locked:
            dobj.lock(self._val)

82 83 84 85 86 87 88
    # prevent implicit conversion to bool
    def __nonzero__(self):
        raise TypeError("Field does not support implicit conversion to bool")

    def __bool__(self):
        raise TypeError("Field does not support implicit conversion to bool")

89 90
    @staticmethod
    def full(domain, val, dtype=None):
Martin Reinecke's avatar
Martin Reinecke committed
91 92 93 94 95 96 97 98 99 100 101 102 103 104
        """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
        """
105 106 107 108 109 110 111 112
        if not np.isscalar(val):
            raise TypeError("val must be a scalar")
        return Field(DomainTuple.make(domain), val, dtype)

    @staticmethod
    def empty(domain, dtype=None):
        return Field(DomainTuple.make(domain), None, dtype)

113
    @staticmethod
114
    def from_global_data(domain, arr, sum_up=False):
Martin Reinecke's avatar
Martin Reinecke committed
115 116 117 118 119 120 121 122 123
        """Returns a Field constructed from `domain` and `arr`.

        Parameters
        ----------
        domain : DomainTuple, tuple of Domain, or Domain
            the domain of the new Field
        arr : numpy.ndarray
            The data content to be used for the new Field.
            Its shape must match the shape of `domain`.
124 125 126 127 128
        sum_up : bool, optional
            If True, the contents of `arr` are summed up over all MPI tasks
            (if any), and the sum is used as data content.
            If False, the contens of `arr` are used directly, and must be
            identical on all MPI tasks.
Martin Reinecke's avatar
Martin Reinecke committed
129
        """
130
        return Field(domain, dobj.from_global_data(arr, sum_up))
131

Martin Reinecke's avatar
Martin Reinecke committed
132 133 134 135 136
    @staticmethod
    def from_local_data(domain, arr):
        domain = DomainTuple.make(domain)
        return Field(domain, dobj.from_local_data(domain.shape, arr))

137
    def to_global_data(self):
Martin Reinecke's avatar
Martin Reinecke committed
138 139 140 141 142 143 144 145 146 147 148 149
        """Returns an array containing the full data of the field.

        Returns
        -------
        numpy.ndarray : array containing all field entries.
            Its shape is identical to `self.shape`.

        Notes
        -----
        Do not write to the returned array! Depending on whether MPI is
        active or not, this may or may not change the field's data content.
        """
150 151
        return dobj.to_global_data(self._val)

Martin Reinecke's avatar
Martin Reinecke committed
152 153
    @property
    def local_data(self):
Martin Reinecke's avatar
Martin Reinecke committed
154 155 156 157 158 159 160 161 162 163
        """numpy.ndarray : locally residing field data

        Returns a handle to the part of the array data residing on the local
        task (or to the entore array if MPI is not active).

        Notes
        -----
        If the field is not locked, the array data can be modified.
        Use with care!
        """
Martin Reinecke's avatar
Martin Reinecke committed
164 165
        return dobj.local_data(self._val)

Martin Reinecke's avatar
Martin Reinecke committed
166
    def cast_domain(self, new_domain):
Martin Reinecke's avatar
Martin Reinecke committed
167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183
        """Returns a field with the same data, but a different domain

        Parameters
        ----------
        new_domain : Domain, tuple of Domain, or DomainTuple
            The domain for the returned field. Must be shape-compatible to
            `self`.

        Returns
        -------
        Field
            Field living on `new_domain`, but with the same data as `self`.

        Notes
        -----
        No copy is made. If needed, use an additional copy() invocation.
        """
Martin Reinecke's avatar
Martin Reinecke committed
184 185
        return Field(new_domain, self._val)

Martin Reinecke's avatar
Martin Reinecke committed
186
    @staticmethod
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
187
    def _infer_domain(domain, val=None):
188
        if domain is None:
189
            if isinstance(val, Field):
Martin Reinecke's avatar
Martin Reinecke committed
190
                return val._domain
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
191
            if np.isscalar(val):
192
                return DomainTuple.make(())  # empty domain tuple
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
193
            raise TypeError("could not infer domain from value")
Martin Reinecke's avatar
Martin Reinecke committed
194
        return DomainTuple.make(domain)
195

Martin Reinecke's avatar
Martin Reinecke committed
196 197
    @staticmethod
    def _infer_dtype(dtype, val):
Martin Reinecke's avatar
Martin Reinecke committed
198 199 200 201
        if dtype is not None:
            return dtype
        if val is None:
            raise ValueError("could not infer dtype")
Martin Reinecke's avatar
Martin Reinecke committed
202 203
        if isinstance(val, Field):
            return val.dtype
Martin Reinecke's avatar
Martin Reinecke committed
204
        return np.result_type(val)
205

Martin Reinecke's avatar
Martin Reinecke committed
206 207
    @staticmethod
    def from_random(random_type, domain, dtype=np.float64, **kwargs):
208 209 210 211
        """ Draws a random field with the given parameters.

        Parameters
        ----------
Martin Reinecke's avatar
Martin Reinecke committed
212 213
        random_type : 'pm1', 'normal', or 'uniform'
            The random distribution to use.
Martin Reinecke's avatar
Martin Reinecke committed
214
        domain : DomainTuple
215 216 217
            The domain of the output random field
        dtype : type
            The datatype of the output random field
Theo Steininger's avatar
Theo Steininger committed
218

219 220
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
221
        Field
Martin Reinecke's avatar
Martin Reinecke committed
222
            The newly created Field.
223
        """
Martin Reinecke's avatar
Martin Reinecke committed
224
        domain = DomainTuple.make(domain)
Martin Reinecke's avatar
Martin Reinecke committed
225 226 227
        return Field(domain=domain,
                     val=dobj.from_random(random_type, dtype=dtype,
                                          shape=domain.shape, **kwargs))
228

Martin Reinecke's avatar
Martin Reinecke committed
229
    def fill(self, fill_value):
Martin Reinecke's avatar
Martin Reinecke committed
230 231 232 233 234 235 236
        """Fill `self` uniformly with `fill_value`

        Parameters
        ----------
        fill_value: float or complex or int
            The value to fill the field with.
        """
Martin Reinecke's avatar
Martin Reinecke committed
237
        self._val.fill(fill_value)
238
        return self
Martin Reinecke's avatar
Martin Reinecke committed
239

Martin Reinecke's avatar
Martin Reinecke committed
240
    def lock(self):
Martin Reinecke's avatar
Martin Reinecke committed
241 242 243 244 245 246 247 248 249 250 251 252 253 254 255
        """Write-protect the data content of `self`.

        After this call, it will no longer be possible to change the data
        entries of `self`. This is convenient if, for example, a
        DiagonalOperator wants to ensure that its diagonal cannot be modified
        inadvertently, without making copies.

        Notes
        -----
        This will not only prohibit modifications to the entries of `self`, but
        also to the entries of any other Field or numpy array pointing to the
        same data. If an unlocked instance is needed, use copy().

        The fact that there is no `unlock()` method is deliberate.
        """
Martin Reinecke's avatar
Martin Reinecke committed
256
        dobj.lock(self._val)
257
        return self
Martin Reinecke's avatar
Martin Reinecke committed
258 259 260

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

Theo Steininger's avatar
Theo Steininger committed
264 265
    @property
    def val(self):
Martin Reinecke's avatar
Martin Reinecke committed
266 267 268 269 270 271
        """dobj.data_object : the data object storing the field's entries

        Notes
        -----
        This property is intended for low-level, internal use only. Do not use
        from outside of NIFTy's core; there should be better alternatives.
272
        """
Martin Reinecke's avatar
Martin Reinecke committed
273
        return self._val
csongor's avatar
csongor committed
274

Martin Reinecke's avatar
Martin Reinecke committed
275 276
    @property
    def dtype(self):
Martin Reinecke's avatar
Martin Reinecke committed
277
        """type : the data type of the field's entries"""
Martin Reinecke's avatar
Martin Reinecke committed
278 279
        return self._val.dtype

Martin Reinecke's avatar
Martin Reinecke committed
280 281
    @property
    def domain(self):
Martin Reinecke's avatar
Martin Reinecke committed
282
        """DomainTuple : the field's domain"""
Martin Reinecke's avatar
Martin Reinecke committed
283 284
        return self._domain

285 286
    @property
    def shape(self):
Martin Reinecke's avatar
Martin Reinecke committed
287
        """tuple of int : the concatenated shapes of all sub-domains"""
Martin Reinecke's avatar
Martin Reinecke committed
288
        return self._domain.shape
csongor's avatar
csongor committed
289

290
    @property
Martin Reinecke's avatar
Martin Reinecke committed
291
    def size(self):
Martin Reinecke's avatar
Martin Reinecke committed
292
        """int : total number of pixels in the field"""
Martin Reinecke's avatar
Martin Reinecke committed
293
        return self._domain.size
csongor's avatar
csongor committed
294

Theo Steininger's avatar
Theo Steininger committed
295 296
    @property
    def real(self):
Martin Reinecke's avatar
Martin Reinecke committed
297
        """Field : The real part of the field"""
298
        if not np.issubdtype(self.dtype, np.complexfloating):
299
            return self
Martin Reinecke's avatar
Martin Reinecke committed
300
        return Field(self._domain, self.val.real)
Theo Steininger's avatar
Theo Steininger committed
301 302 303

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

Martin Reinecke's avatar
Martin Reinecke committed
309
    def copy(self):
310
        """ Returns a full copy of the Field.
Theo Steininger's avatar
Theo Steininger committed
311

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

315 316
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
317
        Field
Martin Reinecke's avatar
Martin Reinecke committed
318
            An identical, but unlocked copy of 'self'.
319
        """
Martin Reinecke's avatar
Martin Reinecke committed
320
        return Field(val=self, copy=True)
csongor's avatar
csongor committed
321

322 323 324 325 326 327 328 329 330 331 332
    def empty_copy(self):
        """ Returns a Field with identical domain and data type, but
        uninitialized data.

        Returns
        -------
        Field
            A copy of 'self', with uninitialized data.
        """
        return Field(self._domain, dtype=self.dtype)

Martin Reinecke's avatar
Martin Reinecke committed
333 334 335 336 337 338 339 340 341
    def locked_copy(self):
        """ Returns a read-only version of the Field.

        If `self` is locked, returns `self`. Otherwise returns a locked copy
        of `self`.

        Returns
        -------
        Field
Martin Reinecke's avatar
Martin Reinecke committed
342
            A read-only version of `self`.
Martin Reinecke's avatar
Martin Reinecke committed
343
        """
Martin Reinecke's avatar
Martin Reinecke committed
344
        return self if self.locked else Field(val=self, copy=True, locked=True)
Martin Reinecke's avatar
Martin Reinecke committed
345

346
    def scalar_weight(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
347 348 349 350 351 352 353 354 355 356 357 358 359 360
        """Returns the uniform volume element for a sub-domain of `self`.

        Parameters
        ----------
        spaces : int, tuple of int or None
            indices of the sub-domains of the field's domain to be considered.
            If `None`, the entire domain is used.

        Returns
        -------
        float or None
            if the requested sub-domain has a uniform volume element, it is
            returned. Otherwise, `None` is returned.
        """
361
        if np.isscalar(spaces):
Martin Reinecke's avatar
Martin Reinecke committed
362
            return self._domain[spaces].scalar_dvol
363 364

        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
365
            spaces = range(len(self._domain))
Martin Reinecke's avatar
Martin Reinecke committed
366
        res = 1.
367
        for i in spaces:
Martin Reinecke's avatar
Martin Reinecke committed
368
            tmp = self._domain[i].scalar_dvol
369 370 371 372 373
            if tmp is None:
                return None
            res *= tmp
        return res

Martin Reinecke's avatar
Martin Reinecke committed
374
    def total_volume(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
375 376 377 378 379 380 381 382 383 384 385 386 387
        """Returns the total volume of a sub-domain of `self`.

        Parameters
        ----------
        spaces : int, tuple of int or None
            indices of the sub-domains of the field's domain to be considered.
            If `None`, the entire domain is used.

        Returns
        -------
        float
            the total volume of the requested sub-domain.
        """
Martin Reinecke's avatar
Martin Reinecke committed
388
        if np.isscalar(spaces):
Martin Reinecke's avatar
Martin Reinecke committed
389
            return self._domain[spaces].total_volume
Martin Reinecke's avatar
Martin Reinecke committed
390 391 392 393 394

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

398
    def weight(self, power=1, spaces=None, out=None):
Theo Steininger's avatar
Theo Steininger committed
399
        """ Weights the pixels of `self` with their invidual pixel-volume.
400 401 402 403

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

Martin Reinecke's avatar
Martin Reinecke committed
406 407 408
        spaces : None, int or tuple of int
            Determines on which sub-domain the operation takes place.
            If None, the entire domain is used.
Theo Steininger's avatar
Theo Steininger committed
409

410 411 412 413 414
        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"!

415 416
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
417
        Field
Theo Steininger's avatar
Theo Steininger committed
418
            The weighted field.
419
        """
420 421 422 423 424
        if out is None:
            out = self.copy()
        else:
            if out is not self:
                out.copy_content_from(self)
csongor's avatar
csongor committed
425

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

428 429
        fct = 1.
        for ind in spaces:
Martin Reinecke's avatar
Martin Reinecke committed
430
            wgt = self._domain[ind].dvol
431 432 433
            if np.isscalar(wgt):
                fct *= wgt
            else:
Martin Reinecke's avatar
Martin Reinecke committed
434
                new_shape = np.ones(len(self.shape), dtype=np.int)
Martin Reinecke's avatar
Martin Reinecke committed
435 436
                new_shape[self._domain.axes[ind][0]:
                          self._domain.axes[ind][-1]+1] = wgt.shape
437
                wgt = wgt.reshape(new_shape)
Martin Reinecke's avatar
Martin Reinecke committed
438 439
                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
440
                    wgt = dobj.local_data(dobj.from_global_data(wgt))
Martin Reinecke's avatar
Martin Reinecke committed
441
                out.local_data[()] *= wgt**power
442
        fct = fct**power
Martin Reinecke's avatar
Martin Reinecke committed
443
        if fct != 1.:
444
            out *= fct
445

446
        return out
csongor's avatar
csongor committed
447

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

451 452 453
        Parameters
        ----------
        x : Field
454
            x must live on the same domain as `self`.
Theo Steininger's avatar
Theo Steininger committed
455

Martin Reinecke's avatar
Martin Reinecke committed
456
        spaces : None, int or tuple of int (default: None)
457 458
            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
459

460 461
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
462
        float, complex, either scalar (for full dot products)
463
                              or Field (for partial dot products)
464
        """
465
        if not isinstance(x, Field):
466 467
            raise TypeError("The dot-partner must be an instance of " +
                            "the NIFTy field class")
Theo Steininger's avatar
Theo Steininger committed
468

Martin Reinecke's avatar
Martin Reinecke committed
469
        if x._domain != self._domain:
470
            raise ValueError("Domain mismatch")
Theo Steininger's avatar
Theo Steininger committed
471

Martin Reinecke's avatar
Martin Reinecke committed
472
        ndom = len(self._domain)
473 474 475
        spaces = utilities.parse_spaces(spaces, ndom)

        if len(spaces) == ndom:
Martin Reinecke's avatar
Martin Reinecke committed
476
            return dobj.vdot(self.val, x.val)
477 478
        # 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
479
        return (self.conjugate()*x).sum(spaces=spaces)
Theo Steininger's avatar
Theo Steininger committed
480

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

Theo Steininger's avatar
Theo Steininger committed
484 485
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
486
        float
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
487
            The L2-norm of the field values.
csongor's avatar
csongor committed
488
        """
489
        return np.sqrt(abs(self.vdot(x=self)))
csongor's avatar
csongor committed
490

Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
491
    def conjugate(self):
Martin Reinecke's avatar
cleanup  
Martin Reinecke committed
492
        """ Returns the complex conjugate of the field.
Theo Steininger's avatar
Theo Steininger committed
493

494 495
        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
496 497
        Field
            The complex conjugated field.
csongor's avatar
csongor committed
498
        """
Martin Reinecke's avatar
Martin Reinecke committed
499
        return Field(self._domain, self.val.conjugate())
csongor's avatar
csongor committed
500

Theo Steininger's avatar
Theo Steininger committed
501
    # ---General unary/contraction methods---
502

Theo Steininger's avatar
Theo Steininger committed
503 504
    def __pos__(self):
        return self.copy()
505

Theo Steininger's avatar
Theo Steininger committed
506
    def __neg__(self):
507
        return Field(self._domain, -self.val)
csongor's avatar
csongor committed
508

Theo Steininger's avatar
Theo Steininger committed
509
    def __abs__(self):
510
        return Field(self._domain, abs(self.val))
csongor's avatar
csongor committed
511

512
    def _contraction_helper(self, op, spaces):
Theo Steininger's avatar
Theo Steininger committed
513
        if spaces is None:
514
            return getattr(self.val, op)()
515

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

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

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

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

Theo Steininger's avatar
Theo Steininger committed
526 527 528
        # 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
529
        else:
Martin Reinecke's avatar
Martin Reinecke committed
530
            return_domain = tuple(dom
Martin Reinecke's avatar
Martin Reinecke committed
531
                                  for i, dom in enumerate(self._domain)
Theo Steininger's avatar
Theo Steininger committed
532
                                  if i not in spaces)
533

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

536
    def sum(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
537 538 539 540 541 542 543 544 545 546 547 548 549 550
        """Sums up over the sub-domains given by `spaces`.

        Parameters
        ----------
        spaces : None, int or tuple of int (default: None)
            The summation is only carried out over the sub-domains in this
            tuple. If None, it is carried out over all sub-domains.

        Returns
        -------
        Field or scalar
            The result of the summation. If it is carried out over the entire
            domain, this is a scalar, otherwise a Field.
        """
551
        return self._contraction_helper('sum', spaces)
csongor's avatar
csongor committed
552

553
    def integrate(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570
        """Integrates over the sub-domains given by `spaces`.

        Integration is performed by summing over `self` multiplied by its
        volume factors.

        Parameters
        ----------
        spaces : None, int or tuple of int (default: None)
            The summation is only carried out over the sub-domains in this
            tuple. If None, it is carried out over all sub-domains.

        Returns
        -------
        Field or scalar
            The result of the integration. If it is carried out over the
            entire domain, this is a scalar, otherwise a Field.
        """
Martin Reinecke's avatar
Martin Reinecke committed
571 572 573 574 575
        swgt = self.scalar_weight(spaces)
        if swgt is not None:
            res = self.sum(spaces)
            res *= swgt
            return res
576 577 578
        tmp = self.weight(1, spaces=spaces)
        return tmp.sum(spaces)

579
    def prod(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
580 581 582 583 584 585 586 587 588 589 590 591 592 593
        """Computes the product over the sub-domains given by `spaces`.

        Parameters
        ----------
        spaces : None, int or tuple of int (default: None)
            The operation is only carried out over the sub-domains in this
            tuple. If None, it is carried out over all sub-domains.

        Returns
        -------
        Field or scalar
            The result of the product. If it is carried out over the entire
            domain, this is a scalar, otherwise a Field.
        """
594
        return self._contraction_helper('prod', spaces)
csongor's avatar
csongor committed
595

596 597
    def all(self, spaces=None):
        return self._contraction_helper('all', spaces)
csongor's avatar
csongor committed
598

599 600
    def any(self, spaces=None):
        return self._contraction_helper('any', spaces)
csongor's avatar
csongor committed
601

602
    def min(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
603 604 605 606 607 608 609 610 611 612 613 614 615 616
        """Determines the minimum over the sub-domains given by `spaces`.

        Parameters
        ----------
        spaces : None, int or tuple of int (default: None)
            The operation is only carried out over the sub-domains in this
            tuple. If None, it is carried out over all sub-domains.

        Returns
        -------
        Field or scalar
            The result of the operation. If it is carried out over the entire
            domain, this is a scalar, otherwise a Field.
        """
617
        return self._contraction_helper('min', spaces)
csongor's avatar
csongor committed
618

619
    def max(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
620 621 622 623 624 625 626 627 628 629 630 631 632 633
        """Determines the maximum over the sub-domains given by `spaces`.

        Parameters
        ----------
        spaces : None, int or tuple of int (default: None)
            The operation is only carried out over the sub-domains in this
            tuple. If None, it is carried out over all sub-domains.

        Returns
        -------
        Field or scalar
            The result of the operation. If it is carried out over the entire
            domain, this is a scalar, otherwise a Field.
        """
634
        return self._contraction_helper('max', spaces)
csongor's avatar
csongor committed
635

636
    def mean(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653
        """Determines the mean over the sub-domains given by `spaces`.

        ``x.mean(spaces)`` is equivalent to
        ``x.integrate(spaces)/x.total_volume(spaces)``.

        Parameters
        ----------
        spaces : None, int or tuple of int (default: None)
            The operation is only carried out over the sub-domains in this
            tuple. If None, it is carried out over all sub-domains.

        Returns
        -------
        Field or scalar
            The result of the operation. If it is carried out over the entire
            domain, this is a scalar, otherwise a Field.
        """
654 655
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('mean', spaces)
Martin Reinecke's avatar
Martin Reinecke committed
656
        # MR FIXME: not very efficient
657 658
        # MR FIXME: do we need "spaces" here?
        tmp = self.weight(1, spaces)
Martin Reinecke's avatar
Martin Reinecke committed
659
        return tmp.sum(spaces)*(1./tmp.total_volume(spaces))
csongor's avatar
csongor committed
660

661
    def var(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
662 663 664 665 666 667 668 669 670 671 672 673 674 675
        """Determines the variance over the sub-domains given by `spaces`.

        Parameters
        ----------
        spaces : None, int or tuple of int (default: None)
            The operation is only carried out over the sub-domains in this
            tuple. If None, it is carried out over all sub-domains.

        Returns
        -------
        Field or scalar
            The result of the operation. If it is carried out over the entire
            domain, this is a scalar, otherwise a Field.
        """
676 677
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('var', spaces)
Martin Reinecke's avatar
Martin Reinecke committed
678 679 680
        # MR FIXME: not very efficient or accurate
        m1 = self.mean(spaces)
        if np.issubdtype(self.dtype, np.complexfloating):
681
            sq = abs(self-m1)**2
Martin Reinecke's avatar
Martin Reinecke committed
682
        else:
683 684
            sq = (self-m1)**2
        return sq.mean(spaces)
csongor's avatar
csongor committed
685

686
    def std(self, spaces=None):
Martin Reinecke's avatar
Martin Reinecke committed
687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703
        """Determines the standard deviation over the sub-domains given by
        `spaces`.

        ``x.std(spaces)`` is equivalent to ``sqrt(x.var(spaces))``.

        Parameters
        ----------
        spaces : None, int or tuple of int (default: None)
            The operation is only carried out over the sub-domains in this
            tuple. If None, it is carried out over all sub-domains.

        Returns
        -------
        Field or scalar
            The result of the operation. If it is carried out over the entire
            domain, this is a scalar, otherwise a Field.
        """
704
        from .sugar import sqrt
705 706 707
        if self.scalar_weight(spaces) is not None:
            return self._contraction_helper('std', spaces)
        return sqrt(self.var(spaces))
csongor's avatar
csongor committed
708

709 710 711
    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
712
        if other._domain != self._domain:
713
            raise ValueError("domains are incompatible.")
Martin Reinecke's avatar
Martin Reinecke committed
714
        self.local_data[()] = other.local_data[()]
715

Theo Steininger's avatar
Theo Steininger committed
716
    def __repr__(self):
Martin Reinecke's avatar
Martin Reinecke committed
717
        return "<nifty4.Field>"
Theo Steininger's avatar
Theo Steininger committed
718 719

    def __str__(self):
Martin Reinecke's avatar
Martin Reinecke committed
720
        return "nifty4.Field instance\n- domain      = " + \
721
               self._domain.__str__() + \
722
               "\n- val         = " + repr(self.val)
Martin Reinecke's avatar
cleanup  
Martin Reinecke committed
723

724 725 726 727 728 729 730 731 732
    def equivalent(self, other):
        if self is other:
            return True
        if not isinstance(other, Field):
            return False
        if self._domain != other._domain:
            return False
        return (self._val == other._val).all()

733 734 735 736 737 738 739 740 741 742
for op in ["__add__", "__radd__", "__iadd__",
           "__sub__", "__rsub__", "__isub__",
           "__mul__", "__rmul__", "__imul__",
           "__div__", "__rdiv__", "__idiv__",
           "__truediv__", "__rtruediv__", "__itruediv__",
           "__floordiv__", "__rfloordiv__", "__ifloordiv__",
           "__pow__", "__rpow__", "__ipow__",
           "__lt__", "__le__", "__gt__", "__ge__", "__eq__", "__ne__"]:
    def func(op):
        def func2(self, other):
Martin Reinecke's avatar
cleanup  
Martin Reinecke committed
743 744 745 746 747 748 749 750 751 752 753 754
            # if other is a field, make sure that the domains match
            if isinstance(other, Field):
                if other._domain != self._domain:
                    raise ValueError("domains are incompatible.")
                tval = getattr(self.val, op)(other.val)
                return self if tval is self.val else Field(self._domain, tval)

            if np.isscalar(other) or isinstance(other, dobj.data_object):
                tval = getattr(self.val, op)(other)
                return self if tval is self.val else Field(self._domain, tval)

            return NotImplemented
755 756
        return func2
    setattr(Field, op, func(op))