5_strategies/evolutionary_strategies.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<img src=\"https://github.com/d9w/evolution/raw/master/imgs/logo.png\" width=\"20%\" align=\"right\" style=\"margin:0px 20px\">\n",
    "\n",
    "\n",
    "# Evolutionary Algorithms\n",
    "\n",
    "## Evolutionary Strategies\n",
    "\n",
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by-sa/4.0/\"><img alt=\"Creative Commons License\" align=\"left\" src=\"https://i.creativecommons.org/l/by-sa/4.0/80x15.png\" /></a>&nbsp;| Dennis G. Wilson | <a href=\"https://d9w.github.io/evolution/\">https://d9w.github.io/evolution/</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import ticker, cm\n",
    "import matplotlib.colors as colors\n",
    "import numpy as np\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# A simple $(1 + \\lambda)$ ES"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The simplest Evolutionary Strategy is the $(1+1)$ ES. In this algorithm, the best individual generated is replaced by a new individual when its **fitness** is superior to the current best. For continuous optimization, we can use a representation of floating point numbers which are modified at each step using a Normal distribution:\n",
    "\n",
    "    Initialize x randomly in ℝ\n",
    "    while not terminate\n",
    "        x' = x + 𝑁(0, σ)\n",
    "        if f(x′) < f(x)\n",
    "            x = x'\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Evolutionary Strategies are commonly used for continuous optimization. We will use standard [test functions for continuous optimization](https://en.wikipedia.org/wiki/Test_functions_for_optimization), namely the Himmelblau and Rosenbrock functions. They are defined as follows:\n",
    "\n",
    "Himmelblau: $f(x,y) = (x^2+y-11)^2 + (x+y^2-7)^2$\n",
    "\n",
    "Rosenbrock: $f(x) = \\sum_{i=1}^{n-1}(B(x_{i+1}-x_i^2)^2+(a-x_i)^2)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def himmelblau(x, y):\n",
    "    return (x**2 + y - 11)**2 + (x + y**2 - 7)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def rosenbrock(x, y, a=1, B=100):\n",
    "      return (a-x)**2 + B*((y-x**2))**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We can visualize these functions for a better understanding of the search space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "X = np.arange(-5, 5, 0.1)\n",
    "Y = np.arange(-5, 5, 0.1)\n",
    "X, Y = np.meshgrid(X, Y)\n",
    "Z = himmelblau(X, Y)\n",
    "#Z = rosenbrock(X, Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8, 6))\n",
    "cs = plt.contour(X, Y, Z, levels=100, cmap='Spectral',\n",
    "                 norm=colors.Normalize(vmin=Z.min(), vmax=Z.max()), alpha=0.4)\n",
    "fig.colorbar(cs)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "For the (1+1) ES, we start with a random point $x$ and then sample one point around it using a Normal distribution.|"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "x = np.random.rand(2) * 8 - 4\n",
    "x_orig = np.copy(x)\n",
    "print(\"position: \", x, \"fitness :\", himmelblau(x[0], x[1]))\n",
    "σ = 0.3\n",
    "x_t = np.random.randn(2)*σ + x\n",
    "print(\"new position: \", x_t, \"new fitness :\", himmelblau(x_t[0], x_t[1]))\n",
    "fig = plt.figure(figsize=(8, 6))\n",
    "cs = plt.contour(X, Y, Z, levels=100, cmap='Spectral', norm=colors.Normalize(vmin=Z.min(), vmax=Z.max()), alpha=0.4)\n",
    "plt.annotate('x\\'', (x_t))\n",
    "plt.scatter(x_t[0], x_t[1], c='b')\n",
    "plt.annotate('x', (x))\n",
    "plt.scatter(x[0], x[1], c='r')\n",
    "fig.colorbar(cs);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "While this greedy algorithm could allow us to move around the search space and eventually find an optimum, we can use multiple samples of this distribution to better inform the movement in the search space. We therefore instead sample multiple points following a Normal distribution around $x$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "x = np.copy(x_orig)\n",
    "print(\"position: \", x, \"fitness :\", himmelblau(x[0], x[1]))\n",
    "x_t = np.array([x + np.random.normal(size=(2,))*σ for i in range(20)])\n",
    "x_best = x\n",
    "fit_best = himmelblau(x[0], x[1])\n",
    "for i in range(20):\n",
    "    f = himmelblau(x_t[i, 0], x_t[i, 1])\n",
    "    if f < fit_best:\n",
    "        fit_best = f\n",
    "        x_best = x_t[i, 0], x_t[i, 1]\n",
    "print(\"best position: \", x_best, \"best fitness :\", fit_best)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8, 6))\n",
    "cs = plt.contour(X, Y, Z, levels=100, cmap='Spectral', norm=colors.Normalize(vmin=Z.min(), vmax=Z.max()), alpha=0.4)\n",
    "plt.scatter(x_t[:, 0], x_t[:, 1], c='b')\n",
    "plt.scatter(x[0], x[1], c='r')\n",
    "plt.scatter(x_best[0], x_best[1], c='g')\n",
    "fig.colorbar(cs);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Having this population of points will enable us to explore the space around $x$ before moving on to the next step. However, these points might overlap or cover spaces that have already been explored. Generating too many points could slow down search, requiring more evaluations of the objective function. In many applications, that is costly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Performing random optimization but sampling more than 1 point leads to the $(1+\\lambda)$ Evolutionary Strategy. \n",
    "\n",
    "    Initialize x randomly in ℝ\n",
    "    while not terminate\n",
    "        x_p = x\n",
    "        for i in [1,λ]\n",
    "            x_i = x_p + 𝑁(0, σ)\n",
    "            if f(x_i) < f(x)\n",
    "                x = x_i\n",
    "        x_p = x\n",
    "    return x_p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In the scope of evolutionary algorithms, we'll refer to $x$ as a parent and all $x'$ points as offspring. Each point can also be referred to as an individual in a population, and $f(x)$ is called the fitness of the individual. Each iteration of the algorithm is called a generation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "The $(\\mu/\\rho,\\lambda)$ or $(\\mu/\\rho+\\lambda)$ notation signifies the configuration of parents and offspring. In this notation, $\\mu$ is the number of parents, $\\rho$ is the number of parents involved in creating the offspring, and $\\lambda$ is the number of offspring. $(\\mu/\\rho+\\lambda)$ means the parents can be included in the next population, whereas $(\\mu/\\rho,\\lambda)$ means the parents are not included. $(1+\\lambda)$ therefore means that 1 parent is involved in creating a population of $\\lambda$ offspring and can be included in the next generation (ie, if $f(x) < f(x') \\forall x'$, $x$ does not change).\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def oneplus_lambda(x, fitness, gens=100, lam=20):\n",
    "    x_best = x\n",
    "    f_best = fitness(x)\n",
    "    fits = np.zeros(gens)\n",
    "    for g in range(gens):\n",
    "        N = np.random.normal(size=(lam, len(x)))\n",
    "        for i in range(lam):\n",
    "            ind = x + N[i, :]\n",
    "            f = fitness(ind)\n",
    "            if f < f_best:\n",
    "                f_best = f\n",
    "                x_best = ind\n",
    "        x = x_best\n",
    "        fits[g] = f_best\n",
    "    return fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's run the $(1+\\lambda)$ ES on the himmelblau function and record the evolution over multiple trials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x : himmelblau(x[0], x[1])\n",
    "runs = []\n",
    "for i in range(5):\n",
    "    x = np.random.rand(2)*10-5\n",
    "    runs += [oneplus_lambda(x, f)]\n",
    "runs = np.array(runs)\n",
    "runs_mean = np.mean(runs, 0)\n",
    "runs_std = np.std(runs, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8, 6))\n",
    "plt.fill_between(range(len(runs[0])), runs_mean+0.5*runs_std, runs_mean-0.5*runs_std, alpha=0.5)\n",
    "plt.plot(range(len(runs[0])), runs_mean)\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <h3>Exercise 1</h3>\n",
    "    Study the impact of the $\\lambda$ parameter by modifying it and re-running the optimization. What is the best $\\lambda$ value for the Himmelblau problem?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <h3>Exercise 2</h3>\n",
    "    Change to the Rosenbrock problem. What is the difficulty of this search space, compared to Himmelblau? Can you find a $\\lambda$ parameter adapted to this search space?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Approximating the gradient"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "In the $(1+\\lambda)$ ES, we move from the best point in a population to the best point in the next randomly sampled population. However, with the population information, we can do better than that; we can move in the direction of the gradient approximated by the fitness values of the population. This can also better use the full population information, combining all fitness information instead of selecting the best fitness value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "x = np.copy(x_orig)\n",
    "x_t = np.array([x + np.random.normal(size=(2,))*0.5 for i in range(10)])\n",
    "fits = himmelblau(x_t[:, 0], x_t[:, 1])\n",
    "fig = plt.figure(figsize=(8, 6))\n",
    "norm = colors.Normalize(vmin=Z.min(), vmax=Z.max())\n",
    "cs = plt.contour(X, Y, Z, levels=100, cmap='Spectral', alpha=0.4, norm=norm)\n",
    "plt.quiver(x_t[:, 0], x_t[:, 1], x_t[:, 0]-x[0], x_t[:, 1]-x[1],\n",
    "           fits, scale=fits/60, scale_units='xy', cmap='Spectral', norm=norm)\n",
    "plt.scatter(x[0], x[1], c='k')\n",
    "fig.colorbar(cs)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's consider our objective function $f$ and the distribution of individuals $\\pi(x|\\theta)$ for a specific center $\\theta$. The goal is to maximize the expected fitness under the search distribution, in other words:\n",
    "\n",
    "$J(\\theta) = E_\\theta[f(x)] = \\int f(x)\\pi(x|\\theta)dx$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "This can be rewritten as:\n",
    "\n",
    "$E_\\theta[f(x)\\nabla_\\theta \\log \\pi(x|\\theta)]$\n",
    "\n",
    "which can be approximated as\n",
    "\n",
    "$\\nabla_\\theta J(\\theta) ~= \\frac{1}{\\lambda}\\sum_{k=1}^\\lambda f(x_k) \\nabla_\\theta log(\\pi(x_k|\\theta)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Instead of directly using the fitness values $f$, we will normalize all fitness values by the average $\\mu(f(x))$ and standard deviation $\\sigma(f(x))$ of the population.\n",
    "\n",
    "$A = \\frac{f(x) - \\mu(f(x))}{\\sigma(f(x))}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We'll then define a vector at each point, $x_i - x$. Note that since $x_i = x + N(0, 1)$, we could store just this original random sampled point from $N(0, 1)$ as our vector.\n",
    "\n",
    "$N_i = x_i - x$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We multiply $A$ by all individual vectors which scales the vector magnitude by the relative fitness. Finally, we sum all vectors and divide by $\\lambda$ to get the weighted average.\n",
    "\n",
    "$\\frac{1}{\\lambda}\\sum_i A_i N_i$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "This is equivalent to the dot product and is used as our gradient approximation. Since we're minimizing, we'll use the negative value to do gradient descent.\n",
    "\n",
    "$\\nabla f \\approx -\\frac{A \\cdot N}{\\lambda}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "A = (fits - np.mean(fits)) / np.std(fits)\n",
    "N = x_t - x\n",
    "G = -np.dot(A, N) / 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8, 6))\n",
    "norm = colors.Normalize(vmin=Z.min(), vmax=Z.max())\n",
    "cs = plt.contour(X, Y, Z, levels=100, cmap='Spectral', alpha=0.4, norm=norm)\n",
    "plt.quiver(x_t[:, 0], x_t[:, 1], x_t[:, 0]-x[0], x_t[:, 1]-x[1],\n",
    "           fits, scale=fits/100, scale_units='xy', cmap='Spectral', norm=norm)\n",
    "plt.quiver(x[0], x[1], G[0], G[1])\n",
    "plt.scatter(x[0], x[1], c='k')\n",
    "fig.colorbar(cs)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We have our approximate gradient direction, but how far should we move in it? We'll define a learning rate variable $\\alpha$ as the fixed magnitude of movement for now, so our update of $x$ is then\n",
    "\n",
    "$x = x - \\alpha\\frac{A \\cdot N}{\\lambda}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's put this all together as an Evolutionary Strategy. It is a $(\\mu, \\lambda)$ ES as multiple individuals $\\mu$ inform the next generation; however, these individuals are not kept for the next generation.\n",
    "\n",
    "    Initialize x randomly in ℝ\n",
    "    while not terminate\n",
    "        for i in [1,λ]\n",
    "            N_i = 𝑁(0, 1)\n",
    "            F_i = f(x + N_i)\n",
    "        A = (F−𝜇(F))/𝜎(F)\n",
    "        x = x - 𝛼(A⋅N)/𝜆"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def mu_lambda(x, fitness, gens=200, lam=10, alpha=0.2, verbose=False):\n",
    "    x_best = x\n",
    "    f_best = fitness(x)\n",
    "    fits = np.zeros(gens)\n",
    "    for g in range(gens):\n",
    "        N = np.random.normal(size=(lam, len(x)))\n",
    "        F = np.zeros(lam)\n",
    "        for i in range(lam):\n",
    "            ind = x + N[i, :]\n",
    "            F[i] = fitness(ind)\n",
    "            if F[i] < f_best:\n",
    "                f_best = F[i]\n",
    "                x_best = ind\n",
    "                if verbose:\n",
    "                    print(g, \" \", f_best)\n",
    "        fits[g] = f_best\n",
    "        mu_f = np.mean(F)\n",
    "        std_f = np.std(F)\n",
    "        A = F\n",
    "        if std_f != 0:\n",
    "            A = (F - mu_f) / std_f\n",
    "        x = x - alpha * np.dot(A, N) / lam\n",
    "    return fits, x_best"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "runs = []\n",
    "for i in range(5):\n",
    "    x = np.random.randn(2)\n",
    "    fits, _ = mu_lambda(x, f)\n",
    "    runs += [fits]\n",
    "runs = np.array(runs)\n",
    "runs_mean = np.mean(runs, 0)\n",
    "runs_std = np.std(runs, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8, 6))\n",
    "plt.fill_between(range(len(runs[0])), runs_mean+0.5*runs_std, runs_mean-0.5*runs_std, alpha=0.5)\n",
    "plt.plot(range(len(runs[0])), runs_mean)\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <h3>Exercise 3</h3>\n",
    "    \n",
    "Study the impact of the $\\alpha$ parameter by modifying it and re-running the optimization. Compare the best parameters found with the results from the $(1+\\lambda)$ ES.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <h3>Exercise 4</h3>\n",
    "    \n",
    "Study the performance of the $(\\mu, \\lambda)$ ES on the Rastrigin function. Does it perform better or worse than the  $(1+\\lambda)$ ES?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Rank-based update"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We can make two observations about our current $(\\mu, \\lambda)$ Evolutionary Strategy:\n",
    "\n",
    "1. Not all $\\lambda$ points may be useful for the centroid update\n",
    "2. Exact gradient information may be noisier than relative *rank* of individual fitness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's illustrate this point in the Himmelblau function. Instead of using all $\\lambda$ points, we'll use a smaller value of $\\mu$, for example $\\frac{1}{2}\\lambda$. We'll also create a list of weights which exponentially decreases from the first to last weight."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "mu = 5\n",
    "w = np.log(mu + 1/2) - np.log(np.arange(1, mu+1))\n",
    "w /= np.sum(w)\n",
    "w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "As before, we create a random starting point and then add Gaussian noise to create offspring."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "x = np.copy(x_orig)\n",
    "N = np.random.normal(size=(10,2))*0.5\n",
    "x_t = x + N\n",
    "fits = himmelblau(x_t[:, 0], x_t[:, 1])\n",
    "fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Now, however, instead of directly using fitness values, we'll select the best individuals and update using their weighted values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "sorted_ids = np.argsort(fits)\n",
    "weighted = N[sorted_ids[:mu]] * w.reshape(mu, 1)\n",
    "weighted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8, 6))\n",
    "norm = colors.Normalize(vmin=Z.min(), vmax=Z.max())\n",
    "cs = plt.contour(X, Y, Z, levels=100, cmap='Spectral', alpha=0.4, norm=norm)\n",
    "plt.quiver(x_t[:, 0], x_t[:, 1], x_t[:, 0]-x[0], x_t[:, 1]-x[1],\n",
    "           fits, scale=fits/100, scale_units='xy', cmap='Spectral', norm=norm, alpha=0.4)\n",
    "plt.quiver(x[0]+weighted[:, 0], x[1]+weighted[:, 1], weighted[:, 0], weighted[:, 1],\n",
    "           w, scale_units='xy', cmap='Spectral', norm=norm)\n",
    "plt.scatter(x[0], x[1], c='k')\n",
    "fig.colorbar(cs)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The first point is clear: using all points in the centroid update could pull it away from nearby optima. The second point is a bit harder to understand, but motivation is this: any transformation to the search space which **preserves rank** for the top $\\mu$ individuals does not change the update. This invariance property is helpful to search because it increases the predictive power of the fitness by inducing problem equivalence classes.\n",
    "\n",
    "Hansen, N., Ros, R., Mauny, N., Schoenauer, M., & Auger, A. (2011). Impacts of invariance in search: When CMA-ES and PSO face ill-conditioned and non-separable problems. Applied Soft Computing, 11(8), 5755-5769."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Incorporating these two changes to the previous $(\\mu, \\lambda)$ ES brings us an ES referred to as Canonical ES:\n",
    "\n",
    "<img src=\"https://github.com/d9w/evolution/raw/master/imgs/canonical.png\" width=\"40%\" height=\"auto\">\n",
    "\n",
    "Chrabaszcz, P., Loshchilov, I., & Hutter, F. (2018, July). Back to basics: benchmarking canonical evolution strategies for playing Atari. In Proceedings of the 27th International Joint Conference on Artificial Intelligence (pp. 1419-1426). [code](https://github.com/PatrykChrabaszcz/Canonical_ES_Atari)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<img src=\"https://github.com/d9w/evolution/raw/master/imgs/canonical_results.png\" width=\"80%\" height=\"auto\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <h3>Exercise 5</h3>\n",
    "    \n",
    "Consider the Rosenbrock function. How do you expect the gradient rankings to compare to the gradient information? Will they be very similar or different? Plot a sample to investigate.\n",
    "</div>"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}