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

  import matplotlib.pyplot as plt
  import numpy as np

  print ( '' )
  print ( 'fifteen_test():' )
  print ( '  TSP for a set of fifteen cities.' )

  position = np.loadtxt ( 'fifteen_position.txt' )
  city_num = position.shape[0]
  names = [ \
    'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O' ]
#
#  Plot the cities.
#  Probably DON'T want to see this.
#
  if ( False ):
    from graphviz import Graph
    dot = Graph ( format = 'png' )
    for n in names:
      dot.node ( n )
    for i in range ( 0, city_num ):
      for j in range ( i + 1, city_num ):
        dot.edge ( names[i], names[j] )
    dot.render ( 'fifteen.dot', view = True )
    filename = 'fifteen.dot.png'
    print ( '  Graphics saved as "%s"' % ( filename ) )
#
#  Compute distance matrix.
#
  distance = np.zeros ( [ city_num, city_num ] )
  for i in range ( 0, city_num ):
    for j in range ( 0, city_num ):
      distance[i,j] = np.round ( np.linalg.norm ( position[i,:] - position[j,:] ) )

  print ( '' )
  print ( '  Distance table:' )
  print ( distance )
  input ( 'Press RETURN to continue...' )
#
#  Plot the cities with distances.  Probably DON'T want to see this.
#
  if ( False ):
    from graphviz import Graph
    dot = Graph ( format = 'png' )
    for n in names:
      dot.node ( n )
    for i in range ( 0, city_num ):
      for j in range ( i + 1, city_num ):
        dot.edge ( names[i], names[j], str ( distance[i,j] ) )
    dot.render ( 'fifteen_distance.dot', view = True )
    filename = 'fifteen_distance.dot.png'
    print ( '  Graphics saved as "%s"' % ( filename ) )
#
#  Greedy estimate.
#
  print ( '' )
  print ( '  Estimate TSP solution using greedy algorithm.' )

  mileage_best = np.inf
  p_best = []

  for start in range ( 0, city_num ):

    p = path_greedy ( distance, start )
    mileage = path_length ( distance, p )

    if ( mileage < mileage_best ):
      p_best = p.copy()
      mileage_best = mileage
      print ( '  Greedy tour starting at city', start, ' costs ', mileage )

  print ( '' )
  print ( '  Greedy chose the following itinerary:' )
  print ( '' )
  print ( '  Step  From  To        Distance' )
  print ( '' )

  for i1 in range ( 0, city_num ):
    i2 = ( ( i1 + 1 ) % city_num )
    print ( '  %4d    %s  %s  %14.6g' \
      % ( i1, names[p_best[i1]], names[p_best[i2]], distance[p_best[i1], p_best[i2] ] ) )

  print ( '  ----   -- --  --------------' )
  print ( '  Total:        %14.6g' % ( mileage_best ) )
#
#  Investigate path length distribution by random sampling.
#
  input ( 'Press RETURN to continue...' )

  from numpy.random import default_rng

  rng = default_rng ( )

  sample_num = 1000000
  mileage = np.zeros ( sample_num )

  print ( '' )
  print ( '  Generate', sample_num, 'random permutations' )
  print ( '  and look at how the mileages are distributed.' )

  mileage_best = np.inf
  p_best = []

  for sample in range ( 0, sample_num ):
#
#  Choose a random starting route.
#
    p = rng.permutation ( np.arange ( city_num ) )
#
#  Compute the cost.
#
    mileage[sample] = path_length ( distance, p )

    if ( mileage[sample] < mileage_best ):
      p_best = p.copy ( )
      mileage_best = mileage[sample]
      print ( '  Sample %d: mileage = %g' % ( sample, mileage[sample] ) )
#
#  Report.
#
  print ( '' )
  print ( '  Number of samples tried was', sample_num )
  print ( '  The best itinerary found by random samples:' )
  print ( '' )
  print ( '  Step  From  To        Distance' )
  print ( '' )
  for i1 in range ( 0, city_num ):
    i2 = ( ( i1 + 1 ) % city_num )
    print ( '  %4d    %s  %s  %14.6g' \
      % ( i1, names[p_best[i1]], names[p_best[i2]], distance[p_best[i1], p_best[i2] ] ) )

  print ( '  ----   -- --  --------------' )
  print ( '  Total:        %14.6g' % ( mileage_best ) )

  plt.clf ( )
  plt.hist ( mileage, bins = 21 )
  plt.grid ( True )
  plt.title ( 'Distribution of mileages with random sampling' )
  plt.xlabel ( 'mileage' )
  plt.ylabel ( 'frequency' )
  filename = 'fifteen_sampling.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( )
  plt.close ( )

  return

def path_length ( distance, p ):

#*****************************************************************************80
#
## path_length() evaluates the cost of a round trip.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    06 March 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    real DISTANCE(N,N), the city to city distance table.
#
#    integer P(N), a permutation of 0:N-1, the route.
#
#  Output:
#
#    real mileage: the cost of the route.
#
  mileage = 0.0
  c0 = p[-1]
  for c1 in p:
    mileage = mileage + distance[c0,c1]
    c0 = c1

  return mileage

def path_greedy ( distance, start ):

#*****************************************************************************80
#
## path_greedy() finds a greedy route for a given start.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    06 March 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    real DISTANCE(N,N), the city to city distance table.
#
#    integer START, the starting city.
#
#  Output:
#
#    integer P(N), a greedy route that starts at START.
#
  import numpy as np

  n = distance.shape[0]

  p = np.zeros ( n, dtype = int )
  p[0] = start

  d = distance.copy ( )
  d[:,start] = np.inf

  for i in range ( 0, n ):
    d[i,i] = np.inf

  c1 = start
  for j in range ( 1, n ):
    c2 = np.argmin ( d[c1,:] )
    p[j] = c2
    d[:,c2] = np.inf
    c1 = c2

  return p

fifteen_test ( )