Using Gibbs sampling to segment an image

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.

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)

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

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

-
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.])

-
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.])

-
You can then multiply those values by other values:

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

-
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.

StackExchange Code Review Q#19329, answer score: 6