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

    # ---Factory methods---
241

242
    @classmethod
243
    def from_random(cls, random_type, domain=None, dtype=None,
244
                    distribution_strategy=None, **kwargs):
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
        """ 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

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

        # 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
284
        # extract the distributed_data_object from f and apply the appropriate
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
310
        # 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
311
        else:
312
313
            raise KeyError(
                "unsupported random key '" + str(random_type) + "'.")
csongor's avatar
csongor committed
314

315
        return random_arguments
csongor's avatar
csongor committed
316

317
318
    # ---Powerspectral methods---

319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
    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
371
        # check if all spaces in `self.domain` are either harmonic or
372
373
374
        # power_space instances
        for sp in self.domain:
            if not sp.harmonic and not isinstance(sp, PowerSpace):
Theo Steininger's avatar
Theo Steininger committed
375
                self.logger.info(
376
                    "Field has a space in `domain` which is neither "
377
378
379
                    "harmonic nor a PowerSpace.")

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

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

        space_index = spaces[0]

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

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

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

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

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

421
        if real_signal:
422
            hermitian_part, anti_hermitian_part = \
423
                harmonic_domain.hermitian_decomposition(
424
425
426
427
428
429
430
431
432
433
434
435
436
437
                                            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(
438
439
440
441
442
443
444
445
446
                                            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

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

452
453
        result_field = self.copy_empty(
                   domain=result_domain,
454
                   dtype=result_dtype,
455
                   distribution_strategy=power_spectrum.distribution_strategy)
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
483
        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']:
484
            raise ValueError("pindex's distribution strategy must be "
485
486
487
488
489
490
                             "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(
491
                    "A slicing distributor shall not be reshaped to "
492
493
494
495
496
497
498
499
500
501
502
503
504
                    "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

505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
    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
552

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

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

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

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

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

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

        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)

611
        if real_signal:
612
            for power_space_index in spaces:
613
614
                harmonic_domain = result_domain[power_space_index]
                result_val_list = [harmonic_domain.hermitian_decomposition(
615
616
617
618
619
                                    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)]
620
621
622
623
624
625
626

        # 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]
627
        else:
628
629
630
631
632
633
            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]
634
635
636

        # weight the random fields with the power spectrum
        # therefore get the pindex from the power space
637
        pindex = power_space.pindex
638
639
640
641
642
643
644
        # 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:
645
            self.logger.warn(
646
                "The distribution_stragey of pindex does not fit the "
647
648
649
650
651
652
653
654
655
656
                "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
657
658
        local_rescaler = spec[local_blow_up]
        return local_rescaler
659

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

Theo Steininger's avatar
Theo Steininger committed
662
    def set_val(self, new_val=None, copy=False):
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
        """ 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

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

688
    def get_val(self, copy=False):
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
        """ 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

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

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

Theo Steininger's avatar
Theo Steininger committed
715
716
    @property
    def val(self):
717
718
719
720
721
722
723
724
725
726
727
728
        """ Retruns actual distributed_data_object associated with this Field.
        
        Returns
        -------
        out : distributed_data_object
            The output object.

        See Also
        --------
        get_val

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

Theo Steininger's avatar
Theo Steininger committed
731
732
    @val.setter
    def val(self, new_val):
733
734
735
736
737
738
739
740
741
742
743
744
745
        """ 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

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

748
749
    @property
    def shape(self):
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
        """ 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

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

771
        return global_shape
csongor's avatar
csongor committed
772

773
774
    @property
    def dim(self):
775
776
777
778
779
780
781
782
783
784
785
786
787
788
        """ 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

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

795
796
    @property
    def dof(self):
Theo Steininger's avatar
Theo Steininger committed
797
798
799
800
801
802
803
804
        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)
805
        try:
Theo Steininger's avatar
Theo Steininger committed
806
            return reduce(lambda x, y: x * y, volume_tuple)
807
        except TypeError:
Theo Steininger's avatar
Theo Steininger committed
808
            return 0.
809

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

csongor's avatar
csongor committed
812
    def cast(self, x=None, dtype=None):
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
        """ 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
834
835
        if dtype is None:
            dtype = self.dtype
836
837
        else:
            dtype = np.dtype(dtype)
838

839
840
        casted_x = x

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

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

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

851
        return casted_x
csongor's avatar
csongor committed
852

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

        if dtype is None:
            dtype = self.dtype

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

867
    def copy(self, domain=None, dtype=None, distribution_strategy=None):
868
        """ Returns a full copy of the Field.
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
        
        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
896
        copied_val = self.get_val(copy=True)
897
898
899
900
        new_field = self.copy_empty(
                                domain=domain,
                                dtype=dtype,
                                distribution_strategy=distribution_strategy)
Theo Steininger's avatar
Theo Steininger committed
901
902
        new_field.set_val(new_val=copied_val, copy=False)
        return new_field
csongor's avatar
csongor committed
903

904
    def copy_empty(self, domain=None, dtype=None, distribution_strategy=None):
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
        """ 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
934
935
        if domain is None:
            domain = self.domain
csongor's avatar
csongor committed
936
        else:
Theo Steininger's avatar
Theo Steininger committed
937
            domain = self._parse_domain(domain)
csongor's avatar
csongor committed
938

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

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

Theo Steininger's avatar
Theo Steininger committed
947
948
949
950
951
952
953
954
955
956
        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
957
                distribution_strategy == self.distribution_strategy):
Theo Steininger's avatar
Theo Steininger committed
958
959
960
961
            new_field = self._fast_copy_empty()
        else:
            new_field = Field(domain=domain,
                              dtype=dtype,
962
                              distribution_strategy=distribution_strategy)
Theo Steininger's avatar
Theo Steininger committed
963
        return new_field
csongor's avatar
csongor committed
964

Theo Steininger's avatar
Theo Steininger committed
965
966
967
968
969
970
971
    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():
972
            if key != '_val':
Theo Steininger's avatar
Theo Steininger committed
973
974
975
976
977
978
                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):
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
        """ 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.

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

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

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

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

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

1021
    def dot(self, x=None, spaces=None, bare=False):
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
        """ 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
        
        """
1045
1046
1047
        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
1048

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

1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
        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
1070

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

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

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

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

    def conjugate(self, inplace=False):
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
        """ 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
1102
1103
1104
1105
1106
1107
1108

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

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

        return work_field

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

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

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

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

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

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

        Returns an absolute valued copy of `self`.
        """
Theo Steininger's avatar
Theo Steininger committed
1139
1140
1141
1142
        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
1143

1144
    def _contraction_helper(self, op, spaces):
Theo Steininger's avatar
Theo Steininger committed
1145
1146
1147
1148
1149
        # 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
1150

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

        try:
Theo Steininger's avatar
Theo Steininger committed
1154
            axes_list = reduce(lambda x, y: x+y, axes_list)
1155
        except TypeError:
Theo Steininger's avatar
Theo Steininger committed
1156
            axes_list = ()
csongor's avatar
csongor committed
1157

Theo Steininger's avatar
Theo Steininger committed
1158
1159
1160
        # perform the contraction on the d2o
        data = self.get_val(copy=False)
        data = getattr(data, op)(axis=axes_list)
csongor's avatar
csongor committed
1161

Theo Steininger's avatar
Theo Steininger committed
1162
1163
1164
        # 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
1165
        else:
Theo Steininger's avatar
Theo Steininger committed
1166
1167
1168
            return_domain = tuple(self.domain[i]
                                  for i in xrange(len(self.domain))
                                  if i not in spaces)
1169

Theo Steininger's avatar
Theo Steininger committed
1170
1171
1172
1173
            return_field = Field(domain=return_domain,
                                 val=data,
                                 copy=False)
            return return_field
csongor's avatar
csongor committed
1174

1175
1176
    def sum(self, spaces=None):
        return self._contraction_helper('sum', spaces)
csongor's avatar
csongor committed
1177

1178
1179
    def prod(self, spaces=None):
        return self._contraction_helper('prod', spaces)
csongor's avatar
csongor committed
1180

1181
1182
    def all(self, spaces=None):
        return self._contraction_helper('all', spaces)
csongor's avatar
csongor committed
1183

1184
1185
    def any(self, spaces=None):
        return self._contraction_helper('any', spaces)
csongor's avatar
csongor committed
1186

1187
1188
    def min(self, spaces=None):
        return self._contraction_helper('min', spaces)
csongor's avatar
csongor committed
1189

1190
1191
    def nanmin(self, spaces=None):
        return self._contraction_helper('nanmin', spaces)
csongor's avatar
csongor committed
1192

1193
1194
    def max(self, spaces=None):
        return self._contraction_helper('max', spaces)
csongor's avatar
csongor committed
1195

1196
1197
    def nanmax(self, spaces=None):
        return self._contraction_helper('nanmax', spaces)
csongor's avatar
csongor committed
1198

1199
1200
    def mean(self, spaces=None):
        return self._contraction_helper('mean', spaces)
csongor's avatar
csongor committed
1201

1202
1203
    def var(self, spaces=None):
        return self._contraction_helper('var', spaces)
csongor's avatar
csongor committed
1204

1205
1206
    def std(self, spaces=None):
        return self._contraction_helper('std', spaces)
csongor's avatar
csongor committed
1207

Theo Steininger's avatar
Theo Steininger committed
1208
    # ---General binary methods---
csongor's avatar
csongor committed
1209

Theo Steininger's avatar
Theo Steininger committed
1210
    def _binary_helper(self, other, op, inplace=False):
csongor's avatar
csongor committed
1211
        # if other is a field, make sure that the domains match
1212
        if isinstance(other, Field):
Theo Steininger's avatar
Theo Steininger committed
1213
1214
1215
1216
1217
            try:
                assert len(other.domain) == len(self.domain)
                for index in xrange(len(self.domain)):
                    assert other.domain[index] == self.domain[index]
            except AssertionError:
1218
1219
                raise ValueError(
                    "domains are incompatible.")
Theo Steininger's avatar
Theo Steininger committed
1220
            other = other.get_val(copy=False)
csongor's avatar
csongor committed
1221

Theo Steininger's avatar
Theo Steininger committed
1222
1223
        self_val = self.get_val(copy=False)
        return_val = getattr(self_val, op)(other)
csongor's avatar
csongor committed
1224
1225
1226
1227

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

Theo Steininger's avatar
Theo Steininger committed
1230
        working_field.set_val(return_val, copy=False)
csongor's avatar
csongor committed
1231
1232
1233
        return working_field

    def __add__(self, other):
1234
1235
1236
1237
1238
1239
        """ x.__add__(y) <==> x+y

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

1242
    def __radd__(self, other):
1243
1244
1245
1246
1247
1248
        """ x.__radd__(y) <==> y+x

        See Also
        --------
        _builtin_helper
        """
Theo Steininger's avatar
Theo Steininger committed
1249
        return self._binary_helper(other, op='__radd__')
csongor's avatar
csongor committed
1250
1251

    def __iadd__(self, other):
1252
1253
1254
1255
1256
1257
        """ x.__iadd__(y) <==> x+=y

        See Also
        --------
        _builtin_helper
        """
Theo Steininger's avatar
Theo Steininger committed
1258
        return self._binary_helper(other, op='__iadd__', inplace=True)
csongor's avatar
csongor committed
1259
1260

    def __sub__(self, other):
1261
1262
1263
1264
1265
1266
        """ x.__sub__(y) <==> x-y

        See Also
        --------
        _builtin_helper
        """
Theo Steininger's avatar
Theo Steininger committed
1267
        return self._binary_helper(other, op='__sub__')
csongor's avatar
csongor committed
1268
1269

    def __rsub__(self, other):
1270
1271
1272
1273
1274
1275
        """ x.__rsub__(y) <==> y-x

        See Also
        --------
        _builtin_helper
        """
Theo Steininger's avatar
Theo Steininger committed
1276
        return self._binary_helper(other, op='__rsub__')
csongor's avatar
csongor committed
1277
1278

    def __isub__(self, other):
1279
1280
1281
1282
1283
1284
        """ x.__isub__(y) <==> x-=y

        See Also
        --------
        _builtin_helper
        """
Theo Steininger's avatar
Theo Steininger committed
1285
        return self._binary_helper(other, op='__isub__', inplace=True)
csongor's avatar
csongor committed
1286
1287

    def __mul__(self, other):
1288
1289
1290
1291
1292
1293
        """ x.__mul__(y) <==> x*y

        See Also
        --------
        _builtin_helper
        """
Theo Steininger's avatar
Theo Steininger committed
1294
        return self._binary_helper(other, op='__mul__')
1295

1296
    def __rmul__(self, other):
1297
1298
1299
1300
1301
1302
        """ x.__rmul__(y) <==> y*x

        See Also
        --------
        _builtin_helper
        """
Theo Steininger's avatar
Theo Steininger committed
1303
        return self._binary_helper(other, op='__rmul__')
csongor's avatar
csongor committed
1304
1305

    def __imul__(self, other):
1306
1307
1308
1309
1310
1311
        """ x.__imul__(y) <==> x*=y

        See Also
        --------
        _builtin_helper
        """
Theo Steininger's avatar
Theo Steininger committed
1312
        return self._binary_helper(other, op='__imul__', inplace=True)
csongor's avatar
csongor committed
1313
1314

    def __div__(self, other):
1315
1316
1317
1318
1319
1320
        """ x.__div__(y) <==> x/y

        See Also
        --------
        _builtin_helper
        """
Theo Steininger's avatar
Theo Steininger committed
1321
        return self._binary_helper(other, op='__div__')
csongor's avatar
csongor committed
1322
1323

    def __rdiv__(self, other):
1324
1325
1326
1327
1328
1329
        """ x.__rdiv__(y) <==> y/x

        See Also
        --------
        _builtin_helper
        """
Theo Steininger's avatar
Theo Steininger committed
1330
        return self._binary_helper(other, op='__rdiv__')
csongor's avatar
csongor committed
1331
1332

    def __idiv__(self, other):
1333
1334
1335
1336
1337
1338
        """ x.__idiv__(y) <==> x/=y

        See Also
        --------
        _builtin_helper
        """
Theo Steininger's avatar
Theo Steininger committed