function tsp_descent ( filename )

%*****************************************************************************80
%
%% MAIN is the main program for TSP_DESCENT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 October 2018
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILENAME, the name of the file containing the distance
%    matrix information.  This is an optional input.  If it is not provided,
%    the program will prompt for it.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TSP_DESCENT:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Solve small traveling salesman problems using\n' );
  fprintf ( 1, '  a descent method, generating reversals and\n' );
  fprintf ( 1, '  transpositions of a random starting tour.\n' );
%
%  Get the distance table.
%
  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The program needs to load a file containing\n' );
    fprintf ( 1, '  a city-to-city distance matrix D.\n' );  
    filename = input ( '  Enter the file name, in quotes:  ' );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Distance matrix filename is "%s"\n', filename );
  end

  distance = load ( filename );
%
%  Approve the distance table.
%
  [ m, n ] = size ( distance );

  if ( m ~= n )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TSP_DESCENT - Fatal error!\n' );
    fprintf ( 1, '  The distance matrix D must be square.\n' );
    fprintf ( 1, '  Your matrix has M = %d, N = %d\n', m, n );
    return
  end

  if ( n < 4 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TSP_DESCENT - Fatal error!\n' );
    fprintf ( 1, '  This problem is too small!\n' );
    fprintf ( 1, '  The number of cities N must be at least 4.\n' );
    fprintf ( 1, '  Your matrix has N = %d\n', n );
    return
  end

  v = diag ( distance );
  test = norm ( v );

  if ( 0.0 < test )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TSP_DESCENT - Fatal error!\n' );
    fprintf ( 1, '  The distance matrix D must have zero diagaonal.\n' );
    fprintf ( 1, '  Your matrix has ||diag(D)|| = %g\n', test );
    return
  end

  test = norm ( distance - distance' );

  if ( 0.0 < test )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TSP_DESCENT - Fatal error!\n' );
    fprintf ( 1, '  The distance matrix D must be symmetric.\n' );
    fprintf ( 1, '  Your matrix has ||D-D''|| = %g\n', test );
    return
  end
%
%  Print the distance matrix.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The city-to-city distance matrix D:\n' )
  fprintf ( 1, '\n' );

  disp ( distance )
%
%  Choose a random starting route.
%
  p = randperm ( n );
%
%  Save the cost.
%
  cost = path_cost ( n, distance, p );
%
%  Consider VARIATION_NUM variations.
%
  variation_num = 1000;
  transpose_num = 0;
  reversal_num = 0;

  for variation = 1 : variation_num
%
%  Variation 1: transpose.
%
%  Pick two distinct non-neighboring city indices.
%
    c = randperm ( n, 2 );
    c = sort ( c );
    i1 = c(1);
    i2 = c(2);
    if ( i2 == i1 + 1 )
      i2 = i2 + 1;
      if ( n < i2 )
        i2 = i1;
        i1 = 1;
      end
    end

    transpose_num = transpose_num + 1;
    p2 = [ p(1:i1), p(i2), p(i2-1:-1:i1+1), p(i2+1:n) ];
    cost2 = path_cost ( n, distance, p2 );

    if ( cost2 < cost )
      p = p2;
      cost = cost2;
      fprintf ( '  Transpose %d: cost = %g\n', transpose_num, cost );
    end 
%
%  Variation 2, reverse order.
%
%  Pick two distinct city indices.
%
    c = randperm ( n, 2 );
    c = sort ( c );
    i1 = c(1);
    i2 = c(2);

    reversal_num = reversal_num + 1;
    p2 = [ p(1:i1-1), p(i2:-1:i1), p(i2+1:n) ];
    cost2 = path_cost ( n, distance, p2 );

    if ( cost2 < cost )
      p = p2;
      cost = cost2;
      fprintf ( '  Reversal %d: cost = %g\n', reversal_num, cost );
    end 

  end
%
%  Report.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of variations tried was %d\n', 2 * variation_num );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The best itinerary found:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Step  From  To        Distance\n' );
  fprintf ( 1, '\n' );
  for i1 = 1 : n
    i2 = mod ( i1, n ) + 1;
    fprintf ( 1, '  %4d    %2d  %2d  %14.6g\n', ...
      i2, p(i1), p(i2), distance ( p(i1), p(i2) ) );
  end
  fprintf ( 1, '  ----    --  --  --------------\n' );
  fprintf ( 1, '  cost:           %14.6g\n', cost );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TSP_DESCENT\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function cost = path_cost ( n, distance, p )

%*****************************************************************************80
%
%% PATH_COST evaluates the cost of a round trip.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 October 2018
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of cities.
%
%    Input, real DISTANCE(N,N), the city to city distance table.
%
%    Input, integer P(N), a permutation of 1:N, the route.
%
%    Output, real COST, the cost of the route.
%
  cost = 0.0;
  i1 = n;
  for i2 = 1 : n
    cost = cost + distance ( p(i1), p(i2) );
    i1 = i2;
  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end
