smooth_power routine extended.

nifty_power.py

@@ -97,7 +97,7 @@ def weight_power(domain,spec,power=1,pindex=None,pundex=None):
 
 ##-----------------------------------------------------------------------------
 
-def smooth_power(power,kindex,exclude=1,sigma=-1):
+def smooth_power(power,kindex,mode="2s",exclude=1,sigma=-1):
     """
     Smoothes a power spectrum via convolution with a Gaussian kernel.
 
@@ -109,20 +109,65 @@ def smooth_power(power,kindex,exclude=1,sigma=-1):
     kindex : ndarray
         The array specifying the coordinate indices in conjugate space.
 
+    mode : string
+        Specifices the smoothing mode (default: "2s") :
+
+        - "ff" (smoothing in the harmonic basis using fast Fourier
+            transformations)
+        - "bf" (smoothing in the position basis by brute force)
+        - "2s" (smoothing in the position basis restricted to a 2-`sigma`
+            interval)
+
+
     exclude : scalar
         Excludes the first power spectrum entries from smoothing, indicated by
-        the given integer scalar (default=1, the monopol is not smoothed).
+        the given integer scalar (default = 1, the monopol is not smoothed).
 
-    smooth_length : scalar
-        FWHM of Gaussian convolution kernel.
+    sigma : scalar
+        FWHM of Gaussian convolution kernel (default = -1, `sigma` is set
+        automatically).
 
     Returns
     -------
     smoothpower : ndarray
         The smoothed power spectrum.
 
+    Raises
+    ------
+    KeyError
+        If `mode` is unsupported.
+
     """
-    return gs.smooth_power(power,kindex,exclude=exclude,smooth_length=sigma)
+    if(mode=="2s"):
+        return gs.smooth_power_2s(power,kindex,exclude=exclude,smooth_length=sigma)
+    elif(mode=="ff"):
+        return gs.smooth_power(power,kindex,exclude=exclude,smooth_length=sigma)
+    elif(mode=="bf"):
+        return gs.smooth_power_bf(power,kindex,exclude=exclude,smooth_length=sigma)
+    else:
+        raise KeyError(about._errors.cstring("ERROR: unsupported mode '"+str(mode)+"'."))
 
 ##-----------------------------------------------------------------------------
 
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

smoothing.py

@@ -85,6 +85,98 @@ def smooth_power(power, k, exclude=1, smooth_length=None):
         return power2
 
 
+def smooth_power_bf(power, k, exclude=1, smooth_length=None):
+
+    if smooth_length == 0:
+        # No smoothing requested, just return the input array.
+        return power
+
+    if (exclude > 0):
+        k = k[exclude:]
+        excluded_power = np.copy(power[:exclude])
+        power = power[exclude:]
+
+    if (smooth_length is None) or (smooth_length < 0):
+        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])]
+    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])]
+
+    p_smooth = np.empty(mpower.shape)
+    for i in xrange(len(p_smooth)):
+        C = np.exp(-(mk-mk[i])**2/(2.*smooth_length**2))*mdk
+        p_smooth[i] = np.sum(C*mpower)/np.sum(C)
+
+    p_smooth = p_smooth[nmirror - 1:-nmirror + 1]
+
+#    dk = 0.5*(k[2:] - k[:-2])
+#    dk = np.r_[0.5*(k[1]-k[0]),dk]
+#    dk = np.r_[dk,0.5*(k[-1]-k[-2])]
+#    if (smooth_length is None) or (smooth_length < 0):
+#        smooth_length = k[1]-k[0]
+#
+#    p_smooth = np.empty(power.shape)
+#    for i in xrange(len(p_smooth)):
+#        C = np.exp(-(k-k[i])**2/(2.*smooth_length**2))*dk
+#        p_smooth[i] = np.sum(C*power)/np.sum(C)
+
+    if (exclude > 0):
+        p_smooth = np.r_[excluded_power,p_smooth]
+    return p_smooth
+
+
+def smooth_power_2s(power, k, exclude=1, smooth_length=None):
+
+    if smooth_length == 0:
+        # No smoothing requested, just return the input array.
+        return power
+
+    if (exclude > 0):
+        k = k[exclude:]
+        excluded_power = np.copy(power[:exclude])
+        power = power[exclude:]
+
+    if (smooth_length is None) or (smooth_length < 0):
+        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])]
+    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])]
+
+    p_smooth = np.empty(mpower.shape)
+    for i in xrange(len(p_smooth)):
+        l = i-int(2*smooth_length/mdk[i])-1
+        l = max(l,0)
+        u = i+int(2*smooth_length/mdk[i])+2
+        u = min(u,len(p_smooth))
+        C = np.exp(-(mk[l:u]-mk[i])**2/(2.*smooth_length**2))*mdk[l:u]
+        p_smooth[i] = np.sum(C*mpower[l:u])/np.sum(C)
+
+    p_smooth = p_smooth[nmirror - 1:-nmirror + 1]
+
+#    dk = 0.5*(k[2:] - k[:-2])
+#    dk = np.r_[0.5*(k[1]-k[0]),dk]
+#    dk = np.r_[dk,0.5*(k[-1]-k[-2])]
+#    if (smooth_length is None) or (smooth_length < 0):
+#        smooth_length = k[1]-k[0]
+#
+#    p_smooth = np.empty(power.shape)
+#    for i in xrange(len(p_smooth)):
+#        l = i-int(2*smooth_length/dk[i])-1
+#        l = max(l,0)
+#        u = i+int(2*smooth_length/dk[i])+2
+#        u = min(u,len(p_smooth))
+#        C = np.exp(-(k[l:u]-k[i])**2/(2.*smooth_length**2))*dk[l:u]
+#        p_smooth[i] = np.sum(C*power[l:u])/np.sum(C)
+
+    if (exclude > 0):
+        p_smooth = np.r_[excluded_power,p_smooth]
+    return p_smooth
+
+
 def smooth_field(val, fourier, zero_center, enforce_hermitian_symmetry, vol, \
     smooth_length=0.):
         """