field.py 45.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# NIFTy
# Copyright (C) 2017  Theo Steininger
#
# Author: Theo Steininger
#
# 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/>.

csongor's avatar
csongor committed
19
20
21
from __future__ import division
import numpy as np

Theo Steininger's avatar
Theo Steininger committed
22
23
from keepers import Versionable,\
                    Loggable
Jait Dixit's avatar
Jait Dixit committed
24

25
from d2o import distributed_data_object,\
26
    STRATEGIES as DISTRIBUTION_STRATEGIES
csongor's avatar
csongor committed
27

28
from nifty.config import nifty_configuration as gc
csongor's avatar
csongor committed
29

30
from nifty.domain_object import DomainObject
31

32
from nifty.spaces.power_space import PowerSpace
csongor's avatar
csongor committed
33

csongor's avatar
csongor committed
34
import nifty.nifty_utilities as utilities
35
36
from nifty.random import Random

csongor's avatar
csongor committed
37

Jait Dixit's avatar
Jait Dixit committed
38
class Field(Loggable, Versionable, object):
Theo Steininger's avatar
Theo Steininger committed
39
40
41
    """ The discrete representation of a continuous field over multiple spaces.

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

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

52
53
54
55
    val : scalar, numpy.ndarray, distributed_data_object, Field
        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
56

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

60
61
62
63
64
65
    distribution_strategy: optional[{'fftw', 'equal', 'not', 'freeform'}]
        Specifies which distributor will be created and used.
        'fftw'      uses the distribution strategy of pyfftw,
        'equal'     tries to  distribute the data as uniform as possible
        'not'       does not distribute the data at all
        'freeform'  distribute the data according to the given local data/shape
Theo Steininger's avatar
Theo Steininger committed
66

67
68
69
70
71
    copy: boolean

    Attributes
    ----------
    val : distributed_data_object
Theo Steininger's avatar
Theo Steininger committed
72

73
74
75
76
77
78
79
    domain : DomainObject
        See Parameters.
    domain_axes : tuple of tuples
        Enumerates the axes of the Field
    dtype : type
        Contains the datatype stored in the Field.
    distribution_strategy : string
Theo Steininger's avatar
Theo Steininger committed
80
81
        Name of the used distribution_strategy.

82
83
84
85
86
87
88
    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
89

90
91
92
93
94
95
96
97
98
99
100
    Examples
    --------
    >>> a = Field(RGSpace([4,5]),val=2)
    >>> a.val
    <distributed_data_object>
    array([[2, 2, 2, 2, 2],
           [2, 2, 2, 2, 2],
           [2, 2, 2, 2, 2],
           [2, 2, 2, 2, 2]])
    >>> a.dtype
    dtype('int64')
Theo Steininger's avatar
Theo Steininger committed
101

102
103
104
105
106
    See Also
    --------
    distributed_data_object

    """
Theo Steininger's avatar
Theo Steininger committed
107
    # ---Initialization methods---
108

109
    def __init__(self, domain=None, val=None, dtype=None,
110
                 distribution_strategy=None, copy=False):
csongor's avatar
csongor committed
111

112
        self.domain = self._parse_domain(domain=domain, val=val)
113
        self.domain_axes = self._get_axes_tuple(self.domain)
csongor's avatar
csongor committed
114

Theo Steininger's avatar
Theo Steininger committed
115
        self.dtype = self._infer_dtype(dtype=dtype,
116
                                       val=val)
117

118
119
120
        self.distribution_strategy = self._parse_distribution_strategy(
                                distribution_strategy=distribution_strategy,
                                val=val)
csongor's avatar
csongor committed
121

122
123
124
125
        if val is None:
            self._val = None
        else:
            self.set_val(new_val=val, copy=copy)
csongor's avatar
csongor committed
126

127
    def _parse_domain(self, domain, val=None):
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
        """ Returns a tuple of DomainObjects for nomenclature unification.

        Parameters
        ----------
        domain : all supported NIFTY spaces
            The domain over which the Field lives.
        val : a NIFTY Field instance
            Can be used to make Field infere it's domain by adopting val's
            domain.

        Returns
        -------
        out : tuple
            The output object. A tuple with one or multiple DomainObjects.
        """
143
        if domain is None:
144
145
146
147
            if isinstance(val, Field):
                domain = val.domain
            else:
                domain = ()
148
        elif isinstance(domain, DomainObject):
149
            domain = (domain,)
150
151
152
        elif not isinstance(domain, tuple):
            domain = tuple(domain)

csongor's avatar
csongor committed
153
        for d in domain:
154
            if not isinstance(d, DomainObject):
155
156
                raise TypeError(
                    "Given domain contains something that is not a "
157
                    "DomainObject instance.")
csongor's avatar
csongor committed
158
159
        return domain

Theo Steininger's avatar
Theo Steininger committed
160
    def _get_axes_tuple(self, things_with_shape, start=0):
161
162
163
164
165
166
167
168
169
170
        """ Enumerates all axes of the domain.

        This function is used in the greater context of the 'spaces' keyword.

        Parameters
        ----------
        things_with_shape : indexable list of objects with .shape property
            Normal input is a domain/ tuple of domains.
        start : int
            Sets the integer number for the first axis
Theo Steininger's avatar
Theo Steininger committed
171

172
173
174
175
        Returns
        -------
        out : tuple
            Incremental numeration of all axes.
Theo Steininger's avatar
Theo Steininger committed
176

177
178
        Note
        ----
Theo Steininger's avatar
Theo Steininger committed
179

180
181
182
183
        The 'spaces' keyword is used in operators in order to carry out
        operations only on a certain subspace if the domain of the Field is
        a product space.
        """
Theo Steininger's avatar
Theo Steininger committed
184
185
186
187
188
189
190
191
192
        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)
193

194
    def _infer_dtype(self, dtype, val):
195
196
197
198
199
200
201
202
203
        """ Inferes the datatype of the Field

        Parameters
        ----------
        dtype : type
            Can be None
        val : list of arrays
            If the dtype is None, Fields tries to infere the datatype from the
            values given to it at initialization.
Theo Steininger's avatar
Theo Steininger committed
204

205
206
207
208
        Returns
        -------
        out : np.dtype
        """
csongor's avatar
csongor committed
209
        if dtype is None:
210
            try:
211
                dtype = val.dtype
212
            except AttributeError:
Theo Steininger's avatar
Theo Steininger committed
213
214
215
                try:
                    if val is None:
                        raise TypeError
216
                    dtype = np.result_type(val)
Theo Steininger's avatar
Theo Steininger committed
217
                except(TypeError):
218
                    dtype = np.dtype(gc['default_field_dtype'])
Theo Steininger's avatar
Theo Steininger committed
219
        else:
220
            dtype = np.dtype(dtype)
221

Theo Steininger's avatar
Theo Steininger committed
222
        return dtype
223

224
225
    def _parse_distribution_strategy(self, distribution_strategy, val):
        if distribution_strategy is None:
226
            if isinstance(val, distributed_data_object):
227
                distribution_strategy = val.distribution_strategy
228
            elif isinstance(val, Field):
229
                distribution_strategy = val.distribution_strategy
230
            else:
231
                self.logger.debug("distribution_strategy set to default!")
232
                distribution_strategy = gc['default_distribution_strategy']
233
        elif distribution_strategy not in DISTRIBUTION_STRATEGIES['global']:
234
235
236
            raise ValueError(
                    "distribution_strategy must be a global-type "
                    "strategy.")
237
        return distribution_strategy
238
239

    # ---Factory methods---
240

241
    @classmethod
242
    def from_random(cls, random_type, domain=None, dtype=None,
243
                    distribution_strategy=None, **kwargs):
244
245
246
247
248
        """ Draws a random field with the given parameters.

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

250
251
252
        random_type : String
            'pm1', 'normal', 'uniform' are the supported arguments for this
            method.
Theo Steininger's avatar
Theo Steininger committed
253

254
255
        domain : DomainObject
            The domain of the output random field
Theo Steininger's avatar
Theo Steininger committed
256

257
258
        dtype : type
            The datatype of the output random field
Theo Steininger's avatar
Theo Steininger committed
259

260
261
        distribution_strategy : all supported distribution strategies
            The distribution strategy of the output random field
Theo Steininger's avatar
Theo Steininger committed
262

263
264
265
266
267
268
269
270
271
272
        Returns
        -------
        out : Field
            The output object.

        See Also
        --------
        _parse_random_arguments, power_synthesise

        """
273
        # create a initially empty field
274
        f = cls(domain=domain, dtype=dtype,
275
                distribution_strategy=distribution_strategy)
276
277
278
279
280
281
282

        # 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
Martin Reinecke committed
283
        # extract the distributed_data_object from f and apply the appropriate
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
        # random number generator to it
        sample = f.get_val(copy=False)
        generator_function = getattr(Random, random_type)
        sample.apply_generator(
            lambda shape: generator_function(dtype=f.dtype,
                                             shape=shape,
                                             **random_arguments))
        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
310
        else:
311
312
            raise KeyError(
                "unsupported random key '" + str(random_type) + "'.")
csongor's avatar
csongor committed
313

314
        return random_arguments
csongor's avatar
csongor committed
315

316
317
    # ---Powerspectral methods---

318
319
320
    def power_analyze(self, spaces=None, log=False, nbin=None, binbounds=None,
                      real_signal=True):
        """ Computes the powerspectrum of the Field
Theo Steininger's avatar
Theo Steininger committed
321
322

        Creates a PowerSpace with the given attributes and computes the
323
324
325
326
327
328
329
330
331
332
        power spectrum as a field over this PowerSpace.
        It's important to note that this can only be done if the subspace to
        be analyzed is in harmonic space.

        Parameters
        ----------
        spaces : int, *optional*
            The subspace which you want to have the powerspectrum of.
            {default : None}
                if spaces==None : Tries to synthesize for the whole domain
Theo Steininger's avatar
Theo Steininger committed
333

334
335
336
        log : boolean, *optional*
            True if the output PowerSpace should have log binning.
            {default : False}
Theo Steininger's avatar
Theo Steininger committed
337

338
339
340
341
        nbin : int, None, *optional*
            The number of bins the resulting PowerSpace shall have.
            {default : None}
                if nbin==None : maximum number of bins is used
Theo Steininger's avatar
Theo Steininger committed
342

343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
        binbounds : array-like, None, *optional*
            Inner bounds of the bins, if specifield
            {default : None}
                if binbounds==None : bins are inferred. Overwrites nbins and log
        real_signal : boolean, *optional*
            Whether the analysed signal-space Field is real or complex.
            For a real field a complex power spectrum comes out.
            For a compex field all power is put in a real power spectrum.
            {default : True}
        Raise
        -----
        ValueError
            Raised if
                *len(spaces) is either 0 or >1
                *len(domain) is not 1 with spaces=None
                *the analyzed space is not harmonic
Theo Steininger's avatar
Theo Steininger committed
359

360
361
        Returns
        -------
Theo Steininger's avatar
Theo Steininger committed
362
        out : Field
363
364
365
366
367
368
369
            The output object. It's domain is a PowerSpace and it contains
            the power spectrum of 'self's field.

        See Also
        --------
        power_synthesize, PowerSpace
        """
Theo Steininger's avatar
Theo Steininger committed
370
        # check if all spaces in `self.domain` are either harmonic or
371
372
373
        # power_space instances
        for sp in self.domain:
            if not sp.harmonic and not isinstance(sp, PowerSpace):
Theo Steininger's avatar
Theo Steininger committed
374
                self.logger.info(
375
                    "Field has a space in `domain` which is neither "
376
377
378
                    "harmonic nor a PowerSpace.")

        # check if the `spaces` input is valid
379
380
381
382
383
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
        if spaces is None:
            if len(self.domain) == 1:
                spaces = (0,)
            else:
384
385
386
                raise ValueError(
                    "Field has multiple spaces as domain "
                    "but `spaces` is None.")
387
388

        if len(spaces) == 0:
389
390
            raise ValueError(
                "No space for analysis specified.")
391
        elif len(spaces) > 1:
392
393
            raise ValueError(
                "Conversion of only one space at a time is allowed.")
394
395
396
397

        space_index = spaces[0]

        if not self.domain[space_index].harmonic:
398
399
            raise ValueError(
                "The analyzed space must be harmonic.")
400

401
402
403
404
405
406
        # 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.

407
408
409
410
        distribution_strategy = \
            self.val.get_axes_local_distribution_strategy(
                self.domain_axes[space_index])

411
412
        harmonic_domain = self.domain[space_index]
        power_domain = PowerSpace(harmonic_domain=harmonic_domain,
413
                                  distribution_strategy=distribution_strategy,
414
                                  log=log, nbin=nbin, binbounds=binbounds)
415

416
        # extract pindex and rho from power_domain
417
418
        pindex = power_domain.pindex
        rho = power_domain.rho
419

420
        if real_signal:
421
            hermitian_part, anti_hermitian_part = \
422
                harmonic_domain.hermitian_decomposition(
423
424
425
426
427
428
429
430
431
432
433
434
435
436
                                            self.val,
                                            axes=self.domain_axes[space_index])

            [hermitian_power, anti_hermitian_power] = \
                [self._calculate_power_spectrum(
                                            x=part,
                                            pindex=pindex,
                                            rho=rho,
                                            axes=self.domain_axes[space_index])
                 for part in [hermitian_part, anti_hermitian_part]]

            power_spectrum = hermitian_power + 1j * anti_hermitian_power
        else:
            power_spectrum = self._calculate_power_spectrum(
437
438
439
440
441
442
443
444
445
                                            x=self.val,
                                            pindex=pindex,
                                            rho=rho,
                                            axes=self.domain_axes[space_index])

        # create the result field and put power_spectrum into it
        result_domain = list(self.domain)
        result_domain[space_index] = power_domain

446
        if real_signal:
447
448
449
450
            result_dtype = np.complex
        else:
            result_dtype = np.float

451
452
        result_field = self.copy_empty(
                   domain=result_domain,
453
                   dtype=result_dtype,
454
                   distribution_strategy=power_spectrum.distribution_strategy)
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
        result_field.set_val(new_val=power_spectrum, copy=False)

        return result_field

    def _calculate_power_spectrum(self, x, pindex, rho, axes=None):
        fieldabs = abs(x)
        fieldabs **= 2

        if axes is not None:
            pindex = self._shape_up_pindex(
                                    pindex=pindex,
                                    target_shape=x.shape,
                                    target_strategy=x.distribution_strategy,
                                    axes=axes)
        power_spectrum = pindex.bincount(weights=fieldabs,
                                         axis=axes)
        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

        power_spectrum **= 0.5
        return power_spectrum

    def _shape_up_pindex(self, pindex, target_shape, target_strategy, axes):
        if pindex.distribution_strategy not in \
                DISTRIBUTION_STRATEGIES['global']:
483
            raise ValueError("pindex's distribution strategy must be "
484
485
486
487
488
489
                             "global-type")

        if pindex.distribution_strategy in DISTRIBUTION_STRATEGIES['slicing']:
            if ((0 not in axes) or
                    (target_strategy is not pindex.distribution_strategy)):
                raise ValueError(
490
                    "A slicing distributor shall not be reshaped to "
491
492
493
494
495
496
497
498
499
500
501
502
503
                    "something non-sliced.")

        semiscaled_shape = [1, ] * len(target_shape)
        for i in axes:
            semiscaled_shape[i] = target_shape[i]
        local_data = pindex.get_local_data(copy=False)
        semiscaled_local_data = local_data.reshape(semiscaled_shape)
        result_obj = pindex.copy_empty(global_shape=target_shape,
                                       distribution_strategy=target_strategy)
        result_obj.set_full_data(semiscaled_local_data, copy=False)

        return result_obj

504
505
506
    def power_synthesize(self, spaces=None, real_power=True, real_signal=True,
                         mean=None, std=None):
        """Yields a random field in harmonic space with this power spectrum.
Theo Steininger's avatar
Theo Steininger committed
507

508
509
        This method draws a Gaussian random field in the harmic partner domain.
        The drawn field has this field as its power spectrum.
Theo Steininger's avatar
Theo Steininger committed
510

511
512
513
        Notes
        -----
        For this the domain must be a PowerSpace.
Theo Steininger's avatar
Theo Steininger committed
514

515
516
517
518
        Parameters
        ----------
        spaces : {tuple, int, None} *optional*
            Specifies the subspace in which the power will be synthesized in
Theo Steininger's avatar
Theo Steininger committed
519
            case of a product space.
520
521
            {default : None}
                if spaces==None : Tries to synthesize for the whole domain
Theo Steininger's avatar
Theo Steininger committed
522

523
524
525
        real_power : boolean *optional*
            Determines whether the power spectrum is real or complex
            {default : True}
Theo Steininger's avatar
Theo Steininger committed
526

527
528
529
530
531
532
533
534
535
536
        real_signal : boolean *optional*
            True will result in a purely real signal-space field.
            This means that the created field is symmetric wrt. the origin
            after complex conjugation.
            {default : True}
        mean : {float, None} *optional*
            The mean of the noise field the powerspectrum will be multiplied on.
            {default : None}
                if mean==None : mean will be set to 0
        std : float *optional*
Theo Steininger's avatar
Theo Steininger committed
537
            The standard deviation of the noise field the powerspectrum will be
538
539
540
            multiplied on.
            {default : None}
                if std==None : std will be set to 1
Theo Steininger's avatar
Theo Steininger committed
541

542
543
544
545
546
547
548
549
550
        Returns
        -------
        out : Field
            The output object. A random field created with the power spectrum
            stored in 'self'

        See Also
        --------
        power_analyze
551

552
        """
553
554
555
        # check if the `spaces` input is valid
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))

Theo Steininger's avatar
Theo Steininger committed
556
557
558
        if spaces is None:
            spaces = range(len(self.domain))

559
560
561
562
563
        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.")
564
565
566

        # create the result domain
        result_domain = list(self.domain)
567
568
        for power_space_index in spaces:
            power_space = self.domain[power_space_index]
569
570
            harmonic_domain = power_space.harmonic_domain
            result_domain[power_space_index] = harmonic_domain
571
572
573

        # create random samples: one or two, depending on whether the
        # power spectrum is real or complex
574
        if real_power:
575
            result_list = [None]
576
577
        else:
            result_list = [None, None]
578

579
580
        result_list = [self.__class__.from_random(
                             'normal',
581
582
583
                             mean=mean,
                             std=std,
                             domain=result_domain,
584
                             dtype=np.complex,
585
                             distribution_strategy=self.distribution_strategy)
586
587
588
589
590
591
                       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
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609

        spec = self.val.get_full_data()
        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
        result_val_list[0].apply_scalar_function(
                                            lambda x: x * local_rescaler.real,
                                            inplace=True)

        if not real_power:
            result_val_list[1].apply_scalar_function(
                                            lambda x: x * local_rescaler.imag,
                                            inplace=True)

610
        if real_signal:
611
            for power_space_index in spaces:
612
613
                harmonic_domain = result_domain[power_space_index]
                result_val_list = [harmonic_domain.hermitian_decomposition(
614
615
616
617
618
                                    result_val,
                                    axes=result.domain_axes[power_space_index],
                                    preserve_gaussian_variance=True)[0]
                                   for (result, result_val)
                                   in zip(result_list, result_val_list)]
619
620
621
622
623
624
625

        # 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]
626
        else:
627
628
629
630
631
632
            result = result_list[0] + 1j*result_list[1]

        return result

    def _spec_to_rescaler(self, spec, result_list, power_space_index):
        power_space = self.domain[power_space_index]
633
634
635

        # weight the random fields with the power spectrum
        # therefore get the pindex from the power space
636
        pindex = power_space.pindex
637
638
639
640
641
642
643
        # take the local data from pindex. This data must be compatible to the
        # local data of the field given the slice of the PowerSpace
        local_distribution_strategy = \
            result_list[0].val.get_axes_local_distribution_strategy(
                result_list[0].domain_axes[power_space_index])

        if pindex.distribution_strategy is not local_distribution_strategy:
644
            self.logger.warn(
645
                "The distribution_stragey of pindex does not fit the "
646
647
648
649
650
651
652
653
654
655
                "slice_local distribution strategy of the synthesized field.")

        # 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
        local_pindex = pindex.get_local_data(copy=False)

        local_blow_up = [slice(None)]*len(self.shape)
        local_blow_up[self.domain_axes[power_space_index][0]] = local_pindex
        # here, the power_spectrum is distributed into the new shape
656
657
        local_rescaler = spec[local_blow_up]
        return local_rescaler
658

Theo Steininger's avatar
Theo Steininger committed
659
    # ---Properties---
660

Theo Steininger's avatar
Theo Steininger committed
661
    def set_val(self, new_val=None, copy=False):
662
663
664
665
        """ Let's one set the values of the.

        Parameters
        ----------
Theo Steininger's avatar
Theo Steininger committed
666
        new_val : number, numpy.array, distributed_data_object,
667
668
669
670
                Field, None, *optional*
            The values to be stored in the field.
            {default : None}
                if new_val==None : sets the values to 0.
Theo Steininger's avatar
Theo Steininger committed
671

672
673
674
675
676
677
678
679
680
        copy : boolean, *optional*
            True if this field holds a copy of new_val, False if
            it holds the same object
            {default : False}
        See Also
        --------
        val

        """
681
682
        new_val = self.cast(new_val)
        if copy:
Theo Steininger's avatar
Theo Steininger committed
683
684
            new_val = new_val.copy()
        self._val = new_val
685
        return self
csongor's avatar
csongor committed
686

687
    def get_val(self, copy=False):
688
689
690
691
692
        """ Acceses the values stored in the Field.

        Parameters
        ----------
        copy : boolean
Theo Steininger's avatar
Theo Steininger committed
693
            True makes the method retrun a COPY of the Field's underlying
694
            distributed_data_object.
Theo Steininger's avatar
Theo Steininger committed
695

696
697
698
699
700
701
702
703
704
705
        Returns
        -------
        out : distributed_data_object
            The output object.

        See Also
        --------
        val

        """
706
707
708
        if self._val is None:
            self.set_val(None)

709
        if copy:
Theo Steininger's avatar
Theo Steininger committed
710
            return self._val.copy()
711
        else:
Theo Steininger's avatar
Theo Steininger committed
712
            return self._val
csongor's avatar
csongor committed
713

Theo Steininger's avatar
Theo Steininger committed
714
715
    @property
    def val(self):
716
        """ Retruns actual distributed_data_object associated with this Field.
Theo Steininger's avatar
Theo Steininger committed
717

718
719
720
721
722
723
724
725
726
727
        Returns
        -------
        out : distributed_data_object
            The output object.

        See Also
        --------
        get_val

        """
728
        return self.get_val(copy=False)
csongor's avatar
csongor committed
729

Theo Steininger's avatar
Theo Steininger committed
730
731
    @val.setter
    def val(self, new_val):
732
        """ Sets the values in the d2o of the Field.
Theo Steininger's avatar
Theo Steininger committed
733

734
735
736
737
738
739
740
741
742
743
744
        Parameters
        ----------
        new_val : number, numpy.array, distributed_data_object, Field
            If an array is provided it needs to have the same shape as the
            domain of the Field.

        See Also
        --------
        get_val

        """
745
        self.set_val(new_val=new_val, copy=False)
csongor's avatar
csongor committed
746

747
748
    @property
    def shape(self):
749
        """ Returns the shape of the Field/ it's domain.
Theo Steininger's avatar
Theo Steininger committed
750

751
        All axes lengths written down seperately in a tuple.
Theo Steininger's avatar
Theo Steininger committed
752

753
754
755
756
757
758
759
760
761
762
763
        Returns
        -------
        out : tuple
            The output object. The tuple contains the dimansions of the spaces
            in domain.

        See Also
        --------
        dim

        """
764
        shape_tuple = tuple(sp.shape for sp in self.domain)
765
766
767
768
        try:
            global_shape = reduce(lambda x, y: x + y, shape_tuple)
        except TypeError:
            global_shape = ()
csongor's avatar
csongor committed
769

770
        return global_shape
csongor's avatar
csongor committed
771

772
773
    @property
    def dim(self):
774
        """ Returns the dimension of the Field/it's domain.
Theo Steininger's avatar
Theo Steininger committed
775

776
        Multiplies all values from shape.
Theo Steininger's avatar
Theo Steininger committed
777

778
779
780
781
782
783
784
785
786
787
        Returns
        -------
        out : int
            The dimension of the Field.

        See Also
        --------
        shape

        """
788
        dim_tuple = tuple(sp.dim for sp in self.domain)
Theo Steininger's avatar
Theo Steininger committed
789
790
791
792
        try:
            return reduce(lambda x, y: x * y, dim_tuple)
        except TypeError:
            return 0
csongor's avatar
csongor committed
793

794
795
    @property
    def dof(self):
Theo Steininger's avatar
Theo Steininger committed
796
797
798
799
800
801
802
803
        dof = self.dim
        if issubclass(self.dtype.type, np.complexfloating):
            dof *= 2
        return dof

    @property
    def total_volume(self):
        volume_tuple = tuple(sp.total_volume for sp in self.domain)
804
        try:
Theo Steininger's avatar
Theo Steininger committed
805
            return reduce(lambda x, y: x * y, volume_tuple)
806
        except TypeError:
Theo Steininger's avatar
Theo Steininger committed
807
            return 0.
808

Theo Steininger's avatar
Theo Steininger committed
809
    # ---Special unary/binary operations---
810

csongor's avatar
csongor committed
811
    def cast(self, x=None, dtype=None):
812
        """ Transforms x to a d2o with the same shape as the domain of 'self'
Theo Steininger's avatar
Theo Steininger committed
813

814
815
816
817
818
        Parameters
        ----------
        x : number, d2o, Field, array_like
            The input that shall be casted on a d2o of the same shape like the
            domain.
Theo Steininger's avatar
Theo Steininger committed
819

820
821
        dtype : type
            The datatype the output shall have.
Theo Steininger's avatar
Theo Steininger committed
822

823
824
825
826
827
828
829
830
831
832
        Returns
        -------
        out : distributed_data_object
            The output object.

        See Also
        --------
        _actual_cast

        """
csongor's avatar
csongor committed
833
834
        if dtype is None:
            dtype = self.dtype
835
836
        else:
            dtype = np.dtype(dtype)
837

838
839
        casted_x = x

840
        for ind, sp in enumerate(self.domain):
841
            casted_x = sp.pre_cast(casted_x,
842
843
844
                                   axes=self.domain_axes[ind])

        casted_x = self._actual_cast(casted_x, dtype=dtype)
845
846

        for ind, sp in enumerate(self.domain):
847
848
            casted_x = sp.post_cast(casted_x,
                                    axes=self.domain_axes[ind])
849

850
        return casted_x
csongor's avatar
csongor committed
851

Theo Steininger's avatar
Theo Steininger committed
852
    def _actual_cast(self, x, dtype=None):
853
        if isinstance(x, Field):
csongor's avatar
csongor committed
854
855
856
857
858
            x = x.get_val()

        if dtype is None:
            dtype = self.dtype

859
        return_x = distributed_data_object(
860
861
862
                            global_shape=self.shape,
                            dtype=dtype,
                            distribution_strategy=self.distribution_strategy)
863
864
        return_x.set_full_data(x, copy=False)
        return return_x
Theo Steininger's avatar
Theo Steininger committed
865

866
    def copy(self, domain=None, dtype=None, distribution_strategy=None):
867
        """ Returns a full copy of the Field.
Theo Steininger's avatar
Theo Steininger committed
868

869
870
871
872
873
874
875
876
877
        If no keyword arguments are given, the returned object will be an
        identical copy of the original Field. By explicit specification one is
        able to define the domain, the dtype and the distribution_strategy of
        the returned Field.

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

879
880
        dtype : type
            The new dtype the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
881

882
        distribution_strategy : all supported distribution strategies
Theo Steininger's avatar
Theo Steininger committed
883
884
            The new distribution strategy the Field shall have.

885
886
887
888
889
890
891
892
893
894
        Returns
        -------
        out : Field
            The output object. An identical copy of 'self'.

        See Also
        --------
        copy_empty

        """
Theo Steininger's avatar
Theo Steininger committed
895
        copied_val = self.get_val(copy=True)
896
897
898
899
        new_field = self.copy_empty(
                                domain=domain,
                                dtype=dtype,
                                distribution_strategy=distribution_strategy)
Theo Steininger's avatar
Theo Steininger committed
900
901
        new_field.set_val(new_val=copied_val, copy=False)
        return new_field
csongor's avatar
csongor committed
902

903
    def copy_empty(self, domain=None, dtype=None, distribution_strategy=None):
904
905
906
907
        """ Returns an empty copy of the Field.

        If no keyword arguments are given, the returned object will be an
        identical copy of the original Field containing random data. By
Theo Steininger's avatar
Theo Steininger committed
908
909
910
        explicit specification one is able to define the domain, the dtype and
        the distribution_strategy of the returned Field.

911
912
913
914
        Parameters
        ----------
        domain : DomainObject
            The new domain the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
915

916
917
        dtype : type
            The new dtype the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
918

919
920
        distribution_strategy : all supported distribution strategies
            The distribution strategy the new Field should have.
Theo Steininger's avatar
Theo Steininger committed
921

922
923
924
925
926
927
928
929
930
931
932
        Returns
        -------
        out : Field
            The output object. Contains random data.

        See Also
        --------
        copy
        _fast_copy_empty

        """
Theo Steininger's avatar
Theo Steininger committed
933
934
        if domain is None:
            domain = self.domain
csongor's avatar
csongor committed
935
        else:
Theo Steininger's avatar
Theo Steininger committed
936
            domain = self._parse_domain(domain)
csongor's avatar
csongor committed
937

Theo Steininger's avatar
Theo Steininger committed
938
939
940
941
        if dtype is None:
            dtype = self.dtype
        else:
            dtype = np.dtype(dtype)
csongor's avatar
csongor committed
942

943
944
        if distribution_strategy is None:
            distribution_strategy = self.distribution_strategy
csongor's avatar
csongor committed
945

Theo Steininger's avatar
Theo Steininger committed
946
947
948
949
950
951
952
953
954
955
        fast_copyable = True
        try:
            for i in xrange(len(self.domain)):
                if self.domain[i] is not domain[i]:
                    fast_copyable = False
                    break
        except IndexError:
            fast_copyable = False

        if (fast_copyable and dtype == self.dtype and
956
                distribution_strategy == self.distribution_strategy):
Theo Steininger's avatar
Theo Steininger committed
957
958
959
960
            new_field = self._fast_copy_empty()
        else:
            new_field = Field(domain=domain,
                              dtype=dtype,
961
                              distribution_strategy=distribution_strategy)
Theo Steininger's avatar
Theo Steininger committed
962
        return new_field
csongor's avatar
csongor committed
963

Theo Steininger's avatar
Theo Steininger committed
964
965
966
967
968
969
970
    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
        for key, value in self.__dict__.items():
971
            if key != '_val':
Theo Steininger's avatar
Theo Steininger committed
972
973
974
975
976
977
                new_field.__dict__[key] = value
            else:
                new_field.__dict__[key] = self.val.copy_empty()
        return new_field

    def weight(self, power=1, inplace=False, spaces=None):
978
979
980
981
982
983
984
985
        """ Devides every entry in 'self' by the dim of 'self'

        Parameters
        ----------
        power : number
            Here one can set the power to which the dimension is taken before
            division. power=2 will make the method devide every entry in 'self'
            by the square of the dimension.
Theo Steininger's avatar
Theo Steininger committed
986

987
988
        inplace : boolean
            For True the values in 'self' will be changed to the weighted ones.
Theo Steininger's avatar
Theo Steininger committed
989

990
991
        spaces : int
            Determines on what subspace the operation takes place.
Theo Steininger's avatar
Theo Steininger committed
992

993
994
995
996
997
998
        Returns
        -------
        out : Field
            The output object.

        """
999
        if inplace:
csongor's avatar
csongor committed
1000
1001
1002
1003
            new_field = self
        else:
            new_field = self.copy_empty()

1004
        new_val = self.get_val(copy=False)
csongor's avatar
csongor committed
1005

1006
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
csongor's avatar
csongor committed
1007
        if spaces is None:
Theo Steininger's avatar
Theo Steininger committed
1008
            spaces = range(len(self.domain))
csongor's avatar
csongor committed
1009

1010
        for ind, sp in enumerate(self.domain):
Theo Steininger's avatar
Theo Steininger committed
1011
1012
1013
1014
1015
            if ind in spaces:
                new_val = sp.weight(new_val,
                                    power=power,
                                    axes=self.domain_axes[ind],
                                    inplace=inplace)
1016
1017

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

1020
    def dot(self, x=None, spaces=None, bare=False):
1021
        """ Computes the dot product of 'self' with x.
Theo Steininger's avatar
Theo Steininger committed
1022

1023
        For a 1D Field this is the scalar product.
Theo Steininger's avatar
Theo Steininger committed
1024

1025
1026
1027
1028
        Parameters
        ----------
        x : Field
            Must have the same shape as 'self'
Theo Steininger's avatar
Theo Steininger committed
1029

1030
        spaces : int
Theo Steininger's avatar
Theo Steininger committed
1031
1032


1033
1034
1035
1036
1037
        bare : boolean
            bare=True operation will compute the sum over the pointwise product
            of 'self' and x.
            With bare=False this number will be devided by the dimension of the
            space over which the dotproduct is comupted.
Theo Steininger's avatar
Theo Steininger committed
1038

1039
1040
1041
        Returns
        -------
        out : float, complex
Theo Steininger's avatar
Theo Steininger committed
1042

1043
        """
1044
1045
1046
        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
1047

Martin Reinecke's avatar
Martin Reinecke committed
1048
        # Compute the dot respecting the fact of discrete/continuous spaces
Theo Steininger's avatar
Theo Steininger committed
1049
1050
1051
1052
1053
        if bare:
            y = self
        else:
            y = self.weight(power=1)

1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
        if spaces is None:
            x_val = x.get_val(copy=False)
            y_val = y.get_val(copy=False)
            result = (x_val.conjugate() * y_val).sum()
            return result
        else:
            # create a diagonal operator which is capable of taking care of the
            # axes-matching
            from nifty.operators.diagonal_operator import DiagonalOperator
            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
1069

1070
    def norm(self, q=2):
1071
        """ Computes the Lq-norm of the field values.
csongor's avatar
csongor committed
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083

            Parameters
            ----------
            q : scalar
                Parameter q of the Lq-norm (default: 2).

            Returns
            -------
            norm : scalar
                The Lq-norm of the field values.

        """
1084
        if q == 2:
1085
            return (self.dot(x=self)) ** (1 / 2)
csongor's avatar
csongor committed
1086
        else:
1087
            return self.dot(x=self ** (q - 1)) ** (1 / q)
csongor's avatar
csongor committed
1088
1089

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

1092
1093
1094
1095
        Parameters
        ----------
        inplace : boolean
            Decides whether self or a copied version of self shall be used
Theo Steininger's avatar
Theo Steininger committed
1096

1097
1098
1099
1100
        Returns
        -------
        cc : field
            The complex conjugated field.
csongor's avatar
csongor committed
1101
1102
1103
1104
1105
1106
1107

        """
        if inplace:
            work_field = self
        else:
            work_field = self.copy_empty()

1108
        new_val = self.get_val(copy=False)
Theo Steininger's avatar
Theo Steininger committed
1109
        new_val = new_val.conjugate()
1110
        work_field.set_val(new_val=new_val, copy=False)
csongor's avatar
csongor committed
1111
1112
1113

        return work_field

Theo Steininger's avatar
Theo Steininger committed
1114
    # ---General unary/contraction methods---
1115

Theo Steininger's avatar
Theo Steininger committed
1116
    def __pos__(self):
1117
1118
1119
1120
        """ x.__pos__() <==> +x

        Returns a (positive) copy of `self`.
        """
Theo Steininger's avatar
Theo Steininger committed
1121
        return self.copy()
1122

Theo Steininger's avatar
Theo Steininger committed
1123
    def __neg__(self):
1124
1125
1126
1127
        """ x.__neg__() <==> -x

        Returns a negative copy of `self`.
        """
Theo Steininger's avatar
Theo Steininger committed
1128
1129
1130
        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
1131
1132
        return return_field

Theo Steininger's avatar
Theo Steininger committed
1133
    def __abs__(self):
1134
1135
1136
1137
        """ x.__abs__() <==> abs(x)

        Returns an absolute valued copy of `self`.
        """
Theo Steininger's avatar
Theo Steininger committed
1138
1139
1140
1141
        return_field = self.copy_empty()
        new_val = abs(self.get_val(copy=False))
        return_field.set_val(new_val, copy=False)
        return return_field
csongor's avatar
csongor committed
1142

1143
    def _contraction_helper(self, op, spaces):
Theo Steininger's avatar
Theo Steininger committed
1144
1145
1146
1147
1148
        # build a list of all axes
        if spaces is None:
            spaces = xrange(len(self.domain))
        else:
            spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
csongor's avatar
csongor committed
1149

1150
        axes_list = tuple(self.domain_axes[sp_index] for sp_index in spaces)
1151
1152

        try:
Theo Steininger's avatar
Theo Steininger committed
1153
            axes_list = reduce(lambda x, y: x+y, axes_list)