WIP - testing apply_along_axis functionality

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

@@ -50,8 +50,7 @@ def GaussianKernel(mpower, mk, mu,smooth_length):
     return np.sum(C * mpower) / np.sum(C)
 
 
-def smoothie(power, startindex, endindex, k, kernelfunction, exclude=1,
-                                                          smooth_length=None):
+def smoothie(power, startindex, endindex, k, exclude=1, smooth_length=None):
 
     if smooth_length == 0:
         # No smoothing requested, just return the input array.
@@ -70,7 +69,7 @@ def smoothie(power, startindex, endindex, k, kernelfunction, exclude=1,
     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],
+        p_smooth[i-startindex] = GaussianKernel(power[l:u], k[l:u], k[i],
                                           smooth_length)
 
     if (exclude > 0):
@@ -87,6 +86,6 @@ def smooth_something(datablock, axis=0, startindex=None, endindex=None,
     print kernelfunction
     return np.apply_along_axis(smoothie, axis, datablock,
                                startindex=startindex, endindex=endindex, k=k,
-                               smooth_length=sigma, kernelfunction=kernelfunction)
+                               smooth_length=sigma)
 
     
\ No newline at end of file

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

 import pyximport; pyximport.install()
 import extended
+import util
 
 import numpy as np
 
@@ -25,26 +26,33 @@ mirrorsize=7
 startindex=mirrorsize/2
 endindex=n-mirrorsize/2
 print power, k, power.shape
-smooth = extended.smooth_something(datablock=power, axis=(2),
-                                   startindex=startindex, endindex=endindex,
-                                   kernelfunction=extended.GaussianKernel, k=k,
-                                   sigma=sigma)
-
+# smooth = extended.smooth_something(datablock=power, axis=(2),
+#                                    startindex=startindex, endindex=endindex,
+#                                    kernelfunction=extended.GaussianKernel, k=k,
+#                                    sigma=sigma)
+smooth = util.apply_along_axis(extended.smoothie, (2), power,
+                               startindex=startindex, endindex=endindex, k=k,
+                               smooth_length=sigma)
 print "Smoooooth", smooth
 
-doublesmooth = extended.smooth_something(datablock=smooth, axis=(1),
-                                   startindex=startindex, endindex=endindex,
-                                         kernelfunction=extended.GaussianKernel,
-                                         k=k, sigma=sigma)
+# doublesmooth = extended.smooth_something(datablock=smooth, axis=(1),
+#                                    startindex=startindex, endindex=endindex,
+#                                          kernelfunction=extended.GaussianKernel,
+#                                          k=k, sigma=sigma)
+doublesmooth = util.apply_along_axis(extended.smoothie, (1), smooth,
+                               startindex=startindex, endindex=endindex, k=k,
+                               smooth_length=sigma)
 
 print "DoubleSmooth", doublesmooth
 
-tripplesmooth = extended.smooth_something(datablock=doublesmooth, axis=(0),
-                                   startindex=startindex, endindex=endindex,
-                                         kernelfunction=extended.GaussianKernel,
-                                         k=k, sigma=sigma)
-
+# tripplesmooth = extended.smooth_something(datablock=doublesmooth, axis=(0),
+#                                    startindex=startindex, endindex=endindex,
+#                                          kernelfunction=extended.GaussianKernel,
+#                                          k=k, sigma=sigma)
+tripplesmooth = util.apply_along_axis(extended.smoothie, (0), doublesmooth,
+                               startindex=startindex, endindex=endindex, k=k,
+                               smooth_length=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
+print "smooth.len == power.len" , tripplesmooth.shape, power.shape, power.shape==smooth.shape
\ No newline at end of file

nifty/operators/smoothing_operator/Cython/util.py

+import numpy as np
+
+def apply_along_axis(func1d, axis, arr, *args, **kwargs):
+    """
+    Apply a function to 1-D slices along the given axis.
+
+    Execute `func1d(a, *args)` where `func1d` operates on 1-D arrays and `a`
+    is a 1-D slice of `arr` along `axis`.
+
+    Parameters
+    ----------
+    func1d : function
+        This function should accept 1-D arrays. It is applied to 1-D
+        slices of `arr` along the specified axis.
+    axis : integer
+        Axis along which `arr` is sliced.
+    arr : ndarray
+        Input array.
+    args : any
+        Additional arguments to `func1d`.
+    kwargs: any
+        Additional named arguments to `func1d`.
+
+        .. versionadded:: 1.9.0
+
+
+    Returns
+    -------
+    apply_along_axis : ndarray
+        The output array. The shape of `outarr` is identical to the shape of
+        `arr`, except along the `axis` dimension, where the length of `outarr`
+        is equal to the size of the return value of `func1d`.  If `func1d`
+        returns a scalar `outarr` will have one fewer dimensions than `arr`.
+    """
+    arr = np.asarray(arr)
+    nd = arr.ndim
+    if axis < 0:
+        axis += nd
+    if (axis >= nd):
+        raise ValueError("axis must be less than arr.ndim; axis=%d, rank=%d."
+            % (axis, nd))
+    ind = [0]*(nd-1)
+    i = np.zeros(nd, 'O')
+    indlist = list(range(nd))
+    indlist.remove(axis)
+    i[axis] = slice(None, None)
+    outshape = np.asarray(arr.shape).take(indlist)
+    i.put(indlist, ind)
+    res = func1d(arr[tuple(i.tolist())], *args, **kwargs)
+
+    Ntot = np.product(outshape)
+    holdshape = outshape
+    outshape = list(arr.shape)
+    outshape[axis] = len(res)
+    outarr = np.zeros(outshape, np.asarray(res).dtype)
+    outarr[tuple(i.tolist())] = res
+    k = 1
+    while k < Ntot:
+        # increment the index
+        ind[-1] += 1
+        n = -1
+        while (ind[n] >= holdshape[n]) and (n > (1-nd)):
+            ind[n-1] += 1
+            ind[n] = 0
+            n -= 1
+        i.put(indlist, ind)
+        res = func1d(arr[tuple(i.tolist())], *args, **kwargs)
+        outarr[tuple(i.tolist())] = res
+        k += 1
+    return outarr
\ No newline at end of file