GA_basic <- function(lower, upper, fn, popSize = 100, maxIter = 100,
                     mutation_prob = 0.05,
                     initPopulation, crossover, mutation){
  #Create initial population
  population <- initPopulation(popSize, length(upper), lower, upper)
  for(iter in maxIter){
    #Calculate fitness of the initial population
    fitness <- apply(population, 1, fn)
    
    #Selection (for minimization) (Tournament)
    parents <- order(fitness)[1:2]
    par1 <- population[parents[1],]
    par2 <- population[parents[2],]
    
    #Apply crossover
    children <- crossover(par1, par2)
    
    #Apply mutation
    for (i in 1:popSize) {
      if (runif(1) < mutation_prob) {  # mutation probability
        population[i,] <- mutation(population[i,], lower, upper)
      }
    }
    
    #Update population
    toChange <- order(fitness)[(popSize-1):popSize]
    population[toChange[1],] <- children[1,]
    population[toChange[2],] <- children[2,]
  }
  return(list(population = population, sol = population[order(fitness)[1],], 
              fit = fitness, best_fit = fitness[order(fitness)[1]]))
}


initPop <- function(number, size, lower, upper){
  data <- runif(number*size, lower, upper)
  matrix(data, nrow = number, ncol = size, byrow = T)
}

my_crossover <- function(par1, par2) {
  n <- length(par1)
  sel <- c(rep(F, ceiling(n/2)), rep(T, floor(n/2)))
  sel <- sample(sel, n)
  child1 <- par1
  child2 <- par2
  child1[sel] <- par2[sel]
  child2[sel] <- par1[sel]
  matrix(c(child1, child2), nrow = 2, ncol = n, byrow = T)
}

#Apply a mutation in 1 dimension
my_mutation <- function(sol, lower, upper){
  sel <- sample(1:length(lower), 1)
  sol[sel] <- runif(1, lower[sel], upper[sel])
  sol
}



#run
GA_basic(rep(-5.12, 4), rep(5.12, 4), Rastrigin, popSize = 50, maxIter = 1000,
         mutation_prob = 0.05,
         initPopulation = initPop,
         crossover = my_crossover,
         mutation = my_mutation)