Efficient ways to iterate over several numpy arrays and process current and previous elements?
<p>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 href="http://stackoverflow.com/questions/38207528/what-are-the-efficient-ways-to-loop-over-vectors-along-a-specified-axis-in-numpy">a comment here</a>). 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.</p>

<p>Let's say there are N (N is small, usually 4, might be up to 7) 1-D numpy arrays of <code>float128</code> in a list <code>X</code>, 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).</p>

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

<p>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:</p>

<pre><code>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) )
</code></pre>

<p>This is clearly very inefficient and besides not a Pythonic way.</p>

<p>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:</p>

<pre><code># 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?
</code></pre>

<p>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:</p>

<pre><code>res = []
inds, = numpy.where((secaxis[:-1] &lt; 0) * (secaxis[1:] &gt; 0) +
                   (secaxis[:-1] &gt; 0) * (secaxis[1:] &lt; 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) )
</code></pre>

<p>But this is seemingly done in 7 + 2*(N - 1) passes, moreover, I'm not sure about <code>secaxis[inds]</code> 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?).</p>

<p>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:</p>

<pre><code>def filt(x):
  return (x[0] &lt; 0 and x[1] &gt; 0) or (x[0] &gt; 0 and x[1] &lt; 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) )
</code></pre>

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

<p>So, I have following questions:</p>

<ol>
<li>Of all these methods is the third one indeed the best? If so, what can be done to impove the last one?</li>
<li>Are there any other, better ones yet? </li>
<li>Out of sheer curiosity, is there a way to solve the problem using nditer?</li>
<li>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 <code>secaxis</code> array into RAM, keep others on disk and use third method?</li>
<li>(bonus question) List of equal in length 1-D numpy arrays comes from loading N <code>.npy</code> 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?</li>
</ol>

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

<pre><code>import numpy as np

a = np.random.randn(100000)

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

def method2(a):
    inds, = np.where((a[:-1] &lt; 0) * (a[1:] &gt; 0) +
                       (a[:-1] &gt; 0) * (a[1:] &lt; 0))
    return inds + 1

def method3(a):
    m = a &lt; 0
    p = a &gt; 0
    return np.where((m[:-1] &amp; p[1:]) | (p[:-1] &amp; m[1:]))[0] + 1

def method4(a):
    return np.where(np.diff(a &gt;= 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)
</code></pre>

<p>the <code>%timeit</code> result:</p>

<pre><code>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
</code></pre>

python arrays performance numpy itertools