field.py 51.2 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):
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
    """ Combines D2O-objects with meta-information NIFTY needs for computations.
        
    In NIFTY, Fields are used to store dataarrays and carry all the needed
    metainformation (i.e. the domain) for operators to be able to work on them.
    In addition Field has functions to analyse the power spectrum of it's
    content or create a random field of it.
    
    Parameters
    ----------
    domain : DomainObject
        One of the space types NIFTY supports. RGSpace, GLSpace, HPSpace,
        LMSpace or PowerSpace. It might also be a FieldArray, which is 
        an unstructured domain.
    
    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.
    
    dtype : type
        A numpy.type. Most common are int, float and complex.
    
    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
        
    copy: boolean

    Attributes
    ----------
    val : distributed_data_object
            
    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
        Name of the used distribution_strategy.  
    
    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
    
    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')
    
    See Also
    --------
    distributed_data_object

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

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

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

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

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

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

128
    def _parse_domain(self, domain, val=None):
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
        """ 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.
        """
144
        if domain is None:
145
146
147
148
            if isinstance(val, Field):
                domain = val.domain
            else:
                domain = ()
149
        elif isinstance(domain, DomainObject):
150
            domain = (domain,)
151
152
153
        elif not isinstance(domain, tuple):
            domain = tuple(domain)

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

Theo Steininger's avatar
Theo Steininger committed
161
    def _get_axes_tuple(self, things_with_shape, start=0):
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
        """ 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
        
        Returns
        -------
        out : tuple
            Incremental numeration of all axes.
        
        Note
        ----
        
        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
185
186
187
188
189
190
191
192
193
        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)
194

195
    def _infer_dtype(self, dtype, val):
196
197
198
199
200
201
202
203
204
205
206
207
208
209
        """ 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.
        
        Returns
        -------
        out : np.dtype
        """
csongor's avatar
csongor committed
210
        if dtype is None:
211
            try:
212
                dtype = val.dtype
213
            except AttributeError:
Theo Steininger's avatar
Theo Steininger committed
214
215
216
                try:
                    if val is None:
                        raise TypeError
217
                    dtype = np.result_type(val)
Theo Steininger's avatar
Theo Steininger committed
218
                except(TypeError):
219
                    dtype = np.dtype(gc['default_field_dtype'])
Theo Steininger's avatar
Theo Steininger committed
220
        else:
221
            dtype = np.dtype(dtype)
222

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

225
    def _parse_distribution_strategy(self, distribution_strategy, val):
226
227
228
229
230
231
232
233
234
235
236
237
238
        """ Sets a distribution strategy which the underlying D2O then uses.

        Parameters
        ----------
        distribution_strategy : all supported distribution strategies
            The distribution strategy the new d2o should have.
        val : Field or D2O
            Can infere the distribution strategy from another given object.
        
        Returns
        -------
        out : distribution_strategy
        """
239
        if distribution_strategy is None:
240
            if isinstance(val, distributed_data_object):
241
                distribution_strategy = val.distribution_strategy
242
            elif isinstance(val, Field):
243
                distribution_strategy = val.distribution_strategy
244
            else:
245
                self.logger.debug("distribution_strategy set to default!")
246
                distribution_strategy = gc['default_distribution_strategy']
247
        elif distribution_strategy not in DISTRIBUTION_STRATEGIES['global']:
248
249
250
            raise ValueError(
                    "distribution_strategy must be a global-type "
                    "strategy.")
251
        return distribution_strategy
252
253

    # ---Factory methods---
254

255
    @classmethod
256
    def from_random(cls, random_type, domain=None, dtype=None,
257
                    distribution_strategy=None, **kwargs):
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
        """ Draws a random field with the given parameters.

        Parameters
        ----------
        cls : class
        
        random_type : String
            'pm1', 'normal', 'uniform' are the supported arguments for this
            method.
        
        domain : DomainObject
            The domain of the output random field
        
        dtype : type
            The datatype of the output random field
        
        distribution_strategy : all supported distribution strategies
            The distribution strategy of the output random field
        
        Returns
        -------
        out : Field
            The output object.

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

        """
287
        # create a initially empty field
288
        f = cls(domain=domain, dtype=dtype,
289
                distribution_strategy=distribution_strategy)
290
291
292
293
294
295
296

        # 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
297
        # extract the distributed_data_object from f and apply the appropriate
298
299
300
301
302
303
304
305
306
307
308
        # 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):
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
        """ Processes the arguments given to a unified nomenclature.

        Parameters
        ----------
        random_type : String
            'pm1', 'normal', 'uniform'
        
        f : not used! FIXME: delete or use!
        
        Raise
        -----
        KeyError
            Raised if
                *the given randomtype is not 'pm1', 'normal' or 'uniform'
        
        Returns
        -------
        out : The generated random arguments.

        See Also
        --------
        from_random
331

332
        """
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
        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
348
        else:
349
350
            raise KeyError(
                "unsupported random key '" + str(random_type) + "'.")
csongor's avatar
csongor committed
351

352
        return random_arguments
csongor's avatar
csongor committed
353

354
355
    # ---Powerspectral methods---

356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
    def power_analyze(self, spaces=None, log=False, nbin=None, binbounds=None,
                      real_signal=True):
        """ Computes the powerspectrum of the Field
        
        Creates a PowerSpace with the given attributes and computes the 
        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
                
        log : boolean, *optional*
            True if the output PowerSpace should have log binning.
            {default : False}
        
        nbin : int, None, *optional*
            The number of bins the resulting PowerSpace shall have.
            {default : None}
                if nbin==None : maximum number of bins is used
            
        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
        
        Returns
        -------
        out : Field 
            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
408
        # check if all spaces in `self.domain` are either harmonic or
409
410
411
        # power_space instances
        for sp in self.domain:
            if not sp.harmonic and not isinstance(sp, PowerSpace):
Theo Steininger's avatar
Theo Steininger committed
412
                self.logger.info(
413
                    "Field has a space in `domain` which is neither "
414
415
416
                    "harmonic nor a PowerSpace.")

        # check if the `spaces` input is valid
417
418
419
420
421
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
        if spaces is None:
            if len(self.domain) == 1:
                spaces = (0,)
            else:
422
423
424
                raise ValueError(
                    "Field has multiple spaces as domain "
                    "but `spaces` is None.")
425
426

        if len(spaces) == 0:
427
428
            raise ValueError(
                "No space for analysis specified.")
429
        elif len(spaces) > 1:
430
431
            raise ValueError(
                "Conversion of only one space at a time is allowed.")
432
433
434
435

        space_index = spaces[0]

        if not self.domain[space_index].harmonic:
436
437
            raise ValueError(
                "The analyzed space must be harmonic.")
438

439
440
441
442
443
444
        # 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.

445
446
447
448
        distribution_strategy = \
            self.val.get_axes_local_distribution_strategy(
                self.domain_axes[space_index])

449
450
        harmonic_domain = self.domain[space_index]
        power_domain = PowerSpace(harmonic_domain=harmonic_domain,
451
                                  distribution_strategy=distribution_strategy,
452
                                  log=log, nbin=nbin, binbounds=binbounds)
453

454
        # extract pindex and rho from power_domain
455
456
        pindex = power_domain.pindex
        rho = power_domain.rho
457

458
        if real_signal:
459
            hermitian_part, anti_hermitian_part = \
460
                harmonic_domain.hermitian_decomposition(
461
462
463
464
465
466
467
468
469
470
471
472
473
474
                                            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(
475
476
477
478
479
480
481
482
483
                                            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

484
        if real_signal:
485
486
487
488
            result_dtype = np.complex
        else:
            result_dtype = np.float

489
490
        result_field = self.copy_empty(
                   domain=result_domain,
491
                   dtype=result_dtype,
492
                   distribution_strategy=power_spectrum.distribution_strategy)
493
494
495
496
497
        result_field.set_val(new_val=power_spectrum, copy=False)

        return result_field

    def _calculate_power_spectrum(self, x, pindex, rho, axes=None):
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
        """ Calculates the powerspectrum of the given Field x.

        Parameters
        ----------
        x : Field
            The Field of which the power spectrum shall be calculated
        
        pindex : distributed_data_object
        
        rho : numpy.array
        
        axes : tuple
            Used if only a subspace of the whole field is analysed.
            Coresponds to the 'spaces' keyword in power_analyse.
        
        Returns
        -------
        out : distributed_data_object
            The output object. Contains the power spectrum.

        See Also
        --------
        power_analyze

        """
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
        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):
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
572
573
574
575
576
577
        """ Makes the pindex array have the right size.
        
        If the 'spaces' keyword is used in power_analyze this method also
        shapes the pindex array the right way.
        
        Parameters
        ----------
        pindex : distributed_data_object
            
        
        target_shape : tuple
            
        
        target_strategy : a global distribution_strategy
            
        
        axes : tuple
            
        Raise
        -----
        ValueError
            Raised if
                *A wrong distribution strategy of the pindex is provided.
        
        Returns
        -------
        out : distributed_data_object
            The output object. 

        See Also
        --------
        _calculate_power_spectrum

        """
578
579
        if pindex.distribution_strategy not in \
                DISTRIBUTION_STRATEGIES['global']:
580
            raise ValueError("pindex's distribution strategy must be "
581
582
583
584
585
586
                             "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(
587
                    "A slicing distributor shall not be reshaped to "
588
589
590
591
592
593
594
595
596
597
598
599
600
                    "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

601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
    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.
        
        This method draws a Gaussian random field in the harmic partner domain.
        The drawn field has this field as its power spectrum.
        
        Notes
        -----
        For this the domain must be a PowerSpace.
        
        Parameters
        ----------
        spaces : {tuple, int, None} *optional*
            Specifies the subspace in which the power will be synthesized in
            case of a product space. 
            {default : None}
                if spaces==None : Tries to synthesize for the whole domain
        
        real_power : boolean *optional*
            Determines whether the power spectrum is real or complex
            {default : True}
            
        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*
            The standard deviation of the noise field the powerspectrum will be 
            multiplied on.
            {default : None}
                if std==None : std will be set to 1
        
        Returns
        -------
        out : Field
            The output object. A random field created with the power spectrum
            stored in 'self'

        See Also
        --------
        power_analyze
648

649
        """
650
651
652
        # check if the `spaces` input is valid
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))

Theo Steininger's avatar
Theo Steininger committed
653
654
655
        if spaces is None:
            spaces = range(len(self.domain))

656
657
658
659
660
        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.")
661
662
663

        # create the result domain
        result_domain = list(self.domain)
664
665
        for power_space_index in spaces:
            power_space = self.domain[power_space_index]
666
667
            harmonic_domain = power_space.harmonic_domain
            result_domain[power_space_index] = harmonic_domain
668
669
670

        # create random samples: one or two, depending on whether the
        # power spectrum is real or complex
671
        if real_power:
672
            result_list = [None]
673
674
        else:
            result_list = [None, None]
675

676
677
        result_list = [self.__class__.from_random(
                             'normal',
678
679
680
                             mean=mean,
                             std=std,
                             domain=result_domain,
681
                             dtype=np.complex,
682
                             distribution_strategy=self.distribution_strategy)
683
684
685
686
687
688
                       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
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706

        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)

707
        if real_signal:
708
            for power_space_index in spaces:
709
710
                harmonic_domain = result_domain[power_space_index]
                result_val_list = [harmonic_domain.hermitian_decomposition(
711
712
713
714
715
                                    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)]
716
717
718
719
720
721
722

        # 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]
723
        else:
724
725
726
727
728
            result = result_list[0] + 1j*result_list[1]

        return result

    def _spec_to_rescaler(self, spec, result_list, power_space_index):
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
        """ Transforms a 1D power spectrum to a 2D Field.
        
        This happens in order to make it applicable to a random Field by mere 
        pointwise multiplication.

        Parameters
        ----------
        spec : numpy.array
            Has the power spectrum stored in it.
        
        result_list : Field
            Used to infere the local_distribution_strategy
        
        power_space_index :
            Basically the 'spaces' keyword
        
        Returns
        -------
        out : Field
            The output object. The power spectrum in 2D.

        See Also
        --------
        power_synthesize

        """
755
        power_space = self.domain[power_space_index]
756
757
758

        # weight the random fields with the power spectrum
        # therefore get the pindex from the power space
759
        pindex = power_space.pindex
760
761
762
763
764
765
766
        # 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:
767
            self.logger.warn(
768
                "The distribution_stragey of pindex does not fit the "
769
770
771
772
773
774
775
776
777
778
                "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
779
780
        local_rescaler = spec[local_blow_up]
        return local_rescaler
781

Theo Steininger's avatar
Theo Steininger committed
782
    # ---Properties---
783

Theo Steininger's avatar
Theo Steininger committed
784
    def set_val(self, new_val=None, copy=False):
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
        """ Let's one set the values of the.

        Parameters
        ----------
        new_val : number, numpy.array, distributed_data_object, 
                Field, None, *optional*
            The values to be stored in the field.
            {default : None}
                if new_val==None : sets the values to 0.
        
        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

        """
804
805
        new_val = self.cast(new_val)
        if copy:
Theo Steininger's avatar
Theo Steininger committed
806
807
            new_val = new_val.copy()
        self._val = new_val
808
        return self
csongor's avatar
csongor committed
809

810
    def get_val(self, copy=False):
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
        """ Acceses the values stored in the Field.

        Parameters
        ----------
        copy : boolean
            True makes the method retrun a COPY of the Field's underlying 
            distributed_data_object.
        
        Returns
        -------
        out : distributed_data_object
            The output object.

        See Also
        --------
        val

        """
829
830
831
        if self._val is None:
            self.set_val(None)

832
        if copy:
Theo Steininger's avatar
Theo Steininger committed
833
            return self._val.copy()
834
        else:
Theo Steininger's avatar
Theo Steininger committed
835
            return self._val
csongor's avatar
csongor committed
836

Theo Steininger's avatar
Theo Steininger committed
837
838
    @property
    def val(self):
839
840
841
842
843
844
845
846
847
848
849
850
        """ Retruns actual distributed_data_object associated with this Field.
        
        Returns
        -------
        out : distributed_data_object
            The output object.

        See Also
        --------
        get_val

        """
851
        return self.get_val(copy=False)
csongor's avatar
csongor committed
852

Theo Steininger's avatar
Theo Steininger committed
853
854
    @val.setter
    def val(self, new_val):
855
856
857
858
859
860
861
862
863
864
865
866
867
        """ Sets the values in the d2o of the Field.
        
        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

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

870
871
    @property
    def shape(self):
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
        """ Returns the shape of the Field/ it's domain.
        
        All axes lengths written down seperately in a tuple.
        
        Returns
        -------
        out : tuple
            The output object. The tuple contains the dimansions of the spaces
            in domain.

        See Also
        --------
        dim

        """
887
        shape_tuple = tuple(sp.shape for sp in self.domain)
888
889
890
891
        try:
            global_shape = reduce(lambda x, y: x + y, shape_tuple)
        except TypeError:
            global_shape = ()
csongor's avatar
csongor committed
892

893
        return global_shape
csongor's avatar
csongor committed
894

895
896
    @property
    def dim(self):
897
898
899
900
901
902
903
904
905
906
907
908
909
910
        """ Returns the dimension of the Field/it's domain.
        
        Multiplies all values from shape.
        
        Returns
        -------
        out : int
            The dimension of the Field.

        See Also
        --------
        shape

        """
911
        dim_tuple = tuple(sp.dim for sp in self.domain)
Theo Steininger's avatar
Theo Steininger committed
912
913
914
915
        try:
            return reduce(lambda x, y: x * y, dim_tuple)
        except TypeError:
            return 0
csongor's avatar
csongor committed
916

917
918
    @property
    def dof(self):
Theo Steininger's avatar
Theo Steininger committed
919
920
921
922
923
924
925
926
        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)
927
        try:
Theo Steininger's avatar
Theo Steininger committed
928
            return reduce(lambda x, y: x * y, volume_tuple)
929
        except TypeError:
Theo Steininger's avatar
Theo Steininger committed
930
            return 0.
931

Theo Steininger's avatar
Theo Steininger committed
932
    # ---Special unary/binary operations---
933

csongor's avatar
csongor committed
934
    def cast(self, x=None, dtype=None):
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
        """ Transforms x to a d2o with the same shape as the domain of 'self'
        
        Parameters
        ----------
        x : number, d2o, Field, array_like
            The input that shall be casted on a d2o of the same shape like the
            domain.
            
        dtype : type
            The datatype the output shall have.
        
        Returns
        -------
        out : distributed_data_object
            The output object.

        See Also
        --------
        _actual_cast

        """
csongor's avatar
csongor committed
956
957
        if dtype is None:
            dtype = self.dtype
958
959
        else:
            dtype = np.dtype(dtype)
960

961
962
        casted_x = x

963
        for ind, sp in enumerate(self.domain):
964
            casted_x = sp.pre_cast(casted_x,
965
966
967
                                   axes=self.domain_axes[ind])

        casted_x = self._actual_cast(casted_x, dtype=dtype)
968
969

        for ind, sp in enumerate(self.domain):
970
971
            casted_x = sp.post_cast(casted_x,
                                    axes=self.domain_axes[ind])
972

973
        return casted_x
csongor's avatar
csongor committed
974

Theo Steininger's avatar
Theo Steininger committed
975
    def _actual_cast(self, x, dtype=None):
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
        """ Takes an x and makes a d2o out of it.
        
        Here the actual d2o is created for the cast() method.

        Parameters
        ----------
        x : array_like, Field
            The shape of x and 'self' must allready match.
        dtype : type
            The datatpye the returned d2o shall have.
        
        Returns
        -------
        out : distributed_data_object
            The output object.

        See Also
        --------
        cast

        """
997
        if isinstance(x, Field):
csongor's avatar
csongor committed
998
999
1000
1001
1002
            x = x.get_val()

        if dtype is None:
            dtype = self.dtype

1003
        return_x = distributed_data_object(
1004
1005
1006
                            global_shape=self.shape,
                            dtype=dtype,
                            distribution_strategy=self.distribution_strategy)
1007
1008
        return_x.set_full_data(x, copy=False)
        return return_x
Theo Steininger's avatar
Theo Steininger committed
1009

1010
    def copy(self, domain=None, dtype=None, distribution_strategy=None):
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
        """ Retruns a full copy of the Field.
        
        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.
        
        dtype : type
            The new dtype the Field shall have.
            
        distribution_strategy : all supported distribution strategies
            The new distribution strategy the Field shall have.        
        
        Returns
        -------
        out : Field
            The output object. An identical copy of 'self'.

        See Also
        --------
        copy_empty

        """
Theo Steininger's avatar
Theo Steininger committed
1039
        copied_val = self.get_val(copy=True)
1040
1041
1042
1043
        new_field = self.copy_empty(
                                domain=domain,
                                dtype=dtype,
                                distribution_strategy=distribution_strategy)
Theo Steininger's avatar
Theo Steininger committed
1044
1045
        new_field.set_val(new_val=copied_val, copy=False)
        return new_field
csongor's avatar
csongor committed
1046

1047
    def copy_empty(self, domain=None, dtype=None, distribution_strategy=None):
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
        """ 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
        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.
            
        dtype : type
            The new dtype the Field shall have.
            
        distribution_strategy : all supported distribution strategies
            The distribution strategy the new Field should have.
        
        Returns
        -------
        out : Field
            The output object. Contains random data.

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

        """
Theo Steininger's avatar
Theo Steininger committed
1077
1078
        if domain is None:
            domain = self.domain
csongor's avatar
csongor committed
1079
        else:
Theo Steininger's avatar
Theo Steininger committed
1080
            domain = self._parse_domain(domain)
csongor's avatar
csongor committed
1081

Theo Steininger's avatar
Theo Steininger committed
1082
1083
1084
1085
        if dtype is None:
            dtype = self.dtype
        else:
            dtype = np.dtype(dtype)
csongor's avatar
csongor committed
1086

1087
1088
        if distribution_strategy is None:
            distribution_strategy = self.distribution_strategy
csongor's avatar
csongor committed
1089

Theo Steininger's avatar
Theo Steininger committed
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
        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
1100
                distribution_strategy == self.distribution_strategy):
Theo Steininger's avatar
Theo Steininger committed
1101
1102
1103
1104
            new_field = self._fast_copy_empty()
        else:
            new_field = Field(domain=domain,
                              dtype=dtype,
1105
                              distribution_strategy=distribution_strategy)
Theo Steininger's avatar
Theo Steininger committed
1106
        return new_field
csongor's avatar
csongor committed
1107

Theo Steininger's avatar
Theo Steininger committed
1108
    def _fast_copy_empty(self):
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
        """ Creates an empty Field with the same properties as 'self'.
        
        Returns
        -------
        out : Field
            The output object.

        See Also
        --------
        fast_copy

        """
Theo Steininger's avatar
Theo Steininger committed
1121
1122
1123
1124
1125
1126
        # 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():
1127
            if key != '_val':
Theo Steininger's avatar
Theo Steininger committed
1128
1129
1130
1131
1132
1133
                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):
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
        """ 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.
            
        inplace : boolean
            For True the values in 'self' will be changed to the weighted ones.
        
        spaces : int
            Determines on what subspace the operation takes place.
        
        Returns
        -------
        out : Field
            The output object.

        """
1155
        if inplace:
csongor's avatar
csongor committed
1156
1157
1158
1159
            new_field = self
        else:
            new_field = self.copy_empty()

1160
        new_val = self.get_val(copy=False)
csongor's avatar
csongor committed
1161

1162
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
csongor's avatar
csongor committed
1163
        if spaces is None:
Theo Steininger's avatar
Theo Steininger committed
1164
            spaces = range(len(self.domain))
csongor's avatar
csongor committed
1165

1166
        for ind, sp in enumerate(self.domain):
Theo Steininger's avatar
Theo Steininger committed
1167
1168
1169
1170
1171
            if ind in spaces:
                new_val = sp.weight(new_val,
                                    power=power,
                                    axes=self.domain_axes[ind],
                                    inplace=inplace)
1172
1173

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

1176
    def dot(self, x=None, spaces=None, bare=False):
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
        """ Computes the dot product of 'self' with x.
        
        For a 1D Field this is the scalar product.
        
        Parameters
        ----------
        x : Field
            Must have the same shape as 'self'
        
        spaces : int
            
        
        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.
        
        Returns
        -------
        out : float, complex
        
        """
1200
1201
1202
        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
1203

Martin Reinecke's avatar
Martin Reinecke committed
1204
        # Compute the dot respecting the fact of discrete/continuous spaces
Theo Steininger's avatar
Theo Steininger committed
1205
1206
1207
1208
1209
        if bare:
            y = self
        else:
            y = self.weight(power=1)

1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
        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
1225

1226
    def norm(self, q=2):
1227
        """ Computes the Lq-norm of the field values.
csongor's avatar
csongor committed
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239

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

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

        """
1240
        if q == 2:
1241
            return (self.dot(x=self)) ** (1 / 2)
csongor's avatar
csongor committed
1242
        else:
1243
            return self.dot(x=self ** (q - 1)) ** (1 / q)
csongor's avatar
csongor committed
1244
1245

    def conjugate(self, inplace=False):
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
        """ Retruns the complex conjugate of the field.
        
        Parameters
        ----------
        inplace : boolean
            Decides whether self or a copied version of self shall be used
            
        Returns
        -------
        cc : field
            The complex conjugated field.
csongor's avatar
csongor committed
1257
1258
1259
1260
1261
1262
1263

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

1264
        new_val = self.get_val(copy=False)
Theo Steininger's avatar
Theo Steininger committed
1265
        new_val = new_val.conjugate()
1266
        work_field.set_val(new_val=new_val, copy=False)
csongor's avatar
csongor committed
1267
1268
1269

        return work_field

Theo Steininger's avatar
Theo Steininger committed
1270
    # ---General unary/contraction methods---
1271

Theo Steininger's avatar
Theo Steininger committed
1272
    def __pos__(self):
1273
1274
1275
1276
        """ x.__pos__() <==> +x

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

Theo Steininger's avatar
Theo Steininger committed
1279
    def __neg__(self):
1280
1281
1282
1283
        """ x.__neg__() <==> -x

        Returns a negative copy of `self`.
        """
Theo Steininger's avatar
Theo Steininger committed
1284
1285
1286
        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
1287
1288
        return return_field

Theo Steininger's avatar
Theo Steininger committed
1289
    def __abs__(self):
1290
1291
1292
1293
        """ x.__abs__() <==> abs(x)

        Returns an absolute valued copy of `self`.
        """
Theo Steininger's avatar
Theo Steininger committed
1294
1295
1296
1297
        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
1298

1299
    def _contraction_helper(self, op, spaces):
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
        """ Used for contraction operations like sum, norms, mean, std...

        Parameters
        ----------
        op : callable

        spaces : 
            The subspace on which the contraction should be carried out.
        
        Returns
        -------
        out : float
            The result of the operation performed on the values stored in 'self'

        """
Theo Steininger's avatar
Theo Steininger committed
1315
1316
1317
1318
1319
        # 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
1320

1321
        axes_list = tuple(self.domain_axes[sp_index] for sp_index in spaces)
1322
1323

        try:
Theo Steininger's avatar
Theo Steininger committed
1324
            axes_list = reduce(lambda x, y: x+y, axes_list)
1325
        except TypeError:
Theo Steininger's avatar
Theo Steininger committed
1326
            axes_list = ()
csongor's avatar
csongor committed
1327

Theo Steininger's avatar
Theo Steininger committed
1328
1329
1330
        # perform the contraction on the d2o
        data = self.get_val(copy=False)
        data = getattr(data, op)(axis=axes_list)
csongor's avatar
csongor committed
1331

Theo Steininger's avatar
Theo Steininger committed
1332
1333
1334
        # check if the result is scalar or if a result_field must be constr.
        if np.isscalar(data):
            return data
csongor's avatar
csongor committed
1335
        else:
Theo Steininger's avatar
Theo Steininger committed
1336
1337
1338
            return_domain = tuple(self.domain[i]
                                  for i in xrange(len(self.domain))
                                  if i not in spaces)
1339

Theo Steininger's avatar
Theo Steininger committed
1340
1341
1342
1343
            return_field = Field(domain=return_domain,
                                 val=data,
                                 copy=False)
            return return_field
csongor's avatar
csongor committed
1344

1345
1346
    def sum(self, spaces=None):
        return self._contraction_helper('sum', spaces)
csongor's avatar
csongor committed
1347

1348
1349
    def prod(self, spaces=None):
        return self._contraction_helper('prod', spaces)
csongor's avatar
csongor committed
1350

1351
1352
    def all(self, spaces=None):
        return self._contraction_helper('all', spaces)
csongor's avatar
csongor committed
1353

1354
1355
    def any(self, spaces=None):
        return self._contraction_helper('any', spaces)
csongor's avatar
csongor committed
1356

1357
1358
    def min(self, spaces=None):
        return self._contraction_helper('min', spaces)
csongor's avatar
csongor committed
1359

1360
1361
    def nanmin(self, spaces=None):
        return self._contraction_helper('nanmin', spaces)
csongor's avatar
csongor committed
1362

1363
1364
    def max(self, spaces=None):
        return self._contraction_helper('max', spaces)
csongor's avatar
csongor committed
1365

1366
1367
    def nanmax(self, spaces=None):
        return self._contraction_helper('nanmax', spaces)
csongor's avatar
csongor committed
1368

1369
1370
    def mean(self, spaces=None):
        return self._contraction_helper('mean', spaces)
csongor's avatar
csongor committed
1371

1372
1373
    def var(self, spaces=None):
        return self._contraction_helper('var', spaces)
csongor's avatar
csongor committed
1374

1375
1376
    def std(self, spaces=None):
        return self._contraction_helper('std', spaces)
csongor's avatar
csongor committed
1377

Theo Steininger's avatar
Theo Steininger committed
1378
    # ---General binary methods---
csongor's avatar
csongor committed
1379

Theo Steininger's avatar
Theo Steininger committed
1380
    def _binary_helper(self, other, op, inplace=False):
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
        """ Used for various binary operations like +, -, *, /, **, *=, +=,...

        _binary_helper checks whether `other` has the same domain as 'self'.

        Parameters
        ----------
        other : scalar, array-like
        
        op : callable
        
        inplace : boolean
            If the result shall be saved in the data array of `self`. Used for
            +=, -=, etc...
        Returns
        -------
        out : Field
            The Field containing the computation's result.
            Equals `self` if `inplace is True`.

        """
csongor's avatar
csongor committed
1401
        # if other is a field, make sure that the domains match
1402
        if isinstance(other, Field):
Theo Steininger's avatar
Theo Steininger committed
1403
1404
1405
1406
1407
            try:
                assert len(other.domain) == len(self.domain)
                for index in xrange(len(self.domain)):
                    assert other.domain[index] == self.domain[index]
            except AssertionError:
1408
1409
                raise ValueError(
                    "domains are incompatible.")
Theo Steininger's avatar
Theo Steininger committed
1410
            other = other.get_val(copy=False)
csongor's avatar
csongor committed
1411

Theo Steininger's avatar
Theo Steininger committed
1412
1413
        self_val = self.get_val(copy=False)
        return_val = getattr(self_val, op)(other)
csongor's avatar
csongor committed
1414
1415
1416
1417

        if inplace:
            working_field = self
        else:
1418
            working_field = self.copy_empty(dtype=return_val.dtype)
csongor's avatar
csongor committed
1419

Theo Steininger's avatar
Theo Steininger committed
1420
        working_field.set_val(return_val, copy=False)
csongor's avatar
csongor committed
1421
1422
1423
        return working_field

    def __add__(self, other):
1424
1425
1426
1427
1428
1429
        """ x.__add__(y) <==> x+y

        See Also
        --------
        _binary_helper
        """
Theo Steininger's avatar
Theo Steininger committed
1430
        return self._binary_helper(other, op='__add__')
1431

1432
    def __radd__(self, other):
1433
1434
1435
1436