2
\$\begingroup\$

I have a function written in Python 2.7 that I call several times in a program. This function is rather slow, and I'd like to rewrite it in NumPy. However, I don't know how to do it.

I basically have a canvas of width * height pixels, and I'm making something rather complicated inside. I'd like to know how to replace all the parts with:

for i in range(x):
    do x[i]
    ...

by NumPy-style operations.

I know that I can start by making an array of zeros with np.zeros, but how can I implement the equations / operations in NumPy?

def Circular_cumsum_spot(width,height,centre,radius,spectrum,red_or_violet):
    """makes a "light spot" by cumsums
    from the center to the periphery

    width, height = size of canvas
    x, y = center of the spot
    radius = magnitude of the spot in pixels
    spectrum = light source
    red = to red, violet = to violet

    returns a numpy array"""

    x,y = centre[0], centre[1]

    data = []
    X = Y = Z = 0

    for i in range(height):
        for j in range(width):

            distance = radius - np.sqrt((i-x)**2+(j-y)**2)

            if distance >= radius:
                cumsum = 0
                X = Y = Z = 0
            else:
                cumsum = round((distance/radius) * (len(spectrum)))

                if red_or_violet == "red":


                    for k in range(cumsum):

                        X += spectrum[k][0] 
                        Y += spectrum[k][1]
                        Z += spectrum[k][2]

                else:

                    for k in range(cumsum):

                        l = len(spectrum)-k-1

                        X += spectrum[l][0] 
                        Y += spectrum[l][1] 
                        Z += spectrum[l][2] 

            data += [(X,Y,Z)]
            X = Y = Z = 0
    data = np.array(data)

    divisor = np.max(data[:,1])

    if divisor == 0:
        data = data
    else:
        data = np.divide(data, divisor)

    return data
    print "spot created is", red_or_violet
\$\endgroup\$
  • 1
    \$\begingroup\$ Your code has a lot of loops at the Python level. With some hard work trying to convert the loops into ufunc numpy calls, you could probably achieve a few multiples faster. However, if you convert this code to Cython, and set types on your variables, you can realistically expect to get it around 150X faster (15000% faster). The loops should all be explicit Cython loops. \$\endgroup\$ – Caleb Hattingh Oct 14 '14 at 5:47
  • \$\begingroup\$ That sounds great, how do I use Cython? I never heard of it, thanks \$\endgroup\$ – adrienlucca.wordpress.com Oct 14 '14 at 8:51
  • 1
    \$\begingroup\$ docs.cython.org \$\endgroup\$ – Caleb Hattingh Oct 14 '14 at 21:39
2
\$\begingroup\$

don't know how to create the XYZ list within numpy, but some other parts can be converted to numpy (see below). Didn't test the code, but the help pages for those numpy functions should help you, in case this doesn't run ;) Instead of the double for-loop one can use np.nditer which might speed things up, but I haven't really used this much myself.

Otherwise, you want to remove for loops and just operate on the whole matrix or vector using numpy, e.g. calculate the distances in a single line instead of looping over it using two python loops. Using the axis=0 you can also remove the for loop for the spectrum variable.

def Circular_cumsum_spot(width,height,centre,radius,spectrum,red_or_violet):
    """makes a "light spot" by cumsums
    from the center to the periphery

    width, height = size of canvas
    x, y = center of the spot
    radius = magnitude of the spot in pixels
    spectrum = light source
    red = to red, violet = to violet

    returns a numpy array"""

    x,y = centre[0], centre[1]

    data = []
    X = Y = Z = 0

    # calc an 2d array of distances
    Xi, Yj = np.mgrid(0:height, 0:width)
    distance = np.hypot(Xi-y, Yj-x)

    # create a mask
    mask = distance < radius

    # calculate all cumsum values
    cumsum = np.zeros_like(distance)
    cumsum[mask] = np.round(distance[mask]/radius)*len(spectrum)

    # not sure, if this is already a numpy array
    spectrum = np.array(spectrum)

    for i in range(height):
        for j in range(width):
            if distance[i, j] >= radius:
                X = Y = Z = 0
            else:
                if red_or_violet == "red":
                    X, Y, Z = np.sum(spectrum[:cumsum[i, j],:], axis=0)
                else:  
                    X, Y, Z = np.sum(spectrum[-cumsum[i, j]:,:], axis=0)
            data += [(X,Y,Z)]
    data = np.array(data)

    divisor = data[:,1].max()

    if divisor != 0:
        data = data/divisor

    return data
    print "spot created is", red_or_violet
\$\endgroup\$
  • \$\begingroup\$ Hi, thanks for the effort! There are many useful suggestions, even if the format is not what the CR site's expecting... I get a weird error though, cumsum[mask] = np.round((distance/radius)*len(spectrum)) returns: ValueError: boolean index array should have 1 dimension. I quite don't get what you did there, any idea? \$\endgroup\$ – adrienlucca.wordpress.com Dec 1 '14 at 16:23
  • \$\begingroup\$ mask will be an array of the same length as distance. For each element in distance it will contain either true or false depending on the distance[i] < radius condition. By using cumsum[mask] you only set those values that have a true in the mask. I guess There was an extra ")" in the code above. Will fix it. Don't have time to test it though :( \$\endgroup\$ – Arun Dec 2 '14 at 6:02
  • \$\begingroup\$ Hmm, did test a bit afterall and found some other problems... Edited the code a bit more. The problem you ran into was, because on the left side, I used the mask, and on the right side I used the full array, so the sizes of those arrays where different. At least that's my best guess ;) I think I also had x and y flipped when creating the distance array and there was probably an extra radius that shouldn't be there? \$\endgroup\$ – Arun Dec 2 '14 at 6:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.