#!/usr/bin/env python
from __future__ import division

import random
from math import log

# How many (binary) genes does each organism have?
genecount = 128 
# How many generations to interate for?
num_gens = 20
# Chance that an individual gene will mutate this generation.
# Note that there's a 50% chance the mutation will revert to its current value
mut_rate = 0.2
# Number of organisms in the population
pop_size = 100
# The lower, the larger the runs of genes provided by a single parent
# (For each gene, this is the chance of switching to the other parent as gene source)
crossover_rate = 0.1
# How succsesful are less fit members of the population at finding mates to have sex with?
sex_lambda = 6.0
# Show detailed symbols for each generation, or just on/off?
show_detailed = True
# Display population after each generation
show_population = True
# Display initial population
show_initial = True
# Display final population
show_final = True
# Show frequency summaries
show_summaries = True
# Start from a blank slate of no genes expressed, vs. a random set of genes
tabula_rasa = False
# Bad gene rate
bad_rate = 2
# Neutral gene rate
neutral_rate = 4
# Good gene rate
good_rate = 1

filter_choices = tuple(([-1] * bad_rate) + ([0] * neutral_rate) + ([1] * good_rate))

def make_sample():
	"""Generate an organism for our population"""
	if not tabula_rasa:
		return [random.randint(0,1) for i in xrange(genecount)]
	else:
		return [0] * genecount

def sex(a,b):
	"""Produce the offspring of two parents.
	   One parent is chosen at random, and start copying genes. At each step,
		 there is a chance of swapping to the other parent as a gene source."""
	source = random.randint(0,1)
	new = []
	for genes in zip(a,b):
		# 10% chance of swapping for each gene
		if random.random() < crossover_rate:
			source ^= 1
		new.append(genes[source])
	return new

def mutate(ind):
	"""Randomly mutate some fraction of the genes. See note on mut_rate earlier.
		 Only point mutations, since more productive mutations like translations,
		 duplication, etc. don't really make sense in this model."""
	return [random.randint(0,1) if random.random() < mut_rate else x for x in ind]

# -1 = selected against (i.e. bad mutation), 1 == selected for (good mutation)
# 0 = neutral.
filter = [random.choice(filter_choices) for i in xrange(genecount)]

def fitness(individual):
	"""Calculate the fitness of the individual organism. Compute the dot product
	   of the individual's genes and the environment."""
	return sum(map(lambda (x,y): x * y, zip(filter, individual)))

def symbol(x, y):
	"""Map filter/indivdual gene pairings to symbols to easily grasp successes of the organism.
		Change the '  ' in the 4th line to something else to see which neutral genes
    are being expressed."""
	if x == -1:
		return '+x'[y]
	if x == 0:
		return '  '[y]
	if x == 1:
		return '-o'[y]

def choose_partners(pop):
	"""Pick a pair of organisms from the population to have sex. For simplicity's 
     sake, they're hermaphroditic. Mates are chosen based on an exponential 
		 distribution. The value of sex_lambda determines the half-way point. For
     instance, if sex_lambda == 6.0, then organisms in the top 6th of fitness 
     would be a partner half of the time."""
	a = pop_size 
	while a >= pop_size:
		a = int(random.expovariate(sex_lambda / pop_size))
	b = pop_size
	while b >= pop_size or a == b:
	  b = int(random.expovariate(6.0 / pop_size))
	return pop[a],pop[b]
	
def H(x,r):
	if x:
		return -log(r,2)
	else:
		return -log(1-r,2)

def info(ind, genefreq):
	return sum((H(x, float(y) / pop_size) for (x,y) in zip(ind, genefreq)))

# Calculate our initial population, sorting it according to its fitness
init_gen = sorted([make_sample() for i in xrange(pop_size)], key=fitness, reverse=True)

# Print out our filter for reference (although it can be deduced fairly easily
# from the symbols we use to display the organisms)
print ''.join(['x o'[x+1] for x in filter])


# Frequency details
genefreq = [0] * genecount
for ind in init_gen:
	genefreq = [x+y for (x,y) in zip(genefreq, ind)]

if show_initial:
	# Print out the initial population
	for ind in init_gen:
		print ''.join([symbol(x,y) for (x,y) in zip(filter, ind)]), fitness(ind), info(ind, genefreq)

if show_summaries:
	print ''.join(['0123456789+'[int(x * 10.0 / pop_size)] for x in genefreq])

next_gen = init_gen

# Iterate over the remaining generations
for gen in xrange(num_gens):
	# Get the next generation by having the current one have sex with itself.
	next_gen = [sex(*choose_partners(next_gen)) for i in xrange(pop_size)]
	# Mutate the organisms
	next_gen = map(mutate, next_gen)
	# And sort them by fitness.
	next_gen = sorted(next_gen, key=fitness,reverse=True)

	# Nice little horizontal rule to make it easier to tell generations from
  # one another.
	print '-'* genecount

	# Frequency details
	genefreq = [0] * genecount
	for ind in next_gen:
		genefreq = [x+y for (x,y) in zip(genefreq, ind)]

	if show_population:
		# Display the population of this generation
		for ind in next_gen:
			print ''.join([symbol(x,y) for (x,y) in zip(filter, ind)]), fitness(ind), info(ind, genefreq)
	

	if show_summaries:
		print ''.join(['0123456789+'[int(x * 10.0 / pop_size)] for x in genefreq])

if show_final and not show_population:
	for ind in next_gen:
		print ''.join([symbol(x,y) for (x,y) in zip(filter, ind)]), fitness(ind)

totalrate = bad_rate + neutral_rate + good_rate
blind_success = 0
success = 0
for freq, actual in zip(genefreq, filter):
	blind_success += (bad_rate, neutral_rate, good_rate)[actual + 1] / totalrate
	if freq / pop_size < 0.35:
		guess = -1
	elif freq / pop_size < 0.7:
		guess = 0
	else:
		guess = 1
	success += (guess == actual) and 1 or 0

print blind_success, success