function [ cost, order ] = tsp_moler ( D )

%*****************************************************************************80
%
%% tsp_moler() seeks an optimal path for the traveling salesperson problem.
%
%  Discussion:
%
%    The goal is to form a closed circuit of a number of cities that 
%    requires the shortest total distance.
%
%  Licensing:
%
%    Copyright 1984-2018 The MathWorks, Inc.
%
%  Modified:
%
%    15 June 2026
%
%  Author:
%
%    Original MATLAB version by Cleve Moler.
%    This version by John Burkardt
%
%  Reference:
%
%    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 D(n,n): distances between cities.
%
%  Output:
%
%   real cost: the length of the itinerary in miles.
%
%   integer order(n+1): the cities listed in the order in which they will be
%   visited, including a return to the starting city.
%
  n = size ( D, 1 );
%
%  Choose an initial path at random.
%
  order = randperm ( n );
  cost = path_length ( order, D );

  for iter = 1 : 10000
%
%  Try to reverse a portion of the path.
%
    pt1 = floor ( n * rand() ) + 1;
    pt2 = floor ( n * rand() ) + 1;

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

    q = 1 : n;
    q(lo:hi) = q(hi:-1:lo);
    order_new = order(q);
%
%  Compute the cost of the reversal.
%
    cost_new = path_length ( order_new, D );
%
%  If the reversal lowers the cost, make it the best so far.
%
    if ( cost_new < cost )
      order = order_new;
      cost = cost_new;
      iterp = iter;
    end
%
%  Try a single point insertion.
%
    pt1 = floor ( n * rand() ) + 1;
    pt2 = floor ( ( n - 1 ) * rand() ) + 1;

    q = 1 : n;
    q(pt1) = [];
    q = [ q(1:pt2), pt1, q((pt2+1):(n-1)) ];
    order_new = order(q);
%
%  Compute the cost of the insertion.
%
    cost_new = path_length ( order_new, D );
%
%  If the insertion lowers the cost, make it the best so far.
%
    if ( cost_new < cost )
      order = order_new;
      cost = cost_new;
      iterp = iter;
    end

  end
%
%  Close the permutation.
%
  order(end+1) = order(1);

  return
end
function cost = path_length ( order, 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.
%
%    This is a vectorized distance calculation.
%
%  Licensing:
%
%    Copyright 1984-2018 The MathWorks, Inc.
%
%  Modified:
%
%    13 March 2023
%
%  Author:
%
%    Original MATLAB version by Cleve Moler.
%    This version by John Burkardt
%
%  Reference:
%
%    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:
%
%    integer order(n): the sequence of cities in the current path.
%
%    real D(n,n): the city-to-city distance matrix.
%
%  Output:
%
%    real cost: the current length of the path.
%
%  Local:
%
%    integer q(n): the pair ( order(i), q(i) ) describes the i-th leg of the path.
%   
  n = length ( order );
  q = order( [ n, 1:(n-1) ] );
%
%  Use column indexing into D to create a vector of
%  lengths of each leg which can then be summed.
%
  cost = sum ( D(q + (order-1)*n));

  return
end
