### Alex's Genetic Programming Algorithm
### ------------------------------------

### For documentation see Genetic Programming III (Koza et al 1999)

### Copyright (C) 2001 Alex Wilding
### If you modify this code, either for your own purposes or otherwise,
### PLEASE PLEASE PLEASE PLEASE PLEASE PLEASE let me know at least what
### you have done, or better, the source.

### This source is licensed under the GNU General Public License.
### The GPL can be accessed at: http://www.gnu.org/licenses/gpl.html

### This program is free software; you can redistribute it and/or
### modify it under the terms of the GNU General Public License
### as published by the Free Software Foundation; either version 2
### of the License, or any later version.

### This program is distributed in the hope that it will be useful,
### but WITHOUT ANY WARRANTY; without even the implied warranty of
### MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
### GNU General Public License for more details.

### You should have received a copy of the GNU General Public License
### along with this program; if not, write to the Free Software
### Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

import math
import random, whrandom
import string
from copy import deepcopy
import sys

####################### RUN PARAMETERS

popcount = 2000
matesubset = 400
genmax = -1
initialdepth = 2



def GetFitness(Ind, verbosefile=None):      ## MODIFY ME FOR DIFFERENT PROBLEMS
  AbsoluteMaxNodes = 35
  parsimony = float(float(max(0,AbsoluteMaxNodes-len(GetIndividualTree(Ind))))/AbsoluteMaxNodes)   # between 0 and 1

  #print parsimony
  # target remains stationary, while mx moves in direction determined by result. If result is positive, move right. If negative, move left.
#  for i in range(100): # run the simulation for 100 cycles.
#    a = ev(Ind.Events[0]);
#    # move a maximum of 10 units.
#    v = v + max(-2,min(2,v))
#    v = v * 0.9;
#    x = x + v;

  # Initialisation
  x = (2*random.random()-1)*10
  y = (2*random.random()-1)*10
  xv = 0
  yv = 0
  
  for count in range(50):
    Ind.Parameters()[0] = x;
    Ind.Parameters()[1] = y;
    Ind.Parameters()[2] = xv;
    Ind.Parameters()[3] = yv;
    if verbosefile != None:
      verbosefile.write(str(x)+","+str(y)+","+str(xv)+","+str(yv)+",");
    angle = ev(Ind.Events[0]);
    xv = xv + math.cos(angle)*0.4
    yv = yv + math.sin(angle)*0.4
    xv = xv * 0.9;
    yv = yv * 0.9;
    x += xv
    y += yv
  if verbosefile != None:
    verbosefile.write("\n");  
  result = math.sqrt(x*x+y*y);
  ideal = 0;
  maximum = 10
  parsimonyweight = 15
  fitness =max(maximum-abs(result-ideal),0);
  
  # regulate to non-negative only.
  fitness = float(max(0,fitness));  # fitness is 100-end distance from target.
  return [(100-parsimonyweight)*(fitness/maximum)+parsimony*parsimonyweight,fitness,parsimony]




pass #################################### ** Init:
TNone = -1
TInteger = 0
TBoolean = 1
TFloat = 2
pop = []
newpop = []
pass #################################### ** ABSTRACT NODES
class Node:
  def __init__(self, fIndividual):
    self.OutType = TNone;
    self.Parent = None;
    self.AcceptsTypes = []
    self.Children = []
    self.Individual = fIndividual;
    self.Name = "Abstract"
    #print "Created"
  def __str__(self):
    return str(ev(self))
  def __eval__(self):
    return "Abstract Base Class"
  def Accepts(self,ftype):
    for q in self.AcceptsTypes:
      if q == ftype:
        return 1
    return 0
  def __tree__(self, indent):
    # printing of node tree
    r =  indent+self.Name+": "+str(self)+"\n";
    for child in self.Children:
      r = r + child.__tree__(indent+"  ");
    return r;
  def __random__(self, depth=1):
    return None;
#

class LinkNode(Node):
  def __init__(self, fIndividual):
    Node.__init__(self, fIndividual);
    self.Children = []
    self.Name = "Abstract Link"
  def __eval__(self):
    return "Abstract Link"
  def GetRandomChildren(self,depth = 1, NumberOfChildren = 1):
    self.Children = [];
    for i in range(NumberOfChildren):
      self.Children = self.Children + [FindMeAChild(self.AcceptsTypes,(depth==0))(self.Individual)]
    for child in self.Children:
      randomize(child,depth-1)

class UnaryNode(LinkNode):
  def __init__(self, fIndividual, fChild = None):
    LinkNode.__init__(self, fIndividual);
    self.Children = [fChild]
    self.Name = "Abstract Unary"
  def __eval__(self):
    return "Abstract Unary"
  def __random__(self, depth = 1):
    self.GetRandomChildren(depth, 1);

class BinaryNode(LinkNode):
  def __init__(self, fIndividual, fChildA = None, fChildB = None):
    LinkNode.__init__(self, fIndividual);
    self.Children = [fChildA, fChildB]
    self.Name = "Abstract Binary"
  def __eval__(self):
    return "Abstract Binary"
  def __random__(self, depth = 1):
    self.GetRandomChildren(depth, 2);

class ThreewayNode(LinkNode):
  def __init__(self, fIndividual, fChildA = None, fChildB = None, fChildC = None):
    LinkNode.__init__(self, fIndividual);
    self.Children = [fChildA, fChildB, fChildC]
    self.Name = "Abstract Threeway"
  def __eval__(self):
    return "Abstract Threeway"
  def __random__(self, depth = 1):
    self.GetRandomChildren(depth, 3);

class ComparativeBinaryNode(UnaryNode):
  def __init__(self, fIndividual, fChildA = None, fChildB = None):
    BinaryNode.__init__(self, fIndividual,fChildA, fChildB);
    self.OutType = TBoolean;
    self.AcceptsTypes = [TFloat, TBoolean, TInteger]
    self.Name = "Abstract Comparative Binary Node"
#    if isinstance(self.Children[0],Node):
#      if not self.Accepts(self.Children[0].OutType):
#        print "ERROR - Binary comparative attempted with wrong type"
#        # raise an error or something?

class LiteralNode(Node):
  def __init__(self, fIndividual):
    Node.__init__(self, fIndividual);
    self.Name = "Abstract Literal"
  def __eval__(self):
    return "Abstract Literal"


pass #################################### ** SPECIFIC NODES
pass ################ (Literals)
class IntegerLiteralNode(LiteralNode):
  def __init__(self, fIndividual):
    LiteralNode.__init__(self, fIndividual)
    self.val = 0;
    self.OutType = TInteger
    self.Name = "Integer Literal"
  def __eval__(self):
    return self.val;
  def __random__(self, depth=1):
    self.val = int(random.gauss(0,10.0)*random.gauss(0,10.0))

class BooleanLiteralNode(LiteralNode):
  def __init__(self, fIndividual):
    LiteralNode.__init__(self, fIndividual)
    self.val = 0;
    self.OutType = TBoolean;
    self.Name = "Boolean Literal"
  def __str__(self):
    if self.val == 0:
      return "False"
    else:
      return "True";
  def __eval__(self):
    return self.val;
  def __random__(self, depth=1):
    self.val = random.choice(range(2)); # chooses 0 or 1 randomly.

class FloatLiteralNode(LiteralNode):
  def __init__(self, fIndividual):
    LiteralNode.__init__(self, fIndividual)
    self.val = float(0);
    self.OutType = TFloat;
    self.Name = "Float Literal"
  def __eval__(self):
    return self.val;
  def __random__(self, depth = 1):
    self.val = (random.gauss(0,1.0))

class ParameterLiteralNode(LiteralNode):
  def __init__(self, fIndividual):
    LiteralNode.__init__(self, fIndividual)
    self.OutType = TFloat;
    self.Name = "Parameter Literal"
    self.ParamID = -1;
  def __str__(self):
    return "P["+str(self.ParamID)+"]="+str(ev(self))
  def __eval__(self):
    if self.ParamID > -1:
      return self.Individual.Parameters()[self.ParamID];
    else:
      return 0
  def __random__(self, depth = 1):
#    try:
      self.ParamID = random.choice(range(len(self.Individual.Parameters()))) # select a random parameter to reference.
#    except:
#      self.ParamID = -1;

Literals = [IntegerLiteralNode, BooleanLiteralNode, FloatLiteralNode, ParameterLiteralNode]
pass ################ (Unaries)
class SineUnaryNode(UnaryNode):
  def __init__(self, fIndividual, fChild = None):
    UnaryNode.__init__(self, fIndividual,fChild);
    self.OutType = TFloat;
    self.AcceptsTypes = [TFloat, TBoolean, TInteger]
    self.Name = "Sine Unary"
    if isinstance(self.Children[0],Node):
      if not self.Accepts(self.Children[0].OutType):
        print "ERROR - Sine Unary Create attempted with wrong type"
        # raise an error or something?
  def __str__(self):
    return str(ev(self))
  def __eval__(self):
    return math.sin(ev(self.Children[0]));

class CosineUnaryNode(UnaryNode):
  def __init__(self, fIndividual, fChild = None):
    UnaryNode.__init__(self, fIndividual,fChild);
    self.OutType = TFloat;
    self.AcceptsTypes = [TFloat, TBoolean, TInteger]
    self.Name = "Cosine Unary"
  def __str__(self):
    return str(ev(self))
  def __eval__(self):
    return math.cos(ev(self.Children[0]));

class AbsUnaryNode(UnaryNode):
  def __init__(self, fIndividual, fChild = None):
    UnaryNode.__init__(self, fIndividual,fChild);
    self.OutType = TFloat;
    self.AcceptsTypes = [TFloat, TInteger]
    self.Name = "Abs Unary"
  def __str__(self):
    return str(ev(self))
  def __eval__(self):
    return abs(ev(self.Children[0]));

class NegateUnaryNode(UnaryNode):
  def __init__(self, fIndividual, fChild = None):
    UnaryNode.__init__(self, fIndividual,fChild);
    self.OutType = TFloat;
    self.AcceptsTypes = [TFloat, TInteger, TBoolean]
    self.Name = "Negate Unary"
  def __str__(self):
    return str(ev(self))
  def __eval__(self):
    return -(ev(self.Children[0]));

class SquareRootUnaryNode(UnaryNode):
  def __init__(self, fIndividual, fChild = None):
    UnaryNode.__init__(self, fIndividual,fChild);
    self.OutType = TFloat;
    self.AcceptsTypes = [TFloat, TInteger]
    self.Name = "Square Root Unary"
  def __str__(self):
    return str(ev(self))
  def __eval__(self):
    return math.sqrt(abs(ev(self.Children[0])));

class SquareUnaryNode(UnaryNode):
  def __init__(self, fIndividual, fChild = None):
    UnaryNode.__init__(self, fIndividual,fChild);
    self.OutType = TFloat;
    self.AcceptsTypes = [TFloat, TInteger]
    self.Name = "Square Unary"
  def __str__(self):
    return str(ev(self))
  def __eval__(self):
    return (ev(self.Children[0]))**2;

class PassThroughUnaryNode(UnaryNode):
  def __init__(self, fIndividual, fChild = None):
    UnaryNode.__init__(self, fIndividual,fChild);
    self.OutType = TFloat;
    self.AcceptsTypes = [TFloat, TInteger]
    self.Name = "PassThrough Unary"
  def __str__(self):
    return str(ev(self))
  def __eval__(self):
    return ev(self.Children[0]);

Unaries = [SineUnaryNode, CosineUnaryNode, AbsUnaryNode, NegateUnaryNode, SquareRootUnaryNode, SquareUnaryNode];
pass ################ (Binaries)
class AddBinaryNode(BinaryNode):
  def __init__(self, fIndividual, fChildA = None, fChildB = None):
    BinaryNode.__init__(self, fIndividual)
    self.AcceptsTypes = [TFloat, TInteger, TBoolean]
    self.OutType = TFloat;
    self.Name = "Add Binary"
  def __eval__(self):
    return ev(self.Children[0])+ev(self.Children[1]);

class ArcTan2BinaryNode(BinaryNode):
  def __init__(self, fIndividual, fChildA = None, fChildB = None):
    BinaryNode.__init__(self, fIndividual)
    self.AcceptsTypes = [TFloat, TInteger, TBoolean]
    self.OutType = TFloat;
    self.Name = "ArcTan2 Binary"
  def __eval__(self):
    return math.atan2(ev(self.Children[1]),ev(self.Children[1]));

class SubtractBinaryNode(BinaryNode):
  def __init__(self, fIndividual, fChildA = None, fChildB = None):
    BinaryNode.__init__(self, fIndividual)
    self.AcceptsTypes = [TFloat, TInteger, TBoolean]
    self.OutType = TFloat;
    self.Name = "Subtract Binary"
  def __eval__(self):
    return ev(self.Children[0])-ev(self.Children[1]);

class MultiplyBinaryNode(BinaryNode):
  def __init__(self, fIndividual, fChildA = None, fChildB = None):
    BinaryNode.__init__(self, fIndividual)
    self.AcceptsTypes = [TFloat, TInteger, TBoolean]
    self.OutType = TFloat;
    self.Name = "Multiply Binary"
  def __eval__(self):
    return ev(self.Children[0])*ev(self.Children[1]);

class DivideBinaryNode(BinaryNode):
  def __init__(self, fIndividual, fChildA = None, fChildB = None):
    BinaryNode.__init__(self, fIndividual)
    self.AcceptsTypes = [TFloat, TInteger]
    self.OutType = TFloat;
    self.Name = "Divide Binary"
  def __eval__(self):
    try:
      return ev(self.Children[0])/ev(self.Children[1]);
    except:
      return 0

class LessEqualBinaryNode(ComparativeBinaryNode):
  def __init__(self, fIndividual, fChildA = None, fChildB = None):
    ComparativeBinaryNode.__init__(self, fIndividual,fChildA, fChildB);
    self.Name = "Lesser/Equal Binary Node"
  def __eval__(self):
    return (ev(self.Children[0]) <= ev(self.Children[1]))

class GreaterEqualBinaryNode(ComparativeBinaryNode):
  def __init__(self, fIndividual, fChildA = None, fChildB = None):
    ComparativeBinaryNode.__init__(self, fIndividual,fChildA, fChildB);
    self.Name = "Lesser/Equal Binary Node"
  def __eval__(self):
    return (ev(self.Children[0]) >= ev(self.Children[1]))

class PowerBinaryNode(BinaryNode):
  def __init__(self, fIndividual, fChildA = None, fChildB = None):
    BinaryNode.__init__(self, fIndividual,fChildA, fChildB);
    self.Accepts = [TFloat, TInteger]
    self.Name = "Raise To Power Binary Node"
  def __eval__(self):
    try:
      return math.pow(abs(ev(self.Children[0])),ev(self.Children[1]))
    except:
      return 0

Binaries = [Arctan2BinaryNode, AddBinaryNode, SubtractBinaryNode, MultiplyBinaryNode, DivideBinaryNode, PowerBinaryNode]
pass ################ (Threeways)
class IsNearThreewayNode(ThreewayNode): # Target, Try, Radius
  def __init__(self, fIndividual, fChildA = None, fChildB = None, fChildC = None):
    ThreewayNode.__init__(self, fIndividual,fChildA, fChildB, fChildC);
    self.AcceptsTypes = [TFloat, TInteger, TBoolean]
    self.OutType = TBoolean
    self.Name = "Is Near Threeway Node"
  def __eval__(self):
    return (abs(ev(self.Children[0])-ev(self.Children[1])) <= abs(ev(self.Children[2])))

Threeways = [IsNearThreewayNode]
pass #################################### ** Helper Functions
Nodes = Literals + Unaries + Binaries #+ Threeways;
NodeTypes = {}
# Make a dictionary of every result type.
for noddy in Nodes:
  NodeTypes.update({ noddy :noddy(None).OutType})

def ev(N):
  return N.__eval__()

def randomize(N, depth = 1): # randomizes an instance rather than creating an instance
  return N.__random__(depth)

def FindMeAChild(Types,Terminal): # terminal says we must not return a linking node. Returns a *type*
  ShortList = [];
  TypeList = list(Types)
  if Terminal:
    LongList = Literals;
  else:
    LongList = Nodes;
  for n in LongList:
    if TypeList.count(NodeTypes[n]) > 0:
      ShortList.append(n);
  return whrandom.choice(ShortList)

# Generate complete list of node types available.
pass #################################### ** GA Classes

stepsize = (popcount / 70)
if stepsize == 0:
  stepsize = 1;

class Individual:
  def __init__(self):
    self.Events = [PassThroughUnaryNode(self)];
    self.Fitness = [0,0,0,0]
    self.MyParams = []
  def __random__(self,depth = 11):
    for R in self.Events:
      randomize(R,depth)
  def Parameters(self):
    return self.MyParams

def GetIndividualTree(Ind):
  res = []
  for e in Ind.Events:
    res = res + GetTree(e)
  return res

def GetTree(Root): # return all nodes under a tree as a list. (for randomly picking a node)
  res = Root.Children
  for C in Root.Children:
    res = res + GetTree(C)
  return res

def Mate(ParentA, ParentB): # sexually mate the two parents, and return a child individual
  # Mate two parents to produce an offspring.
  # Deep copies ParentA to Child, then finds an appropriate match of a subtree between Child and ParentB,
  # and moves the subtree from ParentB to Child. If no match can be found, None is returned as the individuals
  # are incompatible.
  
  # CHANGE ME:
  # Alter this routine so that the pieces swapped are in roughly the same place in the tree.
  # A node is chosen randomly, and the route through the tree (which children are selected)
  # is calculated. This route is applied to the other tree, so a node of equivalent position
  # is swapped, instead of something that will clearly destroy genetic information.
  # Q: What will odd positioning lead to - ie. if a node can't even nearly reach a position?
  Child = Individual()
  Child.Events = deepcopy(ParentA.Events)
  Rind = random.randint(0,len(Child.Events)-1)
  SourceNode = random.choice(GetTree(ParentB.Events[Rind]))
  TargetNode = random.choice(GetTree(Child.Events[Rind]))
  TargetNode = deepcopy(SourceNode);
  for e in Child.Events:
    for n in GetTree(e):
      n.Individual = Child
  return Child

def Mutate(Ind): # mutate a subtree of this individual
  # select a subtree and randomize it. This could easily kill the organism, but is needed to stop the gene pool
  # getting stagnant.
  MuteNode = random.choice(GetTree(random.choice(Ind.Events)))
  randomize(MuteNode,initialdepth)

# All objects are in place. GAtastic!

def ChooseRandom(Population): # choose a random individual based on fitness
  # Each member of Population has a non-negative fitness. Find the total.
  # Get a random number 0..total. Walk through the population and stop
  # when the accumulated fitness > the random number.
##  total = 0
##  for i in pop:
##    total = total + i.Fitness[0]
##  slicepos = random.random()*total
##  tf = 0
##  index = 0
##  #print "Slicepos = "+str(slicepos)
##  while tf < slicepos:
##    tf = tf + pop[index].Fitness[0];
##    index = index + 1;
##  index = index - 1
##  #print "Chosen "+str(index)+" for parent - fitness" + str(pop[index].Fitness);
##  return pop[index]
  # operate on a small subset of nodes.
  c = []
  for i in range(matesubset):
    c = c + [random.randint(0,len(pop)-1)]
  total = 0
  for i in c:
    total = total + pop[i].Fitness[0]
  slicepos = random.random()*total
  tf = 0
  index = 0
  #print "Slicepos = "+str(slicepos)
  while tf < slicepos:
    tf = tf + pop[c[index]].Fitness[0];
    index = index + 1;
  index = index - 1
  #print "Chosen "+str(index)+" for parent - fitness" + str(pop[index].Fitness);
  return pop[c[index]]


def pdot(n): # print a little dot to indicate we've done something. Progress indicator
  if n % stepsize == 0:
    sys.stdout.write(".");


def NewPopulation(depth): # randomly initialise a new population up to popcount.
  global pop#, GlobalParameters
  while len(pop) < popcount:
    pdot(len(pop))
    I = Individual()
    I.MyParams = RandomParams()
    for e in I.Events:
      randomize(e,depth)
    pop = pop + [I]

def MateAll(): # perform mating on the entire population, kill the parents and replace them with children
  global pop, newpop
  newpop = [];
  #print "Mating:"
  while len(newpop) < popcount: # produce new population by sexual crossover.
 #   try:
      newpop = newpop + [Mate(ChooseRandom(pop),ChooseRandom(pop))]
      pdot(len(newpop))
 #   except:
      #doesn't matter - just try again!
 #     pass
  #print "Done Mating"
  pop = newpop;
  newpop = [];
  print ""
  # now we have a new population in wpop.

def GetAllFitness(): #
  j = 0
  for p in pop:
    try:
      j = j + 1
      pdot(j);
      totfit = 0
      for i in range(10):
        p.MyParams = RandomParams()  # randomize parameters for each test
        p.Fitness = GetFitness(p)
        totfit = totfit + p.Fitness[0]
      p.Fitness[0] = totfit / 10.0;
    except: # Stop it from breeding if it excepts. Cruel :)
      print "Error in evaluation"
      p.Fitness = [0,0,0]

# create the population
#GlobalParameters = [0];
def RandomParams():
  return [0,0,0,0]

print "Generating initial population"
NewPopulation(initialdepth)

def ClearFile(filename):
  fitlog = open(filename,"w");
  fitlog.close();

def WriteToFile(filename, line):
  fitlog = open(filename,"a");
  fitlog.write(line)
  fitlog.close();

ClearFile("fitness.csv")
ClearFile("samplelog.csv")
ClearFile("Programs.txt")

gen = 0
while (gen < genmax) or (genmax == -1):
  # make the global parameters
  #GlobalParameters = [random.random()*3];
  #print GlobalParameters
  print "\nMutating:"
  j = 0
  for p in pop: # mutate a few and randomize the old parameters
    p.MyParams = RandomParams()
    j= j+1
    pdot(j)
    try:
      if random.random() < 0.15: # but not too many
        Mutate(p);
    except:
      pass # who cares?!
  # get the fitness ready for reproduction 
  print "\nFinding Fitness:"
  GetAllFitness()
  qmax = pop[0]
  qmin = pop[0]
  s = 0
  ens = 0
  tot = 0
  print "\nBucketing:"
  bucketcount = 20;
  buckets = [];
  for i in range(bucketcount): buckets.append(0);
  for p in pop:
    buckets[int(bucketcount*(p.Fitness[0]+1)/100.0)] = buckets[int(bucketcount*(p.Fitness[0]+1)/100.0)]+1
  s = ""
  for q in buckets:
    s = s + str(q)+","
  s = s + "\n"
  WriteToFile("fitness.csv",s);
  s = 0
  print "\nGetting Best of Generation:"
  for p in pop: # Get the Best Of Generation
    #p.RootNode.__tree__("  ")
    s = s + 1
    try:
      tot = tot + p.Fitness[0]
      if p.Fitness[0] < qmin.Fitness[0]:
        qmin = p
      if p.Fitness[0] > qmax.Fitness[0]:
        qmax = p
        ens = s
    except:
      print "¶ Error on "+str(s)
      print "p.Fitness = "+ str(p.Fitness)
      print "---------end error report"
  s = "\nNode "+str(ens)+":"+str(qmax.Fitness)
  print "\n"

  besttree = qmax.Events[0].__tree__("")
  print besttree

  WriteToFile("Programs.txt","\n\n"+s+"\n\n"+besttree);
  
  runlog = open("samplelog.csv","a")
  GetFitness(qmax,runlog);
  runlog.close();
  
  print s
  # do the reproduction
  print "\nMating for generation "+str(gen)+":"
  MateAll()
  gen = gen + 1

#Dummy = Individual()
#Dummy.Parameters = [2.1,11.44,-2];gpa.py
#randomize(Dummy.RootNode,14)
#print "---"
#Dummy.RootNode.__tree__("  ")
#Mutate(Dummy)
#print "---"
#Dummy.RootNode.__tree__("  ")

#q = FloatLiteralNode(Dummy);
#randomize(q,30)
#p = Unaries[0](Dummy,q)
##print ev(p)
#randomize(p,11)
#print ev(p)
#p.__tree__("|")
#print FindMeAChild(SineUnaryNode(Dummy).AcceptsTypes,1)