AIND-Simulated_Annealing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# I. Introduction\n",
    "\n",
    "![Simulated Annealing](SA_animation.gif)\n",
    "As illustrated in the lectures, simulated annealing is a probablistic technique used for finding an approximate solution to an optimization problem.  In this exercise you will check your understanding by implementing [simulated annealing](https://en.wikipedia.org/wiki/Simulated_annealing) to solve the [Traveling Salesman Problem](https://en.wikipedia.org/wiki/Travelling_salesman_problem) (TSP) between US state capitals.  Briefly, the TSP is an optimization problem that seeks to find the shortest path passing through every city exactly once.  In our example the TSP path is defined to start and end in the same city (so the path is a closed loop).\n",
    "\n",
    "Image Source: [Simulated Annealing - By Kingpin13 (Own work) [CC0], via Wikimedia Commons (Attribution not required)](https://commons.wikimedia.org/wiki/File:Hill_Climbing_with_Simulated_Annealing.gif)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "Students should read through the code, then:\n",
    "\n",
    "  0. Implement the `simulated_annealing()` main loop function in Section II\n",
    "  0. Complete the `TravelingSalesmanProblem` class by implementing the `successors()` and `get_value()` methods in section III\n",
    "  0. Complete the `schedule()` function to define the temperature schedule in Section IV\n",
    "  0. Use the completed algorithm and problem description to experiment with simulated annealing to solve larger TSP instances on the map of US capitals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import json\n",
    "import copy\n",
    "\n",
    "import numpy as np  # contains helpful math functions like numpy.exp()\n",
    "import numpy.random as random  # see numpy.random module\n",
    "# import random  # alternative to numpy.random module\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.image as mpimg\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\"\"\"Read input data and define helper functions for visualization.\"\"\"\n",
    "\n",
    "# Map services and data available from U.S. Geological Survey, National Geospatial Program.\n",
    "# Please go to http://www.usgs.gov/visual-id/credit_usgs.html for further information\n",
    "map = mpimg.imread(\"map.png\")  # US States & Capitals map\n",
    "\n",
    "# List of 30 US state capitals and corresponding coordinates on the map\n",
    "with open('capitals.json', 'r') as capitals_file:\n",
    "    capitals = json.load(capitals_file)\n",
    "capitals_list = list(capitals.items())\n",
    "\n",
    "def show_path(path, starting_city, w=12, h=8):\n",
    "    \"\"\"Plot a TSP path overlaid on a map of the US States & their capitals.\"\"\"\n",
    "    x, y = list(zip(*path))\n",
    "    _, (x0, y0) = starting_city\n",
    "    plt.imshow(map)\n",
    "    plt.plot(x0, y0, 'y*', markersize=15)  # y* = yellow star for starting point\n",
    "    plt.plot(x + x[:1], y + y[:1])  # include the starting point at the end of path\n",
    "    plt.axis(\"off\")\n",
    "    fig = plt.gcf()\n",
    "    fig.set_size_inches([w, h])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## II. Simulated Annealing -- Main Loop\n",
    "\n",
    "The main loop of simulated annealing repeatedly generates successors in the neighborhood of the current state and considers moving there according to an acceptance probability distribution parameterized by a cooling schedule.  See the [simulated-annealing function](https://github.com/aimacode/aima-pseudocode/blob/master/md/Simulated-Annealing.md) pseudocode from the AIMA textbook online at github.  Note that our Problem class is already a \"node\", so the MAKE-NODE line is not required."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def simulated_annealing(problem, schedule):\n",
    "    \"\"\"The simulated annealing algorithm, a version of stochastic hill climbing\n",
    "    where some downhill moves are allowed. Downhill moves are accepted readily\n",
    "    early in the annealing schedule and then less often as time goes on. The\n",
    "    schedule input determines the value of the temperature T as a function of\n",
    "    time. [Norvig, AIMA Chapter 3]\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    problem : Problem\n",
    "        An optimization problem, already initialized to a random starting state.\n",
    "        The Problem class interface must implement a callable method\n",
    "        \"successors()\" which returns states in the neighborhood of the current\n",
    "        state, and a callable function \"get_value()\" which returns a fitness\n",
    "        score for the state. (See the `TravelingSalesmanProblem` class below\n",
    "        for details.)\n",
    "\n",
    "    schedule : callable\n",
    "        A function mapping time to \"temperature\". \"Time\" is equivalent in this\n",
    "        case to the number of loop iterations.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    Problem\n",
    "        An approximate solution state of the optimization problem\n",
    "        \n",
    "    Notes\n",
    "    -----\n",
    "        (1) DO NOT include the MAKE-NODE line from the AIMA pseudocode\n",
    "\n",
    "        (2) Modify the termination condition to return when the temperature\n",
    "        falls below some reasonable minimum value (e.g., 1e-10) rather than\n",
    "        testing for exact equality to zero\n",
    "        \n",
    "    See Also\n",
    "    --------\n",
    "    AIMA simulated_annealing() pseudocode\n",
    "        https://github.com/aimacode/aima-pseudocode/blob/master/md/Simulated-Annealing.md\n",
    "    \"\"\"     \n",
    "    raise NotImplementedError\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## III. Representing the Problem\n",
    "\n",
    "In order to use simulated annealing we need to build a representation of the problem domain.  The choice of representation can have a significant impact on the performance of simulated annealing and other optimization techniques.  Since the TSP deals with a close loop that visits each city in a list once, we will represent each city by a tuple containing the city name and its position specified by an (x,y) location on a grid.  The _state_ will then consist of an ordered sequence (a list) of the cities; the path is defined as the sequence generated by traveling from each city in the list to the next in order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class TravelingSalesmanProblem:\n",
    "    \"\"\"Representation of a traveling salesman optimization problem.  The goal\n",
    "    is to find the shortest path that visits every city in a closed loop path.\n",
    "    \n",
    "    Students should only need to implement or modify the successors() and\n",
    "    get_values() methods.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    cities : list\n",
    "        A list of cities specified by a tuple containing the name and the x, y\n",
    "        location of the city on a grid. e.g., (\"Atlanta\", (585.6, 376.8))\n",
    "    \n",
    "    Attributes\n",
    "    ----------\n",
    "    names\n",
    "    coords\n",
    "    path : list\n",
    "        The current path between cities as specified by the order of the city\n",
    "        tuples in the list.\n",
    "    \"\"\"\n",
    "    def __init__(self, cities):\n",
    "        self.path = copy.deepcopy(cities)\n",
    "    \n",
    "    def copy(self):\n",
    "        \"\"\"Return a copy of the current board state.\"\"\"\n",
    "        new_tsp = TravelingSalesmanProblem(self.path)\n",
    "        return new_tsp\n",
    "    \n",
    "    @property\n",
    "    def names(self):\n",
    "        \"\"\"Strip and return only the city name from each element of the\n",
    "        path list. For example,\n",
    "            [(\"Atlanta\", (585.6, 376.8)), ...] -> [\"Atlanta\", ...]\n",
    "        \"\"\"\n",
    "        names, _ = zip(*self.path)\n",
    "        return names\n",
    "    \n",
    "    @property\n",
    "    def coords(self):\n",
    "        \"\"\"Strip the city name from each element of the path list and return\n",
    "        a list of tuples containing only pairs of xy coordinates for the\n",
    "        cities. For example,\n",
    "            [(\"Atlanta\", (585.6, 376.8)), ...] -> [(585.6, 376.8), ...]\n",
    "        \"\"\"\n",
    "        _, coords = zip(*self.path)\n",
    "        return coords\n",
    "    \n",
    "    def successors(self):\n",
    "        \"\"\"Return a list of states in the neighborhood of the current state by\n",
    "        switching the order in which any adjacent pair of cities is visited.\n",
    "        \n",
    "        For example, if the current list of cities (i.e., the path) is [A, B, C, D]\n",
    "        then the neighbors will include [A, B, D, C], [A, C, B, D], [B, A, C, D],\n",
    "        and [D, B, C, A]. (The order of successors does not matter.)\n",
    "        \n",
    "        In general, a path of N cities will have N neighbors (note that path wraps\n",
    "        around the end of the list between the first and last cities).\n",
    "        \"\"\"\n",
    "        raise NotImplementedError\n",
    "\n",
    "    def get_value(self):\n",
    "        \"\"\"Calculate the total length of the closed-circuit path of the current\n",
    "        state by summing the distance between every pair of adjacent cities.  Since\n",
    "        the default simulated annealing algorithm seeks to maximize the objective\n",
    "        function, return -1x the path length. (Multiplying by -1 makes the smallest\n",
    "        path the smallest negative number, which is the maximum value.)\n",
    "        \n",
    "        Notes\n",
    "        -----\n",
    "            (1) Remember to include the edge from the last city back to the\n",
    "            first city\n",
    "            \n",
    "            (2) Remember to multiply the path length by -1 so that simulated\n",
    "            annealing finds the shortest path\n",
    "        \"\"\"\n",
    "        raise NotImplementedError"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Testing TravelingSalesmanProblem\n",
    "The following tests should validate the class constructor and functionality of the `successors()` and `get_value()` methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Construct an instance of the TravelingSalesmanProblem\n",
    "test_cities = [('DC', (11, 1)), ('SF', (0, 0)), ('PHX', (2, -3)), ('LA', (0, -4))]\n",
    "tsp = TravelingSalesmanProblem(test_cities)\n",
    "assert(tsp.path == test_cities)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Test the successors() method -- no output means the test passed\n",
    "successor_paths = [x.path for x in tsp.successors()]\n",
    "assert(all(x in [[('LA', (0, -4)), ('SF', (0, 0)), ('PHX', (2, -3)), ('DC', (11, 1))],\n",
    "                 [('SF', (0, 0)), ('DC', (11, 1)), ('PHX', (2, -3)), ('LA', (0, -4))],\n",
    "                 [('DC', (11, 1)), ('PHX', (2, -3)), ('SF', (0, 0)), ('LA', (0, -4))],\n",
    "                 [('DC', (11, 1)), ('SF', (0, 0)), ('LA', (0, -4)), ('PHX', (2, -3))]]\n",
    "          for x in successor_paths))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Test the get_value() method -- no output means the test passed\n",
    "assert(np.allclose(tsp.get_value(), -28.97, atol=1e-3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IV. Define the Temperature Schedule\n",
    "\n",
    "The most common temperature schedule is simple exponential decay:\n",
    "$T(t) = \\alpha^t T_0$\n",
    "\n",
    "(Note that this is equivalent to the incremental form $T_{i+1} = \\alpha T_i$, but implementing that form is slightly more complicated because you need to preserve state between calls.)\n",
    "\n",
    "In most cases, the valid range for temperature $T_0$ can be very high (e.g., 1e8 or higher), and the _decay parameter_ $\\alpha$ should be close to, but less than 1.0 (e.g., 0.95 or 0.99).  Think about the ways these parameters effect the simulated annealing function.  Try experimenting with both parameters to see how it changes runtime and the quality of solutions.\n",
    "\n",
    "You can also experiment with other schedule functions -- linear, quadratic, etc.  Think about the ways that changing the form of the temperature schedule changes the behavior and results of the simulated annealing function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# These are presented as globals so that the signature of schedule()\n",
    "# matches what is shown in the AIMA textbook; you could alternatively\n",
    "# define them within the schedule function, use a closure to limit\n",
    "# their scope, or define an object if you would prefer not to use\n",
    "# global variables\n",
    "alpha = 0.95\n",
    "temperature=1e4\n",
    "\n",
    "def schedule(time):\n",
    "    raise NotImplementedError"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Testing the Temperature Schedule\n",
    "The following tests should validate the temperature schedule function and perform a simple test of the simulated annealing function to solve a small TSP test case"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# test the schedule() function -- no output means that the tests passed\n",
    "assert(np.allclose(alpha, 0.95, atol=1e-3))\n",
    "assert(np.allclose(schedule(0), temperature, atol=1e-3))\n",
    "assert(np.allclose(schedule(10), 5987.3694, atol=1e-3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Failure implies that the initial path of the test case has been changed\n",
    "assert(tsp.path == [('DC', (11, 1)), ('SF', (0, 0)), ('PHX', (2, -3)), ('LA', (0, -4))])\n",
    "result = simulated_annealing(tsp, schedule)\n",
    "print(\"Initial score: {}\\nStarting Path: {!s}\".format(tsp.get_value(), tsp.path))\n",
    "print(\"Final score: {}\\nFinal Path: {!s}\".format(result.get_value(), result.path))\n",
    "assert(tsp.path != result.path)\n",
    "assert(result.get_value() > tsp.get_value())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## V. Run Simulated Annealing on a Larger TSP\n",
    "Now we are ready to solve a TSP on a bigger problem instance by finding a shortest-path circuit through several of the US state capitals.\n",
    "\n",
    "You can increase the `num_cities` parameter up to 30 to experiment with increasingly larger domains.  Try running the solver repeatedly -- how stable are the results?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Create the problem instance and plot the initial state\n",
    "num_cities = 10\n",
    "capitals_tsp = TravelingSalesmanProblem(capitals_list[:num_cities])\n",
    "starting_city = capitals_list[0]\n",
    "print(\"Initial path value: {:.2f}\".format(-capitals_tsp.get_value()))\n",
    "print(capitals_list[:num_cities])  # The start/end point is indicated with a yellow star\n",
    "show_path(capitals_tsp.coords, starting_city)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# set the decay rate and initial temperature parameters, then run simulated annealing to solve the TSP\n",
    "alpha = 0.95\n",
    "temperature=1e6\n",
    "result = simulated_annealing(capitals_tsp, schedule)\n",
    "print(\"Final path length: {:.2f}\".format(-result.get_value()))\n",
    "print(result.path)\n",
    "show_path(result.coords, starting_city)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Experiments (Optional)\n",
    "Here are some ideas for additional experiments with various settings and parameters once you've completed the lab.\n",
    "\n",
    "- Change the number of cities in the final map (between 10 and 30).  How are your results affected?  Why?\n",
    "- Change the alpha and temperature parameters.  How do they affect the results?\n",
    "- Use a different schedule function (something other than exponential decay).  Is the algorithm still effective?\n",
    "- Use a different successors function; e.g., generate successors of a state by swapping _any_ pair of cities in the path, rather than only adjacent cities.  Try defining your own successor function.  What effect does the change have?\n",
    "- Use a different distance metric for get_value (e.g., we used the L2-norm (Euclidean distance), try the L1-norm (manhattan distance) or L$\\infty$-norm (uniform norm)\n",
    "\n",
    "Share and discuss your results with others in the forums!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:aind]",
   "language": "python",
   "name": "conda-env-aind-py"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}