{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "349d0370",
   "metadata": {},
   "source": [
    "# The shortest path problem using QUBO\n",
    "This notebook implements a Perceval code for finding the shortest path on a _directed, weighted graph_ using the Quadratic Unconstrained Binary\n",
    "Optimization (QUBO) model. It is mainly implementing the algorithm in https://arxiv.org/pdf/2112.09766.pdf and represents the work done during the 2022 LOQCathon. \n",
    "\n",
    "Authors: Beata Zjawin, Benjamin Pointard, Nathan Claudet, Noé Delorme, Rina Ismailati.\n",
    "\n",
    "The task is to find the shortest path from start to finish (such that the sum of the weights is minimized) on a simple 5-edge graph:\n",
    "\n",
    "![Graph_5_edges](../_static/img/graph_5edges.png)\n",
    "\n",
    "Start with importing the necessary packages.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ab150578",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "import perceval as pcvl\n",
    "from perceval.components import PS, BS, GenericInterferometer\n",
    "import numpy as np\n",
    "from scipy.optimize import minimize\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcdefaults()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "7c25b79a-d88c-4bf2-9398-93c1434d06cd",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Implementation in Perceval\n",
    "To find the shortest path, first we need to compute the objective function of the weighted graph. The constraints specified by the shortest path problem can be written in the form of quadratic penalty functions, as required by the QUBO formalism [1]:\n",
    "\n",
    "\\begin{align*}\n",
    "H_s&=\\Big(\\sum_{j} x_{s,j} - \\sum_{k} x_{k,s} - 1\\Big)^2, \\\\\n",
    "H_f&=\\Big(\\sum_{j} x_{f,j} - \\sum_{k} x_{k,f} + 1\\Big)^2, \\\\\n",
    "H_i&=\\Big(\\sum_{j} x_{i,j} - \\sum_{k} x_{k,i}\\Big)^2. \\\\\n",
    "\\end{align*}\n",
    "\n",
    "Here, $s$ and $f$ correspond to the start and finish node, respectively. Together with the original objective function:\n",
    "\n",
    "\\begin{equation*}\n",
    "H_c= \\sum_{(i,j)\\in edges} c_{i,j} x^{2}_{i,j},\n",
    "\\end{equation*}\n",
    "\n",
    "the constraints form the final QUBO objective function:\n",
    "\n",
    "\\begin{equation*}\n",
    "H= \\alpha \\Big(H_s + H_f + \\sum_{i \\notin \\{s,f\\}} H_i\\Big) + H_c,\n",
    "\\end{equation*}\n",
    "\n",
    "where $\\alpha$ is a scaling factor satisfying $\\alpha > \\sum_{(i,j)\\in edges} c_{i,j}$.\n",
    "\n",
    "\n",
    "\n",
    "For our simple example, the objective function is given by (see the Appendix to generate Hamiltonians from arbitrary adjacency matrices):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "dcde3aa3-da33-4d7d-9846-7efa216550c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "H1 = [[2., 32., -32., -32.,  32.,   0.],\n",
    "      [0.,  1.,  32.,   0., -32., -32.],\n",
    "      [0.,  0.,  35.,  32., -64., -32.],\n",
    "      [0.,  0.,   0.,   2., -32.,  32.],\n",
    "      [0.,  0.,   0.,   0.,  35.,  32.],\n",
    "      [0.,  0.,   0.,   0.,   0.,   4.]]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "6d0d9f59",
   "metadata": {},
   "source": [
    "Next, we implement a variational algorithm to find the optimal path using the QUBO model.\n",
    "\n",
    "Following the approach introduced in Ref. [2], the output of the simulation, i.e. the Fock states $|n> = |n_1,...,n_M>$, can be mapped into bit strings $b=(b_1,...b_M)$ encoding possible paths on the graph by applying the parity function:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\rho_j: b_i^{(j)} = \\text{mod}[n_i,2] \\oplus j.\n",
    "\\end{equation*}\n",
    "\n",
    "For each evaluation, we need to test for $j=0$,  $j=1$, $n=M$ and $n=M−1$, as it is not possible to know a priori which configuration is the optimal one. Here, the number of modes, $M$, corresponds to the number of possible paths and $n$ is the number of photons.\n",
    "\n",
    "In our example, $b=(b_{sa},b_{sb},b_{ab},b_{af},b_{ba},b_{bf})$, with $b_i \\in \\{0,1\\}$, where subscripts correspond to graph edges.\n",
    "\n",
    "In the cell below, we set up the functions necessary for the simulation, such as implementation of the parity function and sampling of the circuit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3441d4b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def parify_samples(samples, j):\n",
    "    \"\"\"Apply the parity function to the samples.\"\"\"\n",
    "    def _parity(output, j):\n",
    "        m = len(output)\n",
    "        parity = [0] * m\n",
    "        for i in range(m):\n",
    "            parity[i] = (output[i] + j) % 2\n",
    "        return pcvl.BasicState(parity)\n",
    "\n",
    "    new_samples = pcvl.BSCount()\n",
    "    for sample in samples:\n",
    "        new_sample = _parity(sample, j)\n",
    "        if new_sample in new_samples:\n",
    "            new_samples[new_sample] += samples[sample]\n",
    "        else:\n",
    "            new_samples[new_sample] = samples[sample]\n",
    "    return new_samples\n",
    "\n",
    "def set_parameters_circuit(parameters_circuit, values):\n",
    "    \"\"\"Set values of circuit parameters.\"\"\"\n",
    "    for idx, p in enumerate(parameters_circuit):\n",
    "        parameters_circuit[idx].set_value(values[idx])\n",
    "\n",
    "def compute_samples(circuit, input_state, nb_samples, j):\n",
    "    \"\"\"Sample from the circuit.\"\"\"\n",
    "    experiment = pcvl.Experiment(circuit)\n",
    "    experiment.with_input(input_state)\n",
    "    experiment.min_detected_photons_filter(0)\n",
    "    computer = pcvl.SimulatedComputer(\"SLOS\")\n",
    "    factory = pcvl.ExecutionFactory(computer, experiment)\n",
    "    with computer.acquire():\n",
    "        samples = factory.sample_count(nb_samples)['results']\n",
    "    return parify_samples(samples, j)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "bc3125ec-1fe5-4df4-a5bf-0cb9a7de5df0",
   "metadata": {},
   "source": [
    "The central part of solving the shortest path problem with Perceval is initializing a universal circuit and optimizing its parameters. We will work with the following generic interferometer:\n",
    "\n",
    "![generic_6photons_interferometer](../_static/img/generic_6photons_interferometer.svg)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "2bb3985b-3105-4d5b-893f-8d89f41d0051",
   "metadata": {},
   "source": [
    "We sample the circuit for all configurations of $j$ and $n$. In this example, we set all the initial parameters to $\\pi$, which is close to the optimal solution $-$ a good initial guess can dramatically improve simulation time. The loss function is given by\n",
    "\n",
    "$$E(j,\\theta,\\psi) = \\sum_{b^{(j)}} \\beta_{b^{(j)}} <b_1^{(j)},...,b_M^{(j)}|H|b_1^{(j)},...,b_M^{(j)}>$$\n",
    "\n",
    "where $\\theta$ and $\\psi$ stand for the parameters of the linear interferometer: beam-splitters and phase-shifters, respectively. For optimization, we use the Powell minimisation algorithm from the scipy.optimize package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d0d7c887",
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_configuration(circuit, nb_modes, j, n, H, nb_samples):\n",
    "    \"\"\"output the samples for a given configuration of j and n (includes minimisation of the loss function)\"\"\"\n",
    "    parameters_circuit = circuit.get_parameters()\n",
    "\n",
    "    input_state = pcvl.BasicState([1]*nb_modes)\n",
    "    if n!=nb_modes:\n",
    "        input_state = pcvl.BasicState([1]*(nb_modes-1)+[0])\n",
    "\n",
    "    def loss(parameters):\n",
    "        set_parameters_circuit(parameters_circuit, parameters)\n",
    "        samples = compute_samples(circuit, input_state, nb_samples, j)\n",
    "        E = 0\n",
    "        for sample in samples:\n",
    "            b = np.array([sample[i] for i in range(len(sample))])\n",
    "            b_prime = np.dot(H, b)\n",
    "            E += samples[sample]/nb_samples*np.dot(b.conjugate(), b_prime)\n",
    "            \n",
    "        return E.real\n",
    "\n",
    "    # init_parameters = [2*(math.pi)*random.random() for _ in parameters_circuit] # initialize with random initial parameters\n",
    "    init_parameters = [math.pi for _ in parameters_circuit] # initialize with a good guess\n",
    "    best_parameters = minimize(loss, init_parameters, method='Powell', bounds=[(0,2*math.pi) for _ in init_parameters]).x    \n",
    "    set_parameters_circuit(parameters_circuit, best_parameters)\n",
    "    samples = compute_samples(circuit, input_state, nb_samples, j)\n",
    "    return samples\n",
    "\n",
    "def shortest_path(H, nb_samples):\n",
    "    \"\"\"run the universal circuit and optimize the parameters\"\"\"\n",
    "    nb_modes = len(H)\n",
    "    circuit = GenericInterferometer(\n",
    "        nb_modes,\n",
    "        lambda i: BS(theta=pcvl.P(f\"theta{i}\"),phi_tr=pcvl.P(f\"phi_tr{i}\")),\n",
    "        phase_shifter_fun_gen=lambda i: PS(phi=pcvl.P(f\"phi{i}\")))\n",
    "    js = [0,1]\n",
    "    ns = [nb_modes, nb_modes-1]\n",
    "    configuration_samples = []\n",
    "    for j in js:\n",
    "        for n in ns:\n",
    "            current_sample = test_configuration(circuit, nb_modes, j, n, H, nb_samples)\n",
    "            print(f\"Configuration for (j,n)=({j},{n})\")\n",
    "            print(current_sample)\n",
    "            configuration_samples.append(([j,n],current_sample))\n",
    "            \n",
    "    return configuration_samples"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "75ab2dba-e47b-4d53-a92f-c13c077aeb69",
   "metadata": {},
   "source": [
    "Run the algorithm on 10000 samples (this may take a few minutes):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7f440a3e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Configuration for (j,n)=(0,6)\n",
      "{\n",
      "  |1,1,0,0,0,0>: 721\n",
      "  |1,1,1,0,0,1>: 14\n",
      "  |0,0,0,0,0,0>: 2637\n",
      "  |1,0,1,0,0,0>: 50\n",
      "  |1,0,0,0,1,0>: 503\n",
      "  |0,1,1,1,0,1>: 6\n",
      "  |1,0,0,1,0,0>: 693\n",
      "  |1,0,0,0,0,1>: 1053\n",
      "  |0,1,1,0,0,0>: 18\n",
      "  |0,0,0,1,0,1>: 595\n",
      "  |1,1,1,1,0,0>: 8\n",
      "  |0,1,1,1,1,0>: 8\n",
      "  |0,1,0,0,1,0>: 583\n",
      "  |1,1,0,1,0,1>: 279\n",
      "  |0,0,1,1,1,1>: 2\n",
      "  |0,0,1,0,0,1>: 18\n",
      "  |0,1,0,1,0,0>: 624\n",
      "  |0,1,0,0,0,1>: 521\n",
      "  |0,0,0,0,1,1>: 489\n",
      "  |0,0,1,1,0,0>: 23\n",
      "  |0,0,1,0,1,0>: 24\n",
      "  |0,0,0,1,1,0>: 313\n",
      "  |1,1,0,1,1,0>: 281\n",
      "  |1,0,0,1,1,1>: 49\n",
      "  |1,1,0,0,1,1>: 78\n",
      "  |1,0,1,1,1,0>: 5\n",
      "  |0,1,0,1,1,1>: 383\n",
      "  |1,1,1,0,1,0>: 7\n",
      "  |1,0,1,0,1,1>: 2\n",
      "  |1,1,1,1,1,1>: 3\n",
      "  |1,0,1,1,0,1>: 7\n",
      "  |0,1,1,0,1,1>: 3\n",
      "}\n",
      "Configuration for (j,n)=(0,5)\n",
      "{\n",
      "  |1,1,1,0,0,0>: 2\n",
      "  |1,1,1,1,1,0>: 1\n",
      "  |0,1,0,0,0,0>: 6575\n",
      "  |1,1,0,1,0,0>: 4\n",
      "  |1,1,0,0,0,1>: 5\n",
      "  |0,0,1,0,0,0>: 2\n",
      "  |1,1,0,0,1,0>: 8\n",
      "  |0,1,0,1,1,0>: 3397\n",
      "  |0,0,0,1,0,0>: 1\n",
      "  |0,1,1,0,0,1>: 3\n",
      "  |0,1,0,0,1,1>: 1\n",
      "  |0,1,1,0,1,0>: 1\n",
      "}\n",
      "Configuration for (j,n)=(1,6)\n",
      "{\n",
      "  |0,0,0,1,1,0>: 4\n",
      "  |1,0,0,0,0,1>: 11\n",
      "  |1,1,1,1,0,0>: 1\n",
      "  |1,1,0,1,1,0>: 1\n",
      "  |1,0,0,1,1,1>: 1\n",
      "  |1,0,0,1,0,0>: 9972\n",
      "  |1,0,1,1,0,1>: 6\n",
      "  |1,1,0,0,0,0>: 2\n",
      "  |1,0,1,0,0,0>: 1\n",
      "  |0,0,1,1,0,0>: 1\n",
      "}\n",
      "Configuration for (j,n)=(1,5)\n",
      "{\n",
      "  |1,0,0,1,1,0>: 1\n",
      "  |0,1,0,1,0,1>: 1\n",
      "  |1,1,0,1,0,0>: 2\n",
      "  |0,1,0,0,1,1>: 1\n",
      "  |0,1,0,1,1,0>: 9995\n",
      "}\n"
     ]
    }
   ],
   "source": [
    "nb_samples = 10000\n",
    "data = shortest_path(H1, nb_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "308e2867",
   "metadata": {},
   "source": [
    "Plot the results for better visualization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f29aa6de",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_samples(samples):\n",
    "    \n",
    "    values = list(samples.values())\n",
    "\n",
    "    keys = samples.keys()\n",
    "\n",
    "    key_list = []\n",
    "    for x in keys:\n",
    "        s = \"\".join(str(c) for c in list(x))\n",
    "        key_list.append(s)\n",
    "\n",
    "    y_pos = np.arange(len(key_list))\n",
    "    barlist = plt.bar(y_pos, values, align='center', alpha=0.8)\n",
    "    index = key_list.index('100100')\n",
    "    barlist[index].set_color('m')\n",
    "    plt.yscale('log')\n",
    "    plt.xticks(y_pos, key_list)\n",
    "    plt.xticks(rotation = 70)\n",
    "    plt.ylabel('Sample values')\n",
    "    plt.title('Sampling the Fock States')\n",
    "\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f13d9a95",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAHcCAYAAADfvQDCAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAAV71JREFUeJzt3Xl4TGf/BvB7spMQS4gtEjuxJJaIpSQ0pLYWVaq0pIpWLJXS0r5tqtVq67W1YqkiVWqtndIKGku8CLHUGiKNEBEhkYQsM9/fH70yP2OCRJYzc3J/rst1dZ4ZM/ccp5k7Z57zHI2ICIiIiIhUyELpAERERETFhUWHiIiIVItFh4iIiFSLRYeIiIhUi0WHiIiIVItFh4iIiFSLRYeIiIhUi0WHiIiIVItFh4iIiFSLRYeoFNNoNPj888/1t0NDQ6HRaHDt2jXFMj3u2rVr0Gg0+O9//6t0lCKn0WgwduxYpWMQqRqLDlEhnTlzBgMGDICrqyvs7OxQs2ZNdOvWDT/88IPS0czKzp07DUqXEjQaTZ5/qlWrpmiuvGRlZWHevHlo2bIlypcvjwoVKqBp06YYNWoULly4oH/c4cOH8fnnn+PevXvP/VoLFixAaGho4UMTKcBK6QBE5uzw4cPo0qULateujZEjR6JatWqIi4vDkSNHMG/ePIwbN07piAXy5ptv4vXXX4etrW2Jv/bOnTsREhKieNnp1q0b3nrrLYOxMmXKKJTmyV599VX8/vvvGDx4MEaOHIns7GxcuHAB27dvR4cOHdC4cWMA/+6j06ZNw/Dhw1GhQoXneq0FCxbAyckJw4cPL7o3QFRCWHSICuGrr76Co6Mjjh07ZvQhkpiYqEyoQrC0tISlpaXSMRTVsGFDDB06VOkYT3Xs2DFs374dX331FT7++GOD++bPn1+oozdEasOvrogK4cqVK2jatGmevylXrVrV4Pby5cvRtWtXVK1aFba2tnB3d8fChQuN/p6bmxt69+6N/fv3o02bNihTpgyaN2+O/fv3AwA2btyI5s2bw87ODq1bt8bJkycN/v7w4cPh4OCAq1evwt/fH/b29qhRowa++OILiMhT309ec3Ry8xw8eBBt27aFnZ0d6tatixUrVhj9/dOnT8PHxwdlypRBrVq1MH36dCxfvvyZ836GDx+OkJAQAIZfHz3uxx9/RL169WBrawsvLy8cO3bM6DEXLlzAgAEDUKlSJdjZ2aFNmzbYunXrU993QSQmJmLEiBFwdnaGnZ0dPDw88PPPPxs9TqfTYd68efp/qypVquCll17C8ePHn/r806dPh4WFxVO/+rxy5QoAoGPHjkb3WVpaonLlygCAzz//HJMnTwYA1KlTR79dc/8t8rNPurm54e+//8Zff/2l//u+vr76++/du4f3338fLi4usLW1Rf369fHtt99Cp9MZPM+aNWvQunVrlCtXDuXLl0fz5s0xb968p24LoqLAIzpEheDq6oqIiAicPXsWzZo1e+pjFy5ciKZNm+Lll1+GlZUVtm3bhjFjxkCn0yEwMNDgsdHR0XjjjTcwevRoDB06FP/973/Rp08fLFq0CB9//DHGjBkDAJgxYwYGDhyIixcvwsLi/39v0Wq1eOmll9CuXTt899132LVrF4KDg5GTk4MvvviiwO8zOjoaAwYMwIgRIzBs2DAsW7YMw4cPR+vWrdG0aVMAQHx8PLp06QKNRoOpU6fC3t4eP/30U76+Bhs9ejRu3LiBP//8E7/88kuej/n1119x//59jB49GhqNBt999x369++Pq1evwtraGgDw999/o2PHjqhZsyamTJkCe3t7rFu3Dn379sVvv/2Gfv36PTPLw4cPkZSUZDBWrlw52Nra4sGDB/D19UV0dDTGjh2LOnXqYP369Rg+fDju3buHCRMm6P/OiBEjEBoaih49euCdd95BTk4ODhw4gCNHjqBNmzZ5vvZ//vMffP3111i8eDFGjhz5xIyurq4AgFWrVqFjx46wssr7R3n//v1x6dIlrF69GnPmzIGTkxMAoEqVKgDyt0/OnTsX48aNg4ODAz755BMAgLOzMwAgIyMDPj4+iI+Px+jRo1G7dm0cPnwYU6dOxc2bNzF37lwAwJ9//onBgwfjxRdfxLfffgsAOH/+PA4dOmSwzYiKhRDRc/vjjz/E0tJSLC0tpX379vLhhx/K7t27JSsry+ixGRkZRmP+/v5St25dgzFXV1cBIIcPH9aP7d69WwBImTJlJDY2Vj++ePFiASD79u3Tjw0bNkwAyLhx4/RjOp1OevXqJTY2NnL79m39OAAJDg7W316+fLkAkJiYGKM84eHh+rHExESxtbWVDz74QD82btw40Wg0cvLkSf3YnTt3pFKlSkbPmZfAwEDJ60dSTEyMAJDKlStLcnKyfnzLli0CQLZt26Yfe/HFF6V58+by8OFDg/feoUMHadCgwVNfX+Tf7ZHXn+XLl4uIyNy5cwWArFy5Uv93srKypH379uLg4CCpqakiIrJ3714BIOPHjzd6DZ1OZ/B6gYGBIiLywQcfiIWFhYSGhj4zp06nEx8fHwEgzs7OMnjwYAkJCTHYN3LNnDnzids/v/tk06ZNxcfHx+ixX375pdjb28ulS5cMxqdMmSKWlpbyzz//iIjIhAkTpHz58pKTk/PM90ZU1PjVFVEhdOvWDREREXj55Zdx6tQpfPfdd/D390fNmjWNvi55dEJrSkoKkpKS4OPjg6tXryIlJcXgse7u7mjfvr3+tre3NwCga9euqF27ttH41atXjbI9etpy7mnMWVlZ2LNnT4Hfp7u7Ozp16qS/XaVKFTRq1MjgdXft2oX27dvD09NTP1apUiUMGTKkwK+Xl0GDBqFixYr627l5cjMkJydj7969GDhwIO7fv4+kpCQkJSXhzp078Pf3x+XLlxEfH//M13nllVfw559/Gvzx9/cH8O+E6WrVqmHw4MH6x1tbW2P8+PFIS0vDX3/9BQD47bffoNFoEBwcbPT8j38lJyIYO3Ys5s2bh5UrV2LYsGHPzKjRaLB7925Mnz4dFStWxOrVqxEYGAhXV1cMGjQo33N0CrJP5mX9+vXo1KkTKlasqN/eSUlJ8PPzg1arRXh4OACgQoUKSE9Px59//pmvXERFiV9dERWSl5cXNm7ciKysLJw6dQqbNm3CnDlzMGDAAERFRcHd3R0AcOjQIQQHByMiIgIZGRkGz5GSkgJHR0f97UfLDAD9fS4uLnmO371712DcwsICdevWNRhr2LAhADzXGjmP5wGAihUrGrxubGysQTnLVb9+/QK/Xn4y5Jae3AzR0dEQEXz66af49NNP83yOxMRE1KxZ86mvU6tWLfj5+eV5X2xsLBo0aGDwNSEANGnSRH8/8O8cmho1aqBSpUrPeFfAihUrkJaWhoULFxoUqGextbXFJ598gk8++QQ3b97EX3/9hXnz5mHdunWwtrbGypUrn/kcBdkn83L58mWcPn1a/1XY43In5I8ZMwbr1q1Djx49ULNmTXTv3h0DBw7ESy+9lM93S/T8WHSIioiNjQ28vLzg5eWFhg0bIiAgAOvXr0dwcDCuXLmCF198EY0bN8bs2bPh4uICGxsb7Ny5E3PmzDGauPmkM5+eNC7PmGRcWEq9bkEy5G7DSZMm6Y/APK6oSldR6tixI6KiojB//nwMHDgwX+XocdWrV8frr7+OV199FU2bNsW6desQGhr6xLk7AAq8T+ZFp9OhW7du+PDDD/O8P7dcV61aFVFRUdi9ezd+//13/P7771i+fDneeuutPCdyExUlFh2iYpA72fTmzZsAgG3btiEzMxNbt241ODKxb9++Ynl9nU6Hq1ev6j9oAODSpUsA/j2Lpji4uroiOjraaDyvsbzkdZZVQeQewbK2tn7iEZnCcnV1xenTp6HT6QyO6uQu0Jc7SbhevXrYvXs3kpOTn1lc6tevj++++w6+vr546aWXEBYWhnLlyj1XPmtra7Ro0QKXL19GUlISqlWr9sTtWpB98knPUa9ePaSlpeVre9vY2KBPnz7o06cPdDodxowZg8WLF+PTTz81yQJK6sE5OkSFsG/fvjyPauzcuRMA0KhRIwD/fzTi0cempKRg+fLlxZZt/vz5+v8WEcyfPx/W1tZ48cUXi+X1/P39ERERgaioKP1YcnIyVq1ala+/b29vDwDPvQZM1apV4evri8WLF+sL5qNu3779XM/7qJ49eyIhIQFr167Vj+Xk5OCHH36Ag4MDfHx8APy7mJ+IYNq0aUbPkdf+0qJFC+zcuRPnz59Hnz598ODBg6fmuHz5Mv755x+j8Xv37iEiIgIVK1bUf530pO1akH3S3t4+z3+XgQMHIiIiArt3784zS05ODgDgzp07BvdZWFigRYsWAIDMzMwnvU2iIsEjOkSFMG7cOGRkZKBfv35o3LgxsrKycPjwYaxduxZubm4ICAgAAHTv3l3/G+3o0aORlpaGJUuWoGrVqnl+KBeWnZ0ddu3ahWHDhsHb2xu///47duzYgY8//viJ8ykK68MPP8TKlSvRrVs3jBs3Tn96ee3atZGcnPzMIzatW7cGAIwfPx7+/v6wtLTE66+/XqAMISEheOGFF9C8eXOMHDkSdevWxa1btxAREYHr16/j1KlTz/3+AGDUqFFYvHgxhg8fjsjISLi5uWHDhg04dOgQ5s6dqz8S06VLF7z55pv4/vvvcfnyZbz00kvQ6XQ4cOAAunTpkuf1rdq1a4ctW7agZ8+eGDBgADZv3qw/bf5xp06dwhtvvIEePXqgU6dOqFSpEuLj4/Hzzz/jxo0bmDt3rr7I5G7XTz75BK+//jqsra3Rp0+fAu2TrVu3xsKFCzF9+nTUr18fVatWRdeuXTF58mRs3boVvXv31i83kJ6ejjNnzmDDhg24du0anJyc8M477yA5ORldu3ZFrVq1EBsbix9++AGenp76+U1ExUahs72IVOH333+Xt99+Wxo3biwODg5iY2Mj9evXl3HjxsmtW7cMHrt161Zp0aKF2NnZiZubm3z77beybNmyPE/n7tWrl9Fr4ZFTkXPlnno9c+ZM/diwYcPE3t5erly5It27d5eyZcuKs7OzBAcHi1arNXrO/JxenlceHx8fo1OOT548KZ06dRJbW1upVauWzJgxQ77//nsBIAkJCU/ajCIikpOTI+PGjZMqVaqIRqPRn2qe13t8Un4RkStXrshbb70l1apVE2tra6lZs6b07t1bNmzY8NTXz32+x7fx427duiUBAQHi5OQkNjY20rx5c/3p54+/n5kzZ0rjxo3FxsZGqlSpIj169JDIyMinvt6WLVvEyspKBg0aZPTv9WiGb775Rnx8fKR69epiZWUlFStWlK5du+b5Pr/88kupWbOmWFhYGPz75nefTEhIkF69ekm5cuUEgMG/+/3792Xq1KlSv359sbGxEScnJ+nQoYP897//1S+zsGHDBunevbtUrVpVbGxspHbt2jJ69Gi5efPmU7c1UVHQiJTgbEIiKnbDhw/Hhg0bkJaWpnQUAMD777+PxYsXIy0trdRfXoKISh7n6BBRkXl8bsmdO3fwyy+/4IUXXmDJISJFcI4OERWZ9u3bw9fXF02aNMGtW7ewdOlSpKamPnFdGyKi4saiQ0RFpmfPntiwYQN+/PFHaDQatGrVCkuXLkXnzp2VjkZEpRTn6BAREZFqcY4OERERqRaLDhEREalWqZ+jo9PpcOPGDZQrV67QS9ATERFRyRAR3L9/HzVq1DC60O6jSn3RuXHjhtEVoYmIiMg8xMXFoVatWk+8v9QXndwl2+Pi4lC+fHmF0xAREVF+pKamwsXF5ZkXwS31RSf366ry5cuz6BAREZmZZ0074WRkIiIiUi0WHSIiIlItFh0iIiJSLRYdIiIiUi0WHSIiIlItFh0iIiJSLRYdIiIiUi0WHSIiIlIt1RSdjIwMuLq6YtKkSUpHISIiIhOhmqLz1VdfoV27dkrHICIiIhOiiqJz+fJlXLhwAT169FA6ChEREZkQxYtOeHg4+vTpgxo1akCj0WDz5s1GjwkJCYGbmxvs7Ozg7e2No0ePGtw/adIkzJgxo4QSExERkblQvOikp6fDw8MDISEhed6/du1aBAUFITg4GCdOnICHhwf8/f2RmJgIANiyZQsaNmyIhg0blmRsIiIiMgMaERGlQ+TSaDTYtGkT+vbtqx/z9vaGl5cX5s+fDwDQ6XRwcXHBuHHjMGXKFEydOhUrV66EpaUl0tLSkJ2djQ8++ACfffZZnq+RmZmJzMxM/e3cy7ynpKTw6uVEpUBmQiZy7uUoHcOAVQUr2FazVToGkVlJTU2Fo6PjMz+/rUowU4FlZWUhMjISU6dO1Y9ZWFjAz88PERERAIAZM2bov7YKDQ3F2bNnn1hych8/bdq04g1ORCYpMyETp186jZy7JlZ0Klqhxa4WLDtExcCki05SUhK0Wi2cnZ0Nxp2dnXHhwoXnes6pU6ciKChIfzv3iA4RqV/OvRzk3M2BxkYDC1vFv7kHAOgydci5m4OcezksOkTFwKSLTkENHz78mY+xtbWFrS1/mBCVZha2FrCwM42iAwDaLK3SEYhUy3T+T8+Dk5MTLC0tcevWLYPxW7duoVq1agqlIiIiInNh0kXHxsYGrVu3RlhYmH5Mp9MhLCwM7du3L9Rzh4SEwN3dHV5eXoWNSURERCZK8a+u0tLSEB0drb8dExODqKgoVKpUCbVr10ZQUBCGDRuGNm3aoG3btpg7dy7S09MREBBQqNcNDAxEYGCgftY2ERERqY/iRef48ePo0qWL/nbuROFhw4YhNDQUgwYNwu3bt/HZZ58hISEBnp6e2LVrl9EEZSIiIqLHKV50fH198aylfMaOHYuxY8eWUCIiIiJSC5Oeo1OcOEeHiIhI/Upt0QkMDMS5c+dw7NgxpaMQERFRMSm1RYeIiIjUj0WHiIiIVItFh4iIiFSr1BYdTkYmIiJSv1JbdDgZmYiISP1KbdEhIiIi9WPRISIiItVi0SEiIiLVKrVFh5ORiYiI1K/UFh1ORiYiIlK/Ult0iIiISP1YdIiIiEi1WHSIiIhItVh0iIiISLVYdIiIiEi1Sm3R4enlRERE6ldqiw5PLyciIlK/Ult0iIiISP1YdIiIiEi1WHSIiIhItVh0iIiISLVYdIiIiEi1WHSIiIhItUpt0eE6OkREROpXaosO19EhIiJSv1JbdIiIiEj9WHSIiIhItVh0iIiISLVYdIiIiEi1WHSIiIhItVh0iIiISLVYdIiIiEi1WHSIiIhItVh0iIiISLVKbdHhJSCIiIjUr9QWHV4CgoiISP1KbdEhIiIi9WPRISIiItVi0SEiIiLVYtEhIiIi1WLRISIiItVi0SEiIiLVYtEhIiIi1WLRISIiItVi0SEiIiLVYtEhIiIi1WLRISIiItVi0SEiIiLVYtEhIiIi1Sq1RSckJATu7u7w8vJSOgoREREVk1JbdAIDA3Hu3DkcO3ZM6ShERERUTEpt0SEiIiL1Y9EhIiIi1WLRISIiItVi0SEiIiLVYtEhIiIi1WLRISIiItVi0SEiIiLVYtEhIiIi1WLRISIiItVi0SEiIiLVYtEhIiIi1WLRISIiItVi0SEiIiLVYtEhIiIi1WLRISIiItVi0SEiIiLVYtEhIiIi1TL7onPv3j20adMGnp6eaNasGZYsWaJ0JCIiIjIRVkoHKKxy5cohPDwcZcuWRXp6Opo1a4b+/fujcuXKSkcjIiIihZn9ER1LS0uULVsWAJCZmQkRgYgonIqIiIhMgeJFJzw8HH369EGNGjWg0WiwefNmo8eEhITAzc0NdnZ28Pb2xtGjRw3uv3fvHjw8PFCrVi1MnjwZTk5OJZSeiIiITJniRSc9PR0eHh4ICQnJ8/61a9ciKCgIwcHBOHHiBDw8PODv74/ExET9YypUqIBTp04hJiYGv/76K27dulVS8YmIiMiEKV50evTogenTp6Nfv3553j979myMHDkSAQEBcHd3x6JFi1C2bFksW7bM6LHOzs7w8PDAgQMHnvh6mZmZSE1NNfhDRERE6qR40XmarKwsREZGws/PTz9mYWEBPz8/REREAABu3bqF+/fvAwBSUlIQHh6ORo0aPfE5Z8yYAUdHR/0fFxeX4n0TREREpBiTLjpJSUnQarVwdnY2GHd2dkZCQgIAIDY2Fp06dYKHhwc6deqEcePGoXnz5k98zqlTpyIlJUX/Jy4urljfAxERESnH7E8vb9u2LaKiovL9eFtbW9ja2hZfICIiIjIZJn1Ex8nJCZaWlkaTi2/duoVq1aoplIqIiIjMhUkXHRsbG7Ru3RphYWH6MZ1Oh7CwMLRv375Qzx0SEgJ3d3d4eXkVNiYRERGZKMW/ukpLS0N0dLT+dkxMDKKiolCpUiXUrl0bQUFBGDZsGNq0aYO2bdti7ty5SE9PR0BAQKFeNzAwEIGBgUhNTYWjo2Nh3wYRERGZIMWLzvHjx9GlSxf97aCgIADAsGHDEBoaikGDBuH27dv47LPPkJCQAE9PT+zatctogjIRERHR4xQvOr6+vs+8ZMPYsWMxduzYEkpEREREamHSc3SKE+foEBERqV+pLTqBgYE4d+4cjh07pnQUIiIiKialtugQERGR+rHoEBERkWqx6BAREZFqldqiw8nIRERE6ldqiw4nIxMREalfqS06REREpH4sOkRERKRaLDpERESkWqW26HAyMhERkfqV2qLDychERETqV2qLDhEREakfiw4RERGpFosOERERqRaLDhEREakWiw4RERGpVqktOjy9nIiISP1KbdHh6eVERETqV2qLDhEREakfiw4RERGpFosOERERqRaLDhEREakWiw4RERGpFosOERERqVapLTpcR4eIiEj9Sm3R4To6RERE6lfgovPzzz9jx44d+tsffvghKlSogA4dOiA2NrZIwxEREREVRoGLztdff40yZcoAACIiIhASEoLvvvsOTk5OmDhxYpEHJCIiInpeVgX9C3Fxcahfvz4AYPPmzXj11VcxatQodOzYEb6+vkWdj4iIiOi5FfiIjoODA+7cuQMA+OOPP9CtWzcAgJ2dHR48eFC06YiIiIgKocBHdLp164Z33nkHLVu2xKVLl9CzZ08AwN9//w03N7eizkdERET03Ap8RCckJATt27fH7du38dtvv6Fy5coAgMjISAwePLjIAxIRERE9rwIf0alQoQLmz59vND5t2rQiCURERERUVJ5rHZ0DBw5g6NCh6NChA+Lj4wEAv/zyCw4ePFik4YiIiIgKo8BF57fffoO/vz/KlCmDEydOIDMzEwCQkpKCr7/+usgDEhERET2vAhed6dOnY9GiRViyZAmsra314x07dsSJEyeKNFxx4iUgiIiI1K/ARefixYvo3Lmz0bijoyPu3btXFJlKBC8BQUREpH4FLjrVqlVDdHS00fjBgwdRt27dIglFREREVBQKXHRGjhyJCRMm4H//+x80Gg1u3LiBVatWYdKkSXjvvfeKIyMRERHRcynw6eVTpkyBTqfDiy++iIyMDHTu3Bm2traYNGkSxo0bVxwZiYiIiJ5LgYuORqPBJ598gsmTJyM6OhppaWlwd3eHg4NDceQjIiIiem4FLjq5bGxs4O7uXpRZiIiIiIpUgYtOly5doNFonnj/3r17CxWIiIiIqKgUuOh4enoa3M7OzkZUVBTOnj2LYcOGFVUuIiIiokIrcNGZM2dOnuOff/450tLSCh2IiIiIqKg817Wu8jJ06FAsW7asqJ6OiIiIqNCKrOhERETAzs6uqJ6OiIiIqNAK/NVV//79DW6LCG7evInjx4/j008/LbJgRERERIVV4KLj6OhocNvCwgKNGjXCF198ge7duxdZMCIiIqLCKnDRWb58eXHkICIiIipyRTZHx9yEhITA3d0dXl5eSkchIiKiYpKvIzoVK1Z86iKBj0pOTi5UoJISGBiIwMBApKamGn0dR0REROqQr6Izd+7cYo5BREREVPTyVXS44jERERGZo+e+qCcAPHz4EFlZWQZj5cuXL1QgIiIioqJS4MnI6enpGDt2LKpWrQp7e3tUrFjR4A8RERGRqShw0fnwww+xd+9eLFy4ELa2tvjpp58wbdo01KhRAytWrCiOjERERETPpcBfXW3btg0rVqyAr68vAgIC0KlTJ9SvXx+urq5YtWoVhgwZUhw5iYiIiAqswEd0kpOTUbduXQD/zsfJPZ38hRdeQHh4eNGmIyIiIiqEAhedunXrIiYmBgDQuHFjrFu3DsC/R3oqVKhQpOGIiIiICqPARScgIACnTp0CAEyZMgUhISGws7PDxIkTMXny5CIPSERERPS8CjxHZ+LEifr/9vPzw4ULFxAZGYn69eujRYsWRRqOiIiIqDAKXHTi4uLg4uKiv+3q6gpXV9ciDUVERERUFAr81ZWbmxt8fHywZMkS3L17tzgyERERERWJAhed48ePo23btvjiiy9QvXp19O3bFxs2bEBmZmZx5CMiIiJ6bgUuOi1btsTMmTPxzz//4Pfff0eVKlUwatQoODs74+233y6OjERERETPpcBFJ5dGo0GXLl2wZMkS7NmzB3Xq1MHPP/9clNmIiIiICuW5i87169fx3XffwdPTE23btoWDgwNCQkKKMhsRERFRoRT4rKvFixfj119/xaFDh9C4cWMMGTIEW7Zs4ZlXREREZHIKfERn+vTp8Pb2RmRkJM6ePYupU6cqWnLi4uLg6+sLd3d3tGjRAuvXr1csCxEREZmWAh/R+eeff6DRaIojy3OxsrLC3Llz4enpiYSEBLRu3Ro9e/aEvb290tGIiIhIYQUuOqZUcgCgevXqqF69OgCgWrVqcHJyQnJyMosOERERPf9k5KISHh6OPn36oEaNGtBoNNi8ebPRY0JCQuDm5gY7Ozt4e3vj6NGjeT5XZGQktFqtwcrNREREVHopXnTS09Ph4eHxxDO21q5di6CgIAQHB+PEiRPw8PCAv78/EhMTDR6XnJyMt956Cz/++GNJxCYiIiIzUOCvropajx490KNHjyfeP3v2bIwcORIBAQEAgEWLFmHHjh1YtmwZpkyZAgDIzMxE3759MWXKFHTo0OGpr5eZmWmwinNqamoRvAsiIiIyRc91RCcnJwd79uzB4sWLcf/+fQDAjRs3kJaWVqThsrKyEBkZCT8/P/2YhYUF/Pz8EBERAQAQEQwfPhxdu3bFm2+++cznnDFjBhwdHfV/+DUXERGRehW46MTGxqJ58+Z45ZVXEBgYiNu3bwMAvv32W0yaNKlIwyUlJUGr1cLZ2dlg3NnZGQkJCQCAQ4cOYe3atdi8eTM8PT3h6emJM2fOPPE5p06dipSUFP2fuLi4Is1MREREpqPAX11NmDABbdq0walTp1C5cmX9eL9+/TBy5MgiDZcfL7zwAnQ6Xb4fb2trC1tb22JMRERERKaiwEXnwIEDOHz4MGxsbAzG3dzcEB8fX2TBAMDJyQmWlpa4deuWwfitW7dQrVq1In0tIiIiUp8Cf3Wl0+mg1WqNxq9fv45y5coVSahcNjY2aN26NcLCwgxePywsDO3bty/Uc4eEhMDd3R1eXl6FjUlEREQmqsBFp3v37pg7d67+tkajQVpaGoKDg9GzZ88CB0hLS0NUVBSioqIAADExMYiKisI///wDAAgKCsKSJUvw888/4/z583jvvfeQnp6uPwvreQUGBuLcuXM4duxYoZ6HiIiITFeBv7qaNWsW/P394e7ujocPH+KNN97A5cuX4eTkhNWrVxc4wPHjx9GlSxf97aCgIADAsGHDEBoaikGDBuH27dv47LPPkJCQAE9PT+zatctogjIRERHR4zQiIgX9Szk5OVizZg1Onz6NtLQ0tGrVCkOGDEGZMmWKI2OxSk1NhaOjI1JSUlC+fHml4xBRMUq/kI7T/qdhWc4SFnaKr5cKANA91EF7X4sWu1vAvjEvXUOUX/n9/H6uBQOtrKwwdOjQ5w5nCkJCQhASEpLnfCMiIiJSh3wVna1bt+b7CV9++eXnDlOSAgMDERgYqG+EREREpD75Kjp9+/bN15NpNBoeISEiIiKTka+iU5AF+YiIiIhMhWnMxiMiIiIqBs9VdMLCwtC7d2/Uq1cP9erVQ+/evbFnz56izlasuGAgERGR+hW46CxYsAAvvfQSypUrhwkTJmDChAkoX748evbsiZCQkOLIWCy4YCAREZH6Ffj08q+//hpz5szB2LFj9WPjx49Hx44d8fXXXyMwMLBIAxIRERE9rwIf0bl37x5eeuklo/Hu3bsjJSWlSEIRERERFYUCF52XX34ZmzZtMhrfsmULevfuXSShiIiIiIpCgb+6cnd3x1dffYX9+/frryB+5MgRHDp0CB988AG+//57/WPHjx9fdEmLGFdGJiIiUr8CX+uqTp06+XtijQZXr159rlAlide6Iio9eK0rIvUotmtdxcTEFCoYERERUUkxjV9piIiIiIpBgY/oiAg2bNiAffv2ITEx0ejyEBs3biyycERERESFUeCi8/7772Px4sXo0qULnJ2dodFoiiMXERERUaEVuOj88ssv2LhxI3r27FkceYiIiIiKTIHn6Dg6OqJu3brFkaVE8VpXRERE6lfgovP5559j2rRpePDgQXHkKTG81hUREZH6Ffirq4EDB2L16tWoWrUq3NzcYG1tbXD/iRMniiwcERERUWEUuOgMGzYMkZGRGDp0KCcjExERkUkrcNHZsWMHdu/ejRdeeKE48hAREREVmQLP0XFxceGlEoiIiMgsFLjozJo1Cx9++CGuXbtWDHGIiIiIik6Bv7oaOnQoMjIyUK9ePZQtW9ZoMnJycnKRhSMiIiIqjAIXnblz5xZDjJIXEhKCkJAQaLVapaMQERFRMdGIiCgdQkn5vcw7EZm/9AvpOO1/GpblLGFhZxrXNNY91EF7X4sWu1vAvrG90nGIzEZ+P78LfETnUQ8fPkRWVpbBGMsCERERmYoC/0qTnp6OsWPHomrVqrC3t0fFihUN/hARERGZigIXnQ8//BB79+7FwoULYWtri59++gnTpk1DjRo1sGLFiuLISERERPRcCvzV1bZt27BixQr4+voiICAAnTp1Qv369eHq6opVq1ZhyJAhxZGTiIiIqMAKfEQnOTlZf/Xy8uXL608nf+GFFxAeHl606YiIiIgKocBFp27duoiJiQEANG7cGOvWrQPw75GeChUqFGk4IiIiosIocNEJCAjAqVOnAABTpkxBSEgI7OzsMHHiREyePLnIAxIRERE9rwLP0Zk4caL+v/38/HD+/HmcOHEC9evXR4sWLYo0HBEREVFhFGodHQBwc3ODm5tbEUQhIiIiKlr5/uoqIiIC27dvNxhbsWIF6tSpg6pVq2LUqFHIzMws8oDFJSQkBO7u7vDy8lI6ChERERWTfBedL774An///bf+9pkzZzBixAj4+flhypQp2LZtG2bMmFEsIYtDYGAgzp07h2PHjikdhYiIiIpJvotOVFQUXnzxRf3tNWvWwNvbG0uWLEFQUBC+//57/RlYRERERKYg30Xn7t27cHZ21t/+66+/0KNHD/1tLy8vxMXFFW06IiIiokLId9FxdnbWr5+TlZWFEydOoF27dvr779+/D2tr66JPSERERPSc8l10evbsiSlTpuDAgQOYOnUqypYti06dOunvP336NOrVq1csIYmIiIieR75PL//yyy/Rv39/+Pj4wMHBAT///DNsbGz09y9btgzdu3cvlpBEREREzyPfRcfJyQnh4eFISUmBg4MDLC0tDe5fv349HBwcijwgERER0fMq8IKBjo6OeY5XqlSp0GGIiIiIilKBr3VFREREZC5YdIiIiEi1WHSIiIhItVh0iIiISLVYdIiIiEi1WHSIiIhItUpt0QkJCYG7uzu8vLyUjkJERETFpNQWncDAQJw7dw7Hjh1TOgoREREVk1JbdIiIiEj9WHSIiIhItVh0iIiISLVYdIiIiEi1WHSIiIhItVh0iIiISLVYdIiIiEi1WHSIiIhItVh0iIiISLVYdIiIiEi1WHSIiIhItVh0iIiISLVYdIiIiEi1WHSIiIhItVh0iIiISLVYdIiIiEi1WHSIiIhItVRRdPr164eKFStiwIABSkchIiIiE6KKojNhwgSsWLFC6RhERERkYqyUDlAUfH19sX//fqVjqEafHw4qHcHItnEvKB2BiIjMkOJHdMLDw9GnTx/UqFEDGo0GmzdvNnpMSEgI3NzcYGdnB29vbxw9erTkgxIREZHZUbzopKenw8PDAyEhIXnev3btWgQFBSE4OBgnTpyAh4cH/P39kZiYWMJJiYiIyNwo/tVVjx490KNHjyfeP3v2bIwcORIBAQEAgEWLFmHHjh1YtmwZpkyZUuDXy8zMRGZmpv52ampqwUMTERGRWVD8iM7TZGVlITIyEn5+fvoxCwsL+Pn5ISIi4rmec8aMGXB0dNT/cXFxKaq4REREZGJMuugkJSVBq9XC2dnZYNzZ2RkJCQn6235+fnjttdewc+dO1KpV66klaOrUqUhJSdH/iYuLK7b8REREpCzFv7oqCnv27Mn3Y21tbWFra1uMaYiIiMhUmPQRHScnJ1haWuLWrVsG47du3UK1atUUSkVERETmwqSLjo2NDVq3bo2wsDD9mE6nQ1hYGNq3b1+o5w4JCYG7uzu8vLwKG5OIiIhMlOJfXaWlpSE6Olp/OyYmBlFRUahUqRJq166NoKAgDBs2DG3atEHbtm0xd+5cpKen68/Cel6BgYEIDAxEamoqHB0dC/s2iIiIyAQpXnSOHz+OLl266G8HBQUBAIYNG4bQ0FAMGjQIt2/fxmeffYaEhAR4enpi165dRhOUiYiIiB6neNHx9fWFiDz1MWPHjsXYsWNLKBERERGphUnP0SlOnKNDRESkfqW26AQGBuLcuXM4duyY0lGIiIiomJTaokNERETqx6JDREREqsWiQ0RERKpVaosOJyMTERGpX6ktOpyMTEREpH6ltugQERGR+rHoEBERkWqx6BAREZFqsegQERGRail+rSulhISEICQkBFqtVukoRETP1OeHg0pHMLJt3AtKRyB6plJ7RIdnXREREalfqS06REREpH4sOkRERKRaLDpERESkWiw6REREpFqltujwWldERETqV2qLDs+6IiIiUr9SW3SIiIhI/Vh0iIiISLVYdIiIiEi1WHSIiIhItVh0iIiISLVYdIiIiEi1ePXyYrx6Oa82TEREpKxSe0SH6+gQERGpX6ktOkRERKR+LDpERESkWiw6REREpFosOkRERKRaLDpERESkWiw6REREpFosOkRERKRaLDpERESkWiw6REREpFqltuiEhITA3d0dXl5eSkchIiKiYlJqiw4vAUFERKR+pbboEBERkfqx6BAREZFqsegQERGRarHoEBERkWqx6BAREZFqsegQERGRarHoEBERkWqx6BAREZFqsegQERGRarHoEBERkWqx6BAREZFqsegQERGRarHoEBERkWpZKR1AKSEhIQgJCYFWq1U6ChGRavX54aDSEYxsG/eC0hGoBJXaIzqBgYE4d+4cjh07pnQUIiIiKialtugQERGR+rHoEBERkWqx6BAREZFqsegQERGRarHoEBERkWqx6BAREZFqsegQERGRarHoEBERkWqx6BAREZFqsegQERGRarHoEBERkWqx6BAREZFqsegQERGRarHoEBERkWqx6BAREZFqsegQERGRarHoEBERkWqpouhs374djRo1QoMGDfDTTz8pHYeIiIhMhJXSAQorJycHQUFB2LdvHxwdHdG6dWv069cPlStXVjoaERERKczsj+gcPXoUTZs2Rc2aNeHg4IAePXrgjz/+UDoWERERmQDFi054eDj69OmDGjVqQKPRYPPmzUaPCQkJgZubG+zs7ODt7Y2jR4/q77tx4wZq1qypv12zZk3Ex8eXRHQiIiIycYoXnfT0dHh4eCAkJCTP+9euXYugoCAEBwfjxIkT8PDwgL+/PxITE0s4KREREZkbxYtOjx49MH36dPTr1y/P+2fPno2RI0ciICAA7u7uWLRoEcqWLYtly5YBAGrUqGFwBCc+Ph41atR44utlZmYiNTXV4A8RERGpk0lPRs7KykJkZCSmTp2qH7OwsICfnx8iIiIAAG3btsXZs2cRHx8PR0dH/P777/j000+f+JwzZszAtGnTij07lbw+PxxUOoKRbeNeeOZjmLvo5Cc3UX6Y6/5trrmLk+JHdJ4mKSkJWq0Wzs7OBuPOzs5ISEgAAFhZWWHWrFno0qULPD098cEHHzz1jKupU6ciJSVF/ycuLq5Y3wMREREpx6SP6OTXyy+/jJdffjlfj7W1tYWtrW0xJyIiIiJTYNJHdJycnGBpaYlbt24ZjN+6dQvVqlVTKBURERGZC5MuOjY2NmjdujXCwsL0YzqdDmFhYWjfvn2hnjskJATu7u7w8vIqbEwiIiIyUYp/dZWWlobo6Gj97ZiYGERFRaFSpUqoXbs2goKCMGzYMLRp0wZt27bF3LlzkZ6ejoCAgEK9bmBgIAIDA5GamgpHR8fCvg0iIiIyQYoXnePHj6NLly7620FBQQCAYcOGITQ0FIMGDcLt27fx2WefISEhAZ6enti1a5fRBGUiIiKixyledHx9fSEiT33M2LFjMXbs2BJKRERERGph0nN0ihPn6BAREalfqS06gYGBOHfuHI4dO6Z0FCIiIiompbboEBERkfqx6BAREZFqsegQERGRapXaosPJyEREROpXaosOJyMTERGpX6ktOkRERKR+ii8YqLTcxQpTU1OL/LmzH6QX+XMWVn7eJ3MXHeYuWc/KnZ6WjnRdOiy0FrDQmsbveTqtDjqdDqlpqdCmap/4OHPc3gBzFyU15y7M8z5r0WGNPOsRKnf9+nW4uLgoHYOIiIieQ1xcHGrVqvXE+0t90dHpdLhx4wbKlSsHjUajdJw8paamwsXFBXFxcShfvrzScfKNuUsWc5cs5i5ZzF2yzCG3iOD+/fuoUaMGLCyefIS21H91ZWFh8dQmaErKly9vsjvc0zB3yWLuksXcJYu5S5ap53Z0dHzmY0zjS2oiIiKiYsCiQ0RERKrFomMGbG1tERwcDFtbW6WjFAhzlyzmLlnMXbKYu2SZa+68lPrJyERERKRePKJDREREqsWiQ0RERKrFokNERESqxaJDREREqsWiQ0RERKrFokOUTyICnU6ndIwCM9fc5owns5Ysc93e5prb3PD0cjNz6dIlxMXFITY2Fk5OTujcuTMqVKigdCzV02q1sLS01N/W6XTQaDQme320XOaYW0RMOl9BcHuXLHPY3nkx9dzm/rnDomNG5syZg+XLl+P8+fNo1KgR7OzsoNVq0alTJ7z55pvw8vJS1Q8tUxAfH4/ffvsNly5dQnx8PPz8/DBo0CA4OTkpHe2pzDX3o8xpX9ZqtQgLC0NycjISExPRqFEjdOnSBTY2NkpHyzdu7+JnjrnV8LnDomMm7t27BxcXF3zzzTd45513cPnyZZw7dw6RkZGIjIyEVqvF119/jfbt2ysd1YhOp0NGRgYcHByUjlIg9+/fxyuvvIJLly6hTZs2sLGxwaFDh5CQkIDu3bvj448/RqdOnZSOacRcc2dlZeHgwYNo0qQJqlevbnCfKf8gTU9Px5gxY7BlyxbY2tqibt26yMjIQJkyZdCzZ08MHToUdevWNbn3wO1dsswxtzl/7hgQMgs//vijeHp6Go1nZWXJoUOHpEePHlKxYkW5cuWKAumebtGiReLu7i4zZ86UixcvilarNXpMWlqaRERESE5OjgIJ8/bNN99Iy5Yt5fbt2yIicvfuXYmJiZFff/1V/P39xcPDQ8LCwhROacxccy9YsEAqVqwob731lixatEiOHj0qd+/eNXhMUlKSLF26VDIzM5UJmYdvv/1W3N3d5ciRIyIicvz4cVmxYoWMGTNG2rdvL6+//rrR+zAF3N4lyxxzm/PnzqNYdMzEpk2bpGHDhhIREZHn/RkZGeLt7S0//vhjCSd7ttatW0uDBg3ExcVFrK2tpWvXrrJ8+XKJj4/XP2bZsmXy4osvKpjS2CuvvCLjx483GtdqtRITEyOvvPKKNG7cWJKSkhRI92Tmmrtz587SuXNn6d69u9SsWVM8PT1lzJgx8uuvv8q5c+ckMzNTli5dKq6urkpHNfDCCy/Id999ZzSempoqW7dulVq1akmvXr3yLPhK4vYuWeaY25w/dx7Fs67MRNeuXVGjRg189913OHbsGHJycgzuL1OmDKytrREfH69QwrwlJibCwsIC06dPx5UrV7Bjxw5UrVoV48ePR5MmTfDGG29gx44d+OGHH9CwYUOl4xro3bs3tm7ditjYWINxCwsLuLm5Yd68ebCzs8PJkycVSpg3c8x9584dAEBAQAB2796NQ4cO4bXXXsPRo0cxefJkjB49GtOmTcOXX36Jfv36KZz2/+Xk5KBNmzbYsWMH7t69a3BfuXLl0KdPH4SGhuL69es4d+6cQimNcXuXLHPNba6fO0aUblqUfwcPHpRmzZqJvb29vPnmm7Jt2zb5+++/JTIyUpYsWSLly5eX6OhopWMaiI2NlWnTpsnOnTsNxpOTk2XFihXi6+sr1tbWotFo5Nq1awqlzNv169elU6dO0qFDB1m+fLlcvnxZsrKy9PdfuXJF7Ozs5OrVqwqmNGaOue/fvy8bNmyQP/74w+i+Y8eOyfjx46Vu3bomuZ9ERERIo0aNZMqUKXkewr927ZrY29tLbGysAunyxu1d8sw1d3h4uDRt2lTs7e1l6NChZvG58zhORjZDq1atwsKFCxEREYFq1aqhXLlyyMzMxMSJEzF+/Hil4xm5c+cOrK2tUb58+Twn2k2YMAEHDhzAiRMnFEr4ZMePH8eMGTNw5swZ1K1bF23btoWTkxMyMzMRFhaGjIwMhIeHKx3TyMmTJ/HVV1/pc3t5eZl87uzsbGg0GlhZWUGr1QKAwanx06ZNw6ZNmxAVFaVQwrxptVosWbIEn3zyCaytrTFw4EC8/PLLqFKlCq5du4bt27fj+PHjJnUEDQAyMzMBALa2ttBqtRARWFlZ6e/n9i5a5po7V2hoKBYvXoz//e9/qF69usl/7hhQtmdRfmm1WsnOzjYYS0lJkd9++0127NghN27cUChZwWRnZxt8B/3gwQOpU6eOfPXVVwqmerbdu3fLm2++KV5eXuLt7S2urq7y0UcfmfwkvB07dsiQIUOkbdu2ZpX7UVqtVtLS0sTNzU2mTZumdJwnyszMlG+++Ubc3d1Fo9FIw4YNxdnZWfr27aufgGoOuL2Ll7nkfvjwoZw/f17++usvg/F79+7J2rVrZfv27WbzucMjOiYuPT0d9vb2BmO5q9xaWJj+FKuHDx8iNjYWOp0OTZo00Y/LvxPhkZaWhp9++gljxoyBnZ2dgkn/X3Z2No4fP459+/ahdu3a8PDwQPPmzQH8++9x7do1NGnSxGQX+Mrrt/P09HRcv34dDRo0MNncOp3uqfv08ePH0aJFC5NacyQxMRFnz56FnZ0dOnTooB9PSEjA4cOHUa9ePTRp0sSkMosILl++rP/NvG3btihfvrz+fq1WCwsLC0RGRnJ7FxFzy/3HH39gxowZiI2NhYjgzp076NKlCwIDA9G9e3el4xWcgiWL8mH8+PGybNkyOXXqlNy/f9/gPp1OJzk5OZKamqpQuqdbvny5uLi4iLu7uzRt2lSaNWsmH330kVy6dMngcY/OHzEFw4cPFycnJ2nVqpVUqlRJLCwspGXLlvL9998rHe2p/vnnH4PbWq1WMjMzjU7Z1+l0JRnrmRITE43GTC1jXr799lupXLmyNG3aVJycnKRChQoybNgwiYyMVDraU02aNEkqVqwozZo1EwcHB7Gzs5NevXoZzaMzNea6vc0xd/Xq1WXChAmyZs0a2bVrlyxatEj8/PykTJky0r59ezl8+LDSEQuERceErVu3TjQajbi6ukq7du1k0qRJsnHjRomOjtava/HgwQPp2rWrnDhxQuG0hlavXi2urq4SHBws69atk6VLl0pQUJC0bNlSGjRoIMHBwSa1Nkeun3/+WerVqydhYWFy584dycrKkqNHj0pAQIDY29tLgwYNJDw8XOmYRrZt2yb169eXwMBAWb9+vX4NnVw5OTmSnp5ucut0/PHHH+Ln5yczZ86U8PBwSUlJMbhfp9PJgwcPjMaVtnLlSqlTp44sWLBA/vrrL9m7d6/MnDlTvLy8pGzZsvLmm2/KzZs3lY5pZMWKFVKvXj1Zv369XLhwQa5cuSLr16+XXr16iZWVlfj4+Mj58+eVjmnEXLe3OeZet26duLq6GvwCqtPp5N69e7J7927p2bOndOvWTe7cuaNgyoJh0TFhI0eOlICAADl06JBMmTJFmjZtKi4uLtKtWzf56quvZO/evfLjjz+Kra2t0lGN+Pj4yEcffWQwdv/+fTl58qR8+umn4urqKrNnz1Yo3ZO9+uqrMmbMGP3tR+cTxcbGSu/evaVPnz5KRHuq3r17S506daRnz57i7e0t3bt3lylTpsiePXv0hXLz5s1iYWGhcFJD/v7+UqFCBWnbtq288MILMmLECFm0aJGcPHlSPydtx44d4ubmpnBSQ35+fvLBBx8YjGm1WklISJDly5dL8+bNjfZ/U9CzZ095//3387zvr7/+ko4dO8rbb79dwqmezVy3tznmXrVqlbRt2/aJRebIkSPi4uIia9asKeFkz8/0J3mUUlqtFq6urqhYsSI6dOiAGTNm4OzZs1i+fDlq1qypn9cyceJEDBo0SOm4BnJyclClShVYW1sbjDs4OMDT0xNffPEFhg4ditWrVyMxMVGhlHlr3bq1wVkmFhYWyM7ORnZ2NmrXro3x48fjwoULCAsLUy7kYzIyMpCUlIRJkyZhzpw5ePfdd+Hq6oojR47gk08+waBBgzBz5kzMnj0b/fv3VzquXnp6OpKSkjB79mwsWrQI/v7+iI2NxeLFizFp0iRMnjwZa9aswaxZs+Dl5aV0XD2dToc6deogNTXVYNzCwgLOzs4YPnw4RowYgR07diA6OlqhlMZEBE2bNkVMTIzBuE6ng4igc+fOCAwMxKFDh3D8+HGFUhoz1+1trrl9fHxw9epVjBw5En///bd+Tmgub29veHp6IjIyUqGEz0HppkVPlpCQoJ/P8vg8lgcPHsj8+fNFo9GY5He9CxYsEGtra/n555/z/M3g6tWrUrVqVTl37pwC6Z7s5MmTUq5cOenZs6ccPHjQ6P6HDx9KpUqV5Pjx4wqky9udO3fks88+k0WLFunHdDqdREZGypw5c2Tw4MHi7e0tGo1Gjh49qmBSQzdv3pRJkybJsmXL9GNarVb27Nkj77//vnTs2FE8PDxEo9HI//73PwWTGlu/fr1oNBr54osvJCYmxuj+xMREcXJyktOnT5d8uKfYu3evaDQaGTlypJw6dcro/tTUVKlcubKcPHmy5MM9hblub3PNfeDAAWnXrp306tVLZs2aJfv379ev73PgwAFxdHSUQ4cOKZwy/3jWlZnRarXQaDSwsLDAsmXLMH78eKSlpSkdy0hOTg6mTJmC33//Hb6+vnjllVdQv359ODs7w9raGgsXLsS3336LGzduKB3VyJEjR/DZZ58hKSkJDRo0QPv27dG9e3dotVrMmjULhw4dwuXLl5WOaSQ7OxvW1tbIyckxOOMqMzMTH374ITZt2oR//vlHwYTGMjIyICKwt7fX58+Vnp6O8ePHY+/evUZHIUzB999/j9DQUDRq1Ag+Pj5o1qwZmjRpAltbW8yaNQs//fQT4uLilI5pZOPGjfjvf/8LBwcHtGzZEq1bt0b79u1hY2ODGTNmYPPmzbh27ZrSMY18//33WL58ORo1agRfX1+z2d7mlju3Euzfvx8//vgjDh8+jEqVKqFixYq4cuUKLCws0K1bN/z4448KJ80/Fh0TlXuK55NOAxYRzJw5E+np6Zg2bVoJp3u63NOEU1JSEBoaivnz5yMmJgatWrVCrVq1cPjwYdSqVQujR4/GyJEjlY5rIDf72bNnsX37dhw9ehQ3b97E2bNnkZmZid69e2P06NHw9/dXOqqePOFqx4+ert2hQwd4eXlh3rx5JR2vwB7d95s3bw4fHx/Mnz9f6VhGHj58iC1btmDp0qW4ePEiqlevDq1Wi7///hutWrXCe++9hyFDhigd04CIQKvVIjw8HGvXrsWpU6eg0Whw/fp1xMfHw8/PD++9955JXfYhdz9OT0/Hjh07sGzZMpw7dw7Vq1eHTqcz2e1tjrl1Op3R8hM3btzAjh07cO3aNbi4uMDNzQ1+fn4Gv0yZOhYdE5f7z/O0wmNqa6KICFJTU+Ho6Kgfi4qKwrp165CSkoImTZqgc+fOaNasmUmtBfT4kRDg3/UvYmJiYG1tDWtra9SrVw9ly5ZVKGHetFqtwQrCj8vMzMQ333yDESNGoFatWiWY7Omelfvhw4eYNGkSJk2aBDc3t5IL9gx5fRhcvHgRYWFhePjwIVxdXdGmTRu4uroqmNJYXtv7n3/+wYkTJ6DT6eDk5ISmTZuicuXKCiXMW0ZGBtLT01GlShX9WHR0NP78809kZGTAzc3NJLe3ueYG/t1XtFotrKysTOpn9PNi0TFBK1asQOPGjdGiRQuDRfTy+gFravbt24fly5fjwoULSE1NhZ+fHwYNGoROnTopHS3fcnJyoNPpTGbxrvwyp4UkH2WuuXP3E2tra5P+f/Jx5pR7w4YNCA0NxcmTJ6HT6dChQwf069cPffv2hYODg9Lxnsgcc3/zzTdo0aIFfHx8DBapzc7OBgCjk0vMiXn9ZCkFDh48iICAAHzyySeYNGkSfv75Z1y4cAEA9IfzMzMz8cUXX5jc/JZDhw4hMDAQsbGxePXVV9G3b18cPHgQvr6+8PT0xKZNmwD8/1EqU3HkyBF06tQJv/zyC7KysmBlZaUvOVlZWfqVhpOSkkwqe1RUFIYMGYLt27cjJycHFhYW+rKQm9kUnT9/HlOnTsXBgwf1X1OZY+7c/USj0SA7O9voys6m4vH95NHcWVlZJps7PDwcH374IcqUKYM5c+bgP//5D5KTkzFs2DA0b94cP/30k9IR82SOuQ8ePIiPP/4YX375JV5//XV8+eWXOHLkCADoj2Y/ePAA48ePN7l5fvlSolOf6ZkmTJggXl5eEhQUJD4+PtKyZUvp0aOHfPTRR7Jp0ya5fv26REREiEajMVopWWn9+/eXESNGGIxptVo5duyYDBkyROrVqycbN25UKN2TvfXWW2JtbS2urq5SqVIlGTBggOzevdvgMQcPHhR/f3+j640p6a233hI7Oztp0aKFtG7dWiZOnGi0YumhQ4dk0KBBRqsjK+mtt94Se3t76dSpk7z66qsya9YsOXPmjMFjDh8+LKNGjTKpFZLNObc57icDBgyQkSNHGo3fvn1bJk2aJFWqVJE5c+aUfLBnMMfckydPlk6dOsmsWbNk2LBh0qlTJ2nfvr0MHDhQfvjhB7l48aIcOXJENBqNya7E/zTmM5uolLhz5w46duyIWbNmITMzE7t27cKOHTuwf/9+7N27F3Xr1sXZs2fRtWtXkzsEmpSUBA8PD/3t3Ml4bdq0QUhICEaNGoUZM2agc+fOJjUPICYmBp9++in8/Pxw7Ngx7NixA4MHD0bZsmXRv39/jBo1CmvWrMHNmzdNagLehQsXMGnSJDRr1gzHjx/H8ePHsXPnTjg7O6NHjx547bXXsGrVKpw9e/apc2FK2unTp/Hee++hSpUqiIyMxG+//YZNmzahXr166NKlC7p164YVK1bgwIEDJvXVirnmNtf9JDMzE7a2tvrbWVlZsLCwgJOTE2bOnAkRwdKlS/Hqq6/CxcVFwaSGzDF3UlISGjVqhKCgIOTk5ODIkSMIDw/HiRMn8Ouvv2Ljxo24fPky/P39Ua5cOaXjFpzSTYsMnTlzJs9rzty4cUOWLl0qr776qmg0GtmxY4cC6Z5uzpw5UrNmTYmOjjYYz/3t9urVq9KgQQOJiopSIl6e4uPjZcSIEbJkyRIR+Xe9ovj4eAkLC5PPPvtM2rZtKxUrVhSNRiNbt25VOO3/u3r1qvTq1UsWL14sIiLp6ely8uRJ+emnn2TkyJHi7e0tjRs3Fo1GI1u2bFE47f+7dOmS+Pj4yPLly0VE9MvKf/LJJ9K7d2/x9vYWHx8f0Wg0snnzZmXDPsJcc5vrfiLy7+UTqlSpYrT2U+5q5bdv35Y6deqY1BW/Rcwz982bN2Xfvn1G43fu3JHt27fLxIkTTfZzJz9YdExQ7uKAWq1WsrOzDS5DsG3bNnF0dFQo2dMlJiZKt27dpFGjRhIcHCwHDhwwOMy5ceNGcXBwUDBh3pKSkowuiCny7+KA165dk0mTJpnkNk9ISJArV64YjScnJ8vhw4flzTffNMncly9flosXLxqNx8fHy4YNG6R79+5SoUIFBZI9nbnmNsf9RKfTyf3792XQoEFSqVIleeutt2Tjxo0G12pbs2aNyf08yc39+uuvm1XuR+l0OtFqtQZfv27dulXs7e0VTFU4POvKTMi/pRQDBgxAamoq9uzZo3SkPF26dAkLFy7EwYMHYWNjAxcXF5QtWxbp6ek4d+4cXnrpJcycOVPpmE8lj52y37dvX1hZWWHDhg0Kpnq63Em8j3611rdvX9jZ2WHNmjUKJnu6J+W2t7fHqlWrFEz2dOaaG8h7/zbV/SQtLQ2hoaHYunUrkpKSYGNjg3LlykFEEB8fj4EDB5rcOmIAcP/+fSxfvhy///47kpKSYGlpaRa586LT6TBmzBgkJydj3bp1Ssd5Liw6JiI7Oxvnzp3Djh07UL58ebRs2RJubm5wdnaGlZWVfg2MnJwcpKSkmNQcl7ycOXMG27dvx4ULF3D37l1kZGTg/fffR9euXU1uHZqnSUtLw7hx4zBhwgR4enoqHcfI46sJA//+YEpNTUX//v3x7bffmtR1op5GRHDnzh14e3tjxYoV6Nixo9KR8sVccwPAvXv3zGI/uXHjBsLDw3H+/HnExcUhMzMTgYGBaN26tcF8GFNz8eJFHD58GNeuXcP169fx8OFDk8z96Ir7T7r//v37qFChQskGKyIsOiZi0qRJWL16NapWrYrk5GTExcWhbt26GDx4MCZMmAAnJyelIz7RjRs3sHr1akRERKB+/frw9PSEt7c36tSpA61Wi4yMDJOewPas9Yken1yotOjoaMyfPx8nTpxAw4YNUbduXTRv3hxeXl6oVq2a/nEPHjxAmTJlFExqKHc/+d///ocGDRqgadOmaNKkCRo0aAAHBweDlWQfXcdDaeaaO9ezPsQyMjJM6peP3P07MjIS9evXR8OGDdGuXTt4e3ubVM4nyT1d/9GjfY+uUm5KTp8+jRYtWhiMPWtVfnPEomMCzp07h3bt2mHNmjVo2bIlnJ2dERcXh2XLlmHp0qVIS0vD/PnzMXToUJNbCfnatWsYOHAgkpOT0apVK5w+fRq3bt1CjRo10KNHD3z88ceoVKmS0jGNJCYm4siRI+jVq5fBmSaP/0B6+PAh7OzsTGa7X716Fb1790a5cuXQrl07/P3330hMTIS1tTVatGiBMWPGGPxmbiq589pPEhMTUatWLfTo0QOTJ082KPPMXTh37tzBxYsX0aFDB/2YiOj379yMuUeKTSV3Xvv3rVu3YGlpCQ8PD4wbNw5t2rRROqaRvLa3TqfTr29lZWWl3/6mcmZbdHQ0mjRpAm9vb7z44ot47bXX0KxZM4PHZGdn4+TJk/D09DS7BVQNlNBcIHqK6dOnS+fOnfW3H13LIi0tTSZMmCDNmzeXxMREJeI91ejRo6VXr14SFxenH4uJiZHg4GCpUqWKVKtWTXbt2qVgwrwFBgaKRqMRJycnGTZsmNGVeHU6ncTExMjMmTPl4cOHCqU09u6770qfPn0MJjfGx8fLggULpGXLllK+fHlZuXKlcgGf4Fn7SfXq1Y3WLjIF5pp7/PjxotFopGHDhjJ58mS5cOGCwf1arVZiY2Nl3bp1JrU21LP273Llypnk/v2s7Z2Tk6Pf3qayVtHnn38utWvXlnfffVc6dOggjRo1En9/f5k/f77Ex8eLiEhcXJxoNBqD/d8cseiYgN9++00aN24ssbGx+rHs7GzJzMwUkX9PbW3RooUsXLhQqYhP1KFDB5k1a5aI/Hu22KM/NLVarbzyyivSt29fERGTWkTN29tbgoKC5Pvvv5eOHTuKpaWl1K5dW6ZMmaI/Q+Xzzz+XevXqKZzUkL+/v0ybNk1E/v3h+fgPzXfffVc6deokDx8+NKntba77ibnmbtWqlbz99tsyZcoUfUFo1aqVzJ49W5KTk0VEZNq0aVKnTh2Fkxoy1/3bHLf34MGD5f3335f4+Hg5evSo/PDDDzJ06FBp2bKluLu7y5AhQ6Rv377i7u6udNRCY9ExAUlJSdK4cWNxd3eXDRs25HkEoUWLFvq1MEzJZ599Jm3atDHInJWVJRkZGSIiEhYWJvXr1zdaU0JJ169flwEDBujXzklNTZXjx4/Lf/7zH/2aIq1atZJy5cqZ3Aqms2fPljp16hicLpyZmanf/lFRUVKnTh3566+/lIqYJ3PcT0TMM/e1a9fE399fVqxYIZmZmXLp0iVZv369jBw5UurVqyfly5cXf39/qVy5ssyePVvpuAbMcf82x+2dnZ0tK1eulBkzZhiMJyQkyJ9//ilff/219O3bVzQajf7npDlj0TER8fHxMmjQIGnRooX07NlTgoODZf/+/RITEyNBQUFSuXJlSUtLUzqmkWPHjkm1atWkTZs2sm3bNqP7L168KLa2tpKenq5Aurylp6fL1q1bjRbs0mq1kpSUJGFhYdK7d2+xtLTUf6CZiitXroinp6fUrVtXQkNDje4/e/asWFtbm9T2FjHP/UTEPHOnpKRIaGio7N+/32D83r17EhUVJUuXLpVOnTpx/y4i5ry9c+Wu3fao1atXi0ajMalt/bw4GdmEJCYmYufOndizZw9iY2Nx8eJFJCUlwdfXF++88w7eeOMNpSPmKTo6Gh999BGOHz+OypUro2PHjujZsycuXryI1atXw8XFxaTXoJE8JmK++eabiI2NRXh4uEKpnuz+/fuYMmUK1qxZg5ycHHTr1g0vvfQSzp49i/3796N58+b45ZdflI5pxFz3E3PNDfy7b+degPRRgwYNQmJiIvbt26dQsicz1/0bMJ/t/aSzwB69OPDkyZNx7Ngx7N+/v+QDFjEWHYXdunULMTExsLW1RZkyZVC3bl1YWFjgypUryMjIgL29PZycnEzyzKVHpaenIywsDHv37sWxY8dw5swZVK5cGSNGjMDQoUPh5uamdES9Z51O/uDBA7zyyit477330K9fvxJO93S5P6AePnyIM2fOIDw8HHv37kVkZCTq1q2LIUOGoH///qhevbrSUfNkTvvJo8w1dy555IyrBw8ewMfHB1OmTMGrr76qdDQD5r5/5zL17Z37uWNjYwMRgZubm8HabCKCLVu2oGbNmia9vlJ+segoaMmSJVi+fDlOnDgBKysrNGrUCE2aNMGLL76Il19+2eQXBdy5cyfu3r0LrVYLFxcXtG3bFvb29sjIyIClpSXu379v0uv/PEl2djaOHz+O9u3bKx3lmR5d8yIlJQWOjo5KRzJirvuJuefOyclBlSpV4O3tbfCzJDMzE3v27EGvXr0UTJk/5rR/m8v2fvxzx93dHY0bN0bHjh3Rq1cv1KpVS+mIRY5FRyF37txBgwYNEBgYiJEjRyI1NRU7d+5EWFgYLl++jGbNmmHevHmoU6eOyaxxkev+/ft499138eeffyInJwfVq1eHvb09KleujO7du+O1117T/89iSgtl5eTkIDk5GVWrVlU6SoGYa25z3U/UkrtGjRpwcHBA5cqV4evri4EDB8LV1VXpmEbUsn+bw/Z+2udOdHQ0mjdvjjlz5qBOnTrIyckx+grObJX8tCASEZk3b554e3vned/evXvFy8tL3N3dDdaTMBXTp0+X5s2bS3h4uIj8e8X1RYsWyZAhQ6RFixby2muvyb179xROaWzOnDlSoUIFGTt2rISHh+c5yS4lJUW2b9+uP7XfFOQ3986dO/OcVKgUc91P1Jjbw8NDBg4caJK51bh/m+r2NufPncJg0VHIggULpGnTpnL+/HkREXnw4IHBh+v58+elYcOGsm7dOqUiPlHHjh1l7ty5RuNarVZ2794ttWvX1q8tYkratm0rHTp0EC8vL7GwsJDGjRtLcHCwnDlzRr9ex4IFC574g0Ap5prbXPcT5i5Z3L9Ljjl/7hSGaRxzLYVee+01WFhY4IcfftBfZsDGxgY6nQ4A0LhxY1SuXBmxsbEKJzWUnZ2Npk2bYtOmTbhz5w6Afw89536X3r17d4SEhCA6Ohpnz55VOO3/u337NmxsbPDee+/h6NGjOHv2LPr164fQ0FB4enrCx8cHixYtwoIFC+Dt7a10XD1zzW2u+wlzlyzu3yXLXD93Ck3pplUaabVa0el08ttvv0mtWrWkfPnyMnLkSDlx4oSIiNy4cUN+/fVXcXBwkJiYGGXD5iEiIkLq168v//nPfyQpKcno/ri4OLG3t5fr168rkC5vN27ckNmzZxst15+TkyPh4eEyfPhwcXR0NLnlzs01t4h57icizF2SuH+XHHP/3CkMFh0FPXz4UP7++29ZsGCB+Pv7i729vTg4OEijRo2kbt268umnnyod0YhOp5OsrCxZvHixVK5cWSpUqCCjRo2Sffv2ydWrV2Xjxo0yfPhwad26tdJRjWRkZOgX7Mpr+fgPPvhAWrZsWdKxnskcc5vrfsLcJY/7d8kyx8+dwuJZVyUsKSkJa9euxcyZM1G5cmVUqlQJFStWRNu2bdGyZUtkZGTg6tWr6NGjBxo0aGBSZ1s97t69ewgNDcWvv/6KqKgoODo6ws7ODq1atcLUqVPRrl07pSPm28OHD+Hp6YmAgAB89NFHSsfJN3PIba77CXMrj/t30VDT587zYNEpYW+//TZOnTqFHj16wMHBAXfu3EF0dDTi4+Ph6uqKadOmwd3dXemYeXrw4AHKlCljMCYiePDgAdLS0nDmzBk4ODiY1HfpQN6583rMunXrMHjwYNjY2JRQsqdTU25z3U+Yu/hw/y455vy5UxRYdEqQiMDBwQE7d+6Ej4+Pfiw6OhoHDhzATz/9hOTkZGzYsAHNmjVTOK2xDz74AB07dkTr1q1RrVo12NraGj3m7t27qFixokmt/ZOf3Pfu3UOFChVKPtxTqDm3ue4nzF10uH+XDHP/3CkSJftNWel29uxZadasmRw7dizP+zMyMqRFixYSHBxcssHyYdWqVaLRaMTa2lrq1KkjEydOlL1790pCQoJ+bYuUlBR55ZVX5PTp0wqn/X9Pyn3r1i3Jzs4WEZG0tDTp06ePnDlzRuG0/09tuc11P2Hu4sH9u+SY8+dOUWHRKUEZGRnStWtX6dy5s1y9ejXPiXezZs0yyQlsI0aMkPfee0+uXLki06dPFzc3N9FoNNKqVSuZMWOGnDhxQpYtWyZWVlZKRzXA3CWLuUsWc5csc8xtzp87RYVFp4QdPnxYPD09pWPHjrJy5Uq5ceOG/oyDhw8fymuvvSZvvPGGwikNZWdny1dffSVTp041GD916pSMGjVKHB0dxcHBQaytrSUgIEChlMaYu2Qxd8li7pJlrrlFzPNzpyix6Cjg9OnT8tprr4mdnZ04OTlJ37595d1335U6deqIl5eXnDp1SumIRu7evSsXLlwQEZHMzEyj3wpWrlwpGo1GoqKilIj3RMxdspi7ZDF3yTLX3CLm+blTVFRyxS7z0rx5c6xbtw6JiYnYvn07Nm/ejOTkZAQEBGDAgAFo0qSJ0hGNVKhQQT8pMPfsB51OBxGBpaUlMjIyYGdnBw8PDwVTGmPuksXcJYu5S5a55gbM83OnqLDoKKhq1ap4++238fbbb5vUVZDz69G89+/fx7Rp0xRMk3/MXbKYu2Qxd8kyt9zm/rnzPHh6ORWJ7OxsWFpamt3/NMxdspi7ZDF3yTLX3GrHokNERESqxdpJREREqsWiQ0RERKrFokNERESqxaJDREREqsWiQ0RERKrFokNERESqxaJDREREqsWiQ0RERKrFokNERESq9X9Eph74jAXEQgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_samples(data[2][1]) # plot the samples for (j,n) = 1,6, i.e. the configuration which contains both a physical and minimum-energy solution"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "cd31e406-cf82-4ef2-8df0-a17cf06a8d15",
   "metadata": {},
   "source": [
    "We can see that the best path is $b=(1,0,0,1,0,0)$, that is, the path $saf$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbb8e055-474a-4456-823f-9348e06c6176",
   "metadata": {},
   "source": [
    "## Bibliography\n",
    "[1] Thomas Krauss and Joey McCollum. “Solving the Network Shortest Path Problem on a Quantum Annealer”. In: IEEE Transactions on Quantum Engineering 1 (2020), pp. 1–12. doi: 10.1109/TQE.2020.3021921.\n",
    "\n",
    "[2] Kamil Bradler and Hugo Wallner. Certain properties and applications of\n",
    "shallow bosonic circuits. 2021. doi: 10.48550/ARXIV.2112.09766. url:\n",
    "https://arxiv.org/abs/2112.09766."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "3065aa56-30e8-4aee-8b6c-f8aff1615c01",
   "metadata": {},
   "source": [
    "## Appendix\n",
    "If you want to try more complicated graphs, this code generates the objective function for any weighted, directed graph:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "dcfdc370",
   "metadata": {},
   "outputs": [],
   "source": [
    "def alpha_graph(graph): \n",
    "    \"\"\"calculate alpha as sum of all the costs + 1\"\"\"\n",
    "    alpha = 0 \n",
    "    for i in range(0,len(graph)):\n",
    "        for j in range(0,len(graph)):\n",
    "            if graph[i][j] != 0 : alpha+=graph[i][j]\n",
    "    return alpha+1\n",
    "\n",
    "def graph_to_s(graph, s, weight):\n",
    "    \"\"\"calculate Hs\"\"\"\n",
    "    ham = np.zeros((weight, weight))\n",
    "    edges_coeff = [0]*weight\n",
    "    constant_coeff = -1\n",
    "    count = 0\n",
    "    for i in range(0, len(graph)):\n",
    "        for j in range(0, len(graph)):\n",
    "            if graph[i][j] != 0 :\n",
    "                if i == s:\n",
    "                    edges_coeff[count] = 1\n",
    "                if j == s:\n",
    "                    edges_coeff[count] = -1\n",
    "                count += 1\n",
    "    \n",
    "    for i in range(0, weight):\n",
    "        for j in range(i+1, weight):\n",
    "            ham[i][j] = edges_coeff[i]*edges_coeff[j]*2\n",
    "            \n",
    "    for i in range(0, weight):\n",
    "        ham[i][i] += edges_coeff[i]*edges_coeff[i] + constant_coeff*edges_coeff[i]*2            \n",
    "\n",
    "    return ham\n",
    "\n",
    "def graph_to_f(graph, t, weight):\n",
    "    \"\"\"calculate Hf\"\"\"\n",
    "    ham = np.zeros((weight, weight))\n",
    "    edges_coeff = [0]*weight\n",
    "    constant_coeff = +1\n",
    "    count = 0\n",
    "    for i in range(0,len(graph)):\n",
    "        for j in range(0,len(graph)):\n",
    "            if graph[i][j] != 0 :\n",
    "                if i == t:\n",
    "                    edges_coeff[count] = 1\n",
    "                if j == t:\n",
    "                    edges_coeff[count] = -1\n",
    "                count += 1\n",
    "    \n",
    "    for i in range(0, weight):\n",
    "        for j in range(i+1, weight):\n",
    "            ham[i][j] = edges_coeff[i]*edges_coeff[j]*2\n",
    "            \n",
    "    for i in range(0,weight):\n",
    "        ham[i][i] += edges_coeff[i]*edges_coeff[i] + constant_coeff*edges_coeff[i]*2              \n",
    "\n",
    "    return ham\n",
    "\n",
    "def graph_to_i(graph, our_i, weight):\n",
    "    \"\"\"calculate Hi\"\"\"\n",
    "    ham = np.zeros((weight, weight))\n",
    "    edges_coeff = [0]*weight\n",
    "    count = 0\n",
    "    for i in range(0, len(graph)):\n",
    "        for j in range(0, len(graph)):\n",
    "            if graph[i][j] != 0 :\n",
    "                if i==our_i:\n",
    "                    edges_coeff[count] = 1\n",
    "                if j==our_i:\n",
    "                    edges_coeff[count] = -1\n",
    "                count += 1\n",
    "    \n",
    "    for i in range(0, weight):\n",
    "        for j in range(i+1, weight):\n",
    "            ham[i][j] = edges_coeff[i]*edges_coeff[j]*2 \n",
    "            \n",
    "    for i in range(0, weight):\n",
    "        ham[i][i] += edges_coeff[i]*edges_coeff[i]\n",
    "\n",
    "    return ham\n",
    "\n",
    "def graph_to_c(graph, weight):\n",
    "    \"\"\"calculate Hc\"\"\"\n",
    "    ham = np.zeros((weight, weight))\n",
    "    count = 0\n",
    "    for i in range(0,len(graph)):\n",
    "        for j in range(0,len(graph)):\n",
    "            if graph[i][j] != 0 :\n",
    "                ham[count][count] = graph[i][j]\n",
    "                count += 1\n",
    "\n",
    "    return ham\n",
    "\n",
    "def graph_to_hamiltonian(graph,s,f): \n",
    "    \"\"\"returns the graph adjacency matrix as a Hamiltonian\"\"\"\n",
    "    weight = np.count_nonzero(graph)\n",
    "    ham =  np.zeros((weight, weight))\n",
    "    alpha = alpha_graph(graph)\n",
    "    ham += alpha*graph_to_s(graph, s, weight)\n",
    "    ham += alpha*graph_to_f(graph, f, weight)\n",
    "    for i in range(0,weight):\n",
    "        if i!=s and i!=f:\n",
    "            ham += alpha*graph_to_i(graph, i, weight)\n",
    "    ham += graph_to_c(graph, weight)\n",
    "\n",
    "    return ham"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "7bbfb18f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  2.  32. -32. -32.  32.   0.]\n",
      " [  0.   1.  32.   0. -32. -32.]\n",
      " [  0.   0.  35.  32. -64. -32.]\n",
      " [  0.   0.   0.   2. -32.  32.]\n",
      " [  0.   0.   0.   0.  35.  32.]\n",
      " [  0.   0.   0.   0.   0.   4.]]\n"
     ]
    }
   ],
   "source": [
    "\"\"\"example: our 5-edge graph\"\"\"\n",
    "# start,a,b,finish -> 0,1,2,3\n",
    "G1 = np.array([[0,2,1,0], [0,0,3,2], [0,3,0,4], [0,0,0,0]]) # weights associated to all possible paths\n",
    "G1_s = 0 # start index\n",
    "G1_f = 3 # finish index\n",
    "H_example = graph_to_hamiltonian(G1, G1_s, G1_f)\n",
    "print(H_example)"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}