nm <- function(x0, h=1, eps=0.000001, a=1, b=0.5, g=2, d=0.5) {
# A simple implementation of the Nelder-Mead method 
# for optimization of the Rosenbrock function (a demonstration)
#
# x0 ... the n-dimensional initial solution
# h ... initial size of the simplex
# eps ... the terminating measure of the size of the simplex
# a,b,g,d ... the "alpha, beta, gamma, delta" parameters of the method

# The Rosenbrock function to be minimized (the algorithm should work for any f)
f <- function(x) {
  (1 - x[1])^2 + 100 * (x[2] - x[1]^2)^2
}

# A measure of the size of the simplex defined by the matrix X
sz <- function(X, n) {
  (abs(det(X[2:(n + 1),] - matrix(rep(X[1,], n), byrow=TRUE, nrow=n))))^(1 / n)
}

# Plot the image of the Rosenbrock function
yy <- xx <- seq(-10, 10, by=0.1)
zz <- yy%o%xx
for(i in 1:201) for(j in 1:201) zz[i,j] <- f(c(xx[i], yy[j]))
image(xx, yy, zz, xlab="x", ylab="y", col=topo.colors(15), breaks=exp(0:15) - 1.5)
points(1, 1, pch=19, col="red")

# Construct the initial simplex
n <- length(x0)
X <- matrix(rep(x0, n + 1), byrow=TRUE, nrow=n + 1)
for(i in 2:(n + 1)) X[i, i - 1] <- x0[i - 1] + h

# The main cycle
val <- rep(0, n + 1)
while(sz(X, n) > eps) {

  # Plot the actual simplex
  lines(X[1:2, 1], X[1:2, 2], lwd=2)
  lines(X[2:3, 1], X[2:3, 2], lwd=2)
  lines(X[c(1, 3), 1], X[c(1, 3), 2], lwd=2)
  
  # Delay the execution to be able to see the animation
  Sys.sleep(0.3)

  for(i in 1:(n + 1)) val[i] <- f(X[i,])   # This can be implemented more efficiently (most of these values have already been computed).
 
  ord <- order(val)  # This can also be implemented somewhat more efficiently (we do not need the full ordering).
  val <- val[ord]
  X <- X[ord,]       
  xl <- X[1,]
  xs <- X[n,]
  xh <- X[n + 1,]
  fl <- val[1]
  fs <- val[n]
  fh <- val[n + 1]
 
  cen <- apply(X[1:n,], 2, mean)           # Compute the centroid.
  xr <- cen + a * (cen - xh); fr <- f(xr)    # Compute the reflection.
 
  if((fl <= fr) & (fr < fs)) {
    X[n + 1,] <- xr
    next
  }
  if(fr < fl) {
    xe <- cen + g * (xr - cen); fe <- f(xe) # Compute the expanded point.
    if(fe < fr) {
      X[n + 1,] <- xe
    } else {
      X[n + 1,] <- xr
    }
    next
  }
  if(fr >= fs) {
    if(fr < fh) {
      xc <- cen + b * (xr - cen); fc <- f(xc) # Compute the outer contraction point.
      if(fc <= fr) {
        X[n + 1,] <- xc
        next
      }
    } else {
      xc <- cen + b * (xh - cen); fc <- f(xc) # Compute the inner contraction point.
      if(fc < fh) {
        X[n + 1,] <- xc
        next
      }
    }
  }
  
  for(i in 2:(n+1)) X[i,]<-xl+d*(X[i,]-xl) # Shrink.
}

return(X[1,])
}

# Demo:
nm(c(-7, 5), h=1, eps=0.000001, a=1, b=0.5, g=2, d=0.5)
nm(c(-7, -7), h=1, eps=0.000001, a=1, b=0.5, g=2, d=0.5)
nm(c(7, -5), h=1, eps=0.000001, a=1, b=0.5, g=2, d=0.5)

# Odskusame funkciu optim na Rosenbrockovi
f <- function(x){(1 - x[1])^2 + 100 * (x[2] - x[1]^2)^2}

optim(c(7, -5), f, method="Nelder-Mead")
res <- optim(c(7, -5), f, method="Nelder-Mead")
optim(res$par, f, method="Nelder-Mead")

optim(c(7, -5), f, method="CG")
optim(c(7, -5), f, method="CG", control=list(maxit=10000))

optim(c(7, -5), f, method="BFGS")