complexify.py

# -*- coding: utf-8 -*-
import random

class Lineage:
  """ Model one independantly evolving lineage """
  def __init__(self, complexity= 0, min = None):
    self.complexity = complexity
    self.min = min
  
  def __repr__(self):
    return "Lineage(%s, %s, %s)" % (self.complexity, self.min)
  
  def generation(self, report = False):
    """ do the evolution (get more or less complex)"""
    self.complexity += random.normalvariate(0, 0.01)
    if min:
      if self.complexity < self.min:
        self.complexity = self.min
    if report:
      return "complexity:%s" % self.complexity
  
  def run(self, generations = 100, log=None):
    """ Run a simulation with one lineage over several generations"""
    for generation in range(generations):
      self.generation()
      if log:
        log.write("%s,%s\n" % (generation, self.complexity))
    print "complexity after %s generations: %s" % (generation, self.complexity)
    
    
class Simulation:
  """ Run a simulation with many lineages
  
  'Generations' sets the number of generations to run the simulation, "species"
  is the equlibrium  number of species to maintain in each population (the point
  at which the speciation rate and the extinction rate are balanced) and
  "turnover" sets both the extinction and speciation rate at that balance. The
  population of evoling lineages is stored as a tuple list ([id, Lineage())...]

  In order to get multiple species into the population early on, the speciation
  rate is elevated (5x) before the equilibrium point is met and extinction
  doesn't kick in untill half the equilibrium population is in place. 
  """
  
  def __init__(self, generations= 500, species = 50, turnover = 0.005, min= None):
    self.generations = generations
    self.species = species
    self.turnover = turnover
    self.min = min
    
    L = Lineage(min=self.min)
    self.lineage_count = 1
    self.population = [(1,L)]  
  
  def _speciate(self, rate):
    """ Create a new lineage that inherits its parents complexity """
    new_species = []
    for i, L in self.population:
      if random.random() < rate:
        self.lineage_count += 1
        new_species.append((self.lineage_count, Lineage(complexity=L.complexity, min=self.min)) )
    self.population += new_species
    
  def _cladogenesis(self):
    """ Kill of some lineges, duplicate others """
    if len(self.population) > self.species/2:
      self.population = [L for L in self.population if
                         random.random() > self.turnover]
    
    if len(self.population) < self.species:
      self._speciate(self.turnover * 5)
    else:
      self._speciate(self.turnover)
    
  def _csv_line(self, tup, generation):
    """ Converts information in a (lineage ID, lineage) tuple to string"""
    return "%s, %s, %s\n" % (tup[0], generation, tup[1].complexity)
  

  def generation(self, report= False):
    """ Evolve for one generation with extinction/speciation and complexity """
    self._cladogenesis()
    for i, L in self.population:
      L.generation()
    if report:
      print [L for i, L in self.population]
        
  def run(self, log, header = True):
    """ Run a simulation over many generations """
    if header:
      log.write("lineage, generation, complexity\n")
    for generation in range(self.generations):
      log.writelines([ self._csv_line(t, generation) for t in self.population])
      self.generation()
 
def main():
  s = Simulation(generations = 300, species = 20)
  s.run(open("test.csv", "w"))

if __name__ == "__main__":
   main()

plots.r

library(ggplot2)
theme_set( theme_bw() )

test <- read.csv("test.csv")
#take a look
p <- ggplot(test, aes(generation, complexity, group=as.factor(replicate)))
p + geom_line()

# which lineages get past 0.5 complexity?

max.complexity <- tapply(test$complexity, as.factor(test$lineage), max)
colour <- ifelse(max.complexity > 0.5, "red", "grey")
p2 <- ggplot(test, aes(generation, complexity, colour=as.factor(replicate)))
p2 + geom_line(size=0.2, alpha=0.3) + scale_colour_manual(values=colour)