WIP - iteration and testing

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

+import numpy as np
+
+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
+	
+    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]
+
+    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)
+
+#    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):
+
+    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
+	
+#    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]
+
+#    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)
+    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)
+
+    
\ No newline at end of file

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

+import numpy as np
+
+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
+	
+    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]))]
+
+    print "mpower", mpower
+    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
\ No newline at end of file

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

+import pyximport; pyximport.install()
+import extended
+
+import numpy as np
+
+
+print "///////////////////////////////////////First thing ////////////////////////"
+
+k=np.sqrt(np.arange(600))
+power = np.ones(600)
+power[300]=1000
+print power, k
+smooth = extended.smooth_power_2s(power, k)
+
+print "Smoooooth", smooth
+
+print "///////////////////////////////////////Final thing ////////////////////////"
+print "smooth.len == power.len" , len(smooth), len(power), len(power)==len(smooth)
\ No newline at end of file

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

+import pyximport; pyximport.install()
+import extended
+
+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
+print power, k, power.shape
+smooth = extended.smooth_something(datablock=power, axis=(2), startindex=None, endindex=None, kernelfunction=extended.smoothie, k=ksq)
+
+print "Smoooooth", smooth
+
+print "///////////////////////////////////////Final thing ////////////////////////"
+print "smooth.len == power.len" , smooth.shape, power.shape, power.shape==smooth.shape
\ No newline at end of file