mutate <- function(x){
    to.mut <- runif(1);
    ifelse(to.mut <= 1/genome_size, ifelse(x==0,1,0), x);
}

reproduce <- function(x) {
  b <- c();
  for (i in 1:pop_size) {
    b<- rbind(b, sapply(x[i,], FUN = mutate));
  }
  return(b); 
}

gametogenesis<-function(x){
  genepool<-list();
  for (i in 1:gametes) {
    genepool[[i]]<-(reproduce(x))
  }
  nextgen<-c();
  for (i in 1:pop_size) { 
   nextgen<-rbind(nextgen, genepool[[sample(1:gametes,1,replace = T)]][sample(1:pop_size,1,replace = T),]) 
  }
  return(nextgen);
}


# MAIN
pop_size <- 10;
genome_size <- 20;
mut<- 1/genome_size;
gametes<-10;
gen<-200;
# create population
lemons <- c();
for (i in 1:pop_size) {
  lemons<-rbind(lemons, (rbinom(genome_size,1,0.5)));
}
# create fitness record
fitness <- data.frame();
fitness <- data.frame(gen = 0, min = min(rowSums(lemons)/genome_size), 
                      max = max(rowSums(lemons)/genome_size), 
                      mean = mean(rowSums(lemons)/genome_size));
counter<-0;
for (j in 1:gen){
  lemons <- gametogenesis(lemons);
  counter <- counter+1;
  fitness<-rbind(fitness,data.frame(gen = counter, min = min(rowSums(lemons)/genome_size), 
                               max = max(rowSums(lemons)/genome_size), 
                               mean = mean(rowSums(lemons)/genome_size)))
}
#PLOT 

plot(fitness[,1], fitness[,4], type = "l", ylim = c(0,1), main = "neutral", ylab = "fitness", xlab = "generation");
lines(fitness[,1], fitness[,2], col="red");
lines(fitness[,1], fitness[,3], col="green");