library(ggplot2)
library(gganimate)
library(grDevices)

rand.pupulation <- function(n) {
    return(matrix(runif(2*n), ncol=2)) 
}

alpha.fixed <- function(alpha) {
    return(function() alpha)
}

alpha.runif <- function() {
    return(function() runif(1))
}

recom.singlearithm <- function(alpha_gen) {

    return(function(parents) {
        alpha <- alpha_gen()
        gene_index <- sample(1:ncol(parents), 1)

        gene_p1 <- parents[1, gene_index]
        gene_p2 <- parents[2, gene_index]

        children <- parents

        children[1,gene_index] <- alpha * gene_p2 + (1 - alpha) * gene_p1
        children[2,gene_index] <- alpha * gene_p1 + (1 - alpha) * gene_p2

        return(children)
    })
}

recom.wholearithm <- function(alpha_gen) {
    return(function(parents) {
        alpha <- alpha_gen()

        children <- parents

        children[1,] <- alpha * children[1,] + (1 - alpha) * children[2,]
        children[2,] <- alpha * children[2,] + (1 - alpha) * children[1,]

        return(children)
    })
}

next.population <- function(population, recom) {
    next_population <- matrix(, nrow=0, ncol=2)

    parent_pairs <- matrix(sample(1:nrow(population), nrow(population)), ncol=2)

    for(i in 1:nrow(parent_pairs)) {
        parents <- parent_pairs[i,]
        next_population <- rbind(next_population, recom(population[parents,]))
    }

    return(next_population)
}

experiment <- function(population, gen, recom) {
    df <- data.frame( x = population[,1]
                     ,y = population[,2]
                     ,generation = 0)

    for(g in 1:gen) {
        population <- next.population(population, recom)
        
        df <- rbind(df, data.frame( x = population[,1]
                                   ,y = population[,2]
                                   ,generation = g))
    }

    return(df)
}

plot.experiment <- function(df, filename) {
    pdf(file=filename, onefile=TRUE)

    for(g in unique(df$generation)) {
        tmp_df <- df[df$generation == g,]

        p <- ggplot(data=tmp_df, aes(x=x, y=y)) +
            geom_point() +
            labs(title=sprintf("generation: %d", g)) +
            xlim(0, 1) +
            ylim(0, 1)

        hull <- chull(tmp_df$x, tmp_df$y)

        p <- p + geom_polygon(data=tmp_df[hull,], alpha=0.25)

        print(p)
    }

    dev.off()
}

all.experiments <- function() {
    population <- rand.pupulation(20)
    
    df <- experiment(population, 10, recom.singlearithm(alpha.fixed(0.5)))
    plot.experiment(df, "single_fixed_alpha.pdf")

    df <- experiment(population, 10, recom.singlearithm(alpha.runif()))
    plot.experiment(df, "single_rand_alpha.pdf")

    df <- experiment(population, 10, recom.wholearithm(alpha.fixed(0.5)))
    plot.experiment(df, "whole_fixed_alpha.pdf")

    df <- experiment(population, 10, recom.wholearithm(alpha.runif()))
    plot.experiment(df, "whole_rand_alpha.pdf")
}