Efficient ways to iterate over several numpy arrays and process current and previous elements?

Question

I've read a lot about different techniques for iterating over numpy arrays recently and it seems that consensus is not to iterate at all (for instance, see a comment here). There are several similar questions on SO, but my case is a bit different as I have to combine "iterating" (or not iterating) and accessing previous values.

Let's say there are N (N is small, usually 4, might be up to 7) 1-D numpy arrays of float128 in a list X, all arrays are of the same size. To give you a little insight, these are data from PDE integration, each array stands for one function, and I would like to apply a Poincare section. Unfortunately, the algorithm should be both memory- and time-efficient since these arrays are sometimes ~1Gb each, and there are only 4Gb of RAM on board (I've just learnt about memmap'ing of numpy arrays and now consider using them instead of regular ones).

One of these arrays is used for "filtering" the others, so I start with secaxis = X.pop(idx). Now I have to locate pairs of indices where (secaxis[i-1] > 0 and secaxis[i] < 0) or (secaxis[i-1] < 0 and secaxis[i] > 0) and then apply simple algebraic transformations to remaining arrays, X (and save results). Worth mentioning, data shouldn't be wasted during this operation.

There are multiple ways for doing that, but none of them seem efficient (and elegant enough) to me. One is a C-like approach, where you just iterate in a for-loop:

import array # better than lists
res = [ array.array('d') for _ in X ]
for i in xrange(1,secaxis.size):
  if condition: # see above
    co = -secaxis[i-1]/secaxis[i]
    for j in xrange(N):
      res[j].append( (X[j][i-1] + co*X[j][i])/(1+co) )

This is clearly very inefficient and besides not a Pythonic way.

Another way is to use numpy.nditer, but I haven't figured out yet how one accesses the previous value, though it allows iterating over several arrays at once:

# without secaxis = X.pop(idx)
it = numpy.nditer(X)
for vec in it:
  # vec[idx] is current value, how do you get the previous (or next) one?

Third possibility is to first find sought indices with efficient numpy slices, and then use them for bulk multiplication/addition. I prefer this one for now:

res = []
inds, = numpy.where((secaxis[:-1] < 0) * (secaxis[1:] > 0) +
                   (secaxis[:-1] > 0) * (secaxis[1:] < 0))
coefs = -secaxis[inds] / secaxis[inds+1] # array of coefficients
for f in X: # loop is done only N-1 times, that is, 3 to 6
    res.append( (f[inds] + coefs*f[inds+1]) / (1+coefs) )

But this is seemingly done in 7 + 2*(N - 1) passes, moreover, I'm not sure about secaxis[inds] type of addressing (it is not slicing and generally it has to find all elements by indices just like in the first method, doesn't it?).

Finally, I've also tried using itertools and it resulted in monstrous and obscure structures, which might stem from the fact that I'm not very familiar with functional programming:

def filt(x):
  return (x[0] < 0 and x[1] > 0) or (x[0] > 0 and x[1] < 0)
import array
from itertools import izip, tee, ifilter
res = [ array.array('d') for _ in X ] 
iters = [iter(x) for x in X]   # N-1 iterators in a list
prev, curr = tee(izip(*iters)) # 2 similar iterators, each of which
                               # consists of N-1 iterators
next(curr, None) # one of them is now for current value
seciter = tee(iter(secaxis))
next(seciter[1], None)
for x in ifilter(filt, izip(seciter[0], seciter[1], prev, curr)):
  co = - x[0]/x[1]
  for r, p, c in zip(res, x[2], x[3]):
    r.append( (p+co*c) / (1+co) )

Not only this looks very ugly, it also takes an awful lot of time to complete.

So, I have following questions:

1. Of all these methods is the third one indeed the best? If so, what can be done to impove the last one?
2. Are there any other, better ones yet?
3. Out of sheer curiosity, is there a way to solve the problem using nditer?
4. Finally, will I be better off using memmap versions of numpy arrays, or will it probably slow things down a lot? Maybe I should only load secaxis array into RAM, keep others on disk and use third method?
5. (bonus question) List of equal in length 1-D numpy arrays comes from loading N .npy files whose sizes aren't known beforehand (but N is). Would it be more efficient to read one array, then allocate memory for one 2-D numpy array (slight memory overhead here) and read remaining into that 2-D array?

Answer 1

The numpy.where() version is fast enough, you can speedup it a little by method3(). If the > condition can change to >=, you can also use method4().

import numpy as np

a = np.random.randn(100000)

def method1(a):
    idx = []
    for i in range(1, len(a)):
        if (a[i-1] > 0 and a[i] < 0) or (a[i-1] < 0 and a[i] > 0):
            idx.append(i)
    return idx

def method2(a):
    inds, = np.where((a[:-1] < 0) * (a[1:] > 0) +
                       (a[:-1] > 0) * (a[1:] < 0))
    return inds + 1

def method3(a):
    m = a < 0
    p = a > 0
    return np.where((m[:-1] & p[1:]) | (p[:-1] & m[1:]))[0] + 1

def method4(a):
    return np.where(np.diff(a >= 0))[0] + 1

assert np.allclose(method1(a), method2(a))
assert np.allclose(method2(a), method3(a))
assert np.allclose(method3(a), method4(a))

%timeit method1(a)
%timeit method2(a)
%timeit method3(a)
%timeit method4(a)

the %timeit result:

1 loop, best of 3: 294 ms per loop
1000 loops, best of 3: 1.52 ms per loop
1000 loops, best of 3: 1.38 ms per loop
1000 loops, best of 3: 1.39 ms per loop