#---------------------------------------------------------------#
# optimisación swarm pso                                        #
library(pso)
library(ggplot2)

                                       # función a optimizar
d=2                                    # 2 dimensiones
f_ <- function(x1,x2){                 # función de Rastrigin
        d * 10+x1^2 + x2^2-10*(cos(2*pi*x1)+cos(2*pi*x2))
        }

f <- function(x = c(0,0)){             # su versión vectorial
        f_(x[1],x[2])
        } 
f(c(5,5))                              # ejemplo de uso

                                       # los límites
min_x1 <- - 5.2                        # clásicos de Rastrigin
min_x2 <- - 5.2
max_x1 <-   5.2
max_x2 <-   5.2


maxit=100   # num iteración max
s=40        # número de partículas

control =list(trace=1,      # modo traza
              maxit=maxit,
              REPORT=1,     # num de informes
              trace.stats=1,
              s = s
              )

set.seed(15)
                        # optimización
p <- psoptim(rep(NA,d), # tantos NA como d
             fn=f,
             lower=c(min_x1,min_x2),
             upper=c(max_x1,max_x2),
             control=control)

qplot(p$stats$it,p$stats$error,
      xlab= "iteration",ylab = "mínimo")

cat("Coord óptima:",round(p$par,2),
    "valeour de la función correspondiente",
     round(p$value,2),"\n")

#---------------------------------------------------------------#
x1_ <- seq(min_x1, max_x1, by=0.2) # grid x1
x2_ <- seq(min_x2, max_x2, by=0.2) # grid x2
y_  <- outer(x1_, x2_, f_)         # un valor por columna


paleta <- colorRamppaleta(c("red","yellow"))(100)
diag <- persp(x1_,                            # perspectiva 3D
              x2_,
              y_,
              phi=-10,                        # ángulos de presentación
              theta=60,
              col=paleta[round((y_+max(y_))/(max(y_)-min(y_))* 55)],          
              ticktype = "detailed",         # eje detallado
              shade=.001,
              cex.lab=1.5 )             

u <- as.vector(p$par[1])
v <- as.vector(p$par[2])
w <- as.vector(p$value)

nube <- trans3d(u,                       # nube de 1 puntos
                 v,
                 v,
                 pmat = diag)

points(nube,                            # inserción en diag
       pch = 1,      # símbolo
       cex = 8,      # tamaño
       lwd = 5,      # anchura línea
       col = 2       # color
       )

