Genetic Programming.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Genetic Programming"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "This chapter introduces some very basic concepts about Genetic Programming. As usual, we will require a couple of imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import random\n",
    "import sys\n",
    "import re\n",
    "import matplotlib.pyplot as plt\n",
    "from statistics import mean\n",
    "\n",
    "# For presenting as slides\n",
    "#plt.rcParams['figure.figsize'] = [12, 8]\n",
    "#plt.rcParams.update({'font.size': 22})\n",
    "#plt.rcParams['lines.linewidth'] = 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "## Classic Genetic Programming"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The main difference between regular Genetic Algorithms and Genetic Programming lies in the representation: GAs tend to operate on the genotype, while GPs operate on the phenotype (programs); for this, the programs are traditionally represented as trees. Furthermore, GPs have a more open-ended nature in that the number of elements used in a solution as well as their interconnections must be open to evolution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In a tree representation, programs are typically generated using two sets: A set of terminals (items of arity 0) and functions (items of arity > 0). The input to a function must be the result of any other function that can be defined, which leads to the important property of _type closure_: Each function must be able to handle all values it might ever receive as input, all terminals must be allowable inputs for all functions, and the output from any function must be a permitted input to any other function. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The first step to evolving a suitable program for a given problem is thus to define the sets of functions and terminals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def add(x, y): return x + y\n",
    "def sub(x, y): return x - y\n",
    "def mul(x, y): return x * y\n",
    "def div(x, y): return 1 if not y else x/y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "FUNCTIONS = [add, sub, mul, div]\n",
    "TERMINALS = ['A', 'B', 'C', 'D']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Initialisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Our usual approach to define the presentation is by introducing code that creates random instances. Let's do the same here, except that this time we need to use a tree representation. For this, we define a suitable datastructure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "class Node:\n",
    "    \n",
    "    def __init__(self, data=None, left=None, right=None):\n",
    "        self.data = data\n",
    "        self.left = left\n",
    "        self.right = right\n",
    "\n",
    "    def __str__(self):\n",
    "        if self.data in FUNCTIONS:\n",
    "            return self.data.__name__\n",
    "        else:\n",
    "            return str(self.data)\n",
    "        \n",
    "    def to_dot(self, dot):\n",
    "        dot.node(str(id(self)), str(self))\n",
    "        \n",
    "        if self.left:\n",
    "            left_dot = self.left.to_dot(dot)\n",
    "            dot.edge(str(id(self)), str(id(self.left)))\n",
    "\n",
    "        if self.right:\n",
    "            right_dot = self.right.to_dot(dot)\n",
    "            dot.edge(str(id(self)), str(id(self.right)))\n",
    "\n",
    "        return dot      "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The `to_dot` function is a helper function we will use in order to represent trees visually using GraphViz (which will need to be installed on your machine if you want to run this Jupyter notebook)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from graphviz import Digraph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def print_tree(tree):\n",
    "    dot = Digraph()\n",
    "    tree.to_dot(dot)\n",
    "    return dot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Full Initialisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "A basic approach to creating random trees is to grow full trees. We decide on a maximum depth for our tree, and then add non-terminal nodes until the depth limit is reached. Once we have reached it, we only add terminals as leaves. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "MAX_DEPTH = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def get_random_terminal():\n",
    "    return random.choice(TERMINALS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def get_random_function():\n",
    "    return random.choice(FUNCTIONS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def fill_random_tree(depth=0):\n",
    "    tree = Node()\n",
    "    if depth < MAX_DEPTH:\n",
    "        tree.data = get_random_function()\n",
    "        tree.left  = fill_random_tree(depth + 1)\n",
    "        tree.right = fill_random_tree(depth + 1)\n",
    "    else:\n",
    "        tree.data = get_random_terminal()\n",
    "    return tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "root = fill_random_tree()\n",
    "print_tree(root)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Grow Initialisation\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "An alternative approach to filling entire trees up to their maximum depth is to grow trees of various sizes. For this we add random nodes while we haven't reached the depth limit. The resulting trees will be less bushy, but the ratio between terminals and non-terminals will bias the average size. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def grow_random_tree(depth=0):\n",
    "    tree = Node()\n",
    "    if depth >= MAX_DEPTH or random.random() > 0.5:\n",
    "        tree.data = get_random_terminal()\n",
    "    else:\n",
    "        tree.data = get_random_function()\n",
    "        tree.left  = grow_random_tree(depth + 1)\n",
    "        tree.right = grow_random_tree(depth + 1)\n",
    "    return tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "root = grow_random_tree()\n",
    "print_tree(root)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "If the ratio of terminals to functions is high, then we might end up with many trees of size 1. To avoid this, we can define a minimum depth for the trees we want to generate, and only consider terminals once the minimum size has been reached:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "MIN_DEPTH = 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def grow_random_tree(depth=0):\n",
    "    tree = Node()\n",
    "    if depth < MIN_DEPTH:\n",
    "        tree.data = get_random_function()\n",
    "    elif depth >= MAX_DEPTH:\n",
    "        tree.data = get_random_terminal()\n",
    "    else:\n",
    "        if random.random() > 0.5:\n",
    "            tree.data = get_random_terminal()\n",
    "        else:\n",
    "            tree.data = get_random_function()\n",
    "                \n",
    "    if tree.data in FUNCTIONS:\n",
    "        tree.left  = grow_random_tree(depth + 1)\n",
    "        tree.right = grow_random_tree(depth + 1)\n",
    "        \n",
    "    return tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "root = grow_random_tree()\n",
    "print_tree(root)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Ramped Half and Half"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "A common practice when creating the initial population is to use a mix of filled and grown trees. For example, half the population may be initialised with the full method, and half with the grow method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "population_size = 40"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def get_initial_population():\n",
    "    population = []\n",
    "    while len(population) < population_size:\n",
    "        if random.random() < 0.5:\n",
    "            population.append(fill_random_tree())\n",
    "        else: \n",
    "            population.append(grow_random_tree())\n",
    "        \n",
    "    return population"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "By defining a helper function that calculates the size of a tree, we can look at the distribution of tree sizes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def size(tree):\n",
    "    if tree.data in FUNCTIONS:\n",
    "        l = size(tree.left) if tree.left else 0\n",
    "        r = size(tree.right) if tree.right else 0\n",
    "        return 1 + l + r\n",
    "    else:\n",
    "        return 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "plt.hist([size(x) for x in get_initial_population()])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    " Sometimes this is also done using a ramped approach, to generate trees of various sizes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def get_initial_population():\n",
    "    population = []\n",
    "    stages = MAX_DEPTH - MIN_DEPTH + 1\n",
    "    for stage in range(MIN_DEPTH, MAX_DEPTH + 1):\n",
    "        for md in range(int(population_size/stages/2)):\n",
    "            population.append(grow_random_tree(MAX_DEPTH - stage))\n",
    "        for md in range(int(population_size/stages/2)):\n",
    "            population.append(fill_random_tree(MAX_DEPTH - stage))\n",
    "    return population"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "plt.hist([size(x) for x in get_initial_population()])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Ephemeral Random Constants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Our terminals only include constants (`A` etc), but in practice we might also need actual values, such as numbers. Including all possible numbers in the set of terminals does not scale well. A common approach is therefore to include _ephemeral random constants_: A single terminal represents the choice of a random value; when the terminal is chosen then it is instantiated with an actual random value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "EPHEMERAL_CONSTANT = \"R\"\n",
    "TERMINALS = ['A', 'B', 'C', 'D', EPHEMERAL_CONSTANT]\n",
    "EPHEMERAL_RANGE = 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def get_random_terminal():\n",
    "    t = random.choice(TERMINALS)\n",
    "    if t == EPHEMERAL_CONSTANT:\n",
    "        t = random.randint(-EPHEMERAL_RANGE, EPHEMERAL_RANGE)\n",
    "    return t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "root = grow_random_tree()\n",
    "print_tree(root)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Search Operators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Since we have a different representation to what we used in the past (which was mainly list-based), we need to adapt our search operators."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Crossover"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "A basic approach to crossing over two trees is to randomly choose two crossover points in the parent trees, and then to cut and swap subtrees below the crossover points. We first need a function with which we can create copies of trees:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def copy(tree):\n",
    "    t = Node()\n",
    "    t.data = tree.data\n",
    "    if tree.left:  t.left = copy(tree.left)\n",
    "    if tree.right: t.right = copy(tree.right)\n",
    "    return t "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We need a couple further helper functions in order to find subtrees at chosen crossover points. First, here's a helper class that allows us to count up or down while (hackily) traversing a tree:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "class Counter:\n",
    "    def __init__(self, num):\n",
    "        self.num = num\n",
    "    \n",
    "    def reduce(self):\n",
    "        self.num -= 1\n",
    "        \n",
    "    def increase(self):\n",
    "        self.num += 1\n",
    "    \n",
    "    def is_target(self):\n",
    "        return self.num == 0\n",
    "    \n",
    "    def greater_zero(self):\n",
    "        return self.num > 0\n",
    "    \n",
    "    def get_value(self):\n",
    "        return self.num"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Using this helper class, we can now traverse the tree and return a copy of the subtree whenever the counter has reached the target:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def get_subtree_at(tree, counter):\n",
    "    counter.reduce()\n",
    "    if counter.is_target():\n",
    "        return copy(tree)\n",
    "    else:\n",
    "        ret = None\n",
    "        if tree.left and counter.greater_zero(): \n",
    "            ret = get_subtree_at(tree.left, counter)\n",
    "        if not ret and tree.right and counter.greater_zero(): \n",
    "            ret = get_subtree_at(tree.right, counter)\n",
    "        return ret"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Similarly, we can replace a target node with an entirely different subtree given a counter that tells us when we've reached the target node:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def insert_subtree_at(tree, subtree, counter):\n",
    "    counter.reduce()\n",
    "    if counter.is_target():\n",
    "        tree.data  = subtree.data\n",
    "        tree.left  = subtree.left\n",
    "        tree.right = subtree.right\n",
    "    else:\n",
    "        if tree.left and counter.greater_zero(): \n",
    "            insert_subtree_at(tree.left, subtree, counter)\n",
    "        if tree.right and counter.greater_zero(): \n",
    "            insert_subtree_at(tree.right, subtree, counter)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Using these two functions we can now define a simple crossover function that picks a random subtree from one parent (`parent1`) and inserts it at a random position in the other parent (`parent2`):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def subtree_crossover(parent1, parent2):\n",
    "    pos1 = random.randint(1, size(parent1))\n",
    "    pos2 = random.randint(1, size(parent2))\n",
    "    \n",
    "    if size(parent1) == 1:\n",
    "        subtree = copy(parent1)\n",
    "    else:\n",
    "        subtree = get_subtree_at(parent1, Counter(pos1))\n",
    "    offspring = copy(parent2)\n",
    "    insert_subtree_at(offspring, subtree, Counter(pos2))\n",
    "    \n",
    "    return offspring"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Let's create some example parents to cross over:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "parent1 = grow_random_tree()\n",
    "print_tree(parent1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "parent2 = grow_random_tree()\n",
    "print_tree(parent2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Now we can produce an offspring by calling our `crossover` function with these two parents:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "offspring = subtree_crossover(parent1, parent2)\n",
    "print_tree(offspring)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Often crossover points are not actually sampled with a uniform random distribution though: Given a tree with branching factor 2 or more, the majority of nodes will be leaves, which crossover will mostly cut a single leaf. A common alternative is to apply a 90% probability of choosing a function node, and only a 10% chance of choosing a terminal node."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Mutation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "\n",
    "Subtree mutation (a.k.a. headless chicken mutation) is a simple mutation where we pick a random subtree, and replace it with a randomly generated subtree:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "P_mutate = 0.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def subtree_mutation(tree):\n",
    "    mutation_point = random.randint(1, size(tree))\n",
    "    random_subtree = grow_random_tree(MAX_DEPTH - 2)\n",
    "    insert_subtree_at(tree, random_subtree, Counter(mutation_point))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "print_tree(parent1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "offspring = copy(parent1)\n",
    "subtree_mutation(offspring)\n",
    "print_tree(offspring)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "An alternative mutation operator is _point mutation_ where, for each node, we replace the node with a certain probability with another node of the same arity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def point_mutation(tree):\n",
    "    if random.random() < P_mutate:\n",
    "        if tree.data in FUNCTIONS: \n",
    "            tree.data = random.choice(FUNCTIONS)\n",
    "        else:\n",
    "            tree.data = get_random_terminal()\n",
    "    \n",
    "    if tree.left:\n",
    "        point_mutation(tree.left)\n",
    "    if tree.right:\n",
    "        point_mutation(tree.right)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "print_tree(parent1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "offspring = copy(parent1)\n",
    "point_mutation(offspring)\n",
    "print_tree(offspring)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "\n",
    "Since selection is independent of the representation, we don't need any adaptations. Since we use a variable size representation, we will include the size as one of the selection criteria."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "tournament_size = 3\n",
    "def tournament_selection(population):\n",
    "    candidates = random.sample(population, tournament_size)        \n",
    "    winner = min(candidates, key=lambda x: (get_fitness(x), size(x)))\n",
    "    return winner"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "An evolution step of our GA now just needs to integrate these operators:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def evolution_step(population):\n",
    "    new_population = []\n",
    "    while len(new_population) < len(population):\n",
    "        parent1 = copy(selection(population))\n",
    "        parent2 = copy(selection(population))\n",
    "\n",
    "        if random.random() < P_xover:\n",
    "            offspring = crossover(parent1, parent2)\n",
    "        else:\n",
    "            offspring = random.choice([parent1, parent2])\n",
    "\n",
    "        if random.random() < 0.5:\n",
    "            point_mutation(offspring)\n",
    "        else:\n",
    "            subtree_mutation(offspring)\n",
    "\n",
    "        new_population.append(offspring)\n",
    "\n",
    "    population.clear()\n",
    "    population.extend(new_population)\n",
    "\n",
    "    best_fitness = min([get_fitness(k) for k in population])\n",
    "    return best_fitness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Similarly, the overall GA is also the same as always:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "fitness_values = []\n",
    "size_values = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def ga():\n",
    "    population = get_initial_population()\n",
    "    best_fitness = sys.maxsize\n",
    "    for p in population:\n",
    "        fitness = get_fitness(p)\n",
    "        if fitness < best_fitness or (fitness == best_fitness and size(p) < size(best_solution)):\n",
    "            best_fitness = fitness\n",
    "            best_solution = copy(p)\n",
    "\n",
    "    iteration = 0\n",
    "    while iteration < max_iterations and best_fitness > 0.000001:\n",
    "        fitness_values.append(best_fitness)\n",
    "        size_values.append(mean([size(x) for x in population]))\n",
    "        print(f\"GA Iteration {iteration}, best fitness: {best_fitness}, average size: {size_values[-1]}\")\n",
    "        iteration += 1\n",
    "        evolution_step(population)\n",
    "\n",
    "        for p in population:\n",
    "            fitness = get_fitness(p)\n",
    "            if fitness < best_fitness or (fitness == best_fitness and size(p) < size(best_solution)):\n",
    "                best_fitness = fitness\n",
    "                best_solution = copy(p)\n",
    "\n",
    "    print(f\"GA solution after {iteration} iterations, best fitness: {best_fitness}\")\n",
    "    return best_solution\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Example Problem: Symbolic Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We will use symbolic regression as a first example problem. Given a set of points, we would like to come up with a symbolic expression that represents a function approximating the points. As example equation, we will use $x^2 + x + 1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We only have one terminal `x` for this problem, but can also include ephemeral constants:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "TERMINALS = [\"x\", EPHEMERAL_CONSTANT]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "As fitness function for how close we are to approximating a set of points we use the sum of absolute errors for all points at different values of `x` in the range [−1.0, +1.0]. Thus, a smaller fitness value is better; a fitness of zero indicates a perfect fit. In order to calculate the error for any given point, we need a helper function that evaluates a tree for a given value of `x`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def evaluate(tree, assignment):\n",
    "    if tree.data in FUNCTIONS:\n",
    "        return tree.data(evaluate(tree.left, assignment), evaluate(tree.right, assignment))\n",
    "    elif tree.data in assignment:\n",
    "        return assignment[tree.data]\n",
    "    else:\n",
    "        return tree.data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The dictionary `assignment` will map all terminals to actual values during the evaluation. We sample a number of points for measuring the error:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "test_data = {}\n",
    "for value in [x/10 for x in range(-10, 10)]:\n",
    "    test_data[value] = value*value + value + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "test_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Now calculating the fitness function reduces to calling the `evaluate` function for each points in our test set, and summing up the differences:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def get_fitness(tree):\n",
    "    fitness = 0.0\n",
    "    \n",
    "    for (x, expected_result) in test_data.items():\n",
    "        assignment = {\"x\": x}\n",
    "        result = evaluate(tree, assignment)\n",
    "        fitness += (result - expected_result) * (result - expected_result)\n",
    "    fitness /= len(test_data)\n",
    "    \n",
    "    return fitness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "With this, we can finally call our genetic algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "max_iterations = 200\n",
    "selection = tournament_selection\n",
    "crossover = subtree_crossover\n",
    "P_xover = 0.7\n",
    "fitness_values = []\n",
    "size_values = []\n",
    "result = ga()\n",
    "print_tree(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Is the solution actually correct? We can plot the results of this function vs the results of the actual target function to compare:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "points = [x/10 for x in range(-10, 11)]\n",
    "real_values = [(x*x + x + 1) for x in points]\n",
    "gp_values = []\n",
    "for x in points:\n",
    "    assignment = {}\n",
    "    assignment[\"x\"] = x\n",
    "    gp_values.append(evaluate(result, assignment))\n",
    "plt.plot(points, gp_values)\n",
    "\n",
    "plt.plot(points, real_values, linestyle='--', label = \"Target\")\n",
    "plt.plot(points, gp_values,  linestyle=':', label = \"GP\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The solution is usually correct; if it is not, try again, and maybe try to adjust the parameters of the GA."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Handling Bloat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The solution is probably a fairly large tree. Indeed the size is sometimes quite problematic as genetic programming tends to suffer from _bloat_. Let's look at the average population size throughout the evolution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "plt.plot(size_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We have already adapted our selection operator to counter bloat somewhat, but mutation and crossover also cause an increase in size. We therefore define alternative mutation and crossover operators."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Shrink mutation replaces a random subtree with a terminal node:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def shrink_mutation(individual):\n",
    "    num_nodes = size(individual)\n",
    "    if num_nodes < 2:\n",
    "        return\n",
    "    \n",
    "    mutation_point = random.randint(2, num_nodes - 1)\n",
    "    node = Node()\n",
    "    node.data = get_random_terminal()\n",
    "    insert_subtree_at(individual, node, Counter(mutation_point))\n",
    "    \n",
    "    return individual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "print_tree(parent1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "offspring = copy(parent1)\n",
    "shrink_mutation(offspring)\n",
    "print_tree(offspring)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Size-fair subtree mutation first picks a random subtree, and then replaces that with a new random subtree that is at most as big as the replaced subtree:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def depth(node):\n",
    "    d = 0\n",
    "    if node.left:\n",
    "        d = max(d, depth(node.left))\n",
    "    if node.right:\n",
    "        d = max(d, depth(node.right))\n",
    "    return d + 1\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def sizefair_subtree_mutation(tree):\n",
    "    mutation_point = random.randint(1, size(tree))\n",
    "    replaced_tree = get_subtree_at(tree, Counter(mutation_point))\n",
    "    random_subtree = grow_random_tree(MAX_DEPTH - depth(replaced_tree) + 1)\n",
    "    insert_subtree_at(tree, random_subtree, Counter(mutation_point))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "print_tree(parent1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "offspring = copy(parent1)\n",
    "sizefair_subtree_mutation(offspring)\n",
    "print_tree(offspring)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Permutation mutation shuffles the order of arguments for functions with arity > 1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def permutation_mutation(tree):\n",
    "    if random.random() < P_mutate:\n",
    "        if tree.data in FUNCTIONS: \n",
    "            tree.left, tree.right = tree.right, tree.left\n",
    "    \n",
    "    if tree.left:\n",
    "        permutation_mutation(tree.left)\n",
    "    if tree.right:\n",
    "        permutation_mutation(tree.right)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "print_tree(parent1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "offspring = copy(parent1)\n",
    "permutation_mutation(offspring)\n",
    "print_tree(offspring)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "There is also a size-fair variant of the crossover operator we defined earlier: We first pick a random crossover point in the first parent. Then, we pick a subtree in the other parent that is not larger than the subtree in the first parent, and perform the crossover with this. We thus define a helper function that gives us all valid positions with subtrees of a maximum size:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def get_smaller_subtree_positions(tree, max_size, positions, counter):\n",
    "    if size(tree) <= max_size:\n",
    "        positions.append(counter.get_value())\n",
    "    counter.increase()\n",
    "    \n",
    "    if tree.left:\n",
    "        get_smaller_subtree_positions(tree.left, max_size, positions, counter)\n",
    "    if tree.right:\n",
    "        get_smaller_subtree_positions(tree.right, max_size, positions, counter)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We apply this helper function by simply randomly picking one of the valid positions during crossover:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def sizefair_crossover(parent1, parent2):\n",
    "    \n",
    "    # First pick a crossover point in parent1\n",
    "    pos1 = 1\n",
    "    if size(parent1) > 1:\n",
    "        pos1 = random.randint(1, size(parent1) - 1)\n",
    "\n",
    "    # Then select a subtree that isn't larger than what is being replaced from parent2\n",
    "    subtree = get_subtree_at(parent1, Counter(pos1))\n",
    "    subtree_size = size(subtree)    \n",
    "    positions = []\n",
    "    get_smaller_subtree_positions(parent2, subtree_size, positions, Counter(1))\n",
    "    pos2 = random.choice(positions)\n",
    "    subtree = get_subtree_at(parent2, Counter(pos2))\n",
    "    \n",
    "    # Then insert that into parent1\n",
    "    offspring = copy(parent1)\n",
    "    insert_subtree_at(offspring, subtree, Counter(pos1))\n",
    "    \n",
    "    assert size(offspring) <= size(parent1)\n",
    "    \n",
    "    return offspring"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "print_tree(parent1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print_tree(parent2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "offspring = sizefair_crossover(parent1, parent2)\n",
    "print_tree(offspring)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Now we just need to do some minor adaptations to make sure that these new operators are actually applied during the evolution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "mutation_operators = [subtree_mutation, permutation_mutation, shrink_mutation, point_mutation, sizefair_subtree_mutation]\n",
    "crossover = sizefair_crossover"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Besides the adapted operators, another common approach is to apply an upper bound on the size, and reject offspring that exceeds that size. This can have negative implications if it is the only means to control bloat, as larger individuals will then become less likely to produce surviving offspring. However, we can use it in conjunction with out optimised operators:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "REJECT_SIZE = 50"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def evolution_step(population):\n",
    "    new_population = []\n",
    "    while len(new_population) < len(population):\n",
    "        parent1 = copy(selection(population))\n",
    "        parent2 = copy(selection(population))\n",
    "\n",
    "        if random.random() < P_xover:\n",
    "            offspring = crossover(parent1, parent2)\n",
    "        else:\n",
    "            offspring = random.choice([parent1, parent2])\n",
    "\n",
    "        if random.random() < 0.2:\n",
    "            mutation_operator = random.choice(mutation_operators)\n",
    "            mutation_operator(offspring)\n",
    "\n",
    "        if size(offspring) < REJECT_SIZE:\n",
    "            new_population.append(offspring)\n",
    "\n",
    "    population.clear()\n",
    "    population.extend(new_population)\n",
    "\n",
    "    best_fitness = min([get_fitness(k) for k in population])\n",
    "    return best_fitness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Thus, finally, let's run the search again for the same target function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "max_iterations = 200\n",
    "population_size = 100\n",
    "selection = tournament_selection\n",
    "crossover = sizefair_crossover\n",
    "P_xover = 0.7\n",
    "fitness_values = []\n",
    "size_values = []\n",
    "result = ga()\n",
    "print_tree(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "To validate whether we are still affected by bloat we can plot the evolution of size again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "plt.plot(size_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We can also compare the resulting function with the desired target function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "points = [x/10 for x in range(-10, 11)]\n",
    "real_values = [(x*x + x + 1) for x in points]\n",
    "gp_values = []\n",
    "for x in points:\n",
    "    assignment = {}\n",
    "    assignment[\"x\"] = x\n",
    "    gp_values.append(evaluate(result, assignment))\n",
    "plt.plot(points, gp_values)\n",
    "\n",
    "plt.plot(points, real_values,  linestyle='--', label = \"Target\")\n",
    "plt.plot(points, gp_values,  linestyle=':', label = \"GP\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Example Application: Fault Localisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Given a test suite with at least one failing and at least one passing test case, the aim of fault localisation is to identify those statements in a program that are likely faulty. Intuitively, statements that are primarily executed by failed test cases are more likely to be faulty than those that are primarily executed by passed test cases. In spectrum-based fault localisation this is done by evaluating the similarity of each statement with the the error vector (i.e. the vector of pass/fail verdicts for each statement). There are many different competing similarity metrics (which we will consider in more details the Software Analysis course). \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The basis for fault localisation is a coverage matrix in which we have coverage information for each statement and each test case. These matrices serve to calculate several factors for each statement:\n",
    "\n",
    "- $e_f$: Number of times the statement has been executed by failing tests.\n",
    "- $e_p$: Number of times the statement has been executed by passing tests.\n",
    "- $n_f$: Number of times the statement has _not_ been executed by failing tests.\n",
    "- $n_p$: Number of times the statement has _not_ been executed by passing tests.\n",
    "\n",
    "Based on these factors, different metrics to calculate the suspiciousness of a program statement have been defined, e.g.:\n",
    "\n",
    "- Tarantula: $\\frac{\\frac{e_f}{e_f+n_f}}{\\frac{e_p}{e_p+n_p} + \\frac{e_f}{e_f+n_f}}$\n",
    "- Barinel: $1 - \\frac{e_p}{e_f + e_p}$\n",
    "- Ochiai: $\\frac{e_f}{\\sqrt{(e_f + n_f) \\cdot (e_f + e_p)}}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The `middle` example function takes three parameters and should return the parameter that is the middle one when ranked by size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def middle(x, y, z):\n",
    "    if y < z:\n",
    "        if x < y:\n",
    "            return y\n",
    "        elif x < z:\n",
    "            return y\n",
    "    else:\n",
    "        if x > y:\n",
    "            return y\n",
    "        elif x > z:\n",
    "            return x\n",
    "    return z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "middle(2, 3, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Our implementation of `middle` is unfortunately buggy, which we can demonstrate by generating some tests and checking the expected values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "tests = []\n",
    "\n",
    "for i in range(10):\n",
    "    x = random.randrange(10)\n",
    "    y = random.randrange(10)\n",
    "    z = random.randrange(10)\n",
    "    m = sorted([x,y,z])[1]\n",
    "    tests.append((x,y,z,m))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Let's check if our tests can trigger a fault, otherwise we will need to generate some more:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "for (x,y,z,m) in tests:\n",
    "    result = middle(x,y,z)\n",
    "    if result != m:\n",
    "        print(f\"Failed test: {x},{y},{z} == {result} but should be {m}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We will consider the program at the level of its lines:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import inspect\n",
    "\n",
    "lines = inspect.getsource(middle).splitlines()\n",
    "for i, line in enumerate(lines):\n",
    "    print(f\"{i}: {line}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In order to apply fault localisation we need to trace test executions and keep track of which lines were executed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "trace = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def trace_line(frame, event, arg):\n",
    "    if event == \"line\":\n",
    "        trace.append(frame.f_lineno)\n",
    "    return trace_line"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def middle_instrumented(x,y,z):\n",
    "    global trace\n",
    "    sys.settrace(trace_line)\n",
    "    trace = []\n",
    "    ret = middle(x,y,z)\n",
    "    sys.settrace(None)\n",
    "    return ret"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "middle_instrumented(1,2,3)\n",
    "trace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "middle_instrumented(3,4,1)\n",
    "trace"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Now we can derive an execution spectrum for the `middle` function using our tests."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def get_spectrum(tests, statements):\n",
    "    matrix = []\n",
    "    for (x,y,z,m) in tests:\n",
    "        row = []\n",
    "        result = middle_instrumented(x,y,z)\n",
    "        for lineno in statements:\n",
    "            if lineno in trace:\n",
    "                row.append(1)\n",
    "            else:\n",
    "                row.append(0)\n",
    "        if result == m:\n",
    "            row.append(1)\n",
    "        else:\n",
    "            row.append(0)\n",
    "        matrix.append(row)\n",
    "    \n",
    "    spectrum = pd.DataFrame(matrix, columns=statements + [\"Passed\"])\n",
    "    return spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "statements = [i for i in range(len(lines))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "middle_spectrum = get_spectrum(tests, statements)\n",
    "middle_spectrum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Let's calculate Barinel using this spectrum (since it's the simplest formula)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def get_suspiciousness(line, spectrum):\n",
    "    \n",
    "    if line not in spectrum.columns:\n",
    "        return 0\n",
    "    \n",
    "    # Number of times the statement has been executed by failing tests.\n",
    "    e_f = spectrum[(spectrum[\"Passed\"] == 0) & (spectrum[line] == 1)].size\n",
    "    \n",
    "    # Number of times the statement has been executed by passing tests.\n",
    "    e_p = spectrum[(spectrum[\"Passed\"] == 1) & (spectrum[line] == 1)].size\n",
    "\n",
    "    if e_p + e_f == 0:\n",
    "        return 0\n",
    "    \n",
    "    suspiciousness = 1 - e_p / (e_p + e_f)\n",
    "        \n",
    "    return suspiciousness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "lines = inspect.getsource(middle).splitlines()\n",
    "for i in range(len(lines)):\n",
    "    line = lines[i]\n",
    "    suspiciousness = get_suspiciousness(i+1, middle_spectrum)\n",
    "    print(\"%2d %.2f: %s\" % (i, suspiciousness, line))\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Indeed the `return y` in line 5 is incorrect."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "To make things more challenging we will parse some example spectra based on some more complex Java classes in the format produced by the GZoltar tool."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def read_matrix(matrix_file, spectra_file, faulty_line):\n",
    "    line_details = []\n",
    "    fh = open(spectra_file)\n",
    "    num_line = 0\n",
    "    for line in fh:\n",
    "        num_line += 1\n",
    "        if num_line == 1:\n",
    "            continue\n",
    "        result = line.rstrip()\n",
    "        line_details.append(result)\n",
    "    fh.close()\n",
    "\n",
    "    num_statements = len(line_details)\n",
    "\n",
    "    ep = {}\n",
    "    ef = {}\n",
    "    np = {}\n",
    "    nf = {}\n",
    "\n",
    "    for i in range(num_statements):\n",
    "        ep[i] = 0\n",
    "        ef[i] = 0\n",
    "        np[i] = 0\n",
    "        nf[i] = 0\n",
    "\n",
    "    fh = open(matrix_file)\n",
    "    num_test = 0\n",
    "    for line in fh:\n",
    "        result = line.split(\" \")\n",
    "\n",
    "        test_result = result[-1].rstrip()\n",
    "        for i in range(num_statements):\n",
    "            if result[i] == '1' and test_result == '+':\n",
    "                ep[i] += 1\n",
    "            elif result[i] == '1' and test_result == '-':\n",
    "                ef[i] += 1\n",
    "            elif result[i] == '0' and test_result == '+':\n",
    "                np[i] += 1\n",
    "            elif result[i] == '0' and test_result == '-':\n",
    "                nf[i] += 1\n",
    "\n",
    "        num_test += 1\n",
    "    fh.close()\n",
    "\n",
    "    return (ep, ef, np, nf, line_details, num_statements, faulty_line)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Our aim is to use GP in order to evolve a similarity function that performs better than the standard metrics. We need some datapoints for our fitness evaluation again; for this we can use some example matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "data_files = [('data/fl/matrix.txt', 'data/fl/spectra.csv', 7),\n",
    "              ('data/fl/matrix-lift.txt', 'data/fl/spectra-lift.csv', 13),\n",
    "              ('data/fl/matrix-complex.txt', 'data/fl/spectra-complex.csv', 39),\n",
    "              ('data/fl/matrix-rational.txt', 'data/fl/spectra-rational.csv', 10)\n",
    "              ]\n",
    "eval_cases = []\n",
    "\n",
    "for (matrix, spectrum, faulty_line) in data_files:\n",
    "    eval_cases.append(read_matrix(matrix, spectrum, faulty_line))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In the following, we will try to use GP in order to derive a formula that produces better results on our dataset than any of these well-established coefficients. This idea is described in more detail in the following paper:\n",
    "\n",
    "Yoo, S. (2012, September). Evolving human competitive spectra-based fault localisation techniques. In International Symposium on Search Based Software Engineering (pp. 244-258). Springer, Berlin, Heidelberg."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Similar to the existing coefficients, the components of our formula shall be simple arithmetic operations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def add(x, y): return x + y\n",
    "def sub(x, y): return x - y\n",
    "def mul(x, y): return x * y\n",
    "def div(x, y): return 1 if not y else x/y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "FUNCTIONS = [add, sub, mul, div]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The terminals will be the factors listed above, as well as the constant 1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "TERMINALS = ['ep', 'ef', 'np', 'nf', 1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In order to evaluate the fitness of a particular formula, we apply it to our four (faulty) example programs, rank the statements by their suspiciousness, and then report the position at which the faulty statement is ranked. A perfect fault localisation means that the faulty statement is listed first. In case of ties, where multiple statements have the same suspiciousness score, we calculate the average position of the statements with the same suspiciousness:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def get_fitness(tree):\n",
    "    positions = []\n",
    "    for (ep, ef, np, nf, line_details, num_statements, faulty_line) in eval_cases:\n",
    "        lines = [line for line in range(num_statements)]\n",
    "        weights = []\n",
    "        for line in lines:\n",
    "            assignment = {\"np\": np[line], \n",
    "                          \"nf\": nf[line], \n",
    "                          \"ep\": ep[line],\n",
    "                          \"ef\": ef[line]}\n",
    "            result = evaluate(tree, assignment)\n",
    "            weights.append(result)\n",
    "        lines, weights = zip(*sorted(zip(lines, weights), key=lambda t: t[1], reverse=True))\n",
    "\n",
    "        start = lines.index(faulty_line)\n",
    "        end = start\n",
    "        faulty_weight = weights[start]\n",
    "        while end < len(weights) and faulty_weight == weights[end]:\n",
    "            end += 1\n",
    "\n",
    "        positions.append((start + end)/2)\n",
    "        \n",
    "    return sum(positions)/len(positions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In order to compare different ranking formulas, we will also look at the top 5 statements for each example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def print_ranking(tree):\n",
    "    for (ep, ef, np, nf, line_details, num_statements, faulty_line) in eval_cases:\n",
    "        lines = [line for line in range(num_statements)]\n",
    "        weights = []\n",
    "        for line in lines:\n",
    "            assignment = {\"np\": np[line], \n",
    "                          \"nf\": nf[line], \n",
    "                          \"ep\": ep[line],\n",
    "                          \"ef\": ef[line]}\n",
    "            result = evaluate(tree, assignment)\n",
    "            weights.append(result)\n",
    "\n",
    "        lines, weights = zip(*sorted(zip(lines, weights), key=lambda t: t[1], reverse=True))\n",
    "\n",
    "        print(f\"Faulty line: {faulty_line}: {lines[:5]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Before evolving our own suspiciousness formula, let's express some standard formulas using our tree datastructure:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Tarantula: $\\frac{\\frac{e_f}{e_f+n_f}}{\\frac{e_p}{e_p+n_p} + \\frac{e_f}{e_f+n_f}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "tarantula = Node(div,\n",
    "                 Node(div, Node('ef'), Node(add, Node('ef'), Node('nf'))),\n",
    "                 Node(add, Node(div, Node('ep'), Node(add, Node('ep'), Node('np'))),\n",
    "                           Node(div, Node('ef'), Node(add, Node('ef'), Node('nf'))))\n",
    "                )\n",
    "print_tree(tarantula)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Barinel: $1 - \\frac{e_p}{e_f + e_p}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "barinel = Node(sub, Node(1), Node(div, Node('ep'), Node(add, Node('ef'), Node('ep'))))\n",
    "print_tree(barinel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Barinel should perform better than Tarantula:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "get_fitness(tarantula)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "get_fitness(barinel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "But can we do better?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "max_iterations = 50\n",
    "population_size = 50\n",
    "selection = tournament_selection\n",
    "crossover = sizefair_crossover\n",
    "P_xover = 0.7\n",
    "fitness_values = []\n",
    "size_values = []\n",
    "result = ga()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "print_tree(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "get_fitness(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Let's now compare the results of the standard metrics and our new one by looking at the top 5 most suspicious statements for each of the coefficients:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "print_ranking(tarantula)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print_ranking(barinel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print_ranking(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The ranking is clearly better overall given these four example spectra. Obviously this isn't a really fair comparison, however, since we are measuring the performance on the training data. For more details of how this generalises, please refer to the original paper."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Strongly-Typed Genetic Programming and Grammatical Evolution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Our implementation of GP is untyped, and thus requires type closure. Sometimes, depending on the problem to be solved, this can be an issue. In _strongly-typed genetic programming_ all nodes have types, and operators need to follow the type rules."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "A common way to implement strongly-typed genetic programming is using _Grammar-Guided Genetic Programming (G3P)_: G3P is based on a grammar that describes the syntax of the target language.  Each individual of the population is a derivation tree of that grammer. (Derivation step = The process of replacing a nonterminal symbol with the symbols on the right-hand side of a production rule). G3P, in principle, works similar to our GP above, with the additional caveat that initialisation, crossover, and mutation need to be adapted in order to use the grammar and to ensure that only grammatically valid derivation trees are used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "An alternative way to use a grammar for evolution is using _Grammatical Evolution_ (GE): Whereas a classic GP uses syntax trees as representation, the representation in grammatical evolution is a simple linear genome. Each individual is a linear structure of codons (numbers) that select production rules to be chosen during derivation from a given grammar. \n",
    "Grammatical evolution clearly separates the search and solution spaces; syntactic correctness is guaranteed by the translation process; the problem of closure is removed, and any search algorithm can be used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Since we are using a list representation, we can use our list helper class that can cache fitness values again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "class L(list):\n",
    "    \"\"\"\n",
    "    A subclass of list that can accept additional attributes.\n",
    "    Should be able to be used just like a regular list.\n",
    "    \"\"\"\n",
    "    def __new__(self, *args, **kwargs):\n",
    "        return super(L, self).__new__(self, args, kwargs)\n",
    "\n",
    "    def __init__(self, *args, **kwargs):\n",
    "        if len(args) == 1 and hasattr(args[0], '__iter__'):\n",
    "            list.__init__(self, args[0])\n",
    "        else:\n",
    "            list.__init__(self, args)\n",
    "        self.__dict__.update(kwargs)\n",
    "\n",
    "    def __call__(self, **kwargs):\n",
    "        self.__dict__.update(kwargs)\n",
    "        return self"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In GE the chromosome size is fixed. For example, let's start by setting the number of codons to 50. Each codon is a number in a given range, which we will encode using `MAX`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "num_codons = 50\n",
    "MAX = 50"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Individuals are simply lists of (randomly generated) codons:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def get_random_individual():\n",
    "    individual = L([random.randint(0, MAX) for _ in range(num_codons)])\n",
    "    evaluate(individual)\n",
    "    return individual"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "This, in turn, means we can use our usual approach for mutating each gene in a chromosome with a probability dependent on the overall number of genes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def mutate(individual):\n",
    "    P_mutate = 1/len(individual)\n",
    "    mutated = L(individual[:])\n",
    "    for pos in range(len(individual)):\n",
    "        if random.random() < P_mutate:\n",
    "            mutated[pos] = random.randint(0, MAX)\n",
    "    evaluate(mutated)\n",
    "    return mutated"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We can use standard crossover operators:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def singlepoint_crossover(parent1, parent2):\n",
    "    pos = random.randint(0, len(parent1))\n",
    "    offspring1 = L(parent1[:pos] + parent2[pos:])\n",
    "    offspring2 = L(parent2[:pos] + parent1[pos:])\n",
    "    return (offspring1, offspring2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "...and we can also use standard parent and survivor selection operators:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "elite_size = 2\n",
    "def elitism(population):\n",
    "    population.sort(key=lambda k: get_fitness(k))\n",
    "    return population[:elite_size]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "tournament_size = 3\n",
    "def tournament_selection(population):\n",
    "    candidates = random.sample(population, tournament_size)        \n",
    "    winner = min(candidates, key=lambda x: get_fitness(x))\n",
    "    return winner"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In order to explain the process of decoding codons to production rules, we need to define a grammar for an example problem. Let's assume we want to generate Python code for a function with two parameters, which returns the larger of the two parameters. Our grammar thus covers linear sequences of code with if-statements, where we restrict each if-statement to a single operation. Operations can be updates of a local variable `__return__` which, at the end of the function, will define the return value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "EXAMPLE_GRAMMAR = {\n",
    "    \"<start>\":\n",
    "        [\"<code>\"],\n",
    "\n",
    "    \"<code>\":\n",
    "        [\"<line>\", \"<code>\\n<line>\"],\n",
    "\n",
    "    \"<line>\":\n",
    "        [\"    <expr>\"],\n",
    "\n",
    "    \"<expr>\":\n",
    "        [\"<if-statement>\",\n",
    "         \"<op>\"],\n",
    "\n",
    "    \"<if-statement>\":\n",
    "        [\"if <condition>:\\n        <op>\"],\n",
    "\n",
    "    \"<condition>\":\n",
    "        [\"<param> <operator> <param>\"],\n",
    "    \"<param>\":\n",
    "        [\"x\", \"y\"],\n",
    "    \"<operator>\":\n",
    "        [\"<\", \">\", \"==\"],\n",
    "    \"<op>\":\n",
    "        [\"__return__ = x\", \"__return__ = y\"]\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "One of the production rules must be explicitly chosen as the start label\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "START_SYMBOL = \"<start>\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In order to determine whether a production rule contains further non-terminals we use a simple regular expression. Anything within `<` and `>` counts as a non-terminal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "RE_NONTERMINAL = re.compile(r'(<[^<> ]*>)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def nonterminals(expansion):\n",
    "    return re.findall(RE_NONTERMINAL, expansion)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Depending on the grammar used, a derivation can be of any size; we need to put some limit on the size of the derivations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "max_expansions = 50"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In order to decode an individual we start with the start rule, and apply derivation. Whenever there is a choice of multiple production rules during the derivation, the number of the rule is determined by the next codon; i.e., we always pick the rule with the number _codon_ $%$ _number of rules_. Non-terminals are resolved from left to right (we need to make sure that the derivation is deterministic). Whenever we use a codon, we move to the next number in the chromosome. If we hit the last codon in a chromosome then we wrap and start over from the first one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def decode(individual, grammar):\n",
    "    pos = 0\n",
    "    term = START_SYMBOL\n",
    "    num_expansions = 0\n",
    "    count = 0\n",
    "    \n",
    "    while len(nonterminals(term)) > 0 and num_expansions < max_expansions:\n",
    "        \n",
    "        symbol_to_expand = nonterminals(term)[0]\n",
    "        expansions = grammar[symbol_to_expand]\n",
    "        if len(expansions) == 1:\n",
    "            expansion = expansions[0]\n",
    "        else:\n",
    "            codon = individual[pos]\n",
    "            pos = (pos + 1) % len(individual)\n",
    "            expansion = expansions[codon % len(expansions)]\n",
    "            \n",
    "        new_term = term.replace(symbol_to_expand, expansion, 1)\n",
    "    \n",
    "        term = new_term\n",
    "        num_expansions += 1\n",
    "\n",
    "    return term"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "print(decode([1,2,3], EXAMPLE_GRAMMAR))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Given the string representation of a program we need to create an executable version in order to measure the fitness. We already compiled code at runtime in the previous chapter, so we essentially do the same again here. For any given string representation of a function, e.g.:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "text = \"\"\"\n",
    "def foo():\n",
    "  print(\"Foo\")\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We call the `compile` function which produces a code object. This code object we insert into the dictionary of modules such that we can call the function using its name, e.g.:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "code = compile(text, filename=\"<GE>\", mode=\"exec\")\n",
    "current_module = sys.modules[__name__]\n",
    "exec(code, current_module.__dict__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "foo()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We use a function that puts the decoded chromosome into a function, and we'll include the function signature and return value for now:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def create_function(individual):\n",
    "    program_text = decode(individual, EXAMPLE_GRAMMAR)\n",
    "    source = \"\"\"def generated_function(x, y):\n",
    "    __return__ = 0\n",
    "%s\n",
    "    return __return__\n",
    "\"\"\" % program_text\n",
    "    code = compile(source, filename=\"<GE>\", mode=\"exec\")\n",
    "    current_module = sys.modules[__name__]\n",
    "    exec(code, current_module.__dict__)\n",
    "    return generated_function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In order to calculate the fitness value we just need to run the function a number of times with some test data that represents the function we want to learn. Let's produce some test data again; let our example function return the maximum of the two input parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "training_data = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "for i in range(20):\n",
    "    num1 = random.randint(0, MAX)\n",
    "    num2 = random.randint(0, MAX)\n",
    "    training_data.append(([num1, num2], max([num1, num2])))\n",
    "\n",
    "training_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "There may be errors during the execution; for example, if we hit the maximum number of expansions (`max_expansions`) then the code may not compile. There may be an exception if there is a division by zero or some other error. In these cases we assign the function the maximum (worst) fitness value and let selective pressure handle this for us:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def evaluate(individual):\n",
    "    fitness = len(training_data)\n",
    "    individual.fitness = fitness\n",
    "    try:\n",
    "        p = create_function(individual)\n",
    "\n",
    "        for (test_data, test_output) in training_data:\n",
    "            result = p(*test_data)\n",
    "            if result == test_output:\n",
    "                fitness -= 1\n",
    "    except:\n",
    "        return len(training_data)\n",
    "\n",
    "    individual.fitness = fitness\n",
    "    return fitness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "x = get_random_individual()\n",
    "x.fitness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Since we call the `evaluate` function directly after generating or mutating individuals, retrieving the fitness function just consists of retrieving the cached fitness value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def get_fitness(individual):\n",
    "    return individual.fitness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In order to use our new individuals we need to modify our GA such that it uses this new representation rather than the tree structure we defined earlier; that is, we need to replace the ramped half and half initialisation with a simple random sampling:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def get_initial_population():\n",
    "    return [get_random_individual() for _ in range(population_size)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Furthermore, since we now have a fixed size representation we can omit the handling of size as a secondary criterion during selection:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def evolution_step(population):\n",
    "    new_population = []\n",
    "    while len(new_population) < len(population):\n",
    "        parent1 = selection(population)\n",
    "        parent2 = selection(population)\n",
    "\n",
    "        if random.random() < P_xover:\n",
    "            offspring1, offspring2 = crossover(parent1, parent2)\n",
    "        else:\n",
    "            offspring1, offspring2 = parent1, parent2\n",
    "\n",
    "        offspring1 = mutate(offspring1)\n",
    "        offspring2 = mutate(offspring2)\n",
    "        \n",
    "        new_population.append(offspring1)\n",
    "        new_population.append(offspring2)\n",
    "\n",
    "    population.clear()\n",
    "    population.extend(new_population)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def ga():\n",
    "    population = get_initial_population()\n",
    "    best_fitness = sys.maxsize\n",
    "    for p in population:\n",
    "        fitness = get_fitness(p)\n",
    "        if fitness < best_fitness:\n",
    "            best_fitness = fitness\n",
    "            best_solution = L(p[:])\n",
    "\n",
    "    iteration = 0\n",
    "    print(f\"GA Iteration {iteration}, best fitness: {best_fitness}\")\n",
    "    while iteration < max_iterations and best_fitness > 0:\n",
    "        fitness_values.append(best_fitness)\n",
    "        iteration += 1\n",
    "        evolution_step(population)\n",
    "\n",
    "        for p in population:\n",
    "            fitness = get_fitness(p)\n",
    "            if fitness < best_fitness:\n",
    "                best_fitness = fitness\n",
    "                best_solution = L(p[:])\n",
    "                print(f\"GA Iteration {iteration}, best fitness: {best_fitness}\")\n",
    "\n",
    "    print(f\"GA solution after {iteration} iterations, best fitness: {best_fitness}\")\n",
    "    return best_solution\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "With this, we can apply the grammatical evolution to produce our own max function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "population_size = 20\n",
    "P_xover = 0.7\n",
    "max_iterations = 100\n",
    "num_codons = 20\n",
    "selection = tournament_selection\n",
    "crossover = singlepoint_crossover\n",
    "fitness_values = []\n",
    "result = ga()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(decode(result, EXAMPLE_GRAMMAR))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "This problem is probably slightly too easy, so let's increase difficulty by expanding to three parameters. We therefore require a corresponding function template:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def create_function(individual):\n",
    "    program_text = decode(individual, EXAMPLE_GRAMMAR)\n",
    "    source = \"\"\"def generated_function(x, y, z):\n",
    "    __return__ = 0\n",
    "%s\n",
    "    return __return__\n",
    "\"\"\" % program_text\n",
    "    code = compile(source, filename=\"<GE>\", mode=\"exec\")\n",
    "    current_module = sys.modules[__name__]\n",
    "    exec(code, current_module.__dict__)\n",
    "    return generated_function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "As well as some test data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "training_data = []\n",
    "for i in range(1000):\n",
    "    num1 = random.randint(0, MAX)\n",
    "    num2 = random.randint(0, MAX)\n",
    "    num3 = random.randint(0, MAX)\n",
    "    \n",
    "    training_data.append(([num1, num2, num3], max([num1, num2, num3])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We'll slightly expand our grammar to also include conjunctions and disjunctions and more relational operators:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "EXAMPLE_GRAMMAR = {\n",
    "    \"<start>\":\n",
    "        [\"<code>\"],\n",
    "\n",
    "    \"<code>\":\n",
    "        [\"<line>\", \"<code>\\n<line>\"],\n",
    "\n",
    "    \"<line>\":\n",
    "        [\"    <expr>\"],\n",
    "\n",
    "    \"<expr>\":\n",
    "        [\"<if-statement>\",\n",
    "         \"<op>\"],\n",
    "\n",
    "    \"<if-statement>\":\n",
    "        [\"if <condition>:\\n        <op>\"],\n",
    "\n",
    "    \"<condition>\":\n",
    "        [\"<param> <operator> <param>\", \"<param> <operator> <param> and <param> <operator> <param>\", \"<param> <operator> <param> or <param> <operator> <param>\"],\n",
    "\n",
    "    \"<param>\":\n",
    "        [\"x\", \"y\", \"z\", \"__return__\"],\n",
    "\n",
    "    \"<operator>\":\n",
    "        [\"<\", \">\", \"<=\", \">=\", \"==\", \"!=\"],\n",
    "\n",
    "    \"<op>\":\n",
    "        [\"__return__ = <param>\"]\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "population_size = 100\n",
    "P_xover = 0.7\n",
    "max_iterations = 500\n",
    "max_expansions = 80\n",
    "num_codons = 80\n",
    "selection = tournament_selection\n",
    "tournament_size = 4\n",
    "crossover = singlepoint_crossover\n",
    "fitness_values = []\n",
    "result = ga()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "print(decode(result, EXAMPLE_GRAMMAR))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "create_function(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "generated_function(1,2,3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Finding the maximum of three numbers appears to be relatively easy, and the search often finds a solution. Let's try and see if we can also reproduce our `middle` program. For this, we need new tests."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "training_data = []\n",
    "for i in range(1000):\n",
    "    num1 = random.randint(0, MAX)\n",
    "    num2 = random.randint(0, MAX)\n",
    "    num3 = random.randint(0, MAX)\n",
    "    \n",
    "    training_data.append(([num1, num2, num3], sorted([num1, num2, num3])[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "population_size = 100\n",
    "P_xover = 0.7\n",
    "max_iterations = 500\n",
    "max_expansions = 80\n",
    "num_codons = 80\n",
    "selection = tournament_selection\n",
    "tournament_size = 10\n",
    "crossover = singlepoint_crossover\n",
    "fitness_values = []\n",
    "result = ga()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(decode(result, EXAMPLE_GRAMMAR))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "This appears to be a somewhat more challenging problem, and the search tends to take much longer to find an acceptable solution (probably longer than is reasonable to wait in this Jupyter notebook)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "More practical information about Genetic Programming can be found in the free ebook [A Field Guide to Genetic Programming](http://www0.cs.ucl.ac.uk/staff/W.Langdon/ftp/papers/poli08_fieldguide.pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Genetic Improvement"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The programs we tried to evolve in this chapter have all been comparatively tiny. Indeed GPs are not generally used to evolve large programs, as they cannot easily come up with complex system designs and architectures. While evolving large programs from scratch may currently not be possible, it is possible to improve existing programs. This is the aim of program repair and _Genetic Improvement_ (GI). There are two central hypotheses that motivate Genetic Improvement:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The _Competent Programmer Hypothesis_ assumes that programmers write programs that are close to being correct, and program faults are syntactically small."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The _Plastic Surgery Hypothesis_ states that the content of new code can often be assembled out of fragments of code that already exist. Changes to a codebase contain snippets that already exist in the codebase at the time of the change, and these snippets can be efficiently found and exploited."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Given these two hypotheses, the idea of Genetic Improvement is to use GP to slightly modify a large, existing, code base, rather than to evolve something from scratch. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The representation is typically based either in patches or modifications of the abstract syntax tree (AST) of the program. Thus, in contrast to the past examples where we implicitly defined our encoding through the `get_random_individual` function, in this case the representation is already defined by the AST of the programs we are trying to evolve.  For example, let's consider the AST of the `middle` function we considered earlier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import ast"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "middle_ast = ast.parse(inspect.getsource(middle))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "ast.dump(middle_ast)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We know from earlier that the `middle` function contains a fault. The classic application of Genetic Improvement is _Automated Program Repair_: Given a program and a test suite in which there is at least one failing test, try to evolve the program such that it passes all test cases. This idea was initially described in the following seminal paper:\n",
    "\n",
    "Le Goues, C., Nguyen, T., Forrest, S., & Weimer, W. (2011). Genprog: A generic method for automatic software repair. IEEE Transactions on Software Engineering, 38(1), 54-72."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Unlike the GP approaches we defined earlier, the starting point for GI is not a population of random individuals, but the given program we want to improve. Therefore, we will not start by defining a `get_random_individual` function this time, but first introduce mutations on an AST, and then the random individuals will be a result of mutating the existing program. As the mutation will be dependent on executions, we first consider the fitness function. The search in automated program repair is guided by tests: Given a set of tests, we want to maximise the number of passing tests. Let's define some passing and failing tests for our defective `middle` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "PASSING_TESTS = []\n",
    "FAILING_TESTS = []\n",
    "\n",
    "num_tests = 20\n",
    "\n",
    "while len(PASSING_TESTS) + len(FAILING_TESTS) < 2 * num_tests:\n",
    "    x = random.randrange(10)\n",
    "    y = random.randrange(10)\n",
    "    z = random.randrange(10)\n",
    "    m = middle(x, y, z) # Actual value\n",
    "    expected = sorted([x, y, z])[1] # Expected value\n",
    "    if m == expected:\n",
    "        if len(PASSING_TESTS) < num_tests:\n",
    "            PASSING_TESTS.append((x, y, z, m))\n",
    "    else:\n",
    "        if len(FAILING_TESTS) < num_tests:\n",
    "            FAILING_TESTS.append((x, y, z, expected))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "PASSING_TESTS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "FAILING_TESTS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The fitness function is a weighted combination of the passed tests -- intuitively, it is more progress to fix a failing test than to break a passing test, so we assign a higher weight to making failing tests pass."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "WEIGHT_PASSING = 1\n",
    "WEIGHT_FAILING = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "To calculate the fitness value, we need to compile the AST for a given individual, execute all our tests, and count how many passing tests there are."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def get_fitness(individual):\n",
    "    global PASSING_TESTS, FAILING_TESTS\n",
    "    original_middle = middle\n",
    "\n",
    "    code = compile(individual, '<fitness>', 'exec')\n",
    "    exec(code, globals())\n",
    "\n",
    "    passing_passed = 0\n",
    "    failing_passed = 0\n",
    "\n",
    "    for x, y, z, m in PASSING_TESTS:\n",
    "        if middle(x, y, z) == m:\n",
    "            passing_passed += 1\n",
    "\n",
    "    for x, y, z, m in FAILING_TESTS:\n",
    "        if middle(x, y, z) == m:\n",
    "            failing_passed += 1\n",
    "\n",
    "    fitness = WEIGHT_PASSING * passing_passed + WEIGHT_FAILING * failing_passed\n",
    "        \n",
    "    globals()['middle'] = original_middle\n",
    "    return fitness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "get_fitness(middle_ast)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "When mutating individuals we consider the level of statements. We therefore need to extract all statements from a program in order to decide where to mutate, and how to mutate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "class StatementVisitor(ast.NodeVisitor):\n",
    "\n",
    "    def __init__(self):\n",
    "        self.statements = []\n",
    "        self.statements_seen = set()\n",
    "        super().__init__()\n",
    "\n",
    "    def add_statements(self, node, attr):\n",
    "        elems = getattr(node, attr, [])\n",
    "        if not isinstance(elems, list):\n",
    "            elems = [elems]\n",
    "\n",
    "        for stmt in elems:\n",
    "            if stmt in self.statements_seen:\n",
    "                continue\n",
    "\n",
    "            self.statements.append(stmt)\n",
    "            self.statements_seen.add(stmt)\n",
    "\n",
    "    def visit_node(self, node):\n",
    "        # Any node other than the ones listed below\n",
    "        self.add_statements(node, 'body')\n",
    "        self.add_statements(node, 'orelse')\n",
    "\n",
    "    def visit_Module(self, node: ast.Module):\n",
    "        # Module children are defs, classes and globals - don't add\n",
    "        super().generic_visit(node)\n",
    "        \n",
    "    def generic_visit(self, node):\n",
    "        self.visit_node(node)\n",
    "        super().generic_visit(node)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def all_statements(tree, tp = None):\n",
    "    visitor = StatementVisitor()\n",
    "    visitor.visit(tree)\n",
    "    statements = visitor.statements\n",
    "    if tp is not None:\n",
    "        statements = [s for s in statements if isinstance(s, tp)]\n",
    "\n",
    "    return statements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "all_statements(middle_ast)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "This allows us to choose random nodes in our AST."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import astor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "random_node = random.choice(all_statements(middle_ast))\n",
    "print(astor.to_source(random_node))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "GenProg defines three types of mutations: AST-Nodes can be _swapped_, _inserted_, or _deleted_. There are two important aspects to this mutation:\n",
    "\n",
    "- First, the location at which an insertion, swap, or deletion happens depends on the suspiciousness of the nodes as determined using fault localisation. This way, the search should focus on code that is likely faulty.\n",
    "\n",
    "- Second, for the chosen location the mutations do not use random code, but try to re-use the existing code, as the assumption is that the relevant code is likely already contained in the program.\n",
    "\n",
    "In order to implement these mutations, we define a helper function that replaces AST nodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import copy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "class NodeReplacer(ast.NodeTransformer):\n",
    "\n",
    "    def __init__(self, target_node, mutant_node):\n",
    "        self.target_node = target_node\n",
    "        self.mutant_node = mutant_node\n",
    "\n",
    "    def visit(self, node):\n",
    "        if node == self.target_node:\n",
    "            return self.mutant_node\n",
    "        else:\n",
    "            super().visit(node)\n",
    "            return node"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "For the insertion mutation we need a helper function that inserts a new node after a chosen target `node`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def mutate_insert(node, statements):\n",
    "    new_node = copy.deepcopy(random.choice(statements))\n",
    "\n",
    "    if isinstance(node, ast.Return):\n",
    "        if isinstance(new_node, ast.Return):\n",
    "            return new_node\n",
    "        else:\n",
    "            return [new_node, node]\n",
    "\n",
    "    return [node, new_node]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Swapping is very similar, except that we drop the original node."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def mutate_swap(node, statements):     \n",
    "    replacement_node = copy.deepcopy(random.choice(statements))\n",
    "    return replacement_node"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Finally, for deletion we remove a node by setting it to `None`, unless that would make the program syntactically invalid, in which case in Python we need a `pass` statement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def mutate_delete(node, statements):\n",
    "    branches = [attr for attr in ['body', 'orelse', 'finalbody'] if hasattr(node, attr) and getattr(node, attr)]\n",
    "    if branches:\n",
    "        # Replace `if P: S` by `S`\n",
    "        branch = random.choice(branches)\n",
    "        new_node = getattr(node, branch)\n",
    "        return new_node\n",
    "\n",
    "    if isinstance(node, ast.stmt):\n",
    "        # Avoid empty bodies; make this a `pass` statement\n",
    "        new_node = ast.Pass()\n",
    "        ast.copy_location(new_node, node)\n",
    "        return new_node\n",
    "\n",
    "    return None  # Just delete"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The final ingredient is the selection of nodes, which we make dependent on the weights as defined by the fault localisation. Note that our implementation redundantly executes all tests, this could of course be optimised."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def select_node(individual):\n",
    "    tests = PASSING_TESTS + FAILING_TESTS\n",
    "    statements = all_statements(individual)\n",
    "    spectrum = get_spectrum(tests, statements)\n",
    "\n",
    "    if len(statements) < 3:\n",
    "        return random.choice(statements)\n",
    "    \n",
    "    # Crude hack for probabilistic choice\n",
    "    candidates = random.sample(statements, 3)\n",
    "    return max(candidates, key=lambda x: get_suspiciousness(x.lineno-1, spectrum))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Now we have all ingredients in place, and can implement the overall mutation operator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def mutate(individual):\n",
    "    mutant = copy.deepcopy(individual)\n",
    "    statements = all_statements(mutant)\n",
    "\n",
    "    node = select_node(mutant)\n",
    "\n",
    "    mutation = random.choice([mutate_insert, mutate_swap, mutate_delete])\n",
    "    mutated_node = mutation(node, statements)\n",
    "    \n",
    "    mutator = NodeReplacer(node, mutated_node)\n",
    "    mutant = mutator.visit(mutant)\n",
    "\n",
    "    ast.fix_missing_locations(mutant)\n",
    "    return mutant "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Our mutation operator only picks a single node and applies a mutation. In the original GenProg implementation, nodes are probabilistically mutated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "print(astor.to_source(mutate(middle_ast)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "This is also exactly how the random individuals are created: Not by creating random ASTs, but by applying a mutation to the original program."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def get_random_individual():\n",
    "    node = mutate(middle_ast)\n",
    "    return node"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(astor.to_source(get_random_individual()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We need to make sure that our initial population uses this function, rather than the ramped half and half we previously defined."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def get_initial_population():\n",
    "    return [copy.deepcopy(middle_ast)] + [get_random_individual() for _ in range(population_size - 1)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Given our past implementations and the mutation operators we just defined, the crossover is quite simple: We select nodes to split, and then re-attach the programs. Note that this is once again a simplification compared to GenProg; for example, GenProg only considers 'suspicious' statements as possible crossover locations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def crossover(parent1, parent2):\n",
    "    statements1 = all_statements(parent1)\n",
    "    xover_point_1 = random.choice(statements1)\n",
    "\n",
    "    statements2 = all_statements(parent2)\n",
    "    xover_point_2 = random.choice(statements2)\n",
    "    \n",
    "    mutator = NodeReplacer(xover_point_1, copy.deepcopy(xover_point_2))\n",
    "    offspring1 = mutator.visit(parent1)\n",
    "\n",
    "    mutator = NodeReplacer(xover_point_2, copy.deepcopy(xover_point_1))\n",
    "    offspring2 = mutator.visit(parent2)\n",
    "\n",
    "    return offspring1, offspring2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "parent1 = get_random_individual()\n",
    "parent2 = get_random_individual()\n",
    "print(astor.to_source(parent1))\n",
    "print(astor.to_source(parent2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "offspring1, offspring2 = crossover(parent1, parent2)\n",
    "print(astor.to_source(offspring1))\n",
    "print(astor.to_source(offspring2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "While the selection is independent of the representation, we do need to take into account that we now have a maximisation problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def get_size(individual):\n",
    "    return sum([1 for _ in ast.walk(individual)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "tournament_size = 3\n",
    "def tournament_selection(population):\n",
    "    candidates = random.sample(population, tournament_size)\n",
    "    winner = max(candidates, key=lambda x: get_fitness(x))\n",
    "    return winner"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "GenProg makes several adaptations to the overall genetic algorithm; we'll just adapt our existing standard implementation to make use of the new operators."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def ga():\n",
    "    population = get_initial_population()\n",
    "    best_fitness = 0\n",
    "    for p in population:\n",
    "        fitness = get_fitness(p)\n",
    "        if fitness > best_fitness:\n",
    "            best_fitness = fitness\n",
    "            best_solution = copy.deepcopy(p)\n",
    "\n",
    "    iteration = 0\n",
    "    while iteration < max_iterations and \\\n",
    "            best_fitness < WEIGHT_PASSING * len(PASSING_TESTS) + WEIGHT_FAILING * len(FAILING_TESTS):\n",
    "        fitness_values.append(best_fitness)\n",
    "        print(f\"GA Iteration {iteration}, best fitness: {best_fitness}\")\n",
    "        iteration += 1\n",
    "        evolution_step(population, best_solution)\n",
    "\n",
    "        for p in population:\n",
    "            fitness = get_fitness(p)\n",
    "            if fitness > best_fitness:\n",
    "                best_fitness = fitness\n",
    "                best_solution = copy.deepcopy(p)\n",
    "\n",
    "    print(f\"GA solution after {iteration} iterations, best fitness: {best_fitness}\")\n",
    "    return best_solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Since our mutation operator will always change exactly one node, we need to add a probability for mutation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "P_mutation = 0.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def evolution_step(population, elite):\n",
    "    new_population = [elite, get_random_individual()]\n",
    "    viable_population = [p for p in population if get_fitness(p) > 0]\n",
    "    while len(new_population) < len(population):\n",
    "        parent1 = copy.deepcopy(selection(viable_population))\n",
    "        parent2 = copy.deepcopy(selection(viable_population))\n",
    "\n",
    "        offspring1, offspring2 = crossover(copy.deepcopy(parent1), copy.deepcopy(parent2))\n",
    "\n",
    "        if random.random() < P_mutation:\n",
    "            offspring1 = mutate(offspring1)\n",
    "            offspring2 = mutate(offspring2)\n",
    "\n",
    "        new_population.append(offspring1)\n",
    "        new_population.append(offspring2)\n",
    "        new_population.append(parent1)\n",
    "        new_population.append(parent2)\n",
    "\n",
    "    population.clear()\n",
    "    population.extend(new_population)\n",
    "\n",
    "    best_fitness = max([get_fitness(k) for k in population])\n",
    "    return best_fitness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Can we fix the `middle` program?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "P_xover = 0.0\n",
    "P_mutation = 0.4\n",
    "max_iterations = 50\n",
    "population_size = 40\n",
    "tournament_size = 5\n",
    "fitness_values = []\n",
    "result = ga()\n",
    "\n",
    "print(astor.to_source(result))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "So how well does the attempted fix work? Let's print out all the failing tests from our test set, if there are any."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "code = compile(result, '<fitness>', 'exec')\n",
    "exec(code, globals())\n",
    "\n",
    "for x, y, z, m in PASSING_TESTS:\n",
    "    if middle(x, y, z) != m:\n",
    "        print(f\"Failed {x}, {y}, {z} == {m}\")\n",
    "\n",
    "for x, y, z, m in FAILING_TESTS:\n",
    "    if middle(x, y, z) != m:\n",
    "        print(f\"Failed {x}, {y}, {z} == {m}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We've taken several shortcuts in this implementation; the original GenProg implementation has some differences in terms of the genetic operators and the search algorithm. A more elaborate discussion and implementation of the can be found in Andreas Zeller's Debugging Book: https://www.debuggingbook.org/"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.5"
  },
  "rise": {
   "enable_chalkboard": true,
   "scroll": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}