From 1e805c57ae65f0b8691663fad109cbba8d59bb92 Mon Sep 17 00:00:00 2001
From: Marco Selig <mselig@ncg-02.MPA-Garching.MPG.DE>
Date: Thu, 11 Sep 2014 16:10:46 +0200
Subject: [PATCH] intermediate update; healpy update hotfix.
---
__init__.py | 1 +
nifty_core.py | 14 +--
nifty_power_conversion.py | 253 ++++++++++++++++++++++++++++++++++++++
nifty_tools.py | 18 +--
4 files changed, 271 insertions(+), 15 deletions(-)
create mode 100644 nifty_power_conversion.py
diff --git a/__init__.py b/__init__.py
index 5dd278c84..564ff2ece 100644
--- a/__init__.py
+++ b/__init__.py
@@ -26,6 +26,7 @@ from nifty_power import *
from nifty_tools import *
from nifty_explicit import *
+from nifty_power_conversion import *
##-----------------------------------------------------------------------------
diff --git a/nifty_core.py b/nifty_core.py
index 72ab4ea6b..338f35c54 100644
--- a/nifty_core.py
+++ b/nifty_core.py
@@ -514,7 +514,7 @@ class _about(object): ## nifty support class for global settings
"""
## version
- self._version = "0.8.1"
+ self._version = "0.8.4"
## switches and notifications
self._errors = notification(default=True,ccode=notification._code)
@@ -5137,7 +5137,7 @@ class hp_space(space):
if(self.discrete):
x = self.calc_weight(x,power=-0.5)
## transform
- Tx = hp.map2alm(x.astype(np.float64),lmax=codomain.para[0],mmax=codomain.para[1],iter=kwargs.get("iter",self.niter),pol=True,use_weights=False,regression=True,datapath=None)
+ Tx = hp.map2alm(x.astype(np.float64),lmax=codomain.para[0],mmax=codomain.para[1],iter=kwargs.get("iter",self.niter),pol=True,use_weights=False,datapath=None)
else:
raise ValueError(about._errors.cstring("ERROR: unsupported transformation."))
@@ -5181,7 +5181,7 @@ class hp_space(space):
elif(sigma<0):
raise ValueError(about._errors.cstring("ERROR: invalid sigma."))
## smooth
- return hp.smoothing(x,fwhm=0.0,sigma=sigma,invert=False,pol=True,iter=kwargs.get("iter",self.niter),lmax=3*self.para[0]-1,mmax=3*self.para[0]-1,use_weights=False,regression=True,datapath=None)
+ return hp.smoothing(x,fwhm=0.0,sigma=sigma,invert=False,pol=True,iter=kwargs.get("iter",self.niter),lmax=3*self.para[0]-1,mmax=3*self.para[0]-1,use_weights=False,datapath=None)
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -5211,7 +5211,7 @@ class hp_space(space):
if(self.discrete):
x = self.calc_weight(x,power=-0.5)
## power spectrum
- return hp.anafast(x,map2=None,nspec=None,lmax=3*self.para[0]-1,mmax=3*self.para[0]-1,iter=kwargs.get("iter",self.niter),alm=False,pol=True,use_weights=False,regression=True,datapath=None)
+ return hp.anafast(x,map2=None,nspec=None,lmax=3*self.para[0]-1,mmax=3*self.para[0]-1,iter=kwargs.get("iter",self.niter),alm=False,pol=True,use_weights=False,datapath=None)
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -8755,7 +8755,7 @@ class diagonal_operator(operator):
if(np.any(self.val==0)):
if(pseudo):
x_ = field(self.domain,val=None,target=x.target)
- x_.val = x.val*np.where(self.val==0,0,1/self.val) ## bypasses self.domain.enforce_values
+ x_.val = np.ma.filled(x.val/np.ma.masked_where(self.val==0,self.val,copy=False),fill_value=0) ## bypasses self.domain.enforce_values
return x_
else:
raise AttributeError(about._errors.cstring("ERROR: singular operator."))
@@ -8768,7 +8768,7 @@ class diagonal_operator(operator):
if(np.any(self.val==0)):
if(pseudo):
x_ = field(self.domain,val=None,target=x.target)
- x_.val = x.val*np.where(self.val==0,0,1/np.conjugate(self.val)) ## bypasses self.domain.enforce_values
+ x_.val = np.ma.filled(x.val/np.ma.masked_where(self.val==0,np.conjugate(self.val),copy=False),fill_value=0) ## bypasses self.domain.enforce_values
return x_
else:
raise AttributeError(about._errors.cstring("ERROR: singular operator."))
@@ -8781,7 +8781,7 @@ class diagonal_operator(operator):
if(np.any(self.val==0)):
if(pseudo):
x_ = field(self.domain,val=None,target=x.target)
- x_.val = x.val*np.where(self.val==0,0,np.conjugate(1/self.val)) ## bypasses self.domain.enforce_values
+ x_.val = np.ma.filled(x.val/np.conjugate(np.ma.masked_where(self.val==0,self.val,copy=False)),fill_value=0) ## bypasses self.domain.enforce_values
return x_
else:
raise AttributeError(about._errors.cstring("ERROR: singular operator."))
diff --git a/nifty_power_conversion.py b/nifty_power_conversion.py
new file mode 100644
index 000000000..999e18f32
--- /dev/null
+++ b/nifty_power_conversion.py
@@ -0,0 +1,253 @@
+from nifty import *
+
+def power_backward_conversion_rg(k_space,p,mean=None,bare=True):
+ """
+ This function is designed to convert a theoretical/statistical power
+ spectrum of a log-normal field to the theoretical power spectrum of
+ the underlying Gaussian field.
+ The function only works for power spectra defined for rg_spaces
+
+ Parameters
+ ----------
+ k_space : nifty.rg_space,
+ a regular grid space with the attribute `Fourier = True`
+ p : np.array,
+ the power spectrum of the log-normal field.
+ Needs to have the same number of entries as
+ `k_space.get_power_indices()[0]`
+ mean : float, *optional*
+ specifies the mean of the log-normal field. If `mean` is not
+ specified the function will use the monopole of the power spectrum.
+ If it is specified the function will NOT use the monopole of the
+ spectrum (default: None).
+ WARNING: a mean that is too low can violate positive definiteness
+ of the log-normal field. In this case the function produces an
+ error.
+ bare : bool, *optional*
+ whether `p` is the bare power spectrum or not (default: True).
+
+ Returns
+ -------
+ mean : float,
+ the recovered mean of the underlying Gaussian distribution.
+ p1 : np.array,
+ the power spectrum of the underlying Gaussian field, where the
+ monopole has been set to zero. Eventual monopole power has been
+ shifted to the mean.
+
+ References
+ ----------
+ .. [#] M. Greiner and T.A. Ensslin, "Log-transforming the matter power spectrum";
+ `arXiv:1312.1354 <http://arxiv.org/abs/1312.1354>`_
+ """
+ pindex = k_space.get_power_indices()[2]
+ V = k_space.vol.prod()**(-1)
+
+ mono_ind = np.where(pindex==0)
+
+ spec = power_operator(k_space,spec=p,bare=bare).get_power(bare=False)
+
+ if(mean is None):
+ mean = 0.
+ else:
+ spec[0] = 0.
+
+ pf = field(k_space,val=spec[pindex]).transform()+mean**2
+
+ if(np.any(pf.val<0.)):
+ raise ValueError(about._errors.cstring("ERROR: spectrum or mean incompatible with positive definiteness.\n Try increasing the mean."))
+ return None
+
+ p1 = sqrt(log(pf).power())
+
+ p1[0] = (log(pf)).transform()[mono_ind][0]
+
+ p2 = 0.5*V*log(k_space.calc_weight(spec[pindex],1).sum()+mean**2)
+
+ logmean = 1/V * (p1[0]-p2)
+
+ p1[0] = 0.
+
+ if(np.any(p1<0.)):
+ raise ValueError(about._errors.cstring("ERROR: spectrum or mean incompatible with positive definiteness.\n Try increasing the mean."))
+ return None
+
+ if(bare==True):
+ return logmean.real,power_operator(k_space,spec=p1,bare=False).get_power(bare=True).real
+ else:
+ return logmean.real,p1.real
+
+def power_forward_conversion_rg(k_space,p,mean=0,bare=True):
+ """
+ This function is designed to convert a theoretical/statistical power
+ spectrum of a Gaussian field to the theoretical power spectrum of
+ the exponentiated field.
+ The function only works for power spectra defined for rg_spaces
+
+ Parameters
+ ----------
+ k_space : nifty.rg_space,
+ a regular grid space with the attribute `Fourier = True`
+ p : np.array,
+ the power spectrum of the Gaussian field.
+ Needs to have the same number of entries as
+ `k_space.get_power_indices()[0]`
+ mean : float, *optional*
+ specifies the mean of the Gaussian field (default: 0).
+ bare : bool, *optional*
+ whether `p` is the bare power spectrum or not (default: True).
+
+ Returns
+ -------
+ p1 : np.array,
+ the power spectrum of the exponentiated Gaussian field.
+
+ References
+ ----------
+ .. [#] M. Greiner and T.A. Ensslin, "Log-transforming the matter power spectrum";
+ `arXiv:1312.1354 <http://arxiv.org/abs/1312.1354>`_
+ """
+
+ pindex = k_space.get_power_indices()[2]
+
+ spec = power_operator(k_space,spec=p,bare=bare).get_power(bare=False)
+
+ S_x = field(k_space,val=spec[pindex]).transform()
+
+ S_0 = k_space.calc_weight(spec[pindex],1).sum()
+
+ pf = exp(S_x+S_0+2*mean)
+
+ p1 = sqrt(pf.power())
+
+ if(bare==True):
+ return power_operator(k_space,spec=p1,bare=False).get_power(bare=True).real
+ else:
+ return p1.real
+
+
+
+def power_backward_conversion_lm(k_space,p,mean=None):
+ """
+ This function is designed to convert a theoretical/statistical power
+ spectrum of a log-normal field to the theoretical power spectrum of
+ the underlying Gaussian field.
+ The function only works for power spectra defined for lm_spaces
+
+ Parameters
+ ----------
+ k_space : nifty.rg_space,
+ a regular grid space with the attribute `Fourier = True`
+ p : np.array,
+ the power spectrum of the log-normal field.
+ Needs to have the same number of entries as
+ `k_space.get_power_indices()[0]`
+ mean : float, *optional*
+ specifies the mean of the log-normal field. If `mean` is not
+ specified the function will use the monopole of the power spectrum.
+ If it is specified the function will NOT use the monopole of the
+ spectrum. (default: None)
+ WARNING: a mean that is too low can violate positive definiteness
+ of the log-normal field. In this case the function produces an
+ error.
+
+ Returns
+ -------
+ mean : float,
+ the recovered mean of the underlying Gaussian distribution.
+ p1 : np.array,
+ the power spectrum of the underlying Gaussian field, where the
+ monopole has been set to zero. Eventual monopole power has been
+ shifted to the mean.
+
+ References
+ ----------
+ .. [#] M. Greiner and T.A. Ensslin, "Log-transforming the matter power spectrum";
+ `arXiv:1312.1354 <http://arxiv.org/abs/1312.1354>`_
+ """
+
+ p = np.copy(p)
+ if(mean is not None):
+ p[0] = 4*pi*mean**2
+
+ klen = k_space.get_power_indices()[0]
+ C_0_Omega = field(k_space,val=0)
+ C_0_Omega.val[:len(klen)] = p*sqrt(2*klen+1)/sqrt(4*pi)
+ C_0_Omega = C_0_Omega.transform()
+
+ if(np.any(C_0_Omega.val<0.)):
+ raise ValueError(about._errors.cstring("ERROR: spectrum or mean incompatible with positive definiteness.\n Try increasing the mean."))
+ return None
+
+ lC = log(C_0_Omega)
+
+ Z = lC.transform()
+
+ spec = Z.val[:len(klen)]
+
+ mean = (spec[0]-0.5*sqrt(4*pi)*log((p*(2*klen+1)/(4*pi)).sum()))/sqrt(4*pi)
+
+ spec[0] = 0.
+
+ spec = spec*sqrt(4*pi)/sqrt(2*klen+1)
+
+ spec = np.real(spec)
+
+ if(np.any(spec<0.)):
+ spec = spec*(spec>0.)
+ about.warnings.cprint("WARNING: negative modes set to zero.")
+
+ return mean.real,spec
+
+
+def power_forward_conversion_lm(k_space,p,mean=0):
+ """
+ This function is designed to convert a theoretical/statistical power
+ spectrum of a Gaussian field to the theoretical power spectrum of
+ the exponentiated field.
+ The function only works for power spectra defined for lm_spaces
+
+ Parameters
+ ----------
+ k_space : nifty.rg_space,
+ a regular grid space with the attribute `Fourier = True`
+ p : np.array,
+ the power spectrum of the Gaussian field.
+ Needs to have the same number of entries as
+ `k_space.get_power_indices()[0]`
+ m : float, *optional*
+ specifies the mean of the Gaussian field (default: 0).
+
+ Returns
+ -------
+ p1 : np.array,
+ the power spectrum of the exponentiated Gaussian field.
+
+ References
+ ----------
+ .. [#] M. Greiner and T.A. Ensslin, "Log-transforming the matter power spectrum";
+ `arXiv:1312.1354 <http://arxiv.org/abs/1312.1354>`_
+ """
+ m = mean
+ klen = k_space.get_power_indices()[0]
+ C_0_Omega = field(k_space,val=0)
+ C_0_Omega.val[:len(klen)] = p*sqrt(2*klen+1)/sqrt(4*pi)
+ C_0_Omega = C_0_Omega.transform()
+
+ C_0_0 = (p*(2*klen+1)/(4*pi)).sum()
+
+ exC = exp(C_0_Omega+C_0_0+2*m)
+
+ Z = exC.transform()
+
+ spec = Z.val[:len(klen)]
+
+ spec = spec*sqrt(4*pi)/sqrt(2*klen+1)
+
+ spec = np.real(spec)
+
+ if(np.any(spec<0.)):
+ spec = spec*(spec>0.)
+ about.warnings.cprint("WARNING: negative modes set to zero.")
+
+ return spec
\ No newline at end of file
diff --git a/nifty_tools.py b/nifty_tools.py
index db13db114..a91d83881 100644
--- a/nifty_tools.py
+++ b/nifty_tools.py
@@ -155,8 +155,8 @@ class invertible_operator(operator):
x : field
Valid input field.
force : bool
- Indicates wheter to return a field instead of ``None`` is
- forced incase the conjugate gradient fails.
+ Indicates wheter to return a field instead of ``None`` in case
+ the conjugate gradient fails.
Returns
-------
@@ -213,8 +213,8 @@ class invertible_operator(operator):
x : field
Valid input field.
force : bool
- Indicates wheter to return a field instead of ``None`` is
- forced incase the conjugate gradient fails.
+ Indicates wheter to return a field instead of ``None`` in case
+ the conjugate gradient fails.
Returns
-------
@@ -430,10 +430,10 @@ class propagator_operator(operator):
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def _inverse_multiply_1(self,x,**kwargs): ## > applies A1 + A2 in self.codomain
- return self._A1(x)+self._A2(x.transform(self.domain)).transform(self.codomain)
+ return self._A1(x,pseudo=True)+self._A2(x.transform(self.domain)).transform(self.codomain)
def _inverse_multiply_2(self,x,**kwargs): ## > applies A1 + A2 in self.domain
- return self._A1(x.transform(self.codomain)).transform(self.domain)+self._A2(x)
+ return self._A1(x.transform(self.codomain),pseudo=True).transform(self.domain)+self._A2(x)
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -475,8 +475,8 @@ class propagator_operator(operator):
Scalars are interpreted as constant arrays, and an array will
be interpreted as a field on the domain of the operator.
force : bool
- Indicates wheter to return a field instead of ``None`` is
- forced incase the conjugate gradient fails.
+ Indicates wheter to return a field instead of ``None`` in case
+ the conjugate gradient fails.
Returns
-------
@@ -836,6 +836,8 @@ class conjugate_gradient(object):
s = self.W(r)
gamma_ = gamma
gamma = r.dot(s)
+ if(np.signbit(np.real(gamma))):
+ about.warnings.cprint("WARNING: positive definiteness of W violated.")
beta = max(0,gamma/gamma_) ## positive definite
d = s+beta*d ## conjugated gradient
--
GitLab