function [ len, p ] = tsp_moler ( D )

%*****************************************************************************80
%
%% tsp_moler() seeks an optimal path for the traveling salesperson problem.
%
%  Discussion:
%
%    This is a functional form of an old MATLAB demo called "travel".
%
%    It is a pretty good, but certainly not the best, 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:
%
%    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
%
%  Usage:
%
%   [ len, p ] = traveler ( D )
%
%  Input:
%
%   real D(n,n): distances between cities.
%
%  Output:
%
%   real 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.
%
  n = size ( D, 1 );
%
%  Choose an initial path at random.
%
  p = randperm ( n );
  len = path_length ( p, 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);
    pnew = p(q);
%
%  If better, take it.
%
    lennew = path_length ( pnew, D );
    if ( lennew < len )
      p = pnew;
      len = lennew;
      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)) ];
    pnew = p(q);
%
%  If better, take it.
%
    lennew = path_length ( pnew, D );
    if ( lennew < len )
      p = pnew;
      len = lennew;
      iterp = iter;
    end

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

  return
end

