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>I'll need to read your post in more detail, but will start with some general observations (from previous iteration questions).</p>

<p>There isn't an efficient way of iterating over arrays in Python, though there are things that slow things down.  I like to distinguish between the iteration mechanism (<code>nditer</code>, <code>for x in A:</code>) and the action (<code>alist.append(...)</code>, <code>x[i+1] += 1</code>).  The big time consumer is usually the action, done many times, not the iteration mechanism itself.</p>

<p>Letting <code>numpy</code> do the iteration in compiled code is the fastest.</p>

<pre><code> xdiff = x[1:] - x[:-1]
</code></pre>

<p>is much faster than</p>

<pre><code> xdiff = np.zeros(x.shape[0]-1)
 for i in range(x.shape[0]:
     xdiff[i] = x[i+1] - x[i]
</code></pre>

<p>The <code>np.nditer</code> isn't any faster.</p>

<p><code>nditer</code> is recommended as a general iteration tool in compiled code.  But its main value lies in handling broadcasting and coordinating the iteration over several arrays (input/output).  And you need to use buffering and <code>c</code> like code to get the best speed from <code>nditer</code> (I'll look up a recent SO question).</p>

<p><a href="http://stackoverflow.com/a/39058906/901925">http://stackoverflow.com/a/39058906/901925</a></p>

<p>Don't use <code>nditer</code> without studying the relevant <code>iteration</code> tutorial page (the one that ends with a <code>cython</code> example).</p>

<p>=========================</p>

<p>Just judging from experience, this approach will be fastest.  Yes it's going to iterate over <code>secaxis</code> a number of times, but those are all done in compiled code, and will be much faster than any iteration in Python.  And the <code>for f in X:</code> iteration is just a few times.</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: 
    res.append( (f[inds] + coefs*f[inds+1]) / (1+coefs) )
</code></pre>

<p><code>@HYRY</code> has explored alternatives for making the <code>where</code> step faster.  But as you can see the differences aren't that big.  Other possible tweaks</p>

<pre><code>inds1 = inds+1
coefs = -secaxis[inds] / secaxis[inds1]
coefs1 = coefs+1
for f in X:
    res.append(( f[inds] + coefs*f[inds1]) / coefs1)
</code></pre>

<p>If <code>X</code> was an array, <code>res</code> could be an array as well.</p>

<pre><code>res = (X[:,inds] + coefs*X[:,inds1])/coefs1
</code></pre>

<p>But for small <code>N</code> I suspect the list <code>res</code> is just as good.  Don't need to make the arrays any bigger than necessary.  The tweaks are minor, just trying to avoid recalculating things.</p>

<p>=================</p>

<p>This use of <code>np.where</code> is just <code>np.nonzero</code>.  That actually makes two passes of the array, once with <code>np.count_nonzero</code> to determine how many values it will return, and create the return structure (list of arrays of now known length).  And a second loop to fill in those indices.  So multiple iterations are fine if it keeps action simple.</p>

python arrays performance numpy itertools