Lab9_Simulated-Annealing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "np.random.seed(2017)  # for reproducibility\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Προσομοιωμένη Ανόπτηση\n",
    "\n",
    "Στην άσκηση αυτή θα χρησιμοποιήσετε τη μέθοδο της προσομοιωμένης ανόπτησης\n",
    "([simulated annealing](https://en.wikipedia.org/wiki/Simulated_annealing)) για\n",
    "να βρείτε το ελάχιστο μιας συνάρτησης $V: \\mathbb{R} \\to \\mathbb{R}$. \n",
    "Σαν παράδειγμα ξεκινάμε με την πολυωνυμική συνάρτηση:\n",
    "\n",
    "$$ V(x) = \\frac{x^4}{4} + \\frac{7 x^3}{15} + \\frac{4 x^2}{5} + \\frac{4x}{5} + 2.$$\n",
    "\n",
    "Η γραφική της παράσταση κατασκευάζεται στο παρακάτω κελί.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def V(x):\n",
    "    \"\"\"Energy (objective) function to minimize\"\"\"\n",
    "    return np.polyval((0.25, 7/15, -0.8, -0.8, 2.0), x)\n",
    "\n",
    "x = np.linspace(-3, 2, 100)\n",
    "\n",
    "plt.figure(figsize = (9, 6))\n",
    "plt.plot(x, V(x))\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"V(x)\")\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Όπως παρατηρείτε η συνάρτηση έχει ένα ολικό ελάχιστο στο $x= -2$ και ένα τοπικό ελάχιστο\n",
    "στο $x = 1$.\n",
    "\n",
    "Ο αλγόριθμος της προσομοιωμένης ανόπτησης δουλεύει ως εξής. Αρχικά ένα \"σωματίδιο\" είναι\n",
    "τοποθετημένο τυχαία στο διάστημα $(-3.5, +3.5)$ και έχει θερμοκρασία $T = 128$. \n",
    "Σε κάθε επανάληψη επιχειρεί μια μετατόπιση που ακολουθεί ομοιόμορφη κατανομή στο διάστημα \n",
    "$(−\\delta, +\\delta)$ και πραγματοποιεί τη μετατόπιση σύμφωνα με τον αλγόριθμο *Metropolis-Hastings*.\n",
    "Όπως και στο προηγούμενο εργαστήριο, η θερμοκρασία ελέγχει πόσο \"ευέλικτο\" είναι το σωματίδιο μας,\n",
    "δηλαδή πόσο εύκολα αποδέχεται μετακινήσεις $x\\to y$ που φαίνονται χειρότερες ($V(y)>V(x)$) από την\n",
    "παρούσα θέση. \n",
    "\n",
    "Αρχικά αφήνουμε το σωματίδιο να κινείται ελεύθερα και να εξερευνεί το χώρο (μεγάλη θερμοκρασία)\n",
    "αλλά στη συνέχεια, κάθε $n$ βήματα, το περιορίζουμε (ρίχνουμε τη θερμοκρασία). \n",
    "Η ιδέα είναι ότι αφού δώσουμε στο σωματίδιο την ευκαιρία να εξερευνήσει το χώρο για κάποιο διάστημα, καθώς το ψύχουμε\n",
    "αυτό θα προτιμάει να επισκέπτεται κοιλάδες παρά λόφους. Σκοπός μας είναι σιγά σιγά να παγιδεύσουμε\n",
    "το σωματίδιο στην βαθύτερη κοιλάδα (κοντά στο ελάχιστο της συνάρτησης $V$.) Η διαδικασία αυτή (της σταδιακής ψύξης) συνεχίζεται μέχρι η θερμοκρασία να φτάσει στην $T_{min}$\n",
    "\n",
    "*Σημείωση*: Η συνάρτηση `np.polyval(p, x)` παίρνει σαν είσοδο ένα διάνυμα `p` και ένα μία τιμή `x`\n",
    "και υπολογίζει το πολυώνυμο $ p \\cdot {(x^d, x^{d-1}, ..., x^0)}^T$ όπου `d` είναι η τάξη του `p`. Όταν το `x`\n",
    "είναι διάνυσμα (όπως στην περίπτωση μας), υπολογίζει την τιμή πολυωνύμου στο σημείο `x`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Αλγόριθμος\n",
    "\n",
    "Στο επόμενο κελί κώδικα σας δίνουμε μια απλή υλοποίηση της ιδέας αυτής."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def sim_anneal(f, x_low, x_high,\n",
    "               T=2.0, T_min=2**(-12), \n",
    "               cool_factor=0.9, cool_sch=100,\n",
    "               delta = 1.0, shrink_factor=0.5, shrink_thres=0.25,\n",
    "               trace = False):\n",
    "    \n",
    "    # Initialize\n",
    "    x = np.random.uniform(x_low, x_high)\n",
    "    f_cur = f(x)  # current value of f\n",
    "    steps = accepted = 0\n",
    "    if trace: path = []\n",
    "\n",
    "    # Main loop\n",
    "    while T > T_min:\n",
    "        \n",
    "        # proposed new position\n",
    "        x_new = x + np.random.uniform(-delta, delta)\n",
    "        if x_new > x_high or x_new < x_low:\n",
    "            continue  # skip this proposal\n",
    "        f_new = f(x_new)\n",
    "\n",
    "        # accept/reject jump\n",
    "        if np.log(np.random.rand()) < -(f_new - f_cur) / T :\n",
    "            x, f_cur = x_new, f_new\n",
    "            accepted += 1            \n",
    "                \n",
    "        # Cooling\n",
    "        if steps == cool_sch:  # time for the next cooling\n",
    "            T *= cool_factor  # cool down\n",
    "            # Adjust Search Space\n",
    "            if accepted / steps < shrink_thres:\n",
    "                delta *= shrink_factor\n",
    "            steps = accepted = 0\n",
    "        \n",
    "        if trace: path.append(x)\n",
    "        steps += 1        \n",
    "        \n",
    "    \n",
    "    if trace: return path\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Οι παράμετροι της συνάρτησης `sim_anneal` είναι:\n",
    "\n",
    "* **`f`**: η συνάρτηση που θέλουμε να ελαχιστοποιήσουμε\n",
    "* **`x_low, x_high`**: Το όρια μέσα στα οποία ψάχνουμε να βρούμε το ελάχιστο\n",
    "* **`T, T_min`**: Η αρχική και η τελική θερμοκρασία της \"ανόπτησης\"\n",
    "* **`cool_factor`**: Πόσο μειώνεται η θερμοκρασία σε κάθε κύκλο ψύξης (ποσοστό της προηγούμενης)\n",
    "* **`cool_sch`**: Κάθε πόσα βήματα να μειώνεται η θερμοκρασία\n",
    "* **`delta`**: Ακτίνα από την οποία θα προτείνουμε άλμα\n",
    "* **`shrink_factor`**: Πόσο θα μειώσουμε το `delta` αν δούμε ότι η μέθοδος έχει κολλήσει\n",
    "* **`shirnk_thres`**: Ποσοστό το επιτυχημένων αλμάτων κάτω από το οποίο θεωρούμε ότι το σωματίδιο έχει κολλήσει.\n",
    "* **`trace`**: Αν θέλουμε να μας γυρίσει όλη την διαδρομή ή μόνο την τελική θέση\n",
    "\n",
    "Επομένως για να βρούμε το ελάχιστο της συνάρτησης $V$ που ορίσαμε προηγουμένως στο\n",
    "διάστημα $(-3.5, 3.5)$ τρέχουμε την συνάρτηση με τις εξής παραμέτρους.\n",
    "\n",
    "*Παρατήρηση: Ζητάμε όλο το μονοπάτι του σωματιδίου για να το πλοτάρουμε αργότερα.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Parameters\n",
    "T_start, T_end = 2.0, 2**(-12)\n",
    "cool_sch, cool_factor = 100, 0.9\n",
    "\n",
    "path = sim_anneal(V, -3.5, 3.5, delta = 2.,\n",
    "                  T = T_start, T_min = T_end,\n",
    "                  cool_factor = cool_factor, \n",
    "                  cool_sch = cool_sch,\n",
    "                  trace = True)  # keep the default for the rest\n",
    "\n",
    "path[-1]  # last position"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Στο επόμενο κελί, πλοτάρουμε την θέση του σωματιδίου και την θερμοκρασία σαν συνάρτηση του χρόνου"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Compute Temperature\n",
    "Ncool = np.log2(T_end / T_start) / np.log2(cool_factor)\n",
    "Ncool = np.ceil(Ncool) + 1  # number of coolings\n",
    "T = T_start * cool_factor ** np.arange(Ncool)\n",
    "T = np.repeat(T, 2)[:-1]  # step function\n",
    "steps = np.arange(0, len(path) + 1, cool_sch)\n",
    "steps = np.repeat(steps, 2)[1:]  # step function\n",
    "\n",
    "# Plot\n",
    "fig, ax = plt.subplots()\n",
    "fig.set_size_inches((9, 6))\n",
    "\n",
    "# Plot trace\n",
    "ax.plot(path)\n",
    "ax.set_xlabel('Step')\n",
    "ax.set_ylabel('Current Position')\n",
    "\n",
    "# Plot temp\n",
    "ax2 = ax.twinx()\n",
    "ax2.plot(steps, T, 'r--')\n",
    "ax2.set_ylabel('Temperature', color='red')\n",
    "ax2.tick_params('y', colors='red')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Το επόμενο κελί κώδικα οπτικοποιεί τη διαδικασία."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib qt\n",
    "%run animate_simanneal.py -minima 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Το πρόγραμμα `animate_simanneal.py` μπορεί να δεχτεί για ορίσματα:\n",
    "\n",
    "- **`minima`**  επιλέγει ένα πολυώνυμο με {0,1,2,3} ακρότατα. \n",
    "Για μηδέν είναι το παράδειγμα που δόθηκε\n",
    "- **`T0`** Αρχική θερμοκρασία (default 128)Initial Temperature\n",
    "- **`T1`** Τελική θερμοκρασία (default 0.01)\n",
    "- **`CF`** Ποσοστό απώλειας θερμοκρασία σε κάθε βήμα (default 0.5)\n",
    "- **`CS`** Περίοδος ψύξης (default 100)\n",
    "- **`d`** ακτίνα που άλματος (default 3.0)\n",
    "- **`--delay`** χρόνος μεταξύ frames σε ms (default 100)\n",
    "\n",
    "Δοκιμάστε να δώσετε άλλα δεδομένα και δείτε πως αλλάζει η συμπεριφορά της ανόπτησης.\n",
    "\n",
    "Αν η εντολή σας βγάλει σφάλμα στο notebook δοκιμάστε να την τρέξετε\n",
    "από το terminal ως\n",
    "\n",
    "`python animate_simanneal.py -minima 0`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Τρέχτε το παρακάτω κελί ώστε τα διαγράμματα να εμφανίζονται μέσα στο notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Παραδοτέο\n",
    "\n",
    "Σε ένα καινούργιο jupyter notebook, χρησιμοποιήστε την μέθοδο της προσομοιωμένης ανόπτησης για την συνάρτηση:\n",
    "\n",
    "$$ V(x) = \\frac{13}{90} x^4 - \\frac{3}{10} x^3 - \\frac{101}{45} x^2 \\frac{13}{15}x + 4$$\n",
    "\n",
    "Α) Τρέξτε τον αλγόριθμο `200` φορές και υπολογίστε το ποσοστό των πειραμάτων στα οποία \n",
    "ο αλγόριθμος βρήκε το ολικό ελάχιστο και το ποσοστό στο οποίο βρήκε κάποιο τοπικό ελάχιστο\n",
    "(με ακρίβεια 0.01)\n",
    "\n",
    "Β) Τρέχτε τον αλγόριθμο `200` φορές για τιμές του `cool_factor = [0.25, 0.5 0.75, 0.95]` και εκτιμήστε για κάθε τιμή του \n",
    "την πιθανότητα ο αλγόριθμος να υπολογίσει το ολικό ελάχιστο. \n",
    "Κάντε διάγραμμα των τιμών του `cool_factor` και της πιθανότητας εύρεση του ολικού ελάχιστου."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Λύση (μισή)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Evaluation Parameters\n",
    "global_min = -2   # the actual optimal value\n",
    "local_min = 1     # the local minimum\n",
    "precision = 0.1   # the desired precision\n",
    "n_exper = 500     # number of experiments\n",
    "n_correct = 0     # number of right guesses\n",
    "n_incorrect = 0   # number of wrong guess (local minimum)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n_correct = n_incorrect = 0\n",
    "print(\"# Iterations: \", end=\"\")\n",
    "for n in range(n_exper):\n",
    "    if n % 100 == 0: print(\"%d\" % n, end=\"...\")\n",
    "    x = sim_anneal(V, -3.5, 3.5, delta = 2.0)\n",
    "    if np.isclose(x, global_min, rtol=0., atol=precision):\n",
    "        n_correct += 1\n",
    "    elif np.isclose(x, local_min, rtol=0., atol=precision):\n",
    "        n_incorrect += 1\n",
    "print(\"Done!\\n\")\n",
    "print(\"The Algorithm found the global minimum %.2f%% of the times\" % (100 * n_correct / n_exper))\n",
    "print(\"The Algorithm found the local minimum %.2f%% of the times\" % (100 * n_incorrect / n_exper))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Πλανόδιος Πωλητής\n",
    "\n",
    "Στην άσκηση αυτή θα χρησιμοποιήσουμε τον αλγόριθμο της προσομοιωμένης ανόπτησης\n",
    "για να λύσουμε το πρόβλημα του [πλανόδιου πωλητή](https://en.wikipedia.org/wiki/Travelling_salesman_problem).\n",
    "\n",
    "Το πρόβλημα έχει ως εξής:\n",
    "\n",
    "> Ένας πωλητής θέλει να επισκεφτεί $Ν$ πόλεις και να γυρίσει στη βάση του.\n",
    "Ποια είναι η συντομότερη διαδρομή που ενώνει όλες τις πόλεις;\n",
    "\n",
    "Στη συγκεκριμένη άσκηση θα ασχοληθούμε με πόλεις της Ευρώπης.\n",
    "Τρέξτε το παρακάτω κελί κώδικα. Θα πρέπει να σας κατεβάσει το αρχείο [eu.csv](http://math.ntua.gr/~loulakis/info/python_codes_files/eu.csv) που περιέχει \n",
    "το γεωργαφικό μήκος και πλάτος 27 ευρωπαϊκών πόλεων. Ο κώδικας οργανώνει την πληροφορία ως εξής:\n",
    "\n",
    "* **`city_coord`** ένα `dictionary` που για κάθε πόλη έχει το γεωγραφικό μήκος και πλάτος\n",
    "* **`Ncities`** ο αριθμός των πόλεων (27)\n",
    "* **`city_names`** τα ονόματα των πόλεων\n",
    "\n",
    "Το κελί τυπώνει τις συντεταγμένες της Αθήνας για επαλήθευση."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import csv\n",
    "import urllib.request\n",
    "\n",
    "url = \"http://math.ntua.gr/~loulakis/info/python_codes_files/eu.csv\"\n",
    "with urllib.request.urlopen(url) as urlstream:\n",
    "    csvfile = csv.reader(urlstream.read().decode('utf-8').splitlines(), delimiter=\",\")\n",
    "    city_coord = {city: tuple(map(float, (longitude, latitude))) \n",
    "                  for longitude, latitude, city, _ in csvfile}\n",
    "\n",
    "Ncities = len(city_coord)  # number of cities\n",
    "city_names = list(city_coord.keys())\n",
    "city_coord[\"Athens\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Στο επόμενο κελί ορίζουμε την συνάρτηση `geodesic_dist` η οποία υπολογίζει την\n",
    "απόσταση μεταξύ 2 πόλεων δοθέντων των συντεταγμένων τους (credit to [Michael Dunn](http://stackoverflow.com/a/4913653/2167530)) με βάση τον τύπο [Haversine](https://en.wikipedia.org/wiki/Haversine_formula)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def geodesic_dist(city1, city2):\n",
    "    \"\"\"\n",
    "    Calculate the great circle distance between two points \n",
    "    on the earth (specified in decimal degrees)\n",
    "    \"\"\"\n",
    "    if city1 == city2:\n",
    "        return 0.\n",
    "    lon1, lat1 = np.deg2rad(city1)\n",
    "    lon2, lat2 = np.deg2rad(city2)\n",
    "    \n",
    "    # haversine formula \n",
    "    dlon = lon2 - lon1 \n",
    "    dlat = lat2 - lat1 \n",
    "    a = np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2\n",
    "    c = 2 * np.arcsin(np.sqrt(a)) \n",
    "    r = 6371 # Radius of earth in kilometers 6371. Use 3956 for miles\n",
    "    return c * r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Για ευκολία και ταχύτητα θα υπολογίσουμε όλες της αποστάσεις εκ των προτέρων και θα τις αποθηκεύσουμε σε ένα πίνακα"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Distance between every pair of cities\n",
    "dist = np.array([\n",
    "    [geodesic_dist(coord1, coord2) for coord2 in city_coord.values()]\n",
    "    for coord1 in city_coord.values()\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Σαν σύμβαση θεωρούμε ότι η αρχική σειρά των πόλεων είναι αυτή με την οποία\n",
    "μπορούμε να κάνουμε iterate το `city_coord` δηλαδή:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for i, city in enumerate(city_coord): print(i, city)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Αλγόριθμος\n",
    "\n",
    "Επειδή ο δρόμος που θα ακολουθήσει ο πωλητής είναι κυκλικός, δεν έχει σημασία\n",
    "από που θα ξεκινήσει (μόνο ποιες πόλεις θέλει να επισκεφθεί έχει σημασία).\n",
    "Ως εκ τούτου, \"λύση\" στο πρόβλημα μας είναι μια λίστα με τη σειρά με την οποία\n",
    "κάποιος πρέπει να επισκεφθεί τις πόλεις, δεδομένου ότι μετά την τελευταία πόλη\n",
    "θα επιστρέψει στην αρχική χωρίς να επισκεφθεί άλλη πόλη.\n",
    "\n",
    "Επομένως, σκοπός μας είναι να χρησιμοποιήσουμε την τεχνική της προσομοιωμένης ανόπτησης για να\n",
    "\"ψάξουμε\" τον χώρο των μεταθέσεων της λίστας των πόλεων ώστε να\n",
    "βρούμε την μετάθεση που δίνει την συνολική απόσταση κοντά στην ελάχιστη δυνατή.\n",
    "\n",
    "Για να τρέξουμε τον αλγόριθμο της ανόπτησης χρειαζόμαστε 2 στοιχεία:\n",
    "\n",
    "1. Μια συνάρτηση *ενέργειας*\n",
    "2. Έναν τρόπο να μεταβαίνουμε σε γειτονικές καταστάσεις\n",
    "\n",
    "Ως συνάρτηση ενέργειας θα πάρουμε την συνολική απόσταση που διανύει ο πωλητής \n",
    "(αυτή που θέλουμε να ελαχιστοποιήσουμε). \n",
    "Στο παρακάτω κελί ορίζουμε την συνάρτηση (είσοδος η σειρά με την οποία ο πωλητής επισκέπτεται τις πόλεις):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def total_dist(cities):\n",
    "    \"\"\"The total distance covered by the salesman\"\"\"\n",
    "    total = sum(dist[city1][city2] for city1, city2 in zip(cities[:-1], cities[1:]))\n",
    "    total += dist[cities[-1]][cities[0]] # return to start\n",
    "    return total"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ως γειτονικές μεταθέσεις, ορίζουμε τις μεταθέσεις που διαφέρουν μόνο σε ένα ζευγάρι πόλεων. \n",
    "Επομένως για να κινηθούμε στον χώρο το μεταθέσεων, θα διαλέγουμε τυχαία 2 πόλεις και θα αλλάζουμε\n",
    "τη θέση τους στη λίστα.\n",
    "\n",
    "Για να υπολογίσουμε την αλλαγή στην συνολική απόσταση για μια τέτοια αλλαγή,\n",
    "μπορούμε να προσθέσουμε όλες τις αποστασεις από την αρχή αλλά προφανώς κάτι τέτοιο \n",
    "δεν είναι αποδοτικό. Γι' αυτό το λόγο ορίζουμε την συνάρτηση `dist_diff`.\n",
    "Στην επόμενη ενότητα εξηγούμε την λογική της συνάρτησης, αν δεν την καταλαβαίνετε\n",
    "μπορείτε απλά να την χρησιμοποιήσετε σαν μαύρο κουτί.\n",
    "\n",
    "### Κόστος αλλαγής σειράς επίσκεψης\n",
    "\n",
    "Η αλλαγή της σειράς 2 πόλεων είναι μια τοπική αλλαγή επομένως αρκεί να σκεφτούμε ποιές είναι\n",
    "συνέπειες αυτές της αλλαγής. *Πόσες αποστάσεις επηρεάζονται συνολικά;*\n",
    "\n",
    "Θεωρείστε 2 τριάδες πόλεων τυχαία μέσα στη λίστα:\n",
    "\n",
    "* $A \\to B \\to C$ και \n",
    "* $X \\to Y \\to Z$\n",
    "\n",
    "και και χωρίς βλάβη στη γενικότητα θεωρείστε ότι η πόλη $B$ εμφανίζεται\n",
    "πριν από την $Y$ στη λίστα.\n",
    "\n",
    "Υπάρχουν 2 περιπτώσεις\n",
    "\n",
    "1. Γενική περίπτωση $B$ και $Y$ δεν συνορεύουν (δηλαδή $X \\not\\equiv B$)\n",
    "2. Ειδική περιπτωση $B$ και $Y$ συνορεύουν (δηλαδή $X \\equiv B$)\n",
    "\n",
    "Στο παρακάτω διάγραμμα φαίνονται οι συνέπειες που θα είχε μια αλλαγή μεταξύ $B$ και $Y$\n",
    "στις 2 αυτές περιπτώσεις\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![TSP Swap Cases](http://math.ntua.gr/~loulakis/info/python_codes_files/tsp_swap.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Στο παρακάτω κελί ορίζουμε την συνάρτηση."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def dist_diff(ABC, XYZ, tour):\n",
    "    \"\"\"\n",
    "    Differnce of total distance traveled if we swap cities B, Y.\n",
    "    Input: 2 triplets of nodes A, B, C and X, Y, Z.\n",
    "    \"\"\"\n",
    "    # order\n",
    "    A, B, C = ABC\n",
    "    X, Y, Z = XYZ\n",
    "    # cities\n",
    "    a, b, c = tour[ABC]\n",
    "    x, y, z = tour[XYZ]\n",
    "    \n",
    "    BpY = (B - Y) % len(tour) # to understand the order\n",
    "    dx = 0\n",
    "    \n",
    "    if BpY != -1 % len(tour):  #  B is not X\n",
    "        dx += dist[x, b] - dist[x, y] + dist[y, c] - dist[b, c]\n",
    "    if BpY != 1:  # Y is not A\n",
    "        dx += dist[a, y] - dist[a, b] + dist[b, z] - dist[y, z]\n",
    "\n",
    "    return dx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Υλοποίηση\n",
    "\n",
    "Στο παρακάτω κελί υλοποιούμε τον αλγόριθμο της προσομοιωμένης\n",
    "ανόπτοσηση για το πρόβλημα το πλανόδιου πωλητή.\n",
    "\n",
    "Αρχικά θέτουμε τις παραμέτρους του αλγορίθμου και μερικές μεταβλητές γαι να οπτικοποιήσουμε τα αποτελέσματα. Ο κώδικας είναι σχεδόν ίδιος με την συνάρτηση που δόθηκε προηγουμένως. Μια σημαντική διαφορά είναι ότι τώρα αποθηκεύουμε τη βέλτιστη λύση που έχουμε βρει μέχρι κάθε στιγμή. Αυτό το κάνουμε γιατί με μικρό υπολογιστικό κόστος μπορεί να βρούμε μια καλύτερη διαδρομή την οποία ο αλγόριθμος εξερεύνησε στην πορεία αλλά δεν τερμάτισε σε αυτήν.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Annealing Parameters\n",
    "Temp = 10 ** 7\n",
    "Temp_end = 0.02\n",
    "cool_factor = 0.95\n",
    "cool_schedule = 100\n",
    "\n",
    "# Bookeeping parameters\n",
    "tour = np.arange(Ncities)  # initial permutation\n",
    "tour_length = total_dist(tour)  # initial energy level\n",
    "\n",
    "shortest_tour = np.copy(tour)  # optimal permutation (to be updated)\n",
    "shortest_tour_length = tour_length  # shortest distance\n",
    "\n",
    "# To plot later\n",
    "first_guess = tour_length\n",
    "running_lengths = []\n",
    "optimal_lengths = []\n",
    "\n",
    "# Main loop\n",
    "while Temp > Temp_end:\n",
    "    # Constant Temperature Era\n",
    "    for step in range(cool_schedule):\n",
    "        # Choose 2 cities to swap\n",
    "        B, Y = np.random.choice(Ncities, size=2, replace=False)  # index of cities to swap\n",
    "        ABC = np.array([B-1, B, B+1]) % Ncities  # neighbors of B\n",
    "        XYZ = np.array([Y-1, Y, Y+1]) % Ncities  # neighbors of Y\n",
    "\n",
    "        # Compute the energy diffence  \n",
    "        diff = dist_diff(ABC, XYZ, tour)\n",
    "        \n",
    "        # Decide whether to change state.\n",
    "        if np.log(np.random.rand()) < -diff / Temp:\n",
    "            tour[B], tour[Y] = tour[Y], tour[B]\n",
    "            tour_length += diff\n",
    "            # Keep track of the optimal tour observed\n",
    "            if tour_length < shortest_tour_length:\n",
    "                shortest_tour = np.copy(tour)\n",
    "                shortest_tour_length = tour_length\n",
    "\n",
    "        running_lengths.append(tour_length)\n",
    "        optimal_lengths.append(shortest_tour_length)\n",
    "\n",
    "    Temp *= cool_factor  # cool down\n",
    "\n",
    "print('Improvement:\\t%8.0f%%'  % (100 * (first_guess - shortest_tour_length) / first_guess))\n",
    "print('Initial distance:\\t%.0f km'  % first_guess)\n",
    "print('Optimal distance:\\t%.0f km'  % shortest_tour_length)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Αποτελέσματα\n",
    "\n",
    "Παρακάτω πλοτάρουμε την εξέλιξη της τρέχουσας βέλτιστης λύσης κατά την διάρκεια του αλγορίθμου."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "running_lengths = np.array(running_lengths)\n",
    "optimal_lengths = np.array(optimal_lengths)\n",
    "\n",
    "plt.figure(figsize = (10, 7))\n",
    "plt.plot(np.arange(len(running_lengths)), running_lengths, label=\"candidate solution\")\n",
    "plt.plot(np.arange(len(optimal_lengths)), optimal_lengths, 'r-', linewidth=2, label=\"optimal solution\")\n",
    "\n",
    "plt.xlabel(\"Steps (Temperature)\")\n",
    "plt.ylabel(\"Tour Length\")\n",
    "plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "plt.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Τέλος, με το επόμενο κελί κώδικα δημιουργούμε μια οπτική αναπαράσταση \n",
    "της βέλτιστης διαδρομής.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 9))  # increase image size\n",
    "\n",
    "# Get Tour Coordinates\n",
    "coord = np.array([city_coord[city] for city in city_names])\n",
    "cities = shortest_tour.copy()\n",
    "cities = np.append(cities, cities[0])\n",
    "coord = coord[cities,:]  # rearrange\n",
    "\n",
    "# Scale for better readability\n",
    "scale_factor = 100\n",
    "coord = coord * scale_factor\n",
    "\n",
    "# Plot Tour\n",
    "plt.plot(coord[:, 1], coord[:,0], 'o-', mfc='r')\n",
    "plt.xlabel(\"Longitude\")\n",
    "plt.ylabel(\"Latitude\")\n",
    "plt.grid()\n",
    "\n",
    "# Manage spacing\n",
    "spacing = np.fabs(min(coord[:,1]) - max(coord[:,1])) * .1\n",
    "plt.xlim(min(coord[:,1]) - spacing, max(coord[:,1]) + spacing * 3)\n",
    "plt.ylim(min(coord[:,0]) - spacing, max(coord[:,0]) + spacing)\n",
    "\n",
    "# Add labels - remove for more clear picture\n",
    "if len(cities)<=30:\n",
    "    for i, city in enumerate(cities[:-1]):\n",
    "        plt.text(coord[i,1], coord[i,0],\n",
    "                 '%s' % city_names[city],\n",
    "                 withdash = True,\n",
    "                 dashdirection = 1,\n",
    "                 dashlength = 30,\n",
    "                 rotation = 0,\n",
    "                 dashrotation = 20,\n",
    "                 dashpush = 10)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Παραδοτέο\n",
    "\n",
    "Σε συνέχεια του προηγούμενο notebook:\n",
    "\n",
    "Α) Δοκιμάστε να βελτιώσετε την λύση πειράζοντας της παραμέτρους του κώδικα.\n",
    "\n",
    "B) Τρέχτε τον κώδικα του πλανόδιου πωλητή για 35 πόλεις τις Ευρώπης.\n",
    "(Τα δεδομένα είναι σε αυτή τη διεύθυνσης\n",
    "\n",
    "`http://math.ntua.gr/~loulakis/info/python_codes_files/europe.csv`\n",
    "\n",
    "Αλλάχτε το `url` στο αντίστοιχο κελί κώδικα απλά.\n",
    "\n",
    "Και για τα 2 υποερωτήματα, πλοτάρετε τα 2 γραφήματα (εξέλιξη της τρέχουσας βέλτιστης λύσης και προτεινόμενη διαδρομή) όπως σε αυτό το notebook.\n",
    "\n",
    "Γ) Σε ένα κελί κειμένου, απαντήστε στις παρακάτω ερωτήσεις:\n",
    "\n",
    "1. Πόσα είναι τα δυνατά μονοπάτια που περνούν από όλες τις 35 πόλεις μία μόνο φορά; \n",
    "\n",
    "2. Αν σε έναν υπολογιστή ο χρόνος που απαιτείται για τον υπολογισμό του συνολικού μήκους ενός τέτοιου μονοπατιού ήταν 1μsec, πόσο χρόνο θα χρειαζόμασταν για να βρούμε το βέλτιστο μονοπάτι με έναν εξαντλητικό αλγόριθμο (brute force) που θα υπολόγιζε τη συνολική απόσταση που θα διανύσει ο πλανόδιος πωλητής κατά μήκος κάθε μονοπατιού και θα επέστρεφε το μονοπάτι με την ελάχιστη συνολική απόσταση;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}