etsp <- function(coor, x=NULL, sann=TRUE, method=2, tau=10000, N=100000, graph=10000)
{
   # Demo pre (znahodneny pazravy) Hill climbing a pre simulovane zihanie
   # Aplikovane na "Euklidovsky" problem obchodneho cestujuceho (Euclidean TSP)
   
   # coor ... n krat 2 matica udavajuca koordinaty n "miest" v rovine
   # x ... pociatocne riesenie; NULL znamena plne nahodne
   # sann... TRUE znamena simulovane zihanie, inak Hill climbing
   # method ... metoda generovania kandidatov (pripustne hodnoty: 1,2)
   # tau ... pociatocna teplota pre simulovane zihanie
   # N ... celkovy pocet iteracii
   # graph ... ako casto sa ma vykreslovat graf; ak NULL, tak nevykresluj nic
   #
   # Poznamky: Pokles teploty je linearny; moze sa volit rozne.
   
   
   Lng <- function(x, n, D) {
      # x je permutacia 1,...,n urcujuca Hamiltonovsku kruznicu v grafe
      # D je matica vzajomnych vzdialenosti
      # Vystupom je dlzka cesty urcenej x-om
      
      l <- D[x[1], x[n]]
      for (i in 1:(n - 1))
         l <- l + D[x[i], x[i + 1]]
      
      return(l)
   }
   
   # Vypocet matice vzajomnych vzdialenosti
   n <- nrow(coor); D <- matrix(0, ncol = n, nrow = n)
   for (i in 1:(n - 1)) {
      for (j in i:n) {
         D[i, j] <- D[j, i] <- sqrt(sum((coor[i, ] - coor[j, ])^2))
      }
   }
   
   # Vektor, v ktorom si pamatame dlzky ciest (len kvoli grafu).
   lengths <- rep(0, N) 
   
   # Vygenerovanie nahodneho pociatocneho riesenia x, inicializacia xmin, Lmin
   if (is.null(x)) x <- sample(1:n)
   Lx <- Lng(x, n, D)
   Lmin <- Lx
   xmin <- x
   
   if (!is.null(graph)) windows()
   
   # Hlavny cyklus
   for (i in 1:N) {
      lengths[i] <- Lx
      
      # Vygenerovanie kandidata
      ind <- sort(sample(1:n, 2))
      
      if (method == 1) {
         y <- x
         y[ind[1]] <- x[ind[2]]
         y[ind[2]] <- x[ind[1]]
      } else {
         reversed <- rev(x[ind[1]:ind[2]])
         if (ind[1] == 1) {
            left <- c()
         } else {
            left <- x[1:(ind[1] - 1)]
         }
         if (ind[2] == n) {
            right <- c()
         } else {
            right <- x[(ind[2] + 1):n]
         }
         y <- c(left, reversed, right)
      }
      
      Ly <- Lng(y, n, D)
      if (Ly < Lmin) {
         Lmin <- Ly
         xmin <- y
      }
      
      if (sann) {
         if (runif(1) < exp((Lx - Ly) / (tau / i))) {
            x <- y; Lx <- Ly
         }
      } else {
         if (Ly < Lx) {
            x <- y; Lx <- Ly
         }
      }
      
      # Vykreslenie momentalnej cesty
      if (!is.null(graph) && i %% graph == 0) {
         plot(coor[, 1], coor[, 2], type = "n", cex = 2, col = "red",
              main = c(i, round(Lmin, 3), round(Lx, 3)), xlab = "", ylab = "")
         grid()
         for (i in 1:(n - 1))
            lines(c(coor[x[i], 1], coor[x[i + 1], 1]), c(coor[x[i], 2],coor[x[i + 1], 2]))
         lines(c(coor[x[1], 1], coor[x[n], 1]), c(coor[x[1], 2], coor[x[n],2]))
         points(coor[, 1], coor[, 2], cex = 2, col = "red", pch = 19)
      }
   }
   # Koniec vypoctu
   
   # Vykreslenie najlepsej najdenej cesty a casoveho profilu ucelovej funkcie 
   if (!is.null(graph)) {
      plot(coor[, 1], coor[, 2], type = "n", cex = 2, col = "red",
           main = c(i, round(Lmin, 3), round(Lx, 3)), xlab = "", ylab = "")
      grid()
      for (i in 1:(n - 1))
         lines(c(coor[xmin[i], 1], coor[xmin[i + 1], 1]),
               c(coor[xmin[i], 2], coor[xmin[i + 1], 2]), col = "blue", lwd = 2)
      lines(c(coor[xmin[1], 1], coor[xmin[n], 1]),
            c(coor[xmin[1], 2], coor[xmin[n], 2]), col = "blue", lwd = 2)
      points(coor[, 1], coor[, 2], cex = 2, col = "red", pch = 19)
      windows(); plot(lengths, type = "l")
   }
   
   return(list(x = xmin, L = Lmin))
}

# Demo

cities <- cbind(runif(30), runif(30))
etsp(cities, x = NULL, sann = TRUE, method = 2, tau = 1000, N = 100000, graph = 10000)

t <- 2 * pi * seq(1/12, 1, by = 1/12)
cities <- matrix(0, nrow = 60, ncol = 2)
for (i in 0:11) {
   cities[(i * 5 + 1):(i * 5 + 5), 1] <- cos(t[i + 1]) * (1:5)
   cities[(i * 5 + 1):(i * 5 + 5), 2] <- sin(t[i + 1]) * (1:5)
}

# Pazravy znahodneny Hill climbing
etsp(cities, x = NULL, sann = FALSE, method = 1, tau = 30000, N = 300000, graph = 10000)

# Skusime multistart
L.best <- Inf
for (i in 1:30) {
   res <- etsp(cities, x = NULL, sann = FALSE, method = 1, tau = 30000, N = 10000, graph = NULL)
   if (res$L < L.best) {
      L.best <- res$L
      print(L.best)
   }
}

# Simulovane zihanie
etsp(cities, x = NULL, sann = TRUE, method = 1, tau = 30000, N = 300000, graph = 10000)
etsp(cities, x = NULL, sann = TRUE, method = 2, tau = 30000, N = 300000, graph = 10000)