Python in Nashville

Yell123 at Yell123 at
Sat Jun 20 13:12:48 CEST 2009

I maintain a math hobby - magic square web page  

A very bright fellow in England sent me a algorithm written in python to  
solve a problem I was working on.  I can't seem to get it to work ...  syntax 
The segment below gives the explanation and the  code:
If you are interested in the problem this  solves ... see section 1 on my 
"Topographical Model For the Magic  Square"
Sorry about the delay  in replying; I've been having email problems.]

> I read a article by  your wife on magic squares. I have been fooling
> around with magic square  models. I read your bio and see you did some 
> with random walk  stuff. I have been wanting a algorithm ... ( perhaps a
> Greedy algorithm)  that would calculate the amount of retained water for 
> order of magic  square. ( see topographical model setion I on web link
> below)

Hm,  interesting problem. Let's see. (Please forgive me if some or all
of what I'm  about to say is totally familiar and/or obvious to you.)

If I've  understood you right, your problem is as follows. Call the number
written in  each cell its "height". Define the "water height" of a cell to 
the  smallest N such that there's a path from that cell to the outside of
the  square using only cells of height <= N. So, when you pour an  infinite
amount of water over the square and let it settle, each cell ends  up
with water up to its water height, if that's bigger than its height.  Then
the "capacity" of the square is the sum of (water height - height)  over
cells for which that's positive, and you want to calculate the  capacity
of a given square.

So, the following seems like a pretty good  algorithm. It's a bit like
Dijkstra's shortest-path algorithm. It needs a  data structure
called a "priority queue", which means a container that  stores
objects along with numbers called their "priorities", and  has
efficient ways of doing (1) add an object to the queue and (2)
pull  out the object with highest priority. We'll keep track of an
upper bound on  the water height for each cell, and use the queue
to store cells for which we  might be able to reduce the upper
bound. No, actually we'll use it to store  cells that might enable
us to reduce their neighbours' upper  bounds.

1. For each cell, set bound(cell) = n^2. (This is a very  boring
upper bound on the water height of the cell.)
2. For each cell on  the edge of the square:
2a. Set bound(cell) = height(cell).
2b. Put that  cell into the queue with priority -bound(cell).
(So lower cells have higher  priority. Actually, many priority queue
implementations already use the  "smaller number = higher priority"
3. While the queue isn't  empty:
3a. Pull out the highest-priority cell from the queue.
3b. For each  neighbour of that cell:
3b1. Let b =  max(bound(cell),height(neighbour)).
3b2. If if bound(neighbour) >  b:
3b2a. Set bound(neighbour) = b.
3b2b. Add the neighbour to the queue  with priority -b.

When this is finished, for every cell bound(cell)  should equal
(not just be an upper bound for) the cell's water height,  because:

1. bound(cell) is always an upper bound on the cell's water  height,
because the only way it gets decreased is when we find that the  cell
has a neighbour whose water height is known to be <= its new  value.
(So there's a path from the cell to the edge, using nothing  higher
than that value, beginning with that neighbour.)

2. For cells  along the edge, bound(cell) equals the cell's water height,
because for these  cells height = water height, and right at the start
we set bound =  height.

3. Suppose that when the algorithm finishes bound(cell) is  actually
bigger than cell's water height w. Then there is a path from the  cell
to the edge using nothing of height > w. Follow that path until you  find
a cell for which bound <= w. (There must always be one,  since
bound = height for edge cells and the last cell on our path is  an
edge cell and has height <= w.) Then we've found two adjacent
cells  of height <= w, one of which has bound <= w and one of
which as bound  > w. But at some point the former must have been
put into the queue with  bound <= w; when it came out, its neighbour
would have been examined and  had its bound decreased to <= w.

>From 1 and 3, we  find that at the end of the algorithm the bounds
equal the water  heights.

Now we can just sum max(bound-height,0) over all cells, and  get
the total capacity of the square.

Here's some Python code that  does it (rather inefficiently):

def process(square):
n =  len(square) # square is a list of lists
bounds = [n*[n*n] for i in  range(n)]
queue = [] # lousy implementation, just for algorithm  testing
# head of queue is last element
# big numbers for  high priority
for i in range(n-1):
bounds[i][0] =  square[i][0]
bounds[n-1][i] = square[n-1][i]
bounds[i+1][n-1] =  square[i+1][n-1]
bounds[0][i+1] =  square[0][i+1]
while queue:
(b,i,j) =  queue.pop()
for (di,dj) in  [(-1,0),(1,0),(0,-1),(0,1)]:
u=i+di; v=j+dj
if 0<=u<n and 0<=v<n:
b1 = max(square[u][v],-b)
if b1 <  bounds[u][v]:
bounds[u][v] = b1
return  bounds

def capacity(square):
bounds = process(square)
n = len(square)
s = 0
for i in range(n):
for j in range(n):
delta =  bounds[i][j]-square[i][j]
if delta>0: s +=  delta
return s

Testing on one of your  examples:

>>>  capacity([[6,26,49,34,35,13,12],[5,48,1,28,30,36,27],

(Your  page says 321, but I think 320 is correct. The sum of the heights
in the  submerged region is 394, not 393.)

How efficient is this? Well, it's not  quite obvious to me that a cell
can't be put into the queue more than once,  but I think that should
be rare. So let's suppose that typically each cell  goes in once. The
queue should never be bigger than the number of cells aside  from
complications that arise when a cell is inserted more than once;
it  certainly shouldn't be bigger than, say, 4 times the number of
cells. And  each priority queue operation should take no worse
than some modest multiple  of log(queue_size) basic operations.
So the total time to process the whole  square shouldn't be worse
than some small multiple of n^2 log n. That seems  like about the
best you could reasonably expect. (I have another approach in  mind
that might reduce that to n^2 f(n) where f grows even more  slowly
when n is very large, but at the cost of greater complexity when
n  is not very large. I bet you'll only be doing this for n not very  

Note: the idiotic "sort the list every time" priority queue  implementation
in the Python code above is emphatically *not* suitable for  production
use. Find a proper priority queue implementation in your language  of

If some one could help me get this working ... I would be  grateful. 
I have loaded Python Ver 2.6 ... but when I paste this code in it  keeps 
giving me errors. 
Craig Knecht

**************A Good Credit Score is 700 or Above. See yours in just 2 easy 
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <>

More information about the Python-list mailing list