nifty_lm.py 84.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
# NIFTY (Numerical Information Field Theory) has been developed at the
# Max-Planck-Institute for Astrophysics.
#
# Copyright (C) 2015 Max-Planck-Society
#
# Author: Marco Selig
# Project homepage: <http://www.mpa-garching.mpg.de/ift/nifty/>
#
# 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/>.
Marco Selig's avatar
Marco Selig committed
21
22
23
24
25
26
27
28
29
30
31
32
33
34

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

    NIFTY submodule for grids on the two-sphere.

"""
from __future__ import division
35

Marco Selig's avatar
Marco Selig committed
36
37
import os
import numpy as np
38
from numpy import pi
Marco Selig's avatar
Marco Selig committed
39
40
41
import pylab as pl
from matplotlib.colors import LogNorm as ln
from matplotlib.ticker import LogFormatter as lf
42
43
44
45
46
47
48

from nifty.nifty_core import space,\
                             point_space,\
                             field
from keepers import about,\
                    global_configuration as gc,\
                    global_dependency_injector as gdi
Ultimanet's avatar
Ultimanet committed
49
from nifty.nifty_paradict import lm_space_paradict,\
50
51
52
53
                                 gl_space_paradict,\
                                 hp_space_paradict
from nifty.nifty_power_indices import lm_power_indices

Ultimanet's avatar
Ultimanet committed
54
from nifty.nifty_random import random
55

56
57
58
59
gl = gdi['libsharp_wrapper_gl']
hp = gdi['healpy']

if gl is None and gc['lm2gl']:
Marco Selig's avatar
Marco Selig committed
60
    about.warnings.cprint("WARNING: global setting 'about.lm2gl' corrected.")
61
    gc['lm2gl'] = False
Marco Selig's avatar
Marco Selig committed
62

63
LM_DISTRIBUTION_STRATEGIES = []
Marco Selig's avatar
Marco Selig committed
64
65


66
class lm_space(point_space):
Marco Selig's avatar
Marco Selig committed
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
    """
        ..       __
        ..     /  /
        ..    /  /    __ ____ ___
        ..   /  /   /   _    _   |
        ..  /  /_  /  / /  / /  /
        ..  \___/ /__/ /__/ /__/  space class

        NIFTY subclass for spherical harmonics components, for representations
        of fields on the two-sphere.

        Parameters
        ----------
        lmax : int
            Maximum :math:`\ell`-value up to which the spherical harmonics
            coefficients are to be used.
        mmax : int, *optional*
            Maximum :math:`m`-value up to which the spherical harmonics
            coefficients are to be used (default: `lmax`).
86
        dtype : numpy.dtype, *optional*
Marco Selig's avatar
Marco Selig committed
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
            Data type of the field values (default: numpy.complex128).

        See Also
        --------
        hp_space : A class for the HEALPix discretization of the sphere [#]_.
        gl_space : A class for the Gauss-Legendre discretization of the
            sphere [#]_.

        Notes
        -----
        Hermitian symmetry, i.e. :math:`a_{\ell -m} = \overline{a}_{\ell m}` is
        always assumed for the spherical harmonics components, i.e. only fields
        on the two-sphere with real-valued representations in position space
        can be handled.

        References
        ----------
        .. [#] K.M. Gorski et al., 2005, "HEALPix: A Framework for
               High-Resolution Discretization and Fast Analysis of Data
               Distributed on the Sphere", *ApJ* 622..759G.
        .. [#] M. Reinecke and D. Sverre Seljebotn, 2013, "Libsharp - spherical
               harmonic transforms revisited";
               `arXiv:1303.4945 <http://www.arxiv.org/abs/1303.4945>`_

        Attributes
        ----------
        para : numpy.ndarray
            One-dimensional array containing the two numbers `lmax` and
            `mmax`.
116
        dtype : numpy.dtype
Marco Selig's avatar
Marco Selig committed
117
118
119
120
121
122
123
            Data type of the field values.
        discrete : bool
            Parameter captioning the fact that an :py:class:`lm_space` is
            always discrete.
        vol : numpy.ndarray
            Pixel volume of the :py:class:`lm_space`, which is always 1.
    """
124
125
126

    def __init__(self, lmax, mmax=None, dtype=np.dtype('complex128'),
                 datamodel='np'):
Marco Selig's avatar
Marco Selig committed
127
128
129
130
131
132
133
134
135
136
137
        """
            Sets the attributes for an lm_space class instance.

            Parameters
            ----------
            lmax : int
                Maximum :math:`\ell`-value up to which the spherical harmonics
                coefficients are to be used.
            mmax : int, *optional*
                Maximum :math:`m`-value up to which the spherical harmonics
                coefficients are to be used (default: `lmax`).
138
            dtype : numpy.dtype, *optional*
Marco Selig's avatar
Marco Selig committed
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
                Data type of the field values (default: numpy.complex128).

            Returns
            -------
            None.

            Raises
            ------
            ImportError
                If neither the libsharp_wrapper_gl nor the healpy module are
                available.
            ValueError
                If input `nside` is invaild.

        """
154

155
156
157
158
        # check imports
        if 'libsharp_wrapper_gl' not in gdi and 'healpy' not in gdi:
            raise ImportError(about._errors.cstring(
                "ERROR: neither libsharp_wrapper_gl nor healpy available."))
159

160
        self.paradict = lm_space_paradict(lmax=lmax, mmax=mmax)
Marco Selig's avatar
Marco Selig committed
161

162
163
164
        # check data type
        dtype = np.dtype(dtype)
        if dtype not in [np.dtype('complex64'), np.dtype('complex128')]:
Marco Selig's avatar
Marco Selig committed
165
            about.warnings.cprint("WARNING: data type set to default.")
166
167
            dtype = np.dtype('complex128')
        self.dtype = dtype
168

169
        # set datamodel
170
171
172
173
        if datamodel not in ['np']:
            about.warnings.cprint("WARNING: datamodel set to default.")
            self.datamodel = 'np'
        else:
174
175
            self.datamodel = datamodel

Marco Selig's avatar
Marco Selig committed
176
        self.discrete = True
177
        self.harmonic = True
178
179
180
181
182
183
184
185
        self.distances = (np.float(1),)

        self.power_indices = lm_power_indices(
                    lmax=self.paradict['lmax'],
                    dim=self.get_dim(),
                    comm=self.comm,
                    datamodel=self.datamodel,
                    allowed_distribution_strategies=LM_DISTRIBUTION_STRATEGIES)
Marco Selig's avatar
Marco Selig committed
186

187
188
    @property
    def para(self):
189
        temp = np.array([self.paradict['lmax'],
190
191
                         self.paradict['mmax']], dtype=int)
        return temp
192

193
194
195
196
197
    @para.setter
    def para(self, x):
        self.paradict['lmax'] = x[0]
        self.paradict['mmax'] = x[1]

198
    def copy(self):
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
        return lm_space(lmax=self.paradict['lmax'],
                        mmax=self.paradict['mmax'],
                        dtype=self.dtype)

#    def _enforce_shape(self, x):
#        """
#            Shapes an array of valid field values correctly, according to the
#            specifications of the space instance.
#
#            Parameters
#            ----------
#            x : numpy.ndarray
#                Array containing the field values to be put into shape.
#
#            Returns
#            -------
#            y : numpy.ndarray
#                Correctly shaped array.
#        """
#        about.warnings.cprint("WARNING: enforce_shape is deprecated!")
#
#        x = np.array(x)
#        if(np.size(x) != self.get_dim(split=False)):
#            raise ValueError(about._errors.cstring("ERROR: dimension mismatch ( " +
#                                                   str(np.size(x)) + " <> " + str(self.get_dim(split=False)) + " )."))
#        elif(not np.all(np.array(np.shape(x))==self.get_dim(split=True))):
#            about.warnings.cprint("WARNING: reshaping forced.")
#
#        return x.reshape(self.get_dim(split=True), order='C')

#    def lmax(self):
#        """
#            Returns the maximum quantum number :math:`\ell`.
#
#            Returns
#            -------
#            lmax : int
#                Maximum quantum number :math:`\ell`.
#        """
#        return self.paradict['lmax']
#
#    def mmax(self):
#        """
#            Returns the maximum quantum number :math:`m`.
#
#            Returns
#            -------
#            mmax : int
#                Maximum quantum number :math:`m`.
#
#        """
#        return self.paradict['mmax']
Marco Selig's avatar
Marco Selig committed
251

252
    def get_shape(self):
Ultima's avatar
Ultima committed
253
254
        mmax = self.paradict['mmax']
        lmax = self.paradict['lmax']
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
        return (np.int((mmax + 1) * (lmax + 1) - ((lmax + 1) * lmax) // 2),)

# -> inherited
#    def get_dim(self, split=False):
#        """
#            Computes the dimension of the space, i.e.\  the number of spherical
#            harmonics components that are stored.
#
#            Parameters
#            ----------
#            split : bool, *optional*
#                Whether to return the dimension as an array with one component
#                or as a scalar (default: False).
#
#            Returns
#            -------
#            dim : {int, numpy.ndarray}
#                Number of spherical harmonics components.
#
#            Notes
#            -----
#            Due to the symmetry assumption, only the components with
#            non-negative :math:`m` are stored and only these components are
#            counted here.
#        """
#        # dim = (mmax+1)*(lmax-mmax/2+1)
#        if(split):
#            return self.get_shape()
#            # return
#            # np.array([(self.para[0]+1)*(self.para[1]+1)-(self.para[1]+1)*self.para[1]//2],dtype=np.int)
#        else:
#            return np.prod(self.get_shape())
#            # return
#            # (self.para[0]+1)*(self.para[1]+1)-(self.para[1]+1)*self.para[1]//2

    def get_dof(self, split=False):
Marco Selig's avatar
Marco Selig committed
291
292
293
294
295
296
297
298
299
300
301
302
303
304
        """
            Computes the number of degrees of freedom of the space, taking into
            account symmetry constraints and complex-valuedness.

            Returns
            -------
            dof : int
                Number of degrees of freedom of the space.

            Notes
            -----
            The number of degrees of freedom is reduced due to the hermitian
            symmetry, which is assumed for the spherical harmonics components.
        """
305
306
307
308
309
310
311
312
        # dof = 2*dim-(lmax+1) = (lmax+1)*(2*mmax+1)*(mmax+1)*mmax
        mmax = self.paradict['lmax']
        lmax = self.paradict['lmax']
        dof = np.int((lmax + 1) * (2 * mmax + 1) - (mmax + 1) * mmax)
        if split:
            return (dof, )
        else:
            return dof
Marco Selig's avatar
Marco Selig committed
313

314
    def get_meta_volume(self, split=False):
Marco Selig's avatar
Marco Selig committed
315
        """
316
            Calculates the meta volumes.
Marco Selig's avatar
Marco Selig committed
317

318
319
320
321
            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.
Marco Selig's avatar
Marco Selig committed
322
323
324

            Parameters
            ----------
325
326
327
            total : bool, *optional*
                Whether to return the total meta volume of the space or the
                individual ones of each field component (default: False).
Marco Selig's avatar
Marco Selig committed
328
329
330

            Returns
            -------
331
332
            mol : {numpy.ndarray, float}
                Meta volume of the field components or the complete space.
Marco Selig's avatar
Marco Selig committed
333

334
335
336
337
338
            Notes
            -----
            The spherical harmonics components with :math:`m=0` have meta
            volume 1, the ones with :math:`m>0` have meta volume 2, sinnce they
            each determine another component with negative :math:`m`.
Marco Selig's avatar
Marco Selig committed
339
        """
340
341
342
343
344
345
        if not split:
            return np.float(self.get_dof())
        else:
            mol = self.cast(1, dtype=np.float)
            mol[self.paradict['lmax'] + 1:] = 2  # redundant: (l,m) and (l,-m)
            return mol
Marco Selig's avatar
Marco Selig committed
346

347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
    # TODO: Extend to binning/log
    def enforce_power(self, spec):
        size = self.paradict['lmax'] + 1
        kindex = np.arange(size, dtype=np.array(self.distances).dtype)
        return self._enforce_power_helper(spec=spec,
                                          size=size,
                                          kindex=kindex)

    def _getlm(self):  # > compute all (l,m)
        index = np.arange(self.get_dim())
        n = 2 * self.paradict['lmax'] + 1
        m = np.ceil(
            (n - np.sqrt(n**2 - 8 * (index - self.paradict['lmax']))) / 2
                    ).astype(np.int)
        l = index - self.paradict['lmax'] * m + m * (m - 1) // 2
        return l, m

    def _enforce_values(self, x, extend=True):
Marco Selig's avatar
Marco Selig committed
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
        """
            Computes valid field values from a given object, taking into
            account data types, size, and hermitian symmetry.

            Parameters
            ----------
            x : {float, numpy.ndarray, nifty.field}
                Object to be transformed into an array of valid field values.

            Returns
            -------
            x : numpy.ndarray
                Array containing the valid field values.

            Other parameters
            ----------------
            extend : bool, *optional*
                Whether a scalar is extented to a constant array or not
                (default: True).
        """
385
386
387
388
389
        if(isinstance(x, field)):
            if(self == x.domain):
                if(self.dtype is not x.domain.dtype):
                    raise TypeError(about._errors.cstring("ERROR: inequal data types ( '" + str(
                        np.result_type(self.dtype)) + "' <> '" + str(np.result_type(x.domain.dtype)) + "' )."))
Marco Selig's avatar
Marco Selig committed
390
391
392
                else:
                    x = np.copy(x.val)
            else:
393
394
                raise ValueError(about._errors.cstring(
                    "ERROR: inequal domains."))
Marco Selig's avatar
Marco Selig committed
395
        else:
396
            if(np.size(x) == 1):
Marco Selig's avatar
Marco Selig committed
397
                if(extend):
398
399
                    x = self.dtype(
                        x) * np.ones(self.get_dim(split=True), dtype=self.dtype, order='C')
Marco Selig's avatar
Marco Selig committed
400
401
                else:
                    if(np.isscalar(x)):
402
                        x = np.array([x], dtype=self.dtype)
Marco Selig's avatar
Marco Selig committed
403
                    else:
404
                        x = np.array(x, dtype=self.dtype)
Marco Selig's avatar
Marco Selig committed
405
            else:
406
                x = self._enforce_shape(np.array(x, dtype=self.dtype))
Marco Selig's avatar
Marco Selig committed
407

408
        if(np.size(x) != 1)and(np.any(x.imag[:self.para[0] + 1] != 0)):
Marco Selig's avatar
Marco Selig committed
409
            about.warnings.cprint("WARNING: forbidden values reset.")
410
411
412
            x.real[:self.para[0] + 1] = np.absolute(x[:self.para[0] + 1]) * (np.sign(x.real[:self.para[0] + 1]) + (
                np.sign(x.real[:self.para[0] + 1]) == 0) * np.sign(x.imag[:self.para[0] + 1])).astype(np.int)
            x.imag[:self.para[0] + 1] = 0  # x.imag[l,m==0] = 0
Marco Selig's avatar
Marco Selig committed
413

414
        # check finiteness
Marco Selig's avatar
Marco Selig committed
415
416
417
418
419
        if(not np.all(np.isfinite(x))):
            about.warnings.cprint("WARNING: infinite value(s).")

        return x

420
    def get_random_values(self, **kwargs):
Marco Selig's avatar
Marco Selig committed
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
        """
            Generates random field values according to the specifications given
            by the parameters, taking into account complex-valuedness and
            hermitian symmetry.

            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.array, nifty.field, function}, *optional*
                Power spectrum (default: 1).
            vmin : float, *optional*
                Lower limit for a uniform distribution (default: 0).
            vmax : float, *optional*
                Upper limit for a uniform distribution (default: 1).
        """
457
        arg = random.parse_arguments(self, **kwargs)
Marco Selig's avatar
Marco Selig committed
458
459

        if(arg is None):
460
            return np.zeros(self.get_dim(split=True), dtype=self.dtype, order='C')
Marco Selig's avatar
Marco Selig committed
461

462
463
        elif(arg[0] == "pm1"):
            x = random.pm1(dtype=self.dtype, shape=self.get_dim(split=True))
Marco Selig's avatar
Marco Selig committed
464

465
466
467
        elif(arg[0] == "gau"):
            x = random.gau(dtype=self.dtype, shape=self.get_dim(
                split=True), mean=None, dev=arg[2], var=arg[3])
Marco Selig's avatar
Marco Selig committed
468

469
470
471
472
473
        elif(arg[0] == "syn"):
            lmax = self.para[0]  # lmax
            if(self.dtype == np.dtype('complex64')):
                if 'libsharp_wrapper_gl' in gdi:  # default
                    x = gl.synalm_f(arg[1], lmax=lmax, mmax=lmax)
Marco Selig's avatar
Marco Selig committed
474
                else:
475
476
                    x = hp.synalm(arg[1].astype(np.complex128), lmax=lmax, mmax=lmax).astype(
                        np.complex64)  # FIXME: `verbose` kwarg
Marco Selig's avatar
Marco Selig committed
477
            else:
478
479
480
                if 'healpy' in gdi:  # default
                    # FIXME: `verbose` kwarg
                    x = hp.synalm(arg[1], lmax=lmax, mmax=lmax)
Marco Selig's avatar
Marco Selig committed
481
                else:
482
                    x = gl.synalm(arg[1], lmax=lmax, mmax=lmax)
Marco Selig's avatar
Marco Selig committed
483
484
485

            return x

486
487
488
        elif(arg[0] == "uni"):
            x = random.uni(dtype=self.dtype, shape=self.get_dim(
                split=True), vmin=arg[1], vmax=arg[2])
Marco Selig's avatar
Marco Selig committed
489
490

        else:
491
492
            raise KeyError(about._errors.cstring(
                "ERROR: unsupported random key '" + str(arg[0]) + "'."))
Marco Selig's avatar
Marco Selig committed
493

494
495
496
497
        if(np.any(x.imag[:self.para[0] + 1] != 0)):
            x.real[:self.para[0] + 1] = np.absolute(x[:self.para[0] + 1]) * (np.sign(x.real[:self.para[0] + 1]) + (
                np.sign(x.real[:self.para[0] + 1]) == 0) * np.sign(x.imag[:self.para[0] + 1])).astype(np.int)
            x.imag[:self.para[0] + 1] = 0  # x.imag[l,m==0] = 0
Marco Selig's avatar
Marco Selig committed
498
499
500

        return x

501
    def check_codomain(self, codomain):
Marco Selig's avatar
Marco Selig committed
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
        """
            Checks whether a given codomain is compatible to the
            :py:class:`lm_space` or not.

            Parameters
            ----------
            codomain : nifty.space
                Space to be checked for compatibility.

            Returns
            -------
            check : bool
                Whether or not the given codomain is compatible to the space.

            Notes
            -----
            Compatible codomains are instances of :py:class:`lm_space`,
            :py:class:`gl_space`, and :py:class:`hp_space`.
        """
521
522
        if codomain is None:
            return False
523

524
        if(not isinstance(codomain, space)):
Marco Selig's avatar
Marco Selig committed
525
526
            raise TypeError(about._errors.cstring("ERROR: invalid input."))

527
528
529
        if self.datamodel is not codomain.datamodel:
            return False

530
        if(self == codomain):
Marco Selig's avatar
Marco Selig committed
531
532
            return True

533
534
535
536
        elif(isinstance(codomain, gl_space)):
            # lmax==mmax                         nlat==lmax+1
            # nlon==2*lmax+1
            if(self.para[0] == self.para[1])and(codomain.para[0] == self.para[0] + 1)and(codomain.para[1] == 2 * self.para[0] + 1):
Marco Selig's avatar
Marco Selig committed
537
538
539
540
                return True
            else:
                about.warnings.cprint("WARNING: unrecommended codomain.")

541
542
543
        elif(isinstance(codomain, hp_space)):
            # lmax==mmax                        3*nside-1==lmax
            if(self.para[0] == self.para[1])and(3 * codomain.para[0] - 1 == self.para[0]):
Marco Selig's avatar
Marco Selig committed
544
545
546
547
548
549
                return True
            else:
                about.warnings.cprint("WARNING: unrecommended codomain.")

        return False

550
    def get_codomain(self, coname=None, **kwargs):
Marco Selig's avatar
Marco Selig committed
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
        """
            Generates a compatible codomain to which transformations are
            reasonable, i.e.\  a pixelization of the two-sphere.

            Parameters
            ----------
            coname : string, *optional*
                String specifying a desired codomain (default: None).

            Returns
            -------
            codomain : nifty.space
                A compatible codomain.

            Notes
            -----
            Possible arguments for `coname` are ``'gl'`` in which case a Gauss-
            Legendre pixelization [#]_ of the sphere is generated, and ``'hp'``
            in which case a HEALPix pixelization [#]_ is generated.

            References
            ----------
            .. [#] K.M. Gorski et al., 2005, "HEALPix: A Framework for
                   High-Resolution Discretization and Fast Analysis of Data
                   Distributed on the Sphere", *ApJ* 622..759G.
            .. [#] M. Reinecke and D. Sverre Seljebotn, 2013, "Libsharp - spherical
                   harmonic transforms revisited";
                   `arXiv:1303.4945 <http://www.arxiv.org/abs/1303.4945>`_

        """
581
582
583
584
        if(coname == "gl")or(coname is None)and(about.lm2gl.status):  # order matters
            if(self.dtype == np.dtype('complex64')):
                # nlat,nlon = lmax+1,2*lmax+1
                return gl_space(self.para[0] + 1, nlon=2 * self.para[0] + 1, dtype=np.float32)
Marco Selig's avatar
Marco Selig committed
585
            else:
586
587
                # nlat,nlon = lmax+1,2*lmax+1
                return gl_space(self.para[0] + 1, nlon=2 * self.para[0] + 1, dtype=np.float64)
Marco Selig's avatar
Marco Selig committed
588

589
590
        elif(coname == "hp")or(coname is None)and(not about.lm2gl.status):
            return hp_space((self.para[0] + 1) // 3)  # nside = (lmax+1)/3
Marco Selig's avatar
Marco Selig committed
591
592

        else:
593
594
            raise ValueError(about._errors.cstring(
                "ERROR: unsupported or incompatible space '" + coname + "'."))
Marco Selig's avatar
Marco Selig committed
595
596
597



598
599
600
601
602
603
604
    def _dotlm(self, x, y):  # > compute inner product
        dot = np.sum(x.real[:self.para[0] + 1] *
                     y.real[:self.para[0] + 1], axis=None, dtype=None, out=None)
        dot += 2 * np.sum(x.real[self.para[0] + 1:] *
                          y.real[:self.para[0] + 1:], axis=None, dtype=None, out=None)
        dot += 2 * np.sum(x.imag[self.para[0] + 1:] *
                          y.imag[:self.para[0] + 1:], axis=None, dtype=None, out=None)
Marco Selig's avatar
Marco Selig committed
605
606
        return dot

607
    def calc_dot(self, x, y):
Marco Selig's avatar
Marco Selig committed
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
        """
            Computes the discrete inner product of two given arrays of field
            values.

            Parameters
            ----------
            x : numpy.ndarray
                First array
            y : numpy.ndarray
                Second array

            Returns
            -------
            dot : scalar
                Inner product of the two arrays.
        """
624
625
626
627
628
629
        x = self._enforce_shape(np.array(x, dtype=self.dtype))
        y = self._enforce_shape(np.array(y, dtype=self.dtype))
        # inner product
        if 'libsharp_wrapper_gl' in gdi:  # default
            if(self.dtype == np.dtype('complex64')):
                return gl.dotlm_f(x, y, lmax=self.para[0], mmax=self.para[1])
Marco Selig's avatar
Marco Selig committed
630
            else:
631
                return gl.dotlm(x, y, lmax=self.para[0], mmax=self.para[1])
Marco Selig's avatar
Marco Selig committed
632
        else:
633
            return self._dotlm(x, y)
Ultima's avatar
Ultima committed
634

635
    def calc_transform(self, x, codomain=None, **kwargs):
Marco Selig's avatar
Marco Selig committed
636
637
638
639
640
641
642
643
        """
            Computes the transform of a given array of field values.

            Parameters
            ----------
            x : numpy.ndarray
                Array to be transformed.
            codomain : nifty.space, *optional*
644
                codomain space to which the transformation shall map
Marco Selig's avatar
Marco Selig committed
645
646
647
648
649
650
651
                (default: self).

            Returns
            -------
            Tx : numpy.ndarray
                Transformed array
        """
652
        x = self._enforce_shape(np.array(x, dtype=self.dtype))
Marco Selig's avatar
Marco Selig committed
653
654

        if(codomain is None):
655
            return x  # T == id
Marco Selig's avatar
Marco Selig committed
656

657
658
        # check codomain
        self.check_codomain(codomain)  # a bit pointless
Marco Selig's avatar
Marco Selig committed
659

660
661
        if(self == codomain):
            return x  # T == id
Marco Selig's avatar
Marco Selig committed
662

663
664
665
666
667
        elif(isinstance(codomain, gl_space)):
            # transform
            if(self.dtype == np.dtype('complex64')):
                Tx = gl.alm2map_f(x, nlat=codomain.para[0], nlon=codomain.para[
                                  1], lmax=self.para[0], mmax=self.para[1], cl=False)
Marco Selig's avatar
Marco Selig committed
668
            else:
669
670
671
                Tx = gl.alm2map(x, nlat=codomain.para[0], nlon=codomain.para[
                                1], lmax=self.para[0], mmax=self.para[1], cl=False)
            # weight if discrete
Marco Selig's avatar
Marco Selig committed
672
            if(codomain.discrete):
673
                Tx = codomain.calc_weight(Tx, power=0.5)
Marco Selig's avatar
Marco Selig committed
674

675
676
677
678
679
        elif(isinstance(codomain, hp_space)):
            # transform
            Tx = hp.alm2map(x.astype(np.complex128), codomain.para[0], lmax=self.para[0], mmax=self.para[
                            1], pixwin=False, fwhm=0.0, sigma=None, invert=False, pol=True, inplace=False)  # FIXME: `verbose` kwarg
            # weight if discrete
Marco Selig's avatar
Marco Selig committed
680
            if(codomain.discrete):
681
                Tx = codomain.calc_weight(Tx, power=0.5)
Marco Selig's avatar
Marco Selig committed
682
683

        else:
684
685
            raise ValueError(about._errors.cstring(
                "ERROR: unsupported transformation."))
Marco Selig's avatar
Marco Selig committed
686

687
        return Tx.astype(codomain.dtype)
Marco Selig's avatar
Marco Selig committed
688

689
    def calc_smooth(self, x, sigma=0, **kwargs):
Marco Selig's avatar
Marco Selig committed
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
        """
            Smoothes an array of field values by convolution with a Gaussian
            kernel in position space.

            Parameters
            ----------
            x : numpy.ndarray
                Array of field values to be smoothed.
            sigma : float, *optional*
                Standard deviation of the Gaussian kernel, specified in units
                of length in position space; for testing: a sigma of -1 will be
                reset to a reasonable value (default: 0).

            Returns
            -------
            Gx : numpy.ndarray
                Smoothed array.
        """
708
709
710
        x = self._enforce_shape(np.array(x, dtype=self.dtype))
        # check sigma
        if(sigma == 0):
Marco Selig's avatar
Marco Selig committed
711
            return x
712
        elif(sigma == -1):
Marco Selig's avatar
Marco Selig committed
713
            about.infos.cprint("INFO: invalid sigma reset.")
714
715
            sigma = 4.5 / (self.para[0] + 1)  # sqrt(2)*pi/(lmax+1)
        elif(sigma < 0):
Marco Selig's avatar
Marco Selig committed
716
            raise ValueError(about._errors.cstring("ERROR: invalid sigma."))
717
718
719
720
        # smooth
        if 'healpy' in gdi:  # default
            # no overwrite
            return hp.smoothalm(x, fwhm=0.0, sigma=sigma, invert=False, pol=True, mmax=self.para[1], verbose=False, inplace=False)
Marco Selig's avatar
Marco Selig committed
721
        else:
722
723
            # no overwrite
            return gl.smoothalm(x, lmax=self.para[0], mmax=self.para[1], fwhm=0.0, sigma=sigma, overwrite=False)
Marco Selig's avatar
Marco Selig committed
724

725
    def calc_power(self, x, **kwargs):
Marco Selig's avatar
Marco Selig committed
726
727
728
729
730
731
732
733
734
735
736
737
738
739
        """
            Computes the power of an array of field values.

            Parameters
            ----------
            x : numpy.ndarray
                Array containing the field values of which the power is to be
                calculated.

            Returns
            -------
            spec : numpy.ndarray
                Power contained in the input array.
        """
740
741
742
743
744
        x = self._enforce_shape(np.array(x, dtype=self.dtype))
        # power spectrum
        if(self.dtype == np.dtype('complex64')):
            if 'libsharp_wrapper_gl' in gdi:  # default
                return gl.anaalm_f(x, lmax=self.para[0], mmax=self.para[1])
Marco Selig's avatar
Marco Selig committed
745
            else:
746
                return hp.alm2cl(x.astype(np.complex128), alms2=None, lmax=self.para[0], mmax=self.para[1], lmax_out=self.para[0], nspec=None).astype(np.float32)
Marco Selig's avatar
Marco Selig committed
747
        else:
748
749
            if 'healpy' in gdi:  # default
                return hp.alm2cl(x, alms2=None, lmax=self.para[0], mmax=self.para[1], lmax_out=self.para[0], nspec=None)
Marco Selig's avatar
Marco Selig committed
750
            else:
751
                return gl.anaalm(x, lmax=self.para[0], mmax=self.para[1])
Marco Selig's avatar
Marco Selig committed
752

753
    def get_plot(self, x, title="", vmin=None, vmax=None, power=True, norm=None, cmap=None, cbar=True, other=None, legend=False, mono=True, **kwargs):
Marco Selig's avatar
Marco Selig committed
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
        """
            Creates a plot of field values according to the specifications
            given by the parameters.

            Parameters
            ----------
            x : numpy.ndarray
                Array containing the field values.

            Returns
            -------
            None

            Other parameters
            ----------------
            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)``).
            power : bool, *optional*
                Whether to plot the power contained in the field or the field
                values themselves (default: True).
            unit : string, *optional*
                Unit of the field values (default: "").
            norm : string, *optional*
                Scaling of the field values before plotting (default: None).
            cmap : matplotlib.colors.LinearSegmentedColormap, *optional*
                Color map to be used for two-dimensional plots (default: None).
            cbar : bool, *optional*
                Whether to show the color bar or not (default: True).
            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).
            mono : bool, *optional*
                Whether to plot the monopole or not (default: True).
            save : string, *optional*
                Valid file name where the figure is to be stored, by default
                the figure is not saved (default: False).

        """
798
        if(not pl.isinteractive())and(not bool(kwargs.get("save", False))):
Marco Selig's avatar
Marco Selig committed
799
800
801
802
803
            about.warnings.cprint("WARNING: interactive mode off.")

        if(power):
            x = self.calc_power(x)

804
805
806
            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])
Marco Selig's avatar
Marco Selig committed
807

808
            xaxes = np.arange(self.para[0] + 1, dtype=np.int)
Marco Selig's avatar
Marco Selig committed
809
            if(vmin is None):
810
811
                vmin = np.min(x[:mono].tolist(
                ) + (xaxes * (2 * xaxes + 1) * x)[1:].tolist(), axis=None, out=None)
Marco Selig's avatar
Marco Selig committed
812
            if(vmax is None):
813
814
815
816
                vmax = np.max(x[:mono].tolist(
                ) + (xaxes * (2 * xaxes + 1) * x)[1:].tolist(), axis=None, out=None)
            ax0.loglog(xaxes[1:], (xaxes * (2 * xaxes + 1) * x)[1:], color=[0.0,
                                                                            0.5, 0.0], label="graph 0", linestyle='-', linewidth=2.0, zorder=1)
Marco Selig's avatar
Marco Selig committed
817
            if(mono):
818
819
                ax0.scatter(0.5 * (xaxes[1] + xaxes[2]), x[0], s=20, color=[0.0, 0.5, 0.0], marker='o',
                            cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, verts=None, zorder=1)
Marco Selig's avatar
Marco Selig committed
820
821

            if(other is not None):
822
                if(isinstance(other, tuple)):
Marco Selig's avatar
Marco Selig committed
823
824
                    other = list(other)
                    for ii in xrange(len(other)):
825
                        if(isinstance(other[ii], field)):
Marco Selig's avatar
Marco Selig committed
826
827
828
                            other[ii] = other[ii].power(**kwargs)
                        else:
                            other[ii] = self.enforce_power(other[ii])
829
                elif(isinstance(other, field)):
Marco Selig's avatar
Marco Selig committed
830
831
832
                    other = [other.power(**kwargs)]
                else:
                    other = [self.enforce_power(other)]
833
                imax = max(1, len(other) - 1)
Marco Selig's avatar
Marco Selig committed
834
                for ii in xrange(len(other)):
835
836
                    ax0.loglog(xaxes[1:], (xaxes * (2 * xaxes + 1) * other[ii])[1:], 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)], label="graph " + str(ii + 1), linestyle='-', linewidth=1.0, zorder=-ii)
Marco Selig's avatar
Marco Selig committed
837
                    if(mono):
838
839
                        ax0.scatter(0.5 * (xaxes[1] + xaxes[2]), other[ii][0], s=20, 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', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, verts=None, zorder=-ii)
Marco Selig's avatar
Marco Selig committed
840
841
842
                if(legend):
                    ax0.legend()

843
            ax0.set_xlim(xaxes[1], xaxes[-1])
Marco Selig's avatar
Marco Selig committed
844
            ax0.set_xlabel(r"$\ell$")
845
            ax0.set_ylim(vmin, vmax)
Marco Selig's avatar
Marco Selig committed
846
847
848
849
            ax0.set_ylabel(r"$\ell(2\ell+1) C_\ell$")
            ax0.set_title(title)

        else:
850
            x = self._enforce_shape(np.array(x))
Marco Selig's avatar
Marco Selig committed
851
852
853
            if(np.iscomplexobj(x)):
                if(title):
                    title += " "
854
855
856
857
858
859
                if(bool(kwargs.get("save", False))):
                    save_ = os.path.splitext(
                        os.path.basename(str(kwargs.get("save"))))
                    kwargs.update(save=save_[0] + "_absolute" + save_[1])
                self.get_plot(np.absolute(x), title=title + "(absolute)", vmin=vmin, vmax=vmax,
                              power=False, norm=norm, cmap=cmap, cbar=cbar, other=None, legend=False, **kwargs)
Marco Selig's avatar
Marco Selig committed
860
861
862
863
#                self.get_plot(np.real(x),title=title+"(real part)",vmin=vmin,vmax=vmax,power=False,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs)
#                self.get_plot(np.imag(x),title=title+"(imaginary part)",vmin=vmin,vmax=vmax,power=False,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs)
                if(cmap is None):
                    cmap = pl.cm.hsv_r
864
865
866
867
868
                if(bool(kwargs.get("save", False))):
                    kwargs.update(save=save_[0] + "_phase" + save_[1])
                self.get_plot(np.angle(x, deg=False), title=title + "(phase)", vmin=-3.1416, vmax=3.1416, power=False,
                              norm=None, cmap=cmap, cbar=cbar, other=None, legend=False, **kwargs)  # values in [-pi,pi]
                return None  # leave method
Marco Selig's avatar
Marco Selig committed
869
870
            else:
                if(vmin is None):
871
                    vmin = np.min(x, axis=None, out=None)
Marco Selig's avatar
Marco Selig committed
872
                if(vmax is None):
873
874
875
876
877
878
879
880
881
882
                    vmax = np.max(x, axis=None, out=None)
                if(norm == "log")and(vmin <= 0):
                    raise ValueError(about._errors.cstring(
                        "ERROR: nonpositive value(s)."))

                # not a number
                xmesh = np.nan * \
                    np.empty(self.para[::-1] + 1, dtype=np.float16, order='C')
                xmesh[4, 1] = None
                xmesh[1, 4] = None
Marco Selig's avatar
Marco Selig committed
883
                lm = 0
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
                for mm in xrange(self.para[1] + 1):
                    xmesh[mm][mm:] = x[lm:lm + self.para[0] + 1 - mm]
                    lm += self.para[0] + 1 - mm

                s_ = np.array([1, self.para[1] / self.para[0]
                               * (1.0 + 0.159 * bool(cbar))])
                fig = pl.figure(num=None, figsize=(
                    6.4 * s_[0], 6.4 * s_[1]), dpi=None, facecolor="none", edgecolor="none", frameon=False, FigureClass=pl.Figure)
                ax0 = fig.add_axes(
                    [0.06 / s_[0], 0.06 / s_[1], 1.0 - 0.12 / s_[0], 1.0 - 0.12 / s_[1]])
                ax0.set_axis_bgcolor([0.0, 0.0, 0.0, 0.0])

                xaxes = np.arange(self.para[0] + 2, dtype=np.int) - 0.5
                yaxes = np.arange(self.para[1] + 2, dtype=np.int) - 0.5
                if(norm == "log"):
                    n_ = ln(vmin=vmin, vmax=vmax)
Marco Selig's avatar
Marco Selig committed
900
901
                else:
                    n_ = None
902
903
904
905
                sub = ax0.pcolormesh(xaxes, yaxes, np.ma.masked_where(np.isnan(
                    xmesh), xmesh), cmap=cmap, norm=n_, vmin=vmin, vmax=vmax, clim=(vmin, vmax))
                ax0.set_xlim(xaxes[0], xaxes[-1])
                ax0.set_xticks([0], minor=False)
Marco Selig's avatar
Marco Selig committed
906
                ax0.set_xlabel(r"$\ell$")
907
908
                ax0.set_ylim(yaxes[0], yaxes[-1])
                ax0.set_yticks([0], minor=False)
Marco Selig's avatar
Marco Selig committed
909
910
911
                ax0.set_ylabel(r"$m$")
                ax0.set_aspect("equal")
                if(cbar):
912
913
914
915
916
917
                    if(norm == "log"):
                        f_ = lf(10, labelOnlyBase=False)
                        b_ = sub.norm.inverse(
                            np.linspace(0, 1, sub.cmap.N + 1))
                        v_ = np.linspace(
                            sub.norm.vmin, sub.norm.vmax, sub.cmap.N)
Marco Selig's avatar
Marco Selig committed
918
919
920
921
                    else:
                        f_ = None
                        b_ = None
                        v_ = None
922
923
                    fig.colorbar(sub, ax=ax0, orientation="horizontal", fraction=0.1, pad=0.05, shrink=0.75, aspect=20, ticks=[
                                 vmin, vmax], format=f_, drawedges=False, boundaries=b_, values=v_)
Marco Selig's avatar
Marco Selig committed
924
925
                ax0.set_title(title)

926
927
928
        if(bool(kwargs.get("save", False))):
            fig.savefig(str(kwargs.get("save")), dpi=None, facecolor="none", edgecolor="none", orientation="portrait",
                        papertype=None, format=None, transparent=False, bbox_inches=None, pad_inches=0.1)
Marco Selig's avatar
Marco Selig committed
929
930
931
932
933
934
935
936
            pl.close(fig)
        else:
            fig.canvas.draw()

    def __repr__(self):
        return "<nifty_lm.lm_space>"

    def __str__(self):
937
        return "nifty_lm.lm_space instance\n- lmax     = " + str(self.para[0]) + "\n- mmax     = " + str(self.para[1]) + "\n- dtype = numpy." + str(np.result_type(self.dtype))
Marco Selig's avatar
Marco Selig committed
938

939
# -----------------------------------------------------------------------------
Marco Selig's avatar
Marco Selig committed
940
941


942
# -----------------------------------------------------------------------------
Marco Selig's avatar
Marco Selig committed
943

944
class gl_space(point_space):
Marco Selig's avatar
Marco Selig committed
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
    """
        ..                 __
        ..               /  /
        ..     ____ __  /  /
        ..   /   _   / /  /
        ..  /  /_/  / /  /_
        ..  \___   /  \___/  space class
        .. /______/

        NIFTY subclass for Gauss-Legendre pixelizations [#]_ of the two-sphere.

        Parameters
        ----------
        nlat : int
            Number of latitudinal bins, or rings.
        nlon : int, *optional*
            Number of longitudinal bins (default: ``2*nlat - 1``).
962
        dtype : numpy.dtype, *optional*
Marco Selig's avatar
Marco Selig committed
963
964
965
966
967
968
969
970
971
972
            Data type of the field values (default: numpy.float64).

        See Also
        --------
        hp_space : A class for the HEALPix discretization of the sphere [#]_.
        lm_space : A class for spherical harmonic components.

        Notes
        -----
        Only real-valued fields on the two-sphere are supported, i.e.
973
        `dtype` has to be either numpy.float64 or numpy.float32.
Marco Selig's avatar
Marco Selig committed
974
975
976
977
978
979
980
981
982
983
984
985
986
987

        References
        ----------
        .. [#] M. Reinecke and D. Sverre Seljebotn, 2013, "Libsharp - spherical
               harmonic transforms revisited";
               `arXiv:1303.4945 <http://www.arxiv.org/abs/1303.4945>`_
        .. [#] K.M. Gorski et al., 2005, "HEALPix: A Framework for
               High-Resolution Discretization and Fast Analysis of Data
               Distributed on the Sphere", *ApJ* 622..759G.

        Attributes
        ----------
        para : numpy.ndarray
            One-dimensional array containing the two numbers `nlat` and `nlon`.
988
        dtype : numpy.dtype
Marco Selig's avatar
Marco Selig committed
989
990
991
992
993
994
995
            Data type of the field values.
        discrete : bool
            Whether or not the underlying space is discrete, always ``False``
            for spherical spaces.
        vol : numpy.ndarray
            An array containing the pixel sizes.
    """
996
997
998

    def __init__(self, nlat, nlon=None, dtype=np.dtype('float'),
                 datamodel='np'):
Marco Selig's avatar
Marco Selig committed
999
1000
1001
1002
1003
1004
1005
1006
1007
        """
            Sets the attributes for a gl_space class instance.

            Parameters
            ----------
            nlat : int
                Number of latitudinal bins, or rings.
            nlon : int, *optional*
                Number of longitudinal bins (default: ``2*nlat - 1``).
1008
            dtype : numpy.dtype, *optional*
Marco Selig's avatar
Marco Selig committed
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
                Data type of the field values (default: numpy.float64).

            Returns
            -------
            None

            Raises
            ------
            ImportError
                If the libsharp_wrapper_gl module is not available.
            ValueError
                If input `nlat` is invaild.

        """
1023
1024
1025
1026
        # check imports
        if 'libsharp_wrapper_gl' not in gdi:
            raise ImportError(about._errors.cstring(
                "ERROR: libsharp_wrapper_gl not available."))
1027
1028

        self.paradict = gl_space_paradict(nlat=nlat, nlon=nlon)
Marco Selig's avatar
Marco Selig committed
1029

1030
1031
1032
        # check data type
        dtype = np.dtype(dtype)
        if dtype not in [np.dtype('float32'), np.dtype('float64')]:
Marco Selig's avatar
Marco Selig committed
1033
            about.warnings.cprint("WARNING: data type set to default.")
1034
1035
            dtype = np.dtype('float')
        self.dtype = dtype
1036

1037
        # set datamodel
1038
1039
1040
1041
        if datamodel not in ['np']:
            about.warnings.cprint("WARNING: datamodel set to default.")
            self.datamodel = 'np'
        else:
1042
            self.datamodel = datamodel
Marco Selig's avatar
Marco Selig committed
1043
1044

        self.discrete = False
1045
        self.harmonic = False
1046
1047
1048
        self.distances = tuple(gl.vol(self.paradict['nlat'],
                                      nlon=self.paradict['nlon']
                                      ).astype(np.float))
1049
1050
1051

    @property
    def para(self):
1052
        temp = np.array([self.paradict['nlat'],
1053
1054
                         self.paradict['nlon']], dtype=int)
        return temp
1055

1056
1057
1058
1059
    @para.setter
    def para(self, x):
        self.paradict['nlat'] = x[0]
        self.paradict['nlon'] = x[1]
1060

1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
    def _enforce_shape(self, x):
        """
            Shapes an array of valid field values correctly, according to the
            specifications of the space instance.

            Parameters
            ----------
            x : numpy.ndarray
                Array containing the field values to be put into shape.

            Returns
            -------
            y : numpy.ndarray
                Correctly shaped array.
        """
        about.warnings.cprint("WARNING: enforce_shape is deprecated!")

        x = np.array(x)
        if(np.size(x) != self.get_dim(split=False)):
            raise ValueError(about._errors.cstring("ERROR: dimension mismatch ( " +
                                                   str(np.size(x)) + " <> " + str(self.get_dim(split=False)) + " )."))
#        elif(not np.all(np.array(np.shape(x))==self.get_dim(split=True))):
#            about.warnings.cprint("WARNING: reshaping forced.")

        return x.reshape(self.get_dim(split=True), order='C')
Marco Selig's avatar
Marco Selig committed
1086

1087
    def copy(self):
1088
1089
1090
1091
        return gl_space(nlat=self.paradict['nlat'],
                        nlon=self.paradict['nlon'],
                        dtype=self.dtype)

Marco Selig's avatar
Marco Selig committed
1092
1093
1094
1095
1096
1097
1098
1099
1100
    def nlat(self):
        """
            Returns the number of latitudinal bins.

            Returns
            -------
            nlat : int
                Number of latitudinal bins, or rings.
        """
1101
        return self.paradict['nlat']
Marco Selig's avatar
Marco Selig committed
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111

    def nlon(self):
        """
            Returns the number of longitudinal bins.

            Returns
            -------
            nlon : int
                Number of longitudinal bins.
        """
1112
        return self.paradict['nlon']
Marco Selig's avatar
Marco Selig committed
1113

1114
    def get_shape(self):
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
        return (np.int((self.paradict['nlat'] * self.paradict['nlon'])),)

# -> inherited
#    def get_dim(self, split=False):
#        """
#            Computes the dimension of the space, i.e.\  the number of pixels.
#
#            Parameters
#            ----------
#            split : bool, *optional*
#                Whether to return the dimension as an array with one component
#                or as a scalar (default: False).
#
#            Returns
#            -------
#            dim : {int, numpy.ndarray}
#                Dimension of the space.
#        """
#        # dim = nlat*nlon
#        if(split):
#            return self.get_shape()
#            # return np.array([self.para[0]*self.para[1]],dtype=np.int)
#        else:
#            return np.prod(self.get_shape())
#            # return self.para[0]*self.para[1]

    def get_dof(self, split=False):
Marco Selig's avatar
Marco Selig committed
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
        """
            Computes the number of degrees of freedom of the space.

            Returns
            -------
            dof : int
                Number of degrees of freedom of the space.

            Notes
            -----
            Since the :py:class:`gl_space` class only supports real-valued
            fields, the number of degrees of freedom is the number of pixels.
        """
1155
1156
        # dof = dim
        return self.get_dim(split=split)
Marco Selig's avatar
Marco Selig committed
1157

1158
1159
1160
1161
1162
1163
1164
    # TODO: Extend to binning/log
    def enforce_power(self, spec):
        size = self.paradict['nlat']
        kindex = np.arange(size, dtype=np.array(self.distances).dtype)
        return self._enforce_power_helper(spec=spec,
                                          size=size,
                                          kindex=kindex)
Marco Selig's avatar
Marco Selig committed
1165

1166
    def get_random_values(self, **kwargs):
Marco Selig's avatar
Marco Selig committed
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
        """
            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.array, nifty.field, function}, *optional*
                Power spectrum (default: 1).
            codomain : nifty.lm_space, *optional*
                A compatible codomain for power indexing (default: None).
            vmin : float, *optional*
                Lower limit for a uniform distribution (default: 0).
            vmax : float, *optional*
                Upper limit for a uniform distribution (default: 1).
        """
1204
        arg = random.parse_arguments(self, **kwargs)
Marco Selig's avatar
Marco Selig committed
1205
1206

        if(arg is None):
1207
            x = np.zeros(self.get_dim(split=True), dtype=self.dtype, order='C')
Marco Selig's avatar
Marco Selig committed
1208

1209
1210
        elif(arg[0] == "pm1"):
            x = random.pm1(dtype=self.dtype, shape=self.get_dim(split=True))
Marco Selig's avatar
Marco Selig committed
1211

1212
1213
1214
        elif(arg[0] == "gau"):
            x = random.gau(dtype=self.dtype, shape=self.get_dim(
                split=True), mean=None, dev=arg[2], var=arg[3])
Marco Selig's avatar
Marco Selig committed
1215

1216
1217
1218
1219
1220
        elif(arg[0] == "syn"):
            lmax = self.para[0] - 1  # nlat-1
            if(self.dtype == np.dtype('float32')):
                x = gl.synfast_f(arg[1], nlat=self.para[0], nlon=self.para[
                                 1], lmax=lmax, mmax=lmax, alm=False)
Marco Selig's avatar
Marco Selig committed
1221
            else:
1222
1223
1224
                x = gl.synfast(arg[1], nlat=self.para[0], nlon=self.para[
                               1], lmax=lmax, mmax=lmax, alm=False)
            # weight if discrete
Marco Selig's avatar
Marco Selig committed
1225
            if(self.discrete):
1226
                x = self.calc_weight(x, power=0.5)
Marco Selig's avatar
Marco Selig committed
1227

1228
1229
1230
        elif(arg[0] == "uni"):
            x = random.uni(dtype=self.dtype, shape=self.get_dim(
                split=True), vmin=arg[1], vmax=arg[2])
Marco Selig's avatar
Marco Selig committed
1231
1232

        else:
1233
1234
            raise KeyError(about._errors.cstring(
                "ERROR: unsupported random key '" + str(arg[0]) + "'."))
Marco Selig's avatar
Marco Selig committed
1235
1236
1237

        return x

1238
    def check_codomain(self, codomain):
Marco Selig's avatar
Marco Selig committed
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
        """
            Checks whether a given codomain is compatible to the space or not.

            Parameters
            ----------
            codomain : nifty.space
                Space to be checked for compatibility.

            Returns
            -------
            check : bool
                Whether or not the given codomain is compatible to the space.

            Notes
            -----
            Compatible codomains are instances of :py:class:`gl_space` and
            :py:class:`lm_space`.
        """
1257
1258
        if codomain is None:
            return False
1259

1260
        if(not isinstance(codomain, space)):
Marco Selig's avatar
Marco Selig committed
1261
1262
            raise TypeError(about._errors.cstring("ERROR: invalid input."))

1263
1264
1265
        if self.datamodel is not codomain.datamodel:
            return False

1266
        if(self == codomain):
Marco Selig's avatar
Marco Selig committed
1267
1268
            return True

1269
1270
1271
1272
        if(isinstance(codomain, lm_space)):
            # nlon==2*lat-1                          lmax==nlat-1
            # lmax==mmax
            if(self.para[1] == 2 * self.para[0] - 1)and(codomain.para[0] == self.para[0] - 1)and(codomain.para[0] == codomain.para[1]):
Marco Selig's avatar
Marco Selig committed
1273
1274
1275
1276
1277
1278
                return True
            else:
                about.warnings.cprint("WARNING: unrecommended codomain.")

        return False

1279
    def get_codomain(self, **kwargs):
Marco Selig's avatar
Marco Selig committed
1280
1281
1282
1283
1284
1285
1286
1287
1288
        """
            Generates a compatible codomain to which transformations are
            reasonable, i.e.\  an instance of the :py:class:`lm_space` class.

            Returns
            -------
            codomain : nifty.lm_space
                A compatible codomain.
        """
1289
1290
1291
        if(self.dtype == np.dtype('float32')):
            # lmax,mmax = nlat-1,nlat-1
            return lm_space(self.para[0] - 1, mmax=self.para[0] - 1, dtype=np.complex64)
Marco Selig's avatar
Marco Selig committed
1292
        else:
1293
1294
            # lmax,mmax = nlat-1,nlat-1
            return lm_space(self.para[0] - 1, mmax=self.para[0] - 1, dtype=np.complex128)
Marco Selig's avatar
Marco Selig committed
1295

1296
    def get_meta_volume(self, split=False):
Marco Selig's avatar
Marco Selig committed
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
        """
            Calculates the meta volumes.

            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.

            Parameters
            ----------
            total : bool, *optional*
                Whether to return the total meta volume of the space or the
                individual ones of each field component (default: False).

            Returns
            -------
            mol : {numpy.ndarray, float}
                Meta volume of the field components or the complete space.

            Notes
            -----
            For Gauss-Legendre pixelizations, the meta volumes are the pixel
            sizes.
        """
1321
1322
        if not split:
            return np.float(4 * pi)
Marco Selig's avatar
Marco Selig committed
1323
        else:
1324
1325
            mol =