'''
Code from 07-12-sui04

(You can skip to sections by searching for `<num> --- i.e. backtic
followed by section number.)

TABLE OF CONTENTS:
`1. Random Variables and Statistics
`2. Simulation and Population Management
`3. Decision-Makers (Agent Brain)
`4. Food
`5. Agents
`6. Simulation Running
`6-1. Parameters
`6-2. Experiment Management

'''
from random import *
import gc
import sys
import os
from copy import *
from types import *
import time
import string
import math
import glob
import re
try:
	import psyco
	psyco.full()
except:
	pass

################################################
#  `1. Random Variables and Statistics         #
################################################

def multinomial(pv):
	total = 0
	for prob in pv:
		if prob > 0: total += prob
	# In the unlikely event total=0, return an index uniform randomly
	if total==0:
		return randint(0, len(pv)-1)
	point = total * random() #Pick a random point along the range [0, total)
	total = 0
	for i,prob in enumerate(pv):
		if prob>0: total += prob
		if total > point:
			return i

def clipNegative(val):
	if val > 0: return val
	return 0
	
def normalizeMultinomial(pv):
	total = 0
	for prob in pv:
		if prob > 0: total += prob
	if total==0:
		return [1.0/len(pv) for i in pv]
	return [clipNegative(prob)/total for prob in pv]

def printMultinomial(pv, out = sys.stdout, showUnNormed = 0, noNewLine = 0):
	if not showUnNormed:
		total = 0
		for prob in pv:
			if prob > 0: total += prob
	print >>out, "PV[",
	for prob in pv:
		if showUnNormed:
			val = prob
		else:
			if prob > 0 and total > 0:
				val = prob / total
			else:
				val = 0
		print >>out, "%f," % val,
	if noNewLine:
		print >>out, "]",
	else:
		print >>out, "]"

# Allows statistics to be updated progressively, without having to store
# each data item
# EXAMPLE:
#   progStats = ProgressiveStats()
#   progStats.update(3)
#   progStats.update(7)
#   print progStats.calcMean()
#   -> '5'
class ProgressiveStats:
	lastValue = 0;
	n = 0
	sumof_x_ = 0
	sumof_xsquared_ = 0
	min = 0
	max = 0
	haveMinMax = 0
	
	def update(self, xVal):
		self.lastValue = xVal
		self.n += 1
		self.sumof_x_ += xVal
		self.sumof_xsquared_ += xVal ** 2
		self.confLevel = 2.58
		if not self.haveMinMax:
			self.min = xVal
			self.max = xVal
			self.haveMinMax = 1
		else:
			if xVal < self.min:
				self.min = xVal
			elif xVal > self.max:
				self.max = xVal

	def calcCount(self):
		return self.n

	def calcTotal(self):
		return self.sumof_x_

	def calcMean(self):
		if self.n > 0:
			return self.sumof_x_/self.n
		return 0

	def calcVar(self):
		if self.n > 0:
			return self.sumof_xsquared_/self.n - (self.sumof_x_/self.n) ** 2
		return 0

	def calcSd(self):
		## Can be negative if the quantities above get too big and thus inaccurate
		return self.calcVar<0 and 0 or sqrt(self.calcVar)

	def calcSe(self):
		if self.n > 0:
			return self.calcSd()/sqrt(self.n)
		return 0

################################################
#  `2. Simulation and Population Management    #
################################################
# A set of agents or food pieces, etc.
# pyClass is the name of the entity class (e.g. Agent)
# EXAMPLE:
#   pop = Population('agents', Agent)
class Population:
	def __init__(self, name, pyClass):
		self.name = name
		self.pyClass = pyClass
		self.members = []
	def size(self):
		return len(self.members)
	def generate(self, number):
		newMembers = []
		for i in xrange(number):
			newMember = self.pyClass()
			newMember.generate()
			self.members.append(newMember)
			newMembers.append(newMember)
		return newMembers

# The n x n grid on which agents, food, etc. live
class Board:
	def __init__(self):
		self.width = pBoardWidth
		self.height = pBoardHeight
		self.board = [[None for y in xrange(pBoardHeight)] for x in xrange(pBoardWidth)]
	def placeAnywhereWithin(self, entity, refEntity, width, height):
		(startx, starty, endx, endy) = (refEntity.x - width/2, refEntity.y - height/2, refEntity.x + width/2, refEntity.y + height/2)
		if startx<0: startx=0
		if starty<0: starty=0
		if endx>=pBoardWidth: endx = self.width - 1
		if endy>=pBoardHeight: endy = self.height - 1
		tries = 0
		while 1:
			(x,y) = (randint(startx, endx), randint(starty, endy))
			tries += 1
			if not self.board[x][y]:
				entity.x = x
				entity.y = y
				self.board[x][y] = entity
				return True
			if tries>(endx-startx)*(endy-starty): break
		return False
		
	def placeAnywhere(self, entity):
		tries = 0
		while 1:
			(x,y) = (randint(0, self.width-1), randint(0,self.height-1))
			tries += 1
			if not self.board[x][y]:
				entity.x = x
				entity.y = y
				self.board[x][y] = entity
				return True
			if tries>10: break
		return False
	
	# This does not search exhaustively
	def getEntity(self, entity, width, height, pyClass, test = 0, testArgs = ()):
		(startx, starty, endx, endy) = (entity.x - width/2, entity.y - height/2, entity.x + width/2, entity.y + height/2)
		if startx<0: startx=0
		if starty<0: starty=0
		if endx>=pBoardWidth: endx = self.width - 1
		if endy>=pBoardHeight: endy = self.height - 1
		tries = 0
		while 1:
			(x,y) = (randint(startx, endx), randint(starty,endy))
			tries += 1
			if self.board[x][y] and isinstance(self.board[x][y], pyClass) and (not test or test(self.board[x][y], *testArgs)):
				return self.board[x][y]
			if tries>(endx-startx)*(endy-starty): break
		return False
		
	def moveBy(self, entity, x, y):
		(newX, newY) = (entity.x + x, entity.y + y)
		# If out of bounds, fail
		if newX < 0 or newY < 0 or newX >= pBoardWidth or newY >= pBoardHeight:
			return False
		# If something's already there, fail
		if self.board[newX][newY]:
			return False
		self.board[newX][newY] = entity
		self.board[entity.x][entity.y] = None
		(entity.x, entity.y) = (newX, newY)
		return True
	def display(self):
		for y in xrange(pBoardHeight):
			for x in xrange(pBoardWidth):
				print self.board[x][y] and self.board[x][y].boardLabel or "_",
			print ""

# A simulation that holds a board and population, and can be run
# EXAMPLE:
#   sim = Simulation()
#   sim.addPopulation('agents', Agent)
#   sim.run(cycles = 100)
#   -> runs for 100 cycles, outputs to stdout
class Simulation:
	cycle = 0
	def __init__(self):
		self.populations = []
		# Shouldn't actually be possible (because different pops can have same class), but for most sims works fine
		self.popClasses = {}
		self.popNames = {}
		self.board = Board()
	def addPopulation(self, name, pyClass):
		self.populations.append(Population(name, pyClass))
		self.popClasses[pyClass.__name__] = self.populations[-1]
		self.popNames[name] = self.populations[-1]
	def getPopulation(self, pyClassOrName):
		if isinstance(pyClassOrName, StringType):
			return self.popNames[pyClassOrName]
		else:
			return self.popClasses[pyClassOrName.__name__]
	def runCycle(self):
		for pop in self.populations:
			for member in pop.members:
				try:
					member.makeObservations(self)
				except:
					pass
			for member in pop.members:
				toRemove = member.run(self)
				if toRemove:
					for en in toRemove:
						remPop = self.getPopulation(en.__class__)
						remPop.members.remove(en)
						self.board.board[en.x][en.y] = None
						(en.x,en.y) = (-1,-1)
					
	def run(self, cycles):
		for c in xrange(cycles):
			self.runCycle()
			self.cycle+=1

###############################################
#  `3. Decision-Makers (Agent Brains)         #
###############################################
class DecisionNode:
	pass

class DecisionBranch(DecisionNode):
	def __init__(self):
		self.childNodes = []
		self.splitValues = [] # One less than the number of nodes
		self.attribute = None # The attribute to make a decision about
	def normalize(self):
		for node in self.childNodes:
			node.normalize()
	def generate(self, withNodes = 0): #Just the 2 child nodes for now
		self.attribute = randint(0, pNumObs-1)
		self.splitValues.append(random())
		#Use the nodes supplied, if given
		if withNodes:
			self.childNodes = withNodes
		else:
			for i in xrange(2):
				newChildNode = random() < pGen_DTreeBranch and DecisionBranch() or DecisionLeaf()
				newChildNode.generate()
				self.childNodes.append(newChildNode)
	def display(self, out = sys.stdout, dtree = 0, indent = 0, showUnNormed = False):
		for i in xrange(indent): print >>out, " ",
		print >>out, "[%s]: %f" % (gObsPosToName[self.attribute], self.splitValues[0])
		for i in xrange(len(self.childNodes)):
			self.childNodes[i].display(out, dtree, indent + 2, showUnNormed = showUnNormed)
	def crossover(self, other):
		# Can't do a crossover with a branch and leaf, so just copy self for the rest of the recursion
		other = other.__class__== DecisionLeaf and self or other
		# Branch may disappear, and become one of the child nodes
		if random() < pMutate_DTreeStructureBranch:
			(copyFrom, otherFrom) = random()<0.5 and (self,other) or (other,self)
			randChild = randint(0, len(self.childNodes)-1)
			return copyFrom.childNodes[randChild].crossover(otherFrom.childNodes[randChild])
		else:
			child = DecisionBranch()
			child.attribute = random() < pMutate_DTreeAttribute and randint(0, pNumObs-1) or self.attribute
			child.splitValues = [self.splitValues[0] + gauss(0, pMutate_DTreeSplitValue)]
			child.splitValues = map(lambda a: a > 1 and 1 or (a < 0 and 0 or a), child.splitValues)
			# At this stage, guaranteed to have 2 branches, with equal number of splits (because they're all 2 for now)
			for i in xrange(len(self.childNodes)):
				(copyFrom, otherFrom) = random()<0.5 and (self,other) or (other,self)
				child.childNodes.append(copyFrom.childNodes[i].crossover(otherFrom.childNodes[i]))
			return child
	def generateWithStructure(self, structure, pos = 0):
		try:
			self.attribute = int(structure[pos]); pos += 1
		except ValueError:
			obsName = structure[pos]
			#sys.exit()
			self.attribute = gObsNameToPos[obsName]; pos += 1
		self.splitValues.append(float(structure[pos])); pos += 1
		# pos+2 should be '('
		pos += 1
		for i in xrange(2):
			if structure[pos]=='L':
				self.childNodes.append(DecisionLeaf())
				self.childNodes[-1].generate()
				pos += 1
			else:
				self.childNodes.append(DecisionBranch())
				pos = self.childNodes[-1].generateWithStructure(structure, pos)
		return pos + 1

class DecisionLeaf(DecisionNode):
	def __init__(self):
		self.pv = []
		self.called = 0
	def normalize(self):
		self.pv = normalizeMultinomial(self.pv)
	def generate(self):
		if (pGen_ActProbs):
			for prob in pGen_ActProbs:
				newProb = prob + gauss(0, pGen_ActProbSd)
				if newProb < 0: newProb = 0
				elif newProb > 1: newProb = 1
				self.pv.append(newProb)
		else:
			for i in xrange(pNumActs):
				self.pv.append(random())
	def display(self, out = sys.stdout, dtree = 0, indent = 0, showUnNormed = False):
		for i in xrange(indent): print >>out, " ",
		printMultinomial(self.pv, out, showUnNormed = showUnNormed, noNewLine = True)
		if dtree:
			if not dtree.dtreeCalled:
				print >>out, " {-}"
			else:
				print >>out, " {%.3f / %d}" % (float(self.called) / dtree.dtreeCalled, self.called)
	def reproduce(self):
		child = DecisionLeaf()
		for i, prob in enumerate(self.pv):
			newProb = prob + gauss(0, pMutateActProb)
			if newProb < 0: newProb = 0
			elif newProb > 1: newProb = 1
			child.pv.append(newProb)
		return child
	def crossover(self, other):
		#Leaf may split into a branch
		if random() < pMutate_DTreeStructureLeaf:
			child = DecisionBranch()
			child.generate([self.reproduce(),self.reproduce()])
			return child
		else:
			return self.reproduce()

# A decision tree that currently allows 2 arms per branch, and holds
# multinomial probability distributions at the leaves
# EXAMPLE:
#   dtree = DecisionTree()
#   dtree.generate()
#   print dtree
#   -> prints randomly generated decision tree
# EXAMPLE 2:
#   dtree = DecisionTree()
#   dtree.generateWithStructure("Age 0.35 ( L Energy 0.1 ( L L ) )")
#   print dtree
#   -> shows decision tree with 'Age' as root branch node, and 'Energy' as one
#      of the child nodes. The probability distributions in the leaves are randomly
#      generated
class DecisionTree:
	def __init__(self):
		self.rootNode = None
		self.dtreeCalled = 0
	def normalize(self):
		self.rootNode.normalize()
	def generate(self):
		self.rootNode = DecisionBranch()
		self.rootNode.generate()
	def display(self, out = sys.stdout, showUnNormed = False):
		print >>out, "{\n'Call Count': ",self.dtreeCalled
		self.rootNode.display(out, self, showUnNormed = showUnNormed)
		print >>out, "}\n"
	def crossover(self, other):
		child = DecisionTree()
		(copyFrom, otherFrom) = random()<0.5 and (self,other) or (other,self)
		child.rootNode = copyFrom.rootNode.crossover(otherFrom.rootNode)
		return child
	def decide(self, obs):
		self.dtreeCalled += 1
		node = self.rootNode
		i = 0
		while i<500:
			i+=1
			if node.__class__ == DecisionBranch:
				if obs[node.attribute] < node.splitValues[0]:
					node = node.childNodes[0]
				else:
					node = node.childNodes[1]
			else:
				node.called += 1
				return multinomial(node.pv)
		return 0
	def generateWithStructure(self, structure):
		structure = string.split(structure)
		self.rootNode = DecisionBranch()
		self.rootNode.generateWithStructure(structure)

###############################################
#  `4. Food                                   #
###############################################

# A piece of food that holds a certain amount of energy that can be consumed
# by agents
class Food:
	energy = 0
	birthCycle = 0
	def generate(self):
		self.energy = gauss(pFoodEnergy['mean'], pFoodEnergy['sd'])
	def run(self, sim):
		if sim.cycle - self.birthCycle > pFoodMaxAge:
			return [self]
		
Food.boardLabel = "F"

###############################################
#  `5. Agents                                 #
###############################################

# The main class of interest. The actions and observations that agents can
# do and make are defined below
# EXAMPLE:
#   agent = Agent()
#   agent.generate()
#   agent.run()
#   -> Runs agent for 1 cycle, which involves observing world, choosing action
#      and then performing it
#   agent.obsAge()
#   -> Returns the agent's observation of its age
class Agent:
	birthCycle = 0
	energy = 0
	lastActCycle = 0
	children = None
	obs = None
	actToDo = None
	def __init__(self):
		self.children = []
	def generate(self):
		self.energy = 2*pInvestment
		self.decisionMaker = DecisionTree()
		if pGen_DMStructure:
			self.decisionMaker.generateWithStructure(pGen_DMStructure)
		else:
			self.decisionMaker.generate()
		self.maxAge = gauss(pGen_MaxAge['mean'], pGen_MaxAge['sd'])
	def makeObservations(self, sim):
		obs = [o(self, sim) for o in pObs]
		self.actToDo = self.decisionMaker.decide(obs)
	def run(self, sim):
		global gBasicStats
		self.energy -= pPerCycleEnergy
		for child in self.children:
			child.energy -= pBurden(self, child, sim)
		# Check for death, or act already performed
		if sim.cycle - self.birthCycle > self.maxAge:
			gBasicStats['deaths'] += 1
			return [self]
		if self.energy <= 0:
			gBasicStats['deaths'] += 1
			return [self]
		if sim.cycle == self.lastActCycle:
			return
		self.lastActCycle = sim.cycle
		gBasicStats['actsDone'] += 1
		return pActs[self.actToDo](self, sim)
	def getAge(self, sim):
		return sim.cycle - self.birthCycle
	####################
	# Actions          #
	####################
	# ALL action functions MUST start with 'act'
	def actRest(self, sim):
		pass
	def actWalk(self, sim):
		while 1:
			(x,y) = (randint(-1,1),randint(-1,1))
			if x!=0 or y!=0: break
		sim.board.moveBy(self, x, y)
	def actEat(self, sim):
		food = sim.board.getEntity(self, pEatWidth, pEatHeight, Food)
		if food:
			self.energy += food.energy
			return [food]
	def actSuicide(self, sim):
		if pSuicideOffspringHarm:
			for child in self.children:
				child.energy -= pSuicideOffspringHarm
		gBasicStats['suicides'] += 1
		gBasicStats['deaths'] += 1
		return [self]
	def requestMate(self, requester, sim):
		if pCantReproduce(self, sim):
			return False
		obs = map(lambda o: o(self, sim), pObs)
		actToDo = self.decisionMaker.decide(obs)
		return pActs[actToDo]==Agent.actReproduce
	def actReproduce(self, sim):
		global gBasicStats
		if pCantReproduce(self, sim):
			return
		otherAgent = sim.board.getEntity(self, pMateWidth, pMateHeight, Agent,
			lambda a,sim: a.lastActCycle!=sim.cycle, [sim])
		if otherAgent:
			if not pMateConsensual or otherAgent.requestMate(self, sim):
				if self.energy > pMateRequiredEnergy and otherAgent.energy > pMateRequiredEnergy:
					child = Agent()
					placedSomewhere = False
					if pBirthsAnywhere:
						### Pick some random agent to put the child next to
						### (That way, will have someone to mate with)
						someOtherAgent = choice(sim.getPopulation('agents').members)
						placedSomewhere = sim.board.placeAnywhereWithin(child, someOtherAgent, pBirthWidth, pBirthHeight)
					else:
						placedSomewhere = sim.board.placeAnywhereWithin(child, self, pBirthWidth, pBirthHeight)
					if placedSomewhere:
						### Child is born
						gBasicStats['births'] += 1
						sim.getPopulation("agents").members.append(child)
						### PI child.energy = 
						child.energy = 2*pInvestment
						self.energy -= pInvestment
						otherAgent.energy -= pInvestment
						self.children.append(child)
						otherAgent.children.append(child)
						child.birthCycle = sim.cycle
						child.lastActCycle = sim.cycle
						child.maxAge = gauss(pGen_MaxAge['mean'], pGen_MaxAge['sd'])
						child.decisionMaker = self.decisionMaker.crossover(otherAgent.decisionMaker)
						return
		self.energy -= pMateFailCost
	############################
	# Observations             #
	############################
	# ALL observation functions MUST start with 'obs'
	def obsEnergy(self, sim):
		return float(self.energy) / pMostEnergy
	def obsAge(self, sim):
		#age = sim.cycle - self.birthCycle
		return float(sim.cycle-self.birthCycle) / pMostAge
	def obsGlobalFood(self, sim):
		#print "global food:", float(len(sim.getPopulation("food").members))/(sim.board.width*sim.board.height)
		return float(len(sim.getPopulation("food").members))/(sim.board.width*sim.board.height)
	def obsGlobalPop(self, sim):
		#print "global pop:", float(len(sim.getPopulation("agents").members))/(sim.board.width*sim.board.height)
		return float(len(sim.getPopulation("agents").members))/(sim.board.width*sim.board.height)
	def obsGlobalRelPop(self, sim):
		if sim.maxPopValue == sim.minPopValue:
			return 0
		r = (self.obsGlobalPop(sim)-sim.minPopValue)/(sim.maxPopValue-sim.minPopValue)
		return r
		
Agent.boardLabel = "A"

# Proto is an almost-agent. Not normally used.
class Proto:
	def __init__(self):
		self.pv = [] # Vector specifying what protos can do
	
	def generate(self):
		for i in xrange(pNumActs):
			self.pv.append(random())
			
	def run(self, sim):
		actToDo = self.actions[multinomial(self.pv)] # Pick an action according to prob vector
		actToDo(self, sim) #Run that action
		
	def reproduce(self, sim):
		newProto = Proto()
		# Reproduce pv
		for i, prob in enumerate(self.pv):
			newProto.pv.append(self.pv[i] + gauss(0, pMutateActProb))
		sim.newMembers.append(newProto)
		
	def doNothing(self, sim):
		#print "nothing"
		pass
		





################################################
#  `6. Simulation Running                      #
################################################


#----------------------------------------------#
#  `6-1. Parameters                            #
#----------------------------------------------#

pBoardWidth = 50
pBoardHeight = 50
pMutateActProb = 0.002
pDTFixedStructure = True
pMutate_DTreeSplitValue = 0
pMutate_DTreeAttribute = 0
pMutate_DTreeStructureLeaf = 0
pMutate_DTreeStructureBranch = 0
pGen_DTreeBranch = 0.35
Proto.actions = [Proto.reproduce, Proto.doNothing]
pObsNames = None
pObs = [Agent.obsEnergy, Agent.obsAge, Agent.obsGlobalFood, Agent.obsGlobalPop]
#pActs = [Agent.actRest, Agent.actWalk, Agent.actEat, Agent.actReproduce, Agent.actSuicide]
pActs = [Agent.actWalk, Agent.actEat, Agent.actReproduce, Agent.actSuicide]
#pGen_ActProbs = [0, 0, 0.5, 0.5, 0]
#pGen_ActProbs = [0, 0, 0.7, 0.3, 0]
pGen_ActProbs = [0.25, 0.25, 0.25, 0.25]
pGen_ActProbSd = 0.25
pGen_MaxAge = {'mean':60, 'sd':5}
pGen_NumAgents = 1000
def pGen_AgentsModFunc(en): en.energy = 400
pGen_NumInitialFood = 300
pFoodMaxAge = 8
#pGen_DMStructure = None
#pGen_DMStructure = "1 0.3 ( L 1 0.5 ( L 3 0.5 ( L 0 0.05 ( L L ) ) ) )"
#pGen_DMStructure = "Age 0.35 ( L Energy 0.1 ( L L ) )"
pGen_DMStructure = '''Age 0.05 (
	Energy 0.025 (
		GlobalFood 0.025 ( L L )
		GlobalFood 0.025 ( L L )
	)
	Age 0.2 (
		Energy 0.04 (
			GlobalFood 0.025 ( L L )
			GlobalFood 0.025 ( L L )
		)
		Energy 0.08 (
			GlobalFood 0.025 ( L L )
			GlobalFood 0.025 ( L L )
		)
	)
)'''

pSimTag = 'T012'
pInvestment = 50
#0.104
pMateRequiredEnergy = 200
pMateConsensual = False
pMateFailCost = 3
pEatWidth = 25
pEatHeight = 25
pMateWidth = 25
pMateHeight = 25
pBirthWidth = 5
pBirthHeight = 5
pMostEnergy = 2400
pMostAge = 140
pEnergyGroupSize = 200
pAgeGroupSize = 4
pPerCycleEnergy = 8
pFoodEnergy  = {'mean':130, 'sd':5}
pFoodPeriod = 60
pFoodMagnitude = 60
pFoodOffset = 60
pSuicideOffspringHarm = 45
#Turn kin selection on/off
pBirthsAnywhere = False
pMinimumAgents = 10
def pCantReproduce(self, sim):
	return False #self.obsAge(sim)>=0.5
def pBurden(self, child, sim):
	return self.obsAge(sim)>=0.5 and 0 or 0
def pConstantFood(cycle):
	return 60
def pCyclicalFood(cycle):
	return int(pFoodOffset + pFoodMagnitude*math.sin(float(cycle)*(2*math.pi)/pFoodPeriod))
def pPeriodicDrought(cycle):
	if cycle/8 % 6 == 5: return 0
	return 70
pFoodPerCycle = pCyclicalFood

## Constants once parameters are defined
pNumActs = len(pActs)
pNumObs = len(pObs)
gObsPosToName = [o.__name__.replace('obs', '') for o in pObs]
gObsNameToPos = {}
for i,obsName in enumerate(gObsPosToName):
	gObsNameToPos[obsName] = i

fieldNames = ['population', 'energy', 'age', 'births', 'deaths', 'suicideProp']
def handleFieldPrint(value):
	if isinstance(value, ProgressiveStats):
		return str(value.calcMean())
	return str(value)

gBasicStats = {
	'births': 0,
	'deaths': 0,
	'suicides': 0,
	'actsDone': 0,
}

# This extends the Simulation class to deal with the current experiment.
# I am thinking I should just merge this back into the first definition.
class Simulation(Simulation):
	# Called at the end of each epoch
	def stats(self, statsFile, energyFile, ageFile):
		# Get the average probability vector
		avgPv = [0 for i in xrange(pNumActs)]
		agentPop = self.getPopulation('agents')
		fields = {
			'energy': ProgressiveStats(),
			'age': ProgressiveStats(),
			'births': gBasicStats['births'],
			'deaths': gBasicStats['deaths'],
			'suicideProp': gBasicStats['suicides']/(gBasicStats['actsDone']+0.00001),
		}
		# Clear the basic stats for the new epoch
		for fn in gBasicStats.iterkeys():
			gBasicStats[fn] = 0
			
		energyDistribution = [0 for i in xrange(0,pMostEnergy,pEnergyGroupSize)]
		ageDistribution = [0 for i in xrange(0,pMostAge,pAgeGroupSize)]
		for agent in agentPop.members:
			fields['energy'].update(agent.energy)
			fields['age'].update(agent.obsAge(self))
			energyGroup = int(agent.energy/pEnergyGroupSize)
			energyGroup = (agent.energy<pMostEnergy and [energyGroup] or [int((pMostEnergy-1)/pEnergyGroupSize)])[0]
			energyDistribution[energyGroup] += 1
			ageGroup = int(agent.getAge(self)/pAgeGroupSize)
			ageGroup = (agent.getAge(self)<pMostAge and [ageGroup] or [int((pMostAge-1)/pAgeGroupSize)])[0]
			ageDistribution[ageGroup] += 1
		
		energyFile.write(str(self.cycle)+'\t')
		for energyCount in energyDistribution:
			energyFile.write(str(energyCount))
			energyFile.write('\t')
		energyFile.write('\n')
		ageFile.write(str(self.cycle)+'\t')
		for ageCount in ageDistribution:
			ageFile.write(str(ageCount))
			ageFile.write('\t')
		ageFile.write('\n')

		statsFile.write(str(self.cycle)+'\t')
		fields['population'] = len(self.populations[0].members)
		for field in fieldNames:
			statsFile.write(handleFieldPrint(fields[field]))
			statsFile.write('\t')
		statsFile.write('\n')
	def generateNewPopulation(self, popName, pyClass, num, modFunc = None):
		self.addPopulation(popName, pyClass)
		self.generatePopulation(popName, num, modFunc)
	def generatePopulation(self, popName, num, modFunc = None):
		pop = self.getPopulation(popName)
		newMembers = pop.generate(num)
		for en in newMembers:
			if modFunc: modFunc(en)
			if not self.board.placeAnywhere(en):
				pop.members.remove(en)
			else:
				en.birthCycle = self.cycle
	def addToPopulation(self, popName, popMembers):
		pop = self.getPopulation(popName)
		for en in popMembers:
			if self.board.placeAnywhere(en):
				en.birthCycle = self.cycle
				pop.members.append(en)
	def generate(self, numAgents = None):
		self.generateNewPopulation('agents', Agent, numAgents==None and pGen_NumAgents or numAgents, pGen_AgentsModFunc)
		self.generateNewPopulation('food', Food, pGen_NumInitialFood)

def mout(msg):
	print msg


#----------------------------------------------#
#  `6-2. Experiment management                 #
#----------------------------------------------#

# Configuration for running the current experiment
# EXAMPLE:
#   exp = Experiment()
#   exp.run()
#   -> runs all the simulations for this experiment
#   exp.runSimulation(['me', 'main/p0/', runId], 5000)
#	-> Runs a single simulation for 5000 cycles, puts the results into the
#      the directory 'main/p0/', and prefixes the output files with the
#      value (normally an int) in runId. ('me' is just to substitute for
#      argv[0], which is normally the name of the script.)
class Experiment:
	def getOutFileName(self, label, runSet = None, runId = None, ext = ".txt"):
		if runSet == None: runSet = self.runSet
		if runId == None: runId = self.runId
		
		# runId might not be a number
		try:
			return "%s%03d-%s.txt" % (runSet, int(runId), label)
		except ValueError:
			return "%s%s-%s.txt" % (runSet, runId, label)
		
	def runSimulation(self, argv = sys.argv, cycles = 1000, agents = None):
		runSet = argv[1]
		runId = str(argv[2])
		self.runSet = runSet
		self.runId = runId

		try:
			os.makedirs(runSet)
		except:
			pass

		s = Simulation()
		if agents:
			s.generate(0)
			s.addToPopulation('agents', agents)
		else:
			s.generate()
		s.maxPopValue = 0
		s.minPopValue = 0
		s.globalPopWindow = []
		s.board.display()
		c = time.clock()
		t = time.time()
		ageFile = open(self.getOutFileName("age"), "w")
		energyFile = open(self.getOutFileName("energy"), "w")
		probFile = open(self.getOutFileName("demog"), "w")
		genomeFile = open(self.getOutFileName("genome"), "w")
		energyFile.write('Cycle\t')
		for energyLevel in xrange(0, pMostEnergy, 200):
			energyFile.write(str(energyLevel)+"\t")
		energyFile.write("\n")
		ageFile.write('Cycle\t')
		for ageLevel in xrange(0, pMostAge, pAgeGroupSize):
			ageFile.write(str(ageLevel)+"\t")
		ageFile.write("\n")
		probFile.write('Cycle\t')
		probFile.write('\t'.join(fieldNames))
		probFile.write('\n')
		for i in xrange(cycles):
			### Duplicate existing agents if there are too few
			### (Duplicates everything, including state (age/energy))
			if len(s.getPopulation('agents').members)<pMinimumAgents:
				mout('!!!!!!!!!!!!!!!!!!!!!!!!! --  Population < '+str(pMinimumAgents))
				errorFile = open(self.getOutFileName("errors"), "a")
				print >>errorFile, "Population <", pMinimumAgents, "in Cycle", s.cycle, "; generating more agents"
				errorFile.close()
				newAgents = []
				for i in xrange(5):
					newAgents.append(deepcopy(choice(s.getPopulation('agents').members)))
					newAgents[-1].energy = pInvestment*2
					newAgents[-1].age = 0
				s.addToPopulation('agents', newAgents)
				# s.generatePopulation('agents', 5)
			if i % 20==0:
				print "Cycle:",s.cycle
				print ""
				s.board.display()
				totalEnergy = 0
				co = 0 # count of agents
				avgDTree = None
				totalMembers = len(s.getPopulation('agents').members)
				for ag in s.getPopulation('agents').members:
					#print "chil: %d" % (len(ag.children))
					totalEnergy += ag.energy
					co += 1
					### Decision tree
					if pDTFixedStructure:
						if not avgDTree:
							avgDTree = deepcopy(ag.decisionMaker)
							avgDTree.normalize()
							avgRawDTree = deepcopy(ag.decisionMaker)
						else:
							avgNode = avgDTree.rootNode
							avgRawNode = avgRawDTree.rootNode
							agentNode = ag.decisionMaker.rootNode
							avgDTree.dtreeCalled += ag.decisionMaker.dtreeCalled
							avgRawDTree.dtreeCalled += ag.decisionMaker.dtreeCalled
							avgToCheck = [avgNode]
							avgRawToCheck = [avgRawNode]
							agentToCheck = [agentNode]
							while len(avgToCheck):
								if isinstance(avgNode, DecisionBranch):
									avgToCheck.extend(avgNode.childNodes)
									avgRawToCheck.extend(avgRawNode.childNodes)
									agentToCheck.extend(agentNode.childNodes)
								else: ### avgNode is leaf
									avgNode.called += agentNode.called
									avgRawNode.called += agentNode.called
									normedAgentPv = normalizeMultinomial(agentNode.pv)
									for i,prob in enumerate(avgNode.pv):
										avgNode.pv[i] += normedAgentPv[i]
										avgRawNode.pv[i] += agentNode.pv[i]
										if co==totalMembers:
											avgNode.pv[i] /= co
											avgRawNode.pv[i] /= co
								avgNode = avgToCheck.pop()
								avgRawNode = avgRawToCheck.pop()
								agentNode = agentToCheck.pop()
				if pDTFixedStructure:
					print "AvgDTree:"
					if avgDTree: avgDTree.display()

				print len(s.getPopulation('agents').members)
				sag = choice(s.getPopulation('agents').members)
				sag.decisionMaker.display()
				print "Avg agent energy: %f" % (float(totalEnergy)/co)
				print "Num agents: %d" % (co)
				
				print >>genomeFile, "{'Cycle':", s.cycle, ","
				if avgDTree:
					print >>genomeFile, "'Average Decision Tree': "
					avgDTree.display(genomeFile)
					print >>genomeFile
					print >>genomeFile, "'[UNNORMALISED] Average Decision Tree': "
					avgRawDTree.display(genomeFile, showUnNormed = True)
					print >>genomeFile
				
				print >>genomeFile, "'Random Decision Tree': "
				sag.decisionMaker.display(genomeFile)
				print >>genomeFile
				print >>genomeFile, "'[UNNORMALISED] Random Decision Tree': "
				sag.decisionMaker.display(genomeFile, showUnNormed = True)
				print >>genomeFile
				print >>genomeFile

			s.generatePopulation('food', pFoodPerCycle(s.cycle))
			s.globalPopWindow.append(Agent.obsGlobalPop(Agent(), s))
			if (len(s.globalPopWindow))>120:
				s.globalPopWindow.pop(0)
			maxPopValue = -1
			minPopValue = 2
			for popValue in s.globalPopWindow:
				if popValue > maxPopValue:
					maxPopValue = popValue
				if popValue < minPopValue:
					minPopValue = popValue
			s.maxPopValue = maxPopValue
			s.minPopValue = minPopValue
			#print "Pop: (%f,%f)" % (s.minPopValue, s.maxPopValue)
			#print "PopDens:", Agent.obsGlobalPop(Agent(), s)
			#print "       PopRel:", Agent.obsGlobalRelPop(Agent(), s)
			s.run(1)
			if len(s.getPopulation('agents').members)==0:
				break
			s.stats(probFile, energyFile, ageFile)
		print time.clock() - c
		print time.time() - t
		return s

	def run(self):
		for runId in xrange(1, 41):
			self.runSimulation(['me', 'main/'+pSimTag+'/', runId], 5000)
			#pBirthsAnywhere = True
			#self.runSimulation(['me', 'main/p1/', runId], 20000)
	
	def blendFiles(self, outputFileName, inputFileNamePatterns):
		if type(inputFileNamePatterns) is not list: inputFileNamePatterns = [inputFileNamePatterns]
		
		outputFile = open(outputFileName, 'w')
		inputFileNames = []
		for inputFileNamePattern in inputFileNamePatterns:
			inputFileNames.extend(glob.glob(inputFileNamePattern))
		
		print "outputFile:", outputFileName
		print "inputFiles:", inputFileNames
		
		inputFiles = [open(inputFileName) for inputFileName in inputFileNames]
		
		numCols = 0
		for f in inputFiles: ### Headers
			header = f.readline()
			if numCols==0:
				numCols = len(header.split('\t'))
		outputFile.write(header)

		moreInput = True
		while moreInput:
			moreInput = False
			avgLine = [0 for i in xrange(numCols)]
			numInputFiles = 0
			for f in inputFiles:
				ln = f.readline()
				if ln:
					moreInput = True
				else:
					break
				avgLine = [float(val1)+val2 for val1, val2 in zip(ln[:-1].split('\t'),avgLine)]
				numInputFiles += 1
			if numInputFiles>0:
				avgLine = [val / numInputFiles for val in avgLine]
				outputFile.write('\t'.join([str(val) for val in avgLine]) + '\n')

	def gatherStats(self):
		self.blendFiles('main/'+pSimTag+'-demog_summary.txt', 'main/'+pSimTag+'/*demog*')

'''
#Debugging
dt1 = DecisionTree()
dt1.generate()
dt1.display()
dt2 = DecisionTree()
dt2.generate()
dt2.display()
child = dt1.crossover(dt2)
child.display()
print child.decide([0.4,0.2,0.3,0.7])
'''

# This is where the program starts running:
if __name__ == '__main__':
	exp = Experiment()
	exp.run()
	exp.gatherStats()