Python, how to optimize this code

Question

I tried to optimize the code below but I cannot figure out how to improve computation speed. I tried Cthon but the performance is like in python.

Is it possible to improve the performance without rewrite everything in C/C++?

Thanks for any help

import numpy as np

heightSequence = 400
widthSequence = 400
nHeights = 80

DOF = np.zeros((heightSequence, widthSequence), dtype = np.float64)
contrast = np.float64(np.random.rand(heightSequence, widthSequence, nHeights))

initDOF = np.zeros([heightSequence, widthSequence], dtype = np.float64)
initContrast = np.zeros([heightSequence, widthSequence, nHeights], dtype = np.float64)
initHeight = np.float64(np.r_[0:nHeights:1.0])
initPixelContrast = np.array(([0 for ii in range(nHeights)]), dtype = np.float64)


# for each row
for row in range(heightSequence):
    # for each col
    for col in range(widthSequence):

        # initialize variables            
        height = initHeight # array ndim = 1
        c = initPixelContrast # array ndim = 1

        # for each height            
        for indexHeight in range(0, nHeights):
            # get contrast profile for current pixel
            tempC = contrast[:, :, indexHeight]
            c[indexHeight] = tempC[row, col]

        # save original contrast            
        # originalC = c
        # originalHeight = height                

        # remove profile before maximum and after minumum contrast
        idxMaxContrast = np.argmax(c)
        c = c[idxMaxContrast:]
        height = height[idxMaxContrast:]

        idxMinContrast = np.argmin(c) + 1
        c = c[0:idxMinContrast]
        height = height[0:idxMinContrast]              

        # remove some refraction
        if (len(c) <= 1) | (np.max(c) <= 0):
            DOF[row, col] = 0                  

        else:

            # linear fitting of profile contrast                                             
            P = np.polyfit(height, c, 1)
            m = P[0]
            q = P[1]

            # remove some refraction               
            if m >= 0:
                DOF[row, col] = 0

            else:
                DOF[row, col] = -q / m

    print 'row=%i/%i' %(row, heightSequence)

# set range of DOF
DOF[DOF < 0] = 0
DOF[DOF > nHeights] = 0

What's the data? Can you add some for reference so that this piece of code runs by itself? You could generate some arrays of the right shape with `np.random.rand` , for example. — YXD, Feb 28 '13 at 11:32
If you **do not** want C or C++ suggestions then do not add such tags. — Öö Tiib, Feb 28 '13 at 11:33
The data is a stack of 80 images... I will update the code to run by itself but I need some time to test if everything is OK xD this is only a small portion of my whole algorithm where is the bottleneck — opensw, Feb 28 '13 at 11:46
That would be good - otherwise it's all guesswork on the part of anyone trying to answer your question and that's not much use to you or us... — YXD, Feb 28 '13 at 11:50
[*Possibly helpful ...*](http://stackoverflow.com/a/4299378/23771) — Mike Dunlavey, Feb 28 '13 at 12:42
I updated the code, now it is running by itself :) on i7 CPU it is very slow... and this is only for one stack of images, usually I have 260 stacks and the resolution is much greater than in this demo xD — opensw, Feb 28 '13 at 13:52

score 4 · Accepted Answer · answered Feb 28 '13 at 11:53

4

By looking at the code it seems that you can get rid of the two outer loops completely, converting the code to a vectorised form. However, the np.polyfit call must then be replaced by some other expression, but the coefficients for a linear fit are easy to find, also in vectorised form. The last if-else can then be turned into a np.where call.

answered Feb 28 '13 at 11:53

Michael Wild

24,977
3
43
43

Thanks for the tip, I will try to follow your advice... but I am pretty new in python, anyway, I will try it :) – opensw Feb 28 '13 at 14:01
1

It's the same advice that applies to e.g. using Matlab: try to get rid of the loops (at least most of the time) – Michael Wild Feb 28 '13 at 14:04
OK, I finished to implement the part about polyfit algorithm (I implemented these formulas http://www.codecogs.com/code/maths/approximation/regression/linear.php in python in very efficient way without the "for row" and "for col" loops). My last problem is to optimize the statements with argmax(c) and argmin(c). Have you an idea how to do this computation outside the loops? – opensw Feb 28 '13 at 21:58
1

Just specify the `axis` along which you want to operate. If you don't specify an `axis`, `numpy` operates on the flattened array (i.e. the global minimum and maximum). In your case, `axis` would be probably be `2`. – Michael Wild Mar 01 '13 at 06:48
Dear Michael Wild, thanks a lot :) it's working like a charm :) now I can analyze one experiment in less than 1 minute, I really appreciate your help – opensw Mar 01 '13 at 11:32

Python, how to optimize this code

1 Answers1