Optimization-Visualization-mmorais2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualization of Optimization Algorithms\n",
    "\n",
    "## CS519: Scientific Visualization - Fall 2020\n",
    "\n",
    "### netid: <mark>mmorais2</mark>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import display, HTML"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Table of Contents\n",
    "1. Abstract\n",
    "1. Introduction\n",
    "1. Test Functions\n",
    "    1. Rosenbrock Function\n",
    "    2. Goldstein-Price Function\n",
    "    3. Bartels-Conn Function\n",
    "    4. Egg Crate Function\n",
    "1. Optimization Algorithms\n",
    "    1. Gradient Descent\n",
    "        1. Gradient Descent: Rosenbrock\n",
    "        1. Gradient Descent: Goldstein-Price\n",
    "    2. BFGS\n",
    "        1. BFGS: Rosenbrock\n",
    "        1. BFGS: Goldstein-Price\n",
    "    3. Simulated Annealing\n",
    "        1. Simulated Annealing: Rosenbrock\n",
    "        1. Simulated Annealing: Goldstein-Price\n",
    "        1. Simulated Annealing: Bartels-Conn\n",
    "        1. Simulated Annealing: Egg Crate\n",
    "    4. Particle Swarm\n",
    "        1. Particle Swarm: Rosenbrock\n",
    "        1. Particle Swarm: Goldstein-Price\n",
    "        1. Particle Swarm: Bartels-Conn\n",
    "        1. Particle Swarm: Egg Crate\n",
    "1. Results\n",
    "1. Discussion\n",
    "1. References\n",
    "1. Appendix: Links to Animations\n",
    "1. Appendix: Source Code\n",
    "    1. Library Imports\n",
    "    1. Library Versions\n",
    "    1. Surface Generation\n",
    "    1. Rosenbrock Function\n",
    "    1. Goldstein-Price Function\n",
    "    1. Bartels-Conn Function\n",
    "    1. Egg Crate Function\n",
    "    1. Gradient Descent\n",
    "    1. BFGS\n",
    "    1. Simulated Annealing\n",
    "    1. Particle Swarm\n",
    "    1. Static 2d simulation result plot\n",
    "    1. Static 3d simulation result plot\n",
    "    1. Animated 2d simulation result plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Abstract"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Optimization algorithms seek to find the best solution $x^{*}$ from a set S such that $f(x^{*}) \\leq f(x)$ for all $x$ in S. For this project we describe and implement a handful of optimization algorithms, evaluate their performance on some well known test functions, and create visualizations to build some intuition and verify their function.  Finally, we perform a comparative analysis between algorithms using repeated trials on these test functions in order to draw broader conclusions. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given a fuction $f: \\mathbb{R}^n \\rightarrow \\mathbb{R}$ and the set $S \\in \\mathbb{R}^n$ the solution to an optimization problem seeks to find $x^{*} \\in S$ such that $f(x^{*}) \\leq f(x)$ for all $x$ in S.  The function $f$ is referred to as the objective function and the set $S$ is referred to the as the constraints.\n",
    "\n",
    "In this project we will examine 4 different algorithms for solving optimization problems.  Each algorithm represents a more general class of problem-solving approaches to optimzation.\n",
    "\n",
    "In order to understand how each algorithm functions as well as appreciate the strengths and weakness, we choose a set of test functions of varying difficulty.  These functions are described in the next section. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The table below summarizes the test functions used in the project.  Each test function is classified according to criteria given for optimization test functions in Jamil and Yang (2013) [[4](#References)].\n",
    "\n",
    "* Continous\n",
    "    * A continuous function has a value $f(x)$ for every $x$ in the domain $S$.\n",
    "* Differentiable\n",
    "    * Gradient-based optimization methods require the objective function to be differentiable within the search domain.  Differentiability criteria are similarly restricted.  \n",
    "    * Examples of functions that are not differentiable are absolute value functions or functions with asymptote(s).\n",
    "* Unimodal vs. Multimodal\n",
    "    * A unimodal function is either monotonically increasing or decreasing.  As a result, the function has no local minima or maxima.\n",
    "    * Some of the optimization algorithms we investigate use only gradient information to decide on a search direction and are guaranteed to find the global minimum of a unimodal function.  Although gradient-based methods can be used on multimodal functions, there is no guarantee that they will find the global minimum.\n",
    "* Separability\n",
    "    * Degree of independence of parameters in the search space.\n",
    "    * Certain kinds of optimization algorithms can exploit separability during their search, but this criteria is unimportant for the algorithms used in this project.\n",
    "\n",
    "| Test Function Name | Criteria |\n",
    "|-|-|\n",
    "| Rosenbrock | Continuous, Differentiable, Non-Separable, Unimodal |\n",
    "| Goldstein-Price | Continuous, Differentiable, Non-Separable, Multimodal |\n",
    "| Bartels-Conn | Continuous, Non-Differentiable, Non-Separable, Multimodal |\n",
    "| Egg Crate | Continuous, Differentiable, Separable, Multimodal |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rosenbrock Function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Rosenbrock function [[3](#References)] is characterized by a parabolic shaped flat valley bounded by steep canyon walls. \n",
    "\n",
    "$$\n",
    "f(x_1, \\cdots, x_D) = \\sum_{i}^{D-1} [ 100(x_{i+1} - x_i^2)^2 + (1 - x_i)^2]\n",
    "$$\n",
    "\n",
    "Global minimum is located at $x^* = f(1, \\cdots, 1)$ and $f(x^*) = 0$.\n",
    "\n",
    "### Classification\n",
    "* Continuous\n",
    "* Differentiable\n",
    "* Non-Separable\n",
    "* Unimodel\n",
    "\n",
    "### Implementation Details\n",
    "The domain of the Rosenbrock function used in this project is restricted to (-2,2) along both $x_1$ and $x_2$ dimensions.\n",
    "\n",
    "source code: [Rosenbrock-Function:-Source-Code](#Rosenbrock-Function:-Source-Code)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rosenbrock: Visualization\n",
    "A 2d filled contour plot of the Rosenbrock function is shown below.\n",
    "* The global minimum is labeled with a diamond.\n",
    "* There are 12 regularly sampled points superimposed on the plot and labeled as $x_0$ that are used to initialize separate trials of each optimization algorithm.\n",
    "\n",
    "<div align=\"middle\">\n",
    "    <img src=\"./sims/rosenbrock-plot2d.png\">\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Goldstein-Price Function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Goldstein-Price [[4](#References)] function has a single global minimum, but many local minima and several orders of magnitude difference in range.\n",
    "\n",
    "$$\n",
    "f(x_1, x_2) = [1 + (x_1 + x_2 + 1)^2 (19 - 14 x_1 + 3 x_1^2 - 14 x_2 + 6 x_1 x_2 + 3 x_2^2)] \\\\ \\times [30 + (2 x_1 - 3 x_2)^2 (18 - 32 x_1 + 12 x_1^2 + 48 x_2 - 36 x_1 x_2 + 27 x_2^2)]\n",
    "$$\n",
    "\n",
    "Global minimum is located at $x^* = f(0, -1)$ and $f(x^*) = 3$.\n",
    "\n",
    "### Classification\n",
    "* Continuous\n",
    "* Differentiable\n",
    "* Non-Separable\n",
    "* Multimodal\n",
    "\n",
    "### Implementation Details\n",
    "The domain of the Goldstein-Price function used in this project is restricted to (-2,2) along both $x_1$ and $x_2$ dimensions.\n",
    "\n",
    "source code: [Goldstein-Price-Function:-Source-Code](#Goldstein-Price-Function:-Source-Code)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Goldstein-Price: Visualization\n",
    "A 2d filled contour plot of the Goldstein-Price function is shown below.\n",
    "* The global minimum is labeled with a diamond.\n",
    "* There are 12 regularly sampled points superimposed on the plot and labeled as $x_0$ that are used to initialize separate trials of each optimization algorithm.\n",
    "\n",
    "<div align=\"middle\">\n",
    "    <img src=\"./sims/goldstein_price-plot2d.png\">\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bartels-Conn Function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Bartels-Conn [[2](#References)] function is characterized by a central elliptical valley with discontinuities appearing at all points along the descent direction. \n",
    "\n",
    "$$\n",
    "f(x_1, x_2) =  |x_1^2 + x_2^2 + x_1 x_2| + |\\sin(x_1)| + |\\cos(x_2)|\n",
    "$$\n",
    "\n",
    "Global minimum is located at $x^* = f(0, 0)$ and $f(x^*) = 1$.\n",
    "\n",
    "### Classification\n",
    "* Continuous\n",
    "* Non-Differentiable\n",
    "* Non-Separable\n",
    "* Multimodal\n",
    "\n",
    "### Implementation Details\n",
    "The domain of the Bartels-Conn function used in this project is restricted to (-5,5) along both $x_1$ and $x_2$ dimensions.\n",
    "\n",
    "source code: [Bartels-Conn-Function:-Source-Code](#Bartels-Conn-Function:-Source-Code)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bartels-Conn: Visualization\n",
    "A 2d filled contour plot of the Bartels-Conn function is shown below.\n",
    "* The global minimum is labeled with a diamond.\n",
    "* There are 12 regularly sampled points superimposed on the plot and labeled as $x_0$ that are used to initialize separate trials of each optimization algorithm.\n",
    "\n",
    "<div align=\"middle\">\n",
    "    <img src=\"./sims/bartels_conn-plot2d.png\">\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Egg Crate Function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Egg Crate function [[2](#References)] is characterized by many local minima arranged in a two-dimensional grid.  The minimum value of each local minima becomes progressively lower as you move towards the origin.\n",
    "\n",
    "$$\n",
    "f(x_1, x_2) = x_1^2 + x_2^2 + 25 (\\sin^2(x_1) + \\sin^2(x_2))\n",
    "$$\n",
    "\n",
    "Global minimum is located at $x^* = f(0, 0)$ and $f(x^*) = 0$.\n",
    "\n",
    "### Classification\n",
    "* Continuous\n",
    "* Differentiable\n",
    "* Separable\n",
    "* Multimodal\n",
    "\n",
    "### Implementation Details\n",
    "The domain of the Egg Crate function used in this project is restricted to (-5,5) along both $x_1$ and $x_2$ dimensions.\n",
    "\n",
    "source code: [Egg-Crate-Function:-Source-Code](#Egg-Crate-Function:-Source-Code)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Egg Crate: Visualization\n",
    "A 2d filled contour plot of the Egg Crate function is shown below.\n",
    "* The global minimum is labeled with a diamond.\n",
    "* There are 12 regularly sampled points superimposed on the plot and labeled as $x_0$ that are used to initialize separate trials of each optimization algorithm.\n",
    "\n",
    "<div align=\"middle\">\n",
    "    <img src=\"./sims/egg_crate-plot2d.png\">\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimization Algorithms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The table below summarizes the algorithms used in the project. Each algorithm is classified according to criteria that are based on the de- scription of algorithms for optimzation provided in Kochenderfer and Wheeler (2019) [4].\n",
    "\n",
    "* Gradient-based\n",
    "    * Gradient-based methods use the derivative of the function to decide on a search direction.  These methods are further subcategorized into whether they use only first-order or second-order information.\n",
    "    * Gradient-based methods are restricted to unimodal or nearly unimodal test functions such as Rosenbrock and Goldstein-Price.\n",
    "* Stochastic\n",
    "    * Stochastic optimization methods incorporate randomness when deciding on a search direction or distance.  Unlike gradient-based methods, these methods are not deterministic and evaluating their behavior is more challenging.\n",
    "    * Although simulated annealing and particle swarm use drastically different search methods (each is described in more detail later), the contrast to make here is between a stochastic method that explores from a single point versus a stochastic method that explores from multiple points in parallel. For these reasons we refer to simulated annealing as a stochastic serial method and particle swarm as a stochastic parallel (or population) method.\n",
    "    * Stochastic methods are not restricted to any test functions.\n",
    "\n",
    "| Method | Approach |\n",
    "|-|-|\n",
    "| Gradient Descent | Gradient-Based, First-Order |\n",
    "| BFGS | Gradient-Based, Second-Order |\n",
    "| Simulated Annealing | Stochastic, Serial |\n",
    "| Particle Swarm | Stochastic, Parallel/Population |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient Descent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Gradient Descent method [[1](#References)] is a first-order iterative optimization algorithm for finding the local minimum of a differentiable function.\n",
    "\n",
    "### Gradient Descent Algorithm\n",
    "1. Start with some initial guess $x_0$ and learning rate $\\alpha$.\n",
    "2. Update $x_k$ in the direction of negative gradient $x_k = x_{k-1} - \\alpha \\nabla f(x_{k-1})$.\n",
    "3. Evaluate the gradient at the new minimum $\\nabla f(x_k)$\n",
    "4. Repeat from step 2 until $\\nabla f(x_k) \\approx 0$\n",
    "\n",
    "### Implementation Details\n",
    "Gradient descent is implemented using the [autograd](https://autograd.readthedocs.io/en/latest/index.html) automatic differentiation library to compute the gradient of each test function.  A constant learning rate $\\alpha$ is used for all iterations of the algorithm.  The algorithm terminates when the norm of the gradient is less than _tol_.\n",
    "\n",
    "source code: [Gradient-Descent:-Source-Code](#Gradient-Descent:-Source-Code)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gradient Descent: Rosenbrock\n",
    "The result of a single trial of gradient descent on the Rosenbrock test function is shown below.  Observe that search direction is perpendicular to the isovalue given by the contour.  Although gradient descent quickly finds the narrow valley containing the global minimum, most of the iterations of the algorithm (nit=8288) are spent slowly traveling along this relatively flat surface.\n",
    "\n",
    "<div align=\"middle\">\n",
    "    <img src=\"./sims/gradient_descent-rosenbrock-plot2d-10.png\">\n",
    "</div>\n",
    "\n",
    "### Gradient Descent: Goldstein Price\n",
    "The result of a single trial of gradient descent on the Goldstein-Price test function is shown below.  Goldstein-Price is multimodal and as a result, gradient descent is not guaranteed to find the global minimum which is what happens on this trial.  Instead of finding the global minimum, the search ends when finding a flat surface inside of a local minimum ($\\min(f)=30$ instead of global $\\min(f)=3$).\n",
    "\n",
    "<div align=\"middle\">\n",
    "<img src=\"./sims/gradient_descent-goldstein_price-plot2d-05.png\">\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## BFGS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Broyden-Fletcher-Goldfarb-Shanno aka BFGS method is a second-order iterative optimization method.\n",
    "\n",
    "### Quasi-Newton Methods\n",
    "The BFGS method [[5](#References)] is referred to as Quasi-Newton in reference to the fact that unlike Newton's method which uses an explicit Hessian matrix, these methods approximate the Hessian. \n",
    "\n",
    "$$\n",
    "x_{k+1} = x_k - \\alpha_k B_k^{-1} \\nabla f(x_k)\n",
    "$$\n",
    "\n",
    "where\n",
    "* $B_k$ is approxmation to Hessian\n",
    "* $\\alpha_k$ is obtained from line search\n",
    "\n",
    "Secant updating methods have superlinear convergence ($1 < r < 2$).\n",
    "* Slower to converge than Newton's method, but cost-per-iteration is less.\n",
    "\n",
    "### BFGS Algorithm\n",
    "1. Start with some initial guess $x_0$ and approximate Hessian $B_0 = I$.\n",
    "2. Solve $B_k s_k = -\\nabla f(x_k)$ for $s_k$ or use a line search (described below) to find $s_k$.\n",
    "3. Compute $x_{k+1} = x_k + s_k$.\n",
    "4. Compute the difference in gradients $y_k = \\nabla f(x_{k+1}) - \\nabla f(x_k)$.\n",
    "5. Update approximate Hessian.\n",
    "$$\n",
    "B_{k+1} = B_k + \\frac{y_k y_k^T}{y_k^T s_k} - \\frac{B_k s_k s_k^T B_k}{s_k^T B_k s_k}\n",
    "$$\n",
    "6. Repeat from step 2 until some stopping criteria is reached.\n",
    "\n",
    "### Line Search\n",
    "A line search is used to find the distance along the descent direction of the next step of the optimization.  The scalar multiple $\\alpha$ along the descent direction $d$ is found by minimizing the function below.\n",
    "\n",
    "$$\n",
    "\\underset{\\alpha}{\\text{minimize}} f(x_k + \\alpha d)\n",
    "$$\n",
    "where\n",
    "* $f(...)$ is the function to minimize\n",
    "* $x_k$ is the current solution\n",
    "* $\\alpha$ is a scalar \n",
    "* $d$ is a vector that describes the descent direction of the function\n",
    "\n",
    "For first-order optimization problems the descent direction is given by the negative gradient $-\\nabla f(x_k)$.\n",
    "\n",
    "For second-order optimization problems the descent direction is given by by the product of the negative gradient and Hessian $-\\nabla f(x_k) H_k$.\n",
    "\n",
    "### Implementation Details\n",
    "BFGS requires the gradient of the function to minimize and similar to gradient descent the [autograd](https://autograd.readthedocs.io/en/latest/index.html) automatic differentiation library is used to compute the gradient of each test function.  We observed that the performance of the line search has a big impact on the quality of the results obtained and as a result, we replaced a simple version of BFGS written in numpy with a version from scipy having a far more sophisticated line search.  The algorithm terminates when the norm of the gradient is less than _tol_.\n",
    "\n",
    "source code: [BFGS:-Source-Code](#BFGS:-Source-Code)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### BFGS: Rosenbrock\n",
    "The result of a single trial of BFGS on the Rosenbrock test function is shown below. Observe that similar to gradient descent the search direction is perpendicular to the isovalue given by the contour.  However unlike gradient descent, the line search in the BFGS algorithm is able to exploit second-order information to take larger steps between iterations.  As a result, the algorithm converges very quickly (nit=30) in comparison to all the other methods tested.\n",
    "\n",
    "<div align=\"middle\">\n",
    "    <img src=\"./sims/bfgs-rosenbrock-plot2d-11.png\">\n",
    "</div>\n",
    "\n",
    "animation: [youtube](https://youtu.be/_HfCZAnnIgI)\n",
    "\n",
    "### BFGS: Goldstein Price\n",
    "The result of a single trial of BFGS on the Goldstein-Price test function is shown below.  Compare this result to the result initialized at the same initial position $x_0$ with gradient descent [Gradient-Descent:-Goldstein-Price](#Gradient-Descent:-Goldstein-Price).  The solution trajectory used by BFGS avoids passing through the local minimum in which gradient descent was caught.  Despite the improvement, BFGS fails to find the global minimum on 6 of the 12 trials of the Goldstein-Price test function.\n",
    "\n",
    "<div align=\"middle\">\n",
    "<img src=\"./sims/bfgs-goldstein_price-plot2d-05.png\">\n",
    "</div>\n",
    "\n",
    "animation: [youtube](https://youtu.be/fIyGMrPIsGk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulated Annealing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simulated annealing [[1](#References)] is a stochastic optimization method based on the natural physical optimization process that occurs when a material is heated to a relatively high temperature and allowed to cool.  At high temperature the atoms in the material more readily break apart and redistribute allowing the material to become more easily deformed and disordered.  As the material cools, the amount of free energy needed for such motion decreases and the material hardens into an ordered crystal structure.\n",
    "\n",
    "In the context of optimization, this process suggests two mechanisms:\n",
    "* A means by which the search continues in the direction of the local minimum or restarts in a new position that might be initially worse than the current local minimum. \n",
    "* A slow decrease in the probability that the algorithm restarts the search in some other position.\n",
    "\n",
    "### Transition Distribution\n",
    "The mean and covariance of the transition distribution is used to select a new position.  The new position is described in terms of an offset from the current position according to a multivariate normal distribution.\n",
    "\n",
    "### Annealing Schedule\n",
    "The annealing schedule describes the probability $p(z)$ that the algorithm restarts the search in some other position.  The initial value and rate of decay are parameters of the algorithm.\n",
    "\n",
    "### Simulated Annealing Algorithm\n",
    "1. Start with some initial guess $x_0$ and set this as the global minimum $f(x_{min}) = f(x_0)$.\n",
    "2. Generate a new position $x_k$ by adding to $x_{k-1}$ an offset randomly chosen from the transition distribution.\n",
    "3. Evaluate the function at the new position and compute the change in the objective function $\\Delta f(x_k) = f(x_k) - f(x_{k-1})$.\n",
    "4. If the objective function is improved $\\Delta f(x_k) < 0$, then move to the new position, else use the annealing schedule to compute the probability that despite the lack of improvement a change in position is still made.\n",
    "5. If the function evaluated at this position is less than the global minimum, then update the global minimum $f(x_{min}) = f(x_k)$.\n",
    "6. Repeat from step 2 until the number of iterations are reached.\n",
    "\n",
    "### Implementation Details\n",
    "Simulated annealing is implemented in numpy with a multivariate normal initialized from a list of means (set to 1. for all trials) and shared covariance matrix (set to the identity matrix for all trials).  Points sampled from this multivariate normal that fall outside of the domain are clipped to the boundary.  More effort could be spent to tune the distribution to each test function, but that wasn't done.  The initial temperature used by the annealing schedule is set to $T_0=1.0$ for all simulations.\n",
    "\n",
    "source code: [Simulated-Annealing:-Source-Code](#Simulated-Annealing:-Source-Code)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulated Annealing: Rosenbrock\n",
    "The result of a single trial of simulated annealing on the Rosenbrock test function is shown below.  The algorithm maintains a history of the absolute minimum that it has found and the location where the algorithm is searching does not necessarily corresond to the algorithm minimum.  To help visualize the algorithm progress, the plots and animation show the evolution of the minimum rather than current location.\n",
    "\n",
    "<div align=\"middle\">\n",
    "    <img src=\"./sims/simulated_annealing-rosenbrock-plot2d-10.png\">\n",
    "</div>\n",
    "\n",
    "animation: [youtube](https://youtu.be/dcKCRAYu-Oo)\n",
    "\n",
    "### Simulated Annealing: Goldstein-Price\n",
    "The result of a single trial of simulated annealing on the Goldstein-Price test function is shown below.  Similar to the previous example, the algorithm does not make rapid progress towards the minimum until about 1000 iterations have passed.\n",
    "\n",
    "<div align=\"middle\">\n",
    "    <img src=\"./sims/simulated_annealing-goldstein_price-plot3d-06.png\">\n",
    "</div>\n",
    "\n",
    "animation: [youtube](https://youtu.be/AgMDXNWJH24)\n",
    "\n",
    "### Simulated Annealing: Bartels-Conn\n",
    "The result of a single trial of simulated annealing on the Bartels-Conn test function is shown below.  Since the Bartels-Conn function is not differentiable, we cannot use a gradient-based solver.  However, we can compare the progress of the simulated annealing algorithm along the downhill direction of the test surface.  In contrast to what we would see from gradient descent or BFGS, the progress downhill zig-zags randomly since there is no gradient information to direct progress.  Nevertheless the simulated annealing algorithm finds the global minimum after 200 iterations.\n",
    "\n",
    "<div align=\"middle\">\n",
    "    <img src=\"./sims/simulated_annealing-bartels_conn-plot2d-12.png\">\n",
    "</div>\n",
    "\n",
    "animation: [youtube](https://youtu.be/KKK3SiV80Ls)\n",
    "\n",
    "### Simulated Annealing: Egg Crate\n",
    "The result of a single trial of simulated annealing on the Egg Crate test function is shown below.  This test function has 8 local minima surrounding a global minimum.  The simulated annealing function spends about 1000 iterations exploring 2 of the surrounding local minima before finding the central global minimum.  After so many iterations, the annealing process will discourage restarts and the algorithm spends another 500 iterations to get within 2 decimal places of accuracy ($\\min(f)=2.06e-02$) of the global minimum.\n",
    "\n",
    "<div align=\"middle\">\n",
    "    <img src=\"./sims/simulated_annealing-egg_crate-plot3d-01.png\">\n",
    "</div>\n",
    "\n",
    "animation: [youtube](https://youtu.be/bfBrm2unoOg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Particle Swarm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Particle swarm [[6](#References)] is a stochastic optimization method based on particles at different positions that simultaneously explore the optimization function and influence each other's search.  \n",
    "\n",
    "Each particle in the swarm is characterized by the following properties:\n",
    "* current position\n",
    "* current velocity\n",
    "* position of the minimum found by this particle\n",
    "\n",
    "The rule that each particle uses to update its' next position is based on the following:\n",
    "* current velocity of the particle\n",
    "* velocity in the direction of the minimum found by this particle so far\n",
    "* velocity in the direction of the global minimum found by all particles so far\n",
    "\n",
    "### Particle Swarm Algorithm\n",
    "1. Initialize a list of $n \\times p$-dimensional particles with a random initial position $(x_{1,0}, \\cdots, x_{n,0})$ and random velocity $(v_{1,0}, \\cdots, v_{n,0})$.\n",
    "2. For each particle save the position and minimum function value found by that particle $(x_{1,\\min}, \\cdots, x_{n,\\min})$ and $(f(x_{1,\\min}), \\cdots, f(x_{n,\\min}))$.\n",
    "3. Save the global position and minimum function value $x_\\min$ and $f(x_\\min)$.\n",
    "4. Initialize a counter $k=1$ used to track the iteration number.\n",
    "5. Update the velocity of each particle for the current iteration $v_{i,k} = \\omega v_{i,k-1} + p_1 r_{i,1} (x_{i,\\min} - x_{i,k-1}) + p_2 r_{i,2} (x_\\min - x_{i,k-1})$ where $\\omega$ is a inertia constant, $p_1, p_2$ are momentum constants, and $r_{i,1}, r_{i,2}$ are per-particle random numbers in the range $[0,1]$.\n",
    "6. Update the position of each particle for the current iteration $x_{i,k} = x_{i,k-1} + v_{i,k}$.\n",
    "7. Evaluate the objective function at each new position and update the per-particle and global minimum function values and position.\n",
    "8. Repeat from step 4 until the number of iterations are reached.\n",
    "\n",
    "### Implementation Details\n",
    "Particle swarm is implemented in numpy. There are a few algorithm hyperparameters which are set to default values for all trials:\n",
    "* $\\omega$ inertia coefficient (default: 1.0)\n",
    "* $p_1$ momentum coefficient towards min position of current particle (default: 1.0)\n",
    "* $p_2$ momentum coefficient towards min position among all particles (default: 1.0)\n",
    "\n",
    "An additional hyperparamter controls the number of particles $n$ used in the swarm.  Each particle will participate in the search during an iteration.  As a result, the effective number of iterations of the algorithm is the number of particles multiplied by the number of iterations.  The effective number of iterations (nit) is shown in the plots (the animations have a per-particle iteration counter that ends at nit/n).\n",
    "\n",
    "Based on experiments adjusting the value of $n$ such that the effective number of iterations is held constant, we observed that the algorithm will always find the global minimum when $n \\geq 5$.  A value of $n=3$ was chosen so that the algorithm doesn't always find the global minimum. \n",
    "\n",
    "source code: [Particle-Swarm:-Source-Code](#Particle-Swarm:-Source-Code)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Particle Swarm: Rosenbrock\n",
    "The result of a single trial of particle swarm on the Rosenbrock test function is shown below.  To help visualize the algorithm progress, the plots and animation show the trajectory of each particle rather than the value of the minimum position among all the particles (note the difference in convention compared to simulated annealing).\n",
    "\n",
    "For this trial $n=3$ particles start from different initial positions that are well dispersed across the test function.  One of the consistent patterns that emerges from the animations is the tendency for particles to converge.  In this trial, all particles converge at the entrance of the narrow valley leading to the global minimum.  Since the particles have only 50 iterations each, there isn't enough time to get closer.\n",
    "\n",
    "<div align=\"middle\">\n",
    "    <img src=\"./sims/particle_swarm-rosenbrock-plot3d-05.png\">\n",
    "</div>\n",
    "\n",
    "animation: [youtube](https://youtu.be/SzwsbCBg-tk)\n",
    "\n",
    "### Particle Swarm: Goldstein-Price\n",
    "The result of a single trial of particle swarm on the Goldstein-Price test function is shown below. Compare the result to the result with BFGS [BFGS:-Goldstein-Price](#BFGS:-Goldstein-Price).  Although none of the particles start as close to the global minimum, the collection of particles are able to influence each other such that they all move in the same direction down the slope. \n",
    "\n",
    "<div align=\"middle\">\n",
    "    <img src=\"./sims/particle_swarm-goldstein_price-plot2d-08.png\">\n",
    "</div>\n",
    "\n",
    "animation: [youtube](https://youtu.be/FIHklbjZrQ0)\n",
    "\n",
    "### Particle Swarm: Bartels-Conn\n",
    "The result of a single trial of particle swarm on the Bartels-Conn test function is shown below.\n",
    "\n",
    "<div align=\"middle\">\n",
    "    <img src=\"./sims/particle_swarm-bartels_conn-plot3d-02.png\">\n",
    "</div>\n",
    "\n",
    "animation: [youtube](https://youtu.be/crtcMyoKOzQ)\n",
    "\n",
    "### Particle Swarm: Egg Crate\n",
    "The result of a single trial of particle swarm on the Egg Crate test function is shown below.  This animation is an excellent demonstration of particles starting from opposite corners of the domain converging and moving together to reach the global minimum.\n",
    "\n",
    "<div align=\"middle\">\n",
    "    <img src=\"./sims/particle_swarm-egg_crate-plot2d-11.png\">\n",
    "</div>\n",
    "\n",
    "animation: [youtube](https://youtu.be/Z0m8CiTAb3M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The table below summarizes the results obtained from running 12 trials of each combination of algorithm and test function.  Each column is defined as follows:\n",
    "* alg\n",
    "    * Name of algorithm.\n",
    "* func\n",
    "    * Name of test function.\n",
    "* ntrials\n",
    "    * Number of trials\n",
    "* nmin\n",
    "    * Number of trials reaching minimum (or near-minimum) threshold.\n",
    "    * For the Rosenbrock, Bartels-Conn, and Egg Crate functions, any result with an absolute error less than 1.0 from the test surface minimum is considered to have reached the minimum.\n",
    "    * For the Goldstein-Price function the same rule applies, but the threshold is 10.0 due to the magnitude of the range across the test function.\n",
    "* mae\n",
    "    * Mean average absolute error of all trials.\n",
    "    * The absolute error is computed by taking the magnitude of the difference between the minimum reported by the algorithm and the known global minimum value of the test function.\n",
    "* mnit\n",
    "    * Mean number of iterations to reach a solution.\n",
    "    * Gradient-based algorithms will terminate at convergence, but the other algoritms are run for a fixed number of iterations.\n",
    "    * Particle swarm uses multiple particles running in parallel. The number of iterations reported for this algorithm in the table reflects the number of iterations of the algorithm multiplied by number of particles.  This makes the number of iterations of this algorithm comparable to a serial algorithm such as simulated annealing.\n",
    "* msec\n",
    "    * Mean elapsed runtime (in seconds) across all trials.\n",
    "    * Running time values are reported using a laptop from early 2015 with 2.9 GHz Dual-Core Intel i5 processor and 8GB of memory. The operating system used is ubuntu 18.04.\n",
    "    * Versions of software used appears in [Library-Versions](#Library-Versions).\n",
    "* algparm\n",
    "    * Algorithm hyperparameters and settings used for all trials.\n",
    "    * Gradient Descent: $\\alpha$ is learning rate and _tol_ is convergence threshold\n",
    "    * BFGS: _tol_ is convergence threshold\n",
    "    * Simulated Annealing: $T_0$ is the initial temperature used in the [Annealing-Schedule](#Annealing-Schedule)\n",
    "    * Particle Swarm: $n$ is the number of particles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>alg</th>\n",
       "      <th>func</th>\n",
       "      <th>ntrials</th>\n",
       "      <th>nmin</th>\n",
       "      <th>mae</th>\n",
       "      <th>mnit</th>\n",
       "      <th>msec</th>\n",
       "      <th>algparm</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>Gradient Descent</td>\n",
       "      <td>Rosenbrock</td>\n",
       "      <td>12</td>\n",
       "      <td>12</td>\n",
       "      <td>1.23e-04</td>\n",
       "      <td>8281</td>\n",
       "      <td>11.58</td>\n",
       "      <td>$\\alpha$=0.001 tol=0.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>Gradient Descent</td>\n",
       "      <td>Goldstein Price</td>\n",
       "      <td>12</td>\n",
       "      <td>3</td>\n",
       "      <td>2.04e+02</td>\n",
       "      <td>2113</td>\n",
       "      <td>25.45</td>\n",
       "      <td>$\\alpha$=1e-05 tol=0.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>BFGS</td>\n",
       "      <td>Rosenbrock</td>\n",
       "      <td>12</td>\n",
       "      <td>12</td>\n",
       "      <td>2.51e-06</td>\n",
       "      <td>28</td>\n",
       "      <td>0.38</td>\n",
       "      <td>tol=0.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>BFGS</td>\n",
       "      <td>Goldstein Price</td>\n",
       "      <td>12</td>\n",
       "      <td>6</td>\n",
       "      <td>1.49e+02</td>\n",
       "      <td>16</td>\n",
       "      <td>0.58</td>\n",
       "      <td>tol=0.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>Simulated Annealing</td>\n",
       "      <td>Rosenbrock</td>\n",
       "      <td>12</td>\n",
       "      <td>12</td>\n",
       "      <td>2.21e-01</td>\n",
       "      <td>200</td>\n",
       "      <td>0.30</td>\n",
       "      <td>$T_0$=1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>Simulated Annealing</td>\n",
       "      <td>Goldstein Price</td>\n",
       "      <td>12</td>\n",
       "      <td>7</td>\n",
       "      <td>3.65e+01</td>\n",
       "      <td>1500</td>\n",
       "      <td>1.23</td>\n",
       "      <td>$T_0$=1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>Simulated Annealing</td>\n",
       "      <td>Bartels Conn</td>\n",
       "      <td>12</td>\n",
       "      <td>11</td>\n",
       "      <td>1.80e+00</td>\n",
       "      <td>200</td>\n",
       "      <td>0.19</td>\n",
       "      <td>$T_0$=1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>Simulated Annealing</td>\n",
       "      <td>Egg Crate</td>\n",
       "      <td>12</td>\n",
       "      <td>6</td>\n",
       "      <td>6.40e+00</td>\n",
       "      <td>1500</td>\n",
       "      <td>1.43</td>\n",
       "      <td>$T_0$=1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>Particle Swarm</td>\n",
       "      <td>Rosenbrock</td>\n",
       "      <td>12</td>\n",
       "      <td>10</td>\n",
       "      <td>1.23e+00</td>\n",
       "      <td>150</td>\n",
       "      <td>0.07</td>\n",
       "      <td>n=3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>Particle Swarm</td>\n",
       "      <td>Goldstein Price</td>\n",
       "      <td>12</td>\n",
       "      <td>9</td>\n",
       "      <td>2.24e+01</td>\n",
       "      <td>750</td>\n",
       "      <td>0.38</td>\n",
       "      <td>n=3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>Particle Swarm</td>\n",
       "      <td>Bartels Conn</td>\n",
       "      <td>12</td>\n",
       "      <td>11</td>\n",
       "      <td>9.21e-01</td>\n",
       "      <td>150</td>\n",
       "      <td>0.15</td>\n",
       "      <td>n=3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>Particle Swarm</td>\n",
       "      <td>Egg Crate</td>\n",
       "      <td>12</td>\n",
       "      <td>6</td>\n",
       "      <td>2.71e+00</td>\n",
       "      <td>750</td>\n",
       "      <td>0.43</td>\n",
       "      <td>n=3</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(HTML('./sims/results-sum.html'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Discussion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusions\n",
    "The following conclusions can be drawn from the comparative results presented.\n",
    "\n",
    "* When using gradient-based solvers, first-order methods such as gradient descent require more iterations to find the global minimum than second-order methods such as BFGS. Use of a dyanamic learning rate with gradient descent would reduce the number of iterations required, but would not change the fundamental conclusion.\n",
    "* Gradient-based solvers achieve higher levels of accuracy than stochastic solvers. In our tests gradient based solvers achieve about 6 decimal digits of precision on the Rosenbrock function, whereas stochastic solvers achieve at most 1 decimal digit of precision on the same test function despite using more iterations.\n",
    "* Stochastic solvers can sometimes find the global minimum despite getting caught in a local minimum. In contrast, a gradient-based solver cannot escape a local minimum. In our tests simulated annealing and particle swarm are able to find the global minimum of the Goldstein-Price function more frequently than BFGS and gradient descent. Increasing the number of iterations given to BFGS and gradient descent would not change this outcome.\n",
    "* Stochastic solvers require more iterations than gradient-based solvers, but less computational effort per iteration. The mean runtime per iteration (msec/mnit) of simulated annealing and particle swarm on the Rosenbrock and Goldstein-Price function is less than BFGS and gradient descent.\n",
    "* Particle swarm is more computationally efficient than simulated annealing. Both particle swarm and simulated annealing reached the benchmark on 36 of 48 trials, but particle swarm required 1800 mean number of iterations (mnit) whereas simulated annealing required 3400 mnit. \n",
    "\n",
    "## Future Work\n",
    "The following should be considered for future work.\n",
    "\n",
    "* All of the test functions are from $f: \\mathbb{R}^2$. Repeating these experiments on higher dimensional test functions could produce different results. Interesting properties to study would be the growth in number of iterations required by the stochastic solvers versus the growth in computational cost of the gradient-based methods on higher dimensional surfaces.\n",
    "* One alternative that wasn't explored by this study is to use repeated trials of stochastic solvers initialized at different initial points rather than increasing the number of iterations from the same initial point."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[1]\n",
    "> Mykel J. Kochenderfer and Tim A. Wheeler. 2019. Algorithms for Optimization. The MIT Press.\n",
    "\n",
    "[2]\n",
    "> Momin Jamil and Xin-She Yang, A literature survey of benchmark functions for global optimization problems, Int. Journal of Mathematical Modelling and Numerical Optimisation, Vol. 4, No. 2, pp. 150–194 (2013). DOI: 10.1504/IJMMNO.2013.055204\n",
    "\n",
    "[3]\n",
    "> H. H. Rosenbrock, “An Automatic Method for Finding the Greatest or least Value of a Function,” Computer Journal, vol. 3, no. 3, pp. 175-184, 1960. [Available Online]: http://comjnl.oxfordjournals.org/content/3/3/175.full.pdf\n",
    "\n",
    "[4]\n",
    "> A. A. Goldstein, J. F. Price, “On Descent from Local Minima,” Mathematics and Comptutaion, vol. 25, no. 115, pp. 569-574, 1971.\n",
    "\n",
    "[5]\n",
    "> Michael T. Heath. 2018. Scientific Computing: An Introductory Survey, Revised Second Edition. SIAM-Society for Industrial and Applied Mathematics, Philadelphia, PA, USA.\n",
    "\n",
    "[6]\n",
    ">  J. Kennedy, R. C. Eberhart, and Y. Shi, Swarm Intelligence. Morgan Kaufmann, 2001."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Appendix: Links to Animations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The table belows lists links to selected 2d animations referred to in the project and uploaded to youtube.\n",
    "\n",
    "| Algorithm | Surface | Trial | Youtube Link |\n",
    "|-|-|-|-|\n",
    "| BFGS | Rosenbrock | 1 | https://youtu.be/PDk9d_65sHs |\n",
    "| BFGS | Rosenbrock | 11 | https://youtu.be/_HfCZAnnIgI |\n",
    "| BFGS | Goldstein-Price | 4 | https://youtu.be/Y01H7iUr6js |\n",
    "| BFGS | Goldstein-Price | 5 | https://youtu.be/fIyGMrPIsGk |\n",
    "| Simulated Annealing | Rosenbrock | 10 | https://youtu.be/dcKCRAYu-Oo |\n",
    "| Simulated Annealing | Goldstein-Price | 6 | https://youtu.be/AgMDXNWJH24 | |\n",
    "| Simulated Annealing | Bartels-Conn | 1 | https://youtu.be/wp2_u-zHY7c |\n",
    "| Simulated Annealing | Bartels-Conn | 12 | https://youtu.be/KKK3SiV80Ls |\n",
    "| Simulated Annealing | Egg Crate | 1 | https://youtu.be/bfBrm2unoOg |\n",
    "| Particle Swarm | Rosenbrock | 5 | https://youtu.be/SzwsbCBg-tk |\n",
    "| Particle Swarm | Goldstein-Price | 4 | https://youtu.be/cyyI9hAzjqg |\n",
    "| Particle Swarm | Goldstein-Price | 7 | https://youtu.be/sQjNwbgXpvc |\n",
    "| Particle Swarm | Goldstein-Price | 8 | https://youtu.be/FIHklbjZrQ0 |\n",
    "| Particle Swarm | Bartels-Conn | 2 | https://youtu.be/crtcMyoKOzQ |\n",
    "| Particle Swarm | Egg Crate | 10 | https://youtu.be/s9MEM_ML3kg |\n",
    "| Particle Swarm | Egg Crate | 11 | https://youtu.be/Z0m8CiTAb3M |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Appendix: Source Code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Source code listings for the following components of the project:\n",
    "* Test functions\n",
    "* Optimization algorithms\n",
    "* Static 2d and 3d plots\n",
    "* Animated 2d plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Library Imports\n",
    "\n",
    "Library imports for optimization algorithms.\n",
    "\n",
    "```python\n",
    "from autograd import grad\n",
    "import numpy as np\n",
    "import scipy.optimize as opt\n",
    "```\n",
    "\n",
    "Note that gradient-based algorithms such as Gradient Descent and BFGS use [autograd](https://autograd.readthedocs.io/en/latest/index.html) to wrap numpy operations with equivalent functions that support automatic differentiation (AD).  Replacing the standard numpy import with the line below is sufficient to be able to differentiate functions which are composed of numpy operations.\n",
    "\n",
    "```python\n",
    "import autograd.numpy as np\n",
    "```\n",
    "\n",
    "Library imports for plotting functions.\n",
    "\n",
    "```python\n",
    "import json\n",
    "import os\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker as ticker\n",
    "from matplotlib.animation import FuncAnimation, FFMpegWriter\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Library Versions\n",
    "Python version `3.6.9` and the packages listed in the `requirements.txt` below were used to produce the results reported in this project. \n",
    "\n",
    "```\n",
    "nbconvert==5.6.1\n",
    "jupyter~=1.0\n",
    "matplotlib~=3.3\n",
    "numpy~=1.19\n",
    "scipy~=1.5\n",
    "scikit-image~=0.17\n",
    "scikit-learn~=0.23\n",
    "pandas~=1.0\n",
    "pyvista~=0.25\n",
    "itkwidgets~=0.32\n",
    "autograd~=1.3\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Surface Generation: Source Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def surface(fx, start=-30, stop=30, num=60):\n",
    "    \"\"\"\n",
    "    surface evaluates fx at regularly spaced grid of points\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    fx : func\n",
    "        fx is a vector valued function that returns a scalar result\n",
    "    start : float\n",
    "        lower bound of the coordinate grid\n",
    "    stop : float\n",
    "        upper bound of the coordinate grid\n",
    "    num : int\n",
    "        number of points along one dimension of the grid\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    array\n",
    "        2D array formed by evaluating fx at each grid point\n",
    "    \"\"\"\n",
    "    x = np.linspace(start=start, stop=stop, num=num)\n",
    "    x1, x2 = np.meshgrid(x, x, indexing='ij')\n",
    "    X = np.vstack((x1.ravel(), x2.ravel()))\n",
    "    z = np.apply_along_axis(fx, 0, X).reshape(num,num)\n",
    "    return x1, x2, z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rosenbrock Function: Source Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rosenbrock(x):\n",
    "    \"\"\"\n",
    "    rosenbrock evaluates Rosenbrock function at vector x\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    x : array\n",
    "        x is a D-dimensional vector, [x1, x2, ..., xD]\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    float\n",
    "        scalar result\n",
    "    \"\"\"\n",
    "    D = len(x)\n",
    "    i, iplus1 = np.arange(0,D-1), np.arange(1,D)\n",
    "    return np.sum(100*(x[iplus1] - x[i]**2)**2 + (1-x[i])**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Goldstein-Price Function: Source Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def goldstein_price(x):\n",
    "    \"\"\"\n",
    "    goldstein_price evaluates Goldstein-Price function at vector x\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    x : array\n",
    "        x is a 2-dimensional vector, [x1, x2]\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    float\n",
    "        scalar result\n",
    "    \"\"\"\n",
    "    a = (x[0] + x[1] + 1)**2\n",
    "    b = 19 - 14*x[0] + 3*x[0]**2 - 14*x[1] + 6*x[0]*x[1] + 3*x[1]**2\n",
    "    c = (2*x[0] - 3*x[1])**2\n",
    "    d = 18 - 32*x[0] + 12*x[0]**2 + 48*x[1] - 36*x[0]*x[1] + 27*x[1]**2\n",
    "    return (1. + a*b) * (30. + c*d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bartels-Conn Function: Source Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bartels_conn(x):\n",
    "    \"\"\"\n",
    "    bartels_conn evaluates Bartels-Conn function at vector x\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    x : array\n",
    "        x is a 2-dimensional vector, [x1, x2]\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    float\n",
    "        scalar result\n",
    "    \"\"\"\n",
    "    a = np.abs(x[0]**2 + x[1]**2 + x[0]*x[1])\n",
    "    b = np.abs(np.sin(x[0]))\n",
    "    c = np.abs(np.cos(x[1]))\n",
    "    return a + b +c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Egg Crate Function: Source Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def egg_crate(x):\n",
    "    \"\"\"\n",
    "    egg_crate evaluates Egg Crate function at vector x\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    x : array\n",
    "        x is a 2-dimensional vector, [x1, x2]\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    float\n",
    "        scalar result\n",
    "    \"\"\"\n",
    "    return x[0]**2 + x[1]**2 + 25.*(np.sin(x[0])**2 + np.sin(x[1])**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient Descent: Source Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gradient_descent(fx, gradfx, x0, alpha, tol, maxiter):\n",
    "    \"\"\"\n",
    "    gradient_descent returns the point xk where fx is minimum\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    fx : function\n",
    "        function to minimize\n",
    "    gradfx : function\n",
    "        gradient of function to minimize\n",
    "    x0 : numpy.ndarray\n",
    "        initial guess for xk\n",
    "    alpha : float\n",
    "        learning rate\n",
    "    tol : float\n",
    "        convergence threshold\n",
    "    maxiter : int\n",
    "        maximum number of iterations\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    numpy.ndarray\n",
    "        point xk where fx is minimum\n",
    "    numpy.ndarray\n",
    "        position and value history\n",
    "        [[x0, fx(x0), gradfx(x0)],\n",
    "         [x1, fx(x1), gradfx(x1)],...]\n",
    "    \"\"\"\n",
    "\n",
    "    xk, fxk, gradfxk = x0, fx(x0), gradfx(x0)\n",
    "\n",
    "    # Save current and minimum position and value to history.\n",
    "    steps = np.zeros((maxiter, (x0.size*2)+1))\n",
    "    steps[0,:] = np.hstack((x0, fxk, gradfxk))\n",
    "\n",
    "    # Repeat up to maximum number of iterations.\n",
    "    for k in range(1,maxiter):\n",
    "\n",
    "        # Stop iteration when gradient is near zero.\n",
    "        if np.linalg.norm(gradfxk) < tol:\n",
    "            steps = steps[:-(maxiter-k),:]\n",
    "            break\n",
    "\n",
    "        # Update xk based on product of learning rate and gradient.\n",
    "        xk = xk - alpha * gradfxk\n",
    "\n",
    "        # Evaluate gradient at new value of xk.\n",
    "        gradfxk = gradfx(xk)\n",
    "\n",
    "        # Evaluate the function at new value of xk.\n",
    "        fxk = fx(xk)\n",
    "\n",
    "        # Save iteration history.\n",
    "        steps[k,:] = np.hstack((xk, fxk, gradfxk))\n",
    "\n",
    "    return xk, steps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## BFGS: Source Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NOTE(mmorais): bfgs is failing on line search, use scipy instead.\n",
    "def scipy_bfgs(fx, gradfx, x0, tol, maxiter):\n",
    "    \"\"\"\n",
    "    scipy_bfgs wraps scipy implementation of bfgs\n",
    "    \"\"\"\n",
    "    # Save current and minimum position and value to history.\n",
    "    steps = np.zeros((maxiter+1, (x0.size*2)+1))\n",
    "    steps[0,:] = np.hstack((x0, fx(x0), gradfx(x0)))\n",
    "    def make_callback():\n",
    "        k = 1\n",
    "        def callback(xk):\n",
    "            nonlocal k\n",
    "            # Save iteration history.\n",
    "            steps[k,:] = np.hstack((xk, fx(xk), gradfx(xk)))\n",
    "            k = k + 1\n",
    "        return callback\n",
    "\n",
    "    # Invoke scipy minimize with BFGS option.\n",
    "    res = opt.minimize(fx, x0, method='BFGS', jac=gradfx, tol=tol,\n",
    "                       options={'maxiter': maxiter},\n",
    "                       callback=make_callback())\n",
    "\n",
    "    # Copy OptimizeResult to equivalent returned from bfgs.\n",
    "    xk = res.x\n",
    "    if res.nit < maxiter:\n",
    "        steps = steps[:-(maxiter-res.nit),:]\n",
    "    return xk, steps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulated Annealing: Source Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulated_annealing(fx, x0, mean, cov, tk, bounds, niter):\n",
    "    \"\"\"\n",
    "    simulated_annealing returns the point xk where fx is minimum\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    fx : function\n",
    "        function to minimize\n",
    "    x0 : numpy.ndarray\n",
    "        initial guess for xk\n",
    "    mean : numpy.ndarray\n",
    "        means of multivariate normal transition distribution\n",
    "    cov : numpy.ndarray\n",
    "        covariance of multivariate normal transition distribution\n",
    "    tk : function\n",
    "        annealing schedule as a function of iteration number\n",
    "    bounds : numpy.ndarray\n",
    "        domain boundaries [x1_min, x1_max, ..., xn_min, xn_max]\n",
    "    niter : int\n",
    "        number of iterations\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    numpy.ndarray\n",
    "        point xk where fx is minimum\n",
    "    numpy.ndarray\n",
    "        current and minimum position and value history\n",
    "        [[x0, fx(x0), xk_min, fx(xk_min)],\n",
    "         [x1, fx(x1), xk_min, fx(xk_min)],...]\n",
    "    \"\"\"\n",
    "\n",
    "    # Initialize solution at x0.\n",
    "    xk, fxk = x0, fx(x0)\n",
    "    xk_min, fxk_min = xk, fxk\n",
    "\n",
    "    # Setup random transition distribution.\n",
    "    mvnorm = partial(np.random.multivariate_normal, mean, cov)\n",
    "\n",
    "    # Save current and minimum position and value to history.\n",
    "    steps = np.zeros((niter, (x0.size+1)*2))\n",
    "    steps[0,:] = np.hstack((x0, fxk, xk_min, fxk_min))\n",
    "\n",
    "    # Perform fixed number of iterations.\n",
    "    for k in range(1,niter):\n",
    "\n",
    "        # Generate a new random point.\n",
    "        xk1 = xk + mvnorm()\n",
    "        xk1 = np.clip(xk1, a_min=bounds[::2], a_max=bounds[1::2])\n",
    "\n",
    "        # Evaluate the function at the new point.\n",
    "        fxk1 = fx(xk1)\n",
    "\n",
    "        # Compute the change in the objective function.\n",
    "        delta_fxk = fxk1 - fxk\n",
    "\n",
    "        # If objective function is improved or escape current position,\n",
    "        # then update xk, fxk with the new position.\n",
    "        if delta_fxk < 0. or np.random.random() < np.exp(-fxk1/tk(k)):\n",
    "            xk, fxk = xk1, fxk1\n",
    "            if fxk1 < fxk_min:\n",
    "                xk_min, fxk_min = xk1, fxk1\n",
    "\n",
    "        # Save iteration history.\n",
    "        steps[k,:] = np.hstack((xk1, fxk1, xk_min, fxk_min))\n",
    "\n",
    "    return xk_min, steps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Particle Swarm: Source Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def particle_swarm(fx, x0s, omega, p1, p2, bounds, niter):\n",
    "    \"\"\"\n",
    "    particle_swarm returns the point xk where fx is minimum\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    fx : function\n",
    "        function to minimize\n",
    "    x0s : numpy.ndarray\n",
    "        initial positions of particles in swarm\n",
    "    omega : float\n",
    "        inertia coefficient\n",
    "    p1 : float\n",
    "        momentum coefficient towards min position of current particle\n",
    "    p2 : float\n",
    "        momentum coefficient towards min position among all particles\n",
    "    bounds : numpy.ndarray\n",
    "        domain boundaries [x1_min, x1_max, ..., xn_min, xn_max]\n",
    "    niter : int\n",
    "        number of iterations\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    numpy.ndarray\n",
    "        point xk where fx is minimum\n",
    "    numpy.ndarray\n",
    "        current and minimum position and value history for each particle\n",
    "        [[x_1,0, fx(x_1,0), xk_1,min, fx(xk_1,min),...,\n",
    "            x_n,0, fx(x_n,0), xk_n,min, fx(xk_n,min)],\n",
    "         [x_1,1, fx(x_1,1), xk_1,min, fx(xk_1,min),...,\n",
    "            x_n,1, fx(x_n,1), xk_n,min, fx(xk_n,min)],\n",
    "    \"\"\"\n",
    "\n",
    "    # Initialize swarm with position, velocity, and min position.\n",
    "    pos = np.copy(x0s)\n",
    "    x0sdelta = np.max(x0s, axis=0) - np.min(x0s, axis=0)\n",
    "    vel = (np.random.random(x0s.shape)-0.5)*x0sdelta\n",
    "    posmin, fxmin = np.copy(x0s), np.apply_along_axis(fx, 1, x0s)\n",
    "\n",
    "    # Global minimum position.\n",
    "    xk_min, fxk_min = posmin[np.argmin(fxmin),:], np.min(fxmin)\n",
    "\n",
    "    # Save position, velocity, and min position by particle to history.\n",
    "    # Also save global min position and value with each particle history.\n",
    "    npart, ndim, nvecs = x0s.shape[0], x0s.shape[1], 4\n",
    "    steps = np.zeros((npart*niter, ndim*nvecs+2))\n",
    "    steps[:npart,:] = np.hstack((pos, vel, posmin, fxmin[:,np.newaxis],\n",
    "                                 np.broadcast_to(xk_min,(npart,ndim)),\n",
    "                                 np.broadcast_to(fxk_min,(npart,1))))\n",
    "\n",
    "    # Perform fixed number of iterations.\n",
    "    for k in range(1,niter):\n",
    "\n",
    "        # Compute new velocity of each particle.\n",
    "        rs = np.random.random((npart,2))\n",
    "        vel = omega*vel + p1*rs[0]*(posmin-pos) + p2*rs[1]*(xk_min-pos)\n",
    "\n",
    "        # Update the position of each particle based on velocity.\n",
    "        pos = pos + vel\n",
    "        pos = np.clip(pos, a_min=bounds[::2], a_max=bounds[1::2])\n",
    "\n",
    "        # Evaluate the objective function at each new position.\n",
    "        fxpart = np.apply_along_axis(fx, 1, pos)\n",
    "\n",
    "        # If objective function is improved,\n",
    "        # then replace particle minimum position and value.\n",
    "        inds = fxpart < fxmin\n",
    "        posmin[inds,:], fxmin[inds] = pos[inds,:], fxpart[inds]\n",
    "\n",
    "        # If global objective function is improved,\n",
    "        # then replace global minimum position and value.\n",
    "        ind = np.argmin(fxmin)\n",
    "        if fxmin[ind] < fxk_min:\n",
    "            xk_min, fxk_min = posmin[ind,:], fxmin[ind]\n",
    "\n",
    "        # Save particle history.\n",
    "        ind0 = k*npart\n",
    "        steps[ind0:ind0+npart,:] = (\n",
    "            np.hstack((pos, vel, posmin, fxmin[:,np.newaxis],\n",
    "                       np.broadcast_to(xk_min,(npart,ndim)),\n",
    "                       np.broadcast_to(fxk_min,(npart,1)))))\n",
    "\n",
    "    return xk_min, steps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Static 2d Simulation Result Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_steps(**params):\n",
    "    \"\"\"Return solution steps based on simulation properties.\"\"\"\n",
    "    savefn = os.path.join(params['base_dirn'],\n",
    "                          params['savefn_fmt'].format(**params))\n",
    "    return np.load(savefn)\n",
    "\n",
    "\n",
    "def load_meta(**params):\n",
    "    \"\"\"Return metafile based on simulation properties.\"\"\"\n",
    "    metafn = os.path.join(params['base_dirn'],\n",
    "                          params['metafn_fmt'].format(**params))\n",
    "    return json.load(open(metafn, 'r'))\n",
    "\n",
    "\n",
    "def plot2d_solutions(**params):\n",
    "    \"\"\"\n",
    "    plot2d_solutions creates 2d solution plot from simulation results\n",
    "    \"\"\"\n",
    "    algstr = params['alg'].replace('_',' ').title()\n",
    "    funcstr = params['func'].replace('_',' ').title()\n",
    "    ngridpts = params.get('ngridpts', 500)\n",
    "    bounds = params['bounds']\n",
    "    trial = params['trial']  # Single trial only.\n",
    "    xkmind = params.get('xkmind', slice(2))\n",
    "    color = params.get('color', 'darkorange')\n",
    "    ticker_locator = params.get('ticker_locator', 'LinearLocator')\n",
    "    colorbar_label = params.get('colorbar_label', 'z')\n",
    "    show_legend = params.get('show_legend', True)\n",
    "\n",
    "    # Imbue title with simulation meta information.\n",
    "    meta = load_meta(**params)\n",
    "    expmin, expxkmin = meta['exp_fxkmin'], meta['exp_xkmin']\n",
    "    expminstr = 'abs $\\\\min(f)$={0:.0f}'.format(expmin)\n",
    "    algstr = algstr if len(algstr) > 4 else algstr.upper()\n",
    "    algmeta = [('nx0','n'),('T0','$T_0$'),\n",
    "               ('alpha','$\\\\alpha$'),('tol','tol')]\n",
    "    algmetastr = ' '.join(['{0}={1}'.format(n2, meta[n1])\n",
    "                           for n1, n2 in algmeta if n1 in meta])\n",
    "    nitstr = 'nit={0:d}'.format(meta['nsteps'][trial-1])\n",
    "    minfx = meta['f(xk)'][trial-1]\n",
    "    minfmt = '.2e' if minfx < 1e-1 else '.1f'\n",
    "    minstr = '$\\\\min(f)$={0:{1}}'.format(minfx, minfmt)\n",
    "    metastrs = [algstr, minstr, nitstr, algmetastr]\n",
    "    titlestr = ' '.join([s for s in metastrs if len(s) > 0])\n",
    "    suptitlestr = 'Solution Trajectories: {0} Function'.format(funcstr)\n",
    "\n",
    "    # Generate surface for filled contour plot.\n",
    "    fx = globals()[params['func']]\n",
    "    start, stop = np.min(bounds[::2]), np.max(bounds[1::2])\n",
    "    x1, x2, z = surface(fx, start, stop, ngridpts)\n",
    "\n",
    "    fig = plt.figure(figsize=(8,6))\n",
    "\n",
    "    # Plot 2d filled contour.\n",
    "    locator = getattr(ticker, ticker_locator)\n",
    "    cs = plt.contourf(x1, x2, z, locator=locator(), cmap='viridis_r',\n",
    "                      alpha=0.7)\n",
    "\n",
    "    # Plot expected minimum.\n",
    "    plt.scatter(expxkmin[0], expxkmin[1], marker='D', c='red', s=30,\n",
    "                label=expminstr)\n",
    "\n",
    "    # Plot initial point.\n",
    "    x0 = np.array(meta['x0'][trial-1]).reshape(-1,2)\n",
    "    plt.scatter(x0[:,0], x0[:,1], marker='X', c='dodgerblue', s=30,\n",
    "                label='$x_0$')\n",
    "\n",
    "    # Plot solution trajectory.\n",
    "    steps = load_steps(**params)\n",
    "    nx0 = meta.get('nx0', 1)  # Multiple particles?\n",
    "    xks = steps[:,xkmind]\n",
    "    xks = np.clip(xks, a_min=bounds[::2], a_max=bounds[1::2])\n",
    "    nxks = 0 if np.isnan(xks).any() else len(xks)//nx0\n",
    "    for p in range(nx0):\n",
    "        p0, pN, pstep = p, nxks, nx0\n",
    "        plt.plot(xks[p0:pN:pstep,0], xks[p0:pN:pstep,1],\n",
    "                 marker='.', ms=5, markevery=0.25,\n",
    "                 ls='-', lw=1, c=color,\n",
    "                 label='$x_k$, trial={:d}'.format(trial))\n",
    "\n",
    "    plt.suptitle(suptitlestr)\n",
    "    plt.title(titlestr)\n",
    "    plt.xlabel('x1')\n",
    "    plt.xlim(bounds[:2])\n",
    "    plt.ylabel('x2')\n",
    "    plt.ylim(bounds[2:])\n",
    "    plt.colorbar(cs, label=colorbar_label)\n",
    "    if show_legend:\n",
    "        plt.legend()\n",
    "    if params.get('plot2dfn_fmt') is not None:\n",
    "        imgn = params['plot2dfn_fmt'].format(**params)\n",
    "        plotfn = os.path.join(params['base_dirn'], imgn)\n",
    "        plt.savefig(plotfn)\n",
    "    else:\n",
    "        plt.show()\n",
    "    plt.close(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Static 3d Simulation Result Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_steps(**params):\n",
    "    \"\"\"Return solution steps based on simulation properties.\"\"\"\n",
    "    savefn = os.path.join(params['base_dirn'],\n",
    "                          params['savefn_fmt'].format(**params))\n",
    "    return np.load(savefn)\n",
    "\n",
    "\n",
    "def load_meta(**params):\n",
    "    \"\"\"Return metafile based on simulation properties.\"\"\"\n",
    "    metafn = os.path.join(params['base_dirn'],\n",
    "                          params['metafn_fmt'].format(**params))\n",
    "    return json.load(open(metafn, 'r'))\n",
    "\n",
    "\n",
    "def plot3d_solutions(**params):\n",
    "    \"\"\"\n",
    "    plot3d_solutions creates 3d solution plot from simulation results\n",
    "    \"\"\"\n",
    "    algstr = params['alg'].replace('_',' ').title()\n",
    "    funcstr = params['func'].replace('_',' ').title()\n",
    "    ngridpts = params.get('ngridpts', 500)\n",
    "    bounds = params['bounds']\n",
    "    elev = params['elev']\n",
    "    azim = params['azim']\n",
    "    trial = params['trial']  # Single trial only.\n",
    "    xkmind = params.get('xkmind', slice(2))\n",
    "    color = params.get('color', 'crimson')\n",
    "    show_legend = params.get('show_legend', True)\n",
    "\n",
    "    # Imbue title with simulation meta information.\n",
    "    meta = load_meta(**params)\n",
    "    expmin, expxkmin = meta['exp_fxkmin'], meta['exp_xkmin']\n",
    "    expminstr = 'abs $\\\\min(f)$={0:.0f}'.format(expmin)\n",
    "    algstr = algstr if len(algstr) > 4 else algstr.upper()\n",
    "    algmeta = [('nx0','n'),('T0','$T_0$'),\n",
    "               ('alpha','$\\\\alpha$'),('tol','tol')]\n",
    "    algmetastr = ' '.join(['{0}={1}'.format(n2, meta[n1])\n",
    "                           for n1, n2 in algmeta if n1 in meta])\n",
    "    nitstr = 'nit={0:d}'.format(meta['nsteps'][trial-1])\n",
    "    minfx = meta['f(xk)'][trial-1]\n",
    "    minfmt = '.2e' if minfx < 1e-1 else '.1f'\n",
    "    minstr = '$\\\\min(f)$={0:{1}}'.format(minfx, minfmt)\n",
    "    metastrs = [algstr, minstr, nitstr, algmetastr]\n",
    "    titlestr = ' '.join([s for s in metastrs if len(s) > 0])\n",
    "    suptitlestr = 'Solution Trajectories: {0} Function'.format(funcstr)\n",
    "\n",
    "    # Generate surface for filled contour plot.\n",
    "    fx = globals()[params['func']]\n",
    "    start, stop = np.min(bounds[::2]), np.max(bounds[1::2])\n",
    "    x1, x2, z = surface(fx, start, stop, ngridpts)\n",
    "\n",
    "    fig = plt.figure(figsize=(10,8))\n",
    "    ax = fig.gca(projection='3d')\n",
    "    ax.view_init(elev=elev, azim=azim)\n",
    "    surf = ax.plot_surface(x1, x2, z, cmap='viridis_r', alpha=0.7)\n",
    "    fig.colorbar(surf, shrink=0.5, aspect=5)\n",
    "\n",
    "    # Plot expected minimum.\n",
    "    ax.scatter3D([expxkmin[0]],[expxkmin[1]], [expmin],\n",
    "                 marker='D', c='black', s=30,\n",
    "                 label=expminstr)\n",
    "\n",
    "    # Plot initial point.\n",
    "    x0 = np.array(meta['x0'][trial-1]).reshape(-1,2)\n",
    "    ax.scatter3D(x0[:,0], x0[:,1], [fx(xk) for xk in x0],\n",
    "                 marker='X', c='dodgerblue', s=30,\n",
    "                 label='$x_0$')\n",
    "\n",
    "    # Plot solution trajectory.\n",
    "    steps = load_steps(**params)\n",
    "    nx0 = meta.get('nx0', 1)  # Multiple particles?\n",
    "    xks = steps[:,xkmind]\n",
    "    xks = np.clip(xks, a_min=bounds[::2], a_max=bounds[1::2])\n",
    "    nxks = 0 if np.isnan(xks).any() else len(xks)//nx0\n",
    "    for p in range(nx0):\n",
    "        p0, pN, pstep = p, nxks, nx0\n",
    "        ax.plot3D(xks[p0:pN:pstep,0],\n",
    "                  xks[p0:pN:pstep,1],\n",
    "                  [fx(xk) for xk in xks[p0:pN:pstep,:]],\n",
    "                  ls='-', lw=1, c=color,\n",
    "                  label='$x_k$, trial={:d}'.format(trial))\n",
    "\n",
    "    plt.suptitle(suptitlestr)\n",
    "    plt.title(titlestr)\n",
    "    plt.xlabel('x1')\n",
    "    plt.xlim(bounds[:2])\n",
    "    plt.ylabel('x2')\n",
    "    plt.ylim(bounds[2:])\n",
    "    if show_legend:\n",
    "        ax.legend()\n",
    "    if params.get('plot3dfn_fmt') is not None:\n",
    "        imgn = params['plot3dfn_fmt'].format(**params)\n",
    "        plotfn = os.path.join(params['base_dirn'], imgn)\n",
    "        plt.savefig(plotfn)\n",
    "    else:\n",
    "        plt.show()\n",
    "    plt.close(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Animated 2d Simulation Result Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_steps(**params):\n",
    "    \"\"\"Return solution steps based on simulation properties.\"\"\"\n",
    "    savefn = os.path.join(params['base_dirn'],\n",
    "                          params['savefn_fmt'].format(**params))\n",
    "    return np.load(savefn)\n",
    "\n",
    "\n",
    "def load_meta(**params):\n",
    "    \"\"\"Return metafile based on simulation properties.\"\"\"\n",
    "    metafn = os.path.join(params['base_dirn'],\n",
    "                          params['metafn_fmt'].format(**params))\n",
    "    return json.load(open(metafn, 'r'))\n",
    "\n",
    "\n",
    "def anim2d_solutions(**params):\n",
    "    \"\"\"\n",
    "    anim2d_solutions creates 2d animation from simulation results\n",
    "    \"\"\"\n",
    "    algstr = params['alg'].replace('_',' ').title()\n",
    "    funcstr = params['func'].replace('_',' ').title()\n",
    "    ngridpts = params.get('ngridpts', 500)\n",
    "    bounds = params['bounds']\n",
    "    trial = params['trial']  # Single trial only.\n",
    "    xkmind = params.get('xkmind', slice(2))\n",
    "    fxkmind = params.get('fxkmind', 2)\n",
    "    color = params.get('color', 'darkorange')\n",
    "    ticker_locator = params.get('ticker_locator', 'LinearLocator')\n",
    "    colorbar_label = params.get('colorbar_label', 'z')\n",
    "    fps = params.get('fps', 30)\n",
    "    bitrate = params.get('bitrate', 1000)\n",
    "    show_legend = params.get('show_legend', True)\n",
    "\n",
    "    # Imbue title with simulation meta information.\n",
    "    meta = load_meta(**params)\n",
    "    expmin, expxkmin = meta['exp_fxkmin'], meta['exp_xkmin']\n",
    "    expminstr = 'abs $\\\\min(f)$={0:.0f}'.format(expmin)\n",
    "    algstr = algstr if len(algstr) > 4 else algstr.upper()\n",
    "    algmeta = [('nx0','n'),('T0','$T_0$'),\n",
    "               ('alpha','$\\\\alpha$'),('tol','tol')]\n",
    "    algmetastr = ' '.join(['{0}={1}'.format(n2, meta[n1])\n",
    "                           for n1, n2 in algmeta if n1 in meta])\n",
    "    nitstr = 'nit={0:d}'.format(meta['nsteps'][trial-1])\n",
    "    minfx = meta['f(xk)'][trial-1]\n",
    "    minfmt = '.2e' if minfx < 1e-1 else '.1f'\n",
    "    minstr = '$\\\\min(f)$={0:{1}}'.format(minfx, minfmt)\n",
    "    metastrs = [algstr, minstr, nitstr, algmetastr]\n",
    "    titlestr = ' '.join([s for s in metastrs if len(s) > 0])\n",
    "    suptitlestr = 'Solution Trajectories: {0} Function'.format(funcstr)\n",
    "\n",
    "    # Generate surface for filled contour plot.\n",
    "    fx = globals()[params['func']]\n",
    "    start, stop = np.min(bounds[::2]), np.max(bounds[1::2])\n",
    "    x1, x2, z = surface(fx, start, stop, ngridpts)\n",
    "\n",
    "    fig = plt.figure(figsize=(8,6))\n",
    "\n",
    "    # Plot 2d filled contour.\n",
    "    locator = getattr(ticker, ticker_locator)\n",
    "    cs = plt.contourf(x1, x2, z, locator=locator(), cmap='viridis_r',\n",
    "                      alpha=0.7)\n",
    "\n",
    "    # Plot expected minimum.\n",
    "    plt.scatter(expxkmin[0], expxkmin[1], marker='D', c='red', s=30,\n",
    "                label=expminstr)\n",
    "\n",
    "    # Plot initial point.\n",
    "    x0 = np.array(meta['x0'][trial-1]).reshape(-1,2)\n",
    "    plt.scatter(x0[:,0], x0[:,1], marker='X', c='dodgerblue', s=30,\n",
    "                label='$x_0$')\n",
    "\n",
    "    # Plot bounds.\n",
    "    plt.suptitle(suptitlestr)\n",
    "    plt.title(titlestr)\n",
    "    plt.xlabel('x1')\n",
    "    plt.xlim(bounds[:2])\n",
    "    plt.ylabel('x2')\n",
    "    plt.ylim(bounds[2:])\n",
    "    plt.colorbar(cs, label=colorbar_label)\n",
    "\n",
    "    # Load solution trajectory.\n",
    "    steps = load_steps(**params)\n",
    "    nx0 = meta.get('nx0', 1)  # Multiple particles?\n",
    "    xks, fxks = steps[:,xkmind], steps[:,fxkmind]\n",
    "    xks = np.clip(xks, a_min=bounds[::2], a_max=bounds[1::2])\n",
    "    nxks = 0 if np.isnan(xks).any() else len(xks)//nx0\n",
    "\n",
    "    txtx = bounds[0] + (bounds[1]-bounds[0])*0.5\n",
    "    txty = bounds[3] - (bounds[3]-bounds[2])*0.025\n",
    "    txt = plt.text(txtx, txty, '', ha='center', va='top')\n",
    "    lns = []\n",
    "    for _ in range(nx0):\n",
    "        ln, = plt.plot([], [],\n",
    "                       ls=(0, (1,1)), lw=2, c=color,\n",
    "                       label='$x_k$, trial={:d}'.format(trial))\n",
    "        lns.append(ln)\n",
    "\n",
    "    def update(ind):\n",
    "        # All of the particles store the same min.\n",
    "        txt.set_text('k={0:d} $f(x_k)$={1:{2}}'.format(\n",
    "                     ind+1, fxks[ind*nx0],\n",
    "                     '.2e' if fxks[ind*nx0] < 1e-1 else '.1f'))\n",
    "        for p in range(nx0):\n",
    "            # Starting from first position up to + incl current position.\n",
    "            p0, pN, pstep = p, (ind*nx0)+p+nx0, nx0\n",
    "            lns[p].set_data(xks[p0:pN:pstep,0], xks[p0:pN:pstep,1])\n",
    "        return lns\n",
    "\n",
    "    if show_legend:\n",
    "        plt.legend()\n",
    "\n",
    "    anim = FuncAnimation(fig, update, frames=nxks, blit=True)\n",
    "    if params.get('anim2dfn_fmt') is not None:\n",
    "        imgn = params['anim2dfn_fmt'].format(**params)\n",
    "        animfn = os.path.join(params['base_dirn'], imgn)\n",
    "        writer = FFMpegWriter(fps=fps, bitrate=bitrate,\n",
    "                              extra_args=['-vcodec', 'libx264'])\n",
    "        anim.save(animfn, writer=writer)\n",
    "    plt.close(fig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}