#! /usr/bin/env python3
#
def tsp_test ( ):

#*****************************************************************************80
#
## tsp_test() tests tsp().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    28 March 2026
#
#  Author:
#
#    John Burkardt
#
  from numpy.random import default_rng
  import matplotlib.pyplot as plt
  import numpy as np
  import platform

  print ( '' )
  print ( 'tsp_test():' )
  print ( '  numpy version:  ' + np.version.version )
  print ( '  python version: ' + platform.python_version ( ) )

  rng = default_rng ( )
#
#  Get the data.
#
  state, lat, lon = state_data ( )

  state_num = len ( state )
#
#  Verify that the data has been read correctly.
#
  print ( '' )
  print ( '  First 5 state, latitudes, longitudes values:' )
  for i in range ( 0, 5 ):
    print ( '  %s  %12.4f  %12.4f' % ( state[i], lat[i], lon[i] ) )
#
#  Compute the distance matrix.
#
  D = distances ( lat, lon )

  print ( '' )
  print ( '  5x5 initial block of distance matrix' )
  print ( D[0:5,0:5] )
#
#  Call the TSP solver.
#
  p_len, p = traveler ( D, rng ) 
#
#  Report the results.
#
  print ( '' )
  print ( '  Length of path = ', p_len )
  print ( '' )
  print ( '  Itinerary:' )
  print ( '  ', end = '' )
  for i in range ( 0, state_num ):
    print ( state[p[i]] + '->', end = '' )
    if ( ( ( i + 1 ) % 15 ) == 0 ):
      print ( '' )
      print ( '  ', end = '')
  print ( state[p[0]] )
#
#  Give the mileages.
#
  print ( '' )
  print ( '  Mileage per leg:' )
  print ( '')
  for leg in range ( 0, state_num ):
    i = p[leg]
    j = p[leg+1]
    print ( '  %s -> %s  %g' % ( state[i], state[j], D[i,j] ) )
#
#  Plot (longitude,latitude) as (X,Y).
#
  plt.clf ( )
  plt.plot ( lon, lat, 'bo', markersize = 10 )
  for leg in range ( 0, state_num ):
    i = p[leg]
    j = p[leg+1]
    x = [ lon[i], lon[j] ]
    y = [ lat[i], lat[j] ]
    plt.plot ( x, y, 'r-' )
  for i in range ( 0, state_num ):
    x = lon[i]
    y = lat[i]
    plt.text ( x, y, state[i] )
  plt.grid ( True )
  plt.xlabel ( '<-- Longitude -->' )
  plt.ylabel ( '<-- Latitude -->' )
  s = ( 'TSP approximate solution with length %g miles' % p_len )
  plt.title ( s )
  filename = 'tsp.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "'+ filename + '"' )
  plt.close ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'tsp_test():' )
  print ( '  Normal end of execution.' )

  return

def state_data ( ):

#*****************************************************************************80
#
## state_data() the names, latitude, and longitude of 48 states.
#
#  Discussion:
#
#    Alaska and Hawaii are omitted.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    28 March 2026
#
#  Author:
#
#    John Burkardt.
#
#  Output:
#
#    string state[48]: the 2-letter abbreviated state name.
#
#    real lat[48], lon[48]: the latitude and longitude of the state
#    capitals, in degrees.
#
  import numpy as np

  state = np.array ( [ \
    'AL', 'AZ', 'AR', 'CA', 'CO', 'CT', 'DE', 'FL', 'GA', 'ID', \
    'IL', 'IN', 'IA', 'KS', 'KY', 'LA', 'ME', 'MD', 'MA', 'MI', \
    'MN', 'MS', 'MO', 'MT', 'NE', 'NV', 'NH', 'NJ', 'NM', 'NY', \
    'NC', 'ND', 'OH', 'OK', 'OR', 'PA', 'RI', 'SC', 'SD', 'TN', \
    'TX', 'UT', 'VT', 'VA', 'WA', 'WV', 'WI', 'WY' ] )

  lat = np.array ( [
    32.361538, 33.448457, 34.736009, 38.555605, 39.7391667, \
    41.767,    39.161921, 30.4518,   33.76,     43.613739,  \
    39.783250, 39.790942, 41.590939, 39.04,     38.197274,  \
    30.45809,  44.323535, 38.972945, 42.2352,   42.7335,    \
    44.95,     32.320,    38.572954, 46.595805, 40.809868,  \
    39.160949, 43.220093, 40.221741, 35.667231, 42.659829,  \
    35.771,    46.8055,   39.962245, 35.482309, 44.931109,  \
    40.269789, 41.82355,  34.000,    44.367966, 36.165,     \
    30.266667, 40.7547,   44.26639,  37.54,     47.042418,  \
    38.349497, 43.074722, 41.145548 ] )

  lon = np.array ( [
    -86.279118, -112.073844,  -92.331122, -121.468926, -104.984167, \
    -72.677,     -75.526755,  -84.27277,   -84.39,     -116.237651, \
    -89.650373,  -86.147685,  -93.620866,  -95.69,      -84.86311,  \
    -91.140229,  -69.765261,  -76.501157,  -71.0275,    -84.5467,   \
    -93.094,     -90.207,     -92.189283, -112.027031,  -96.675345, \
   -119.753877,  -71.549127,  -74.756138, -105.964575,  -73.781339, \
    -78.638,    -100.779004,  -83.000647,  -97.534994, -123.029159, \
    -76.875613,  -71.422132,  -81.035,    -100.336378,  -86.784,    \
    -97.75,     -111.892622,  -72.57194,   -77.46,     -122.893077, \
    -81.633294,  -89.384444, -104.802042 ] )

  return state, lat, lon

def distances ( lat, lon ):

#*****************************************************************************80
#
## distances() computes a distance table from latitude and longitudes.
#
#  Discussion:
#
#    The distances are computed on the surface of a sphere representing
#    the earth.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    14 March 2023
#
#  Author:
#
#    John Burkardt.
#
#  Reference:
#
#    Great circle distance matrix between pairs of cities.
#    https://en.wikipedia.org/wiki/Great-circle_distance.
#
#    Cleve Moler,
#    USA Traveling Salesman Tour,
#    Posted 17 September 2018,
#    blogs.mathworks.com/cleve/2018/09/17/usa-traveling-salesman-tour/?s_tid=srchtitle_traveling#20salesman_1
#
#  Input:
#
#    real lat(*), lon(*): the longitude and latitude of the cities,
#    in degrees.
#
#  Output:
#
#    real D(*,*):  D(i,j) is the distance in miles between cities i and j.
#
#  Local:
#
#    real R: the radius of the earth in miles.
#
  import numpy as np

  R = 3959.0

  n = len ( lon )
  D = np.zeros ( [ n, n ] )
  lat = np.deg2rad ( lat )
  lon = np.deg2rad ( lon )

  for i in range ( 0, n ):
    for j in range ( 0, i ):
      D[i,j] = R * np.arccos ( np.sin ( lat[i] ) * np.sin ( lat[j] ) + \
        np.cos ( lat[i] ) * np.cos ( lat[j] ) * np.cos ( lon[i] - lon[j] ) )
      D[j,i] = D[i,j]

  return D

def path_length ( p, D ):

#*****************************************************************************80
#
## path_length() computes the current path length.
#
#  Discussion:
#
#    This function calculates the total length of the current path p
#    in the traveling salesman problem, using a distance matrix D.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    14 March 2023
#
#  Author:
#
#    John Burkardt.
#
#  Reference:
#
#    Cleve Moler,
#    USA Traveling Salesman Tour,
#    Posted 17 September 2018,
#    logs.mathworks.com/cleve/2018/09/17/usa-traveling-salesman-tour/?s_tid=srchtitle_traveling#20salesman_1
#
#  Usage:
#
#    p_len = path_length ( p, D )
#
#  Input:
#
#    integer p(n): the sequence of cities in the current path.
#
#    real D(n,n): the city-to-city distance matrix.
#
#  Output:
#
#    real p_len: the current length of the path.
#
  n = D.shape[0]

  c2 = p[n-1]
  p_len = 0.0

  for i in range ( 0, n ):
    c1 = c2
    c2 = p[i]
    p_len = p_len + D[c1,c2]

  return p_len

def timestamp ( ):

#*****************************************************************************80
#
## timestamp() prints the date as a timestamp.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    21 August 2019
#
#  Author:
#
#    John Burkardt
#
  import time

  t = time.time ( )
  print ( time.ctime ( t ) )

  return

def traveler ( D, rng ):

#*****************************************************************************80
#
## traveler() seeks an optimal path for the traveling salesperson problem.
#
#  Discussion:
#
#    It is a neuristic algorithm for an approximate solution of the
#    Traveling Salesman Problem (TSP).  
#
#    The goal is to form a closed circuit of a number of cities that 
#    requires the shortest total distance.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    14 March 2023
#
#  Author:
#
#    John Burkardt.
#
#  Reference:
#
#    Cleve Moler,
#    USA Traveling Salesman Tour,
#    Posted 17 September 2018,
#    logs.mathworks.com/cleve/2018/09/17/usa-traveling-salesman-tour/?s_tid=srchtitle_traveling#20salesman_1
#
#  Usage:
#
#   p_len, p = traveler ( D )
#
#  Input:
#
#    real D(n,n): distances between cities.
#
#    rng(): the current random number generator.
#
#  Output:
#
#   real p_len: the length of the itinerary in miles.
#
#   integer p(n+1): the cities listed in the order in which they will be
#   visited, including a return to the starting city.
#
  import numpy as np

  n = D.shape[0]
#
#  Choose an initial path at random.
#
  p = rng.permutation ( n )
  p_len = path_length ( p, D )

  iter_max = 10000

  for iter in range ( 0, iter_max ):
#
#  Try to reverse a portion of the path.
#
    pt1 = rng.integers ( 0, n )
    pt2 = rng.integers ( 0, n )

    lo = min ( pt1, pt2 )
    hi = max ( pt1, pt2 )

    q = np.arange ( n )
    q[lo:hi+1] = np.flip ( q[lo:hi+1] )
    pnew = p[q]
#
#  If better, take it.
#
    lennew = path_length ( pnew, D )
    if ( lennew < p_len ):
      p = pnew.copy()
      p_len = lennew
      iterp = iter
#
#  Try a single point insertion.
#
    pt1 = rng.integers ( 0, n )
    pt2 = rng.integers ( 0, n - 1 )

    q = np.arange ( n )
    q = np.delete ( q, pt1 )
    q = np.insert ( q, pt2+1, pt1 )
    pnew = p[q]
#
#  If better, take it.
#
    lennew = path_length ( pnew, D )
    if ( lennew < p_len ):
      p = pnew.copy()
      p_len = lennew
      iterp = iter
#
#  Make this a round trip by padding a copy of the first city to the end.
#
  p = np.append ( p, p[0] )

  return p_len, p

if ( __name__ == "__main__" ):
  timestamp ( )
  tsp_test ( )
  timestamp ( )