#!/usr/bin/env python
# coding: utf-8

# # Simulated annealing & Traveling Salesman problem

# Simulated annealing is a technique which uses Matropolis to find a global minimum of a
# function of many variables. It is typically used for large scale problems, especially ones where a desired global minimum is hidden among many poorer local minima.
# 
# The class of problems known as NP-complete problems, whose computation time for an exact solution increases with N as $\exp(const.× N)$, were basically unsolvable until recently.
# 
# Only in 1983, the method was invented ( Kirkpatrick, S., Gelatt, C.D., and Vecchi, M.P. 1983, Science, vol. \textbf{220}, 671-680 (1083).) which effectively solved most of these problems for practical purposes.
# 
# Simulated annealing does not "solve" them exactly, but in most cases we just need reasonable good solution even if its not the one with exacly the lowest total energy.
# 
# The archtype example of NP-complete problems, the traveling salesman problem, is easy to solve with this method.
# 
# Other very important application is in designing complex integrated circuits: The arrangement of several hundred thousand circuit elements on a tiny silicon substrate is optimized so as to minimize interference among their connecting wires.
# 
# 
# The simulated annealing uses Metropolis algorithm with slowly decreasing temperature in order that system relaxes into its "ground state".
# 
# The steps consis of
# - Pick initial configuration $X$ and an initial temperature
#   $T_0$. $T_0$ should be much higher than the changes in function $f$
#   to be minimized when typical Monte Carlo step is taken.
# - Loop through decreasing temperature $T_i$
# 
#   - Equilibrate at $T_i$ using Metropolis with selected, allowed
#   elementary changes in the system. The system is equilibrated when
#   the rate of change of $f$ averaged over some number of Monte Carlo
#   steps is small.
#   - Measure the termal average of $f$. If $f$ does not decrease at
#   this temperature $T_i$ compared to previous temperature $T_{i-1}$,
#   exit the loop.
#   - Decrease $T_i$ to $T_{i+1}$. For example, reduces the effective temperature by 10\%.
# 
# 
# The major difficulty (art) in implementation of the algorithm is that
# there is no obvious analogy for the temperature T with respect to a
# free parameter in the combinatorial problem. Furthermore, avoidance of
# entrainment in local minima (quenching) is dependent on the "annealing
# schedule", the choice of initial temperature, how many iterations are
# performed at each temperature, and how much the temperature is
# decremented at each step as cooling proceeds.
# 

# ## The Traveling Salesman problem
# 
# The seller visits $N$ cities ($i=0...N-1$) with given positions
# $R_i$, returning finally to his or her city of origin. 
# **Each city** is to be **visited only once**, and the **route** is to be made **as
# short as possible**. This problem belongs to a class known as
# *NP-complete* problems, whose computation time for an exact solution
# increases with N exponentially.
# 
# The traveling salesman problem also belongs to a class of minimization
# problems for which the objective function E has many local minima. In
# practical cases, it is often enough to be able to choose from these a
# minimum which, even if not absolute, cannot be significantly improved
# upon.
# 
# The annealing method manages to achieve this, while limiting its
# calculations to scale as a small power of $N$.
# 
# 
# As a problem in simulated annealing, the traveling salesman problem is handled
# as follows:
# 
# - 1 **Configuration:** The cities are numbered $i = 0. . .N-1$ and each has coordinate $R_i$. A configuration is a permutation of the number $0 ...N-1$, interpreted as the order in which the cities are visited.
# - 2 **Rearrangements:** An efficient set of moves are:
#     - **(a)** A section of path is removed and then replaced with the same cities running in   the opposite order; or
#     - **(b)** A section of path is removed and then replaced in between two cities on another, randomly chosen, part of the path.
# 
# - 3 **Objective Function:** In the simplest form of the problem, $E$ is taken just as the total length of journey
# 
# 
# The above two mentioned moves are hardest to implement.  The following
# plot explains them by an example
# 
# <img src="reverse.png" width="50%">
# <img src="transportn.png" width="50%">
# 
# For details of implementation, check the example source code.
# 

# In[16]:


import scipy as scp
from numpy import random
from numpy import linalg
from numba import jit

ncity=100

# random coordinates in 2D for n-cities
R = random.random((ncity,2)) # completely randomly positioned cities
city = range(ncity)          # we are visiting them in order in the array


# In[3]:


def Distance(R1,R2):
    return linalg.norm(R1-R2) # sqrt(x^2+y^2)

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


# In[17]:


def Plot(city, R, dist):
    Pt = [R[city[i]] for i in range(len(city))] # list of [(x_i,y_i),...] coordinates
    Pt += [R[city[0]]]   # add the very last/first city so that the curve is connected to itself
    Pt = array(Pt)
    title('Total distance='+str(dist))
    plot(Pt[:,0],Pt[:,1],'o-')
    show()


# In[18]:


from pylab import *
get_ipython().run_line_magic('matplotlib', 'inline')
Plot(city,R, TotalDistance(city,R))


# In[19]:


@jit(nopython=True)
def FindASegment(R):
    nct = len(R) # number of cities
    while True:
        # Two cities n[0] and n[1] chosen at random
        n0 = int(nct*rand())
        n1 = int((nct-1)*rand())
        if n1>=n0 : n1 +=1
        if n1<n0 : (n0,n1) = (n1,n0) # always make sure that n1>n0
        nn = (nct-(n1-n0+1))         # the rest of the cities
        if nn>=3 : break             # the rest is big enough
    n2 = (n0-1) % nct        # the city before first
    n3 = (n1+1) % nct        # the city past second
    return (n0,n1,n2,n3)


# In[20]:


def CostReverse(R, city, n0, n1, n2, n3):
    # cost for reverse move
    de = Distance(R[city[n2]],R[city[n1]])+Distance(R[city[n0]],R[city[n3]])
    de-= Distance(R[city[n2]],R[city[n0]])+Distance(R[city[n1]],R[city[n3]])
    return de

def Reverse(R, city, n0, n1, n2, n3):
    newcity = copy(city)
    for j in range(n1-n0+1):
        newcity[n0+j] = city[n1-j]
    return newcity


# In[21]:


@jit(nopython=True)
def FindTSegment(R):
    (n0,n1,n2,n3) = FindASegment(R)
    nct = len(R)
    nn = nct - (n1-n0+1)  # number for the rest of the cities
    n4 = (n1+1 + int(rand()*(nn-1)) ) % nct # city on the rest of the path
    n5 = (n4+1) % nct
    return (n0,n1,n2,n3,n4,n5)


# In[22]:


def CostTranspose(R, city, n0,n1,n2,n3,n4,n5):
    de = -Distance(R[city[n1]], R[city[n3]])
    de-= Distance(R[city[n0]], R[city[n2]])
    de-= Distance(R[city[n4]], R[city[n5]])
    de+= Distance(R[city[n0]], R[city[n4]])
    de+= Distance(R[city[n1]], R[city[n5]])
    de+= Distance(R[city[n2]], R[city[n3]])
    return de


def Transpose(R, city, n0,n1,n2,n3,n4,n5):
    nct = len(R)
    newcity = []
    # Segment in the range n0,...n1
    for j in range(n1-n0+1):
        newcity.append( city[ (j+n0)%nct ] )
    # is followed by segment n5...n2
    for j in range( (n2-n5)%nct + 1):
        newcity.append( city[ (j+n5)%nct ] )
    # is followed by segement n3..n4
    for j in range( (n4-n3)%nct + 1):
        newcity.append( city[ (j+n3)%nct ] )
    return newcity


# In[25]:


nn = FindTSegment(R)
de = CostTranspose(R, city, *nn) # *nn = nn[0],nn[1],nn[2],nn[3],nn[4],nn[5], i.e., unpacking the list
print(nn)
print(de)
r1 = Transpose(R, city, *nn)
print(r1)


# In[26]:


def TravelingSalesman(city, R, maxSteps, maxAccepted, Tstart, fCool, maxTsteps, Preverse=0.5):
    T = Tstart
    dist = TotalDistance(city,R)
    for t in range(maxTsteps):
        accepted = 0
        for i in range(maxSteps):
            if Preverse > rand():
                # Try reverse
                nn = FindASegment(R)
                de = CostReverse(R, city, *nn)
                if de < 0 or exp(-de/T) > rand():
                    accepted += 1
                    dist += de
                    city = Reverse(R, city, *nn)
            else: 
                # here we transpose
                nn = FindTSegment(R)
                de = CostTranspose(R, city, *nn)
                if de < 0 or exp(-de/T) > rand():
                    accepted += 1
                    dist += de
                    city = Transpose(R, city, *nn)
            if accepted > maxAccepted: 
                break    
        T *= fCool
        Plot(city, R, dist)
        print("T=%10.5f , distance=%10.5f acc.steps=%d" % (T, dist,accepted))
        if accepted == 0:
            break
    Plot(city, R, dist)
    return city 


# In[28]:


from numpy import random

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

random.seed(1)

R = random.random((ncity,2))
city = range(ncity)

ncity = TravelingSalesman(city, R, maxSteps, maxAccepted, Tstart, fCool, maxTsteps)


# In[29]:


Plot(ncity,R, TotalDistance(ncity,R))


# ## Homework
# 
# 1) Implement Traveling Salesman problem.
# 
# 2) Test traveling salesman problem and plot distance $d$ versus $N$, where $N$ is number of cities. What is the optimal curve distance?

# In[ ]: