from random import randint

def create_population(create, pop_size=1000):
    """Create population list"""
    return [create() for i in xrange(pop_size)]

def crossover(combine, population, pop_size):
    """Increase population size by combining"""
    original_size = len(population)
    while len(population) < pop_size:
        population.append(combine(population[randint(0, original_size-1)], population[randint(0, original_size-1)]))
    return population

def mutation(mutate, population, mutation_rate):
    """Mutate portion of population"""
    for i in xrange(int(len(population) * mutation_rate)):
        pos = randint(0, len(population)-1)
        population[pos] = mutate(population[pos])
    return population

def selection(fitness, population):
    """Keep best in population"""
    return sorted(population, key=lambda x: fitness(x), reverse=False)[:len(population)/2]

def main():
    """Goal: identify best approximation of pi"""

    pop_size = 1000       # Number of solutions
    mutation_rate = 0.05  # Fraction of population to mutate
    epochs = 20           # Number of iterations to run for

    # Create randomised population
    population = create_population(lambda: [randint(1,5000), randint(1,5000)], pop_size)

    # Run for number of iterations
    for epoch in xrange(epochs):

        # Remove worst solutions
        population = selection(lambda x: abs(float(x[0])/(0.0001+x[1]) - 3.14159265), population)

        # Show the best so far
        best = float(population[0][0]) / (0.0001+population[0][1])
        error = abs(best - 3.14159265)
        print('Epoch %.03d: %f\t(%d / %d)\terror %f' % (epoch, best, population[0][0], population[0][1], error))

        # Combine solutions to increase population size to original
        population = crossover(lambda x, y: [(x[0]+y[0])/2, (x[1]+y[1])/2], population, pop_size)

        # Mutate some solutions
        population = mutation(lambda x: [x[0]+randint(-1,+1), x[1]+randint(-1,+1)], population, mutation_rate)

if __name__ == '__main__':
    main()