OneMax - Toy Problem

knitr::opts_chunk$set(echo = FALSE)

library(ggplot2)
library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

#### set up mutation 1/genome_size ####
mutate <- function(x) {
  p <- 1/(genome_size)
  bit_flip <- rbinom(1, 1, prob = p)
  x <- (x+bit_flip)%%2
}

#### set up reproduction with mutation ####
gametogeneis <- function(x) {
  gamete <-c()
  for (i in 1:nrow(x)) {
    gamete <- rbind(gamete, sapply(x[i, ], FUN = mutate))
  }
  return(gamete)  
}

#### generate 10 gametes for each individual
#### and randomly select from individual gametes ####
neut_select <- function(x) {
  gamete <- list()
  for (i in 1:gam_per_ind) {
    gamete[[i]] <- (gametogeneis(x))
  }
  selection <- c()
  for (k in 1:pop_size) {
    selection <- rbind(selection,
    gamete[[sample(1:gam_per_ind,1,T)]][sample(1:pop_size,1,T),]) #randomly select gamete for each individual
  }
  return(selection)
} 

#### setup population, pop size, with genome size ####
#### consisting bits (0,1) ####
populate_binom <- function(pop_size, genome_size) {
  N <- c()
  for (i in 1:pop_size) {
    N <- rbind(N,(rbinom(n = genome_size, size = 1, prob = 0.5)))
  } 
  return(N)
}
### alternative ###
populate_sample <- function(pop_size, genome_size) {
  N <- c()
  for (i in 1:pop_size) {
    N <- rbind(N,(sample(x = 0:1, size = genome_size, prob = c(0.5, 0.5), replace = T)))
  } 
  return(N)
}

#### MAIN ####
genome_size <- 20
pop_size <- 100
gam_per_ind <- 10
generation <- 200
pop_test <- populate_sample(pop_size, genome_size)

rowSums(pop_test) # check fitness of initial population

##   [1] 11  8 10 10  7  8 11 11 11  8 10 11 12  9  9 14 11 10 12 10  9 11  6
##  [24] 13 10 10 12 10  5  9 10 11 11 12 11 10 10 15  9  9 10  8  8 11  8  9
##  [47] 14  7 11 11  7 10 13  9 14  9  8 16  9  8  6 12  7  8  9 10 16 10  8
##  [70]  8 11 12 13 13  8 11  9  8  9  7  5  8 14 11 12 10  9 10 10 10 10 12
##  [93] 10 11  9 10 11 11 11  5

out_neutral <- data.frame(gen = 0, mean = mean(rowSums(pop_test))/genome_size, 
           min = min(rowSums(pop_test))/genome_size, 
           max = max(rowSums(pop_test))/genome_size)

##### preview fitness ####
out_neutral

##   gen mean  min max
## 1   0  0.5 0.25 0.8

which(rowSums(pop_test)/genome_size ==  max(rowSums(pop_test)/genome_size))

## [1] 58 67

pop_test[which(rowSums(pop_test)/genome_size ==  max(rowSums(pop_test)/genome_size)),]

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,]    1    0    1    1    1    1    1    0    1     1     0     1     1
## [2,]    1    1    1    1    1    1    1    1    1     1     0     0     1
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## [1,]     1     1     0     1     1     1     1
## [2,]     0     1     1     1     0     1     1

which(rowSums(pop_test)/genome_size ==  min(rowSums(pop_test)/genome_size))

## [1]  29  81 100

pop_test[which(rowSums(pop_test)/genome_size ==  min(rowSums(pop_test)/genome_size)),]

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,]    0    0    1    0    0    0    0    1    1     0     0     1     0
## [2,]    0    0    0    0    0    1    1    0    0     1     0     0     0
## [3,]    1    0    0    1    0    0    0    0    1     0     0     0     0
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## [1,]     0     0     0     1     0     0     0
## [2,]     0     1     0     0     1     0     0
## [3,]     1     0     0     0     1     0     0

counter <- 0 
for (i in 1:generation) {
  pop_test <- neut_select(pop_test)
  counter <- counter + 1
  out_neutral <- rbind(out_neutral, data.frame(gen = counter, 
                       mean = mean(rowSums(pop_test))/genome_size, 
                       min = min(rowSums(pop_test))/genome_size, 
                       max = max(rowSums(pop_test))/genome_size))

}

head(out_neutral)

##   gen   mean  min max
## 1   0 0.5000 0.25 0.8
## 2   1 0.5030 0.25 0.8
## 3   2 0.4985 0.20 0.8
## 4   3 0.4900 0.20 0.8
## 5   4 0.4965 0.25 0.8
## 6   5 0.4905 0.30 0.8

plot_neutral <- ggplot(out_neutral, aes(x = gen)) + 
  geom_line(aes(y = mean, colour = "mean")) + 
  geom_line(aes(y = max, colour = "max")) + 
  geom_line(aes(y = min, colour = "min")) +
  coord_cartesian(ylim = c(0:10/10)) +
  ggtitle("neutral") + xlab("generation") + ylab("fitness")
ggplotly(plot_neutral)