topt.gen <- function(b, N_pop, lambda_mut, N_gen) {
   
   # Konstrukcia t-optimalneho jednoducheho grafu pomocou GA.
   # <=> Konstrukcia D-optimalneho navrhu pre blokovy experiment s blokmi velkosti 2.
   # Pocet vrcholov (treatmentov) je fixne 10.
   #
   # Grafy su reprezentovane postupnostou 45 binarnych hodnot; ich suma je b.
   # "Rank proportional selection" najlepsich 50% jedincov.
   # "Rovnomerne simplexova" rekombinacia a "Vymenna" mutacia.
   
   # b ... pocet hran grafu (blokov experimentu). Musi platit (t-1)<=b<=(t nad 2).
   # lambda_mut ... intenzita mutacie v intervale (0,infty), obvykle 0.1 az 3.
   # N_pop ... velkost populacie, radovo desiatky
   # N_gen ... pocet generacii, radovo tisicky
   
   # Premenne na zaznamenanie priebehu vypoctu
   maxval <- medval <- rep(0, N_gen)
   date_start <- date()
   
   # Premenne na vykreslovanie grafov
   x <- y <- rep(0, 10)
   for (it in 1:10) {
      x[it] <- cos(2 * pi * (it - 1) / 10)
      y[it] <- sin(2 * pi * (it - 1) / 10)
   }
   
   # Vektory potrebne na kodovanie grafu (rx[i], gx[i]), i=1,...,45, su vsetky potencialne hrany
   rx <- c(rep(1, 9), rep(2, 8), rep(3, 7), rep(4, 6), rep(5, 5), 6, 6, 6, 6, 7, 7, 7, 8, 8, 9)
   gx <- c(2:10, 3:10, 4:10, 5:10, 6:10, 7:10, 8, 9, 10, 9, 10, 10)
   
   D <- function(xi)
   {
      # Kriterium "D-optimality", ktore rata pocet kostier zadaneho grafu
      # Funkcia fitness v terminologii GA
      # rt ... vektor indexov lavych vrcholov hran
      # gt ... vektor indexov pravych vrcholov hran
      
      rt <- rx[(1:45)[as.logical(xi)]]
      gt <- gx[(1:45)[as.logical(xi)]]
      b_act <- length(rt)
      M <- matrix(0, ncol = 10, nrow = 10)
      
      # Vypocet Laplacianu grafu
      for (ib in 1:b_act) {
         r <- rt[ib]
         g <- gt[ib]
         M[r, r] <- M[r, r] + 1
         M[g, g] <- M[g, g] + 1
         M[g, r] <- M[r, g] <- M[r, g] - 1
      }
      
      # Vypocet poctu kostier
      return(det(M[1:9, 1:9]))
   }
   
   # Vytvorenie pociatocnej populacie
   pop <- matrix(0, nrow = 2*N_pop, ncol = 45)
   Dvals <- rep(0, 2*N_pop)
   for (j in 1:(2*N_pop)) {
      pop[j, sample(1:45, b)] <- 1
   }
   
   # perc_old je len na priebezne info
   perc_old <- (-1)
   
   # Hlavny cyklus
   for (i in 1:N_gen) {
      
      # Informacia o tom, kolko percent vypoctu je uz dokoncenych.
      perc <- floor(100*i/N_gen) 
      if (perc > perc_old) {
         print(perc)
         perc_old <- perc
      }
      
      # Usporiadanie zlucenej populacie podla kriteria poctu kostier
      for (j in 1:(2*N_pop)) {
         Dvals[j] <- D(pop[j,])
      }
      
      ord <- order(Dvals, decreasing = TRUE)
      pop <- pop[ord,]
      Dvals <- Dvals[ord]
      maxval[i] <- max(Dvals)
      medval[i] <- median(Dvals)
      
      # Vykreslenie grafu "zastupenia genov" (prekreslenie kazdych 50 iteracii)
      # a momentalne najlepsieho grafu
      if (i %% 50 == 0) {
         plot(x, y, type = "n", main = c(i, round(maxval[i])))
         propnt <- apply(pop, 2, mean)
         for (j in 1:45) {
            lines(c(x[rx[j]], x[gx[j]]), c(y[rx[j]], y[gx[j]]),
                  col = rgb(1 - propnt[j], 1, 1 - propnt[j]), lwd = 4)
         }
         
         best <- which.max(Dvals)
         best_r <- rx[(1:45)[as.logical(pop[best,])]]
         best_g <- gx[(1:45)[as.logical(pop[best,])]]
         for (j in 1:b) {
            lines(c(x[best_r[j]], x[best_g[j]]), c(y[best_r[j]], y[best_g[j]]),
                  col = "red", lwd = 2)
         }
         
      }
      
      # Vygenerovanie novych N_pop jedincov z najlepsich N_pop.
      for (j in 1:N_pop) {
         
         # Selekcia dvojice pomocou rank proportional selection z prvych N_pop
         md <- sample(1:N_pop, 2, prob = N_pop:1)
         mum <- pop[md[1],]
         dad <- pop[md[2],]
         
         # "Rovnomerna simplexova" rekombinacia dvojice jedincov na potomka
         diff <- (1:45)[mum != dad]
         child <- mum
         if (length(diff) > 0) {
            child[diff] <- 0
            child[sample(diff, length(diff)/2)] <- 1
         }
         
         # "Vymenna" mutacia potomka
         nex <- rpois(1, lambda = lambda_mut)
         if (nex > 0) {
            for (k in 1:nex) {
               del <- sample((1:45)[as.logical(child)], 1)
               add <- sample((1:45)[as.logical(1 - child)], 1)
               child[del] <- 0
               child[add] <- 1
            }
         }
         
         # Pridanie potomka do zlucenej populacie
         pop[N_pop + j,] <- child
      }
   }
   
   windows()
   plot(maxval, type = "l", col = "red", ylim = c(0, max(maxval)))
   lines(medval, col = "green")
   
   # Najdenie najlepsieho a jeho navrat z procedury
   for (j in 1:(2*N_pop)) Dvals[j] <- D(pop[j,])
   best <- which.max(Dvals)
   best_r <- rx[(1:45)[as.logical(pop[best,])]]
   best_g <- gx[(1:45)[as.logical(pop[best,])]]
   
   return(list(best = rbind(best_r, best_g), val = Dvals[best],
               pars = c(N_pop, lambda_mut, N_gen), 
               date_start = date_start, date_end = date()))
}

# Demo
res <- topt.gen(10, 20, 1, 300)

library(igraph)
windows()
ig <- graph_from_edgelist(matrix(t(res$best), ncol = 2), directed = FALSE)
coords <- layout_(ig, with_gem())
plot.igraph(ig, layout = coords, main = "gem")

res <- topt.gen(15, 40, 2, 1000)

ig <- graph_from_edgelist(matrix(t(res$best), ncol = 2), directed = FALSE)
coords <- layout_(ig, with_gem())
tkplot(ig, layout = coords, main = "gem", vertex.color = "orange", edge.width = 3)