diff --git a/nifty/operators/smoothing_operator/Cython 2/test2.py b/nifty/operators/smoothing_operator/Cython 2/test2.py deleted file mode 100644 index db2507b01417f6ae01768ff3f21c128eb82c1419..0000000000000000000000000000000000000000 --- a/nifty/operators/smoothing_operator/Cython 2/test2.py +++ /dev/null @@ -1,50 +0,0 @@ -import pyximport; pyximport.install() -import extended - -import numpy as np - - -print "///////////////////////////////////////First thing ////////////////////////" - -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=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 diff --git a/nifty/operators/smoothing_operator/Cython 2/extended.py b/nifty/operators/smoothing_operator/Cython/extended.py similarity index 90% rename from nifty/operators/smoothing_operator/Cython 2/extended.py rename to nifty/operators/smoothing_operator/Cython/extended.py index 9c3ef619ebb20f684657c39796b5ab486c62c773..6fbf41aebbebd353007793744a584be522de0981 100644 --- a/nifty/operators/smoothing_operator/Cython 2/extended.py +++ b/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 diff --git a/nifty/operators/smoothing_operator/Cython 2/original.py b/nifty/operators/smoothing_operator/Cython/original.py similarity index 100% rename from nifty/operators/smoothing_operator/Cython 2/original.py rename to nifty/operators/smoothing_operator/Cython/original.py diff --git a/nifty/operators/smoothing_operator/Cython 2/test.py b/nifty/operators/smoothing_operator/Cython/test.py similarity index 100% rename from nifty/operators/smoothing_operator/Cython 2/test.py rename to nifty/operators/smoothing_operator/Cython/test.py diff --git a/nifty/operators/smoothing_operator/Cython/test2.py b/nifty/operators/smoothing_operator/Cython/test2.py new file mode 100644 index 0000000000000000000000000000000000000000..9556ef2bc8553ef651d60b8003feb18682636c28 --- /dev/null +++ b/nifty/operators/smoothing_operator/Cython/test2.py @@ -0,0 +1,58 @@ +import pyximport; pyximport.install() +import extended +import util + +import numpy as np + + +print "///////////////////////////////////////First thing ////////////////////////" + +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=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 = 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 = 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" , tripplesmooth.shape, power.shape, power.shape==smooth.shape \ No newline at end of file diff --git a/nifty/operators/smoothing_operator/Cython/util.py b/nifty/operators/smoothing_operator/Cython/util.py new file mode 100644 index 0000000000000000000000000000000000000000..b8f7b6aa9bb74e1c05283492be27351ef38f16d8 --- /dev/null +++ b/nifty/operators/smoothing_operator/Cython/util.py @@ -0,0 +1,70 @@ +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