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

#*****************************************************************************80
#
## tspsa solves the traveling salesperson problem 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:
#
#    21 March 2020
#
#  Author:
#
#    Original R code by James Howard;
#    Modifications by John Burkardt.
#
#  Reference:
#
#    James Howard,
#    Computational Methods for Numerical Analysis with R,
#    CRC Press, 2017,
#    ISBN13: 978-1-4987-2363-3.
#
{
  source ( "/home/burkardt/public_html/r_src/vecnorm/vecnorm.R" )  

  step <- 1.0 - rate
  n <- nrow ( x )

  xbest <- c(1:n)
  xnext <- c(1:n)
  xcurr <- c(1:n)

  ynext <- 0.0
  for ( i in 2 : n )
  {
    a <- xnext[i-1]
    b <- xnext[i]
    ynext <- ynext + vecnorm ( x[a,] - x[b,] )
  }
  a <- xnext[n]
  b <- xnext[1]
  ynext <- ynext + vecnorm ( x[a,] - x[b,] )

  ybest <- ynext
  ycurr <- ynext

  while ( 1.0 < temp )
  {

    temp <- temp * step
    i <- ceiling ( runif ( 1, 1, n ) )
    xnext <- xcurr
    temporary <- xnext[i]
    xnext[i] <- xnext[i-1]
    xnext[i-1] <- temporary

    ynext <- 0.0
    for ( i in 2 : n )
    {
      a <- xnext[i-1]
      b <- xnext[i]
      ynext <- ynext + vecnorm ( x[a,] - x[b,] )
    }
    a <- xnext[n]
    b <- xnext[1]
    ynext <- ynext + vecnorm ( x[a,] - x[b,] )

    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 ) )
}