2

I have a problem with the speed of my program. I want to calculate the average of four neighbors in a huge array. Here is a part of my code. Do you have any ideas how to change the last line? Or should I use another array?

for a in np.arange(100000):
    for x in np.arange(size):
        for y in np.arange(size):
            if unchangeableflag[x*size+y] == 0:
                vnew[x*size+y] = (v[(x+1)*size+y] + v[(x-1)*size+y] + v[x*size+y+1] + v[x*size+y-1]) / 4.0
unwind
  • 391,730
  • 64
  • 469
  • 606
kame
  • 20,848
  • 33
  • 104
  • 159

3 Answers3

4

You would not need the loop at all. Assuming v, vnew and unchangeableflag are 1-d arrays with size*size entries, you can do

v = v.reshape(size, size)
vnew = vnew.reshape(size, size)
unchangeableflag = unchangeableflag.reshape(size, size)
average = v[1:-1, 2:]
average += v[1:-1, :-2] 
average += v[2:, 1:-1]
average += v[-2:, 1:-1]
average /= 4.0
vnew[1:-1, 1:-1][unchangeableflag[1:-1, 1:-1] == 0] = average

But what are you actually trying to achieve? This looks suspiciously like you could get away with some application of the discrete Laplacian.

(Note that this assumes that v contains floating point numbers. If the dtype of `v´ is sime integral type, you need a slight modification.)

Sven Marnach
  • 574,206
  • 118
  • 941
  • 841
  • I think that i loose always the outermost line. Right? – kame Jan 28 '11 at 18:08
  • 1
    @kame: Your own code has somewhat undefined behaviour on the boundaries. If you define what should be done on the boundaries, it can of course be incorporated. – Sven Marnach Jan 28 '11 at 20:40
  • 2
    @kame - The numpy and scipy convolve functions allow some different approaches to how the boundaries are treated. – tom10 Jan 28 '11 at 20:41
3

You should consider using SciPy's convolution filter or the generic_filter. This is still computationally intensive, but way faster than looping. Normally, when doing this type of averaging, the central element is also included. Note that these solutions also apply to multi-dimensional arrays.

from scipy import ndimage
footprint = scipy.array([[0,0.25,0],[0.25,0,0.25],[0,0.25,0]])
filtered_array = scipy.convolve(array, footprint)

OR

from scipy import ndimage
def myfunction(window):
    return (window[0,1] + window[1,0] + window[1,2] + window[2,1]) / 4
filtered_array = scipy.generic_filter(array, myfunction, size=3)
Benjamin
  • 11,560
  • 13
  • 70
  • 119
  • These solutions are concise, but not fastest. Especially the second solution is quite slow. In the first solution, you probably want to use `ndimage.convolve()` instead of `scipy.convolve`, since the latter is for 1d arrays only. Furthermore, both solution disregard the mysterious `unchangeableflag` in the original post. – Sven Marnach Jan 28 '11 at 16:37
  • convolve1d is for 1d arrays, convolve is the multi-dimentional version. Because of the way I import ndimage, I don't think I need to specify scipy.ndimage.convolve. I'm less concerned with keeping and fixing the original code than with providing an example of a simpler and more effective alternative... – Benjamin Jan 28 '11 at 16:45
  • 4
    `scipy.convolve` and `scipy.ndimage.convolve` are _not_ the same function at all!! Anything in the "bare" `scipy` namespace is just there for historical reasons (i.e. they're just numpy functions, though some have different defaults). _All_ of scipy's functionality is in the other namespaces such as `scipy.ndimage`, `scipy.special`, `scipy.signal`, etc. Other than that, though, I agree, this can certainly be done with a convolution. – Joe Kington Jan 28 '11 at 16:59
  • Sorry if I sounded a bit harsh... I just wanted to make the difference clear... The fact that `import scipy` doesn't actually do anything (other than importing `numpy`) is confusing to a lot of people at first (Or it was to me for quite awhile, anyway). – Joe Kington Jan 28 '11 at 17:06
1

I am unsure but you can remove some invariant.

for a in np.arange(100000):
    for x in np.arange(size):
        for y in np.arange(size):
            t = x*size+y
            if unchangeableflag[t] == 0:
                vnew[t] = (v[t+size] + v[t-size] + v[t+1] + v[t-1]) / 4.0
VGE
  • 4,171
  • 18
  • 17
  • Instead of dividing by 4, you could use shift the bits to the right twice: vnet[t] = (v[t+size] + v[t-size] + v[t+1] + v[t-1]); vnet[t] >>= 2; – tdobek Jan 28 '11 at 15:21
  • @tdobek What? Can you explain this? – kame Jan 28 '11 at 15:27
  • 3
    @tdobek: These are floating point computations. Shifting won't work (nor would it be any faster even if it did). – Sven Marnach Jan 28 '11 at 15:40
  • @Sven Marnach: In this case my idea is void. And it's Python so yeah shifting bits won't necessarily help. Sorry! – tdobek Jan 28 '11 at 16:00