I maintain a math hobby - magic square web page
_http://www.knechtmagicsquare.paulscomputing.com/_
(http://www.knechtmagicsquare.paulscomputing.com/)
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
problems
The segment below gives the explanation and the code:
If you are interested in the problem this solves ... see section 1 on my
website
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
work
with random walk stuff. I have been wanting a algorithm ... ( perhaps a
Greedy algorithm) that would calculate the amount of retained water for
any
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
be
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
convention.)
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.
Contradiction.
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]
queue.append((-bounds[i][0],i,0))
bounds[n-1][i] = square[n-1][i]