field.py 36 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
Theo Steininger's avatar
Theo Steininger committed
13
14
15
16
17
#
# Copyright(C) 2013-2017 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
18

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

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

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

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

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

csongor's avatar
csongor committed
34

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

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

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

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

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

57
58
59
60
    copy: boolean

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

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

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

78
    """
79

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

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

87
88
89
90
        if val is None:
            self._val = None
        else:
            self.set_val(new_val=val, copy=copy)
csongor's avatar
csongor committed
91

92
    def _parse_domain(self, domain, val=None):
93
        if domain is None:
94
95
96
97
            if isinstance(val, Field):
                domain = val.domain
            else:
                domain = ()
98
        elif isinstance(domain, DomainObject):
99
            domain = (domain,)
100
101
102
        elif not isinstance(domain, tuple):
            domain = tuple(domain)

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

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

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

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

137
    # ---Factory methods---
138

139
    @classmethod
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
140
    def from_random(cls, random_type, domain=None, dtype=None, **kwargs):
141
142
143
144
145
        """ Draws a random field with the given parameters.

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

147
148
149
        random_type : String
            'pm1', 'normal', 'uniform' are the supported arguments for this
            method.
Theo Steininger's avatar
Theo Steininger committed
150

151
152
        domain : DomainObject
            The domain of the output random field
Theo Steininger's avatar
Theo Steininger committed
153

154
155
        dtype : type
            The datatype of the output random field
Theo Steininger's avatar
Theo Steininger committed
156

157
158
159
160
161
162
163
        Returns
        -------
        out : Field
            The output object.

        See Also
        --------
164
        power_synthesize
Theo Steininger's avatar
Theo Steininger committed
165

166
167

        """
Theo Steininger's avatar
Theo Steininger committed
168

169
        # create a initially empty field
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
170
        f = cls(domain=domain, dtype=dtype)
171

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
172
        # extract the data from f and apply the appropriate
173
174
175
        # random number generator to it
        sample = f.get_val(copy=False)
        generator_function = getattr(Random, random_type)
176

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
177
178
        sample[:]=generator_function(dtype=f.dtype,
                                             shape=sample.shape,
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
179
                                             **kwargs)
180
181
        return f

182
183
    # ---Powerspectral methods---

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

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

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

213
214
215
216
        Raise
        -----
        ValueError
            Raised if
Theo Steininger's avatar
Theo Steininger committed
217
218
                *len(domain) is != 1 when spaces==None
                *len(spaces) is != 1 if not None
219
                *the analyzed space is not harmonic
Theo Steininger's avatar
Theo Steininger committed
220

221
222
        Returns
        -------
Theo Steininger's avatar
Theo Steininger committed
223
        out : Field
Martin Reinecke's avatar
Martin Reinecke committed
224
            The output object. Its domain is a PowerSpace and it contains
225
226
227
228
229
            the power spectrum of 'self's field.

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

231
        """
Theo Steininger's avatar
Theo Steininger committed
232

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

        # check if the `spaces` input is valid
242
243
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
244
            spaces = list(range(len(self.domain)))
245
246

        if len(spaces) == 0:
247
248
            raise ValueError(
                "No space for analysis specified.")
249

250
251
252
253
254
255
256
257
258
259
260
261
262
        if keep_phase_information:
            parts_val = self._hermitian_decomposition(
                                              domain=self.domain,
                                              val=self.val,
                                              spaces=spaces,
                                              domain_axes=self.domain_axes,
                                              preserve_gaussian_variance=False)
            parts = [self.copy_empty().set_val(part_val, copy=False)
                     for part_val in parts_val]
        else:
            parts = [self]

        parts = [abs(part)**2 for part in parts]
263
264

        for space_index in spaces:
265
266
            parts = [self._single_power_analyze(
                                work_field=part,
267
                                space_index=space_index,
268
269
                                binbounds=binbounds)
                     for part in parts]
270

271
272
273
274
275
276
        if keep_phase_information:
            result_field = parts[0] + 1j*parts[1]
        else:
            result_field = parts[0]

        return result_field
277
278

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

281
        if not work_field.domain[space_index].harmonic:
282
283
            raise ValueError(
                "The analyzed space must be harmonic.")
284

285
286
287
288
289
290
        # 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.

291
        harmonic_domain = work_field.domain[space_index]
292
        power_domain = PowerSpace(harmonic_partner=harmonic_domain,
Theo Steininger's avatar
Theo Steininger committed
293
                                  binbounds=binbounds)
294
295
        power_spectrum = cls._calculate_power_spectrum(
                                field_val=work_field.val,
Martin Reinecke's avatar
Martin Reinecke committed
296
                                pdomain=power_domain,
297
                                axes=work_field.domain_axes[space_index])
298
299

        # create the result field and put power_spectrum into it
300
        result_domain = list(work_field.domain)
301
        result_domain[space_index] = power_domain
302
        result_dtype = power_spectrum.dtype
303

304
        result_field = work_field.copy_empty(
305
                   domain=result_domain,
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
306
                   dtype=result_dtype)
307
308
309
310
        result_field.set_val(new_val=power_spectrum, copy=False)

        return result_field

311
    @classmethod
Martin Reinecke's avatar
Martin Reinecke committed
312
    def _calculate_power_spectrum(cls, field_val, pdomain, axes=None):
313

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

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

        return power_spectrum

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

        return result_obj

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

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

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

374
375
376
377
        Returns
        -------
        out : Field
            The output object. A random field created with the power spectrum
Theo Steininger's avatar
Theo Steininger committed
378
            stored in the `spaces` in `self`.
379

Theo Steininger's avatar
Theo Steininger committed
380
381
382
383
384
385
        Notes
        -----
        For this the spaces specified by `spaces` must be a PowerSpace.
        This expects this field to be the square root of a power spectrum, i.e.
        to have the unit of the field to be sampled.

386
387
388
        See Also
        --------
        power_analyze
Theo Steininger's avatar
Theo Steininger committed
389
390
391
392
393

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

394
        """
Theo Steininger's avatar
Theo Steininger committed
395

396
397
398
        # check if the `spaces` input is valid
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))

Theo Steininger's avatar
Theo Steininger committed
399
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
400
            spaces = list(range(len(self.domain)))
Theo Steininger's avatar
Theo Steininger committed
401

402
403
404
405
406
        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.")
407
408
409

        # create the result domain
        result_domain = list(self.domain)
410
411
        for power_space_index in spaces:
            power_space = self.domain[power_space_index]
412
            harmonic_domain = power_space.harmonic_partner
413
            result_domain[power_space_index] = harmonic_domain
414
415
416

        # create random samples: one or two, depending on whether the
        # power spectrum is real or complex
417
        if real_power:
418
            result_list = [None]
419
420
        else:
            result_list = [None, None]
421

422
423
        result_list = [self.__class__.from_random(
                             'normal',
424
425
426
                             mean=mean,
                             std=std,
                             domain=result_domain,
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
427
                             dtype=np.complex)
428
429
430
431
432
433
                       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
434

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
435
        spec = self.val.copy()
436
437
        spec = np.sqrt(spec)

438
439
440
441
442
443
444
        for power_space_index in spaces:
            spec = self._spec_to_rescaler(spec, result_list, power_space_index)
        local_rescaler = spec

        result_val_list = [x.val for x in result_list]

        # apply the rescaler to the random fields
Martin Reinecke's avatar
Martin Reinecke committed
445
        result_val_list[0] *= local_rescaler.real
446
447

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

450
        if real_signal:
451
            result_val_list = [self._hermitian_decomposition(
452
453
454
455
456
                                            result_domain,
                                            result_val,
                                            spaces,
                                            result_list[0].domain_axes,
                                            preserve_gaussian_variance=True)[0]
457
                               for result_val in result_val_list]
458
459
460
461
462
463
464

        # 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]
465
466
467
            if not issubclass(result_val_list[0].dtype.type,
                              np.complexfloating):
                result = result.real
468
        else:
469
470
471
472
            result = result_list[0] + 1j*result_list[1]

        return result

473
    @staticmethod
474
475
    def _hermitian_decomposition(domain, val, spaces, domain_axes,
                                 preserve_gaussian_variance=False):
476

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
477
478
        h = val.real.copy()
        a = 1j * val.imag.copy()
479
480

        # correct variance
481
        if preserve_gaussian_variance:
Martin Reinecke's avatar
Martin Reinecke committed
482
483
            assert issubclass(val.dtype.type, np.complexfloating),\
                    "complex input field is needed here"
484
485
486
            h *= np.sqrt(2)
            a *= np.sqrt(2)

487
488
        return (h, a)

489
490
    def _spec_to_rescaler(self, spec, result_list, power_space_index):
        power_space = self.domain[power_space_index]
491
492
493

        # weight the random fields with the power spectrum
        # therefore get the pindex from the power space
494
        pindex = power_space.pindex
495
496
497
498

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

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

Theo Steininger's avatar
Theo Steininger committed
510
    # ---Properties---
511

Theo Steininger's avatar
Theo Steininger committed
512
    def set_val(self, new_val=None, copy=False):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
513
        """ Sets the field's data object.
514
515
516

        Parameters
        ----------
Theo Steininger's avatar
Theo Steininger committed
517
        new_val : scalar, array-like, Field, None *optional*
518
519
            The values to be stored in the field.
            {default : None}
Theo Steininger's avatar
Theo Steininger committed
520

521
        copy : boolean, *optional*
Theo Steininger's avatar
Theo Steininger committed
522
523
            If False, Field tries to not copy the input data but use it
            directly.
524
525
526
527
528
529
            {default : False}
        See Also
        --------
        val

        """
Theo Steininger's avatar
Theo Steininger committed
530

531
532
        new_val = self.cast(new_val)
        if copy:
Theo Steininger's avatar
Theo Steininger committed
533
534
            new_val = new_val.copy()
        self._val = new_val
535
        return self
csongor's avatar
csongor committed
536

537
    def get_val(self, copy=False):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
538
        """ Returns the data object associated with this Field.
539
540
541
542

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

546
547
        Returns
        -------
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
548
        out : numpy.ndarray
549
550
551
552
553
554

        See Also
        --------
        val

        """
Theo Steininger's avatar
Theo Steininger committed
555

556
557
558
        if self._val is None:
            self.set_val(None)

559
        if copy:
Theo Steininger's avatar
Theo Steininger committed
560
            return self._val.copy()
561
        else:
Theo Steininger's avatar
Theo Steininger committed
562
            return self._val
csongor's avatar
csongor committed
563

Theo Steininger's avatar
Theo Steininger committed
564
565
    @property
    def val(self):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
566
        """ Returns the data object associated with this Field.
Theo Steininger's avatar
Theo Steininger committed
567

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

        See Also
        --------
        get_val

        """
Theo Steininger's avatar
Theo Steininger committed
577

578
        return self.get_val(copy=False)
csongor's avatar
csongor committed
579

Theo Steininger's avatar
Theo Steininger committed
580
581
    @val.setter
    def val(self, new_val):
582
        self.set_val(new_val=new_val, copy=False)
csongor's avatar
csongor committed
583

584
585
    @property
    def shape(self):
Theo Steininger's avatar
Theo Steininger committed
586
        """ Returns the total shape of the Field's data array.
Theo Steininger's avatar
Theo Steininger committed
587

588
589
590
        Returns
        -------
        out : tuple
Martin Reinecke's avatar
Martin Reinecke committed
591
            The output object. The tuple contains the dimensions of the spaces
592
593
594
595
596
597
598
            in domain.

        See Also
        --------
        dim

        """
Theo Steininger's avatar
Theo Steininger committed
599
600
601
602
603
604
605
606
        if not hasattr(self, '_shape'):
            shape_tuple = tuple(sp.shape for sp in self.domain)
            try:
                global_shape = reduce(lambda x, y: x + y, shape_tuple)
            except TypeError:
                global_shape = ()
            self._shape = global_shape
        return self._shape
csongor's avatar
csongor committed
607

608
609
    @property
    def dim(self):
Theo Steininger's avatar
Theo Steininger committed
610
        """ Returns the total number of pixel-dimensions the field has.
Theo Steininger's avatar
Theo Steininger committed
611

Theo Steininger's avatar
Theo Steininger committed
612
        Effectively, all values from shape are multiplied.
Theo Steininger's avatar
Theo Steininger committed
613

614
615
616
617
618
619
620
621
622
623
        Returns
        -------
        out : int
            The dimension of the Field.

        See Also
        --------
        shape

        """
Theo Steininger's avatar
Theo Steininger committed
624

625
        dim_tuple = tuple(sp.dim for sp in self.domain)
Theo Steininger's avatar
Theo Steininger committed
626
        try:
Martin Reinecke's avatar
Martin Reinecke committed
627
            return int(reduce(lambda x, y: x * y, dim_tuple))
Theo Steininger's avatar
Theo Steininger committed
628
629
        except TypeError:
            return 0
csongor's avatar
csongor committed
630

631
632
    @property
    def dof(self):
Theo Steininger's avatar
Theo Steininger committed
633
634
635
636
637
638
        """ Returns the total number of degrees of freedom the Field has. For
        real Fields this is equal to `self.dim`. For complex Fields it is
        2*`self.dim`.

        """

Theo Steininger's avatar
Theo Steininger committed
639
640
641
642
643
644
645
        dof = self.dim
        if issubclass(self.dtype.type, np.complexfloating):
            dof *= 2
        return dof

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

Theo Steininger's avatar
Theo Steininger committed
649
        volume_tuple = tuple(sp.total_volume for sp in self.domain)
650
        try:
Theo Steininger's avatar
Theo Steininger committed
651
            return reduce(lambda x, y: x * y, volume_tuple)
652
        except TypeError:
Theo Steininger's avatar
Theo Steininger committed
653
            return 0.
654

Theo Steininger's avatar
Theo Steininger committed
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
    @property
    def real(self):
        """ The real part of the field (data is not copied).
        """
        real_part = self.val.real
        result = self.copy_empty(dtype=real_part.dtype)
        result.set_val(new_val=real_part, copy=False)
        return result

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

Theo Steininger's avatar
Theo Steininger committed
673
    # ---Special unary/binary operations---
674

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

678
679
        Parameters
        ----------
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
680
681
682
        x : scalar, numpy.ndarray, Field, array_like
            The input that shall be casted on a numpy.ndarray of the same shape
            like the domain.
Theo Steininger's avatar
Theo Steininger committed
683

684
        dtype : type
Theo Steininger's avatar
Theo Steininger committed
685
            The datatype the output shall have. This can be used to override
Martin Reinecke's avatar
typos    
Martin Reinecke committed
686
            the field's dtype.
Theo Steininger's avatar
Theo Steininger committed
687

688
689
        Returns
        -------
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
690
        out : numpy.ndarray
691
692
693
694
695
696
697
            The output object.

        See Also
        --------
        _actual_cast

        """
csongor's avatar
csongor committed
698
699
        if dtype is None:
            dtype = self.dtype
700
701
        else:
            dtype = np.dtype(dtype)
702

703
704
        casted_x = x

705
        for ind, sp in enumerate(self.domain):
706
            casted_x = sp.pre_cast(casted_x,
707
708
709
                                   axes=self.domain_axes[ind])

        casted_x = self._actual_cast(casted_x, dtype=dtype)
710
711

        for ind, sp in enumerate(self.domain):
712
713
            casted_x = sp.post_cast(casted_x,
                                    axes=self.domain_axes[ind])
714

715
        return casted_x
csongor's avatar
csongor committed
716

Theo Steininger's avatar
Theo Steininger committed
717
    def _actual_cast(self, x, dtype=None):
718
        if isinstance(x, Field):
csongor's avatar
csongor committed
719
720
721
722
            x = x.get_val()

        if dtype is None:
            dtype = self.dtype
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
723
        if x is not None:
Martin Reinecke's avatar
more    
Martin Reinecke committed
724
725
            if np.isscalar(x):
                return np.full(self.shape,x, dtype=dtype)
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
726
727
728
            return np.asarray(x, dtype=dtype).reshape(self.shape)
        else:
            return np.empty(self.shape, dtype=dtype)
csongor's avatar
csongor committed
729

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

733
734
        If no keyword arguments are given, the returned object will be an
        identical copy of the original Field. By explicit specification one is
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
735
        able to define the domain and the dtype of the returned Field.
736
737
738
739
740

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

742
743
        dtype : type
            The new dtype the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
744

745
746
747
748
749
750
751
752
753
754
        Returns
        -------
        out : Field
            The output object. An identical copy of 'self'.

        See Also
        --------
        copy_empty

        """
Theo Steininger's avatar
Theo Steininger committed
755

Theo Steininger's avatar
Theo Steininger committed
756
        copied_val = self.get_val(copy=True)
757
758
        new_field = self.copy_empty(
                                domain=domain,
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
759
                                dtype=dtype)
Theo Steininger's avatar
Theo Steininger committed
760
761
        new_field.set_val(new_val=copied_val, copy=False)
        return new_field
csongor's avatar
csongor committed
762

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
763
    def copy_empty(self, domain=None, dtype=None):
764
765
766
        """ Returns an empty copy of the Field.

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

772
773
774
775
        Parameters
        ----------
        domain : DomainObject
            The new domain the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
776

777
778
        dtype : type
            The new dtype the Field shall have.
Theo Steininger's avatar
Theo Steininger committed
779

780
781
782
        Returns
        -------
        out : Field
Theo Steininger's avatar
Theo Steininger committed
783
            The output object.
784
785
786
787
788
789

        See Also
        --------
        copy

        """
Theo Steininger's avatar
Theo Steininger committed
790

Theo Steininger's avatar
Theo Steininger committed
791
792
        if domain is None:
            domain = self.domain
csongor's avatar
csongor committed
793
        else:
Theo Steininger's avatar
Theo Steininger committed
794
            domain = self._parse_domain(domain)
csongor's avatar
csongor committed
795

Theo Steininger's avatar
Theo Steininger committed
796
797
798
799
        if dtype is None:
            dtype = self.dtype
        else:
            dtype = np.dtype(dtype)
csongor's avatar
csongor committed
800

Theo Steininger's avatar
Theo Steininger committed
801
802
        fast_copyable = True
        try:
Martin Reinecke's avatar
Martin Reinecke committed
803
            for i in range(len(self.domain)):
Theo Steininger's avatar
Theo Steininger committed
804
805
806
807
808
809
                if self.domain[i] is not domain[i]:
                    fast_copyable = False
                    break
        except IndexError:
            fast_copyable = False

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
810
        if (fast_copyable and dtype == self.dtype):
Theo Steininger's avatar
Theo Steininger committed
811
812
            new_field = self._fast_copy_empty()
        else:
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
813
            new_field = Field(domain=domain, dtype=dtype)
Theo Steininger's avatar
Theo Steininger committed
814
        return new_field
csongor's avatar
csongor committed
815

Theo Steininger's avatar
Theo Steininger committed
816
817
818
819
820
821
    def _fast_copy_empty(self):
        # make an empty field
        new_field = EmptyField()
        # repair its class
        new_field.__class__ = self.__class__
        # copy domain, codomain and val
Martin Reinecke's avatar
Martin Reinecke committed
822
        for key, value in list(self.__dict__.items()):
823
            if key != '_val':
Theo Steininger's avatar
Theo Steininger committed
824
825
                new_field.__dict__[key] = value
            else:
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
826
                new_field.__dict__[key] = np.empty_like(self.val)
Theo Steininger's avatar
Theo Steininger committed
827
828
829
        return new_field

    def weight(self, power=1, inplace=False, spaces=None):
Theo Steininger's avatar
Theo Steininger committed
830
        """ Weights the pixels of `self` with their invidual pixel-volume.
831
832
833
834

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

837
        inplace : boolean
Theo Steininger's avatar
Theo Steininger committed
838
839
            If True, `self` will be weighted and returned. Otherwise, a copy
            is made.
Theo Steininger's avatar
Theo Steininger committed
840

Theo Steininger's avatar
Theo Steininger committed
841
842
        spaces : tuple of ints
            Determines on which subspace the operation takes place.
Theo Steininger's avatar
Theo Steininger committed
843

844
845
846
        Returns
        -------
        out : Field
Theo Steininger's avatar
Theo Steininger committed
847
            The weighted field.
848
849

        """
850
        if inplace:
csongor's avatar
csongor committed
851
852
853
854
            new_field = self
        else:
            new_field = self.copy_empty()

855
        new_val = self.get_val(copy=False)
csongor's avatar
csongor committed
856

857
        spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
csongor's avatar
csongor committed
858
        if spaces is None:
Martin Reinecke's avatar
Martin Reinecke committed
859
            spaces = list(range(len(self.domain)))
csongor's avatar
csongor committed
860

861
        for ind, sp in enumerate(self.domain):
Theo Steininger's avatar
Theo Steininger committed
862
863
864
865
866
            if ind in spaces:
                new_val = sp.weight(new_val,
                                    power=power,
                                    axes=self.domain_axes[ind],
                                    inplace=inplace)
867
868

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

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

874
875
876
        Parameters
        ----------
        x : Field
Theo Steininger's avatar
Theo Steininger committed
877
            The domain of x must contain `self.domain`
Theo Steininger's avatar
Theo Steininger committed
878

Theo Steininger's avatar
Theo Steininger committed
879
880
881
        spaces : tuple of ints
            If the domain of `self` and `x` are not the same, `spaces` specfies
            the mapping.
Theo Steininger's avatar
Theo Steininger committed
882

883
        bare : boolean
Theo Steininger's avatar
Theo Steininger committed
884
            If true, no volume factors will be included in the computation.
Theo Steininger's avatar
Theo Steininger committed
885

886
887
888
        Returns
        -------
        out : float, complex
Theo Steininger's avatar
Theo Steininger committed
889

890
        """
891
892
893
        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
894

Martin Reinecke's avatar
Martin Reinecke committed
895
        # Compute the dot respecting the fact of discrete/continuous spaces
Theo Steininger's avatar
Theo Steininger committed
896
897
898
899
900
        if bare:
            y = self
        else:
            y = self.weight(power=1)

901
902
903
        if spaces is None:
            x_val = x.get_val(copy=False)
            y_val = y.get_val(copy=False)
Martin Reinecke's avatar
Martin Reinecke committed
904
            result = (y_val.conjugate() * x_val).sum()
905
906
907
908
            return result
        else:
            # create a diagonal operator which is capable of taking care of the
            # axes-matching
909
            from .operators.diagonal_operator import DiagonalOperator
910
911
912
913
914
915
            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
916

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

Theo Steininger's avatar
Theo Steininger committed
920
921
922
        Returns
        -------
        norm : scalar
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
923
            The L2-norm of the field values.
csongor's avatar
csongor committed
924
925

        """
Martin Reinecke's avatar
fixes    
Martin Reinecke committed
926
        return np.sqrt(np.abs(self.vdot(x=self)))
csongor's avatar
csongor committed
927
928

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

931
932
933
        Parameters
        ----------
        inplace : boolean
Theo Steininger's avatar
Theo Steininger committed
934
            Decides whether the conjugation should be performed inplace.
Theo Steininger's avatar
Theo Steininger committed
935

936
937
938
939
        Returns
        -------
        cc : field
            The complex conjugated field.
csongor's avatar
csongor committed
940
941

        """
Theo Steininger's avatar
Theo Steininger committed
942

csongor's avatar
csongor committed
943
944
945
946
947
        if inplace:
            work_field = self
        else:
            work_field = self.copy_empty()

948
        new_val = self.get_val(copy=False)
Theo Steininger's avatar
Theo Steininger committed
949
        new_val = new_val.conjugate()
950
        work_field.set_val(new_val=new_val, copy=False)
csongor's avatar
csongor committed
951
952
953

        return work_field

Theo Steininger's avatar
Theo Steininger committed
954
    # ---General unary/contraction methods---
955

Theo Steininger's avatar
Theo Steininger committed
956
    def __pos__(self):
957
958
959
        """ x.__pos__() <==> +x

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

961
        """