Geodesic

# Geodesic for Last passage percolation
# author:: Alexandre THIERY
# date:: 18 Dec 2011
 
 
# LIBRARY IMPORT
from numpy import *
import pylab as pylab
 
#constants
c_up, c_down, c_right, c_left, c_still = 1, 2, 3, 4, 0
 
#create weights on each edge of Z^2
def create_weights(m,n):
    vertical,horizontal = random.exponential(1,(m,n)), random.exponential(1,(m,n))
    return(vertical,horizontal)
 
 
#compute the geodesic between the origin and (x,y)
def compute_geodesic(path,x,y):
    path_x, path_y = [x], [y]
    i,j = x,y
    while path[i,j] != c_still:
        if path[i,j] == c_up:
            i,j = i+1,j
        elif path[i,j] == c_down:
            i,j = i-1,j
        elif path[i,j] == c_right:
            i,j = i,j+1
        elif path[i,j] == c_left:
            i,j = i,j-1
        path_x.append(i), path_y.append(j)
    return (path_x, path_y)
 
 
 
#board = [-N,+N]^2
N = 100
vertical,horizontal = create_weights(2*N+1,2*N+1)
#geodesic[i,j] = geodesic distance between (i,j) and (0,0)
big_number = (10000 + N*N)
geodesic = ones((2*N+1,2*N+1))*big_number
#path = will store all the geodesics
path = zeros((2*N+1,2*N+1))
path[N,N] = c_still
geodesic[N,N] = 0
 
#update geodesic
iteration_kont = 0
while True:
    iteration_kont += 1
    print 'Iteration Number = ', iteration_kont
    n_changes = 0
    for i in range(1,2*N):
        for j in range(1,2*N):
            current = geodesic[i,j]
            g_right, g_left = geodesic[i,j+1], geodesic[i,j-1]
            g_up, g_down = geodesic[i+1,j], geodesic[i-1,j]
            p_up, p_down = vertical[i,j],vertical[i-1,j]
            p_left, p_right = horizontal[i,j-1],horizontal[i,j]
            if (g_left + p_left < current):
                current, path[i,j] = g_left + p_left, c_left
                geodesic[i,j] = current
                n_changes += 1
            if (g_right + p_right < current):
                current, path[i,j] = g_right + p_right, c_right
                geodesic[i,j] = current
                n_changes += 1
            if (g_down + p_down < current):
                current, path[i,j] = g_down + p_down, c_down
                geodesic[i,j] = current
                n_changes += 1
            if (g_up + p_up < current):
                current, path[i,j] = g_up + p_up, c_up
                geodesic[i,j] = current
                n_changes += 1
    if n_changes == 0: 
        break
 
#Draw all the geodesics
pylab.figure()
for i in arange(1,2*N):
    (path_x, path_y)=compute_geodesic(path,1,i)
    pylab.plot(path_x, path_y,'r-')
    (path_x, path_y)=compute_geodesic(path,2*N-1,i)
    pylab.plot(path_x, path_y,'r-')
    (path_x, path_y)=compute_geodesic(path,i,1)
    pylab.plot(path_x, path_y,'r-')
    (path_x, path_y)=compute_geodesic(path,i,2*N-1)
    pylab.plot(path_x, path_y,'r-')
pylab.plot([N],[N],'bo')
 
#Draw ball of radius alpha*N
pylab.figure()
alpha = 0.3
ball_x, ball_y = [N],[N]
for i in arange(1,2*N):
    for j in  arange(1,2*N):
        if geodesic[i,j] < alpha*N:
            ball_x.append(i), ball_y.append(j)
pylab.plot(ball_x, ball_y,'r.')
pylab.plot([N],[N],'bo')
pylab.show()