# How can I make this code faster using Numpy?

I have been reading some Numpy guides but can't seem to figure it out. my TA told me I should be able to speed up my code by using a numpy array instead of a for loop in the following segment of code.

for neighbor in get_neighbors(estimates,i,j):
pXgivenX_ *= edge_model(True,neighbor)*observation_model(obs,True)
pX_givenX_  *= edge_model(False,neighbor)*observation_model(obs,False)


here is the entire code of the method it is in if it helps:

def gibbs_segmentation(image, burnin, collect_frequency, n_samples):
"""
Uses Gibbs sampling to segment an image into foreground and background.

Inputs
------
image : a numpy array with the image. Should be Nx x Ny x 3
burnin : Number of iterations to run as 'burn-in' before collecting data
collect_frequency : How many samples in between collected samples
n_samples : how many samples to collect in total

Returns
-------
A distribution of the collected samples: a numpy array with a value between
0 and 1 (inclusive) at every pixel.
"""
(Nx, Ny, _) = image.shape
total_iterations = burnin + (collect_frequency * (n_samples - 1))
pixel_indices = list(itertools.product(xrange(Nx),xrange(Ny)))
# The distribution that you will return
distribution = np.zeros( (Nx, Ny) )

# Initialize binary estimates at every pixel randomly. Your code should
# update this array pixel by pixel at each iteration.
estimates = np.random.random( (Nx, Ny) ) > .5

# PreProcessing
preProObs = {}
for (i,j) in pixel_indices:
preProObs[(i,j)] = []
preProObs[(i,j)].append(observation_model(image[i][j],False))
preProObs[(i,j)].append(observation_model(image[i][j],True))

for iteration in xrange(total_iterations):

# Loop over entire grid, using a random order for faster convergence
random.shuffle(pixel_indices)

for (i,j) in pixel_indices:

pXgivenX_ = 1
pX_givenX_ = 1

for neighbor in get_neighbors(estimates,i,j):
pXgivenX_ *= edge_model(True,neighbor)*preProObs[(i,j)][1]
pX_givenX_  *= edge_model(False,neighbor)*preProObs[(i,j)][0]
estimates[i][j] = np.random.random() > pXgivenX_/(pXgivenX_+pX_givenX_)

if iteration > burnin and (iteration-burnin)%collect_frequency == 0:
distribution += estimates

return distribution / n_samples


Basically my TA hinted that I could do it without a for loop using numpy arrays and that it would be faster.

Edit-- I forgot to include the important part:

def edge_model(label1, label2):
"""
Given the values at two pixels, returns the edge potential between
those two pixels.

Hint: there might be a more efficient way to compute this for an array
of values using numpy!
"""
if label1 == label2:
return ALPHA
else:
return 1-ALPHA

-
What is neighbor? What does observation_model(obs,True) return? –  Sam Mussmann Dec 4 '12 at 3:22
neighbor is a boolean corresponding to whether the pixel is currently believed to be in the foreground or background. observation_model(obs,True) returns a number corresponding to the probability of observing True given that you are at obs (where True corresponds to the pixel being in the foreground). –  Coolguy123 Dec 4 '12 at 3:29
add comment

## migrated from stackoverflow.comDec 5 '12 at 13:24

This question came from our site for professional and enthusiast programmers.

## 1 Answer

Since this is homework, I won't give you the exact answer, but here are some hints.

1. == is overloaded in numpy to return an array when you pass in an array. So you can do things like this:

>>> numpy.arange(5) == 3
array([False, False, False,  True, False], dtype=bool)
>>> (numpy.arange(5) == 3) == False
array([ True,  True,  True, False,  True], dtype=bool)

2. You can use boolean arrays to assign to specific locations in an array. For example:

>>> mostly_true = (numpy.arange(5) == 3) == False
>>> empty = numpy.zeros(5)
>>> empty[mostly_true] = 5
>>> empty
array([ 5.,  5.,  5.,  0.,  5.])

3. You can also negate boolean arrays; together, these facts allow you to conditionally assign values to an array. (numpy.where can be used to do something similar.):

>>> empty[~mostly_true] = 1
>>> empty
array([ 5.,  5.,  5.,  1.,  5.])

4. You can then multiply those values by other values:

>>> empty * numpy.arange(5)
array([  0.,   5.,  10.,   3.,  20.])

5. And many different numpy functions (really ufuncs) provide a reduce method that applies the function along the entire array:

>>> results = empty * numpy.arange(5)
>>> numpy.multiply.reduce(results)
0.0


You should be able to completely eliminate that for loop using only the above techniques.

-
If it is True then I want to multiply by ALPHA and if it is false I multiply by 1-ALPHA but how can I do this in one line? –  Coolguy123 Dec 4 '12 at 3:53
Nevermind I looked at np.where –  Coolguy123 Dec 4 '12 at 3:54
add comment