tsp_anneal <- function ( x, temp = 1.0e2, rate = 1.0e-4 )

#*****************************************************************************80
#
## tsp_anneal() solves the traveling salesperson problem (TSP) using simulated annealing.
#
#  Licensing:
#
#    Copyright 2016 James P. Howard, II
#
#    The computer code and data files on this web page are distributed under
#    https://opensource.org/licenses/BSD-2-Clause, the BSD-2-Clause license.
#
#  Modified:
#
#    04 June 2026
#
#  Author:
#
#    Original R code by James Howard;
#    This version by John Burkardt.
#
#  Reference:
#
#    James Howard,
#    Computational Methods for Numerical Analysis with R,
#    CRC Press, 2017,
#    ISBN13: 978-1-4987-2363-3.
#
#  Input:
#
#    real x(n,2): the coordinates of the cities.
#
#    real temp: the annealing temperature.
#
#    real rate: the annealing rate.
#
#  Output:
#
#    integer order(n): the optimal ordering.
#
#    real distance: the length of the optimal round trip.
#
{
  source ( "/home/john/public_html/r_src/vecnorm/vecnorm.R" )  

  step <- 1.0 - rate
#
#  N is the number of cities.
#
  n <- nrow ( x )
#
#  Create n-vectors xbest, xnext, xcurr.
#
# xbest <- c(1:n)
  xbest <- sample ( 1:n )
  xnext <- xbest
  xcurr <- xbest
#
#  Compute the lengths ybest, ynext, ycurr.
#
  ynext <- 0.0
  from = n
  for ( to in 1 : n )
  {
    ynext <- ynext + vecnorm ( x[xnext[from],] - x[xnext[to],] )
    from = to
  }

  ybest <- ynext
  ycurr <- ynext
#
#  As the temperature decreases, repeatedly adjust the current path.
#
  while ( 1.0 < temp )
  {
#
#  Decrease the temperature.
#
    temp <- temp * step
#
#  Pick a path component i to swap.
#  ceiling() implicitly rules out the case i=1.
#
    i <- ceiling ( runif ( 1, 1, n ) )
#
#  Create the next path candidate, by swapping positions i and i-1.
#
    xnext <- xcurr
#
#  Swap positions of cities i and i-1.
#
    temporary <- xnext[i]
    xnext[i] <- xnext[i-1]
    xnext[i-1] <- temporary
#
#  Compute the path length.
#
    ynext <- 0.0
    from = n
    for ( to in 1 : n )
    {
      ynext <- ynext + vecnorm ( x[xnext[from],] - x[xnext[to],] )
      from = to
    }
#
#  Acceptance depends in part on the current "temperature".
#
    accept <- exp ( - ( ynext - ycurr ) / temp )

    if ( ynext < ycurr || runif ( 1 ) < accept )
    {
      xcurr <- xnext
      ycurr <- ynext
    }

    if ( ynext < ybest )
    {
      xbest <- xcurr
      ybest <- ycurr
    }
  }

  return ( list ( order = xbest, distance = ybest ) )
}