from random import random
from collections import Counter

class Creature:

    RECESSIVE = "a"
    DOMINANT  = "A"

    def __init__(self, g1, g2, parent1 = None, parent2 = None):
        self.gene1 = g1
        self.gene2 = g2
        self.parent1 = parent1
        self.parent2 = parent2

    def breedWith(self, other):
        return (self.__class__)(self.randomGene(), other.randomGene(), self, other)

    def randomGene(self):
        if random() > 0.5:
            return self.gene2
        else:
            return self.gene1

    def genomeAsString(self):
        return self.gene1 + self.gene2

    def phenotype(self):
        if self.gene1 == self.RECESSIVE and self.gene2 == self.RECESSIVE:
            return self.RECESSIVE
        else:
            return self.DOMINANT

class AvoidSiblingsCreature(Creature):

    def chooseMate(self, population):
        mate = None
        n = 1.0
        for x in population:
            if not self.isSibling(x):
                if random() <= 1/n:
                    mate = x
                    n += 1
        return mate

    def isSibling(self, other):
        if other.parent1 == self.parent1 != None:
            return True
        if other.parent1 == self.parent2 != None:
            return True
        if other.parent2 == self.parent1 != None:
            return True
        if other.parent2 == self.parent2 != None:
            return True
        return False

class SeekSiblingsCreature(Creature):

    # Uniformly picks a sibling, but uniformly picks anyone if that's not possible
    def chooseMate(self, population):
        mate = None
        secondChoice = None
        n = 1.0
        k = 1.0
        for x in population:
            if self.isSibling(x):
                if random() <= 1/n:
                    mate = x
                    n += 1
            else:
                if random() <= 1/k:
                    secondChoice = x
                    k += 1
        if mate != None:
            return mate
        else:
            return secondChoice

    def isHalfSibling(self, other):
        if other.parent1 == self.parent1 != None:
            return True
        if other.parent1 == self.parent2 != None:
            return True
        if other.parent2 == self.parent1 != None:
            return True
        if other.parent2 == self.parent2 != None:
            return True
        return False

    def isSibling(self, other):
        if other.parent1 == self.parent1 != None:
            if other.parent2 == self.parent2 != None:
                return True
        if other.parent1 == self.parent2 != None:
            if other.parent2 == self.parent1 != None:
                return True
        return False

class Population:

    DEATH_RATE_a = 0.1
    DEATH_RATE_A = 0.05

    BIRTH_RATE = 0.068/3
    CHILDREN_PER_COUPLE = 3

    def __init__(self, population):
        self.population = population

    def breed(self):
        # Death
        toKill = []
        for x in self.population:
            if x.phenotype() == Creature.DOMINANT:
                if random() < self.DEATH_RATE_A:
                    toKill.append(x)
            else:
                if random() < self.DEATH_RATE_a:
                    toKill.append(x)
        for x in toKill:
            self.population.remove(x)

        # Birth
        newborns = []
        for x in self.population:
            if random() < self.BIRTH_RATE:
                y = x.chooseMate(self.population)
                n = 0
                while n < self.CHILDREN_PER_COUPLE:
                    newborns.append(x.breedWith(y))
                    n += 1
        self.population = self.population + newborns

    def stats(self):
        return Counter(map(lambda x: x.genomeAsString(), self.population))

    def size(self):
        return len(self.population)

def randomGene(p = 0.5):
    if random() < p:
        return Creature.DOMINANT
    else:
        return Creature.RECESSIVE

N = 100
initialPopulation = []
while len(initialPopulation) < N:
    g1 = randomGene()
    g2 = randomGene()
    x = SeekSiblingsCreature(g1, g2)
    initialPopulation.append(x)

P = Population(initialPopulation)
maxGenerations = 100
for n in range(maxGenerations):
    print P.size(), P.stats()
    P.breed()

print "Final frequency of recessive phenotype: %f" % (P.stats()["aa"]/float(P.size()))