Evolve continuous mountain car.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evolutionary Robotics\n",
    "\n",
    "In this notebook, we will discuss evolutionary robotics and have you implement an evolutionary algorithm to solve one of the tasks in the AI gym, namely the continuous mountain car task. As you will see, an initial implementation of an evolutionary algorithm for solving a robotics task is easily made. However, obtaining good results may - depending on the task - be hard, and understanding the solution may even be harder still.\n",
    "\n",
    "The figure below shows the typical evolutionary robotics approach. An initial population is randomly generated. Then there is an iterative process of: (1) Evaluating all individuals (genomes) in the population, resulting in a fitness value for each individual, (2) Selecting the individuals that will be allowed to procreate, i.e., form the new generation, and (3) Vary on the genomes of the selected individuals (using cross-over, mutation, etc.). The process typically terminates either after a specified number of generations, or after convergence to an optimal solution. Evaluation involves the conversion of the genome (genotype) to the phenotype (e.g., setting the weights of a neural network to the values in the genome). Then the phenotype is tested out on the task, typically in simulation but in some works also on real robots. In robotics tasks, evaluation is a stochastic process and execution of the task by the robot can take a long time. \n",
    "\n",
    "<img src=\"evolutionary_robotics_process.jpg\" width=\"50%\"></img>\n",
    "*Figure 1:* Depiction of the typical evolutionary robotics approach. Figure from: _Doncieux, S., Bredeche, N., Mouret, J. B., & Eiben, A. E. G. (2015). Evolutionary robotics: what, why, and where to. Frontiers in Robotics and AI, 2, 4._\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MountainCarContinuous-v0\n",
    "\n",
    "In this notebook, you will apply an evolutionary robotics approach to the continuous mountain car task. In this task, the car needs to reach the flag on the right mountain, while observing its position and velocity, and acting by means of accelerating the car left or right. The car cannot go straight up the mountain, but has to build up momentum to succeed. The fitness function rewards reaching the flag, and punishes the use of control actions (using less energy to reach the hill top is better). Please see the general description of the task <A HREF=\"https://gym.openai.com/envs/MountainCarContinuous-v0/\" TARGET=\"_blank\">here</A> and the details of the task <A HREF=\"https://github.com/openai/gym/wiki/MountainCarContinuous-v0\" TARGET=\"_blank\">here</A>. Most importantly, it has two observations, i.e., the position and velocity along the trajectory, and one control input, i.e., the applied force, which has to be in the interval $[-1, 1]$.\n",
    "\n",
    "<img src=\"continuous_mountain_car.png\" width=\"50%\"></img>\n",
    "*Figure 2:* Screenshot of the continuous mountain car task. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below you find code to evaluate an agent a single time in the mountain car environment. Please study the code, and note that the method ```act``` should be replaced in the end with a learned controller."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'CMC_original' object has no attribute 'seed'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[1], line 16\u001b[0m\n\u001b[0;32m     13\u001b[0m         \u001b[38;5;28;01mreturn\u001b[39;00m [\u001b[38;5;241m2.0\u001b[39m\u001b[38;5;241m*\u001b[39mnp\u001b[38;5;241m.\u001b[39mrandom\u001b[38;5;241m.\u001b[39mrand()\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1.0\u001b[39m]\n\u001b[0;32m     15\u001b[0m agent \u001b[38;5;241m=\u001b[39m random_agent()\n\u001b[1;32m---> 16\u001b[0m reward \u001b[38;5;241m=\u001b[39m \u001b[43mrun_cart\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrun_cart_continuous\u001b[49m\u001b[43m(\u001b[49m\u001b[43magent\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43menv\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mrun_cart\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mCMC_original\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mgraphics\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m)\u001b[49m\n\u001b[0;32m     17\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mReward = \u001b[39m\u001b[38;5;124m'\u001b[39m \u001b[38;5;241m+\u001b[39m \u001b[38;5;28mstr\u001b[39m(reward))\n",
      "File \u001b[1;32mC:\\Education\\notebooks\\evolve_cart\\run_cart.py:236\u001b[0m, in \u001b[0;36mrun_cart_continuous\u001b[1;34m(agent, simulation_seed, n_episodes, env, max_steps, graphics)\u001b[0m\n\u001b[0;32m    226\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\" Runs the continous cart problem, with the agent mapping observations to actions \u001b[39;00m\n\u001b[0;32m    227\u001b[0m \u001b[38;5;124;03m    - agent: should implement a method act(observation, reward, done)\u001b[39;00m\n\u001b[0;32m    228\u001b[0m \u001b[38;5;124;03m    - simulation_seed: used to set the random seed for simulation\u001b[39;00m\n\u001b[1;32m   (...)\u001b[0m\n\u001b[0;32m    231\u001b[0m \u001b[38;5;124;03m    - graphics: If True, render() will be called.\u001b[39;00m\n\u001b[0;32m    232\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[0;32m    234\u001b[0m \u001b[38;5;66;03m#gym.make('MountainCarContinuous-v0') # CMC() # cc.Continuous_MountainCarEnv()\u001b[39;00m\n\u001b[0;32m    235\u001b[0m \u001b[38;5;66;03m#env.reset() # seed=simulation_seed\u001b[39;00m\n\u001b[1;32m--> 236\u001b[0m \u001b[43menv\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mseed\u001b[49m(simulation_seed)\n\u001b[0;32m    238\u001b[0m reward \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m0\u001b[39m\n\u001b[0;32m    239\u001b[0m cumulative_reward \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m0\u001b[39m\n",
      "\u001b[1;31mAttributeError\u001b[0m: 'CMC_original' object has no attribute 'seed'"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import os\n",
    "import matplotlib\n",
    "import run_cart\n",
    "import gym\n",
    "import numpy as np\n",
    "\n",
    "class random_agent(object):\n",
    "    \"\"\"Random agent\"\"\"\n",
    "\n",
    "    def act(self, observation, reward, done):\n",
    "        return [2.0*np.random.rand()-1.0]\n",
    "    \n",
    "agent = random_agent()\n",
    "reward = run_cart.run_cart_continuous(agent, env=run_cart.CMC_original(), graphics=True)\n",
    "print('Reward = ' + str(reward))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using a continuous time recurrent neural network as controller\n",
    "\n",
    "Of course, random control is not going to solve the task. Below, we introduce an agent that uses a continuous time recurrent neural network for control. This network was introduced by Randall Beer and used for one of the first times in the following article: Beer, Randall D. The dynamics of active categorical perception in an evolved model agent. Adaptive Behavior 11.4 (2003): 209-243.\n",
    "\n",
    "CTRNNs are different from typical artificial neural networks (ANNs) in that they have an activation state $s$ which leads to a 'neural inertia'. Specifically, the activation formula, expressed as a differential equation, is:\n",
    "\n",
    "\\begin{equation}\n",
    "\\tau_i \\dot{s_i} = -s_i + I_i + \\sum_{j=1}^{N}{w_{j,i} \\sigma (g_j(s_j + \\theta_j))},\n",
    "\\end{equation}\n",
    "\n",
    "where $\\tau_i$ is the time constant of neuron $i$, $s_i$ its activation state, and $I_i$ its external inputs. $\\sigma$ is the activation function, $g_j$ the gain of neuron $j$, and $\\theta_j$ its bias. \n",
    "\n",
    "Below, we use <A HREF=\"https://github.com/madvn/CTRNN\" TARGET=\"_blank\">this CTRNN package</A>. It makes no explicit difference between types of neurons. So a network of size 10 means essentially a fully connected network of 10 neurons, where each neuron has connections to all other neurons and itself. In the code below, we give external inputs to the first two neurons, setting these inputs to the observations in the continuous mountain car task. We read out the last neuron as the output to be used for the car control. The activation function is $\\sigma(s) = \\frac{1}{1+e^{-s}}$.\n",
    "\n",
    "<FONT COLOR=\"red\">Exercise 1.</FONT>\n",
    "1. Run the network multiple times. What values does it output over the different runs? Why?\n",
    "2. Suppose that we want to start optimizing the weights and other parameters to have the car achieve the task. The ```act``` function below then has a fundamental problem. Can you spot what the problem is?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from CTRNN import CTRNN\n",
    "from scipy.sparse import csr_matrix\n",
    "\n",
    "class CTRNN_agent(object):\n",
    "    \n",
    "    \"\"\" Continuous Time Recurrent Neural Network agent. \"\"\"\n",
    "    \n",
    "    n_observations = 2;\n",
    "    n_actions = 1;\n",
    "    \n",
    "    def __init__(self, network_size, weights=[], taus = [], gains = [], biases = []):\n",
    "        self.network_size = network_size;\n",
    "        if(self.network_size < self.n_observations + self.n_actions):\n",
    "            self.network_size = self.n_observations + self.n_actions;\n",
    "        self.cns = CTRNN(self.network_size, step_size=0.1) \n",
    "        if(len(weights) > 0):\n",
    "            # weights must be a matrix size: network_size x network_size\n",
    "            self.cns.weights = csr_matrix(weights)\n",
    "        if(len(biases) > 0):\n",
    "            self.cns.biases = biases\n",
    "        if(len(taus) > 0):\n",
    "            self.cns.taus = taus\n",
    "        if(len(gains) > 0):\n",
    "            self.gains = gains\n",
    "    \n",
    "    def act(self, observation, reward, done):\n",
    "        external_inputs = np.asarray([0.0]*self.network_size)\n",
    "        external_inputs[0:self.n_observations] = observation\n",
    "        self.cns.euler_step(external_inputs)\n",
    "        return self.cns.outputs[-self.n_actions:]\n",
    "\n",
    "# set up a CTRNN agent:\n",
    "n_neurons = 10;\n",
    "weights = np.zeros([n_neurons, n_neurons])\n",
    "taus = np.asarray([0.1]*n_neurons)\n",
    "gains = np.ones([n_neurons,])\n",
    "biases = np.zeros([n_neurons,])\n",
    "agent = CTRNN_agent(n_neurons, weights=weights, taus = taus, gains = gains, biases = biases)\n",
    "\n",
    "# run the agent:\n",
    "reward = run_cart.run_cart_continuous(agent, env=run_cart.CMC_original(), simulation_seed=0, graphics=False)\n",
    "\n",
    "# print the reward:\n",
    "print('Reward = ' + str(reward))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The simplest Evolutionary Algorithm for evolving a CTRNN\n",
    "\n",
    "Below we implement a very simple evolutionary algorithm that optimizes all CTRNN parameters for the mountain car task. \n",
    "\n",
    "### Genotype to Phenotype mapping\n",
    "We start by changing the CTRNN_agent class, so that it can be initialized with a _genome_. An advantage of using metaheuristic search methods such as evolutionary algorithms is that they can optimize any type of parameter. For evolutionary robotics this means that even the body of the robot and its sensory layout can be optimized. For the CTRNN agent it means that the genome does not only represent the weights, but all parameters for the network, including the time constants and activation function gains. In this case, we map the genotype (all numbers in the range $[0,1]$) to the right phenotype in the ```__init__``` function. We also implement an evaluation function. Read the code in the cell below and run it to activate the new definitions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "from CTRNN import CTRNN\n",
    "from scipy.sparse import csr_matrix\n",
    "import run_cart\n",
    "import gym\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "# added unpacking of genome:\n",
    "class CTRNN_agent(object):\n",
    "    \n",
    "    \"\"\" Continuous Time Recurrent Neural Network agent. \"\"\"\n",
    "    \n",
    "    n_observations = 2;\n",
    "    n_actions = 1;\n",
    "    \n",
    "    def __init__(self, network_size, genome = [], weights=[], taus = [], gains = [], biases = []):\n",
    "        \n",
    "        self.network_size = network_size;\n",
    "        if(self.network_size < self.n_observations + self.n_actions):\n",
    "            self.network_size = self.n_observations + self.n_actions;\n",
    "        self.cns = CTRNN(self.network_size, step_size=0.1) \n",
    "        \n",
    "        if(len(genome) == self.network_size*self.network_size+3*self.network_size):\n",
    "            # Get the network parameters from the genome:\n",
    "            weight_range = 3\n",
    "            ind = self.network_size*self.network_size\n",
    "            w = weight_range * (2.0 * (genome[:ind] - 0.5))\n",
    "            weights = np.reshape(w, [self.network_size, self.network_size])\n",
    "            biases = weight_range * (2.0 * (genome[ind:ind+self.network_size] - 0.5))\n",
    "            ind += self.network_size\n",
    "            taus = 0.9 * genome[ind:ind+self.network_size] + 0.05\n",
    "            ind += self.network_size\n",
    "            gains = 2.0 * (genome[ind:ind+self.network_size]-0.5)\n",
    "        \n",
    "        if(len(weights) > 0):\n",
    "            # weights must be a matrix size: network_size x network_size\n",
    "            self.cns.weights = csr_matrix(weights)\n",
    "        if(len(biases) > 0):\n",
    "            self.cns.biases = biases\n",
    "        if(len(taus) > 0):\n",
    "            self.cns.taus = taus\n",
    "        if(len(gains) > 0):\n",
    "            self.gains = gains\n",
    "    \n",
    "    def act(self, observation, reward, done):\n",
    "        external_inputs = np.asarray([0.0]*self.network_size)\n",
    "        external_inputs[0:self.n_observations] = observation\n",
    "        self.cns.euler_step(external_inputs)\n",
    "        output = 2.0 * (self.cns.outputs[-self.n_actions:] - 0.5)\n",
    "        return output\n",
    "\n",
    "def evaluate(genome, seed = 0, n_episodes = 1, graphics = False):\n",
    "    # create the phenotype from the genotype:\n",
    "    agent = CTRNN_agent(n_neurons, genome=genome)\n",
    "    # run the agent:\n",
    "    reward = run_cart.run_cart_continuous(agent, env=run_cart.CMC_original(), simulation_seed=seed, n_episodes=n_episodes, graphics=graphics)\n",
    "    #print('Reward = ' + str(reward))\n",
    "    return reward\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The evolutionary algorithm\n",
    "The cell below shows the code for a simple evolutionary algorithm (EA). \n",
    "\n",
    "<FONT COLOR=\"red\">Exercise 2.</FONT>\n",
    "1. How often is the simulator run in the course of the entire evolution?\n",
    "2. The code includes the statement ```np.random.seed(1)```.  What happens if you run the evolution again? What is the reason to include this in the code?\n",
    "3. Run the EA for seeds $\\{1,2,3, \\ldots, 6\\}$. How often does the evolution succeed? What can be the reasons for this?\n",
    "4. The EA bases the new generation only on the best 3 individuals. It is hence very _elitist_. What is the disadvantage of a very elitist EA? And what is the potential advantage? \n",
    "5. Change some EA parameters to see if the evolution finds solutions easier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters CTRNN:\n",
    "network_size = 10\n",
    "genome_size = (network_size+3)*network_size\n",
    "\n",
    "# Evolutionary algorithm:\n",
    "n_individuals = 30\n",
    "n_generations = 30\n",
    "p_mut = 0.05\n",
    "n_best = 3\n",
    "\n",
    "np.random.seed(1)\n",
    "Population = np.random.rand(n_individuals, genome_size)\n",
    "Reward = np.zeros([n_individuals,])\n",
    "max_fitness = np.zeros([n_generations,])\n",
    "mean_fitness = np.zeros([n_generations,])\n",
    "Best = []\n",
    "fitness_best = []\n",
    "for g in range(n_generations):\n",
    "    \n",
    "    # evaluate:\n",
    "    for i in range(n_individuals):\n",
    "        Reward[i] = evaluate(Population[i, :])\n",
    "    mean_fitness[g] = np.mean(Reward)\n",
    "    max_fitness[g] = np.max(Reward)\n",
    "    print('Generation {}, mean = {} max = {}'.format(g, mean_fitness[g], max_fitness[g]))\n",
    "    # select:\n",
    "    inds = np.argsort(Reward)\n",
    "    inds = inds[-n_best:]\n",
    "    if(len(Best) == 0 or Reward[-1] > fitness_best):\n",
    "        Best = Population[inds[-1], :]\n",
    "        fitness_best = Reward[-1]\n",
    "    # vary:\n",
    "    NewPopulation = np.zeros([n_individuals, genome_size])\n",
    "    for i in range(n_individuals):\n",
    "        ind = inds[i % n_best]\n",
    "        NewPopulation[i,:] = Population[ind, :]\n",
    "        for gene in range(genome_size):\n",
    "            if(np.random.rand() <= p_mut):\n",
    "                NewPopulation[i,gene] = np.random.rand()\n",
    "    Population = NewPopulation\n",
    "\n",
    "print('Best fitness ' + str(fitness_best))\n",
    "print('Genome = ')\n",
    "for gene in range(len(Best)):\n",
    "    if(gene == 0):\n",
    "        print('[' + str(Best[gene]) + ', ', end='');\n",
    "    elif(gene == len(Best)-1):\n",
    "        print(str(Best[gene]) + ']');\n",
    "    else:\n",
    "        print(str(Best[gene]) + ', ', end='');\n",
    "\n",
    "\n",
    "plt.figure();\n",
    "plt.plot(range(n_generations), mean_fitness)\n",
    "plt.plot(range(n_generations), max_fitness)\n",
    "plt.xlabel('Generations')\n",
    "plt.ylabel('Fitness')\n",
    "plt.legend(['Mean fitness', 'Max fitness'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "evaluate(Best, graphics=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<FONT COLOR=\"red\">Exercise 3.</FONT>\n",
    "1. The EA runs evaluate without changing the simulation seed. It may hence be that the successful evolution has overfitted on the initial conditions with seed 0. Change the seed in the first cell below and see if it still works.\n",
    "2. Run the more extensive evaluation in the second cell below. What is your conclusion?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "evaluate(Best, seed = 646028073333, graphics=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_tests = 100\n",
    "fit = np.zeros([n_tests,])\n",
    "for t in range(n_tests):\n",
    "    fit[t] = evaluate(Best, seed = 100+t, graphics=False)\n",
    "\n",
    "plt.figure()\n",
    "plt.boxplot(fit)\n",
    "plt.ylabel('Fitness')\n",
    "plt.xticks([1], ['Fitness best individual'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Changing the fitness function\n",
    "The reward function is a bit _deceptive_. Specifically, agents are rewarded 100 points for achieving the goal, but they are also punished each time step for performing control actions. At the start of learning, when agents do not yet reach the goal, this may lead to the agents learning to stand still. Most learning methods learn best when the progress is rewarded gradually more. When thinking about the search space, the ideal case is a convex function that has a gradient everywhere and only a single maximum. The worst case is a 'needle in a haystack' problem, in which the fitness is 0 everywhere and maximal only at a single point in the search space.\n",
    "\n",
    "Given the low evolvability of the problem setting above, we now change the reward function. In the cell below, we introduce a class CMC (Continuous Mountain Car) that extends the Continuous_MountainCarEnv class. It reimplements the init, reset and step methods to: (i) keep track of the minimal distance until now to the goal position, and (ii) at each time step add ```1. - self.min_distance / self.max_distance``` to the reward. This should favor individuals early on in evolution that are able to go higher up the slope.\n",
    "\n",
    "In the evaluate and test_best functions, we add a parameter ```original_reward``` that can switch between the original reward definition and this new one.\n",
    "\n",
    "Run the cell below to load these new definitions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gym.envs.classic_control as cc\n",
    "import math\n",
    "\n",
    "class CMC(cc.Continuous_MountainCarEnv):\n",
    "    \"\"\" Derived class of Continuous Mountain Car, so that we can change, e.g., the reward function.\n",
    "    \"\"\"\n",
    "    # Based on: https://raw.githubusercontent.com/openai/gym/master/gym/envs/classic_control/continuous_mountain_car.py\n",
    "    \n",
    "    n_steps_per_render = 5\n",
    "    \n",
    "    def __init__(self):\n",
    "        self.render_step = 0\n",
    "        self.figure_handle = []\n",
    "        super(CMC, self).__init__()\n",
    "        self.max_distance = self.max_position - self.min_position\n",
    "        self.min_distance = self.max_distance\n",
    "    \n",
    "    def reset(self):\n",
    "        super(CMC, self).reset()\n",
    "        self.max_distance = self.max_position - self.min_position\n",
    "        self.min_distance = self.max_distance\n",
    "        self.render_step = 0\n",
    "        if(self.figure_handle != []):\n",
    "            plt.close('mountain_car')\n",
    "            self.figure_handle = []\n",
    "        return np.array(self.state)\n",
    "    \n",
    "    def step(self, action):\n",
    "\n",
    "        position = self.state[0]\n",
    "        velocity = self.state[1]\n",
    "        force = min(max(action[0], -1.0), 1.0)\n",
    "\n",
    "        velocity += force*self.power -0.0025 * math.cos(3*position)\n",
    "        if (velocity > self.max_speed): velocity = self.max_speed\n",
    "        if (velocity < -self.max_speed): velocity = -self.max_speed\n",
    "        position += velocity\n",
    "        if (position > self.max_position): position = self.max_position\n",
    "        if (position < self.min_position): position = self.min_position\n",
    "        if (position==self.min_position and velocity<0): velocity = 0\n",
    "\n",
    "        done = bool(position >= self.goal_position and velocity >= self.goal_velocity)\n",
    "\n",
    "        # Now you can change the reward function here:\n",
    "        distance = abs(position - self.goal_position)\n",
    "        if(distance < self.min_distance):\n",
    "            self.min_distance = distance\n",
    "            \n",
    "        reward = 0\n",
    "        if done:\n",
    "            reward = 100.0\n",
    "        reward -= math.pow(action[0],2)*0.1\n",
    "        reward += 1. - self.min_distance / self.max_distance\n",
    "\n",
    "        self.state = np.array([position, velocity])\n",
    "        return self.state, reward, done, {}\n",
    "    \n",
    "    def render(self, mode='human', sleep_time=0.033):\n",
    "        \n",
    "        if os.name == 'nt':\n",
    "            normal_display = True\n",
    "        else:\n",
    "            if 'DISPLAY' in os.environ.keys():\n",
    "                normal_display = True\n",
    "            else:\n",
    "                normal_display = False\n",
    "            \n",
    "        if normal_display:\n",
    "            super(CMC, self).render()\n",
    "        else:\n",
    "            \n",
    "            self.render_step += 1\n",
    "            \n",
    "            if(self.render_step % self.n_steps_per_render == 0):\n",
    "                # first plot the landscape:\n",
    "                step = 0.01\n",
    "                x_coords = np.arange(self.min_position, self.max_position, step)\n",
    "                y_coords = self._height(x_coords)\n",
    "\n",
    "                if(self.figure_handle == []):\n",
    "                    self.figure_handle = plt.figure('mountain_car')\n",
    "                    self.ax = self.figure_handle.add_subplot(111)\n",
    "                    plt.ion()\n",
    "                    #self.figure_handle.show()\n",
    "                    self.figure_handle.canvas.draw()\n",
    "                else:\n",
    "                    plt.figure('mountain_car')\n",
    "\n",
    "                self.ax.clear()\n",
    "                self.ax.plot(x_coords, y_coords)\n",
    "                self.ax.plot(self.state[0], self._height(self.state[0]), 'ro')\n",
    "                self.ax.text(self.goal_position, self._height(self.goal_position)+0.02, 'Goal')        \n",
    "                #        self.figure_handle.canvas.draw()\n",
    "                #        self.figure_handle.show()\n",
    "                display.clear_output(wait=True)\n",
    "                display.display(plt.gcf())\n",
    "                time.sleep(sleep_time)\n",
    "    \n",
    "def evaluate(genome, seed = 0, graphics = False, original_reward=True):\n",
    "    # create the phenotype from the genotype:\n",
    "    agent = CTRNN_agent(n_neurons, genome=genome)\n",
    "    # run the agent:\n",
    "    if(original_reward):\n",
    "        reward = run_cart.run_cart_continuous(agent, env=run_cart.CMC_original(), simulation_seed=seed, graphics=graphics)\n",
    "    else:\n",
    "        reward = run_cart.run_cart_continuous(agent, env=CMC(), simulation_seed=seed, graphics=graphics)\n",
    "    #print('Reward = ' + str(reward))\n",
    "    return reward\n",
    "\n",
    "def test_best(Best, original_reward=True):    \n",
    "    n_tests = 30\n",
    "    fit = np.zeros([n_tests,])\n",
    "    for t in range(n_tests):\n",
    "        fit[t] = evaluate(Best, seed = 100+t, graphics=False, original_reward=True)\n",
    "    \n",
    "    plt.figure()\n",
    "    plt.boxplot(fit)\n",
    "    plt.ylabel('Fitness')\n",
    "    plt.xticks([1], ['Fitness best individual'])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<FONT COLOR=\"red\">Exercise 4.</FONT>\n",
    "1. Run the EA for seeds $\\{1,2,3, \\ldots, 6\\}$. How often does the evolution succeed?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters CTRNN:\n",
    "network_size = 10\n",
    "genome_size = (network_size+3)*network_size\n",
    "\n",
    "# Evolutionary algorithm:\n",
    "n_individuals = 30\n",
    "n_generations = 30\n",
    "p_mut = 0.05\n",
    "n_best = 3\n",
    "\n",
    "np.random.seed(1) \n",
    "original_reward = False\n",
    "Population = np.random.rand(n_individuals, genome_size)\n",
    "Reward = np.zeros([n_individuals,])\n",
    "max_fitness = np.zeros([n_generations,])\n",
    "mean_fitness = np.zeros([n_generations,])\n",
    "Best = []\n",
    "fitness_best = []\n",
    "for g in range(n_generations):\n",
    "    \n",
    "    # evaluate:\n",
    "    for i in range(n_individuals):\n",
    "        Reward[i] = evaluate(Population[i, :], original_reward=original_reward)\n",
    "    mean_fitness[g] = np.mean(Reward)\n",
    "    max_fitness[g] = np.max(Reward)\n",
    "    print('Generation {}, mean = {} max = {}'.format(g, mean_fitness[g], max_fitness[g]))\n",
    "    # select:\n",
    "    inds = np.argsort(Reward)\n",
    "    inds = inds[-n_best:]\n",
    "    if(len(Best) == 0 or Reward[-1] > fitness_best):\n",
    "        Best = Population[inds[-1], :]\n",
    "        fitness_best = Reward[-1]\n",
    "    # vary:\n",
    "    NewPopulation = np.zeros([n_individuals, genome_size])\n",
    "    for i in range(n_individuals):\n",
    "        ind = inds[i % n_best]\n",
    "        NewPopulation[i,:] = Population[ind, :]\n",
    "        for gene in range(genome_size):\n",
    "            if(np.random.rand() <= p_mut):\n",
    "                NewPopulation[i,gene] = np.random.rand()\n",
    "    Population = NewPopulation\n",
    "\n",
    "print('Best fitness ' + str(fitness_best))\n",
    "print('Genome = ')\n",
    "for gene in range(len(Best)):\n",
    "    if(gene == 0):\n",
    "        print('[' + str(Best[gene]) + ', ', end='');\n",
    "    elif(gene == len(Best)-1):\n",
    "        print(str(Best[gene]) + ']');\n",
    "    else:\n",
    "        print(str(Best[gene]) + ', ', end='');\n",
    "\n",
    "plt.figure();\n",
    "plt.plot(range(n_generations), mean_fitness)\n",
    "plt.plot(range(n_generations), max_fitness)\n",
    "plt.xlabel('Generations')\n",
    "plt.ylabel('Fitness')\n",
    "plt.legend(['Mean fitness', 'Max fitness'])\n",
    "\n",
    "evaluate(Best, graphics=True, original_reward=original_reward)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<FONT COLOR=\"red\">Exercise 5.</FONT>\n",
    "1. Try out different seeds for a visual evaluation of the performance in the 1st cell below. What do you observe?\n",
    "2. Run the ```test_best``` function in the 2nd cell below. What do you observe? \n",
    "3. What two conclusions can you draw from this concerning the new fitness function?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "evaluate(Best, seed=5, graphics=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_best(Best, original_reward=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analyzing the robot's behavior\n",
    "\n",
    "This notebook has focused on setting up a simple evolutionary algorithm for solving the continuous mountain car task in the AI gym, highlighting some of the challenges in obtaining a successful controller.\n",
    "\n",
    "When evolution finds a solution to the task, it is often not clear to the designer _how_ the controller has solved the task. Often, evolution finds solutions that exploit the interaction of the robot with the environment. Specifically, evolved controllers exploit sensory-motor coordination, in which the closed loop of motor actions and sensory observations is exploited allowing to solve tasks with less processing. \n",
    "\n",
    "Analyzing how an evolved neural controller works is useful for understanding how well the controller will generalize to different instantiations of the task. Also, it gives insight into the nature of intelligence.  \n",
    "\n",
    "Currently, analysis of the evolved controller is beyond the scope of this notebook. However, we recommend references [2] and [3] for getting ideas on the types of analysis that can be applied."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "[1] Doncieux, S., Bredeche, N., Mouret, J. B., & Eiben, A. E. G. (2015). Evolutionary robotics: what, why, and where to. Frontiers in Robotics and AI, 2, 4. <A HREF=\"https://www.frontiersin.org/articles/10.3389/frobt.2015.00004/full\" TARGET=\"_blank\">(open access)</A>\n",
    "\n",
    "[2] Nolfi, Stefano. Power and the limits of reactive agents. Neurocomputing 42, 1 (2002): 119-145.\n",
    "\n",
    "[3] Beer, Randall D. The dynamics of active categorical perception in an evolved model agent. Adaptive Behavior 11.4 (2003): 209-243."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Answers\n",
    "\n",
    "<FONT COLOR=\"red\">Exercise 1.</FONT>\n",
    "\n",
    "1. Run the network multiple times. What values does it output over the different runs? Why?\n",
    "\n",
    "It always outputs ```Reward = -24.99999999999965```. The reason is that the network is always the same and ```simulation_seed``` is always 0. This seed is used for random number generation in the simulation, and the same number will hence always result in the same initial conditions. \n",
    "\n",
    "2. Suppose that we want to start optimizing the weights and other parameters to have the car achieve the task. The ```act``` function below then has a fundamental problem. Can you spot what the problem is?\n",
    "\n",
    "The ```act``` function returns an output of the neural network as the action. The action for the continuous mountain car task can be in the range $[-1,1]$. However, the output of the neurons in the CTRNNs is determined with a sigmoid function: $\\sigma(s) = \\frac{1}{1+e^{-s}}$, which is in the range $[0,1]$. This means the car would only be able to accelerate to the right.\n",
    "\n",
    "\n",
    "<FONT COLOR=\"red\">Exercise 2.</FONT>\n",
    "1. How often is the simulator run in the course of the entire evolution?\n",
    "\n",
    "The number of simulator runs is: #generations $\\times$ #individuals $\\times$ #episodes, which is in this case $30 \\times 30 \\times 1 = 900$.\n",
    "\n",
    "2. The code includes the statement ```np.random.seed(1)```.  What happens if you run the evolution again? What is the reason to include this in the code?\n",
    "\n",
    "This random seed determines the initial population and subsequent mutations. So running the evolution again will give the same results. This is useful for reproducibility of the results.\n",
    "\n",
    "3. Run the EA for seeds $\\{1,2,3, \\ldots, 6\\}$. How often does the evolution succeed? What can be the reasons for this?\n",
    "\n",
    "The evolution only succeeds for the random seed 6, so once out of 7 times. The current setup has a low _evolvability_. There can be many reasons for the EA failing to find a good solution:\n",
    "* The EA does not search the search space well enough. \n",
    "* The fitness function is not suitable for the way in which the EA searches.\n",
    "* The genome cannot encode for the solution of the task, i.e., there is no solution in the search space. In this case, since we have at least one successful evolution we can exclude this option, but up until seed 5 this could have been a plausible explanation. A less binary version of this statement is that there are few solutions in the (high-dimensional) search space and the genome is ill-suited for encoding the solution.\n",
    "\n",
    "4. The EA bases the new generation only on the best 3 individuals. It is hence very _elitist_. What is the disadvantage of a very elitist EA? And what is the potential advantage?  \n",
    "\n",
    "The disadvantage is that there will be very little genetic diversity in the population. The search can easily get stuck in a local optimum. The potential advantage is that one needs few generations to get a good solution (in this case only 30). Setting the number of selected individuals higher will lead to a slower convergence to a solution, costing more simulation and hence evolution time.\n",
    "\n",
    "5. Change some EA parameters to see if the evolution finds solutions easier.\n",
    "\n",
    "This is for the students to try out. Let us know when you find settings that result in substantially better evolvability.\n",
    "\n",
    "<FONT COLOR=\"red\">Exercise 3.</FONT>\n",
    "1. The EA runs evaluate without changing the simulation seed. It may hence be that the successful evolution has overfitted on the initial conditions with seed 0. Change the seed in the first cell below and see if it still works.\n",
    "\n",
    "It still works for most seeds. \n",
    "\n",
    "2. Run the more extensive evaluation in the second cell below. What is your conclusion?\n",
    "\n",
    "The error plot shows a very high performance (~98.1 on average) with little variance. However, looking at the graphical rendering in the previous question shows that the car always starts in similar conditions, despite a changing random seed. \n",
    "\n",
    "<FONT COLOR=\"red\">Exercise 4.</FONT>\n",
    "1. Run the EA for seeds $\\{1,2,3, \\ldots, 6\\}$. How often does the evolution succeed?\n",
    "\n",
    "Most of the seeds now lead to a successful agent.\n",
    "\n",
    "<FONT COLOR=\"red\">Exercise 5.</FONT>\n",
    "1. Try out different seeds for a visual evaluation of the performance in the 1st cell below. What do you observe?\n",
    "\n",
    "Some of the seeds now do not lead to the car reaching the flag. The agent has a tendency to mostly focus on accelerating to the right.\n",
    "\n",
    "2. Run the ```test_best``` function in the 2nd cell below. What do you observe? \n",
    "\n",
    "The bulk of the results indicates that the car does not reach the flag. Only in a few cases does the car reach the flag. Then still, the original reward is rather low compared with the 98 before.\n",
    "\n",
    "3. What two conclusions can you draw from this concerning the new fitness function?\n",
    "\n",
    "We can conclude that (i) the found solution generalizes less well to different initial conditions, and (ii) the new fitness function leads to strategies that use higher control inputs, and hence a lower fitness with the original reward function.\n",
    "\n",
    "Finding the right fitness / reward function is called _reward shaping_, and it is currently a major problem in both evolutionary robotics and reinforcement learning."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}