[Numpy-discussion] when and where to use numpy arrays vs nested lists

[Mark P. Miller]

Charles R Harris wrote:
> I won't comment on the code itself. Tell us what you want to do and I 
> bet we can speed it up.
> 


Here's a bit of sample code.  It's actually very useful for my purposes. 
  Are there good ways to express these ideas in numpy as opposed to 
using all of the nested loops that I currently employ?

#########
import random
import numpy
uniformRND = random.random #faster than numpy's
normalRND = numpy.random.normal #faster than Python's


def pick(a,b,d1,d2,sd):
     # this function picks two random values with mean = 0 and
     # sd = sd.  Random values are added (in one form or another)
     # to a and b.  All other code ensures that the random values
     # do not generate "index out of range" errors after new1 and
     # new2 are returned (needed when sd is set to large values)

     new1 = 2*d1 + 1
     while new1 < -d1 or new1 > 2*d1:
         new1=int(round(normalRND(0,sd)))

     new2 = 2*d2 + 1
     while new2 < -d2 or new2 > 2*d2:
         new2=int(round(normalRND(0,sd)))

     tmp1 = a + new1
     if 0 <= tmp1 < d1:
         new1 = tmp1
     elif tmp1 < 0:
         new1 = -(a + new1)
     else:
         new1 = a - new1

     tmp2 = b + new2
     if 0 <= tmp2 < d2:
         new2 = tmp2
     elif tmp2 < 0:
         new2 = -(b + new2)
     else:
         new2 = b - new2

     return new1, new2

def main():
     changerate = 0.25
     sd = 2
     d1=50
     d2=100
     d3=10
     ncycles = 10
     array1 = numpy.zeros((d1, d2, d3, 2), int)
     array2 = numpy.zeros((d1, d2, d3, 2), int)

     for zz in xrange(ncycles):
         for aa in xrange(d1):
             for bb in xrange(d2):
                 item1a, item1b = pick(aa,bb,d1,d2,sd)
                 item2a, item2b = pick(aa,bb,d1,d2,sd)

                 for cc in xrange(d3):
                     element1 = array1[item1a,item1b,cc,
			              int(uniformRND() < 0.5)]
                     element2 = array1[item2a,item2b,cc,
				      int(uniformRND() < 0.5)]

                     if uniformRND() < changerate:
                         element1 = int(uniformRND()*100)

                     if uniformRND() < changerate:
                         element2 = int(uniformRND()*100)

                     array2[aa,bb,cc,0]=element1
                     array2[aa,bb,cc,1]=element2

         array1=array2.copy()

if __name__ == '__main__':
     main()