Merge pull request #3 from mselig/develop

version update.

README.rst

@@ -96,7 +96,7 @@ Requirements
 Download
 ........
 
-The latest release is tagged **v0.4.0** and is available as a source package
+The latest release is tagged **v0.4.2** and is available as a source package
 at `<https://github.com/mselig/nifty/tags>`_. The current version can be
 obtained by cloning the repository::
 

__init__.py

@@ -21,8 +21,9 @@
 
 from __future__ import division
 from nifty_core import *
-#from nifty_power import *
-#from nifty_cmaps import *
+from nifty_cmaps import *
+from nifty_power import *
+
 
 
 ##-----------------------------------------------------------------------------

nifty_core.py

@@ -484,7 +484,7 @@ class _about(object): ## nifty support class for global settings
 
         """
         ## version
-        self._version = "0.4.0"
+        self._version = "0.4.2"
 
         ## switches and notifications
         self._errors = notification(default=True,ccode=notification._code)
@@ -2178,6 +2178,9 @@ class rg_space(space):
         Notes
         -----
         Only even numbers of grid points per axis are supported.
+        The basis transformations between position `x` and Fourier mode `k`
+        rely on (inverse) fast Fourier transformations using the
+        :math:`exp(2 \pi i k^\dagger x)`-formulation.
 
         Attributes
         ----------
@@ -8367,10 +8370,6 @@ class operator(object):
 
             Parameters
             ----------
-            bare : bool, *optional*
-                Indicatese whether the diagonal entries are `bare` or not
-                (mandatory for the correct incorporation of volume weights)
-                (default: False)
             domain : space, *optional*
                 space wherein the probes live (default: self.domain)
             target : space, *optional*
@@ -8431,10 +8430,6 @@ class operator(object):
 
             Parameters
             ----------
-            bare : bool, *optional*
-                Indicatese whether the diagonal entries are `bare` or not
-                (mandatory for the correct incorporation of volume weights)
-                (default: False)
             domain : space, *optional*
                 space wherein the probes live (default: self.domain)
             target : space, *optional*
@@ -8605,7 +8600,7 @@ class diagonal_operator(operator):
             self.sym = True
         else:
             self.val = diag
-            about.infos.cprint("INFO: non-self-adjoint complex diagonal operator.")
+#            about.infos.cprint("INFO: non-self-adjoint complex diagonal operator.")
             self.sym = False
 
         ## check whether identity
@@ -8649,7 +8644,7 @@ class diagonal_operator(operator):
             self.sym = True
         else:
             self.val = newdiag
-            about.infos.cprint("INFO: non-self-adjoint complex diagonal operator.")
+#            about.infos.cprint("INFO: non-self-adjoint complex diagonal operator.")
             self.sym = False
 
         ## check whether identity
@@ -9652,7 +9647,7 @@ class projection_operator(operator):
             band = int(band)
             if(band>self.bands()-1)or(band<0):
                 raise TypeError(about._errors.cstring("ERROR: invalid band."))
-            Px = np.zeros(self.domain.dim(split=True),dtype=self.domain.datatype,order='C').flatten(order='C')
+            Px = np.zeros(self.domain.dim(split=False),dtype=self.domain.datatype,order='C')
             Px[self.ind[band]] += x.val.flatten(order='C')[self.ind[band]]
             Px = field(self.domain,val=Px,target=x.target)
             return Px
@@ -9664,7 +9659,7 @@ class projection_operator(operator):
                 bandsup = np.array(bandsup,dtype=np.int)
             if(np.any(bandsup>self.bands()-1))or(np.any(bandsup<0)):
                 raise ValueError(about._errors.cstring("ERROR: invalid input."))
-            Px = np.zeros(self.domain.dim(split=True),dtype=self.domain.datatype,order='C').flatten(order='C')
+            Px = np.zeros(self.domain.dim(split=False),dtype=self.domain.datatype,order='C')
             x_ = x.val.flatten(order='C')
             for bb in bandsup:
                 Px[self.ind[bb]] += x_[self.ind[bb]]

nifty_power.py

@@ -291,7 +291,7 @@ def _calc_inverse(tk,var,kindex,rho,b1,Amem): ## > computes the inverse Hessian
     ## inversion
     return np.linalg.inv(T2+np.diag(b2,k=0)),b2,Amem
 
-def infer_power(m,domain=None,Sk=None,D=None,pindex=None,pundex=None,kindex=None,rho=None,q=1E-42,alpha=1,perception=(1,0),smoothness=False,var=100,bare=True,**kwargs):
+def infer_power(m,domain=None,Sk=None,D=None,pindex=None,pundex=None,kindex=None,rho=None,q=1E-42,alpha=1,perception=(1,0),smoothness=True,var=10,bare=True,**kwargs):
     """
         Infers the power spectrum.
 
@@ -335,9 +335,9 @@ def infer_power(m,domain=None,Sk=None,D=None,pindex=None,pundex=None,kindex=None
             (default: (1,0)).
         smoothness : bool, *optional*
             Indicates whether the smoothness prior is used or not
-            (default: False).
+            (default: True).
         var : {scalar, list, array}, *optional*
-            Variance of the assumed spectral smoothness prior (default: 100).
+            Variance of the assumed spectral smoothness prior (default: 10).
         bare : bool, *optional*
             Indicates whether the power spectrum entries returned are "bare"
             or not (mandatory for the correct incorporation of volume weights)
@@ -495,7 +495,7 @@ def infer_power(m,domain=None,Sk=None,D=None,pindex=None,pundex=None,kindex=None
         trB2 = Sk.pseudo_tr(D,**kwargs) ## probing of the partial traces of D
         if(np.any(trB2<0)):
             about.warnings.cprint("WARNING: nonpositive value(s) in tr[DSk] reset to 0.")
-            trB2 = np.minimum(0,trB2)
+            trB2 = np.maximum(0,trB2)
     ## power spectrum
     numerator = 2*q+trB1+perception[0]*trB2 ## non-bare(!)
     denominator1 = rho+2*(alpha-1+perception[1])
@@ -503,30 +503,45 @@ def infer_power(m,domain=None,Sk=None,D=None,pindex=None,pundex=None,kindex=None
     if(smoothness):
         if(not domain.discrete):
             numerator = weight_power(domain,numerator,power=-1,pindex=pindex,pundex=pundex) ## bare(!)
-        pk = numerator/denominator1 ## bare(!)
 
         ## smoothness prior
-        tk = np.log(pk)
-        Amemory = None
-        var_ = var*1.1 # temporally increasing the variance
-        breakinfo = False
-        while(var_>=var): # slowly lowering the variance
-            absdelta = 1
-            while(absdelta>1E-3): # solving with fixed variance
-                ## solution of A delta = b1 - b2
-                Ainverse,denominator2,Amemory = _calc_inverse(tk,var_,kindex,rho,denominator1,Amemory)
-                delta = np.dot(Ainverse,numerator/pk-denominator2,out=None)
-                if(np.abs(delta).max()>absdelta): # increasing variance when speeding up
-                    var_ *= 1.1
-                absdelta = np.abs(delta).max()
-                tk += min(1,0.1/absdelta)*delta # adaptive step width
-                pk *= np.exp(min(1,0.1/absdelta)*delta) # adaptive step width
-            var_ /= 1.1 # lowering the variance when converged
-            if(var_<var):
-                if(breakinfo): # making sure there's one iteration with the correct variance
-                    break
-                var_ = var
-                breakinfo = True
+        divergence = 1
+        while(divergence):
+            pk = numerator/denominator1 ## bare(!)
+            tk = np.log(pk)
+            Amemory = None
+            var_ = var*1.1 # temporally increasing the variance
+            var_OLD = -1
+            breakinfo = False
+            while(var_>=var): # slowly lowering the variance
+                absdelta = 1
+                while(absdelta>1E-3): # solving with fixed variance
+                    ## solution of A delta = b1 - b2
+                    Ainverse,denominator2,Amemory = _calc_inverse(tk,var_,kindex,rho,denominator1,Amemory)
+                    delta = np.dot(Ainverse,numerator/pk-denominator2,out=None)
+                    if(np.abs(delta).max()>absdelta): # increasing variance when speeding up
+                        var_ *= 1.1
+                    absdelta = np.abs(delta).max()
+                    tk += min(1,0.1/absdelta)*delta # adaptive step width
+                    pk *= np.exp(min(1,0.1/absdelta)*delta) # adaptive step width
+                var_ /= 1.1 # lowering the variance when converged
+                if(var_<var):
+                    if(breakinfo): # making sure there's one iteration with the correct variance
+                        break
+                    var_ = var
+                    breakinfo = True
+                elif(var_==var_OLD):
+                    if(divergence==3):
+                        pot = int(np.log10(var_))
+                        var = int(1+var_*10**-pot)*10**pot
+                        about.warnings.cprint("WARNING: smoothness variance increased ( var = "+str(var)+" ).")
+                        break
+                    else:
+                        divergence += 1
+                else:
+                    var_OLD = var_
+            if(breakinfo):
+                break
 
         ## weight if ...
         if(not domain.discrete)and(not bare):

powerspectrum.py

@@ -216,7 +216,6 @@ def calc_ps_fast(field,axes,dgrid,zerocentered=False,fourier=False,pindex=None,k
             Fourier space have the same length (default=None).
 
     """
-
     ## field absolutes
     if(not fourier):
         foufield = np.fft.fftshift(np.fft.fftn(field))
@@ -243,8 +242,11 @@ def calc_ps_fast(field,axes,dgrid,zerocentered=False,fourier=False,pindex=None,k
                 ps[position] += fieldabs[ii]
                 rho[position] += 1
         else:
+            ## zerocenter pindex
+            if(np.any(zerocentered==False)):
+                pindex = np.fft.fftshift(pindex, axes=shiftaxes(zerocentered,st_to_zero_mode=True))
             ## power spectrum
-            ps = np.zeros(len(set(pindex.flatten())))
+            ps = np.zeros(np.max(pindex)+1)
             rho = np.zeros(ps.size)
             for ii in np.ndindex(pindex.shape):
                 ps[pindex[ii]] += fieldabs[ii]
@@ -263,6 +265,9 @@ def calc_ps_fast(field,axes,dgrid,zerocentered=False,fourier=False,pindex=None,k
             position = np.searchsorted(klength,kdict[ii])
             ps[position] += fieldabs[ii]
     else:
+        ## zerocenter pindex
+        if(np.any(zerocentered==False)):
+            pindex = np.fft.fftshift(pindex, axes=shiftaxes(zerocentered,st_to_zero_mode=True))
         ## power spectrum
         ps = np.zeros(rho.size)
         for ii in np.ndindex(pindex.shape):

setup.py

@@ -23,7 +23,7 @@ from distutils.core import setup
 import os
 
 setup(name="nifty",
-      version="0.4.0",
+      version="0.4.2",
       description="Numerical Information Field Theory",
       author="Marco Selig",
       author_email="mselig@mpa-garching.mpg.de",