# Zaznam prikazov pre testovanie roznych optimalizacnych metod na vypocet 
# geometrickeho (L1) medianu
# Predmet "Stochasticke optimalizacne metody", Radoslav Harman, FMFI UK

sumdist <- function(x, X) 
{
   # x ... m-rozmerny vektor
   # X ... matica rozmeru nxm s polohami n datovych bodov v R^m
   # Vysledok je suma vzdialenosti bodu x ku vsetkym datovym bodom
   
   return(sum(sqrt(rowSums((X - rep(1,nrow(X)) %*% t(x))^2))))
}

# Otestujme si rychlost vypoctu pre maly aj vacsi dataset
X <- matrix(rnorm(20), ncol = 2); x <- rnorm(2)
system.time({for (i in 1:10000) res <- sumdist(x, X)})
X <- matrix(rnorm(2000), ncol = 2); x <- rnorm(2)
system.time({for (i in 1:10000) res <- sumdist(x, X)})

# Poznamky:
# Na jemnejsie casovanie je vhodna kniznica microbenchmark
# Na hlbsiu analyzu casovej narocnosti je dobre urobit "profiling"

plot.sumdist <- function(X) 
{
   # X ... matica rozmeru nxm s polohami n datovych bodov v R^m
   # Nakresli 201x201 filled contour pre funkciu sumdist
   
   min1 <- min(X[,1]); max1 <- max(X[,1])
   min2 <- min(X[,2]); max2 <- max(X[,2])
   x1 <- seq(min1 - 0.1*(max1 - min1), max1 + 0.1*(max1 - min1), length = 201)
   x2 <- seq(min2 - 0.1*(max2 - min2), max2 + 0.1*(max2 - min2), length = 201)
   Z <- matrix(0, ncol = 201, nrow = 201)
   for (i1 in 1:201) {
      for (i2 in 1:201) {
         Z[i1, i2] <- sumdist(c(x1[i1],x2[i2]), X)
      }
   }
   
   library(viridisLite)
   image(x1, x2, Z, col = viridis(100), asp = 1)
   points(X, pch = 19, col = "red")
}

# Nakreslime si ucelovu fciu pre maly aj velky dataset

X <- matrix(rnorm(8), ncol = 2)
plot.sumdist(X)
X <- matrix(rnorm(80), ncol = 2)
plot.sumdist(X)

# Zakladne procedury vo funkcii optim

optim(c(2, 2), sumdist, X = X, method = "Nelder-Mead")
system.time(for (i in 1:1000) res <- optim(c(2, 2), sumdist, 
                                           X = X, method = "Nelder-Mead"))

optim(c(2, 2), sumdist, X = X, method = "BFGS")  # Nemusi skonvergovat
system.time(for (i in 1:1000) res <- optim(c(2, 2), sumdist,
                                           X = X, method = "BFGS"))

optim(c(2, 2), sumdist, X = X, method = "SANN")  # Nehodi sa az tak pre tuto situaciu
system.time(for (i in 1:10) res <- optim(c(2, 2), sumdist,
                                         X = X, method = "SANN"))

# Vyskusajme PSO

library(pso); help(pso)
psoptim(c(2, 2), sumdist, X = X)
system.time(for (i in 1:10) res <- psoptim(c(2, 2), sumdist, X = X))

# Da sa pouzit vela inych metod, napriklad GA alebo DE

# Nakoniec skusme vypocet pomocou "specializovaneho" algoritmu
# https://en.wikipedia.org/wiki/Geometric_median#Computation
library(Gmedian); help(Gmedian)
Weiszfeld(X)
system.time(for (i in 1:1000) res <- Weiszfeld(X))