Merge pull request #1 from mselig/develop

Develop > Master

README.rst

@@ -96,7 +96,7 @@ Requirements
 Download
 ........
 
-The latest release is tagged **v0.3.0** and is available as a source package
+The latest release is tagged **v0.4.0** and is available as a source package
 at `<https://github.com/mselig/nifty/tags>`_. The current version can be
 obtained by cloning the repository::
 
@@ -128,8 +128,9 @@ References
 ..........
 
 .. [1] Selig et al., "NIFTY - Numerical Information Field Theory - a
-    versatile Python library for signal inference", submitted to A&A, 2013;
-    `arXiv:1301.4499 <http://www.arxiv.org/abs/1301.4499>`_
+    versatile Python library for signal inference",
+    `A&A, vol. 554, id. A26 <http://dx.doi.org/10.1051/0004-6361/201321236>`_,
+    2013; `arXiv:1301.4499 <http://www.arxiv.org/abs/1301.4499>`_
 
 Release Notes
 -------------

demos/demo_excaliwir.py

@@ -52,7 +52,7 @@ g = gl_space(48)
 z = s_space = k = k_space = p = d_space = None
 
 ## power spectrum (and more)
-power = powerindex = powerundex = kindex = rho = None
+power = kindex = rho = powerindex = powerundex = None
 
 ## operators
 S = Sk = R = N = Nj = D = None
@@ -111,7 +111,7 @@ def setup(space,s2n=3,nvar=None):
             the noise variance, `nvar` will be calculated according to
             `s2n` if not specified (default: None)
     """
-    global z,s_space,k,k_space,p,d_space,power,powerindex,powerundex,kindex,rho,S,Sk,R,N,Nj,D,s,n,d,j,m
+    global z,s_space,k,k_space,p,d_space,power,kindex,rho,powerindex,powerundex,S,Sk,R,N,Nj,D,s,n,d,j,m
 
     ## signal space
     z = s_space = space
@@ -119,17 +119,15 @@ def setup(space,s2n=3,nvar=None):
     ## conjugate space
     k = k_space = s_space.get_codomain()
     ## the power indices are calculated once and saved
-    powerindex = k_space.get_power_index()
-    powerundex = k_space.get_power_undex()
-    kindex,rho = k_space.get_power_index(irreducible=True)
+    kindex,rho,powerindex,powerundex = k_space.get_power_indices()
 
     ## power spectrum
-    power = [42/(kk+1)**3 for kk in kindex]
+    power = (lambda kk: 42/(kk+1)**3)
 
     ## prior signal covariance operator (power operator)
-    S = power_operator(k_space,spec=power,pindex=powerindex)
+    S = power_operator(k_space,spec=power)
     ## projection operator to the spectral bands
-    Sk = S.get_projection_operator(pindex=powerindex)
+    Sk = S.get_projection_operator()
     ## the Gaussian random field generated from its power operator S
     s = S.get_random_field(domain=s_space,target=k_space)
 
@@ -169,7 +167,7 @@ def run(space=r1,s2n=3,nvar=None,**kwargs):
             the noise variance, `nvar` will be calculated according to
             `s2n` if not specified (default: None)
     """
-    global z,s_space,k,k_space,p,d_space,power,powerindex,powerundex,kindex,rho,S,Sk,R,N,Nj,D,s,n,d,j,m
+    global s_space,k_space,d_space,power,S,Sk,R,N,Nj,D,s,n,d,j,m
 
     ## setting up signal, noise, data and the operators S, N and R
     setup(space,s2n=s2n,nvar=nvar)
@@ -196,7 +194,7 @@ def run(space=r1,s2n=3,nvar=None,**kwargs):
     m.plot(title="reconstructed map",vmin=np.min(s.val),vmax=np.max(s.val),**kwargs)
 
     ## power spectrum
-#    s.plot(title="power spectra",power=True,other=(m,power),mono=False,kindex=kindex)
+#    s.plot(title="power spectra",power=True,other=(m,power),mono=False)
 
     ## uncertainty
 #    uncertainty = D.hat(bare=True,nrun=D.domain.dim()//4,target=k_space)
@@ -234,7 +232,7 @@ def run_critical(space=r2,s2n=3,nvar=None,q=1E-12,alpha=1,perception=[1,0],**kwa
         infer_power
 
     """
-    global z,s_space,k,k_space,p,d_space,power,powerindex,powerundex,kindex,rho,S,Sk,R,N,Nj,D,s,n,d,j,m
+    global s_space,k_space,d_space,power,rho,S,Sk,R,N,Nj,D,s,n,d,j,m
 
     ## setting up signal, noise, data and the operators S, N and R
     setup(space,s2n=s2n,nvar=nvar)
@@ -245,7 +243,7 @@ def run_critical(space=r2,s2n=3,nvar=None,q=1E-12,alpha=1,perception=[1,0],**kwa
     j = R.adjoint_times(N.inverse_times(d))
 
     ## unknown power spectrum, the power operator is given an initial value
-    S.set_power(42,bare=True,pindex=powerindex) ## The answer is 42. I double checked.
+    S.set_power(42,bare=True) ## The answer is 42. I double checked.
     ## the power spectrum is drawn from the first guess power operator
     pk = S.get_power(bare=False) ## non-bare(!)
 
@@ -266,14 +264,14 @@ def run_critical(space=r2,s2n=3,nvar=None,q=1E-12,alpha=1,perception=[1,0],**kwa
         b1 = Sk.pseudo_tr(m) ## == Sk(m).pseudo_dot(m), but faster
         b2 = Sk.pseudo_tr(D,nrun=np.sqrt(Sk.domain.dim())//4) ## probing of the partial traces of D
         pk_new = (2*q+b1+perception[0]*b2)/(rho+2*(alpha-1+perception[1])) ## non-bare(!)
-        pk_new = smooth_power(pk_new,kindex,mode="2s",exclude=min(8,len(power))) ## smoothing
+        pk_new = smooth_power(pk_new,domain=k_space,mode="2s",exclude=min(8,len(rho))) ## smoothing
         ## the power operator is given the new spectrum
-        S.set_power(pk_new,bare=False,pindex=powerindex) ## auto-updates D
+        S.set_power(pk_new,bare=False) ## auto-updates D
 
         ## check convergence
         log_change = np.max(np.abs(log(pk_new/pk)))
         if(log_change<=0.01):
-            note.cprint("NOTE: accuracy reached in iteration %u."%(iteration))
+            note.cprint("NOTE: desired accuracy reached in iteration %u."%(iteration))
             break
         else:
             note.cprint("NOTE: log-change == %4.3f ( > 1%% ) in iteration %u."%(log_change,iteration))
@@ -290,7 +288,7 @@ def run_critical(space=r2,s2n=3,nvar=None,q=1E-12,alpha=1,perception=[1,0],**kwa
     m.plot(title="reconstructed map",vmin=np.min(s.val),vmax=np.max(s.val),**kwargs)
 
     ## power spectrum
-    s.plot(title="power spectra",power=True,other=(S.get_power(pundex=powerundex),power),mono=False,kindex=kindex)
+    s.plot(title="power spectra",power=True,other=(S.get_power(),power),mono=False)
 
     ## uncertainty
 #    uncertainty = D.hat(bare=True,nrun=D.domain.dim()//4,target=k_space)

demos/demo_faraday.py

@@ -91,8 +91,7 @@ def run(projection=False, power=False):
     m1 = m0.transform(l)
     if(projection):
         ## define projection operator
-        powerindex = l.get_power_index()
-        Sk = projection_operator(l, assign=powerindex)
+        Sk = projection_operator(l)
         ## project quadrupole
         m2 = Sk(m0, band=2)
 
@@ -116,7 +115,7 @@ def run(projection=False, power=False):
     m3.plot(title=r"$m$ on a spherical Gauss-Legendre grid", vmin=-4, vmax=4, cmap=ncmap.fm())
     m4.plot(title=r"$m$ on a regular 2D grid", vmin=-4, vmax=4, cmap=ncmap.fm())
     if(power):
-        m5.plot(title=r"(restricted) Fourier power spectrum of $m$", vmin=1E-3, vmax=1E+0, mono=False)
+        m5.plot(title=r"(restricted, binned) Fourier power spectrum of $m$", vmin=1E-3, vmax=1E+0, mono=False, log=True)
 
 ##=============================================================================
 

powerspectrum.py

@@ -68,7 +68,7 @@ def draw_vector_nd(axes,dgrid,ps,symtype=0,fourier=False,zerocentered=False,kpac
         kdict = np.fft.fftshift(nkdict(axes,dgrid,fourier))
         klength = nklength(kdict)
     else:
-        kdict = kpack[1][kpack[0]]
+        kdict = kpack[1][np.fft.ifftshift(kpack[0],axes=shiftaxes(zerocentered,st_to_zero_mode=False))]
         klength = kpack[1]
 
     #output is in position space
@@ -329,6 +329,120 @@ def get_power_index(axes,dgrid,zerocentered,irred=False,fourier=True):
         return ind
 
 
+def get_power_indices(axes,dgrid,zerocentered,fourier=True):
+    """
+        Returns the index of the Fourier grid points in a numpy
+        array, ordered following the zerocentered flag.
+
+        Parameters
+        ----------
+        axes : ndarray
+            An array with the length of each axis.
+
+        dgrid : ndarray
+            An array with the pixel length of each axis.
+
+        zerocentered : bool
+            Whether the output array should be zerocentered, i.e. starting with
+            negative Fourier modes going over the zero mode to positive modes,
+            or not zerocentered, where zero, positive and negative modes are
+            simpy ordered consecutively.
+
+        irred : bool : *optional*
+            If True, the function returns an array of all k-vector lengths and
+            their degeneracy factors. If False, just the power index array is
+            returned.
+
+        fourier : bool : *optional*
+            Whether the output should be in Fourier space or not
+            (default=False).
+
+        Returns
+        -------
+        index, klength, rho : ndarrays
+            Returns the power index array, an array of all k-vector lengths and
+            their degeneracy factors.
+
+    """
+
+    ## kdict, klength
+    if(np.any(zerocentered==False)):
+        kdict = np.fft.fftshift(nkdict(axes,dgrid,fourier),axes=shiftaxes(zerocentered,st_to_zero_mode=True))
+    else:
+        kdict = nkdict(axes,dgrid,fourier)
+    klength = nklength(kdict)
+    ## output
+    ind = np.empty(axes,dtype=np.int)
+    rho = np.zeros(klength.shape,dtype=np.int)
+    for ii in np.ndindex(kdict.shape):
+        ind[ii] = np.searchsorted(klength,kdict[ii])
+        rho[ind[ii]] += 1
+    return ind,klength,rho
+
+
+def bin_power_indices(pindex,kindex,rho,log=False,nbin=None,binbounds=None):
+    """
+        Returns the (re)binned power indices associated with the Fourier grid.
+
+        Parameters
+        ----------
+        pindex : ndarray
+            Index of the Fourier grid points in a numpy.ndarray ordered
+            following the zerocentered flag (default=None).
+        kindex : ndarray
+            Array of all k-vector lengths (default=None).
+        rho : ndarray
+            Degeneracy of the Fourier grid, indicating how many k-vectors in
+            Fourier space have the same length (default=None).
+        log : bool
+            Flag specifying if the binning is performed on logarithmic scale
+            (default: False).
+        nbin : integer
+            Number of used bins (default: None).
+        binbounds : {list, array}
+            Array-like inner boundaries of the used bins (default: None).
+
+        Returns
+        -------
+        pindex, kindex, rho : ndarrays
+            The (re)binned power indices.
+
+    """
+    ## boundaries
+    if(binbounds is not None):
+        binbounds = np.sort(binbounds)
+    ## equal binning
+    else:
+        if(log is None):
+            log = False
+        if(log):
+            k = np.r_[0,np.log(kindex[1:])]
+        else:
+            k = kindex
+        dk = np.max(k[2:]-k[1:-1]) ## minimal dk
+        if(nbin is None):
+            nbin = int((k[-1]-0.5*(k[2]+k[1]))/dk-0.5) ## maximal nbin
+        else:
+            nbin = min(int(nbin),int((k[-1]-0.5*(k[2]+k[1]))/dk+2.5))
+            dk = (k[-1]-0.5*(k[2]+k[1]))/(nbin-2.5)
+        binbounds = np.r_[0.5*(3*k[1]-k[2]),0.5*(k[1]+k[2])+dk*np.arange(nbin-2)]
+        if(log):
+            binbounds = np.exp(binbounds)
+    ## reordering
+    reorder = np.searchsorted(binbounds,kindex)
+    rho_ = np.zeros(len(binbounds)+1,dtype=rho.dtype)
+    kindex_ = np.empty(len(binbounds)+1,dtype=kindex.dtype)
+    for ii in range(len(reorder)):
+        if(rho_[reorder[ii]]==0):
+            kindex_[reorder[ii]] = kindex[ii]
+            rho_[reorder[ii]] += rho[ii]
+        else:
+            kindex_[reorder[ii]] = (kindex_[reorder[ii]]*rho_[reorder[ii]]+kindex[ii]*rho[ii])/(rho_[reorder[ii]]+rho[ii])
+            rho_[reorder[ii]] += rho[ii]
+
+    return reorder[pindex],kindex_,rho_
+
+
 def nhermitianize(field,zerocentered):
 
     """

setup.py

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

smoothing.py

@@ -20,9 +20,10 @@
 ## along with this program. If not, see <http://www.gnu.org/licenses/>.
 
 ## TODO: optimize
+## TODO: doc strings
 
 from __future__ import division
-import gfft as gf
+#import gfft as gf
 import numpy as np
 
 
@@ -63,8 +64,8 @@ def smooth_power(power, k, exclude=1, smooth_length=None):
     pbinned = pbinned / counter
 
     nmirror = int(5 * smooth_length / dk) + 2
-    zpbinned = np.r_[(2 * pbinned[0] - pbinned[1:nmirror][::-1]), pbinned,
-                     (2 * pbinned[-1] - pbinned[-nmirror:-1][::-1])]
+    zpbinned = np.r_[np.exp(2 * np.log(pbinned[0]) - np.log(pbinned[1:nmirror][::-1])), pbinned,
+                     np.exp(2 * np.log(pbinned[-1])- np.log(pbinned[-nmirror:-1][::-1]))]
     zpbinned = np.maximum(0, zpbinned)
 
     tpbinned = np.fft.fftshift(np.fft.fft(zpbinned))
@@ -100,7 +101,7 @@ def smooth_power_bf(power, k, exclude=1, smooth_length=None):
         smooth_length = k[1]-k[0]
 
     nmirror = int(5*smooth_length/(k[1]-k[0]))+2
-    mpower = np.r_[(2*power[0]-power[1:nmirror][::-1]),power,(2*power[-1]-power[-nmirror:-1][::-1])]
+    mpower = np.r_[np.exp(2*np.log(power[0])-np.log(power[1:nmirror][::-1])),power,np.exp(2*np.log(power[-1])-np.log(power[-nmirror:-1][::-1]))]
     mk = np.r_[(2*k[0]-k[1:nmirror][::-1]),k,(2*k[-1]-k[-nmirror:-1][::-1])]
     mdk = np.r_[0.5*(mk[1]-mk[0]),0.5*(mk[2:]-mk[:-2]),0.5*(mk[-1]-mk[-2])]
 
@@ -142,7 +143,7 @@ def smooth_power_2s(power, k, exclude=1, smooth_length=None):
         smooth_length = k[1]-k[0]
 
     nmirror = int(5*smooth_length/(k[1]-k[0]))+2
-    mpower = np.r_[(2*power[0]-power[1:nmirror][::-1]),power,(2*power[-1]-power[-nmirror:-1][::-1])]
+    mpower = np.r_[np.exp(2*np.log(power[0])-np.log(power[1:nmirror][::-1])),power,np.exp(2*np.log(power[-1])-np.log(power[-nmirror:-1][::-1]))]
     mk = np.r_[(2*k[0]-k[1:nmirror][::-1]),k,(2*k[-1]-k[-nmirror:-1][::-1])]
     mdk = np.r_[0.5*(mk[1]-mk[0]),0.5*(mk[2:]-mk[:-2]),0.5*(mk[-1]-mk[-2])]
 
@@ -207,11 +208,16 @@ def smooth_field(val, fourier, zero_center, enforce_hermitian_symmetry, vol, \
             tfield = val
             vol = 1/np.array(val.shape)/vol
         else:
-#
-            tfield = gf.gfft(val, ftmachine='fft', \
-                in_zero_center=zero_center, out_zero_center=True, \
-                enforce_hermitian_symmetry=enforce_hermitian_symmetry, W=-1, \
-                alpha=-1, verbose=False)
+
+            if(zero_center):
+                tfield = np.fft.fftshift(np.fft.fftn(np.fft.fftshift(val)))
+            else:
+                tfield = np.fft.fftshift(np.fft.fftn(val))
+#            # Transform back to the signal space using GFFT.
+#            tfield = gf.gfft(val, ftmachine='fft', \
+#                in_zero_center=zero_center, out_zero_center=True, \
+#                enforce_hermitian_symmetry=enforce_hermitian_symmetry, W=-1, \
+#                alpha=-1, verbose=False)
 
         # Construct the Fourier transformed smoothing kernel
         tkernel = gaussian_kernel(val.shape, vol, smooth_length)
@@ -221,12 +227,16 @@ def smooth_field(val, fourier, zero_center, enforce_hermitian_symmetry, vol, \
         if(fourier):
             sfield = tfield
         else:
-#
-            # Transform back to the signal space using GFFT.
-            sfield = gf.gfft(tfield, ftmachine='ifft', \
-                in_zero_center=True, out_zero_center=zero_center, \
-                enforce_hermitian_symmetry=enforce_hermitian_symmetry, W=-1, \
-                alpha=-1, verbose=False)
+
+            if(zero_center):
+                sfield = np.fft.ifftshift(np.fft.ifftn(np.fft.ifftshift(tfield)))
+            else:
+                sfield = np.fft.ifftn(np.fft.ifftshift(tfield))
+#            # Transform back to the signal space using GFFT.
+#            sfield = gf.gfft(tfield, ftmachine='ifft', \
+#                in_zero_center=True, out_zero_center=zero_center, \
+#                enforce_hermitian_symmetry=enforce_hermitian_symmetry, W=-1, \
+#                alpha=-1, verbose=False)
 
         return sfield