space.py 29.7 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
theos's avatar
theos committed
324
325
326
    def volume_total(self):
        raise NotImplementedError(about._errors.cstring(
            "ERROR: There is no generic volume for the Space base class."))
327
328

    @property
theos's avatar
theos committed
329
330
331
    def volume_pixel(self):
        raise NotImplementedError(about._errors.cstring(
            "ERROR: There is no generic volume for the Space base class."))
Marco Selig's avatar
Marco Selig committed
332

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

338
339
340
341
342
            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
343
344
345

            Parameters
            ----------
346
347
348
            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
349
350
351

            Returns
            -------
352
353
            mol : {numpy.ndarray, float}
                Meta volume of the pixels or the complete space.
Ultimanet's avatar
Ultimanet committed
354
        """
355
356
357
358
359
360
        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)
361

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

365
366
367
368
    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
369
        """
370
371
372
        raise AttributeError(about._errors.cstring(
            "ERROR: the definition of power spectra is ill-defined for " +
            "(unstructured) point spaces."))
Ultimanet's avatar
Ultimanet committed
373

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

379
380
381
        # Now it's about to extract a powerspectrum from spec
        # First of all just extract a numpy array. The shape is cared about
        # later.
382
        dtype = np.dtype('float')
383
384
385
386
        # Case 1: spec is a function
        if callable(spec):
            # Try to plug in the kindex array in the function directly
            try:
387
                spec = np.array(spec(kindex), dtype=dtype)
388
389
390
391
392
            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),
393
                                    dtype=dtype)
394
395
396
397
398
399
400
401
402
403
                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
404
            spec = spec.astype(dtype).flatten()
Marco Selig's avatar
Marco Selig committed
405

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

        # 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
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
440
        # 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
441

442
    def check_codomain(self, codomain):
Ultima's avatar
Ultima committed
443
444
445
446
447
448
449
450
451
452
        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
453
        """
454
            Checks whether a given codomain is compatible to the space or not.
Marco Selig's avatar
Marco Selig committed
455
456
457

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

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

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

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

478
479
480
481
482
    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
483

484
485
486
487
488
            Returns
            -------
            codomain : nifty.point_space
                A compatible codomain.
        """
theos's avatar
theos committed
489
490
        raise NotImplementedError(about._errors.cstring(
            "ERROR: There is no generic codomain for the Space base class."))
Marco Selig's avatar
Marco Selig committed
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
595
#    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
596

597
    def calc_weight(self, x, power=1):
Marco Selig's avatar
Marco Selig committed
598
        """
Ultimanet's avatar
Ultimanet committed
599
600
            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
601
602
603

            Parameters
            ----------
Ultimanet's avatar
Ultimanet committed
604
605
606
607
            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
608
609

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

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

Ultima's avatar
Ultima committed
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
649
    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

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

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

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

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

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

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

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

695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
    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()
715
        except AttributeError:
716
            return np.all(np.isreal(x))
Marco Selig's avatar
Marco Selig committed
717

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

725
726
727
728
729
730
        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
731
        else:
732
            return x.bincount(weights=weights, minlength=minlength)
Marco Selig's avatar
Marco Selig committed
733

734
735
    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
736
        """
737
738
            Creates a plot of field values according to the specifications
            given by the parameters.
Marco Selig's avatar
Marco Selig committed
739
740
741

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

            Returns
            -------
747
            None
748

Ultimanet's avatar
Ultimanet committed
749
            Other parameters
750
            ----------------
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
            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
770
        """
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
        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'))

787
        if (norm == "log") and (vmin <= 0):
788
789
790
791
792
793
794
795
            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
796
        else:
797
798
799
800
            if vmin is None:
                vmin = x.min()
            if vmax is None:
                vmax = x.max()
Ultimanet's avatar
Ultimanet committed
801

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

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

809
810
811
812
813
814
815
816
817
818
819
820
821
822
        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
823

824
825
826
827
828
829
830
831
832
833
834
835
        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
836

837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
        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
854

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