submodule lm implemented.

lm/nifty_lm.py

+## 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/>.
+
+"""
+    ..                  __   ____   __
+    ..                /__/ /   _/ /  /_
+    ..      __ ___    __  /  /_  /   _/  __   __
+    ..    /   _   | /  / /   _/ /  /   /  / /  /
+    ..   /  / /  / /  / /  /   /  /_  /  /_/  /
+    ..  /__/ /__/ /__/ /__/    \___/  \___   /  lm
+    ..                               /______/
+
+    NIFTY submodule for grids on the two-sphere.
+
+"""
+from __future__ import division
+#from nifty import *
+import os
+import numpy as np
+import pylab as pl
+from matplotlib.colors import LogNorm as ln
+from matplotlib.ticker import LogFormatter as lf
+from nifty import pi,                                                        \
+                  about,                                                     \
+                  random,                                                    \
+                  space,                                                     \
+                  field
+#import libsharp_wrapper_gl as gl
+try:
+    import libsharp_wrapper_gl as gl
+except(ImportError):
+    about.infos.cprint("INFO: libsharp_wrapper_gl not available.")
+    _gl_available = False
+    about.warnings.cprint("WARNING: global setting 'about.lm2gl' corrected.")
+    about.lm2gl.off()
+else:
+    _gl_available = True
+#import healpy as hp
+try:
+    import healpy as hp
+except(ImportError):
+    about.infos.cprint("INFO: healpy not available.")
+    _hp_available = False
+else:
+    _hp_available = True
+
+
+##-----------------------------------------------------------------------------
+
+class lm_space(space):
+    """
+        ..       __
+        ..     /  /
+        ..    /  /    __ ____ ___
+        ..   /  /   /   _    _   |
+        ..  /  /_  /  / /  / /  /
+        ..  \___/ /__/ /__/ /__/  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`).
+        datatype : numpy.dtype, *optional*
+            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`.
+        datatype : numpy.dtype
+            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.
+    """
+    def __init__(self,lmax,mmax=None,datatype=None):
+        """
+            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`).
+            datatype : numpy.dtype, *optional*
+                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.
+
+        """
+        ## check imports
+        if(not _gl_available)and(not _hp_available):
+            raise ImportError(about._errors.cstring("ERROR: neither libsharp_wrapper_gl nor healpy available."))
+        ## check parameters
+        if(lmax<1):
+            raise ValueError(about._errors.cstring("ERROR: nonpositive number."))
+        if(lmax%2==0)and(lmax>2): ## exception lmax == 2 (nside == 1)
+            about.warnings.cprint("WARNING: unrecommended parameter ( lmax <> 2*n+1 ).")
+        if(mmax is None):
+            mmax = lmax
+        elif(mmax<1)or(mmax>lmax):
+            about.warnings.cprint("WARNING: parameter set to default.")
+            mmax = lmax
+        if(mmax!=lmax):
+            about.warnings.cprint("WARNING: unrecommended parameter ( mmax <> lmax ).")
+        self.para = np.array([lmax,mmax],dtype=np.int)
+
+        ## check data type
+        if(datatype is None):
+            datatype = np.complex128
+        elif(datatype not in [np.complex64,np.complex128]):
+            about.warnings.cprint("WARNING: data type set to default.")
+            datatype = np.complex128
+        self.datatype = datatype
+
+        self.discrete = True
+        self.vol = np.real(np.array([1],dtype=self.datatype))
+
+    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+    def lmax(self):
+        """
+            Returns the maximum quantum number :math:`\ell`.
+
+            Returns
+            -------
+            lmax : int
+                Maximum quantum number :math:`\ell`.
+        """
+        return self.para[0]
+
+    def mmax(self):
+        """
+            Returns the maximum quantum number :math:`m`.
+
+            Returns
+            -------
+            mmax : int
+                Maximum quantum number :math:`m`.
+
+        """
+        return self.para[1]
+
+    def 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 np.array([(self.para[0]+1)*(self.para[1]+1)-(self.para[1]+1)*self.para[1]//2],dtype=np.int)
+        else:
+            return (self.para[0]+1)*(self.para[1]+1)-(self.para[1]+1)*self.para[1]//2
+
+    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+    def dof(self):
+        """
+            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.
+        """
+        ## dof = 2*dim-(lmax+1) = (lmax+1)*(2*mmax+1)*(mmax+1)*mmax
+        return (self.para[0]+1)*(2*self.para[1]+1)-(self.para[1]+1)*self.para[1]
+
+    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+    def enforce_power(self,spec,**kwargs):
+        """
+            Provides a valid power spectrum array from a given object.
+
+            Parameters
+            ----------
+            spec : {float, list, numpy.ndarray, nifty.field, function}
+                Fiducial power spectrum from which a valid power spectrum is to
+                be calculated. Scalars are interpreted as constant power
+                spectra.
+
+            Returns
+            -------
+            spec : numpy.ndarray
+                Valid power spectrum.
+        """
+        if(isinstance(spec,field)):
+            spec = spec.val.astype(dtype=self.datatype)
+        elif(callable(spec)):
+            try:
+                spec = np.array(spec(np.arange(self.para[0]+1,dtype=self.vol.dtype)),dtype=self.datatype) ## prevent integer division
+            except:
+                raise TypeError(about._errors.cstring("ERROR: invalid power spectra function.")) ## exception in ``spec(kindex)``
+        elif(np.isscalar(spec)):
+            spec = np.array([spec],dtype=self.datatype)
+        else:
+            spec = np.array(spec,dtype=self.datatype)
+
+        ## drop imaginary part
+        spec = np.real(spec)
+        ## check finiteness
+        if(not np.all(np.isfinite(spec))):
+            about.warnings.cprint("WARNING: infinite value(s).")
+        ## check positivity (excluding null)
+        if(np.any(spec<0)):
+            raise ValueError(about._errors.cstring("ERROR: nonpositive value(s)."))
+        elif(np.any(spec==0)):
+            about.warnings.cprint("WARNING: nonpositive value(s).")
+
+        size = self.para[0]+1 ## lmax+1
+        ## extend
+        if(np.size(spec)==1):
+            spec = spec*np.ones(size,dtype=spec.dtype,order='C')
+        ## size check
+        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
+
+    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+    def _getlm(self): ## > compute all (l,m)
+        index = np.arange(self.dim(split=False))
+        n = 2*self.para[0]+1
+        m = np.ceil((n-np.sqrt(n**2-8*(index-self.para[0])))/2).astype(np.int)
+        l = index-self.para[0]*m+m*(m-1)//2
+        return l,m
+
+    def set_power_indices(self,**kwargs):
+        """
+            Sets the (un)indexing objects for spectral indexing internally.
+
+            Parameters
+            ----------
+            None
+
+            Returns
+            -------
+            None
+
+            See Also
+            --------
+            get_power_indices
+
+        """
+        ## check storage
+        if(not hasattr(self,"power_indices")):
+            ## power indices
+#            about.infos.cflush("INFO: setting power indices ...")
+            kindex = np.arange(self.para[0]+1,dtype=np.int)
+            rho = 2*kindex+1
+            if(_hp_available): ## default
+                pindex = hp.Alm.getlm(self.para[0],i=None)[0] ## l of (l,m)
+            else:
+                pindex = self._getlm()[0] ## l of (l,m)
+            pundex = np.unique(pindex,return_index=True,return_inverse=False)[1]
+            ## storage
+            self.power_indices = {"kindex":kindex,"pindex":pindex,"pundex":pundex,"rho":rho} ## alphabetical
+#            about.infos.cprint(" done.")
+
+    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+    def enforce_values(self,x,extend=True):
+        """
+            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).
+        """
+        if(isinstance(x,field)):
+            if(self==x.domain):
+                if(self.datatype is not x.domain.datatype):
+                    raise TypeError(about._errors.cstring("ERROR: inequal data types ( '"+str(np.result_type(self.datatype))+"' <> '"+str(np.result_type(x.domain.datatype))+"' )."))
+                else:
+                    x = np.copy(x.val)
+            else:
+                raise ValueError(about._errors.cstring("ERROR: inequal domains."))
+        else:
+            if(np.size(x)==1):
+                if(extend):
+                    x = self.datatype(x)*np.ones(self.dim(split=True),dtype=self.datatype,order='C')
+                else:
+                    if(np.isscalar(x)):
+                        x = np.array([x],dtype=self.datatype)
+                    else:
+                        x = np.array(x,dtype=self.datatype)
+            else:
+                x = self.enforce_shape(np.array(x,dtype=self.datatype))
+
+        if(np.size(x)!=1)and(np.any(x.imag[:self.para[0]+1]!=0)):
+            about.warnings.cprint("WARNING: forbidden values reset.")
+            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
+
+        ## check finiteness
+        if(not np.all(np.isfinite(x))):
+            about.warnings.cprint("WARNING: infinite value(s).")
+
+        return x
+
+    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+    def get_random_values(self,**kwargs):
+        """
+            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).
+        """
+        arg = random.arguments(self,**kwargs)
+
+        if(arg is None):
+            return np.zeros(self.dim(split=True),dtype=self.datatype,order='C')
+
+        elif(arg[0]=="pm1"):
+            x = random.pm1(datatype=self.datatype,shape=self.dim(split=True))
+
+        elif(arg[0]=="gau"):
+            x = random.gau(datatype=self.datatype,shape=self.dim(split=True),mean=None,dev=arg[2],var=arg[3])
+
+        elif(arg[0]=="syn"):
+            lmax = self.para[0] ## lmax
+            if(self.datatype==np.complex64):
+                if(_gl_available): ## default
+                    x = gl.synalm_f(arg[1],lmax=lmax,mmax=lmax)
+                else:
+                    x = hp.synalm(arg[1].astype(np.complex128),lmax=lmax,mmax=lmax).astype(np.complex64)
+            else:
+                if(_hp_available): ## default
+                    x = hp.synalm(arg[1],lmax=lmax,mmax=lmax)
+                else:
+                    x = gl.synalm(arg[1],lmax=lmax,mmax=lmax)
+
+            return x
+
+        elif(arg[0]=="uni"):
+            x = random.uni(datatype=self.datatype,shape=self.dim(split=True),vmin=arg[1],vmax=arg[2])
+
+        else:
+            raise KeyError(about._errors.cstring("ERROR: unsupported random key '"+str(arg[0])+"'."))
+
+        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
+
+        return x
+
+    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+    def check_codomain(self,codomain):
+        """
+            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`.
+        """
+        if(not isinstance(codomain,space)):
+            raise TypeError(about._errors.cstring("ERROR: invalid input."))
+
+        if(self==codomain):
+            return True
+
+        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):
+                return True
+            else:
+                about.warnings.cprint("WARNING: unrecommended codomain.")
+
+        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]):
+                return True
+            else:
+                about.warnings.cprint("WARNING: unrecommended codomain.")
+
+        return False
+
+    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+    def get_codomain(self,coname=None,**kwargs):
+        """
+            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>`_
+
+        """
+        if(coname=="gl")or(coname is None)and(about.lm2gl.status): ## order matters
+            if(self.datatype==np.complex64):
+                return gl_space(self.para[0]+1,nlon=2*self.para[0]+1,datatype=np.float32) ## nlat,nlon = lmax+1,2*lmax+1
+            else:
+                return gl_space(self.para[0]+1,nlon=2*self.para[0]+1,datatype=np.float64) ## nlat,nlon = lmax+1,2*lmax+1
+
+        elif(coname=="hp")or(coname is None)and(not about.lm2gl.status):
+            return hp_space((self.para[0]+1)//3) ## nside = (lmax+1)/3
+
+        else:
+            raise ValueError(about._errors.cstring("ERROR: unsupported or incompatible space '"+coname+"'."))
+
+    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+    def get_meta_volume(self,total=False):
+        """
+            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
+            -----
+            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`.
+        """
+        if(total):
+            return self.dof()
+        else:
+            mol = np.ones(self.dim(split=True),dtype=self.vol.dtype,order='C')
+            mol[self.para[0]+1:] = 2 ## redundant in (l,m) and (l,-m)
+            return mol
+
+    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+    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)
+        return dot
+
+    def calc_dot(self,x,y):
+        """
+            Computes the discrete inner product of two given arrays of field
+            values.
+
+            Parameters
+            ----------
+            x : numpy.ndarray
+                First array