#!/usr/bin/env python
from scipy import *
from pylab import *

def Distance2(R1,R2):
    return (R1[0]-R2[0])**2+(R1[1]-R2[1])**2

def TotalDistance2(city, R):
    dist = 0
    for i in range(len(city)-1):
        dist += Distance2(R[city[i]], R[city[i+1]])
    dist + Distance2(R[city[-1]],R[city[0]])
    return dist

def Plot(city, R, dist):
    # Plot
    Pt = [R[city[i]] for i in range(len(city))]
    Pt += [R[city[0]]]
    Pt = array(Pt)
    title('Total distance='+str(dist))
    plot(Pt[:,0], Pt[:,1], '-o')
    show()

    
if __name__ == '__main__':

    ncity = 100
    maxTsteps = 100
    Tstart = 0.2
    fCool = 0.9
    maxSteps= 100*ncity
    maxAccepted = 10*ncity

    Preverse = 0.5

    R=[]
    for i in range(ncity):
        R.append( [rand(),rand()] )
    R = array(R)

    # The index table
    city = range(ncity)
    dist = TotalDistance2(city, R)

    n = zeros(6, dtype=int)
    nct = len(city)
    
    T = Tstart
    Plot(city, R, dist)

    accepted = 0

    for i in range(maxSteps):

        while True: # Will find two random cities sufficiently close by
            # Two cities n[0] and n[1] are choosen at random
            n[0] = int(nct*rand())      # select one city
            n[1] = int((nct-1)*rand())  # select another city, but not the same
            if (n[1]>=n[0]): n[1] += 1 
            if (n[1] < n[0]): (n[0],n[1]) = (n[1],n[0]) # swap, because it is convenient to have: n[0]<n[1]
            nn = (n[0]+nct - n[1]-1) % nct # number of cities not on the segment n[0]..n[1]
            if nn>=3: break

        n[2] = (n[0]-1) % nct
        n[3] = (n[1]+1) % nct