field.py 42.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
Theo Steininger's avatar
Theo Steininger committed
13
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
cleanup    
Martin Reinecke committed
26
from keepers import Loggable
Jait Dixit's avatar
Jait Dixit committed
27

Martin Reinecke's avatar
Martin Reinecke committed
28
from .domain_object import DomainObject
29

Martin Reinecke's avatar
Martin Reinecke committed
30
from .spaces.power_space import PowerSpace
csongor's avatar
csongor committed
31

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

csongor's avatar
csongor committed
36

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

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

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

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
51
    val : scalar, numpy.ndarray, Field
52
53
54
        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
55

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

59
60
61
62
    copy: boolean

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

65
66
67
68
69
70
    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
71

72
73
74
75
76
77
78
    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
79

80
    """
81

Theo Steininger's avatar
Theo Steininger committed
82
    # ---Initialization methods---
83

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

Theo Steininger's avatar
Theo Steininger committed
88
        self.dtype = self._infer_dtype(dtype=dtype,
89
                                       val=val)
90

91
92
93
94
        if val is None:
            self._val = None
        else:
            self.set_val(new_val=val, copy=copy)
csongor's avatar
csongor committed
95

96
    def _parse_domain(self, domain, val=None):
97
        if domain is None:
98
99
100
101
            if isinstance(val, Field):
                domain = val.domain
            else:
                domain = ()
102
        elif isinstance(domain, DomainObject):
103
            domain = (domain,)
104
105
106
        elif not isinstance(domain, tuple):
            domain = tuple(domain)

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

Theo Steininger's avatar
Theo Steininger committed
114
115
116
117
118
119
120
121
122
123
    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)
124

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

139
140
        dtype = np.result_type(dtype, np.float)

Theo Steininger's avatar
Theo Steininger committed
141
        return dtype
142

143
    # ---Factory methods---
144

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

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

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

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

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

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

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

172
173

        """
Theo Steininger's avatar
Theo Steininger committed
174

175
        # create a initially empty field
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
176
        f = cls(domain=domain, dtype=dtype)
177
178
179
180
181
182
183

        # now use the processed input in terms of f in order to parse the
        # random arguments
        random_arguments = cls._parse_random_arguments(random_type=random_type,
                                                       f=f,
                                                       **kwargs)

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
184
        # extract the data from f and apply the appropriate
185
186
187
        # random number generator to it
        sample = f.get_val(copy=False)
        generator_function = getattr(Random, random_type)
188

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
189
190
191
        sample[:]=generator_function(dtype=f.dtype,
                                             shape=sample.shape,
                                             **random_arguments)
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
        return f

    @staticmethod
    def _parse_random_arguments(random_type, f, **kwargs):
        if random_type == "pm1":
            random_arguments = {}

        elif random_type == "normal":
            mean = kwargs.get('mean', 0)
            std = kwargs.get('std', 1)
            random_arguments = {'mean': mean,
                                'std': std}

        elif random_type == "uniform":
            low = kwargs.get('low', 0)
            high = kwargs.get('high', 1)
            random_arguments = {'low': low,
                                'high': high}

csongor's avatar
csongor committed
211
        else:
212
213
            raise KeyError(
                "unsupported random key '" + str(random_type) + "'.")
csongor's avatar
csongor committed
214

215
        return random_arguments
csongor's avatar
csongor committed
216

217
218
    # ---Powerspectral methods---

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

Theo Steininger's avatar
Theo Steininger committed
223
224
225
        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
226
        harmonic space. The resulting field has the same units as the initial
Theo Steininger's avatar
Theo Steininger committed
227
        field, corresponding to the square root of the power spectrum.
228
229
230

        Parameters
        ----------
Theo Steininger's avatar
Theo Steininger committed
231
232
233
234
235
        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
236
            if binbounds==None : bins are inferred.
237
238
239
240
241
242
243
244
245
246
        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
247

248
249
250
251
        Raise
        -----
        ValueError
            Raised if
Theo Steininger's avatar
Theo Steininger committed
252
253
                *len(domain) is != 1 when spaces==None
                *len(spaces) is != 1 if not None
254
                *the analyzed space is not harmonic
Theo Steininger's avatar
Theo Steininger committed
255

256
257
        Returns
        -------
Theo Steininger's avatar
Theo Steininger committed
258
        out : Field
Martin Reinecke's avatar
Martin Reinecke committed
259
            The output object. Its domain is a PowerSpace and it contains
260
261
262
263
264
            the power spectrum of 'self's field.

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

266
        """
Theo Steininger's avatar
Theo Steininger committed
267

Theo Steininger's avatar
Theo Steininger committed
268
        # check if all spaces in `self.domain` are either harmonic or
269
270
271
        # power_space instances
        for sp in self.domain:
            if not sp.harmonic and not isinstance(sp, PowerSpace):
Theo Steininger's avatar
Theo Steininger committed
272
                self.logger.info(
273
                    "Field has a space in `domain` which is neither "
274
275
276
                    "harmonic nor a PowerSpace.")

        # check if the `spaces` input is valid
277
278
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
279
            spaces = list(range(len(self.domain)))
280
281

        if len(spaces) == 0:
282
283
            raise ValueError(
                "No space for analysis specified.")
284

285
286
287
288
289
290
291
292
293
294
295
296
297
        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]
298
299

        for space_index in spaces:
300
301
            parts = [self._single_power_analyze(
                                work_field=part,
302
                                space_index=space_index,
303
304
                                binbounds=binbounds)
                     for part in parts]
305

306
307
308
309
310
311
        if keep_phase_information:
            result_field = parts[0] + 1j*parts[1]
        else:
            result_field = parts[0]

        return result_field
312
313

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

316
        if not work_field.domain[space_index].harmonic:
317
318
            raise ValueError(
                "The analyzed space must be harmonic.")
319

320
321
322
323
324
325
        # 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.

326
        harmonic_domain = work_field.domain[space_index]
327
        power_domain = PowerSpace(harmonic_partner=harmonic_domain,
Theo Steininger's avatar
Theo Steininger committed
328
                                  binbounds=binbounds)
329
330
        power_spectrum = cls._calculate_power_spectrum(
                                field_val=work_field.val,
Martin Reinecke's avatar
Martin Reinecke committed
331
                                pdomain=power_domain,
332
                                axes=work_field.domain_axes[space_index])
333
334

        # create the result field and put power_spectrum into it
335
        result_domain = list(work_field.domain)
336
        result_domain[space_index] = power_domain
337
        result_dtype = power_spectrum.dtype
338

339
        result_field = work_field.copy_empty(
340
                   domain=result_domain,
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
341
                   dtype=result_dtype)
342
343
344
345
        result_field.set_val(new_val=power_spectrum, copy=False)

        return result_field

346
    @classmethod
Martin Reinecke's avatar
Martin Reinecke committed
347
    def _calculate_power_spectrum(cls, field_val, pdomain, axes=None):
348

Martin Reinecke's avatar
Martin Reinecke committed
349
350
351
        pindex = pdomain.pindex
        # MR FIXME: how about iterating over slices, instead of replicating
        # pindex? Would save memory and probably isn't slower.
352
        if axes is not None:
353
354
355
356
            pindex = cls._shape_up_pindex(
                            pindex=pindex,
                            target_shape=field_val.shape,
                            axes=axes)
Theo Steininger's avatar
Theo Steininger committed
357

Martin Reinecke's avatar
Martin Reinecke committed
358
        power_spectrum = utilities.bincount_axis(pindex, weights=field_val,
359
                                         axis=axes)
Martin Reinecke's avatar
Martin Reinecke committed
360
        rho = pdomain.rho
361
362
363
364
365
366
367
368
        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

369
    @staticmethod
Martin Reinecke's avatar
Martin Reinecke committed
370
    def _shape_up_pindex(pindex, target_shape, axes):
Theo Steininger's avatar
Theo Steininger committed
371
        semiscaled_local_shape = [1, ] * len(target_shape)
Theo Steininger's avatar
Theo Steininger committed
372
        for i in range(len(axes)):
Martin Reinecke's avatar
Martin Reinecke committed
373
374
            semiscaled_local_shape[axes[i]] = pindex.shape[i]
        local_data = pindex
Theo Steininger's avatar
Theo Steininger committed
375
        semiscaled_local_data = local_data.reshape(semiscaled_local_shape)
Martin Reinecke's avatar
Martin Reinecke committed
376
377
        result_obj = np.empty(target_shape, dtype=pindex.dtype)
        result_obj[:] = semiscaled_local_data
378
379
380

        return result_obj

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

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

388
389
390
        Parameters
        ----------
        spaces : {tuple, int, None} *optional*
Theo Steininger's avatar
Theo Steininger committed
391
392
393
            Specifies the subspace containing all the PowerSpaces which
            should be converted (default : None).
            if spaces==None : Tries to convert the whole domain.
394
        real_power : boolean *optional*
Theo Steininger's avatar
Theo Steininger committed
395
396
            Determines whether the power spectrum is treated as intrinsically
            real or complex (default : True).
397
        real_signal : boolean *optional*
Theo Steininger's avatar
Theo Steininger committed
398
399
400
401
402
403
            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
404
        std : float *optional*
Theo Steininger's avatar
Theo Steininger committed
405
406
407
            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
408

409
410
411
412
        Returns
        -------
        out : Field
            The output object. A random field created with the power spectrum
Theo Steininger's avatar
Theo Steininger committed
413
            stored in the `spaces` in `self`.
414

Theo Steininger's avatar
Theo Steininger committed
415
416
417
418
419
420
        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.

421
422
423
        See Also
        --------
        power_analyze
Theo Steininger's avatar
Theo Steininger committed
424
425
426
427
428

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

429
        """
Theo Steininger's avatar
Theo Steininger committed
430

431
432
433
        # check if the `spaces` input is valid
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))

Theo Steininger's avatar
Theo Steininger committed
434
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
435
            spaces = list(range(len(self.domain)))
Theo Steininger's avatar
Theo Steininger committed
436

437
438
439
440
441
        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.")
442
443
444

        # create the result domain
        result_domain = list(self.domain)
445
446
        for power_space_index in spaces:
            power_space = self.domain[power_space_index]
447
            harmonic_domain = power_space.harmonic_partner
448
            result_domain[power_space_index] = harmonic_domain
449
450
451

        # create random samples: one or two, depending on whether the
        # power spectrum is real or complex
452
        if real_power:
453
            result_list = [None]
454
455
        else:
            result_list = [None, None]
456

457
458
        result_list = [self.__class__.from_random(
                             'normal',
459
460
461
                             mean=mean,
                             std=std,
                             domain=result_domain,
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
462
                             dtype=np.complex)
463
464
465
466
467
468
                       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
469

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
470
        spec = self.val.copy()
471
472
        spec = np.sqrt(spec)

473
474
475
476
477
478
479
        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
480
        result_val_list[0] *= local_rescaler.real
481
482

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

485
        if real_signal:
486
            result_val_list = [self._hermitian_decomposition(
487
488
489
490
491
                                            result_domain,
                                            result_val,
                                            spaces,
                                            result_list[0].domain_axes,
                                            preserve_gaussian_variance=True)[0]
492
                               for result_val in result_val_list]
493
494
495
496
497
498
499

        # 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]
500
501
502
            if not issubclass(result_val_list[0].dtype.type,
                              np.complexfloating):
                result = result.real
503
        else:
504
505
506
507
            result = result_list[0] + 1j*result_list[1]

        return result

508
    @staticmethod
509
510
    def _hermitian_decomposition(domain, val, spaces, domain_axes,
                                 preserve_gaussian_variance=False):
511
512
513
514
515
516

        flipped_val = val
        for space in spaces:
            flipped_val = domain[space].hermitianize_inverter(
                                                    x=flipped_val,
                                                    axes=domain_axes[space])
517
518
        # if no flips at all where performed `h` is a real field.
        # if all spaces use the default implementation of doing nothing when
Theo Steininger's avatar
Theo Steininger committed
519
        # no flips are applied, one can use `is` to infer this case.
520
521

        if flipped_val is val:
Martin Reinecke's avatar
Martin Reinecke committed
522
523
            h = flipped_val.real.copy()
            a = 1j * flipped_val.imag.copy()
524
525
526
527
        else:
            flipped_val = flipped_val.conjugate()
            h = (val + flipped_val)/2.
            a = val - h
528
529

        # correct variance
530
        if preserve_gaussian_variance:
Martin Reinecke's avatar
Martin Reinecke committed
531
532
            assert issubclass(val.dtype.type, np.complexfloating),\
                    "complex input field is needed here"
533
534
535
            h *= np.sqrt(2)
            a *= np.sqrt(2)

536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
#            The code below should not be needed in practice, since it would
#            only ever be called when hermitianizing a purely real field.
#            However it might be of educational use and keep us from forgetting
#            how these things are done ...

#            if not issubclass(val.dtype.type, np.complexfloating):
#                # in principle one must not correct the variance for the fixed
#                # points of the hermitianization. However, for a complex field
#                # the input field loses half of its power at its fixed points
#                # in the `hermitian` part. Hence, here a factor of sqrt(2) is
#                # also necessary!
#                # => The hermitianization can be done on a space level since
#                # either nothing must be done (LMSpace) or ALL points need a
#                # factor of sqrt(2)
#                # => use the preserve_gaussian_variance flag in the
#                # hermitian_decomposition method above.
#
#                # This code is for educational purposes:
#                fixed_points = [domain[i].hermitian_fixed_points()
#                                for i in spaces]
#                fixed_points = [[fp] if fp is None else fp
#                                for fp in fixed_points]
#
#                for product_point in itertools.product(*fixed_points):
#                    slice_object = np.array((slice(None), )*len(val.shape),
#                                            dtype=np.object)
#                    for i, sp in enumerate(spaces):
#                        point_component = product_point[i]
#                        if point_component is None:
#                            point_component = slice(None)
#                        slice_object[list(domain_axes[sp])] = point_component
#
#                    slice_object = tuple(slice_object)
#                    h[slice_object] /= np.sqrt(2)
#                    a[slice_object] /= np.sqrt(2)

572
573
        return (h, a)

574
575
    def _spec_to_rescaler(self, spec, result_list, power_space_index):
        power_space = self.domain[power_space_index]
576
577
578

        # weight the random fields with the power spectrum
        # therefore get the pindex from the power space
579
        pindex = power_space.pindex
580
581
582
583

        # 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
584
        local_pindex = pindex
585

586
587
588
589
590
        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
591
        # here, the power_spectrum is distributed into the new shape
592
593
        local_rescaler = spec[local_blow_up]
        return local_rescaler
594

Theo Steininger's avatar
Theo Steininger committed
595
    # ---Properties---
596

Theo Steininger's avatar
Theo Steininger committed
597
    def set_val(self, new_val=None, copy=False):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
598
        """ Sets the field's data object.
599
600
601

        Parameters
        ----------
Theo Steininger's avatar
Theo Steininger committed
602
        new_val : scalar, array-like, Field, None *optional*
603
604
            The values to be stored in the field.
            {default : None}
Theo Steininger's avatar
Theo Steininger committed
605

606
        copy : boolean, *optional*
Theo Steininger's avatar
Theo Steininger committed
607
608
            If False, Field tries to not copy the input data but use it
            directly.
609
610
611
612
613
614
            {default : False}
        See Also
        --------
        val

        """
Theo Steininger's avatar
Theo Steininger committed
615

616
617
        new_val = self.cast(new_val)
        if copy:
Theo Steininger's avatar
Theo Steininger committed
618
619
            new_val = new_val.copy()
        self._val = new_val
620
        return self
csongor's avatar
csongor committed
621

622
    def get_val(self, copy=False):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
623
        """ Returns the data object associated with this Field.
624
625
626
627

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

631
632
        Returns
        -------
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
633
        out : numpy.ndarray
634
635
636
637
638
639

        See Also
        --------
        val

        """
Theo Steininger's avatar
Theo Steininger committed
640

641
642
643
        if self._val is None:
            self.set_val(None)

644
        if copy:
Theo Steininger's avatar
Theo Steininger committed
645
            return self._val.copy()
646
        else:
Theo Steininger's avatar
Theo Steininger committed
647
            return self._val
csongor's avatar
csongor committed
648

Theo Steininger's avatar
Theo Steininger committed
649
650
    @property
    def val(self):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
651
        """ Returns the data object associated with this Field.
Theo Steininger's avatar
Theo Steininger committed
652

653
654
        Returns
        -------
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
655
        out : numpy.ndarray
656
657
658
659
660
661

        See Also
        --------
        get_val

        """
Theo Steininger's avatar
Theo Steininger committed
662

663
        return self.get_val(copy=False)
csongor's avatar
csongor committed
664

Theo Steininger's avatar
Theo Steininger committed
665
666
    @val.setter
    def val(self, new_val):
667
        self.set_val(new_val=new_val, copy=False)
csongor's avatar
csongor committed
668

669
670
    @property
    def shape(self):
Theo Steininger's avatar
Theo Steininger committed
671
        """ Returns the total shape of the Field's data array.
Theo Steininger's avatar
Theo Steininger committed
672

673
674
675
        Returns
        -------
        out : tuple
Martin Reinecke's avatar
Martin Reinecke committed
676
            The output object. The tuple contains the dimensions of the spaces
677
678
679
680
681
682
683
            in domain.

        See Also
        --------
        dim

        """
Theo Steininger's avatar
Theo Steininger committed
684
685
686
687
688
689
690
691
        if not hasattr(self, '_shape'):
            shape_tuple = tuple(sp.shape for sp in self.domain)
            try:
                global_shape = reduce(lambda x, y: x + y, shape_tuple)
            except TypeError:
                global_shape = ()
            self._shape = global_shape
        return self._shape
csongor's avatar
csongor committed
692

693
694
    @property
    def dim(self):
Theo Steininger's avatar
Theo Steininger committed
695
        """ Returns the total number of pixel-dimensions the field has.
Theo Steininger's avatar
Theo Steininger committed
696

Theo Steininger's avatar
Theo Steininger committed
697
        Effectively, all values from shape are multiplied.
Theo Steininger's avatar
Theo Steininger committed
698

699
700
701
702
703
704
705
706
707
708
        Returns
        -------
        out : int
            The dimension of the Field.

        See Also
        --------
        shape

        """
Theo Steininger's avatar
Theo Steininger committed
709

710
        dim_tuple = tuple(sp.dim for sp in self.domain)
Theo Steininger's avatar
Theo Steininger committed
711
        try:
Martin Reinecke's avatar
Martin Reinecke committed
712
            return int(reduce(lambda x, y: x * y, dim_tuple))
Theo Steininger's avatar
Theo Steininger committed
713
714
        except TypeError:
            return 0
csongor's avatar
csongor committed
715

716
717
    @property
    def dof(self):
Theo Steininger's avatar
Theo Steininger committed
718
719
720
721
722
723
        """ 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
724
725
726
727
728
729
730
        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
731
732
733
        """ Returns the total volume of all spaces in the domain.
        """

Theo Steininger's avatar
Theo Steininger committed
734
        volume_tuple = tuple(sp.total_volume for sp in self.domain)
735
        try:
Theo Steininger's avatar
Theo Steininger committed
736
            return reduce(lambda x, y: x * y, volume_tuple)
737
        except TypeError:
Theo Steininger's avatar
Theo Steininger committed
738
            return 0.
739

Theo Steininger's avatar
Theo Steininger committed
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
    @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
758
    # ---Special unary/binary operations---
759

csongor's avatar
csongor committed
760
    def cast(self, x=None, dtype=None):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
761
        """ Transforms x to an object with the correct dtype and shape.
Theo Steininger's avatar
Theo Steininger committed
762

763
764
        Parameters
        ----------
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
765
766
767
        x : scalar, numpy.ndarray, Field, array_like
            The input that shall be casted on a numpy.ndarray of the same shape
            like the domain.
Theo Steininger's avatar
Theo Steininger committed
768

769
        dtype : type
Theo Steininger's avatar
Theo Steininger committed
770
            The datatype the output shall have. This can be used to override
Martin Reinecke's avatar
typos    
Martin Reinecke committed
771
            the field's dtype.
Theo Steininger's avatar
Theo Steininger committed
772

773
774
        Returns
        -------
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
775
        out : numpy.ndarray
776
777
778
779
780
781
782
            The output object.

        See Also
        --------
        _actual_cast

        """
csongor's avatar
csongor committed
783
784
        if dtype is None:
            dtype = self.dtype
785
786
        else:
            dtype = np.dtype(dtype)
787

788
789
        casted_x = x

790
        for ind, sp in enumerate(self.domain):
791
            casted_x = sp.pre_cast(casted_x,
792
793
794
                                   axes=self.domain_axes[ind])

        casted_x = self._actual_cast(casted_x, dtype=dtype)
795
796

        for ind, sp in enumerate(self.domain):
797
798
            casted_x = sp.post_cast(casted_x,
                                    axes=self.domain_axes[ind])
799

800
        return casted_x
csongor's avatar
csongor committed
801

Theo Steininger's avatar
Theo Steininger committed
802
    def _actual_cast(self, x, dtype=None):
803
        if isinstance(x, Field):
csongor's avatar
csongor committed
804
805
806
807
            x = x.get_val()

        if dtype is None:
            dtype = self.dtype
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
808
        if x is not None:
Martin Reinecke's avatar
more    
Martin Reinecke committed
809
810
            if np.isscalar(x):
                return np.full(self.shape,x, dtype=dtype)
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
811
812
813
            return np.asarray(x, dtype=dtype).reshape(self.shape)
        else:
            return np.empty(self.shape, dtype=dtype)
csongor's avatar
csongor committed
814

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

818
819
        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
820
        able to define the domain and the dtype of the returned Field.
821
822
823
824
825

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

827
828
        dtype : type
            The new dtype the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
829

830
831
832
833
834
835
836
837
838
839
        Returns
        -------
        out : Field
            The output object. An identical copy of 'self'.

        See Also
        --------
        copy_empty

        """
Theo Steininger's avatar
Theo Steininger committed
840

Theo Steininger's avatar
Theo Steininger committed
841
        copied_val = self.get_val(copy=True)
842
843
        new_field = self.copy_empty(
                                domain=domain,
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
844
                                dtype=dtype)
Theo Steininger's avatar
Theo Steininger committed
845
846
        new_field.set_val(new_val=copied_val, copy=False)
        return new_field
csongor's avatar
csongor committed
847

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
848
    def copy_empty(self, domain=None, dtype=None):
849
850
851
        """ 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
852
853
854
        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
855
        to change the domain and the dtype of the returned Field.
Theo Steininger's avatar
Theo Steininger committed
856

857
858
859
860
        Parameters
        ----------
        domain : DomainObject
            The new domain the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
861

862
863
        dtype : type
            The new dtype the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
864

865
866
867
        Returns
        -------
        out : Field
Theo Steininger's avatar
Theo Steininger committed
868
            The output object.
869
870
871
872
873
874

        See Also
        --------
        copy

        """
Theo Steininger's avatar
Theo Steininger committed
875

Theo Steininger's avatar
Theo Steininger committed
876
877
        if domain is None:
            domain = self.domain
csongor's avatar
csongor committed
878
        else:
Theo Steininger's avatar
Theo Steininger committed
879
            domain = self._parse_domain(domain)
csongor's avatar
csongor committed
880

Theo Steininger's avatar
Theo Steininger committed
881
882
883
884
        if dtype is None:
            dtype = self.dtype
        else:
            dtype = np.dtype(dtype)
csongor's avatar
csongor committed
885

Theo Steininger's avatar
Theo Steininger committed
886
887
        fast_copyable = True
        try:
Martin Reinecke's avatar
Martin Reinecke committed
888
            for i in range(len(self.domain)):
Theo Steininger's avatar
Theo Steininger committed
889
890
891
892
893
894
                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
895
        if (fast_copyable and dtype == self.dtype):
Theo Steininger's avatar
Theo Steininger committed
896
897
            new_field = self._fast_copy_empty()
        else:
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
898
            new_field = Field(domain=domain, dtype=dtype)
Theo Steininger's avatar
Theo Steininger committed
899
        return new_field
csongor's avatar
csongor committed
900

Theo Steininger's avatar
Theo Steininger committed
901
902
903
904
905
906
    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
907
        for key, value in list(self.__dict__.items()):
908
            if key != '_val':
Theo Steininger's avatar
Theo Steininger committed
909
910
                new_field.__dict__[key] = value
            else:
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
911
                new_field.__dict__[key] = np.empty_like(self.val)
Theo Steininger's avatar
Theo Steininger committed
912
913
914
        return new_field

    def weight(self, power=1, inplace=False, spaces=None):
Theo Steininger's avatar
Theo Steininger committed
915
        """ Weights the pixels of `self` with their invidual pixel-volume.
916
917
918
919

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

922
        inplace : boolean
Theo Steininger's avatar
Theo Steininger committed
923
924
            If True, `self` will be weighted and returned. Otherwise, a copy
            is made.
Theo Steininger's avatar
Theo Steininger committed
925

Theo Steininger's avatar
Theo Steininger committed
926
927
        spaces : tuple of ints
            Determines on which subspace the operation takes place.
Theo Steininger's avatar
Theo Steininger committed
928

929
930
931
        Returns
        -------
        out : Field
Theo Steininger's avatar
Theo Steininger committed
932
            The weighted field.
933
934

        """
935
        if inplace:
csongor's avatar
csongor committed
936
937
938
939
            new_field = self
        else:
            new_field = self.copy_empty()

940
        new_val = self.get_val(copy=False)
csongor's avatar
csongor committed
941

942
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
csongor's avatar
csongor committed
943
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
944
            spaces = list(range(len(self.domain)))
csongor's avatar
csongor committed
945

946
        for ind, sp in enumerate(self.domain):
Theo Steininger's avatar
Theo Steininger committed
947
948
949
950
951
            if ind in spaces:
                new_val = sp.weight(new_val,
                                    power=power,
                                    axes=self.domain_axes[ind],
                                    inplace=inplace)
952
953

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

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

959
960
961
        Parameters
        ----------
        x : Field
Theo Steininger's avatar
Theo Steininger committed
962
            The domain of x must contain `self.domain`
Theo Steininger's avatar
Theo Steininger committed
963

Theo Steininger's avatar
Theo Steininger committed
964
965
966
        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
967

968
        bare : boolean
Theo Steininger's avatar
Theo Steininger committed
969
            If true, no volume factors will be included in the computation.
Theo Steininger's avatar
Theo Steininger committed
970

971
972
973
        Returns
        -------
        out : float, complex
Theo Steininger's avatar
Theo Steininger committed
974

975
        """
976
977
978
        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
979

Martin Reinecke's avatar
Martin Reinecke committed
980
        # Compute the dot respecting the fact of discrete/continuous spaces
Theo Steininger's avatar
Theo Steininger committed
981
982
983
984
985
        if bare:
            y = self
        else:
            y = self.weight(power=1)

986
987
988
        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
989
            result = (y_val.conjugate() * x_val).sum()
990
991
992
993
            return result
        else:
            # create a diagonal operator which is capable of taking care of the
            # axes-matching
994
            from .operators.diagonal_operator import DiagonalOperator
995
996
997
998
999
1000
            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
1001

Theo Steininger's avatar
Theo Steininger committed
1002
    def norm(self):
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
1003