{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Travelling Salesman Problem\n",
"\n",
" Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Overview\n",
"\n",
"One of the most famous NP-hard problems in combinatorial optimization, the travelling salesman problem (TSP) considers the following question: \"Given a list of cities and the distances between each pair of cities, what is the shortest possible route that visits each city exactly once and returns to the origin city?\" \n",
"\n",
"This question can also be formulated in the language of graph theory. Given a weighted undirected complete graph $G = (V,E)$, where each vertex $i \\in V$ corresponds to city $i$ and the weight $w_{i,j}$ of each edge $(i,j,w_{i,j}) \\in E$ represents the distance between cities $i$ and $j$, the TSP is to find the shortest Hamiltonian cycle in $G$, where a Hamiltonian cycle is a closed loop on a graph in which every vertex is visited exactly once. Note that because $G$ is an undirected graph, weights are symmetric, i.e., $w_{i,j} = w_{j,i}$. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Use QNN to solve TSP\n",
"\n",
"To use QNN to solve travelling salesman problem, we need to first encode the classical problem to quantum. \n",
"The encoding consists of two parts:\n",
"\n",
"1. The route how the salesman visits each city is encoded in quantum states -- ${\\rm qubit}_{i,t} = |1\\rangle$ corresponds to salesman visiting city $i$ at time $t$. \n",
" 1. As an example, if there are two cities $\\{A,B\\}$, visiting $A$ then $B$ will be in state $|1001\\rangle$, as the salesman visits the city $A$ at time $1$ and the city $B$ at time $2$.\n",
" 2. Similarly, $|0110\\rangle$ means visiting $B$ then $A$.\n",
" 3. Note: $|0101\\rangle$ means visiting $A$ and $B$ both at time $2$, so it is infeasible. To avoid such states, a penalty function will be used (see the next section for details.)\n",
"\n",
"2. The total distance is encoded in a loss function: \n",
"\n",
"$$\n",
"L(\\psi(\\vec{\\theta})) = \\langle\\psi(\\vec{\\theta})|H_C|\\psi(\\vec{\\theta})\\rangle,\n",
"\\tag{1}\n",
"$$\n",
"\n",
"where $|\\psi(\\vec{\\theta})\\rangle$ is the output state from a parameterized quantum circuit. \n",
"\n",
"The details about how to encode the classical problem to quantum is given in detail in the next section. \n",
"After optimizing the loss function, we will obtain the optimal quantum state. Then a decoding process will be performed to get the final route."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Encoding the TSP\n",
"\n",
"To transform the TSP into a problem applicable for parameterized quantum circuits, we need to encode the TSP into a Hamiltonian. \n",
"\n",
"We realize the encoding by first constructing an integer programming problem. Suppose there are $n=|V|$ vertices in graph $G$. Then for each vertex $i \\in V$, we define $n$ binary variables $x_{i,t}$, where $t \\in [0,n-1]$, such that\n",
"\n",
"$$\n",
"x_{i, t}=\n",
"\\begin{cases}\n",
"1, & \\text {if in the resulting Hamiltonian cycle, vertex } i \\text { is visited at time } t\\\\\n",
"0, & \\text{otherwise}\n",
"\\end{cases}.\n",
"\\tag{2}\n",
"$$\n",
"\n",
"As there are $n$ vertices, we have $n^2$ variables in total, whose value we denote by a bit string $x=x_{1,1}x_{1,2}\\dots x_{n,n}$. Assume for now that the bit string $x$ represents a Hamiltonian cycle. Then for each edge $(i,j,w_{i,j}) \\in E$, we will have $x_{i,t} = x_{j,t+1}=1$, i.e., $x_{i,t}\\cdot x_{j,t+1}=1$, if and only if the Hamiltonian cycle visits vertex $i$ at time $t$ and vertex $j$ at time $t+1$; otherwise, $x_{i,t}\\cdot x_{j,t+1}$ will be $0$. Therefore the length of a Hamiltonian cycle is\n",
"\n",
"$$\n",
"D(x) = \\sum_{i,j} w_{i,j} \\sum_{t} x_{i,t} x_{j,t+1}.\n",
"\\tag{3}\n",
"$$\n",
"\n",
"For $x$ to represent a valid Hamiltonian cycle, the following constraint needs to be met:\n",
"\n",
"$$\n",
"\\sum_t x_{i,t} = 1 \\quad \\forall i \\in [0,n-1] \\quad \\text{ and } \\quad \\sum_i x_{i,t} = 1 \\quad \\forall t \\in [0,n-1],\n",
"\\tag{4}\n",
"$$\n",
"\n",
"where the first equation guarantees that each vertex is only visited once and the second guarantees that only one vertex is visited at each time $t$. Then the cost function under the constraint can be formulated below, with $A$ being the penalty parameter set to ensure that the constraint is satisfied:\n",
"\n",
"$$\n",
"C(x) = D(x)+ A\\left( \\sum_{i} \\left(1-\\sum_t x_{i,t}\\right)^2 + \\sum_{t} \\left(1-\\sum_i x_{i,t}\\right)^2 \\right).\n",
"\\tag{5}\n",
"$$\n",
"\n",
"Note that as we would like to minimize the length $D(x)$ while ensuring $x$ represents a valid Hamiltonian cycle, we had better set $A$ large, at least larger than the largest weight of edges.\n",
"\n",
"We now need to transform the cost function $C(x)$ into a Hamiltonian to realize the encoding of the TSP. Each variable $x_{i,j}$ has two possible values, $0$ and $1$, corresponding to quantum states $|0\\rangle$ and $|1\\rangle$. **Note that every variable corresponds to a qubit and so $n^2$ qubits are needed for solving the TSP.** Similar as in the Max-Cut problem, we consider the Pauli $Z$ operator as it has two eigenstates, $|0\\rangle$ and $|1\\rangle$. Their corresponding eigenvalues are 1 and -1, respectively.\n",
"\n",
"Now we would like to consider the mapping\n",
"\n",
"$$\n",
"x_{i,t} \\mapsto \\frac{I-Z_{i,t}}{2}, \\tag{6}\n",
"$$\n",
"\n",
"where $Z_{i,t} = I \\otimes I \\otimes \\ldots \\otimes Z \\otimes \\ldots \\otimes I$ with $Z$ operates on the qubit at position $(i,t)$. Under this mapping, if a qubit $(i,t)$ is in state $|1\\rangle$, then $x_{i,t}|1\\rangle = \\frac{I-Z_{i,t}}{2} |1\\rangle = 1 |1\\rangle$, which means vertex $i$ is visited at time $t$. Also, for a qubit $(i,t)$ in state $|0\\rangle$, $x_{i,t} |0\\rangle= \\frac{I-Z_{i,t}}{2} |0\\rangle = 0|0\\rangle$.\n",
"\n",
"Thus using the above mapping, we can transform the cost function $C(x)$ into a Hamiltonian $H_C$ for the system of $n^2$ qubits and realize the quantization of the TSP. Then the ground state of $H_C$ is the optimal solution to the TSP. In the following section, we will show how to use a parametrized quantum circuit to find the ground state, i.e., the eigenvector with the smallest eigenvalue.\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Paddle Quantum Implementation\n",
"\n",
"To investigate the TSP using Paddle Quantum, there are some required packages to import, which are shown below. The ``networkx`` package is the tool to handle graphs."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"ExecuteTime": {
"end_time": "2021-05-17T08:24:17.197426Z",
"start_time": "2021-05-17T08:24:12.896488Z"
}
},
"outputs": [],
"source": [
"# Import related modules from Paddle Quantum and PaddlePaddle\n",
"import paddle\n",
"import paddle_quantum\n",
"from paddle_quantum.ansatz import Circuit\n",
"from paddle_quantum import Hamiltonian\n",
"\n",
"# Functions for Salesman Problem\n",
"from paddle_quantum.QAOA.tsp import tsp_hamiltonian # Get the Hamiltonian for salesman problem\n",
"from paddle_quantum.QAOA.tsp import solve_tsp_brute_force # Solve the salesman problem by brute force\n",
"\n",
"# Create Graph\n",
"import networkx as nx\n",
"\n",
"# Import additional packages needed\n",
"from numpy import pi as PI\n",
"import matplotlib.pyplot as plt\n",
"import random\n",
"import time"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Generate a weighted complete graph"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we generate a weighted complete graph $G$ with four vertices. For the convenience of computation, the vertices here are labeled starting from $0$."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"ExecuteTime": {
"end_time": "2021-05-17T08:24:24.302458Z",
"start_time": "2021-05-17T08:24:24.060967Z"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADnCAYAAAC9roUQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAyu0lEQVR4nO2deXhV5bX/PyskYQpDRJBB5CgiiiOioshgxWoltlarVasF8Vbbam2rVTGt2trbGmwd6yz+CnJrHW+1V0NV1AqCM1IUsUCBoBhRURlzMnHe3x9rh5xzEiAJ+5y998n6PE8ezcnO3ivh5Lvfvd61vkuccxiGYRjZIS/oAAzDMNoTJrqGYRhZxETXMAwji5joGoZhZBETXcMwjCxiomsYhpFFTHQNwzCyiImuYRhGFjHRNQzDyCImuoZhGFnERNcwDCOLmOgahmFkERNdwzCMLGKiaxiGkUVMdA3DMLKIia5hGEYWMdE1DMPIIia6hmEYWcRE1zAMI4uY6BqGYWSR/KADMHKTWGl5HjAY2A/oDBQCtUAcWAasqCgrSQQXoWEEg9g0YMMPPJEdD5QAY4ADgARQD4j34byPfPQp6wPgFaAceNFE2GgPmOgau0SstLwYmAxcAXQDuqIC21IcsAXYCNwMTK8oK/nK7zgNIyyY6BptIlZa3gW4EfgBuqLt4sNpq9AV8APAlIqykiofzmkYocJE12g1sdLyMcAjQDGar/WbOPAVcFZFWcm8DJzfMALDRNdoMbHS8o7ArcD5ZEZs04kDM4DLKspKarJwPcPIOCa6RouIlZYXAc8Dh5EdwW0gDiwETqooK9mcxesaRkYw0TV2iie484ChQKcAQqgGlgKjTXiNqGPNEcYO8VIKzxOc4OJddyjwnBePYUQWE11jZ9yKphSCEtwGOgHDgVsCjsMwdglLLxjbxatSeI7s5nB3Rhw40aoajKhioms0i1eHuxzoH3QszVAJDLE6XiOKWHrB2B5/QOtww0gxMDXoIAyjLdhK12iC19pbSfB53B1RDfS3lmEjathK12iOyWhrb5hJoE0ahhEpbKVrpOC5ha0B+gUdSwuoBAaaO5kRJWyla6QzHnULiwLdgeODDsIwWoOJrpFOCWrPGAW6ABOCDsIwWoOJrpHOGFrnhxskecDYoIMwjNZgOV1jG14+dws+VS1ccGyMM0cMZL89utEhT7jthWXc9uJyP06dTBzoWlFWYm9kIxLYStdIZjCw1a+THTSgBxvidXyyIe7XKZsjgcZtGJHARNdIZj90ppkvXP7YIs6e9jpLKjf6dcrmqEfjNoxIYKJrJNOZ6ORzGxDC5Q1hGDvERNdIppBoiq7ZPRqRwUTXSKYWnc4bJRxgo3yMyGCiayQTJ5qim9GdOsPwk/ygAzBCxTJ8fE+cdcRAjowVc+CAHgCcOGwP9izuzPNLPuX5JZ/6dZl8NG7DiAQmukYyK/Dx6efIWDFnjBi47fNh/XswrH8P1nwV91N089C4DSMSWHOEkUKstHwBcHjQcbSCBRVlJUcEHYRhtBTL6RrpvEJ08roJYG7QQRhGazDRNdIpR1uBo0AVMCvoIAyjNZjoGum8CGwKOoiW4JzbCLwUdByG0RpMdI0UPEPwm9BVZGhJ1FazYf7DdaunnnKaiHQIOh7DaCkmukZzTCfs7428PDa99fdBwBPAv0XkRyJi7cBG6An3H5YRCN6wxwcIb9NBXODeRM2WS4FVwL7APcBqEblWRHoFG55hbB8TXWN7TAHCOmn3S8kv/IVz7k7UYews4G2gN/Bb4EMRuUNE9g4ySMNoDqvTNbZLrLR8NPA84XLxigNfrygrmZ/8oogIMA64CjjZezkBPA780Tm3IKtRGsZ2sJWusV0qykrmATMIT5ohDkxPF1wAp7zsnJsAHAw8iIruWcDbIvKSiHzDE2fDCAwTXWNnXAYsBKoDjqPai+PynR3onFvsnDsf2ButxNgEfA34B7BIRCaKSGEGYzWM7WLpBWOnxErLi4B5zrmhIuLL/LRWUg0sBUZXlJVsbu03i0gP4IfAz4D+3ssfA7cC07x6X8PICia6RosY+POHixJVG1d06N67T15BVj3D4+gK96S2CG4y3ur2e8CVwDDv5Y3AvcDtzrnKXTm/YbQEE12jRYjIj+iQf89u4y+qLxr+jXqRvGyseONozfDlFWUlvhmVi0geutl2Jbr5BlAH/AW4yTm3xK9rGUY6JrrGThGRY4A5QAFw3qCrn1kNPAoUk5nKhjharnaWt5mXMUTkKFR8T6dxj+MZ4I/AK87+QAyfsY00Y4eIyB5o11cBcIdz7iFPCIegDRTV+NcyXOWd7wFgSKYFF8A596Zz7ky03vduVPBPQW8yr4vIGbvYZpwPnAgcgf29GdhK19gBIpIPvIA+gs8HjnfO1SYfEystLwbOB64AugNdaIW4uEQCycvbjOZWbwJmeB1xgSAivYFLgJ8ADZ1tK4CbgRnOudaUz/0E+D06PFOAzcDv0JRJqL0tjMxhomtsFxG5GS3RWgsc7pz7ZHvHxkrL84DxaK50LLpRlQDqaRQd533kO+c61K79T6eaNUsSXYaOOjW/e+9ZntlOKBCRLsBk9Offx3t5HXAHcLdzbt1OTnEI8DpN0y9b0N/L7d7Hzs5j5BgmukaziMhZwCOoaH7NOdeqR31PhPdBH9s7o2PSa9DH92XAitVTT3kdOAo4xTlX7mP4vuGlFk5H875Hei/HgT8DtzjnVm7nW/+BphW2t+pvWDH/FfgNsMaPeI3wY6JrNEFEDgLeQFMFP3XO3ZGh6/wO+BVarvXzTFzDL5LajK8EJngvJ9B89x+dc28nHb4v8B7QkgqPOqAWOBMVaiPHscS+kYLXSPA3VHAfAu7M4OVme//9egav4QtJbcYlNLYZbwW+C7wlIv8UkZM9cS4Fmt18++qrJunqAqArWi1htANspWtsw6tffRL4FvAucIxzLmMbPl6zwpeo6OzpnPs4U9fKBCIyAO1y+yG6iUj//v0/WL169b75+fkF6cffc889vPXWWyxevJgJEybwq1/9ioKCbYfFUTG3ycY5jq10jWRKUcFdD5yeScEF8CohXvY+PSGT18oEzrmPnXNXAXuh7maVF1988QF1dXVNBPe2227j6aef5tprr2X69OmsXLmS+fNTfHvy0OoGI8fJDzoAIxyIyEnAf6PVBec657K14poNlKCbTg9m6Zq+4pzbAPxx5syZ959zzjlrCwoKUnK5y5cv57nnnmPs2LHEYjFEhJkzZ7J5c4rGvgd8ms24jWCwla6BZ/b9MFrWdb1zLpsTdhvyuid46Y3IMnHixO8XFBRsTX/90UcfZeHChdTW1jJixAjuvfdeAIqKihoO2Qxcn71IjSCxlW47x5sr9je0pbccXe1mkw9Qx68BaE5zUZav7xf5wLVofjqFZ599lhtvvJFJkyYxcuRInnnmGQCcc3j2vl9io+TbDZFeWRi7hrfTfg9wGLAS+L5zLqsNCp63QWSqGHbA6WjFRwrr1q2jW7duTJo0CYA+ffrwySefUFFRgYhQV1dXXVdXdz1afma0A0x02zc/AiahO+enOeeCar9tEN0TA7q+H/wYKEp/cffdd+fAAw/koosuory8nIcffpji4mJisRgAmzZt6tS9e/dfisjFXheckeOY6LZTPOew271PL3TOvRtgOC94/x0TkEm6HwzZ3heuu+46BgwYwB133MFhhx3GNddcA0BdXV3tXXfd9WV1dfVg4C50mvGvRWT3LMVsBIDV6bZDRKQvsACdovAn59zPAg4JEVmIpjm+7px7YSeHh5E/oDW7rRkDVLV06dKB+++///FoyVlym/F0tM3Y6nZzDFvptjNEpAD1wu0PzEPdwcJA1PO696HeEnU7OihpkVMDPDB06NAvnXNPACOB49DNzM7AxcAyEXlMRI5s7lxGNDHRbX/ciLqArQW+65zboUhkkajndVcA+wP/i3oCNzvpImkYsUOtLPUTZY5z7hTgIHQK81bUk+FNEXlZREqiXlZnWHqhXSEiZ6P1uPXAcc65JqPMg8IrXfsKdSPbwzn3WcAh7QoD0CeIi7zP0zfIalCjnPN2dBKvzfin6IZnd+/l91Gx/mu6t3HY8JzmBtPoNFeImvtsc5oLk51ntjDRbSekOYdd6pzLpJFNmxCR59H0wveccw8HHY8PNEwhvho1wOmKCs4XwOFofe5OEZHuqID/HBV00Nrm24H7vY64wEnyVC4BxgAHsANPZfRJ+wPgFTSt8mJ7EGET3XaAiPQE3kR32P8CTAzj7C8RuRLdkJrunLsg6Hh8pBA4Cc3bPofm0lv9+/cMgs5BV9EHeS9vQvPJtzvnAvHk9aaHTPbi6obeXGSH35SKQ83dN6ITOqYHOT0k05jo5jheDvAp4Jtot9eoTBvZtBUROQwdt/4xMDCMN4Yw4DW1fAP19v2a93Idaoh+k3NucTbiiJWWd0H3CH6Armj9qDOuQlfADwBTKspKQvle3RVMdHMcEbkGbe1dDxwR5hIk7waxFugNDHPOfRBwSKFHRI5AxfcMGjfGZ6H+vHMydeOKlZaPQSeLRH4idLYx0c1hROQbNPb0lzjnQj+ZQEQeAr4H/Mw596eg44kKIrIPOs/tAhpF8G1UfP/mnKv34zqx0vKOwK3oMNJMiG06cbSS47KKspJmK0KihpWf5Ciec9hf0dzab6IguB5Rr9cNBOfcSufcT1Bv31+jAy+PQGuyl4rIJbvaZhwrLS8C/kn2BBfvOucDL3nXjzy20s1BvD+u+WiH1zPAqdk2smkrXpnUGnRjZbewl0WFFe89MAn4BVq2BVo1cSdwl3Pu89aczxO8ecBQWjb7zW+qgaXA6IqykkibvZvo5hjeJssMYCJasH+Ec259kDG1FhFZgpYbjXPOzQ06nijjTTP+NtpmfJT3cjWNbcb/2dk5vJTCP4HhBCO4DVQD7wDHRznVYOmF3OPHqODG0ZE764MNp0087/3XUgy7iHNuq3Puf4Gj0WnGz6DC+WM07fC4iBy1o3OgOdzDCFZw8a4/HLgl4Dh2CVvpRoex6Bt/APAYaliTguccNgedMHuuc+6v2QzQL0SkBBWHN5xzRwcdT64hIsPQmtrz0PcKwFx0021WcirKq1J4juzlcFtCHDgxqlUNJrrR4Cp0cyQPbZOtBu71Xq+HJs5htzvnfh5IpD4gIkVot1YHYPcAfX5zGhHpj7YZ/5jGNuMleG3Gg65+pgOwHH1PhY1KYEgU63gtvRB+vo0Kbhf08UrQVcdFaPtk7zTnsFfQus3I4pzbDLyGvj+PDzicnMU5V+mcuxoYiK581wDDgD8DK+Or353tnCsOMsYdUAxMDTqItmCiG26Gom27zZX6dEX795ecdtppD6Lph08Il3PYrmB53SzhnNvonLsZrXKYBCzO69i1f8f+Q0d5RkRhpDNwodeCHClMdMNLN3aeSysEdv/LX/5yzkUXXbQVOMM5tzYr0WUeq9fNMs65WufcTOCQ3U/75TREwl5mmEBreCOF5XTDiQBPo45NLdoxrqurqy0oKHgMuBDN+UYar9Tpc/Qxct8wty/nGp5b2BqgX9CxtIBKYGCU3MlspRtOrkSnCGwT3Pr6eqqqqkgkmn9vFRQUFALfQTfT9spCjBnFObcVeMn71Fa72WU8+qQVBboTsby/iW742Af4DZqz3call17K97//faZMmcL69eu3972dUcPo92ictxVlLK8bDCWkvf9CTBdgQtBBtIb8oAMwmtDkDXT11Vezfv16/vCHPzB58mTmzp3Lt771LRKJBHl5Te6b+ejd/wVgb1polB1SGvK6x4tIB2/1a2SeMbTOD3e7lJ12MEcMKqZ/z87Ubk3wr4/Wc8OsD1j+mW+dvHnoJnJksJVu+KgkabjhW2+9xbPPPst9993H4MGD2Weffbj//vs544wzuPnmm9mwYbtDAxzaShtZnHOr0Fbmnqh5i5FhvHzuML/Od85Re7G5pp7/W1TJ5up6vja0DzMvOIqO+b5Kz7BYabkvN4lsYCvd8PEuSf8uRx55JP/4xz/o3r078+bN47nnnmPx4sUsXLiQ6dOn8/TTT3Peec2O2spD6y+jzmy0lOlEdNyQkVkGowMxfeH0e17lnQ+1t2XPnp2ZN+V4+vXozL59ini/cqNfl0mgce/URyIM2Eo3fPwHLYPZ1mnTr59uIo8ePZoFCxbQq1cvTjjhBE466STeeecd6uubtUrNRxsMoo7ldbPLfnhdjn7QILgABd7qdmvC8dkmX/1q6tG4I4GJbjh5HN1B/gIv1dBQ2te3b99tBz322GMMGjSI/PwmDyxb0M241ZkPNeP8E13JHCMiUdlRjzKd8Smfm0yXwg7cdMahAEx7ZSWf+yu6DV2akcBEN7y8DhwI/Kuurq5WHRshLy+PrVu3ctFFF9GjRw9+9rOfpX9fNfAyal4SeTyXtDfRlfu4YKNpFxTis+gWdyngrz84mhGDivnrmx8y9dl/+3l60Hg7+n3STGGiG24+LS4u/sV9992Xv2XLlm0vbtq0iWOOOYYHH3ww/fgE8Ck6MTaXul4aqhhODDSK9kEtPr53BvTszBM/GsVhA3ty98v/4ZdPvufXqZNxQGT8dU10Q4yI9F2/fv2jl156ad4NN9wwC83zup49ezJ58uTm0gpxdNT3pmzHmmEsr5s94vgouv/7o1EM7l3Emq+q6FTQgetOGcZ1pwzj0D17+HUJ0Hjjfp4wk5johhTPOewxtBVz7g033PBt1Pm/kubv6lXA99GRJrnGG+iNZH8R2RNNNRyD5q1fA1aiNpeFQQWYQyxzzvlW1dS3hzZV7lnchQuO3Xvbx5A+vqbn84Flfp4wk5j3QkgRkVuBn6MiOyLJyKYHutE2Cu0aarjLT0VHreciMnz48NnHHHPM+ClTpiwZNGjQ3uiOdScaTbi3oF4NB5AD3hPZRES6o++nsYiMG3j5E6PyCiKTIgV9/3etKCuJhJiZ6IYQETkHneRbBxznnHs1/RDgDOC7aInZ68Dfsxpk5ikGTgC+CXyjrq6uZ21tbUHXrjvsTt0MXIrOiDO2g4j0AkajnVzj0BE42556+55/Gx377htQdG1iQUVZSWSaZ6w5ImSIyMHAA96nlzUjuKCr28e9j1yjEB2a+B00jdINkIKCAgoKCnb4jUARmmKZkckAo4Y3VWRs0sfBaYfUoymcucDcgt32PAX4ERkoHcsACTTuyGCiGyJEpCfwJGriMRO4O9CAguFSdFpGR9pWBjQKTTnkgpF7mxCRvdAVbIPIpjcO1KBPR3NQwXrdObetPCZWWl6H3ryKshLwrlEFzAo6iNZgohsSRCQP+B+0nfFfwI9c+8v9dAaup/lJGS2lFt1wnO9LRCFHtIB7XxpTBWOBQWmHbUF/H3NRoX3LObejEqsX0Y3LKIjuRhotQCOBiW54+BVwCvAVOjo9MiUwPrKVVo75ds5RWVnJ8uXLOe6440CF+xvkqOh6N+dhpIps37TD1qOz8hpEdqFzrsWtvRVlJYlYaflN6MbsrtwAM00VcFOUDMzBRDcUiMjJ6ArPAd/z3LXaI7XonLc9m/tic1aWIkJdXR1XXnklM2bM4MADDyxA0xPXZjrYbOBN0DiMxlTBGKBX2mGf4eVjvY/3kseot5HpwO938RyZJo8I5u+tTjdgRGQftFJBgF87554NOKSgKUc3R5owdepULrvsMtatW7fttZqaGmKxGGPHjuWpp55qeHk/tLQucohIoYiMEpGrRWQW+uTzNnALejPphY7SeQj4IVoi19c5d6Zz7g7n3CIfBJeKspKv0A3dsD5xxYFpXpyRwla6ASIiXYC/oX6xTxP+lUU2KEfbmLunf2HNmjU88cQTVFZWMnjwYCZOnMj+++8PwODBg1m5cmXDoXF0h35elmJuM9603ZE0pgqOoal5ywoaUwVzgYos5funAKc3E08Y+BK4Ougg2oKJbkB4GyD3AoeitbYT/Vih5AAvs5287nnnnceaNWuYNm0aU6dO5eKLL6a4uJhOnTrx3nvvMW3atIZDOxDS7jTPKW0UjSJ7FI0NHg0soVFkX3HOfZzVID0qykqqYqXlZznnng/ZKPY4cFZFWUnVTo8MIdYcERAicglwJ7oZMNI5tzjgkMLEIuCQ9Bfr6+s54IADWL58OaAr32eeeYaePXsycuRI9t57b9C8+HtoHjTwN7eI7EbTRoQOSYc4tFqlIR/7inPu8yyHuV1EJL/3d659r1PssP1D0qUWB6ZXlJVcEnQgbcVENwBEZBS6iskHznHOPRJwSGHj12g1R5NuiN/+9recc845DBkyJP1LCbQ0qho1/VmY6SCbQ0T2oGkjQnKTQT2ao20Q2fmefWXo8NJfj9Ehv6Tv96YmCvvtVy95eUE+QVQD7wDHV5SVRMZVLB0T3SzjdQe9gxrZ3OacuyzgkMLIcDQf26RcqaamhsLCQkSEeDxe37lz563AWjQn/gxaY+rb5IOdISIDSW1EGJoeMtrtldyI4NtUxkwhIrujv9OjgS8K+ux9Zv8L7rgV/flaVdbnE9WomdPoirKS0P/+doSJbhbxnMNeRMt+5gInOOfabefUTngDOIhU4a0GtiYSiQ6rVq3q1K1bt0Tv3r1rROQTYDIZbgf18vCDSa2RjaUdtgV4ldRGhEgZ8IjIIOA5VGBXAyc555bGSsuLUJvNw8ju5locfXI5KeqCCya6WUVEbgN+hjqHHe6c+zTYiEJNJ+AXwDXe5xvRzbGuzrl4IpHo3qFDcmqUKrSF+M9+BeA1IhxAqsj2SztsA00bESJ7IxWRQ4Bn0Z/zXeBk51xlw9djpeUd0fK1yWRHeONozfDlUU4pJGOimyXSnMPGOedyYWhkptkPHdWTh7ak7syAZTVNV54txmtEOJTUnGx6I8LnNG1E8G16bpCIyHGoW113tIrk2865Dc0dGystHw08irrBZUJ842iN8lkVZSWhL/1rDSa6WcBbPbyOvjkvcc61RyOb1iJoK+9RpO7274hq1IegRSVWXrrnCBoFdjRN64Mr0RVsQ07237noiSEiZ6ANF4XAE8D3d5YWiZWWd0F9nC9ENzL9aBmuQm+y04Cro1oWtiNMdDOM5xz2NpoLnAmcn4t/tBngYPRG1Zo/5E3ARUCz1SBerelRpDYipJ9/JY2r2DnAqlz/9/LKF+9Ab3R3AT9rzeo9VlpeDJwPXIHetLrQum7XBCq2G4GbgBlR7DRrKSa6GcTLCf4dNbL5FzCqnRrZtIUY8D5povjhhx8ye/ZsFi9eTElJCSeccELylx3wIJpvbGhEOIbURoT0kqcPSEoXOOfW+P+jhBNvY/B3wC+9l34FlLX1JhMrLc8DxgMno7/vYaig1qOCLui/kUPLJfNobASZBbwUNfOatmCim0FE5DrUyOYrdOROezWyaQt56Mp1m+iuXr2aSZMmMXDgQEaNGsW9997LnXfeyZgxY7Z9U1VV1bquXbvOQIX2cJo2IiwitRHhs8z/KOFDRPKB+9Eb1FbgQufcdD+v4YnwPmhuvjPqj1yD5muXASuiMmLHT0x0M4TnHFbufXqyc+65IOOJKC8BX2v4ZMOGDaxdu5ahQ7UU9oorrqBPnz5cccUV29zHqqurGThwYIMpzlZgAY352PnOuZx9bG0pItIV3QQrQQXwTOdc+Y6/y/ALcxnLACIymEbnsOtMcNtMOUlDJnv06MHQoUOpq9OKrKKiIlasWJFi91hfX7/1wgsvfAI4EejpnBvpnLvKOfeMCe62pocXUcH9AjjeBDe7mOj6jNc6+b80OofdEGhAEUVE5JFHHllRW1vb5GsFBQVs3LiR+fPnc8EFF6R8raioKO+GG2743Dk3OwqdX9lERGJop99ItLzuWOfc64EG1Q4xlzEf8TYm7kNrPZejZTc5vzHgB97vrqERYSww9rzzzhuwadOmZo9/6KGHGDBgAEcc0WQIrKCrXCOJnTU9GNnDRNdfLgHOQ8tfTt9eYbmxrRHhEFIbEXZPPmbr1q3rVq1alTds2LDdGl5zzlFbW8ucOXOYOnUqHTp0YN26dfTq1QvVbQAGok8a67Pwo4Se1jQ9GFnAOWcfPnwAx6LdZg44K+h4wvaBOoaNBK5CjWnW01g+1PBRCTyMjv8ehm70/tI5V+OSuP76613v3r3d2Wef7YYNG+YuvfRSt3nz5uRD1jvnTgr6Zw7DB3AGWjHggMeBTkHH1N4/2u1K1ytnGUxjOUshOqMruZylRakBzznscfTJ4Vbn3KMZCTpCiEgnUhsRRtG0EWEVqRMRVjpPKZKYg26mFYIuEjp16sQJJ5zAueeey/Dhw+nfv3/65fPZzpy19oSI/AT4E5pyuRP4ucuRluUo025KxpIKt0tQl68D2Hnh9geomUk58GJzIuy1kr6EtpC2W+cwESlCGxEazGFG0rQR4d+kNiJ81IJTd6TR7KalxFF7yKWt+J6wk4+aJR2CVh08ib43m9BM08MvganN3NCMAMh50fVaFCejLYrdgK7s3DglGYfa9W0EbkZd67eVHonI7cBPaWfOYV57c8NEhLHACFL3CBy6YZPciNDW382bwJEtOdA5VyUidwNXtvFaYaQXMBt9KuvqvRYHfox24G3DWwTcR2PTww+cczOyFqmxU3JWdD0zjhuBH+C/GccDwJTVU0/5NmoSUgeMdTlcfiMivUnd9DqU1JvXVtScPbkR4UufLj8RuIfm/w1rgNqampqit99+WzZs2PCbCRMm/J4sGplnmEHo77MvTVf7ceAstDSxoenhMWAC1vQQWnIypxsrLR+Dmp4U46/LfcMf/X+5rfVnddzzwO41a94HNQjJKcEVkQGk+sgekHZILboCbcjJvuaca76+a9d5CNgD+A16g+uA3vzeQgVnbvfu3SfX1tb+GHDOuVwR3IOBf6KVGM05rXVG3+cniMhydINyJJp+KHHOvZGlOI1WkFMrXc9g+VbU8SjjBsuJuhqqKxYuLeyzz6Fr7p4cWYNlLwe4N6kiu0/aYXFSJyK86bJv3tMfzWl+hObbt+XYReRU4Cl0hT06y3FlgnGoiO40HZZIJDYffvjhny9atGhvkiY9ZCFGow3kjOgGNUrEORcXkUiNEvFEdn8aUwXjgAFph21Eu5cacrILnHNN28NCgoj0QFd4ALs55zYGGc8uUoKmCVLSKc45tmzZQlFRUcrBiUSCdevWMXLkyCUVFRVfd9b0EGpyQnQ9wZ2HDc1rFq8R4WBSc7K90w77gtSJCIuiVl4kIvPR0rRTnXP/F3Q8bWRf1Aa0a/KL69at47vf/S5XXHEF48ePp2PH1HHo9fX1TkQqOnTocDjWFBJqIp/T9VIKzxOc4OJddyjwXKy0PPDx0N4O9nAaUwWj0bxgMp+QWiP7gYt+y/JsVHS/DkRVdI8kKW0CKrjf+c53GDduHBMmTGj2m/Lz8wVNv8xGSyIjNQyzPRF50UVzuIcRnOA20AkVulvQduCs4TUiHElqI0LXtMMqSBXZFTlYt/k88Gui7b3QjbRNs3g8zpAhQ/jtb38LwN///ncGDRpE165dGTJkSPKhHYEDgb8B30QrSoyQEWnR9aoUzie746B3RGdgcqy0/OFMDtPzSoPSGxE6ph22lNRGhA8zFU+IeBPNRe8nIntF9GdeQNrG2bp16/jwww9xznHJJZewatUq+vTpw4YNG7j22msZMWJE8uGd0ffFFWjJpBEyIpvT9epwl6OPVGGjEhji11A9rxHhWBrzsUfQ9IaZ3oiw1o9rRw0ReQo4FW0K+H8Bh9NWfgf8nKSnlYkTJ/LRRx9x1FFHceONN/Lpp58yY8YM6urquOaaa5o7x0agR3bCNVpDlP10/4DW4YaRYnRKapsQkd4icrqI3CYi7wBfouVDVwFHo/9ub6MdcqcCvZxzhzrnLnXOPd5eBddjtvffrwcaxa5xLeoKtu2mfcMNN9CpUyfmzdMHqD322IMePXoQj2+3ai9XapVzjkiudL3W3kqCz+PuiGqgf0ummopIf1JrZIelHVJH00aEKJdEZQwRGYIaFn0B9Inw5mA+8FwikRiVl5fXCeCDDz5g0qRJjBs3jpKSEn7xi19w+eWXc+6556Z/r0PrmA/McsxGC4jqSncyaTu8ISSB5ptTEGVvEZkkIn8Wkf8AH5NqaRhHTXR+AxwP9HDOjXZqc/icCe4O+Q/aINAL3diMKvXjxo27ZNmyZVJTo8UwBxxwALNnz6ZTp07MmjWLKVOmNCe4oCvkZr9gBE/kVrqeW9ga1AE/7FRu+fe8geuemrofqY0I6baDm0htRHg7zI0IYUdEpqGeG6XOuTaneYJERA4F/tGrV69+7733Xl3fvn3zvHrrHZFA30snok9GRgiJ4kp3PFpWE3oSdTV9Nr87+wv0Ue8+dPWxJ5qj/TtwOboptptzboJzbqpz7lUT3F0m0nldb9LDXKDfF1988fLbb799uIhsQNMG26MO+Bz1MDbBDTFRLBkroWkNaiiRDgX5nWPDe1avXLCW1BrZJRHONUaBF1GBGi0iXZxzvlSRZAMRORP4C+oo9jjw/W9+85s16Cj6+UBRM99WA3yIPkm1503USBBF0R1D6/xwt0vH/DxKTz6AUw7pR1HHfBZ/vIHfzfqAf3203o/TI3l5FA0/efFXLz1wSA42IoQW59wXIrIAfYoYiw5kDD07mfTwLtrwUE6qJ0OV97WT0DIxI+REKr3g5XPTd/bbzHWnDOP8UTHWba7h+SVrOXyvYv7ngqMo7lLg1yXIK+g4eNDVz/h2PqPFRCbF4G2u/h64AxXcXwI/bcb74mX056lErTUTwP3AcZjgRoaorXQH41NrY6+uhZw5YiBbE45zH3iDL7bUUp9wnD58TyYdE+O2F5f7cRnQP4zB6K66kT1mA6WEXHQ9n4z70UqXlkx6eJVGR7iGEVNGhIjUShcdV+JL0fd+e3SjMD+PyvVxvtii+1bvrdGp1MP6dffjEg3Uo3Eb2eVV9NH7YBEJZaWL1879FCq4VcC3WjlaxwQ3gkRNdDvjUz539yKdfLKltlHDq2p1Ed27W7qNwS4hhMcbot3gnKtBNy4BTggyluYQkd3RDb8JaCPH8c65WcFGZWSDqIluIT6J7rrNurrtWtiYYenaUcsgP9/kqzOj0NSMxsgOoczrikgMrUQYiTZyHGujddoPURPdWnx6pFr+2SZq6xP079l526r3kD17AvDBWl/3JBxa0mNknwbRPcGblhE4XtPDa2jKaREwykbrtC+iJrpxfBLddZtreeKdNXTIEx76r6O54+zhfOuQ/myuqefB11b7cYkGHBq3kX3eR83a+xECHwIR+RqNk31fBsbZaJ32R9REdxk+Vlxc//T7zHytgt2LCjlx2B4s/Gg9E//8Bl9u8a8hzCUSHes3rVvp2wmNFuPVRjesdgM1NveaHp4FuqNND99wzm0IMiYjGCLlveDV6W4mQhtTiboaPrr5O3G0NXO+9/Gac26n7mM+cyLwU+Az4HUvjpSJurmIiJwH/A/wrHPu5IBiSG56uANtesjp37uxfSIlugCx0vIFwOFBx9FSaj+rqP7kzz9pzoJyCY0iPJ/Mjc8R4BHgFBo7mbagaY884B10BTYPeIskD9dcQET6oimGOFDsVTVk69qCGpL/0nupFLjRuhPbN1FrjgB4BbXsC8XGyE5IFPaJ3QPcgM4tO9b7OALtrBsGXOgd+6mIvEqjCL/jk/HNsWhZUnLraLJ3xWjUJKUafYJYAbyAllvNRwUrsjjn1orIe+g05FGx0vI5aLPKfujPW4hu0MbR9NWKirKSXV6FtqHpwWgnRFF0y4H/onnjj7BRBcxyzq1Dp9P+H4CIdARGkCrEewCneR8ANSLyFo0i/Kpz7os2xPBdUgW3OQq9D4D9UUGaiJa6bUCnVtwKLG7D9QMlVlqeV3TYN5YW9Bp4cJeho/6KjrBJoE0rQmNXl0P/HvJipeUfoDf3cuDF1oqw1/TwGHqzqwLOtBpco4Eophci5acLDNzZH633GLovjQJ8LHBAM4f+m9SUxPIWPKouST/XjTfeSH5+Pqeeeir77rtvS36OrehKcCj6M4Ueb7rIZOAKl9jaE6Sz5LVq39ihaZiN6Fik6S2cArI7epMaiTY9lFgNrpFM5EQXIFZafjnw3+x8BRckVcA1FWUlt7blm0WkFzrxt2E1fBRNxxOtQ9tdG0T47bScZWd0pbrNwae2tpZbbrmFhx9+mN69e/PCCy8AUFNTw6JFizjwwAPp2rVZ58w6dMT5KW35ebKFN7D0RtTEPIE/75EqNP/9ADBlewNHvaaH59AnhdXASVaDa6QTVdHNqRlpLUFECtFcdvJqeI+0w2rRgZXzgflvvvnm1iOPPPIhtEwJAOccIsLUqVP5/PPPufnmm1m5ciVPPPEEjz32GBs2bODCCy/kqquuai6Mr4Dd/Ph5MkGstHwMumlYTGYqXOLo7+CsirKSeclf8JoenkVrcBcBE6wG12iOqNXpAuAJ2QOEt+kgDkzzS3ABnHO1zrk3nHO3OOe+g6ZX9gUmoVMpFqMr2lHAlcBTf/vb356ura1NmbLR0Jg1Z84cRo0aBcBdd93F2rVree2113jyySepqKjgww8/bC6Mz/z6efwkVlreMVZafje6yuxP5koKO3vnfz5WWn53rLS8I1jTg9E6oriR1sAU4HTCWbP7JXB1Ji/g5XJXeB8zAUSkJ5qSOBY4dsKECWMLCwub3Fg///xzKisrt4nuk08+yeOPP05BQQEHHXQQy5Yto7Kykr322ivlkqighIpYaXkRmvY4jOy9FzqjVQmH5vfsey+6ANg26SGbZWlG9IjkShfAy6udRfhWu3H08TPr9a7OufXOuX84565xzo0fPXp0dXPHvf/++3Ts2JF+/frx1VdfsX79ekaMGLHt6++++y7DhjXxit8M/DNz0bceT3DnoWmXbN98O7vE1iP7nP6rmVLQqRBtejjbBNfYGZEVXQAvrzaD8AhvHN3lnh90IMBQEWk2Yb9y5UqGD9fp5K+88gr777//tq8tXLiQrl270r17d9Ly/QXopl0o8B7tn0crKgLJ7Uteh4L83QbQ/wd3rx7w4z9faV1mRkuItOh6XAYsRDeugqTai+PygONo4Nj0F55//nkmTJjAddddR36+Zpa2bNmSsqqdNWsWo0ePBkgR3Y0bN3YQkYkicpyIhKFq5FY0pRDoZmpeQUfye/Tpk9+jzy1BxmFEh8iLbkVZSQ06lG8pwQlvtXf9k7x4wsB40qYmjx49mh/+8IecffbZzJ07l2uvvZYxY8ZQW1vL9OnTeeihh1i6dCkXXHAB0LjpBjBnzpwCtKX1n8AGEXlTRG4VkTOyPZnBq1I4n/Dk8zsDk2Ol5aODDsQIP5EsGWuOgDZUQFMKC1HB3ZzF6+6MfwGHtuTAxx9/nJkzZ1JVVcWtt97KwQcfnCK4iUSi+sknn3zkjDPO2IJWRxxK0xv2KpK654D3mxmsuMt4dbjL0SqCsFEJDAkin29Eh5wRXdiW57sF7UTKhvDGgenA5SFa4TYwG/ga0MGHc20CxqE3F0SkG9px1VAvfDTQLe17NqJm3Q1C/IZzbsuuBhIrLb8TuIDwrHKTiQMPVJSV/DToQIzwklOi24D3mPcoARTJh4j90NWuHz9/Dep10exQUBHpgBrKJDdu7JV22FYvnoaV8Hzn3JrWBNEem2KM3CPyOd3m8IRwCFo/WY1/doVV3vkeQB8jwyq4oI5ZFwAfoTeJXUl9LGIHU5idc1udc/9yzt3lnPuec24QMBAt6fsTsMA7dATq6fsI8JGIrBaRv4rIJSJymCfeO2Iy4ff/TaD5ZsNolpxc6SbjrY7OB65A22G70LqbTQIV243ATcCMiK1iBNgbXX2OB45D86FxdKNtZ0IXB34LTN2lIESKUP+IhpXwMSS1J3tsotFg/VXgdefcJshNoyOjfZLzotuA90c7HjgZGIt62e7Q4g916JoLzAJeyqE/op5oHnYsOlHiINS3oYDUR/cGn9n9gbV+BuCtaoeRmpLYO+2wBPAuML/HmPM29hj13Z+K5DXrxhMyNgOnVZSVvBB0IEb4aDeim44nwvvQaGbdEc1dJptZt5dfTgFwCCp8JwFHor6zTwMXkyXPBa/0LFmEh+O1qhePv5BuI75JK+0ZgyIB3F5RVhKWmm0jRLRb0TXCj9eEcRRwbP8L77uqoNeA9HREm7n5zEM5dvDuFHctYEvNVt77eD1/eHYp73+y0a9LLKgoKznCr5MZuYOJrhF6vKeSLfhYtfDIhUfz6cZqNlXXc8zgXgzuXcSar6oY/Qff7CXiQNd29LRktJAou4wZ7YfBaMmZb5w97fVt/39g/+6UXzqGfj06k58n1Cd80ckEGvd//DiZkTuY6BpRYD92ULLWViYeM4ghfboxanAvAKa9stIvwQWNdz9MdI00THSNKNCZDEx/nnBQP47eRwW3cn2cBat9rQQUwtk1ZwRMJLaCjXZPIRkQ3bOnvc7Qa//BhTPfZo/unbj73MPZs6dvOiloRYxhpGCia0SBWrR+2hc65ueR50l4TX2COcs+Z0ttPQUd8hi4m2+ulQ4tQTSMFCy9YESBOD6K7vCBPbn97OG8uepLNsTrODK2G907FbBucw2LP97g12Uc4THXN0KEia4RBZbh43v10001rFq3hdFDdqdrYT5fbqnlmXcr+dNLy9lU49t+XT4at2GkYKJrRIEV+JgKW7VuS0rJWIbIQ+M2jBQsp2uEHs/z4oOg42glS6wxwmgOE10jKryCj3ndDJNAjZIMowkmukZUKEdbgaNAFepMZxhNMNE1osKLqN9uFNgIvBR0EEY4MdE1IoGX170J/6aAZIoq4KYc8l42fMZE14gS0wn/ezYPmBF0EEZ4Cfsb2DC24Y1JeoDwNh3EgWkRG+dkZBkTXSNqTEEnMYeRL4Grgw7CCDcmukakqCgrqUKnDIdttRsHzvLiM4ztYqJrRI6KspJ5aN40LMIbB6ZXlJXMDzoQI/yY6BpR5TJgIVAdcBzVXhw2hNJoETYjzYgssdLyImAeMBQf56e1gmpgKTC6oqxkcwDXNyKIrXSNyOIJ3Wh0pZntVEMceAcTXKOVmOgakcYTvK+hNbzZEt64d73jTXCN1mLpBSNniJWWjwYeBYrJzHyyOFqudpa3mWcYrcZWukbO4AnhELSBohr/WoarvPM9AAwxwTV2BVvpGjlJrLS8GDgfuALoDnShdYuMBCq2G1HPhxnWaWb4gYmukdPESsvzgPHAycBYYBgqqPXoxF5BfXodOkklD1iC+uHOAl4y8xrDT0x0jXaFJ8L7APuhed+O6NTeODrTbIVNfDAyiYmuYRhGFrGNNMMwjCxiomsYhpFFTHQNwzCyiImuYRhGFjHRNQzDyCImuoZhGFnERNcwDCOLmOgahmFkERNdwzCMLGKiaxiGkUVMdA3DMLKIia5hGEYWMdE1DMPIIia6hmEYWcRE1zAMI4uY6BqGYWQRE13DMIwsYqJrGIaRRUx0DcMwsoiJrmEYRhb5/+czw0W2lxPWAAAAAElFTkSuQmCC",
"text/plain": [
"