Code Review Stack Exchange is a question and answer site for peer programmer code reviews. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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:

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

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

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)] = []

for iteration in xrange(total_iterations):

    # Loop over entire grid, using a random order for faster convergence

    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

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
        return 1-ALPHA
share|improve this question

migrated from Dec 5 '12 at 13:24

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

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
up vote 5 down vote accepted

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)

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

share|improve this answer
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

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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