install.packages("GA")
library(GA)

#Rastrigin ---------------------------------------------------------------------
#We will try to find the minimum value of Rastrigin function
Rastrigin <- function(x) {
  sum(x^2 - 10 * cos(2 * pi * x)) + 10 * length(x)
}

#Define 30 dimensional Rastrigin function
set.seed(1234)
dimensions <- 30
#Set the domain bounds
lower <- rep(-5.12, dimensions)
upper <- rep(5.12, dimensions)

res <- ga(type = "real-valued",
         fitness = function(x) -Rastrigin(x),
         lower = lower,
         upper = upper,
         popSize = 50,
         maxiter = 1000,
         run = 100,
         pmutation = 0.04)

summary(res)
plot(res)

#We can see the function used for mutation/selection/crossover
defaultControl <- gaControl()
print(defaultControl)

#custom population/selection/crossover/mutation---------------------------------

#Simple random initilization
initializePopulation <- function(object, ...) {
  n <- object@nBits
  if(is.na(n)){
    n <- length(object@lower)
  }
  data <- runif(n*object@popSize, object@lower, object@upper, byrow = T)
  matrix(data, nrow = object@popSize, ncol = n)
}

#A mutation that adds noise to every dimension
mutationUniform <- function(object, parent, ...) {
  n <- object@nBits
  if(is.na(n)){
    n <- length(object@lower)
  }
  mut <- object@population[parent,]
  mut <- mut + runif(n, -1, 1) #Add noise
  mut <- pmin(mut, object@upper) #Keep solution inside bounds
  mut <- pmax(mut, object@lower)
  return(mut)
}

#Apply crossover where you randomly crossover over each dimension
crossover5050 <- function(object, parents, ...) {
  n <- object@nBits
  if(is.na(n)){
    n <- length(object@lower)
  }
  sel <- c(rep(F, ceiling(n/2)), rep(T, floor(n/2)))
  sel <- sample(sel, n)
  par1 <- object@population[parents[1],]
  par2 <- object@population[parents[2],]
  child1 <- par1
  child2 <- par2
  child1[sel] <- par2[sel]
  child2[sel] <- par1[sel]
  children <- matrix(c(child1, child2), nrow = 2, ncol = n, byrow = T)
  return(list(children = children, fitness = c(NA, NA) ))
}


set.seed(1234)
res <- ga(type = "real-valued",
          fitness = function(x) -Rastrigin(x),
          lower = lower,
          upper = upper,
          popSize = 10,
          maxiter = 1000,
          run = 100,
          pmutation = 0.4,
          population = initializePopulation,
          mutation = mutationUniform,
          crossover = crossover5050)

summary(res)
plot(res)

#------------TSP----------------------------------------------------------------
eurodistmat <- as.matrix(eurodist)

#A function for drawing routes
drawRoute <- function(sq, distances){
  #find the 2d coordinates
  loc <- -cmdscale(distances, add = TRUE)$points #PCA
  x <- loc[,1]; y <- loc[,2]
  s <- seq_len(nrow(distances))
  tspinit <- loc[c(sq, sq[1]),] #create a cycle
  #draw the cities and the inital solution
  plot(x, y, type = "n", asp = 1, xlab = "", ylab = "",
       main = "TSP route", axes = FALSE)
  arrows(tspinit[s,1], tspinit[s,2], tspinit[s+1,1], tspinit[s+1,2],
         angle = 10, col = "red")
  # Loop through the sequence and draw arrows with distances
  for (i in s) {
    j <- i + 1
    if(j > length(s)){
      j <- 1
    }
    xi <- tspinit[i,1]; yi <- tspinit[i,2]
    xj <- tspinit[j,1]; yj <- tspinit[j,2]
    arrows(xi, yi, xj, yj, angle = 10, col = "red")
    
    # Calculate and draw the distance on each arrow
    dist <- distances[sq[i], sq[j %% (nrow(distances) + 1)]]
    text((xi + xj) / 2, (yi + yj) / 2, labels = round(dist, 1), cex = 0.8, col = "blue")
  }
  text(x, y, labels(eurodist), cex = 0.8)
}


#Draw the starting solution
drawRoute(1:21, eurodistmat)

#Basic local search for TSP

# Define the objective function
objective_function <- function(solution) {
  # Evaluate the solution and return its objective value
  indexes <- embed(c(solution, solution[1]), 2)
  -sum(eurodistmat[cbind(indexes[,2], indexes[,1])])
}

res <- ga(type = "permutation",
          fitness = objective_function,
          nBits = 21,  # The length of the permutations
          lower = rep(1, 21),
          upper = rep(21, 21),
          popSize = 50,
          maxiter = 1000,
          run = 50)
res
summary(res)
plot(res)
#Drawsolution
drawRoute(res@solution, eurodistmat)

#0-1 Backpack problem-----------------------------------------------------------
# Size of the problem
n <- 50
# Item prices
prices <- runif(n, 10, 40)
# Item weights
weights <- runif(n, 5, 10)
# Maximum weight
maxWeight <- sum(weights)*0.6

# Define the objective function
objective_function <- function(solution) {
  if(sum(solution*weights) <= maxWeight){
    return(sum(solution*prices))
  }else{
    return(-sum(solution*weights)) 
  }
}

res <- ga(type = "binary",
          fitness = objective_function,
          nBits = n,  # The length of the binary vector
          lower = rep(1, 21),
          upper = rep(21, 21),
          popSize = 50,
          maxiter = 1000,
          run = 50)
res

sum(res@solution*prices)
maxWeight
sum(res@solution*weights)

summary(res)
plot(res)