Simulated annealing in R

I am new in R and I have to implement simulated annealing for schaffer function and I did it. However I am not sure about the correctness of the code.

simulated_annealing <- function(schaffer, step, s0, niter) {
s_n <- s0
iter_count <- 0
for (k in 1:niter) {     
    T <- (1 - step)^k
    s_new <- rnorm(2, s_n, 1)
    dif <- func(s_new) - func(s_n)
    if (dif < 0) {
        s_n = s_new 
    } else {
        random <- runif(1,0,1)
        if (random < exp(-dif/T)) {
            s_n <- s_new
        }
    }
    iter_count <- iter_count + 1

    print(sprintf("The number of iteration was: %f",iter_count))
    print(sprintf("Function at s_new: %f", func(s_new)))
    print(sprintf("Function at s_n: %f", func(s_n)))
    print(sprintf("Current point: %f", s_n))
    print(sprintf("Temperature: %f", T))
}
return(s_n)
}

schaffer <- function(xx)
{x1 <- xx[1]
x2 <- xx[2]
fact1 <- (sin(x1^2-x2^2))^2 - 0.5
fact2 <- (1 + 0.001*(x1^2+x2^2))^2
y <- 0.5 + fact1/fact2
return(y)
}

It is how the simulated annealing works:

simulated_annealing(schaffer,0.1,c(0,2),10)

Output:

```
[1] "The number of iteration was: 1.000000"
[1] "Function at s_new: 0.899459"
[1] "Function at s_n: 0.572171"
[1] "Current point: 0.000000" "Current point: 2.000000"
[1] "Temperature: 0.900000"
[1] "The number of iteration was: 2.000000"
[1] "Function at s_new: 0.047084"
[1] "Function at s_n: 0.047084"
[1] "Current point: 1.022082" "Current point: 2.095950"
[1] "Temperature: 0.810000"
[1] "The number of iteration was: 3.000000"
[1] "Function at s_new: 0.153559"
[1] "Function at s_n: 0.047084"
[1] "Current point: 1.022082" "Current point: 2.095950"
[1] "Temperature: 0.729000"
[1] "The number of iteration was: 4.000000"
[1] "Function at s_new: 0.359053"
[1] "Function at s_n: 0.047084"
[1] "Current point: 1.022082" "Current point: 2.095950"
[1] "Temperature: 0.656100"
[1] "The number of iteration was: 5.000000"
[1] "Function at s_new: 0.828496"
[1] "Function at s_n: 0.047084"
[1] "Current point

Here is how I would rewrite your code, comments to follow.

simulated_annealing <- function(func, s0, niter = 10, step = 0.1) {

   # Initialize
   ## s stands for state
   ## f stands for function value
   ## b stands for best
   ## c stands for current
   ## n stands for neighbor
   s_b <- s_c <- s_n <- s0
   f_b <- f_c <- f_n <- func(s_n)
   message("It\tBest\tCurrent\tNeigh\tTemp")
   message(sprintf("%i\t%.4f\t%.4f\t%.4f\t%.4f", 0L, f_b, f_c, f_n, 1))

   for (k in 1:niter) {     
      Temp <- (1 - step)^k
      # consider a random neighbor
      s_n <- rnorm(2, s_c, 1)
      f_n <- func(s_n)
      # update current state
      if (f_n < f_c || runif(1, 0, 1) < exp(-(f_n - f_c) / Temp)) {
         s_c <- s_n
         f_c <- f_n
      }
      # update best state
      if (f_n < f_b) {
         s_b <- s_n
         f_b <- f_n         
      }
      message(sprintf("%i\t%.4f\t%.4f\t%.4f\t%.4f", k, f_b, f_c, f_n, Temp))
   }
   return(list(iterations = niter, best_value = f_b, best_state = s_b))
}

sol <- simulated_annealing(schaffer, s0 = c(0, 2))
#  It   Best    Current Neigh   Temp
#  0    0.5722  0.5722  0.5722  1.0000
#  1    0.1054  0.1054  0.1054  0.9000
#  2    0.1054  0.1922  0.1922  0.8100
#  3    0.1054  0.1922  0.3579  0.7290
#  4    0.1054  0.1922  0.9059  0.6561
#  5    0.1054  0.6428  0.6428  0.5905
#  6    0.1054  0.8893  0.8893  0.5314
#  7    0.1054  0.1597  0.1597  0.4783
#  8    0.1054  0.6489  0.6489  0.4305
#  9    0.1054  0.7535  0.7535  0.3874
#  10   0.1054  0.2236  0.2236  0.3487

sol
# $iterations
# [1] 10
#
# $best_value
# [1] 0.1053682
# 
# $best_state
# [1] -0.3203791  1.7080835

Now the comments, in no particular order:

• when defining your function, you meant to name the first argument func, not schaffer. That was a bug.

• niter and steps are the kind of arguments that could use an optional value so I made it so. And once an argument has an optional value, it is recommended to move it to the end of the argument list, hence function(func, s0, niter = 50, step = 0.1).

• I renamed T to Temp for 1) it is more descriptive and more likely to suggest that this is the annealing's temperature variable, but mostly 2) R provides T as a shortcut for TRUE, although you can overwrite it like you just did here. This is a horrible feature of the R language and you should avoid using T both for TRUE or as your own variable.

• message is probably a more appropriate tool for logging. I also transformed the code so it prints one row per iteration. I hope you'll agree the log is much easier to look at this way.

• Keeping track of the best state is an improvement over the "vanilla" version simulated annealing process which only reports the current state at the last iteration. If you want it that way, then you need to use three states: best, current, neighbor. Which I implemented here.

• In addition to keeping track of the current state and best state, it is recommended to keep track of their function values, so you do not have to make unnecessary function calls to re-compute them at each iteration.

• Like the matlab function you referred to, it is preferable to return a list so you can provide more outputs than just the best state.

• Always, it is important to comment your code to make it easier to understand and maintain.

I hope it helps.

StackExchange Code Review Q#84688, answer score: 16