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 ####
# function only mutates a single bit 
# if x is a vector with n>1 elements
# it will treat the whole vector as a bit
# it will not work 
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]  7  7 13 10  9 11 12 11 12  8  9 13 11 13  8 12 10  6  7 10 11 13 11
##  [24] 10 13  9  9 10 10 11 14 12 12  9 13 10  9 12  8 11  8 10 12 10 13 14
##  [47] 12  6 11  6  7  8 11 10 10  7 16 10  9 11 13 11  6 14 12 10 13 12 11
##  [70] 11  7 10 10 12  8 11 14 12  8  9 11  8  9  9 12  8 13  9  9 10 13 11
##  [93] 12  8 11 13 11 14  9 12

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

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

## [1] 57

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

##  [1] 0 1 0 1 1 1 1 1 1 0 1 1 1 0 1 1 1 1 1 1

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

## [1] 18 48 50 63

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,]    1    0    0    1    1    0    1    0    0     0     0     0     0
## [2,]    0    0    0    1    0    1    0    0    0     1     0     1     0
## [3,]    0    0    0    0    0    0    0    0    1     0     0     1     0
## [4,]    0    0    0    0    0    0    1    0    1     0     1     0     0
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## [1,]     0     1     0     0     0     0     1
## [2,]     0     0     1     0     1     0     0
## [3,]     0     1     1     0     1     0     1
## [4,]     1     1     0     0     0     1     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)

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)