\n",
"\n",
"# The Traveling Salesperson Problem\n",
"\n",
"Consider the [***Traveling Salesperson Problem***](http://en.wikipedia.org/wiki/Traveling_salesman_problem) (abbreviated ***TSP***): \n",
"\n",
"- *Given a **set of cities** and the **distance** between each pair of cities, what is the **shortest** possible **tour** that visits each city exactly once, and returns to the starting city?*\n",
"\n",
"In this notebook we will develop some solutions to the problem, and more generally show *how to think about* solving problems. Versions of the algorithms developed here are used in [serious applications](https://research.googleblog.com/2016/09/the-280-year-old-algorithm-inside.html) that millions of people rely on every day. \n",
"\n",
"|![](http://support.sas.com/documentation/cdl/en/ornoaug/66084/HTML/default/images/map002g.png)|\n",
"|---|\n",
"|[An example tour, from the TSP History page](http://www.math.uwaterloo.ca/tsp/history/pictorial/dfj.html)|\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Implementation of Basic Concepts\n",
"\n",
"I'll go through all the concepts from the definition and make Python implementation choices for each one:\n",
"\n",
"- **Distance:** The straight-line distance between two points in a two-dimensional plane (or between two cities on a map) is called the [**Euclidean distance**](http://en.wikipedia.org/wiki/Euclidean_distance). I'll implement that with a function `distance(A, B)`. This distance measure is **non-negative** and **symmetric** (the distance from `A` to `B` is the same as the distance from `B` to `A`). But not all distance metrics are like this. If you have to follow roads rather than straight lines, the distance is at the whim of the road-builders. Perhaps the distance from `A` to `B` is different than the distance from `B` to `A` because of one-way streets. Or perhaps you are taking plane flights and the total time of going from `A` to `B` to `C` is less than going directly from `A` to `C` (because `B` is a hub with frequent flights). Almost everything in this notebook still holds with a different distance function; I'll point out where it doesn't.\n",
"- **City:** Given the choice of a Euclidean distance function, the only thing we need to know about a city is its position on the (*x*, *y*) plane. We don't need to know the city's name, population, best restaurants, or anything else. I will define `City` so that `City(300, 100)` creates a point with *x*-coordinate 300 and *y*-coordinate 100. \n",
"- **Set of cities:** A set of cities can be represented with a Python set. I'll use [`frozenset`](https://docs.python.org/3/library/stdtypes.html?highlight=frozenset#frozenset), which is a set that can't be changed. \n",
"- **Tour** (also called a **cycle** or **circuit** or [**Hamiltonian path**](https://en.wikipedia.org/wiki/Hamiltonian_path)): A tour that goes from city `A` to `B` to`C` and back to `A` will be represented by the list `[A, B, C]`. \n",
"- **Shortest**: The shortest tour is the one with the smallest total tour length. `shortest(tours)` implements this.\n",
"- **Tour length:** The call `tour_length(tour)` gives the sum of the distances between adjacent cities in the tour (including back to the start).\n",
"- **Valid Tour:** A tour is valid if it visits every city exactly once and returns to the start; in other words if it is a *permutation* of the cities. `valid_tour` implements that.\n",
"\n",
"Three more basic concepts, not explicitly mentioned in the definition:\n",
"\n",
"- **Point:** At first glance, Python does not have a builtin type for a two-dimensional point. But there is one: [complex numbers](https://docs.python.org/3/c-api/complex.html). \n",
"- **Link**: A tour consists of a sequence of **links** between cities. A link can be represented implicitly by the adjacency of two cities in a tour, or explicitly by an `(A, B)` pair of cities.\n",
"- **TSP algorithm**: A function with the signature `tsp(cities: Cities) -> Tour`.\n",
"\n",
"First some imports (don't worry about these; they will be explained later as they are used):"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import functools\n",
"import itertools\n",
"import pathlib\n",
"import random\n",
"import time \n",
"import math\n",
"import re\n",
"import matplotlib.pyplot as plt \n",
"from collections import Counter, defaultdict, namedtuple\n",
"from statistics import mean, median, stdev\n",
"from typing import Set, List, Tuple, Iterable, Dict"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now the implementation of the basic concepts:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"City = complex # e.g. City(300, 100)\n",
"Cities = frozenset # A set of cities\n",
"Tour = list # A list of cities visited, in order\n",
"TSP = callable # A TSP algorithm is a callable function\n",
"Link = Tuple[City, City] # A city-city link\n",
"\n",
"def distance(A: City, B: City) -> float: \n",
" \"Distance between two cities\"\n",
" return abs(A - B)\n",
"\n",
"def shortest(tours: Iterable[Tour]) -> Tour: \n",
" \"The tour with the smallest tour length.\"\n",
" return min(tours, key=tour_length)\n",
"\n",
"def tour_length(tour: Tour) -> float:\n",
" \"The total distances of each link in the tour, including the link from last back to first.\"\n",
" return sum(distance(tour[i], tour[i - 1]) for i in range(len(tour)))\n",
"\n",
"def valid_tour(tour: Tour, cities: Cities) -> bool:\n",
" \"Does `tour` visit every city in `cities` exactly once?\"\n",
" return Counter(tour) == Counter(cities)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Sets of Random Cities\n",
"\n",
"To test TSP algorithms, I'll need sets of cities. I'll define `random_cities(n)` to return a set of `n` cities, sprawled out randomly over the map. \n",
"\n",
"*Note*: I want to be able to do reproducible comparisons of algorithms: to run two or more algorithms on the same set of cities. Therefore `random_cities` has an optional argument, `seed`, with a default value. If you pass in the same seed, you will get back the same set of cities. This is true even after restarting with a different version of Python (it could possibly change with a major revision, as in Python 2 to Python 3). If you want a different set of *n* random cities, pass in a different seed."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def random_cities(n, seed=1234, width=9999, height=6666) -> Cities:\n",
" \"Make a set of n cities, sampled uniformly from a (width x height) rectangle.\"\n",
" random.seed((n, seed)) # To make `random_cities` reproducible\n",
" return Cities(City(random.randrange(width), random.randrange(height))\n",
" for c in range(n))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Exhaustive TSP Search Algorithm: `exhaustive_tsp`"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's start with an algorithm that is *guaranteed* to find the shortest tour, although inefficiently:\n",
"\n",
"- **Exhaustive TSP Search Algorithm**: *Generate all possible tours of the cities, and choose the shortest one.*\n",
"\n",
"My design philosophy is to first write an English description of the algorithm (as above), then write Python code that closely mirrors the English description. I note that the possible tours of a set of cities are the permutations of the cities. So implementing this algorithm is easy. I note that the possible tours of a set of cities are just the permutations of the city, and `permutations` is defined in the `itertools` module, so we get:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def exhaustive_tsp(cities) -> Tour:\n",
" \"Generate all possible tours of the cities and choose the shortest one.\"\n",
" return shortest(possible_tours(cities))\n",
"\n",
"possible_tours = itertools.permutations "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's try it on a random set of 8 cities:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"((4698+4764j),\n",
" (2301+3350j),\n",
" (1276+3305j),\n",
" (4215+1920j),\n",
" (8883+324j),\n",
" (9744+950j),\n",
" (8315+4692j),\n",
" (9081+6286j))"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"exhaustive_tsp(random_cities(8))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Quick, is that the shortest tour? I can't tell, and I bet you can't either. But if I plotted the points on a map, maybe we would be enlightened.\n",
"\n",
"# Visualizing results: `plot_tour`\n",
"\n",
"I'll define `plot_tour` to plot all the cities and links in a tour, highlighting the first city. \n",
"\n",
"- *Vocabulary note:* A **segment** is a portion of a tour that does not loop back to the start. The **segment** `[A, B, C]` has only two links, A-B and B-C, whereas the **tour** `[A, B, C]` has three links, A-B, B-C, and C-A."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"Segment = list # A portion of a tour; it does not loop back to the start.\n",
"\n",
"def plot_tour(tour: Tour, style='bo-', hilite='rs', title=''): \n",
" \"Plot every city and link in the tour, and highlight the start city.\"\n",
" scale = 1 + len(tour) ** 0.5 // 10\n",
" plt.figure(figsize=((3 * scale, 2 * scale)))\n",
" start = tour[0]\n",
" plot_segment([*tour, start], style)\n",
" plot_segment([start], hilite) \n",
" plt.title(title)\n",
" \n",
"def Xs(cities) -> List[float]: \"X coordinates\"; return [c.real for c in cities]\n",
"def Ys(cities) -> List[float]: \"Y coordinates\"; return [c.imag for c in cities]\n",
"\n",
"def plot_segment(segment: Segment, style='bo:'):\n",
" \"Plot every city and link in the segment.\"\n",
" plt.plot(Xs(segment), Ys(segment), style, linewidth=2/3, markersize=4, clip_on=False)\n",
" plt.axis('scaled'); plt.axis('off')"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAO4AAACuCAYAAAAxtjFOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAUD0lEQVR4nO3deXRU5RkG8CeyBVlkCUE5ikIRcCGhCAWk7EuRTUAgSS1toeVQa7FaLQVsTekGQl2AHmxrbdWiMwEyRMQKKKtsIpRIBEKaJhIWZZM9hITJ7R9Pc4aEJEySmfnunXl+53AQl8x7jnny3fst7xdlWZYFEXGUm0wXICJVp+CKOJCCK+JACq6IAym4Ig6k4Io4kIIr4kAKrogDKbgiDqTgijiQgiviQAquiAMpuCIOpOCK3IDHA8THA/Xr83ePx3RFQJSO9YlUzOMBHnkEiIoCLMv3e2oqMHasuboUXJFKxMcDGRkMa4moKCAuDkhPN1aWgitSHssCdu8GevQAvN7r/3l0NHD5cujrKlHb3EeL2M+BA4DLBbjdQPPmQMuWwBdfXD/iduhgrkZAk1Mi+PxzYO5cPhYnJHASas0aYPt2YNGi60NrWUBysrFyWYcelSUSHT8OLF3K0fXECSAxEUhKAu67r/S/d+ECcOedwG23ATk5HGmTk4ExY8zUXULBlYhx9ixniV0u4OBBzhYnJQHdunEkLc+8ecCRI8DChSEt9YYUXAlrly4Bq1YxrDt3AiNHcnTt0weoVavy//byZY6wW7cCd9wRmnr9pckpCTuFhXxHdbmA9euBQYOAKVP4aFy3rv9f57XXgCFD7BdaQCOuhAmvF9i0iWF97z2ge3c+Bo8YAdx8c9W/XmEh0LEjsHYt0K5d4OutKY244liWxcdfl4s7me69l4/B8+cDTZrU7GsvWcI1XDuGFlBwxYEyMhjWlBTO9iYlcbNEbGxgvr7XCzz/PH8Y2JWCK47w3/9yU4TLBdSpw7CuX8+lmkBbtgy45x7g/vsD/7UDRe+4YlvHjvnWWs+e9a21duwYvM8sLgY6d+bEVLduwfucmtKIK7Zy+jQfUV0uIDcXGD8eeOUV4Otfr3itNZBWrQJuvdXeoQU04ooNXLwIvPMOw7pnD/DwwxxZe/UCbgrhplzL4oTUvHlA376h+9zq0IgrRhQUAKtXM6ybN3O99Cc/AQYO5DusCevW8bP79DHz+VWhEVdC5upVTii5XAxtr14cWYcN48Z+0/r3B6ZPBx56yHQlN6YRV4KquJinbFwuIC2NB9CTkoAFC4DGjU1X57NtG3D+PDB0qOlK/KMRVwLOsoBPP2VYly4FWrdmWMeNA2JiTFdXvuHDgcmTefDACTTiSsBkZXGt1e3mNsOkJOCjj4DbbzddWeX27OGZXNNH9apCwZUaOXyYO5hcLp6mSUriDPHdd5uuzH9/+AMwY0ZoZ7BrSo/KUmUnTwLLlzOsR44AEyYwsHFxoVlrDaQDB3jULzMTqO2gYUzBFb+cP8/JJZeLe4XHjGFYe/Z0Xliv9b3vAQ8+CEydarqSqlFwpUKXL/OInMvFWddhwxjWfv2cNTpVJDeXS0AHDwL16pmupmoUXCmlqAj48EOGde1ahjQpicskTvvmvpHHHgPatweeesp0JVWn4AqKi4EtWxjWlSuBLl0Y1lGjgIYNTVcXHMeO8bB9ZibQoIHpaqpOwY1QJQ2/XS5ONLVrx7COHQs0a2a6uuB7+mn2TZ41y3Ql1RMGbypSFdc2/G7WjGHdsYMH0iPFqVNcwtq3z3Ql1afgRoBDh3yH0L1ehnX1aqBtW9OVmbFgAfD97wO33GK6kurTo3KYKtvwOyGBgbVzV4dQOHeOTc/T0+27/dIfGnHDSHkNv19+ufKG35Fm8WJuGHFyaAGNuI6Xnw+8+66v4feIERxZ/Wn4HWny89ng/OOPgVatTFdTMxpxHaik4bfbzcPfAwcCP/xh1Rt+R5pXX+UpIKeHFtCIaw95eZzqLCsmhmfiULrh96pVbLFSk4bfkebKFY62GzYAbdqYrqbmNOKalpfH76iCguv+kRUdjXT3QbyxoTVSU9kyNCkpMA2/I82bb/L1IRxCCyi45p06VW5oASCqoAALnzuFB6a0xq5dvGRZ/OfxALNnc6LOsvgDL1wouDb3j38A6GK6CufxeDirXnIRNQD89Kc81D92rNnaAsFBR4dF/Dd7dunQAvzzb35jrqZAUnANOnsW+OMfTVcRnrKySocW4J8PHjRTT6ApuAYUFwN//zs7RtihLWk4at/++k0nUVGcBwwHCm6I7dzJrhFpaVyaeDw5BoiOLv9fjo52/hYfQ5KTOcJeG17L4t8PB1rHDZETJ4CZM4GtW4GXXirTdNuPdVypOo+H77SZmXzKcbmc0371RhTcICsq4v7Y+fOBadOAJ58Mv04STjBsGGeVv/Ut05UEhh6Vg2jDBnaT2LmT+2N/8QuF1pTJkzmvEC404gZBXh7wzDOc2Vy0COjd23RFcuUKzx/v3cvOF06nETeACgqA3/2O7T779QN27VJo7aJePd61+/bbpisJDAU3ACyLTdY6dWJn//R04Mc/Do8WpuFk0qT/70QLA/rWqqGsLE56nDvHY3YPPGC6IqlIfDyXh/bs4Q33TqYRt5ouXOBk0+DBPLGzZYtC6wSTJ4fHqKvgVpFlAW+9xcdir5fXcXz3u866MCqSJSUBqamcrHIyPSpXQXo612Lr1wfef5/nY8VZmjXjhOHKlZysciqNE344fZqTTePGAT/7GdvGKLTOFQ5rugpuJbxe4M9/Bjp35iH2klvq1DHR2QYO5DbIw4dNV1J9Cm4Ftm5lW9MPPuCt6snJOskTLmrVAiZOZDsbp9LOqTK++AKYPh3497/Z8X7QINMVSTDk5HDfclaWM5+gNOL+X2EhDwJ068Y1vvR0hTactW3LNjYffWS6kupRcMHJpvh4XgK1axcnoOrUMV2VBJuTJ6ki+lE5N5eXGh85wsMAPXuarkhCKT+f14tmZgKNG5uupmoicsTNzweee459docP55E7hTby3HwzG8ovXWq6kqqLqOBaFi9x7tQJOHOGR7ymTNEdO5HMqVsgI2bn1P79wBNPsCPFihVs1CbSvTt/iGdmAh07mq7Gf2E/4p47x/fY4cN5MdbGjQqt+ERFOXPUDdvgFhfzf0anTtw4kZEBJCY6c81OgmviRDaSu3rVdCX+C8tH5U8+4WGAmBhg/XrOHIpUpGVLrt2vXs3JKicIqxH35Ek+Dn/nO8CvfsXrKBVa8YfT1nTDIrhXrwILF/Kn5te+xtni4cNNVyVOMmwYn9ROnDBdiX8cH9yNG9kCdft2YMcONh1XC1Spqjp1OAeyZInpSvzjmJ1TJXedZmXxXpjHHwfWrQMOHOBo26+f6QrF6fbvByZM4ESm3ScxAzbiejzc71u/Pn/3eAL1lX13nWZksAXq3r3A1KlAo0Y8xaPQSiDcey/QsCH3q9tdQEbcspcIl/yemnrjS4SLioBLl3y/Ll68/q+ffRY4erRM4VFcj01Pr2n1Ij5//Su7QL7yiulKKheQ4MbHczQs+5WaNuXRuPLCePky//3atYEGDXy/GjYs/ecGDYA//YndKMqKjubXEQmUc+fYlig7m3uZ7Sogwa1fn4+wZdWpw6Zq5YWyfn3/OyOW94NBI64Ey8SJwNChwKOPmq6kYgF5xy3vEmGAwTx6lP2GO3UC2rQBYmMZ3Kq0My1712nJo3i43HUq9uKELZABCW55wQKA3/8eeO893gK+aBGP01XH2LF8X46L4+NxXBzfq8eMCUT1IqX17Qt8/jl/2VXAloNKLhE+eJBBTU72Bes//2FbmPffB370I7Y6bdo0EJ8qEhy//S3nVX79a9OVlC+k67jHjvE29rff5vvDU08Bt90Wqk8X8V9eHpcZs7PteUtFSEtq1Yojb0YGW4V068b12OzsUFYhcmOtW3Of+/r1pispn5GfJc2aAb/8JR+r77uPF2clJmqGWOzFzpNUttjyWFTEx+fnnwfuuguYMYP3u9h925mEt4ICtnHdt89+czK2CG6J4mJexjRnDjdmzJzJUz4KsJjy+OPA/fcDjz1mupLSbBXcEpYFbNgAzJ0LfPkl76FNSNAN7xJ6u3dzJeSTT0xXUpoN58s4wg4YAKxdy8PNaWls5LV4sbY4Smh16cJbLjIyTFdSmi2De62uXYFly7iRY/du7tKaM4d7SkWCza7N5Gz5qFyZI0eAF18EUlJ4E/yTT7JnkEiwnDrFq1ZzcoC6dU1XQ7Yfccu6/XYGd+9edrp44AHuxMrNNV2ZhKuYGKBHD/YwswvHBbdE8+bcjpaZCdx9N9C/P5vE2e1dRMKD3R6XHfeoXJHCQvYLmjePQZ4xA+jVy3RVEi6uXmUjwu3buQPQtLAJbgmvl7PQc+bwIPTMmTxbqbVgqalZs4BbbuHypGlhF9wSlgV8+CHXgk+f5gg8bpzWgqX6srKAUaPYoND0QODYd9wbiYriHuh164C//IWz0Pfcw78ur1uHyI20bw+0aAFs22a6kjAO7rW6d+cNfStX8h2lQweeUjp/3nRl4jR2maQK20flyuTlAS+8wLtyJ03i9ZuxsaarEie4eJEjb1YWe6iZEhEjblmtWwMLFvAY4U038eqSadOAQ4dMVyZ217AhJzuXLzdbR0QGt0SLFmy3k5kJ3Hkn0KcPd2Pt22e6MrEzO1wQFtHBLdGoEfDMM3z86d2bvbJGj+ZdRCJl9eoFHD/OXmqmROQ77o14vewqOWcO0KQJ14IHDza/BCD2MXcucOECO5maoOBWwrKANWt8/5NmzGCr2Fq1TFcmph07BvTsyYMHJr4f9KhciagoTkRs3MhrUP75T14M9be/AVeumK5OTGrVip0x1q418/kKrp969uQ6cGoqsGkT14JffJHLAxKZTE5S6VG5mnJzuRaclgb84AdcToqJMV2VhFJhIa/V+fTT0P+/14hbTW3a8PF5926eHImP56H+w4dNVyahUrcuMH48O5SGmoJbQy1bcmZx/37eytCrFx+hMjODe9m32MOkScBrr11/xWyw6VE5wAoKgNdfB2bPZofK6lz2Lc7StSsvxO7SJXSfqeAGSVwc8NlnutM3EixezKN+ixaF7jMV3CCp6LLvqCi+G48bp4MN4eLMGV6lk5PDa2BDQe+4QVLeZd9RUWx/cvAgm9wNHsw14a++MlOjBEbTprxT9513QveZCm6QlHfZt2WxJ9aCBbw0+dlngV27gE6dgBEjuMFDZ4SdKdRrunpUDqLKLvu+VlERr3NMSeHl3z178vbC4cOBBg1CX7dUndfLazk3bwbuuCP4n6fg2syVK9xGl5LCnln9+/PepKFDQ/f+JNWTnAzUqcMrZINNwbWx/HzgX/9iiLdsAYYM4Ug8aBC/QcRecnM5b5GVFfxb7BVch7hwAXj3XcDt5nvxiBEMcd++Oq1kJwMGAM89B/TrF9zPUXAd6MwZ7pFOSeFa8ZgxfJx+8MHg/6SXyi1ZwledN98M7ucouA538iQnwdxuPqo98ghH4q5ddfDfhPx8TlJlZgKNGwfvcxTcMHLsGJuYpaSwtcqECRyJ4+IU4lCaOpU/OKdMCd5nKLhhKi8PWLqUI/GlSwxwQgKbwktwffwxT4pt3x68z1BwI0B2Nkdht5vvwImJDHHbtqYrC0+Wxe4Yy5cH7welghth9u1jiFNSeIFVQgIfqUOxaSCSvPACX1fmzQvO11dwI5Rl8ZRSSgofqVu14kg8bhxw662mq3O+48f5npuTE5w1dwVXYFnAzp18lE5N5f3CCQmcoW7e3HR1zjV6NPcwjxoV+K+t4EopxcXcpZWSwrXi+HiGePRoPlqL/1au5MGDtLTAf20FVyp09Spb07rd3Hr5jW/wcXrkSB1+8EdREScAd+1ii6NAUnDFL4WFwAcfcCReu5ZbLRMSgIceYtMAKd/Pf845g6efDuzXVXClyi5fBlav5ki8eTMPPSQmcoN93bqmq7OX/fvZCfKzzwK7CUbBlRq5eBFYtYoj8Y4dPEOcmMhN9rVrm67OHjp04NzBkSPsjJKcXPOGgQquBMy5c2zf4nazSfjo0Xyc/uY3I/fwg8fD2fkSger2qeBKUJw+7Tv8kJ3Nb9LERE5wRdK+6fh4ICMj8N0+FVwJui+/9B1+OHrUd/ihc+fwD3FF3T6jozlXUF0KroTU4cPAsmUcic+fZ4gTE3kLYri5dIlbSc+cKf33NeKKo+Xk+PZNe72+ww/t2pmurObS04Fvf5sTU2lp199o4fGU3zjQXxE6ZSB20LYtMHMmv8mXLuWGhZEjucd3/nzg0CHTFVadZQELF3Ji7uWXgRUrOBEVF8fH47i4mocW0IgrNmNZwN69vsMPsbEciceP56VqdnbyJC8B83qBN94I7k0VCq7YlmVxu2DJ4Yc2bXyHH1q0MF1daevX80DBE0/wEL26PIqAGxi2beNIvGIFD6onJPCRs0kTc3UVFXFDhccDvPUWr5YJBQVXHMfrBTZt4ki8ahXfiUsOPzRqFLo6cnN9E1CLFoX2sxVccbSiIt74kJICrFnDXVoJCdx6GczDD243DxDMnQs8+mjwPqciCq6EjYIChtft5nHEAQM4Eg8ZAtSrF5jPuHQJmDaNhwZcLt6+aIKWgyRsREcDDz/MQGVn869ffx246y5OHK1ZwxG6uvbs4WN5bCybDZgKLaARVyLA+fPsRuF2M3wjR3Ik7t3bv+tbStZmX3oJePVVHl80TcGViPLVV5yVdrt5/emYMQxxjx6+fdMeDzB7Ni/vatuWI3lsbPDXZqtCwZWIdfw414dTUthAfvx4rg9Pn+7bmlhi2TJ2wLQLBVcEPLW0bBkwa9b1p3YCcSgg0BRckWsE6xheoGlWWeQa7dtff0Y4KoqbLOxEwRW5RnKy7/gd4HvXTU42W1dZCq7INcaODc4xvEDTO66IA2nEFXEgBVfEgRRcEQdScEUcSMEVcSAFV8SBFFwRB1JwRRxIwRVxIAVXxIEUXBEHUnBFHOh/S7+OTZd2rNEAAAAASUVORK5CYII=\n",
"text/plain": [
"