space.py 29.5 KB
Newer Older
1
2
# NIFTY (Numerical Information Field Theory) has been developed at the
# Max-Planck-Institute for Astrophysics.
Marco Selig's avatar
Marco Selig committed
3
##
4
# Copyright (C) 2013 Max-Planck-Society
Marco Selig's avatar
Marco Selig committed
5
##
6
7
# Author: Marco Selig
# Project homepage: <http://www.mpa-garching.mpg.de/ift/nifty/>
Marco Selig's avatar
Marco Selig committed
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.
Marco Selig's avatar
Marco Selig committed
13
##
14
15
16
17
# 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.
Marco Selig's avatar
Marco Selig committed
18
##
19
20
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
Marco Selig's avatar
Marco Selig committed
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

"""
    ..                  __   ____   __
    ..                /__/ /   _/ /  /_
    ..      __ ___    __  /  /_  /   _/  __   __
    ..    /   _   | /  / /   _/ /  /   /  / /  /
    ..   /  / /  / /  / /  /   /  /_  /  /_/  /
    ..  /__/ /__/ /__/ /__/    \___/  \___   /  core
    ..                               /______/

    .. The NIFTY project homepage is http://www.mpa-garching.mpg.de/ift/nifty/

    NIFTY [#]_, "Numerical Information Field Theory", is a versatile
    library designed to enable the development of signal inference algorithms
    that operate regardless of the underlying spatial grid and its resolution.
    Its object-oriented framework is written in Python, although it accesses
    libraries written in Cython, C++, and C for efficiency.

    NIFTY offers a toolkit that abstracts discretized representations of
    continuous spaces, fields in these spaces, and operators acting on fields
    into classes. Thereby, the correct normalization of operations on fields is
    taken care of automatically without concerning the user. This allows for an
    abstract formulation and programming of inference algorithms, including
    those derived within information field theory. Thus, NIFTY permits its user
Marco Selig's avatar
Marco Selig committed
45
    to rapidly prototype algorithms in 1D and then apply the developed code in
Marco Selig's avatar
Marco Selig committed
46
47
48
49
50
    higher-dimensional settings of real world problems. The set of spaces on
    which NIFTY operates comprises point sets, n-dimensional regular grids,
    spherical spaces, their harmonic counterparts, and product spaces
    constructed as combinations of those.

51
52
53
54
55
56
57
    References
    ----------
    .. [#] Selig et al., "NIFTY -- Numerical Information Field Theory --
        a versatile Python library for signal inference",
        `A&A, vol. 554, id. A26 <http://dx.doi.org/10.1051/0004-6361/201321236>`_,
        2013; `arXiv:1301.4499 <http://www.arxiv.org/abs/1301.4499>`_

Marco Selig's avatar
Marco Selig committed
58
59
60
61
62
63
    Class & Feature Overview
    ------------------------
    The NIFTY library features three main classes: **spaces** that represent
    certain grids, **fields** that are defined on spaces, and **operators**
    that apply to fields.

64
65
    .. Overview of all (core) classes:
    ..
Marco Selig's avatar
Marco Selig committed
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
    .. - switch
    .. - notification
    .. - _about
    .. - random
    .. - space
    ..     - point_space
    ..     - rg_space
    ..     - lm_space
    ..     - gl_space
    ..     - hp_space
    ..     - nested_space
    .. - field
    .. - operator
    ..     - diagonal_operator
    ..         - power_operator
    ..     - projection_operator
    ..     - vecvec_operator
    ..     - response_operator
    .. - probing
    ..     - trace_probing
    ..     - diagonal_probing

88
89
    Overview of the main classes and functions:

Marco Selig's avatar
Marco Selig committed
90
91
    .. automodule:: nifty

92
93
94
95
96
97
98
99
100
101
102
103
104
105
    - :py:class:`space`
        - :py:class:`point_space`
        - :py:class:`rg_space`
        - :py:class:`lm_space`
        - :py:class:`gl_space`
        - :py:class:`hp_space`
        - :py:class:`nested_space`
    - :py:class:`field`
    - :py:class:`operator`
        - :py:class:`diagonal_operator`
            - :py:class:`power_operator`
        - :py:class:`projection_operator`
        - :py:class:`vecvec_operator`
        - :py:class:`response_operator`
Marco Selig's avatar
Marco Selig committed
106

107
        .. currentmodule:: nifty.nifty_tools
Marco Selig's avatar
Marco Selig committed
108

109
110
        - :py:class:`invertible_operator`
        - :py:class:`propagator_operator`
Marco Selig's avatar
Marco Selig committed
111

112
        .. currentmodule:: nifty.nifty_explicit
Marco Selig's avatar
Marco Selig committed
113

114
        - :py:class:`explicit_operator`
Marco Selig's avatar
Marco Selig committed
115

116
    .. automodule:: nifty
Marco Selig's avatar
Marco Selig committed
117

118
119
120
    - :py:class:`probing`
        - :py:class:`trace_probing`
        - :py:class:`diagonal_probing`
Marco Selig's avatar
Marco Selig committed
121

122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
        .. currentmodule:: nifty.nifty_explicit

        - :py:class:`explicit_probing`

    .. currentmodule:: nifty.nifty_tools

    - :py:class:`conjugate_gradient`
    - :py:class:`steepest_descent`

    .. currentmodule:: nifty.nifty_explicit

    - :py:func:`explicify`

    .. currentmodule:: nifty.nifty_power

    - :py:func:`weight_power`,
      :py:func:`smooth_power`,
      :py:func:`infer_power`,
      :py:func:`interpolate_power`
Marco Selig's avatar
Marco Selig committed
141
142
143
144

"""
from __future__ import division
import numpy as np
Marco Selig's avatar
Marco Selig committed
145
import pylab as pl
146

147
from d2o import STRATEGIES as DISTRIBUTION_STRATEGIES
148

149
from nifty_paradict import space_paradict
Ultimanet's avatar
Ultimanet committed
150

csongor's avatar
csongor committed
151
from nifty.config import about
152

153
POINT_DISTRIBUTION_STRATEGIES = DISTRIBUTION_STRATEGIES['global']
Marco Selig's avatar
Marco Selig committed
154

Ultimanet's avatar
Ultimanet committed
155

theos's avatar
theos committed
156
class Space(object):
Marco Selig's avatar
Marco Selig committed
157
    """
Ultimanet's avatar
Ultimanet committed
158
159
160
161
162
163
164
        ..                            __             __
        ..                          /__/           /  /_
        ..      ______    ______    __   __ ___   /   _/
        ..    /   _   | /   _   | /  / /   _   | /  /
        ..   /  /_/  / /  /_/  / /  / /  / /  / /  /_
        ..  /   ____/  \______/ /__/ /__/ /__/  \___/  space class
        .. /__/
Marco Selig's avatar
Marco Selig committed
165

Ultimanet's avatar
Ultimanet committed
166
        NIFTY subclass for unstructured spaces.
Marco Selig's avatar
Marco Selig committed
167

Ultimanet's avatar
Ultimanet committed
168
169
        Unstructured spaces are lists of values without any geometrical
        information.
Marco Selig's avatar
Marco Selig committed
170
171
172

        Parameters
        ----------
Ultimanet's avatar
Ultimanet committed
173
174
        num : int
            Number of points.
175
        dtype : numpy.dtype, *optional*
Ultimanet's avatar
Ultimanet committed
176
            Data type of the field values (default: None).
Marco Selig's avatar
Marco Selig committed
177

Ultimanet's avatar
Ultimanet committed
178
        Attributes
Marco Selig's avatar
Marco Selig committed
179
        ----------
Ultimanet's avatar
Ultimanet committed
180
181
        para : numpy.ndarray
            Array containing the number of points.
182
        dtype : numpy.dtype
Ultimanet's avatar
Ultimanet committed
183
184
185
186
187
188
            Data type of the field values.
        discrete : bool
            Parameter captioning the fact that a :py:class:`point_space` is
            always discrete.
        vol : numpy.ndarray
            Pixel volume of the :py:class:`point_space`, which is always 1.
Marco Selig's avatar
Marco Selig committed
189
    """
190

191
    def __init__(self, dtype=np.dtype('float'), **kwargs):
Ultimanet's avatar
Ultimanet committed
192
193
        """
            Sets the attributes for a point_space class instance.
Marco Selig's avatar
Marco Selig committed
194

Ultimanet's avatar
Ultimanet committed
195
196
197
198
            Parameters
            ----------
            num : int
                Number of points.
199
            dtype : numpy.dtype, *optional*
Ultimanet's avatar
Ultimanet committed
200
                Data type of the field values (default: numpy.float64).
Marco Selig's avatar
Marco Selig committed
201

Ultimanet's avatar
Ultimanet committed
202
203
204
205
            Returns
            -------
            None.
        """
Ultima's avatar
Ultima committed
206
        self._cache_dict = {'check_codomain': {}}
207
        self.paradict = space_paradict(**kwargs)
208

209
210
        # parse dtype
        dtype = np.dtype(dtype)
Ultima's avatar
Ultima committed
211
212
213
214
215
216
217
218
219
        if dtype not in [np.dtype('bool'),
                         np.dtype('int16'),
                         np.dtype('int32'),
                         np.dtype('int64'),
                         np.dtype('float32'),
                         np.dtype('float64'),
                         np.dtype('complex64'),
                         np.dtype('complex128')]:
            raise ValueError(about._errors.cstring(
220
                             "WARNING: incompatible dtype: " + str(dtype)))
Ultima's avatar
Ultima committed
221
        self.dtype = dtype
222

223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
        self.harmonic = None
        self._distances = None

    @property
    def distances(self):
        return self._distances

    @distances.setter
    def distances(self, distances):
        naxes = len(self.shape)

        if np.isscalar(distances):
            distances = tuple(np.ones(naxes, dtype=np.float) * distances)
        elif not isinstance(distances, tuple):
            distances = tuple(distances)
            if len(distances) != naxes:
                raise ValueError(about._errors.cstring(
                    "ERROR: size mismatch ( " + str(np.size(distances)) +
                    " <> " + str(naxes) + " )."))
        self._distances = distances
Marco Selig's avatar
Marco Selig committed
243

Ultimanet's avatar
Ultimanet committed
244
245
246
247
    @property
    def para(self):
        temp = np.array([self.paradict['num']], dtype=int)
        return temp
248

Ultimanet's avatar
Ultimanet committed
249
250
    @para.setter
    def para(self, x):
Ultima's avatar
Ultima committed
251
        self.paradict['num'] = x[0]
252

Ultima's avatar
Ultima committed
253
254
255
256
    def __hash__(self):
        # Extract the identifying parts from the vars(self) dict.
        result_hash = 0
        for (key, item) in vars(self).items():
Ultima's avatar
Ultima committed
257
258
            if key in ['_cache_dict']:
                continue
Ultima's avatar
Ultima committed
259
260
261
            result_hash ^= item.__hash__() * hash(key)
        return result_hash

262
263
264
265
266
    def _identifier(self):
        # Extract the identifying parts from the vars(self) dict.
        temp = [(ii[0],
                 ((lambda x: x[1].__hash__() if x[0] == 'comm' else x)(ii)))
                for ii in vars(self).iteritems()
Ultima's avatar
Ultima committed
267
                if ii[0] not in ['_cache_dict']
268
269
270
271
                ]
        # Return the sorted identifiers as a tuple.
        return tuple(sorted(temp))

272
    def copy(self):
273
        return Space(dtype=self.dtype, **self.paradict.parameters)
274

Ultimanet's avatar
Ultimanet committed
275
276
    def getitem(self, data, key):
        return data[key]
Marco Selig's avatar
Marco Selig committed
277

Ultimanet's avatar
Ultimanet committed
278
    def setitem(self, data, update, key):
279
        data[key] = update
Marco Selig's avatar
Marco Selig committed
280

Ultimanet's avatar
Ultimanet committed
281
    def apply_scalar_function(self, x, function, inplace=False):
282
        return x.apply_scalar_function(function, inplace=inplace)
283

284
285
    @property
    def shape(self):
286
        return (self.paradict['num'],)
Marco Selig's avatar
Marco Selig committed
287

288
289
    @property
    def dim(self):
Ultimanet's avatar
Ultimanet committed
290
291
        """
            Computes the dimension of the space, i.e.\  the number of points.
Marco Selig's avatar
Marco Selig committed
292

Ultimanet's avatar
Ultimanet committed
293
294
295
296
297
            Parameters
            ----------
            split : bool, *optional*
                Whether to return the dimension as an array with one component
                or as a scalar (default: False).
Marco Selig's avatar
Marco Selig committed
298

Ultimanet's avatar
Ultimanet committed
299
300
301
302
303
            Returns
            -------
            dim : {int, numpy.ndarray}
                Dimension(s) of the space.
        """
304
        return reduce(lambda x, y: x * y, self.shape)
Marco Selig's avatar
Marco Selig committed
305

306
307
    @property
    def dof(self):
Ultimanet's avatar
Ultimanet committed
308
309
310
311
        """
            Computes the number of degrees of freedom of the space, i.e./  the
            number of points for real-valued fields and twice that number for
            complex-valued fields.
Marco Selig's avatar
Marco Selig committed
312

Ultimanet's avatar
Ultimanet committed
313
314
315
316
317
            Returns
            -------
            dof : int
                Number of degrees of freedom of the space.
        """
318
319
320
        dof = self.dim
        if issubclass(self.dtype.type, np.complexfloating):
            dof = dof * 2
Ultima's avatar
Ultima committed
321
        return dof
322

323
    @property
324
325
326
    def vol(self):
        collapsed_distances = [np.sum(x) for x in self.distances]
        return reduce(lambda x, y: x * y, collapsed_distances)
327
328
329
330

    @property
    def vol_split(self):
        return self.distances
Marco Selig's avatar
Marco Selig committed
331

332
333
    @property
    def meta_volume(self):
Marco Selig's avatar
Marco Selig committed
334
        """
335
            Calculates the meta volumes.
Ultimanet's avatar
Ultimanet committed
336

337
338
339
340
341
            The meta volumes are the volumes associated with each component of
            a field, taking into account field components that are not
            explicitly included in the array of field values but are determined
            by symmetry conditions. In the case of an :py:class:`rg_space`, the
            meta volumes are simply the pixel volumes.
Marco Selig's avatar
Marco Selig committed
342
343
344

            Parameters
            ----------
345
346
347
            total : bool, *optional*
                Whether to return the total meta volume of the space or the
                individual ones of each pixel (default: False).
Marco Selig's avatar
Marco Selig committed
348
349
350

            Returns
            -------
351
352
            mol : {numpy.ndarray, float}
                Meta volume of the pixels or the complete space.
Ultimanet's avatar
Ultimanet committed
353
        """
354
355
356
357
358
359
        return self.dim() * self.vol()

    @property
    def meta_volume_split(self):
        mol = self.cast(1, dtype=np.dtype('float'))
        return self.calc_weight(mol, power=1)
360

361
362
363
    def complement_cast(self, x, axis=None):
        return x

364
365
366
367
    def enforce_power(self, spec, **kwargs):
        """
            Raises an error since the power spectrum is ill-defined for point
            spaces.
Marco Selig's avatar
Marco Selig committed
368
        """
369
370
371
        raise AttributeError(about._errors.cstring(
            "ERROR: the definition of power spectra is ill-defined for " +
            "(unstructured) point spaces."))
Ultimanet's avatar
Ultimanet committed
372

373
    def _enforce_power_helper(self, spec, size, kindex):
374
375
376
377
        # TODO: Resolve this import by splitting nifty_core into nifty_space
        # and nifty_point_space
        from nifty_field import field

378
379
380
        # Now it's about to extract a powerspectrum from spec
        # First of all just extract a numpy array. The shape is cared about
        # later.
381
        dtype = np.dtype('float')
382
383
384
385
        # Case 1: spec is a function
        if callable(spec):
            # Try to plug in the kindex array in the function directly
            try:
386
                spec = np.array(spec(kindex), dtype=dtype)
387
388
389
390
391
            except:
                # Second try: Use a vectorized version of the function.
                # This is slower, but better than nothing
                try:
                    spec = np.array(np.vectorize(spec)(kindex),
392
                                    dtype=dtype)
393
394
395
396
397
398
399
400
401
402
                except:
                    raise TypeError(about._errors.cstring(
                        "ERROR: invalid power spectra function."))

        # Case 2: spec is a field:
        elif isinstance(spec, field):
            try:
                spec = spec.val.get_full_data()
            except AttributeError:
                spec = spec.val
403
            spec = spec.astype(dtype).flatten()
Marco Selig's avatar
Marco Selig committed
404

405
406
        # Case 3: spec is a scalar or something else:
        else:
407
            spec = np.array(spec, dtype=dtype).flatten()
408
409
410
411
412

        # Make some sanity checks
        # check finiteness
        if not np.all(np.isfinite(spec)):
            about.warnings.cprint("WARNING: infinite value(s).")
Marco Selig's avatar
Marco Selig committed
413

414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
        # check positivity (excluding null)
        if np.any(spec < 0):
            raise ValueError(about._errors.cstring(
                "ERROR: nonpositive value(s)."))
        if np.any(spec == 0):
            about.warnings.cprint("WARNING: nonpositive value(s).")

        # Set the size parameter
        if size is None:
            size = len(kindex)

        # Fix the size of the spectrum
        # If spec is singlevalued, expand it
        if np.size(spec) == 1:
            spec = spec * np.ones(size, dtype=spec.dtype)
        # If the size does not fit at all, throw an exception
        elif np.size(spec) < size:
            raise ValueError(about._errors.cstring("ERROR: size mismatch ( " +
                                                   str(np.size(spec)) + " < " +
                                                   str(size) + " )."))
        elif np.size(spec) > size:
            about.warnings.cprint("WARNING: power spectrum cut to size ( == " +
                                  str(size) + " ).")
            spec = spec[:size]

        return spec
Ultimanet's avatar
Ultimanet committed
440

441
    def check_codomain(self, codomain):
Ultima's avatar
Ultima committed
442
443
444
445
446
447
448
449
450
451
        check_dict = self._cache_dict['check_codomain']
        temp_id = id(codomain)
        if temp_id in check_dict:
            return check_dict[temp_id]
        else:
            temp_result = self._check_codomain(codomain)
            check_dict[temp_id] = temp_result
            return temp_result

    def _check_codomain(self, codomain):
Marco Selig's avatar
Marco Selig committed
452
        """
453
            Checks whether a given codomain is compatible to the space or not.
Marco Selig's avatar
Marco Selig committed
454
455
456

            Parameters
            ----------
457
458
            codomain : nifty.space
                Space to be checked for compatibility.
Marco Selig's avatar
Marco Selig committed
459
460
461

            Returns
            -------
462
463
            check : bool
                Whether or not the given codomain is compatible to the space.
Marco Selig's avatar
Marco Selig committed
464
        """
465
466
        if codomain is None:
            return False
467

468
        if not isinstance(codomain, Space):
469
            raise TypeError(about._errors.cstring(
Ultima's avatar
Ultima committed
470
                "ERROR: invalid input. The given input is not a nifty space."))
Ultimanet's avatar
Ultimanet committed
471

472
473
474
475
        if codomain == self:
            return True
        else:
            return False
Ultimanet's avatar
Ultimanet committed
476

477
478
479
480
481
    def get_codomain(self, **kwargs):
        """
            Generates a compatible codomain to which transformations are
            reasonable, in this case another instance of
            :py:class:`point_space` with the same properties.
Marco Selig's avatar
Marco Selig committed
482

483
484
485
486
487
            Returns
            -------
            codomain : nifty.point_space
                A compatible codomain.
        """
theos's avatar
theos committed
488
489
        raise NotImplementedError(about._errors.cstring(
            "ERROR: There is no generic codomain for the Space base class."))
Marco Selig's avatar
Marco Selig committed
490

491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
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
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
#    def get_random_values(self, **kwargs):
#        """
#            Generates random field values according to the specifications given
#            by the parameters.
#
#            Returns
#            -------
#            x : numpy.ndarray
#                Valid field values.
#
#            Other parameters
#            ----------------
#            random : string, *optional*
#                Specifies the probability distribution from which the random
#                numbers are to be drawn.
#                Supported distributions are:
#
#                - "pm1" (uniform distribution over {+1,-1} or {+1,+i,-1,-i}
#                - "gau" (normal distribution with zero-mean and a given
#                standard
#                    deviation or variance)
#                - "syn" (synthesizes from a given power spectrum)
#                - "uni" (uniform distribution over [vmin,vmax[)
#
#                (default: None).
#            dev : float, *optional*
#                Standard deviation (default: 1).
#            var : float, *optional*
#                Variance, overriding `dev` if both are specified
#                (default: 1).
#            spec : {scalar, list, numpy.ndarray, nifty.field, function},
#            *optional*
#                Power spectrum (default: 1).
#            pindex : numpy.ndarray, *optional*
#                Indexing array giving the power spectrum index of each band
#                (default: None).
#            kindex : numpy.ndarray, *optional*
#                Scale of each band (default: None).
#            codomain : nifty.space, *optional*
#                A compatible codomain with power indices (default: None).
#            log : bool, *optional*
#                Flag specifying if the spectral binning is performed on
#                logarithmic
#                scale or not; if set, the number of used bins is set
#                automatically (if not given otherwise); by default no binning
#                is done (default: None).
#            nbin : integer, *optional*
#                Number of used spectral bins; if given `log` is set to
#                ``False``;
#                integers below the minimum of 3 induce an automatic setting;
#                by default no binning is done (default: None).
#            binbounds : {list, array}, *optional*
#                User specific inner boundaries of the bins, which are preferred
#                over the above parameters; by default no binning is done
#                (default: None).
#                vmin : {scalar, list, ndarray, field}, *optional*
#                Lower limit of the uniform distribution if ``random == "uni"``
#                (default: 0).
#            vmin : float, *optional*
#                Lower limit for a uniform distribution (default: 0).
#            vmax : float, *optional*
#                Upper limit for a uniform distribution (default: 1).
#        """
#
#        arg = random.parse_arguments(self, **kwargs)
#
#        if arg is None:
#            return self.cast(0)
#
#        # Prepare the empty distributed_data_object
#        sample = distributed_data_object(
#                                    global_shape=self.shape,
#                                    dtype=self.dtype)
#
#        # Case 1: uniform distribution over {-1,+1}/{1,i,-1,-i}
#        if arg['random'] == 'pm1':
#            sample.apply_generator(lambda s: random.pm1(dtype=self.dtype,
#                                                        shape=s))
#
#        # Case 2: normal distribution with zero-mean and a given standard
#        #         deviation or variance
#        elif arg['random'] == 'gau':
#            std = arg['std']
#            if np.isscalar(std) or std is None:
#                processed_std = std
#            else:
#                try:
#                    processed_std = sample.distributor. \
#                        extract_local_data(std)
#                except(AttributeError):
#                    processed_std = std
#
#            sample.apply_generator(lambda s: random.gau(dtype=self.dtype,
#                                                        shape=s,
#                                                        mean=arg['mean'],
#                                                        std=processed_std))
#
#        # Case 3: uniform distribution
#        elif arg['random'] == 'uni':
#            sample.apply_generator(lambda s: random.uni(dtype=self.dtype,
#                                                        shape=s,
#                                                        vmin=arg['vmin'],
#                                                        vmax=arg['vmax']))
#        return sample
595

596
    def calc_weight(self, x, power=1):
Marco Selig's avatar
Marco Selig committed
597
        """
Ultimanet's avatar
Ultimanet committed
598
599
            Weights a given array of field values with the pixel volumes (not
            the meta volumes) to a given power.
Marco Selig's avatar
Marco Selig committed
600
601
602

            Parameters
            ----------
Ultimanet's avatar
Ultimanet committed
603
604
605
606
            x : numpy.ndarray
                Array to be weighted.
            power : float, *optional*
                Power of the pixel volumes to be used (default: 1).
Marco Selig's avatar
Marco Selig committed
607
608

            Returns
Ultimanet's avatar
Ultimanet committed
609
610
611
            -------
            y : numpy.ndarray
                Weighted array.
Marco Selig's avatar
Marco Selig committed
612
        """
613
614
615
        # weight
        return x * self.get_weight(power=power)

616
    def get_weight(self, power=1, split=False):
617
618
        splitted_weight = tuple(np.array(self.distances)**np.array(power))
        if not split:
619
            return reduce(lambda x, y: x * y, splitted_weight)
620
621
        else:
            return splitted_weight
Marco Selig's avatar
Marco Selig committed
622

Ultima's avatar
Ultima committed
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
    def calc_norm(self, x, q=2):
        """
            Computes the Lq-norm of field values.

            Parameters
            ----------
            x : np.ndarray
                The data array
            q : scalar
                Parameter q of the Lq-norm (default: 2).

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

        """
        if q == 2:
            result = self.calc_dot(x, x)
        else:
            y = x**(q - 1)
            result = self.calc_dot(x, y)

        result = result**(1. / q)
        return result

649
    def dot_contraction(self, x, axes):
Ultimanet's avatar
Ultimanet committed
650
651
652
        """
            Computes the discrete inner product of two given arrays of field
            values.
Marco Selig's avatar
Marco Selig committed
653

Ultimanet's avatar
Ultimanet committed
654
655
656
657
658
659
            Parameters
            ----------
            x : numpy.ndarray
                First array
            y : numpy.ndarray
                Second array
Marco Selig's avatar
Marco Selig committed
660

Ultimanet's avatar
Ultimanet committed
661
662
663
664
665
            Returns
            -------
            dot : scalar
                Inner product of the two arrays.
        """
666
        return x.sum(axis=axes)
667

668
    def calc_transform(self, x, codomain=None, **kwargs):
Marco Selig's avatar
Marco Selig committed
669
        """
Ultimanet's avatar
Ultimanet committed
670
            Computes the transform of a given array of field values.
Marco Selig's avatar
Marco Selig committed
671
672
673

            Parameters
            ----------
Ultimanet's avatar
Ultimanet committed
674
675
676
            x : numpy.ndarray
                Array to be transformed.
            codomain : nifty.space, *optional*
677
                codomain space to which the transformation shall map
Ultimanet's avatar
Ultimanet committed
678
                (default: self).
Marco Selig's avatar
Marco Selig committed
679
680
681

            Returns
            -------
Ultimanet's avatar
Ultimanet committed
682
683
            Tx : numpy.ndarray
                Transformed array
684

Ultimanet's avatar
Ultimanet committed
685
            Other parameters
686
            ----------------
687
688
            iter : int, *optional*
                Number of iterations performed in specific transformations.
Marco Selig's avatar
Marco Selig committed
689
        """
690
691
692
        raise AttributeError(about._errors.cstring(
            "ERROR: fourier-transformation is ill-defined for " +
            "(unstructured) point space."))
693

694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
    def calc_smooth(self, x, **kwargs):
        """
            Raises an error since smoothing is ill-defined on an unstructured
            space.
        """
        raise AttributeError(about._errors.cstring(
            "ERROR: smoothing ill-defined for (unstructured) point space."))

    def calc_power(self, x, **kwargs):
        """
            Raises an error since the power spectrum is ill-defined for point
            spaces.
        """
        raise AttributeError(about._errors.cstring(
            "ERROR: power spectra ill-defined for (unstructured) " +
            "point space."))

    def calc_real_Q(self, x):
        try:
            return x.isreal().all()
714
        except AttributeError:
715
            return np.all(np.isreal(x))
Marco Selig's avatar
Marco Selig committed
716

717
    def calc_bincount(self, x, weights=None, minlength=None):
718
719
720
721
722
723
        try:
            complex_weights_Q = issubclass(weights.dtype.type,
                                           np.complexfloating)
        except AttributeError:
            complex_weights_Q = False

724
725
726
727
728
729
        if complex_weights_Q:
            real_bincount = x.bincount(weights=weights.real,
                                       minlength=minlength)
            imag_bincount = x.bincount(weights=weights.imag,
                                       minlength=minlength)
            return real_bincount + imag_bincount
730
        else:
731
            return x.bincount(weights=weights, minlength=minlength)
Marco Selig's avatar
Marco Selig committed
732

733
734
    def get_plot(self, x, title="", vmin=None, vmax=None, unit=None,
                 norm=None, other=None, legend=False, save=None, **kwargs):
Marco Selig's avatar
Marco Selig committed
735
        """
736
737
            Creates a plot of field values according to the specifications
            given by the parameters.
Marco Selig's avatar
Marco Selig committed
738
739
740

            Parameters
            ----------
Ultimanet's avatar
Ultimanet committed
741
            x : numpy.ndarray
742
                Array containing the field values.
Marco Selig's avatar
Marco Selig committed
743
744
745

            Returns
            -------
746
            None
747

Ultimanet's avatar
Ultimanet committed
748
            Other parameters
749
            ----------------
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
            title : string, *optional*
                Title of the plot (default: "").
            vmin : float, *optional*
                Minimum value to be displayed (default: ``min(x)``).
            vmax : float, *optional*
                Maximum value to be displayed (default: ``max(x)``).
            unit : string, *optional*
                Unit of the field values (default: "").
            norm : string, *optional*
                Scaling of the field values before plotting (default: None).
            other : {single object, tuple of objects}, *optional*
                Object or tuple of objects to be added, where objects can be
                scalars, arrays, or fields (default: None).
            legend : bool, *optional*
                Whether to show the legend or not (default: False).
            save : string, *optional*
                Valid file name where the figure is to be stored, by default
                the figure is not saved (default: False).

Ultimanet's avatar
Ultimanet committed
769
        """
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
        if not pl.isinteractive() and save is not None:
            about.warnings.cprint("WARNING: interactive mode off.")

        x = self.cast(x)

        fig = pl.figure(num=None,
                        figsize=(6.4, 4.8),
                        dpi=None,
                        facecolor="none",
                        edgecolor="none",
                        frameon=False,
                        FigureClass=pl.Figure)

        ax0 = fig.add_axes([0.12, 0.12, 0.82, 0.76])
        xaxes = np.arange(self.para[0], dtype=np.dtype('int'))

786
        if (norm == "log") and (vmin <= 0):
787
788
789
790
791
792
793
794
            raise ValueError(about._errors.cstring(
                "ERROR: nonpositive value(s)."))

        if issubclass(self.dtype.type, np.complexfloating):
            if vmin is None:
                vmin = min(x.real.min(), x.imag.min(), abs(x).min())
            if vmax is None:
                vmax = min(x.real.max(), x.imag.max(), abs(x).max())
Ultimanet's avatar
Ultimanet committed
795
        else:
796
797
798
799
            if vmin is None:
                vmin = x.min()
            if vmax is None:
                vmax = x.max()
Ultimanet's avatar
Ultimanet committed
800

801
802
803
        ax0.set_xlim(xaxes[0], xaxes[-1])
        ax0.set_xlabel("index")
        ax0.set_ylim(vmin, vmax)
804

805
806
        if(norm == "log"):
            ax0.set_yscale('log')
Marco Selig's avatar
Marco Selig committed
807

808
809
810
811
812
813
814
815
816
817
818
819
820
821
        if issubclass(self.dtype.type, np.complexfloating):
            ax0.scatter(xaxes, self.unary_operation(x, op='abs'),
                        color=[0.0, 0.5, 0.0], marker='o',
                        label="graph (absolute)", facecolor="none", zorder=1)
            ax0.scatter(xaxes, self.unary_operation(x, op='real'),
                        color=[0.0, 0.5, 0.0], marker='s',
                        label="graph (real part)", facecolor="none", zorder=1)
            ax0.scatter(xaxes, self.unary_operation(x, op='imag'),
                        color=[0.0, 0.5, 0.0], marker='D',
                        label="graph (imaginary part)", facecolor="none",
                        zorder=1)
        else:
            ax0.scatter(xaxes, x, color=[0.0, 0.5, 0.0], marker='o',
                        label="graph 0", zorder=1)
Marco Selig's avatar
Marco Selig committed
822

823
824
825
826
827
828
829
830
831
832
833
834
        if other is not None:
            if not isinstance(other, tuple):
                other = (other, )
            imax = max(1, len(other) - 1)
            for ii in xrange(len(other)):
                ax0.scatter(xaxes, self.dtype(other[ii]),
                            color=[max(0.0, 1.0 - (2 * ii / imax)**2),
                                   0.5 * ((2 * ii - imax) / imax)**2,
                                   max(0.0, 1.0 -
                                       (2 * (ii - imax) / imax)**2)],
                            marker='o', label="'other' graph " + str(ii),
                            zorder=-ii)
Ultimanet's avatar
Ultimanet committed
835

836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
        if legend:
            ax0.legend()

        if unit is not None:
            unit = " [" + unit + "]"
        else:
            unit = ""

        ax0.set_ylabel("values" + unit)
        ax0.set_title(title)

        if save is not None:
            fig.savefig(str(save), dpi=None,
                        facecolor="none", edgecolor="none")
            pl.close(fig)
        else:
            fig.canvas.draw()
Marco Selig's avatar
Marco Selig committed
853

854
    def __repr__(self):
Ultima's avatar
Ultima committed
855
856
        string = ""
        string += str(type(self)) + "\n"
Ultima's avatar
Ultima committed
857
        string += "paradict: " + str(self.paradict) + "\n"
Ultima's avatar
Ultima committed
858
859
860
        string += 'dtype: ' + str(self.dtype) + "\n"
        string += 'distances: ' + str(self.distances) + "\n"
        return string