sa <- function ( f, x, temp = 1.0e4, rate = 1.0e-4 )

#*****************************************************************************80
#
## sa minimizes a multivariate function 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:
#
#    18 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.
#
{
  step <- 1.0 - rate
  n <- length ( x )

  xbest <- x
  ybest <- f(xbest)
  xnext <- x
  ynext <- f(xnext)
  xcurr <- x
  ycurr <- f(xcurr)

  while ( 1.0 < temp )
  {

    temp <- temp * step
    i <- ceiling ( runif ( 1, 0, n ) )
    xnext[i] <- rnorm ( 1, xcurr[i], temp )
    ynext <- f(xnext)
    accept <- exp ( - ( ynext - ycurr ) / temp )

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

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

  }
  x <- xbest

  return ( x )
}