Python/numpy tricky slicing problem

Question

I have a problem with some numpy stuff. I need a numpy array to behave in an unusual manner by returning a slice as a view of the data I have sliced, not a copy. So heres an example of what I want to do:

Say we have a simple array like this:

a = array([1, 0, 0, 0])

I would like to update consecutive entries in the array (moving left to right) with the previous entry from the array, using syntax like this:

a[1:] = a[0:3]

This would get the following result:

a = array([1, 1, 1, 1])

Or something like this:

a[1:] = 2*a[:3]
# a = [1,2,4,8]

To illustrate further I want the following kind of behaviour:

for i in range(len(a)):
    if i == 0 or i+1 == len(a): continue
    a[i+1] = a[i]

Except I want the speed of numpy.

The default behavior of numpy is to take a copy of the slice, so what I actually get is this:

a = array([1, 1, 0, 0])

I already have this array as a subclass of the ndarray, so I can make further changes to it if need be, I just need the slice on the right hand side to be continually updated as it updates the slice on the left hand side.

Am I dreaming or is this magic possible?

Update: This is all because I am trying to use Gauss-Seidel iteration to solve a linear algebra problem, more or less. It is a special case involving harmonic functions, I was trying to avoid going into this because its really not necessary and likely to confuse things further, but here goes.

The algorithm is this:

while not converged:
    for i in range(len(u[:,0])):
        for j in range(len(u[0,:])):
            # skip over boundary entries, i,j == 0 or len(u)
            u[i,j] = 0.25*(u[i-1,j] + u[i+1,j] + u[i, j-1] + u[i,j+1])

Right? But you can do this two ways, Jacobi involves updating each element with its neighbours without considering updates you have already made until the while loop cycles, to do it in loops you would copy the array then update one array from the copied array. However Gauss-Seidel uses information you have already updated for each of the i-1 and j-1 entries, thus no need for a copy, the loop should essentially 'know' since the array has been re-evaluated after each single element update. That is to say, every time we call up an entry like u[i-1,j] or u[i,j-1] the information calculated in the previous loop will be there.

I want to replace this slow and ugly nested loop situation with one nice clean line of code using numpy slicing:

u[1:-1,1:-1] = 0.25(u[:-2,1:-1] + u[2:,1:-1] + u[1:-1,:-2] + u[1:-1,2:])

But the result is Jacobi iteration because when you take a slice: u[:,-2,1:-1] you copy the data, thus the slice is not aware of any updates made. Now numpy still loops right? Its not parallel its just a faster way to loop that looks like a parallel operation in python. I want to exploit this behaviour by sort of hacking numpy to return a pointer instead of a copy when I take a slice. Right? Then every time numpy loops, that slice will 'update' or really just replicate whatever happened in the update. To do this I need slices on both sides of the array to be pointers.

Anyway if there is some really really clever person out there that awesome, but I've pretty much resigned myself to believing the only answer is to loop in C.

Answer 1

You could have a look at np.lib.stride_tricks.

There is some information in these excellent slides: http://mentat.za.net/numpy/numpy_advanced_slides/

with stride_tricks starting at slide 29.

I'm not completely clear on the question though so can't suggest anything more concrete - although I would probably do it in cython or fortran with f2py or with weave. I'm liking fortran more at the moment because by the time you add all the required type annotations in cython I think it ends up looking less clear than the fortran.

There is a comparison of these approaches here:

www. scipy. org/ PerformancePython

(can't post more links as I'm a new user) with an example that looks similar to your case.

Answer 2

Just use a loop. I can't immediately think of any way to make the slice operator behave the way you're saying you want it to, except maybe by subclassing numpy's array and overriding the appropriate method with some sort of Python voodoo... but more importantly, the idea that a[1:] = a[0:3] should copy the first value of a into the next three slots seems completely nonsensical to me. I imagine that it could easily confuse anyone else who looks at your code (at least the first few times).

Answer 3

In the end I came up with the same problem as you. I had to resort to use Jacobi iteration and weaver:

 while (iter_n < max_time_steps):
        expr = "field[1:-1, 1:-1] = (field[2:, 1:-1] "\
                                                      "+ field[:-2, 1:-1]+"\
                                                      "field[1:-1, 2:] +"\
                                                      "field[1:-1, :-2] )/4."                                       

        weave.blitz(expr, check_size=0)

         #Toroidal conditions
        field[:,0] = field[:,self.flow.n_x - 2]
        field[:,self.flow.n_x -1] = field[:,1]

        iter_n = iter_n + 1

It works and is fast, but is not Gauss-Seidel, so convergence can be a bit tricky. The only option of doing Gauss-Seidel as a traditional loop with indexes.

Answer 4

It is not the correct logic. I'll try to use letters to explain it.

Image array = abcd with a,b,c,d as elements.
Now, array[1:] means from the element in position 1 (starting from 0) on.
In this case:bcd and array[0:3] means from the character in position 0 up to the third character (the one in position 3-1) in this case: 'abc'.

Writing something like:
array[1:] = array[0:3]

means: replace bcd with abc

To obtain the output you want, now in python, you should use something like:

a[1:] = a[0]

Answer 5

It must have something to do with assigning a slice. Operators, however, as you may already know, do follow your expected behavior:

>>> a = numpy.array([1,0,0,0])
>>> a[1:]+=a[:3]
>>> a
array([1, 1, 1, 1])

If you already have zeros in your real-world problem where your example does, then this solves it. Otherwise, at added cost, set them to zero either by multiplying by zero or assigning to zero, (whichever is faster)

edit: I had another thought. You may prefer this:

numpy.put(a,[1,2,3],a[:3]) 

Answer 6

i would suggest cython instead of looping in c. there might be some fancy numpy way of getting your example to work using a lot of intermediate steps... but since you know how to write it in c already, just write that quick little bit as a cython function and let cython's magic make the rest of the work easy for you.