Skip to content

GitLab

Commit 57660e7e authored by Marco Selig's avatar Marco Selig
Browse files

submodule lm implemented.

parent b864b70d
Hide whitespace changes
Inline Side-by-side
.gitignore
View file @ 57660e7e
......@@ -22,8 +22,6 @@ rg/*
lm/*
!lm/__init__.py
!lm/nifty_lm.py
!lm/nifty_gl.py
!lm/nifty_hp.py
!lm/nifty_power_conversion_lm.py
demos/*
......
lm/__init__.py
View file @ 57660e7e
......@@ -20,6 +20,25 @@
## along with this program. If not, see <http://www.gnu.org/licenses/>.
from __future__ import division
#from nifty_rg import *
#from nifty_lm import *
try:
import libsharp_wrapper_gl as gl
except(ImportError):
try:
import healpy as hp
except(ImportError):
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 *
nifty_core.py
View file @ 57660e7e
......@@ -140,26 +140,20 @@
:py:func:`interpolate_power`
"""
## standard libraries
from __future__ import division
import os
#import sys
from sys import stdout as so
import numpy as np
import pylab as pl
from matplotlib.colors import LogNorm as ln
from matplotlib.ticker import LogFormatter as lf
from multiprocessing import Pool as mp
from multiprocessing import Value as mv
from multiprocessing import Array as ma
## third party libraries
import healpy as hp
import libsharp_wrapper_gl as gl
pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
__version__ = "1.0.4"
__version__ = "1.0.2"
pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
##-----------------------------------------------------------------------------
......@@ -2102,2080 +2096,2080 @@ class point_space(space):
##-----------------------------------------------------------------------------
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.
"""
## 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 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
pindex = hp.Alm.getlm(self.para[0],i=None)[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):
x = gl.synalm_f(arg[1],lmax=lmax,mmax=lmax)
else:
#x = gl.synalm(arg[1],lmax=lmax,mmax=lmax)
x = hp.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 calc_dot(self,x,y):
"""
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.
"""
x = self.enforce_shape(np.array(x,dtype=self.datatype))
y = self.enforce_shape(np.array(y,dtype=self.datatype))
## inner product
if(self.datatype==np.complex64):
return gl.dotlm_f(x,y,lmax=self.para[0],mmax=self.para[1])
else:
return gl.dotlm(x,y,lmax=self.para[0],mmax=self.para[1])
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def calc_transform(self,x,codomain=None,**kwargs):
"""
Computes the transform of a given array of field values.
Parameters
----------
x : numpy.ndarray