title: Can I avoid Python loop overhead on dynamic programming with numpy?
body: <p>I need help with the Pythonic looping overhead of the following problem:  I'm writing a function that calculates a pixel flow algorithm that's a classic dynamic programming algorithm on a 2D Numpy array.  It requires:</p>
  
  <p>1) visiting all the elements of the array at least once like this:</p>
  
  <pre><code>for x in range(xsize):
      for y in range(ysize):
           updateDistance(x,y)
  </code></pre>
  
  <p>2) potentially following a path of elements based on the values of the neighbors of an element which looks like this</p>
  
  <pre><code>while len(workingList) &gt; 0:
     x,y = workingList.pop()
     #if any neighbors of x,y need calculation, push x,y and neighbors on workingList 
     #else, calculate flow on pixels as a sum of flow on neighboring pixels
  </code></pre>
  
  <p>Unfortunately, I seem to be getting a lot of Pythonic loop overhead on #1 even if the call to updateDistance is <code>pass</code>.  I figure this is a classic enough algorithm that there must be a good approach to it in Python that can avoid some of the looping overhead.  I'm also worried that if I can fix #1 I'm going to get busted on the loop in #2.</p>
  
  <p>Any suggestions about quickly looping through elements in a 2D numpy array and potentially updating chains of elements in that array?</p>
  
  <p><strong>Edit: Flushing out more details of #2</strong></p>
  
  <p>It seems I might be able to vectorize the first loop, perhaps with vectorizing an  <code>np.meshgrid</code> call.</p>
  
  <p>The loop in part is a little complicated, but here's a simplified version.  I'm concerned about both the loop and the indexing into the neighboring elements:</p>
  
  <pre><code>#A is a 2d cost matrix
  workingList = [(x,y)]
  while len(workingList) &gt; 0:
     x,y = workingList.pop()
     neighborsToCalculate = []
     for n in neighborsThatNeedCalculation(x,y): #indexes A to check neighbors of (x,y)
        neighborsToCalculate.append(n)
     if len(neighborstToCalculate) != 0:
         workingList.append((x,y))
         workingList.extend(neighborsToCalculate)
     else:
         for xn,yn in neighbors(x,y):
            A[x,y] += 1+A[xn,yn]
  </code></pre>
  
  <p>This is a classic breadth first search problem.  It would be great if this could be parallelized.  It probably can't in its current form because it follows a path, but I'd be delighted for suggestions.</p>