Merge pull request #10 from mselig/develop

intermediate update.

__init__.py

@@ -26,6 +26,7 @@ from nifty_power import *
 from nifty_tools import *
 from nifty_explicit import *
 
+from nifty_power_conversion import *
 
 
 ##-----------------------------------------------------------------------------

nifty_core.py

@@ -514,7 +514,7 @@ class _about(object): ## nifty support class for global settings
 
         """
         ## version
-        self._version = "0.8.0"
+        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."))

nifty_power_conversion.py

+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

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