3

I have to iterate over two arrays which are 1000x1000 big. I already reduced the resolution to 100x100 to make the iteration faster, but it still takes about 15 minutes for ONE array! So I tried to iterate over both at the same time, for which I found this:

for index, (x,y) in ndenumerate(izip(x_array,y_array)):

but then I get the error:

ValueError: too many values to unpack

Here is my full python code: I hope you can help me make this a lot faster, because this is for my master thesis and in the end I have to run it about a 100 times...

area_length=11
d_circle=(area_length-1)/2

xdis_new=xdis.copy()
ydis_new=ydis.copy()
ie,je=xdis_new.shape


while (np.isnan(np.sum(xdis_new))) and (np.isnan(np.sum(ydis_new))):
xdis_interpolated=xdis_new.copy()
ydis_interpolated=ydis_new.copy()
# itx=np.nditer(xdis_new,flags=['multi_index'])
# for x in itx:
    # print 'next x and y'
for index, (x,y) in ndenumerate(izip(xdis_new,ydis_new)):
    if np.isnan(x):
        print 'index',index[0],index[1]
        print 'interpolate'
        # define indizes of interpolation area
        i1=index[0]-(area_length-1)/2
        if i1<0:
            i1=0
        i2=index[0]+((area_length+1)/2)
        if i2>ie:
            i2=ie
        j1=index[1]-(area_length-1)/2
        if j1<0:
            j1=0
        j2=index[1]+((area_length+1)/2)
        if j2>je:
            j2=je
        # -->
        print 'i1',i1,'','i2',i2
        print 'j1',j1,'','j2',j2

        area_values=xdis_new[i1:i2,j1:j2]
        print area_values

        b=area_values[~np.isnan(area_values)]

        if len(b)>=((area_length-1)/2)*4:

            xi,yi=meshgrid(arange(len(area_values[0,:])),arange(len(area_values[:,0])))

            weight=zeros((len(area_values[0,:]),len(area_values[:,0])))
            d=zeros((len(area_values[0,:]),len(area_values[:,0])))
            weight_fac=zeros((len(area_values[0,:]),len(area_values[:,0])))
            weighted_area=zeros((len(area_values[0,:]),len(area_values[:,0])))

            d=sqrt((xi-xi[(area_length-1)/2,(area_length-1)/2])*(xi-xi[(area_length-1)/2,(area_length-1)/2])+(yi-yi[(area_length-1)/2,(area_length-1)/2])*(yi-yi[(area_length-1)/2,(area_length-1)/2]))
            weight=1/d
            weight[where(d==0)]=0
            weight[where(d>d_circle)]=0
            weight[where(np.isnan(area_values))]=0

            weight_sum=np.sum(weight.flatten())
            weight_fac=weight/weight_sum
            weighted_area=area_values*weight_fac

            print 'weight'
            print weight_fac
            print 'values'
            print area_values
            print 'weighted'
            print weighted_area

            m=nansum(weighted_area)
            xdis_interpolated[index]=m
            print 'm',m

        else:
            print 'insufficient elements'

    if np.isnan(y):
        print 'index',index[0],index[1]
        print 'interpolate'
        # define indizes of interpolation area
        i1=index[0]-(area_length-1)/2
        if i1<0:
            i1=0
        i2=index[0]+((area_length+1)/2)
        if i2>ie:
            i2=ie
        j1=index[1]-(area_length-1)/2
        if j1<0:
            j1=0
        j2=index[1]+((area_length+1)/2)
        if j2>je:
            j2=je
        # -->
        print 'i1',i1,'','i2',i2
        print 'j1',j1,'','j2',j2

        area_values=ydis_new[i1:i2,j1:j2]
        print area_values

        b=area_values[~np.isnan(area_values)]

        if len(b)>=((area_length-1)/2)*4:

            xi,yi=meshgrid(arange(len(area_values[0,:])),arange(len(area_values[:,0])))

            weight=zeros((len(area_values[0,:]),len(area_values[:,0])))
            d=zeros((len(area_values[0,:]),len(area_values[:,0])))
            weight_fac=zeros((len(area_values[0,:]),len(area_values[:,0])))
            weighted_area=zeros((len(area_values[0,:]),len(area_values[:,0])))

            d=sqrt((xi-xi[(area_length-1)/2,(area_length-1)/2])*(xi-xi[(area_length-1)/2,(area_length-1)/2])+(yi-yi[(area_length-1)/2,(area_length-1)/2])*(yi-yi[(area_length-1)/2,(area_length-1)/2]))
            weight=1/d
            weight[where(d==0)]=0
            weight[where(d>d_circle)]=0
            weight[where(np.isnan(area_values))]=0

            weight_sum=np.sum(weight.flatten())
            weight_fac=weight/weight_sum
            weighted_area=area_values*weight_fac

            print 'weight'
            print weight_fac
            print 'values'
            print area_values
            print 'weighted'
            print weighted_area

            m=nansum(weighted_area)
            ydis_interpolated[index]=m
            print 'm',m

        else:
            print 'insufficient elements'

    else:
        print 'no need to interpolate'

xdis_new=xdis_interpolated
ydis_new=ydis_interpolated
Simeon Visser
  • 118,920
  • 18
  • 185
  • 180
Melanie Maza
  • 207
  • 4
  • 13
  • For those not in the know - izip() is from itertools module: `from itertools import izp` – Pierz Sep 10 '13 at 16:34

5 Answers5

2

Some advice:

  • Profile your code to see what is the slowest part. It may not be the iteration but the computations that need to be done each time.
  • Reduce function calls as much as possible. Function calls are not for free in Python.
  • Rewrite the slowest part as a C extension and then call that C function in your Python code (see Extending and Embedding the Python interpreter).
  • This page has some good advice as well.
Simeon Visser
  • 118,920
  • 18
  • 185
  • 180
2

You specifically asked for iterating two arrays in a single loop. Here is a way to do that

l1 = ["abc", "def", "hi"]
l2 = ["ghi", "jkl", "lst"]
for f,s in zip(l1,l2):
    print "%s : %s" %(f,s)

The above is for python 3, you can use izip for python 2

fkl
  • 5,412
  • 4
  • 28
  • 68
1

You may use this as your for loop:

for index, x in ndenumerate((x_array,y_array)):

But it wont help you much, because your computer cant do two things at the same time.

MaxPowers
  • 5,235
  • 2
  • 44
  • 69
0

Comment #1: You don't want to use ndenumerate on the izip iterator, as it'll output you the iterator, which isn't what you want.

Comment #2:

i1=index[0]-(area_length-1)/2
if i1<0:
    i1=0

could be simplified in i1 = min(index[0]-(area_length-1)/2, 0), and you could store your (area_length+/-1)/2 in specific variables.

Idea #1 : try to iterate on flat versions of the arrays, i.e. with something like

for (i, (x, y)) in enumerate(izip(xdis_new.flat,ydis_new.flat)):

You could get the original indices via divmod(i, xdis_new.shape[-1]), as you should be iterating by rows first.

Idea #2 : Iterate only on the nans, i.e. indexing your arrays with np.isnan(xdis_new)|np.isnan(ydis_new), that could save you some iterations

EDIT #1

  • You probably don't need to initialize d, weight_fac and weighted_area in your loop, as you compute them separately.

  • Your weight[where(d>0)] can be simplified in weight[d>0]

  • Do you need weight_fac ? Can't you just compute weight then normalize it in place ? That should save you some temporary arrays.

Pierre GM
  • 19,809
  • 3
  • 56
  • 67
  • Your comments and ideas sound very good! I already thought about flatten the arrays but then I had exactly the problem with getting the original indices, so thank you very much! I will try what you said and then we'll see... =) – Melanie Maza Aug 14 '12 at 09:11
  • I have a problem with the divmod-thing...I put it after the for loop, but I get the error: TypeError: unsupported operand type(s) for divmod(): 'tuple' and 'int' Could you please tell how and where exactly use this function? Thank you! – Melanie Maza Aug 15 '12 at 08:36
  • Ok, I think I have an answer for the divmod func... [link] (http://stackoverflow.com/questions/9482550/argmax-of-numpy-array-returning-non-flat-indices) it seems to work... – Melanie Maza Aug 15 '12 at 09:06
  • You could try something like: `[(i,j,x,y) for ((i,j),x,y) in zip(zip(*np.unravel_index(np.arange(np.multiply(*X.shape)),X.shape)),X.flat,Y.flat)]` – Pierre GM Aug 16 '12 at 09:50
0

Profiling is definitely a good start to identify where all the time spent actually goes.

I usually use the cProfile module, as it requires minimal overhead and gives me more than enough information.

import cProfile
import pstats
cProfile.run('main()', "ProfileData.txt", 'tottime')
p = pstats.Stats('ProfileData.txt')   
p.sort_stats('cumulative').print_stats(100)

I your example you would have to wrap your code into a main() function to be able to use this code snippet at the very end of your file.

Rostyslav Dzinko
  • 39,424
  • 5
  • 49
  • 62
Peter Lustig
  • 941
  • 11
  • 23