WIP - simplifications

nifty/operators/smoothing_operator/Cython 2/extended.py

@@ -41,31 +41,23 @@ def smooth_power_2s(power, k, exclude=1, smooth_length=None):
     print "p_smooth3", p_smooth
     print "p_smooth length ", len(p_smooth)
 
-#    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 smoothie(power, startindex, endindex, k, exclude=1, smooth_length=None):
+def GaussianKernel(mpower, mk, mu,smooth_length):
+    C = np.exp(-(mk - mu) ** 2 / (2. * smooth_length ** 2))
+    return np.sum(C * mpower) / np.sum(C)
+
+
+def smoothie(power, startindex, endindex, k, kernelfunction, exclude=1,
+                                                          smooth_length=None):
 
     if smooth_length == 0:
         # No smoothing requested, just return the input array.
         return power
 
+    excluded_power = []
     if (exclude > 0):
         k = k[exclude:]
         excluded_power = np.copy(power[:exclude])
@@ -74,44 +66,27 @@ def smoothie(power, startindex, endindex, k, exclude=1, smooth_length=None):
     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
-	
-#    print "nmirror", nmirror
-
-    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])]
-#    print "mpower", mpower
-    p_smooth = np.empty(mpower.shape)
-#    print "p_smooth", p_smooth
-    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))
-        #print "i", i, "l", l, "u", u
-        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)
-        # print "p_smooth[",i,"] = ",p_smooth[i]
+    p_smooth = np.empty(endindex-startindex)
+    for i in xrange(startindex, endindex):
+        l = max(i-int(2*smooth_length)-1,0)
+        u = min(i+int(2*smooth_length)+2,len(p_smooth))
+        p_smooth[i-startindex] = kernelfunction(power[l:u], k[l:u], k[i],
+                                          smooth_length)
 
-#    print "p_smooth2", " all ", p_smooth
-    p_smooth = p_smooth[nmirror - 1:-nmirror + 1]
-#    print "p_smooth3", p_smooth
-#    print "p_smooth length ", len(p_smooth)
     if (exclude > 0):
         p_smooth = np.r_[excluded_power,p_smooth]
     return p_smooth
 
-def smooth_something(datablock, axis=None, startindex=None, endindex=None, kernelfunction=lambda x:x, k=None):
-    if axis == None:
-        axis = size(datablock)
+
+def smooth_something(datablock, axis=0, startindex=None, endindex=None,
+                     kernelfunction=lambda x:x, k=None, sigma=None):
     if startindex == None:
         startindex=0
     if endindex == None:
         endindex=len(datablock)
     print kernelfunction
-    return np.apply_along_axis(kernelfunction,axis,datablock, startindex=startindex, endindex=endindex, k=k)
+    return np.apply_along_axis(smoothie, axis, datablock,
+                               startindex=startindex, endindex=endindex, k=k,
+                               smooth_length=sigma, kernelfunction=kernelfunction)
 
     
\ No newline at end of file

nifty/operators/smoothing_operator/Cython 2/test2.py

@@ -6,21 +6,45 @@ import numpy as np
 
 print "///////////////////////////////////////First thing ////////////////////////"
 
-ksq=np.sqrt(np.arange(8))
-k=np.arange(8)
-power = np.ones(512).reshape((8,8,8))
-power[0][4][4]=1000
-power[1][4][4]=1000
-power[2][4][4]=1000
-power[3][4][4]=1000
-power[4][4][4]=1000
-power[5][4][4]=1000
-power[6][4][4]=1000
-power[7][4][4]=1000
+n=8
+
+ksq=np.sqrt(np.arange(n))
+kk=np.arange(n)
+power = np.ones(n**3).reshape((n,n,n))
+# power[0][4][4]=1000
+# power[1][4][4]=1000
+# power[2][4][4]=1000
+# power[3][4][4]=1000
+power[n/2][n/2][n/2]=10000
+# power[5][4][4]=1000
+# power[6][4][4]=1000
+# power[7][4][4]=1000
+k = kk
+sigma=k[1]-k[0]
+mirrorsize=7
+startindex=mirrorsize/2
+endindex=n-mirrorsize/2
 print power, k, power.shape
-smooth = extended.smooth_something(datablock=power, axis=(2), startindex=None, endindex=None, kernelfunction=extended.smoothie, k=ksq)
+smooth = extended.smooth_something(datablock=power, axis=(2),
+                                   startindex=startindex, endindex=endindex,
+                                   kernelfunction=extended.GaussianKernel, k=k,
+                                   sigma=sigma)
 
 print "Smoooooth", smooth
 
+doublesmooth = extended.smooth_something(datablock=smooth, axis=(1),
+                                   startindex=startindex, endindex=endindex,
+                                         kernelfunction=extended.GaussianKernel,
+                                         k=k, sigma=sigma)
+
+print "DoubleSmooth", doublesmooth
+
+tripplesmooth = extended.smooth_something(datablock=doublesmooth, axis=(0),
+                                   startindex=startindex, endindex=endindex,
+                                         kernelfunction=extended.GaussianKernel,
+                                         k=k, sigma=sigma)
+
+print "TrippleSmooth", tripplesmooth
+
 print "///////////////////////////////////////Final thing ////////////////////////"
 print "smooth.len == power.len" , smooth.shape, power.shape, power.shape==smooth.shape
\ No newline at end of file