Merge pull request #13 from mselig/develop

Develop to Master

.gitignore

@@ -12,8 +12,21 @@
 .spyderproject
 build
 
+rg/*
+!rg/__init__.py
+!rg/nifty_rg.py
+!rg/nifty_power_conversion_rg.py
+!rg/gfft_rg.py
+!rg/powerspectrum.py
+
+lm/*
+!lm/__init__.py
+!lm/nifty_lm.py
+!lm/nifty_power_conversion_lm.py
+
 demos/*
 !demos/__init__.py
+!demos/demos_core.py
 !demos/demo_faraday.py
 !demos/demo_faraday_map.npy
 !demos/demo_excaliwir.py

README.rst

@@ -87,19 +87,21 @@ Requirements
         (standard library)
 
 *   `GFFT <https://github.com/mrbell/gfft>`_ (v0.1.0) - Generalized Fast
-    Fourier Transformations for Python
+    Fourier Transformations for Python - **optional**
 
 *   `HEALPy <https://github.com/healpy/healpy>`_ (v1.4 without openmp) - A
-    Python wrapper for `HEALPix <http://sourceforge.net/projects/healpix/>`_
+    Python wrapper for `HEALPix <http://sourceforge.net/projects/healpix/>`_ -
+    **optional**
 
 *   `libsharp-wrapper <https://github.com/mselig/libsharp-wrapper>`_ (v0.1.2
     without openmp) - A Python wrapper for the
-    `libsharp <http://sourceforge.net/projects/libsharp/>`_ library
+    `libsharp <http://sourceforge.net/projects/libsharp/>`_ library -
+    **optional**
 
 Download
 ........
 
-The latest release is tagged **v0.9.0** and is available as a source package
+The latest release is tagged **v1.0.6** and is available as a source package
 at `<https://github.com/mselig/nifty/tags>`_. The current version can be
 obtained by cloning the repository::
 
@@ -109,14 +111,25 @@ obtained by cloning the repository::
 Installation
 ............
 
-NIFTY is installed using Distutils by running the following command::
+*   NIFTY can be installed using `PyPI <https://pypi.python.org/pypi>`_ and
+    **pip** by running the following command::
 
-    python setup.py install
+        pip install nifty
 
-Alternatively, a private or user specific installation can be done by::
+    Alternatively, a private or user specific installation can be done by::
 
-    python setup.py install --user
-    python setup.py install --install-lib=/SOMEWHERE
+        pip install --user nifty
+
+
+*   NIFTY can be installed using **Distutils** by running the following
+    command::
+
+        python setup.py install
+
+    Alternatively, a private or user specific installation can be done by::
+
+        python setup.py install --user
+        python setup.py install --install-lib=/SOMEWHERE
 
 Acknowledgement
 ---------------

__init__.py

@@ -20,43 +20,24 @@
 ## along with this program. If not, see <http://www.gnu.org/licenses/>.
 
 from __future__ import division
-from nifty_core import *
+from nifty_core import * ## imports `about`
 from nifty_cmaps import *
 from nifty_power import *
 from nifty_tools import *
 from nifty_explicit import *
 
-#from nifty_power_conversion import *
+## optional submodule `rg`
+try:
+    from rg import *
+except(ImportError):
+    pass
 
+## optional submodule `lm`
+try:
+    from lm import *
+except(ImportError):
+    pass
 
-##-----------------------------------------------------------------------------
+from demos import *
+from pickling import *
 
-import copy_reg as cr
-from types import MethodType as mt
-
-
-def _pickle_method(method):
-    fct_name = method.im_func.__name__
-    obj = method.im_self
-    cl = method.im_class
-    ## handle mangled function name
-    if(fct_name.startswith("__"))and(not fct_name.endswith("__")):
-        cl_name = cl.__name__.lstrip("_")
-        fct_name = "_" + cl_name + fct_name
-    return _unpickle_method, (fct_name, obj, cl)
-
-
-def _unpickle_method(fct_name, obj, cl):
-    for oo in cl.__mro__:
-        try:
-            fct = oo.__dict__[fct_name]
-        except(KeyError):
-            pass
-        else:
-            break
-    return fct.__get__(obj, cl)
-
-## enable instance methods pickling
-cr.pickle(mt, _pickle_method, _unpickle_method)
-
-##-----------------------------------------------------------------------------
\ No newline at end of file

demos/__init__.py

@@ -19,5 +19,5 @@
 ## You should have received a copy of the GNU General Public License
 ## along with this program. If not, see <http://www.gnu.org/licenses/>.
 
-pass
+from demos_core import *
 

demos/demo_faraday.py

@@ -48,7 +48,7 @@ about.infos.off()
 ##-----------------------------------------------------------------------------
 
 # (global) Faraday map
-m = field(hp_space(128), val=np.load("demo_faraday_map.npy"))
+m = field(hp_space(128), val=np.load(os.path.join(get_demo_dir(),"demo_faraday_map.npy")))
 
 ##-----------------------------------------------------------------------------
 

demos/demos_core.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/>.
+
+"""
+    NIFTY provides a number of demonstrations that you can run either using
+    ``execfile`` provided with (absolute or relative) file names or using
+    ``run -m`` providing the module name. You can retrieve the directory of
+    the NIFTY demos calling :py:func:`get_demo_dir`.
+
+"""
+import os
+import nifty as nt
+
+def get_demo_dir():
+    """
+        Returns the path of the NIFTY demos directory.
+
+    """
+    return os.path.split(nt.demos.__file__)[0]
+

lm/__init__.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/>.
+
+from __future__ import division
+from nifty import about
+#from nifty_lm import *
+try:
+    import libsharp_wrapper_gl as gl
+except(ImportError):
+    try:
+        import healpy as hp
+    except(ImportError):
+        about.infos.cprint("INFO: neither libsharp_wrapper_gl nor healpy available.")
+        pass ## import nothing
+    else:
+        from nifty_lm import lm_space,hp_space ## import lm & hp
+        ## TODO: change about
+else:
+    try:
+        import healpy as hp
+    except(ImportError):
+        from nifty_lm import lm_space,gl_space ## import lm & gl
+    else:
+        from nifty_lm import lm_space,gl_space,hp_space ## import all 3
+from nifty_power_conversion_lm import *
+

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)