field.py 35 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
14
15
16
17
#
# Copyright(C) 2013-2017 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
18

csongor's avatar
csongor committed
19
from __future__ import division
Martin Reinecke's avatar
Martin Reinecke committed
20
21
from builtins import zip
from builtins import range
22

23
import ast
csongor's avatar
csongor committed
24
25
import numpy as np

Martin Reinecke's avatar
Martin Reinecke committed
26
from .domain_object import DomainObject
27

Martin Reinecke's avatar
Martin Reinecke committed
28
from .spaces.power_space import PowerSpace
csongor's avatar
csongor committed
29

Martin Reinecke's avatar
Martin Reinecke committed
30
31
from . import nifty_utilities as utilities
from .random import Random
Martin Reinecke's avatar
Martin Reinecke committed
32
from functools import reduce
33

csongor's avatar
csongor committed
34

Martin Reinecke's avatar
Martin Reinecke committed
35
class Field(object):
Theo Steininger's avatar
Theo Steininger committed
36
37
38
    """ The discrete representation of a continuous field over multiple spaces.

    In NIFTY, Fields are used to store data arrays and carry all the needed
39
    metainformation (i.e. the domain) for operators to be able to work on them.
Theo Steininger's avatar
Theo Steininger committed
40
41
    In addition Field has methods to work with power-spectra.

42
43
44
45
    Parameters
    ----------
    domain : DomainObject
        One of the space types NIFTY supports. RGSpace, GLSpace, HPSpace,
Theo Steininger's avatar
Theo Steininger committed
46
        LMSpace or PowerSpace. It might also be a FieldArray, which is
47
        an unstructured domain.
Theo Steininger's avatar
Theo Steininger committed
48

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
49
    val : scalar, numpy.ndarray, Field
50
51
52
        The values the array should contain after init. A scalar input will
        fill the whole array with this scalar. If an array is provided the
        array's dimensions must match the domain's.
Theo Steininger's avatar
Theo Steininger committed
53

54
55
    dtype : type
        A numpy.type. Most common are int, float and complex.
Theo Steininger's avatar
Theo Steininger committed
56

57
58
59
60
    copy: boolean

    Attributes
    ----------
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
61
    val : numpy.ndarray
Theo Steininger's avatar
Theo Steininger committed
62

63
64
65
66
67
68
    domain : DomainObject
        See Parameters.
    domain_axes : tuple of tuples
        Enumerates the axes of the Field
    dtype : type
        Contains the datatype stored in the Field.
Theo Steininger's avatar
Theo Steininger committed
69

70
71
72
73
74
75
76
    Raise
    -----
    TypeError
        Raised if
            *the given domain contains something that is not a DomainObject
             instance
            *val is an array that has a different dimension than the domain
Theo Steininger's avatar
Theo Steininger committed
77

78
    """
79

Theo Steininger's avatar
Theo Steininger committed
80
    # ---Initialization methods---
81

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
82
    def __init__(self, domain=None, val=None, dtype=None, copy=False):
83
        self.domain = self._parse_domain(domain=domain, val=val)
84
        self.domain_axes = self._get_axes_tuple(self.domain)
85

Martin Reinecke's avatar
Martin Reinecke committed
86
87
88
        shape_tuple = tuple(sp.shape for sp in self.domain)
        if len(shape_tuple)==0:
            global_shape = ()
89
        else:
Martin Reinecke's avatar
Martin Reinecke committed
90
91
92
93
            global_shape = reduce(lambda x, y: x + y, shape_tuple)
        dtype = self._infer_dtype(dtype=dtype, val=val)
        self._val = np.empty(global_shape,dtype=dtype)
        self.set_val(new_val=val, copy=copy)
csongor's avatar
csongor committed
94

95
    def _parse_domain(self, domain, val=None):
96
        if domain is None:
97
98
            if isinstance(val, Field):
                domain = val.domain
Martin Reinecke's avatar
Martin Reinecke committed
99
            elif np.isscalar(val):
100
                domain = ()
Martin Reinecke's avatar
Martin Reinecke committed
101
102
            else:
                raise TypeError("could not infer domain from value")
103
        elif isinstance(domain, DomainObject):
104
            domain = (domain,)
105
106
107
        elif not isinstance(domain, tuple):
            domain = tuple(domain)

csongor's avatar
csongor committed
108
        for d in domain:
109
            if not isinstance(d, DomainObject):
110
111
                raise TypeError(
                    "Given domain contains something that is not a "
112
                    "DomainObject instance.")
csongor's avatar
csongor committed
113
114
        return domain

Theo Steininger's avatar
Theo Steininger committed
115
116
117
118
119
120
121
122
123
124
    def _get_axes_tuple(self, things_with_shape, start=0):
        i = start
        axes_list = []
        for thing in things_with_shape:
            l = []
            for j in range(len(thing.shape)):
                l += [i]
                i += 1
            axes_list += [tuple(l)]
        return tuple(axes_list)
125

126
    def _infer_dtype(self, dtype, val):
csongor's avatar
csongor committed
127
        if dtype is None:
128
            try:
129
                dtype = val.dtype
130
            except AttributeError:
Theo Steininger's avatar
Theo Steininger committed
131
132
133
                try:
                    if val is None:
                        raise TypeError
134
                    dtype = np.result_type(val)
Theo Steininger's avatar
Theo Steininger committed
135
                except(TypeError):
Martin Reinecke's avatar
more    
Martin Reinecke committed
136
                    dtype = np.dtype(np.float64)
Theo Steininger's avatar
Theo Steininger committed
137
        else:
138
            dtype = np.dtype(dtype)
139

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
140
        return np.result_type(dtype, np.float)
141

142
    # ---Factory methods---
143

144
    @classmethod
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
145
    def from_random(cls, random_type, domain=None, dtype=None, **kwargs):
146
147
148
149
150
        """ Draws a random field with the given parameters.

        Parameters
        ----------
        cls : class
Theo Steininger's avatar
Theo Steininger committed
151

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

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

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

162
163
164
165
166
167
168
        Returns
        -------
        out : Field
            The output object.

        See Also
        --------
169
        power_synthesize
Theo Steininger's avatar
Theo Steininger committed
170

171
172

        """
Theo Steininger's avatar
Theo Steininger committed
173

174
        # create a initially empty field
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
175
        f = cls(domain=domain, dtype=dtype)
176

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
177
        # extract the data from f and apply the appropriate
178
179
180
        # random number generator to it
        sample = f.get_val(copy=False)
        generator_function = getattr(Random, random_type)
181

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
182
183
        sample[:]=generator_function(dtype=f.dtype,
                                             shape=sample.shape,
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
184
                                             **kwargs)
185
186
        return f

187
188
    # ---Powerspectral methods---

Martin Reinecke's avatar
Martin Reinecke committed
189
190
    def power_analyze(self, spaces=None, binbounds=None,
                      keep_phase_information=False):
Theo Steininger's avatar
Theo Steininger committed
191
        """ Computes the square root power spectrum for a subspace of `self`.
Theo Steininger's avatar
Theo Steininger committed
192

Theo Steininger's avatar
Theo Steininger committed
193
194
195
        Creates a PowerSpace for the space addressed by `spaces` with the given
        binning and computes the power spectrum as a Field over this
        PowerSpace. This can only be done if the subspace to  be analyzed is a
196
        harmonic space. The resulting field has the same units as the initial
Theo Steininger's avatar
Theo Steininger committed
197
        field, corresponding to the square root of the power spectrum.
198
199
200

        Parameters
        ----------
Theo Steininger's avatar
Theo Steininger committed
201
202
203
204
205
        spaces : int *optional*
            The subspace for which the powerspectrum shall be computed
            (default : None).
        binbounds : array-like *optional*
            Inner bounds of the bins (default : None).
Martin Reinecke's avatar
Martin Reinecke committed
206
            if binbounds==None : bins are inferred.
207
208
209
210
211
212
213
214
215
216
        keep_phase_information : boolean, *optional*
            If False, return a real-valued result containing the power spectrum
            of the input Field.
            If True, return a complex-valued result whose real component
            contains the power spectrum computed from the real part of the
            input Field, and whose imaginary component contains the power
            spectrum computed from the imaginary part of the input Field.
            The absolute value of this result should be identical to the output
            of power_analyze with keep_phase_information=False.
            (default : False).
Theo Steininger's avatar
Theo Steininger committed
217

218
219
220
221
        Raise
        -----
        ValueError
            Raised if
Theo Steininger's avatar
Theo Steininger committed
222
223
                *len(domain) is != 1 when spaces==None
                *len(spaces) is != 1 if not None
224
                *the analyzed space is not harmonic
Theo Steininger's avatar
Theo Steininger committed
225

226
227
        Returns
        -------
Theo Steininger's avatar
Theo Steininger committed
228
        out : Field
Martin Reinecke's avatar
Martin Reinecke committed
229
            The output object. Its domain is a PowerSpace and it contains
230
231
232
233
234
            the power spectrum of 'self's field.

        See Also
        --------
        power_synthesize, PowerSpace
Theo Steininger's avatar
Theo Steininger committed
235

236
        """
Theo Steininger's avatar
Theo Steininger committed
237

Theo Steininger's avatar
Theo Steininger committed
238
        # check if all spaces in `self.domain` are either harmonic or
239
240
241
        # power_space instances
        for sp in self.domain:
            if not sp.harmonic and not isinstance(sp, PowerSpace):
Martin Reinecke's avatar
Martin Reinecke committed
242
                raise TypeError(
243
                    "Field has a space in `domain` which is neither "
244
245
246
                    "harmonic nor a PowerSpace.")

        # check if the `spaces` input is valid
247
248
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
249
            spaces = list(range(len(self.domain)))
250
251

        if len(spaces) == 0:
252
253
            raise ValueError(
                "No space for analysis specified.")
254

255
256
257
258
259
260
261
262
263
264
265
266
267
        if keep_phase_information:
            parts_val = self._hermitian_decomposition(
                                              domain=self.domain,
                                              val=self.val,
                                              spaces=spaces,
                                              domain_axes=self.domain_axes,
                                              preserve_gaussian_variance=False)
            parts = [self.copy_empty().set_val(part_val, copy=False)
                     for part_val in parts_val]
        else:
            parts = [self]

        parts = [abs(part)**2 for part in parts]
268
269

        for space_index in spaces:
270
271
            parts = [self._single_power_analyze(
                                work_field=part,
272
                                space_index=space_index,
273
274
                                binbounds=binbounds)
                     for part in parts]
275

276
277
278
279
280
281
        if keep_phase_information:
            result_field = parts[0] + 1j*parts[1]
        else:
            result_field = parts[0]

        return result_field
282
283

    @classmethod
Martin Reinecke's avatar
Martin Reinecke committed
284
    def _single_power_analyze(cls, work_field, space_index, binbounds):
285

286
        if not work_field.domain[space_index].harmonic:
287
288
            raise ValueError(
                "The analyzed space must be harmonic.")
289

290
291
292
293
294
295
        # Create the target PowerSpace instance:
        # If the associated signal-space field was real, we extract the
        # hermitian and anti-hermitian parts of `self` and put them
        # into the real and imaginary parts of the power spectrum.
        # If it was complex, all the power is put into a real power spectrum.

296
        harmonic_domain = work_field.domain[space_index]
297
        power_domain = PowerSpace(harmonic_partner=harmonic_domain,
Theo Steininger's avatar
Theo Steininger committed
298
                                  binbounds=binbounds)
299
300
        power_spectrum = cls._calculate_power_spectrum(
                                field_val=work_field.val,
Martin Reinecke's avatar
Martin Reinecke committed
301
                                pdomain=power_domain,
302
                                axes=work_field.domain_axes[space_index])
303
304

        # create the result field and put power_spectrum into it
305
        result_domain = list(work_field.domain)
306
        result_domain[space_index] = power_domain
307
        result_dtype = power_spectrum.dtype
308

309
        result_field = work_field.copy_empty(
310
                   domain=result_domain,
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
311
                   dtype=result_dtype)
312
313
314
315
        result_field.set_val(new_val=power_spectrum, copy=False)

        return result_field

316
    @classmethod
Martin Reinecke's avatar
Martin Reinecke committed
317
    def _calculate_power_spectrum(cls, field_val, pdomain, axes=None):
318

Martin Reinecke's avatar
Martin Reinecke committed
319
320
321
        pindex = pdomain.pindex
        # MR FIXME: how about iterating over slices, instead of replicating
        # pindex? Would save memory and probably isn't slower.
322
        if axes is not None:
323
324
325
326
            pindex = cls._shape_up_pindex(
                            pindex=pindex,
                            target_shape=field_val.shape,
                            axes=axes)
Theo Steininger's avatar
Theo Steininger committed
327

Martin Reinecke's avatar
Martin Reinecke committed
328
        power_spectrum = utilities.bincount_axis(pindex, weights=field_val,
329
                                         axis=axes)
Martin Reinecke's avatar
Martin Reinecke committed
330
        rho = pdomain.rho
331
332
333
334
335
336
337
338
        if axes is not None:
            new_rho_shape = [1, ] * len(power_spectrum.shape)
            new_rho_shape[axes[0]] = len(rho)
            rho = rho.reshape(new_rho_shape)
        power_spectrum /= rho

        return power_spectrum

339
    @staticmethod
Martin Reinecke's avatar
Martin Reinecke committed
340
    def _shape_up_pindex(pindex, target_shape, axes):
Theo Steininger's avatar
Theo Steininger committed
341
        semiscaled_local_shape = [1, ] * len(target_shape)
Theo Steininger's avatar
Theo Steininger committed
342
        for i in range(len(axes)):
Martin Reinecke's avatar
Martin Reinecke committed
343
344
            semiscaled_local_shape[axes[i]] = pindex.shape[i]
        local_data = pindex
Theo Steininger's avatar
Theo Steininger committed
345
        semiscaled_local_data = local_data.reshape(semiscaled_local_shape)
Martin Reinecke's avatar
Martin Reinecke committed
346
347
        result_obj = np.empty(target_shape, dtype=pindex.dtype)
        result_obj[:] = semiscaled_local_data
348
349
350

        return result_obj

351
    def power_synthesize(self, spaces=None, real_power=True, real_signal=True,
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
352
                         mean=None, std=None):
Theo Steininger's avatar
Theo Steininger committed
353
        """ Yields a sampled field with `self`**2 as its power spectrum.
Theo Steininger's avatar
Theo Steininger committed
354

Theo Steininger's avatar
Theo Steininger committed
355
        This method draws a Gaussian random field in the harmonic partner
Martin Reinecke's avatar
typos    
Martin Reinecke committed
356
        domain of this field's domains, using this field as power spectrum.
Theo Steininger's avatar
Theo Steininger committed
357

358
359
360
        Parameters
        ----------
        spaces : {tuple, int, None} *optional*
Theo Steininger's avatar
Theo Steininger committed
361
362
363
            Specifies the subspace containing all the PowerSpaces which
            should be converted (default : None).
            if spaces==None : Tries to convert the whole domain.
364
        real_power : boolean *optional*
Theo Steininger's avatar
Theo Steininger committed
365
366
            Determines whether the power spectrum is treated as intrinsically
            real or complex (default : True).
367
        real_signal : boolean *optional*
Theo Steininger's avatar
Theo Steininger committed
368
369
370
371
372
373
            True will result in a purely real signal-space field
            (default : True).
        mean : float *optional*
            The mean of the Gaussian noise field which is used for the Field
            synthetization (default : None).
            if mean==None : mean will be set to 0
374
        std : float *optional*
Theo Steininger's avatar
Theo Steininger committed
375
376
377
            The standard deviation of the Gaussian noise field which is used
            for the Field synthetization (default : None).
            if std==None : std will be set to 1
Theo Steininger's avatar
Theo Steininger committed
378

379
380
381
382
        Returns
        -------
        out : Field
            The output object. A random field created with the power spectrum
Theo Steininger's avatar
Theo Steininger committed
383
            stored in the `spaces` in `self`.
384

Theo Steininger's avatar
Theo Steininger committed
385
386
387
388
389
390
        Notes
        -----
        For this the spaces specified by `spaces` must be a PowerSpace.
        This expects this field to be the square root of a power spectrum, i.e.
        to have the unit of the field to be sampled.

391
392
393
        See Also
        --------
        power_analyze
Theo Steininger's avatar
Theo Steininger committed
394
395
396
397
398

        Raises
        ------
        ValueError : If domain specified by `spaces` is not a PowerSpace.

399
        """
Theo Steininger's avatar
Theo Steininger committed
400

401
402
403
        # check if the `spaces` input is valid
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))

Theo Steininger's avatar
Theo Steininger committed
404
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
405
            spaces = list(range(len(self.domain)))
Theo Steininger's avatar
Theo Steininger committed
406

407
408
409
410
411
        for power_space_index in spaces:
            power_space = self.domain[power_space_index]
            if not isinstance(power_space, PowerSpace):
                raise ValueError("A PowerSpace is needed for field "
                                 "synthetization.")
412
413
414

        # create the result domain
        result_domain = list(self.domain)
415
416
        for power_space_index in spaces:
            power_space = self.domain[power_space_index]
417
            harmonic_domain = power_space.harmonic_partner
418
            result_domain[power_space_index] = harmonic_domain
419
420
421

        # create random samples: one or two, depending on whether the
        # power spectrum is real or complex
422
        if real_power:
423
            result_list = [None]
424
425
        else:
            result_list = [None, None]
426

427
428
        result_list = [self.__class__.from_random(
                             'normal',
429
430
431
                             mean=mean,
                             std=std,
                             domain=result_domain,
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
432
                             dtype=np.complex)
433
434
435
436
437
438
                       for x in result_list]

        # from now on extract the values from the random fields for further
        # processing without killing the fields.
        # if the signal-space field should be real, hermitianize the field
        # components
439

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
440
        spec = self.val.copy()
441
442
        spec = np.sqrt(spec)

443
444
445
446
447
448
449
        for power_space_index in spaces:
            spec = self._spec_to_rescaler(spec, result_list, power_space_index)
        local_rescaler = spec

        result_val_list = [x.val for x in result_list]

        # apply the rescaler to the random fields
Martin Reinecke's avatar
Martin Reinecke committed
450
        result_val_list[0] *= local_rescaler.real
451
452

        if not real_power:
Martin Reinecke's avatar
Martin Reinecke committed
453
            result_val_list[1] *= local_rescaler.imag
454

455
        if real_signal:
456
            result_val_list = [self._hermitian_decomposition(
457
458
459
460
461
                                            result_domain,
                                            result_val,
                                            spaces,
                                            result_list[0].domain_axes,
                                            preserve_gaussian_variance=True)[0]
462
                               for result_val in result_val_list]
463
464
465
466
467
468
469

        # store the result into the fields
        [x.set_val(new_val=y, copy=False) for x, y in
            zip(result_list, result_val_list)]

        if real_power:
            result = result_list[0]
470
471
472
            if not issubclass(result_val_list[0].dtype.type,
                              np.complexfloating):
                result = result.real
473
        else:
474
475
476
477
            result = result_list[0] + 1j*result_list[1]

        return result

478
    @staticmethod
479
480
    def _hermitian_decomposition(domain, val, spaces, domain_axes,
                                 preserve_gaussian_variance=False):
481

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
482
483
        h = val.real.copy()
        a = 1j * val.imag.copy()
484
485

        # correct variance
486
        if preserve_gaussian_variance:
Martin Reinecke's avatar
Martin Reinecke committed
487
488
            assert issubclass(val.dtype.type, np.complexfloating),\
                    "complex input field is needed here"
489
490
491
            h *= np.sqrt(2)
            a *= np.sqrt(2)

492
493
        return (h, a)

494
495
    def _spec_to_rescaler(self, spec, result_list, power_space_index):
        power_space = self.domain[power_space_index]
496
497
498

        # weight the random fields with the power spectrum
        # therefore get the pindex from the power space
499
        pindex = power_space.pindex
500
501
502
503

        # Now use numpy advanced indexing in order to put the entries of the
        # power spectrum into the appropriate places of the pindex array.
        # Do this for every 'pindex-slice' in parallel using the 'slice(None)'s
Martin Reinecke's avatar
Martin Reinecke committed
504
        local_pindex = pindex
505

506
507
508
509
510
        local_blow_up = [slice(None)]*len(spec.shape)
        # it is important to count from behind, since spec potentially grows
        # with every iteration
        index = self.domain_axes[power_space_index][0]-len(self.shape)
        local_blow_up[index] = local_pindex
511
        # here, the power_spectrum is distributed into the new shape
512
513
        local_rescaler = spec[local_blow_up]
        return local_rescaler
514

Theo Steininger's avatar
Theo Steininger committed
515
    # ---Properties---
516

Theo Steininger's avatar
Theo Steininger committed
517
    def set_val(self, new_val=None, copy=False):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
518
        """ Sets the field's data object.
519
520
521

        Parameters
        ----------
Theo Steininger's avatar
Theo Steininger committed
522
        new_val : scalar, array-like, Field, None *optional*
523
524
            The values to be stored in the field.
            {default : None}
Theo Steininger's avatar
Theo Steininger committed
525

526
        copy : boolean, *optional*
Theo Steininger's avatar
Theo Steininger committed
527
528
            If False, Field tries to not copy the input data but use it
            directly.
529
530
531
532
533
534
            {default : False}
        See Also
        --------
        val

        """
Theo Steininger's avatar
Theo Steininger committed
535

Martin Reinecke's avatar
Martin Reinecke committed
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
        if new_val is None:
            pass
        elif isinstance(new_val, Field):
            if self.domain!=new_val.domain:
                raise ValueError("Domain mismatch")
            if copy:
                self._val[()] = new_val.val
            else:
                self._val = new_val.val
        elif (np.isscalar(new_val)):
            self._val[()]=new_val
        elif isinstance(new_val, np.ndarray):
            if copy:
                self._val[()] = new_val
            else:
                if self.shape!=new_val.shape:
                    raise ValueError("Shape mismatch")
                self._val = new_val
        else:
            raise TypeError("unknown source type")
556
        return self
csongor's avatar
csongor committed
557

558
    def get_val(self, copy=False):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
559
        """ Returns the data object associated with this Field.
560
561
562
563

        Parameters
        ----------
        copy : boolean
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
564
            If true, a copy of the Field's underlying data object
Theo Steininger's avatar
Theo Steininger committed
565
            is returned.
Theo Steininger's avatar
Theo Steininger committed
566

567
568
        Returns
        -------
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
569
        out : numpy.ndarray
570
571
572
573
574
575

        See Also
        --------
        val

        """
Theo Steininger's avatar
Theo Steininger committed
576

Martin Reinecke's avatar
Martin Reinecke committed
577
        return self._val.copy() if copy else self._val
csongor's avatar
csongor committed
578

Theo Steininger's avatar
Theo Steininger committed
579
580
    @property
    def val(self):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
581
        """ Returns the data object associated with this Field.
Theo Steininger's avatar
Theo Steininger committed
582

583
584
        Returns
        -------
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
585
        out : numpy.ndarray
586
587
588
589
590
591

        See Also
        --------
        get_val

        """
Theo Steininger's avatar
Theo Steininger committed
592

Martin Reinecke's avatar
Martin Reinecke committed
593
        return self._val
csongor's avatar
csongor committed
594

Theo Steininger's avatar
Theo Steininger committed
595
596
    @val.setter
    def val(self, new_val):
597
        self.set_val(new_val=new_val, copy=False)
csongor's avatar
csongor committed
598

Martin Reinecke's avatar
Martin Reinecke committed
599
600
601
602
    @property
    def dtype(self):
        return self._val.dtype

603
604
    @property
    def shape(self):
Theo Steininger's avatar
Theo Steininger committed
605
        """ Returns the total shape of the Field's data array.
Theo Steininger's avatar
Theo Steininger committed
606

607
608
609
        Returns
        -------
        out : tuple
Martin Reinecke's avatar
Martin Reinecke committed
610
            The output object. The tuple contains the dimensions of the spaces
611
612
613
614
615
616
617
            in domain.

        See Also
        --------
        dim

        """
Martin Reinecke's avatar
Martin Reinecke committed
618
        return self._val.shape
csongor's avatar
csongor committed
619

620
621
    @property
    def dim(self):
Theo Steininger's avatar
Theo Steininger committed
622
        """ Returns the total number of pixel-dimensions the field has.
Theo Steininger's avatar
Theo Steininger committed
623

Theo Steininger's avatar
Theo Steininger committed
624
        Effectively, all values from shape are multiplied.
Theo Steininger's avatar
Theo Steininger committed
625

626
627
628
629
630
631
632
633
634
635
        Returns
        -------
        out : int
            The dimension of the Field.

        See Also
        --------
        shape

        """
Theo Steininger's avatar
Theo Steininger committed
636

637
        dim_tuple = tuple(sp.dim for sp in self.domain)
Theo Steininger's avatar
Theo Steininger committed
638
        try:
Martin Reinecke's avatar
Martin Reinecke committed
639
            return int(reduce(lambda x, y: x * y, dim_tuple))
Theo Steininger's avatar
Theo Steininger committed
640
641
        except TypeError:
            return 0
csongor's avatar
csongor committed
642

643
644
    @property
    def dof(self):
Theo Steininger's avatar
Theo Steininger committed
645
646
647
648
649
650
        """ Returns the total number of degrees of freedom the Field has. For
        real Fields this is equal to `self.dim`. For complex Fields it is
        2*`self.dim`.

        """

Theo Steininger's avatar
Theo Steininger committed
651
652
653
654
655
656
657
        dof = self.dim
        if issubclass(self.dtype.type, np.complexfloating):
            dof *= 2
        return dof

    @property
    def total_volume(self):
Theo Steininger's avatar
Theo Steininger committed
658
659
660
        """ Returns the total volume of all spaces in the domain.
        """

Theo Steininger's avatar
Theo Steininger committed
661
        volume_tuple = tuple(sp.total_volume for sp in self.domain)
662
        try:
Theo Steininger's avatar
Theo Steininger committed
663
            return reduce(lambda x, y: x * y, volume_tuple)
664
        except TypeError:
Theo Steininger's avatar
Theo Steininger committed
665
            return 0.
666

Theo Steininger's avatar
Theo Steininger committed
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
    @property
    def real(self):
        """ The real part of the field (data is not copied).
        """
        real_part = self.val.real
        result = self.copy_empty(dtype=real_part.dtype)
        result.set_val(new_val=real_part, copy=False)
        return result

    @property
    def imag(self):
        """ The imaginary part of the field (data is not copied).
        """
        real_part = self.val.imag
        result = self.copy_empty(dtype=real_part.dtype)
        result.set_val(new_val=real_part, copy=False)
        return result

Theo Steininger's avatar
Theo Steininger committed
685
    # ---Special unary/binary operations---
686

csongor's avatar
csongor committed
687

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
688
    def copy(self, domain=None, dtype=None):
689
        """ Returns a full copy of the Field.
Theo Steininger's avatar
Theo Steininger committed
690

691
692
        If no keyword arguments are given, the returned object will be an
        identical copy of the original Field. By explicit specification one is
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
693
        able to define the domain and the dtype of the returned Field.
694
695
696
697
698

        Parameters
        ----------
        domain : DomainObject
            The new domain the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
699

700
701
        dtype : type
            The new dtype the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
702

703
704
705
706
707
708
709
710
711
712
        Returns
        -------
        out : Field
            The output object. An identical copy of 'self'.

        See Also
        --------
        copy_empty

        """
Theo Steininger's avatar
Theo Steininger committed
713

Theo Steininger's avatar
Theo Steininger committed
714
        copied_val = self.get_val(copy=True)
715
716
        new_field = self.copy_empty(
                                domain=domain,
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
717
                                dtype=dtype)
Theo Steininger's avatar
Theo Steininger committed
718
719
        new_field.set_val(new_val=copied_val, copy=False)
        return new_field
csongor's avatar
csongor committed
720

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
721
    def copy_empty(self, domain=None, dtype=None):
722
723
724
        """ Returns an empty copy of the Field.

        If no keyword arguments are given, the returned object will be an
Theo Steininger's avatar
Theo Steininger committed
725
726
727
        identical copy of the original Field. The memory for the data array
        is only allocated but not actively set to any value
        (c.f. numpy.ndarray.copy_empty). By explicit specification one is able
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
728
        to change the domain and the dtype of the returned Field.
Theo Steininger's avatar
Theo Steininger committed
729

730
731
732
733
        Parameters
        ----------
        domain : DomainObject
            The new domain the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
734

735
736
        dtype : type
            The new dtype the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
737

738
739
740
        Returns
        -------
        out : Field
Theo Steininger's avatar
Theo Steininger committed
741
            The output object.
742
743
744
745
746
747

        See Also
        --------
        copy

        """
Theo Steininger's avatar
Theo Steininger committed
748

Theo Steininger's avatar
Theo Steininger committed
749
750
        if domain is None:
            domain = self.domain
csongor's avatar
csongor committed
751
        else:
Theo Steininger's avatar
Theo Steininger committed
752
            domain = self._parse_domain(domain)
csongor's avatar
csongor committed
753

Theo Steininger's avatar
Theo Steininger committed
754
755
756
757
        if dtype is None:
            dtype = self.dtype
        else:
            dtype = np.dtype(dtype)
csongor's avatar
csongor committed
758

Theo Steininger's avatar
Theo Steininger committed
759
760
        fast_copyable = True
        try:
Martin Reinecke's avatar
Martin Reinecke committed
761
            for i in range(len(self.domain)):
Theo Steininger's avatar
Theo Steininger committed
762
763
764
765
766
767
                if self.domain[i] is not domain[i]:
                    fast_copyable = False
                    break
        except IndexError:
            fast_copyable = False

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
768
        if (fast_copyable and dtype == self.dtype):
Theo Steininger's avatar
Theo Steininger committed
769
770
            new_field = self._fast_copy_empty()
        else:
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
771
            new_field = Field(domain=domain, dtype=dtype)
Theo Steininger's avatar
Theo Steininger committed
772
        return new_field
csongor's avatar
csongor committed
773

Theo Steininger's avatar
Theo Steininger committed
774
775
776
777
778
779
    def _fast_copy_empty(self):
        # make an empty field
        new_field = EmptyField()
        # repair its class
        new_field.__class__ = self.__class__
        # copy domain, codomain and val
Martin Reinecke's avatar
Martin Reinecke committed
780
        for key, value in list(self.__dict__.items()):
781
            if key != '_val':
Theo Steininger's avatar
Theo Steininger committed
782
783
                new_field.__dict__[key] = value
            else:
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
784
                new_field.__dict__[key] = np.empty_like(self.val)
Theo Steininger's avatar
Theo Steininger committed
785
786
787
        return new_field

    def weight(self, power=1, inplace=False, spaces=None):
Theo Steininger's avatar
Theo Steininger committed
788
        """ Weights the pixels of `self` with their invidual pixel-volume.
789
790
791
792

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

795
        inplace : boolean
Theo Steininger's avatar
Theo Steininger committed
796
797
            If True, `self` will be weighted and returned. Otherwise, a copy
            is made.
Theo Steininger's avatar
Theo Steininger committed
798

Theo Steininger's avatar
Theo Steininger committed
799
800
        spaces : tuple of ints
            Determines on which subspace the operation takes place.
Theo Steininger's avatar
Theo Steininger committed
801

802
803
804
        Returns
        -------
        out : Field
Theo Steininger's avatar
Theo Steininger committed
805
            The weighted field.
806
807

        """
808
        if inplace:
csongor's avatar
csongor committed
809
810
811
812
            new_field = self
        else:
            new_field = self.copy_empty()

813
        new_val = self.get_val(copy=False)
csongor's avatar
csongor committed
814

815
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
csongor's avatar
csongor committed
816
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
817
            spaces = list(range(len(self.domain)))
csongor's avatar
csongor committed
818

819
        for ind, sp in enumerate(self.domain):
Theo Steininger's avatar
Theo Steininger committed
820
821
822
823
824
            if ind in spaces:
                new_val = sp.weight(new_val,
                                    power=power,
                                    axes=self.domain_axes[ind],
                                    inplace=inplace)
825
826

        new_field.set_val(new_val=new_val, copy=False)
csongor's avatar
csongor committed
827
828
        return new_field

Martin Reinecke's avatar
Martin Reinecke committed
829
    def vdot(self, x=None, spaces=None, bare=False):
Theo Steininger's avatar
Theo Steininger committed
830
        """ Computes the volume-factor-aware dot product of 'self' with x.
Theo Steininger's avatar
Theo Steininger committed
831

832
833
834
        Parameters
        ----------
        x : Field
Theo Steininger's avatar
Theo Steininger committed
835
            The domain of x must contain `self.domain`
Theo Steininger's avatar
Theo Steininger committed
836

Theo Steininger's avatar
Theo Steininger committed
837
838
839
        spaces : tuple of ints
            If the domain of `self` and `x` are not the same, `spaces` specfies
            the mapping.
Theo Steininger's avatar
Theo Steininger committed
840

841
        bare : boolean
Theo Steininger's avatar
Theo Steininger committed
842
            If true, no volume factors will be included in the computation.
Theo Steininger's avatar
Theo Steininger committed
843

844
845
846
        Returns
        -------
        out : float, complex
Theo Steininger's avatar
Theo Steininger committed
847

848
        """
849
850
851
        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
852

Martin Reinecke's avatar
Martin Reinecke committed
853
        # Compute the dot respecting the fact of discrete/continuous spaces
Theo Steininger's avatar
Theo Steininger committed
854
855
856
857
858
        if bare:
            y = self
        else:
            y = self.weight(power=1)

859
860
861
        if spaces is None:
            x_val = x.get_val(copy=False)
            y_val = y.get_val(copy=False)
Martin Reinecke's avatar
Martin Reinecke committed
862
            result = (y_val.conjugate() * x_val).sum()
863
864
865
866
            return result
        else:
            # create a diagonal operator which is capable of taking care of the
            # axes-matching
867
            from .operators.diagonal_operator import DiagonalOperator
868
869
870
871
872
873
            diagonal = y.val.conjugate()
            diagonalOperator = DiagonalOperator(domain=y.domain,
                                                diagonal=diagonal,
                                                copy=False)
            dotted = diagonalOperator(x, spaces=spaces)
            return dotted.sum(spaces=spaces)
Theo Steininger's avatar
Theo Steininger committed
874

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

Theo Steininger's avatar
Theo Steininger committed
878
879
880
        Returns
        -------
        norm : scalar
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
881
            The L2-norm of the field values.
csongor's avatar
csongor committed
882
883

        """
Martin Reinecke's avatar
fixes    
Martin Reinecke committed
884
        return np.sqrt(np.abs(self.vdot(x=self)))
csongor's avatar
csongor committed
885
886

    def conjugate(self, inplace=False):
887
        """ Retruns the complex conjugate of the field.
Theo Steininger's avatar
Theo Steininger committed
888

889
890
891
        Parameters
        ----------
        inplace : boolean
Theo Steininger's avatar
Theo Steininger committed
892
            Decides whether the conjugation should be performed inplace.
Theo Steininger's avatar
Theo Steininger committed
893

894
895
896
897
        Returns
        -------
        cc : field
            The complex conjugated field.
csongor's avatar
csongor committed
898
899

        """
Theo Steininger's avatar
Theo Steininger committed
900

csongor's avatar
csongor committed
901
902
903
904
905
        if inplace:
            work_field = self
        else:
            work_field = self.copy_empty()

906
        new_val = self.get_val(copy=False)
Theo Steininger's avatar
Theo Steininger committed
907
        new_val = new_val.conjugate()
908
        work_field.set_val(new_val=new_val, copy=False)
csongor's avatar
csongor committed
909
910
911

        return work_field

Theo Steininger's avatar
Theo Steininger committed
912
    # ---General unary/contraction methods---
913

Theo Steininger's avatar
Theo Steininger committed
914
    def __pos__(self):
915
916
917
        """ x.__pos__() <==> +x

        Returns a (positive) copy of `self`.
Theo Steininger's avatar
Theo Steininger committed
918

919
        """
Theo Steininger's avatar
Theo Steininger committed
920

Theo Steininger's avatar
Theo Steininger committed
921
        return self.copy()
922

Theo Steininger's avatar
Theo Steininger committed
923
    def __neg__(self):
924
925
926
        """ x.__neg__() <==> -x

        Returns a negative copy of `self`.
Theo Steininger's avatar
Theo Steininger committed
927

928
        """
Theo Steininger's avatar
Theo Steininger committed
929

Theo Steininger's avatar
Theo Steininger committed
930
931
932
        return_field = self.copy_empty()
        new_val = -self.get_val(copy=False)
        return_field.set_val(new_val, copy=False)
csongor's avatar
csongor committed
933
934
        return return_field

Theo Steininger's avatar
Theo Steininger committed
935
    def __abs__(self):
936
937
938
        """ x.__abs__() <==> abs(x)

        Returns an absolute valued copy of `self`.
Theo Steininger's avatar
Theo Steininger committed
939

940
        """
Theo Steininger's avatar
Theo Steininger committed
941

Theo Steininger's avatar
Theo Steininger committed
942
        new_val = abs(self.get_val(copy=False))
943
        return_field = self.copy_empty(dtype=new_val.dtype)
Theo Steininger's avatar
Theo Steininger committed
944
945
        return_field.set_val(new_val, copy=False)
        return return_field
csongor's avatar
csongor committed
946

947
    def _contraction_helper(self, op, spaces):
Theo Steininger's avatar
Theo Steininger committed
948
949
        # build a list of all axes
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
950
            spaces = range(len(self.domain))
Theo Steininger's avatar
Theo Steininger committed
951
952
        else:
            spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
csongor's avatar
csongor committed
953

954
        axes_list = tuple(self.domain_axes[sp_index] for sp_index in spaces)
955
956

        try:
Theo Steininger's avatar
Theo Steininger committed
957
            axes_list = reduce(lambda x, y: x+y, axes_list)
958
        except TypeError:
Theo Steininger's avatar
Theo Steininger committed
959
            axes_list = ()
csongor's avatar
csongor committed
960

Martin Reinecke's avatar
stage1    
Martin Reinecke committed