{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bigger-pressure",
   "metadata": {},
   "source": [
    "# Dynamic programming, Jupyter notebook  (7/4 - 2021)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "polish-chorus",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "_Bitonic euclidean traveling-salesman problem_\n",
    "\n",
    "Given a set of _n_ points in the plane, we wish to find the shortest\n",
    "closed bitonic tour that connects all _n_ points, that is, tours that\n",
    "start at the leftmost point, go strictly rightward to the rightmost\n",
    "point, and then go strictly leftward back to the starting point.\n",
    "Give a dynamic programming solution assuming no two points have the\n",
    "same _x_-coordinate.\n",
    "\n",
    "Note this is Problem 15-3 from \"Introduction to Algorithms\", Thomas H.\n",
    "Cormen, Charles E. Leiserson, Ronald L. Rivest, Clifford Stein,\n",
    "page 405, The MIT Press, 2009."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "accepting-translator",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running time bitonic_tsp(...) 0.000 sec\n",
      "Most frequent calls to solve:\n",
      " [((19, 19), 1), ((0, 19), 1), ((19, 0), 1), ((18, 0), 1), ((17, 0), 1), ((16, 0), 1), ((15, 0), 1), ((14, 0), 1), ((13, 0), 1), ((12, 0), 1)]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA+hklEQVR4nO2deXhU5fXHPycQdmUTEYEkoOCCKGJYTNSq2CqoFetSbaooaGrV1r0u0bq0VGtdcLf4AxGN1l1AcStaVwRBgYgrRBJAEFAWIYJA3t8f546ZhJmsM3Nn7pzP88wzM++9M/fMTeY75573vOeIcw7DMAwjWGT4bYBhGIYRe0zcDcMwAoiJu2EYRgAxcTcMwwggJu6GYRgBxMTdMAwjgJi4G74jIk5E9oyy7X8ick6ibfKOHdUuw0h2TNyNqIjIxrBbpYj8GPa8IMprDheRZYm2tak09UdERPYTkVdFZI2I7LB4RERyRGS6iKwVkZUicq+INA/bPl5EvvDO81kRXn+J97oNIjJRRFrWYsupIvKZiPwgIp+KyMga23uLyIve9jUicmuN87A57O/8RWPPieEvJu5GVJxz7UI3oBw4Pmys2G/7koytwFPAmCjb7wdWAd2AAcAvgPPDts/3nn9U84UicjRwFTAMyAZ6AzdGOoiIdAceAy4FdgauAB4XkV297S2A14E3gN2AHt7+4VwY9nfeq7YPbSQvJu5GgxGRliIyTkS+8W7jvLG2wMvA7mGe3+4iMlhEZorIOhFZ4XmtLRp57NGeV7rW85Szw7Y5ETlPRL7yjnWfiIi3rZmI3O55ql+LyIXe/s1FZCxwKHCvZ/O9YYc8KtL71cQ594VzbgKwMIrpvYCnnHObnXMrgVeAfmGvv885NwPYHOG1o4AJzrmFzrm1wN+As6Icpwewzjn3slNeAjYBe3jbzwK+cc7d4Zzb5NmzIMp7GSmMibvRGIqAoagHegAwGLjWObcJGI6KR8jz+wbYDlwC7AIcjHqg50d649oQkROAa4DfAF2Ad4Anaux2HDAI2B84FTjaGz/Xs20AMBAYGXqBc67Ie6+Qx3phPd6voYwDThORNp53PRwV+PrQD/XsQ8wHuopI5wj7zgE+E5Ffez9oI4EtQEjAhwJLRORl74fufyLSv8Z73Oxte09EDq+njUaSYeJuNIYC4Cbn3Crn3Go0RHBGtJ2dc3Odcx8457Y555YA/0bDEg3lPOBm59xnzrltwD+AAeHeO3CLc26dc64ceBMVc1Bhvss5t8zzfm+p5zGjvV9DeRsV6Q3AMlSEX6jna9sB68Oehx7vVHNH59x2YDLwOCrqjwN/8H54QT3704C7gd2Bl4ApYVdSV6Jhn+7AeGCaiOyBkXKYuBuNYXegLOx5mTcWERHp603grRSRDago79KI42YDd3khknXA94CgQhRiZdjjClQYQzYvDdsW/rg2or1fvRGRDNRLfw5oi372jsA/6/kWG9H4eYjQ4x8iHOso4FbgcKAF+iP6fyIywNvlR+BdL2zzE3Ab0BnYB8A5N8s594Nzbotz7hHgPWBEPe00kggTd6MxfIMKbYgsbwwgUpnRB4DPgT7OuZ3R0ErE2HUdLEW90A5ht9bOuffr8doVqNcaomeN7fEsj9oJPUf3eqL5HfAw9RfNhWj4K8QBwLfe+9RkAPC2c26Oc67SOfchMAs4ytu+gIZ9Vkfj/laGz5i4G43hCeBaEekiIrsAf6Uq4+JboLOItA/bfyc0HLFRRPYG/tjI4z4IXC0i/QBEpL2InFLP1z4FXCQi3UWkAxp+COdbNBzRKERphXrLiEirULqic24N8DXwR28CtwM6Sbog7PUtvNcLkOm9PvT9nAyMEZF9vddeC0yKYsqHwKEhT11EDkQni0PHegwYKiJHiUgz4GJgDRqn7yAiR3vHbi6a7noY9Z8bMJIJ55zd7FbnDVgCHOU9boXGbFd4t7uBVmH7TgS+A9ah4ZDDUM99IzpxeRMaGgjt74A9oxz3f8A5Yc/PAErQH4ulwMRo74MK4N+9x82BOz27vkYneLcC4m0/GPgSWAvcXdf7RbAzx9s//LYkbPsA77OsRcX0KaBrjc9Z8/WHh22/FP0B2oB6/S3Dti0ECsKeXwgsQsM2pcBlNWz9jbd9g3fcft54F/TH4Qfvb/cB8Eu///fs1rhb6B/bMNIKERkOPOicy65zZ8NIQSwsY6QFItJaREZ44YbuwPXA837bZRjxwjx3Iy0QkTbAW8DeaMbIS8BFzrkNvhpmGHHCxN0wDCOAWFjGMAwjgDSve5f4s8suu7icnBy/zTAMw0gp5s6du8Y51yXStqQQ95ycHObMmeO3GYZhGCmFiJRF22ZhGcMwjABi4m4YhhFATNwNwzACiIm7YRhGADFxNwzDCCAm7oZhpCTFJcXkjMsh48YMcsblUFxibX3DMXFPU+yLYaQyxSXFFE4rpGx9GQ5H2foyCqcV2v9xGCbuaYh9MYxU5+rXi6jYWlFtrGJrBUUzinyyKPkwcU9DimbYF8NIHbZtg48/hn//G0aPhv32g6UbyiPuW74+8ng6khQrVI3EEu0LYF8Mw2+cg7IymDULZs/W+48+gh9/1O277AJDhsAyyWI9Oy7OzGqflWCLkxcT9zQkq30WZevti2H4z9q18OGHVWI+ezasWqXbWrWCgQPhvPNg8GAV9ZwcEIHikrEUTiusdgXaJrMNY4eN9eeDJCEm7mnI2GERvhjN7YthxJctW2D+/CqPfPZs+PJL3SYC++wDI0aoiA8eDP37Q2Zm5Pcq6F8AaIixbF05zSuyGP/7sT+PGybuaUnoC3DWo0Vsa1MO67MYs699MYzY4RwsWlQ9vDJvHvz0k27v1k1F/Kyz9P6gg6B9+9recUcK+hdQ0L+ACRPgnHOg18gYf4gUJymadeTm5jqrCpl4unaFX/8aXn0V+vWDl1/22yIjVVm9urpHPnu2hlwA2raF3Nwqj3zIEOjeXb31WPDDD/pj8dvfwoQJsXnPVEFE5jrnciNtM889jVm/Hjp1glGjYOxYWLoUevb02yoj2amo0OyVcCH/+mvdlpGh4ZSTT64S8333hWbN4mfPTjupsD/5JIwbp88NE/e0ZfNmjYF26ACnngp//ztMngxFlg1phLF9O3z+eXWvfMECHQfIylIRP/98vR84UD31RDNmDEycCE89pY8NE/e0Zd06vW/fHvbYAw4/HB5+GK65JnaXy0bq8c031ePkc+Zo2AP0f2XQILjqKvXIBw+G3Xbz194QBx+sE7ITJpi4hzBxT1PWr9f7Dh30/uyzNTzzzjtw2GG+mWUkkB9+gLlzq4v58uW6LTMTDjgAzjijKrzSt6+GXZIRERX1yy+Hzz5ToU93TNzTlJDnHhL3k06CCy/US1sT9+CxbRt88kn18MrChZrVAnr19otfVE14DhigeeapxBln6FXFhAlw221+W+M/Ju5pSshzD6WftW0Lp50GxcVwzz02KZXKhFZ5hgv53LlVqzw7d1YBP/nkqvBK587+2hwLdt1Vs78eeQT+8Q9o0cJvi/zFxD1Nqem5g9bteOghm5RKNdatq8paCYl5aJVny5Y6yVlYqII+ZAj06hXceZUxY+C552DaNL0aTWdM3NOUSOI+ZAjsvbeGZkzck5OfftJVnuFx8tAqT9BY8/DhVeGV/v3Ty4M9+mjNoZ8wwcTdxD1JKC4ppmhGEeXry8lqn8XYYfFdMVozLAPqzY0eDX/5i6a/7b133A5v1IPQKs9wj/zjj6tWee62mwr4qFF6n5vb8FWeQaNZM131evPNsGwZ9Ojht0X+Uefct4i0EpHZIjJfRBaKyI3e+CQR+VpE5nm3Ad64iMjdIrJIRBaIyMA4f4aUx4/66uvW6RehZk7yGWfo+KRJcTt0WlNbk5TVq+Gll+D66+GYYzQO3rcv/P736om2bAkXXQRPPw3l5Zq2+MILmr46bJgJe4jRo6Gy0v6H6yw/ICICtHXObRSRTOBd4CLgPOBF59wzNfYfAfwJGAEMAe5yzg2p7RjpXn4gZ1xOxCqN2e2zWXLxkrgc84ILdEXfmjU7bvv1r7VS39Kl0Nyu7WJG6Ec8vGBbJm04aPl4vv1vQbVVnvvtVxVaCa3ytL9F/Rk2TFfNLlqUvOmbsaC28gN1fmynbPSeZnq32n4RTgAme6/7AOggIt0aanQ64Ud99fXro3t6o0fDypXwyitxO3zasXkzXPrSjk1StlLBhzsXcdBBcOut8NZb+reZP18nt885B/bf34S9oYwZo+L+v//5bYl/1Os3TUSaicg8YBXwunNulrdprBd6uVNEWnpj3YGlYS9f5o3VfM9CEZkjInNWr17d+E+Q4jgHbbdFrqMez/rq69ZVn0wN59hjNa3s4YfjdvjAs3UrzJypNXuGDdNzvWpz5B/ryp3KefppuOIKXWPQrl1ibQ0iJ56o5zzdComFUy9xd85td84NAHoAg0VkP+BqYG9gENAJuLIhB3bOjXfO5Trncrt06dIwqwPErbfCxiljae7aVBuPd+OBdeuie+6ZmRrnnTpV48BG3VRWaknbO+7QH8dOnSAvD669Fr77TmuvdGmZ+B/xdKV1aygogGefrapOmW40KBrlnFsHvAkc45xb4YVetgAPA4O93ZYD4bUFe3hjRg2eeEJX1J3Wr4CHTxxPVvtscEK7bdmMP3583LNlonnuoOUItm2Dxx6LmwkpjXOaUXT//boYqEsXOPBAuOwyWLxYJ6afflp/HEOif+dxY2mTmdgf8XRmzBgtjlecrn3fnXO13oAuQAfvcWvgHeA4oJs3JsA44Bbv+bHAy974UGB2Xcc46KCDXLrx1lvOtWjh3GGHObd5c9X40Uc7d8AB8T9+jx7OnX127fsMHuzcfvs5V1kZf3tSgSVLnJs40bmCAue6dXNOJd65nj2dO+ss5yZPdm7Zstrf47EFj7nsO7Od3CAu+85s99iCxxJjfJpy4IHODRjgtxXxA5jjouhqfaZpugGPiEgz1NN/yjn3ooi8ISJdPBGfh2bPAExHM2UWARXA2U369Qkgn30GJ5wAvXvD889riluI3Fy45RZdKt66dfxsqC0sE2L0aO1fOXeu2pVurFwJb74Jb7yht9JSHd91VzjyyKpb7971X/EZ6h5kJIYxY7Rm0kcf6UrddKJOcXfOLQAOjDB+ZJT9HXBB000LJitXap/IFi1g+nSNzYYzaJDWyp4/H4YOjY8N27bBxo21h2VAa81cfLGuWE0Hcf/+e81WCYn5p5/qeIcOWhL54otVzPfdN7jL94PG736nobIJE0zcjTiyaRMcd5zW/XjrLa3xUZOQiH74YfzEfcMGva9L3Nu31yXcjz8Ot98e3ysJP9i4Ed59t0rMP/pIAy1t2sChh+rKzyOP1Fh6PDsJGfGjY0f9Hy4u1kqRQfsfrg0T9wSxbZt6wh9/rKsKo3nCu++uy8rjuaYrvFFHXYwerV+MF16A00+Pn02JYMsWTU8MifmsWfp3yczUZg833KBiPnhwetVjCTrnnKMOynPPaQZNumDingCcgz//GV58Ee67D44/Pvq+Ihqaiae412zUURuHHw45ORqaSTVx37ZN5wtCYv7uu7qYKCNDf1wvv1zFPD9fvXUjmPziFzovMmGCibsRY267DR54QAtynX9+3fvn5uoPwcaN8VnQEqkiZDQyMrQQ0403wpIlKvTJSmUllJRUiflbb1W1iNt/f50cPvJIXShkdVjSh4wMvQK99lpNU91jD78tSgwBrrqQHDz5pIr6b3+rlerqQ26uevsffRQfmxoSlgEVd9AmCMmEc1ru9sEH4ZRTNItlwAC49FL44gudTHvySZ3jmD8f7rxTr5pM2NOPs85SkZ840W9LEod57nHknXfgzDPhkEO0Ql19CxiF4vFz5sSn5V1DwjIA2dnq8U6aBNdd528hpvLyKs/8jTeqen726KGT1UceCUccAT171v4+RnrRvbtW2pw0Sa9C06FWj3nuceLzzzWXvVcvmDKlYf0od90VsrLiF3dvSFgmxOjRGpaJRyGm2srgrlql3vcf/gB77qk/NGefDS+/rLHyBx9U7728XL+4Z55pwm5EZswYLZP86qt+W5IY0uD3K/F8+63msmdmRs5lrw+5uZoOGQ9C4r7zzvV/zYknajhj4kT1jmNFzTK4ZevLGPN8IY9OhuWvFPDJJ1W2Hn44/OlPevx+/YJdytWIPccdp2UiJkzQ+j9Bx74eMWbTJo3rrlypk6K9ezfufXJztRZ1PIoerV+vDbAbkrvdurVmyzz7bFVYJxYUzdixDO4WV8Fr24vo1k3nKWbP1uJbU6Zos4r+/U3YjYbTooVe2U2bpg5Y0LGvSAzZvl0FcO5c+M9/NKWxsYReG49J1fqUHojE6NGaSvif/8TOlqg169uX89prWlht0KD0iJEa8WfMGE2RffRRvy2JPybuMcI59SqnTYO779ZuRk3hoIP0Ph6hmdpquddGbq52CIplxkG0crdWBteIB/vso6WYJ0zQ72yQMXGPEXfcoQuULr9cW9g1lY4dNR83HpOqdZX7jUaogfbs2bBwYWxsGXvkWDK2WxlcI3GMGaMJDzNn+m1JfDFxjwFPP62ifsop8M9/xu59c3PjI+6NDcuANvFo3jx2XZrafV1A5Qvj6ZSRjSBkt49/LXsjvTn1VF0cGPQuTSbuTeTdd7UxQ34+TJ4c24m+QYOgrCz23ZAa67mDZhscf7x+1q1bm2bHtm0aU9/rpwJWXrWEyusrWXLxEhN2I660a6eLCp98smoFcxAxcW8CX3yhuexZWQ3PZa8P4YuZYkljY+4hRo+G1bsV0/1fkXPT68vEiXp5fPPNmjZqGIlizBjNbHvySb8tiR8m7o1k1SoYPlzTCV9+GTp3jv0xBg7UOHcsxd059dybsgR/bY9i+HUhq7eW4XCUrS+jcFphgwR+0ya4/nqd3Bo5svG2GEZjGDpUJ1eDHJoxcW8EFRVVuezTpsWvENFOO8Hee8dW3Ddt0pTNpnju1/2vCDKr56ZXbK2gaEZRvd/jjjv0/P3rX9b4wkg8Iuq9f/BBVVOWoGHi3kC2b9eCVB9+qA2uhwyJ7/FivVK1MaUHahItNz1qznoNVq2CW2/VVa95eY23wzCawhlnaHJAUL13E/cG4Jy2WpsyBe66S+Pt8SY3F1as0JoYTaW4pJjBj+fA9Rlcs6pxcXJoem76TTdpj9j6Vsk0jHiw6676HZ48GX76yW9rYo+JewO48064914tKfunPyXmmLGaVA3VcFnxYxmIY822hsfJQ4wdNpY2mY3LTf/qK/j3v+Hcc2GvvRp8aMOIKWPGwJo1Gl4NGnWKu4i0EpHZIjJfRBaKyI3eeC8RmSUii0TkSRFp4Y239J4v8rbnxPkzJIRnntFGuyedpHHiRDFggE7aNlXcI9VwaWicPERB/wLGHz+e5puywTUsN/2aa6BlS51MNQy/+dWvtFx0EEMz9fHctwBHOucOAAYAx4jIUOCfwJ3OuT2BtcAYb/8xwFpv/E5vv5Tm/fd18U5entakSGTRqjZttAJiU+PuTY2T16SgfwFHzFvCoOn1z02fNUt/JC+/XPvEGobfNGumjTxeeQWWLvXbmthSp0w5ZaP3NNO7OeBI4Blv/BFgpPf4BO853vZhIqmbD/Hll1onJpTL7kf39NBK1abUwohHDZesrPp/IZyDK66Arl31CsgwkoWzz9b/z0mT/LYkttTLBxWRZiIyD1gFvA4sBtY557Z5uywDunuPuwNLAbzt64E4ZIHHn1Auu4jmsu+yiz92DBqkccHyxjnZQHxquGRlaTrjli117/vii9qZ6vrrNcXTMJKF3r21R8DEidqHNyjUS9ydc9udcwOAHsBgYO+mHlhECkVkjojMWR3r9fUxoKJCPfZvvolvLnt9CE2qNiU0Uzlfa7h0bha7Gi5ZntO/bFnt+23bBldeCX37wjnnNPpwhhE3xozRTmNvvum3JbGjQVWynXPrRORN4GCgg4g097zzHoDXzZLlQE9gmYg0B9oD30V4r/HAeIDc3NykKr65fbvG2GfP1uYUQ4f6a0///ro8f84cOPnkhr9+wwZt0j04q4CZ1xTEbM4g1M6uvLz2H79Jk+Czz/RcWpkBIxk58URd+zFhAgwb5rc1saE+2TJdRKSD97g18EvgM+BNICQ1o4Ap3uOp3nO87W84l1qVky+7DJ5/HsaN0z+637RsCfvv3/iMmb/9TcMn99wT28ngkOdeW9x90yb461/h4IOT41waRiRat4aCAnjuufh0P/OD+nzVuwFvisgC4EPgdefci8CVwKUisgiNqYeSiSYAnb3xS4GrYm92/Bg3ThcoXXIJ/PnPfltTxaBBjZtU/fxz/UyjR8PgwbG1qUcPva9tLmDcOF2EZWUGjGTnnHNgS99i9ri3aQXxkoU6wzLOuQXAgRHGS9H4e83xzcApMbEuwTz7rC5Q+s1v4Lbb/LamOrm58OCD2le1T5/6vSbUHapt2/isBm3dWlf5RRP31au1vv3IkVoS2TCSmYXNipETCllbWdWsvXBaIUBKlqG2FaoeM2dqnH3IEHjsseRrwNyYlapTpsBrr8GNN6oIx4OePaOL+9/+phPTVmbASAWKZhThmsdmoV8ykGQS5g9ffaVVHnv0gKlT/cllr4t+/bRefH3F/ccfNbTUrx+cf3787IqW675oETzwgF7q7t3k3CrDiD+xXujnN2kv7qtXw4gRVbnsXbr4bVFkmjeHAw+sfzrkv/6lqV333BPfDJWsLPXca84FFBXpRPANN8Tv2IYRS4LWrD2txf3HHzWXfdky9dj33NNvi2onNxc++khTNWtjyRINhZxyChxxRHxtysqCjRurSgmDppA+9ZRmHVmZASNVaEpBvGQkbcV9+3ZNfZo1C4qLNVUv2cnN1dTCL76ofb/LLtMrkURMCofnuoN68H/5i8b4L788/sc3jFgRKojXuZkWxOvWJrWbtaetuF9+ueay33GHZsekAoMG6X1toZn//ldzdYuKqvLQ40nNXPeXXoK33rIyA0ZqUtC/gPdPXQI3VnJjh9Ru1p6W4n7XXZp/fdFF2nwjVejbVzu3R5tU3bpVc/N7905cca6QuJeXV5UZ6NNH67UbRirSp4/WkXr/fb8taRoNKj8QBJ5/XrNITjwRbr/db2saRrNm2jQ7mrjfc48u8586VTNrEkHXrjphW14Ojzyi/SitzICRyohoee/33vPbkqaRVp77Bx9o/9PBgzWXvVkzvy1qOIMGwbx56qWHs3KlZqYMHw7HHZc4ezIyNIX088+tzIARHPLzNUU6CWsa1pu0EfdFizSXvXt3rfLYpk3dr0lGcnNh82ZYuLD6+FVX6fi4cYld5l9cUsw3p+YwZUAG35yaw1GXFFuZASPlCa2oTuXQTFqI+5o1msvuHEyfnry57PUh0krVmTM1JHLppRqXTxShvqxbWmtfVjqUcftXjevLahjJxEEHQYsWqR2aCby4//ijdjgvL9dYdCLFLx7ssYeWJg2J+/bt2qx7993h2msTa0ss+7IaRjLRqpUKvHnuSUplJZxxhnq2xcU6SZLqiKj3HkqHnDgR5s7VFant2iXWlqAt1zaMcPLy1ImqT6exZCTQ4n7FFZq5cdttcNJJflsTO3JzoaREJ1GvuQYOPRROPz3xdgRtubZhhJOfr8L+0Ud+W9I4Aivu99yjC5T+9CdNfQwSFXsWs/WCHLo9mMGa3+cw/C/+TGIGbbm2YYQTutJP1bh7IMV9yhRdoHTCCXDnncFqElFcUsxDKwqhQ9Uk5t8X+DOJGVqund0+dn1ZDSNZ6NpV57hSVdwlGTrg5ebmujmN7SFXg1mztFhW//7a7DZVUx6jkTMuh7L1ZTuMZ7fPZsnFSxJvkGEEmFGj4JVXNASajE6iiMx1zuVG2hYoz33xYs1l79YttXPZa8MmMQ0jceTlwapVqi2pRmDEfc0aXZ25fbvWZY9X5yG/sUlMw0gcqbyYKRDivnmz9uksL9d4e6rnsteGTWIaRuLYd19o3z414+4pK+7FJcXkjNMu5Z3/lsN7G4p59FE45BC/LYsvNolpGIkjI0PrJaWiuNdZFVJEegKTga6AA8Y75+4SkRuAc4FQaZ1rnHPTvddcDYwBtgN/ds69GkujQ8veQ6sjK1qUkXlSIT/tDRB8kSvoX2BibhgJIj8frrtOu4116OC3NfWnPp77NuAy59y+wFDgAhHZ19t2p3NugHcLCfu+wGlAP+AY4H4RiWn9xUjL3rdiy94Nw4g9oXz3mTP9taOh1CnuzrkVzrmPvMc/AJ8B3Wt5yQnAf5xzW5xzXwOLgMGxMDaEZYwYhpEohgzR8uCpNqnaoJi7iOQABwKzvKELRWSBiEwUkY7eWHdgadjLlhHhx0BECkVkjojMWd3AosmWMWIYRqJo2xYGDEi9uHu9xV1E2gHPAhc75zYADwB7AAOAFUCD+ho558Y753Kdc7ldGliD1zJGDMNIJHl5ukCyZpOcZKZe4i4imaiwFzvnngNwzn3rnNvunKsEHqIq9LIc6Bn28h7eWMz4OWOk1W6Ig+wWXSxjxDCMuJGfDxUVsGCB35bUnzrFXUQEmAB85py7I2y8W9huJwKfeI+nAqeJSEsR6QX0AWbHzmSloH8BS0a8RuWNsKTP/SbshmHEjVQsIlYfzz0fOAM4UkTmebcRwK0iUiIiC4AjgEsAnHMLgaeAT4FXgAucc9vjYn2oC/TmzXF5e8MwDICePfWWSpOqdea5O+feBSKVzJley2vGAvEPgLdurfc//hj3QxmGkd7k58O77/ptRf1J2RWqgHnuhmEkjLw8WLZMy5ykAqkt7ua5G4aRIFKtiFhqi7t57oZhJIj999ec91SZVE1tcW/WDDIzTdwNw4g7zZvralXz3BNFq1YWljEMIyHk58P8+bBxo9+W1E3qi3vr1ua5G4aREPLytCHQrFl17+s3qS/u5rkbhpEghg7VXqqpEJoJhrib524YRgLo0AH69UuNSdXUF/fWrc1zNwwjYeTna233ykq/Lamd1Bd389wNw0gg+fmwYQMsXOi3JbWT+uJuE6qGYSSQVCkilvribhOqhmEkkN69oWvX5J9UTX1xN8/dMOJOcUkxOeNyyLgxg5xxORSXFPttkm+IqPdunnu8Mc/dMOJKcUkxhdMKKVtfhsNRtr6MwmmFaS3w+flQWgorV/ptSXSCIe7muRtG3Cj67zVUbK2oNlaxtYKiGUU+WeQ/qVBELPXF3VIhDSM+/PQTTJhA+frINW6jjacDBx4ILVsmd2gm9cXdPHfDiC1btsC//w19+8I555D1Y4uIu2W17xlxPB1o2RIGDTLPPb7YhKphxIbNm+G++2DPPeG882C33WD6dMYWTKRNZptqu7b5CcYu2DWtr5rz8mDu3OQ9Bakv7q1aaSWfrVv9tsQwUpOKCrjrLs3xu/BCyMmB117TZZjDh1OwfwHjjx9PdvtsBCG7fTbjO5xBwaS5cOSRsGaN35/AF/LzVXbmzvXbksjU2UM16Ql1Y9q8WWu7G4ZRPzZtggcegNtug2+/hcMPh+JivZfqbZML+hdQ0L+g+uv3GAkFBerCvvwy7LFHoixPCsIXMx1yiL+2RKJOz11EeorImyLyqYgsFJGLvPFOIvK6iHzl3Xf0xkVE7haRRSKyQEQGxvUThLoxJeu1kWEkGz/8ALfcoh76FVdA//7w1lvw5ptwxBE7CHtUfvMbmDEDvv8eDj44NergxpBddtFpiWSdVK1PWGYbcJlzbl9gKHCBiOwLXAXMcM71AWZ4zwGGA328WyHwQMytDifcczcMIzrr18Pf/66ifvXVVTOCr78Ohx3WuPfMy9P32Gkn/WGYMiWmJic7+fn68Z3z25IdqVPcnXMrnHMfeY9/AD4DugMnAI94uz0CjPQenwBMdsoHQAcR6RZrw3/GPHfDqJ21a+GGGyA7G667ThVp9myYPl097qbSt6/G5/v3V2/+vvua/p4pQl4efPcdfPml35bsSIMmVEUkBzgQmAV0dc6t8DatBLp6j7sDS8Netswbq/lehSIyR0TmrF69uqF2V2FNsg0jMt99B9deq6J+4406+Tl3Lkydql57LNl1V3jjDTj2WJ2UvfLK5K+JGwOSeTFTvcVdRNoBzwIXO+c2hG9zzjmgQRcmzrnxzrlc51xuly5dGvLS6oTCMua5G4ayahVcdZWGX/7xDzjmGG38+dxzMDCOU2Bt28Lzz8Mf/wi33qqTrVu2xO94ScBee0GnTskZd69XtoyIZKLCXuyce84b/lZEujnnVnhhl1Xe+HIgfHVDD28sPpjnbhjKypWa+fLAA/p9+O1voahIWwclimbNNCyTna0/MCtWqOB37Jg4GxJIRoZGtpJR3OuTLSPABOAz59wdYZumAqO8x6OAKWHjZ3pZM0OB9WHhm9hjE6pGurN8OVx0EfTqBXfeCSefDJ9+Co8/nlhhDyGiYZniYo1X5OdDWVni7UgQ+fnw+ecaBUsm6hOWyQfOAI4UkXnebQRwC/BLEfkKOMp7DjAdKAUWAQ8B58fe7DBsQtVIV8rL4YILdPHR/ffD734HX3wBjzyi8QK/+d3vdDHUN99oZ+mPP/bborgQynefOdNfO2pSZ1jGOfcuEC3xdViE/R1wQRPtqj/muRvpxpIlcPPN8PDD+vzsszUE0quXr2ZF5PDDNWYxfLimWz7zDBx9tN9WxZRBg6B5c71IOe44v62pIhjlB8A8dyP4LF4MY8ZAnz4waRKcey4sWqRFvpJR2EP06wcffKA1a449FiZO9NuimNKmjc5TJ1vcPTjibp67EVS+/BJGjdJQy+OPw/nna6eI++6DrCy/rasfu+8Ob78Nw4bpD9T11yfnyp9GkpenSwd++slvS6pIfXG3VEgjQFRrZ/ev3Sk+Lw/22QeeflonTUtLtchX9x2WjiQ/O+0EL76oYaSbboLRowNT8C8/X/3LefP8tqSK1C8cZp67ERBC7exCXY/KKlZQuMsKuOo4Ci6aoAuFUp3MTJgwQVMlb7hBM32eeQZ23tlvy5pEeBGxwYP9tSVE6nvumZmabGribqQ4RTOKdmxnlwlFXUqCIewhRDQsM3GiFis77DAV+RRm9911zVgyrVRNfXEXsSbZRiBIu3Z2Z58NL72kE8UHHwwLF/ptUZPIz1fPPVmmElJf3MG6MRmBIKt95MnRaOOB4Fe/gnfegW3bVB3ffNNvixpNXp4uyF2yxG9LlGCIu3nuRgAYO2wsbVz1abA2mW0YO2ysTxYliAEDNFWye3fNgS8u9tuiRpFsRcSCIe7muRsBoKB/AeOX7E/2psyqdnbHj9+xA1IQycqCd99Vhfz973WRVrLEN+rJfvtpQlCy5LunfrYMmOduBIaCuT9R0Ht42jW9ALS42CuvaIrkNddoeYV77tHlnylAs2ZaZcE891jSqpV57kbq45zmsadZL9JqtGwJjz6q5RQefBBOPFF7vaYI+flQUgIbNtS9b7wJhrhbWMYIAqtWQUWFFgJLZzIyNCxz//3aLeqII7SBdwqQl6c9Sj74wG9LgiLuFpYxgsDixXqf7uIe4o9/1Frwn3yiqZJffOG3RXUyZIj+NiVDaCYY4m6euxEESkv13sS9il//Gv73P9i4Ud3iZJmtjMLOO2sr2WQwMxjibp67EQRKS3VRXk6O35YkF4MHa7H0zp218Nizz/ptUa3k52tYZvt2f+0Ihrib524EgdJSzfUO1UsyqthjD411DBwIp5wC48b5bVFU8vP1QqOkxF87giHu5rkbQaC01EIytbHLLjBjhmbQXHKJ3ior/bZqB8KLiPlJcMTdPHcj1Vm82MS9Llq3hqee0vLH48bBqacmnWOXna2FxPyeVE2N1QF1YWEZI9X58UftNWriXjfNmqmwZ2fDpZdqQZepUzUmnwSIJMfcb3A8959+8n8GwzAaS6jaVDovYGool1yiXvzcuaqmoWyjJCA/H8rK/K1kXKe4i8hEEVklIp+Ejd0gIstFZJ53GxG27WoRWSQiX4hIYjrhhroxbdmSkMMZRsyxNMjGccop8N//wpo1mgv/4Yd+WwQkRxGx+njuk4BjIozf6Zwb4N2mA4jIvsBpQD/vNfeLSLNYGRsVa5JtpDq2gKnxHHKIqmibNnD44drKz2cGDFCf08/QTJ3i7px7G/i+nu93AvAf59wW59zXwCIg/k2nQp67xd2NVKW0FNq2hS5d/LYkNdlrL82F32cfOOEE+Pe/fTUnM1PT85Pdc4/GhSKywAvbdPTGugNLw/ZZ5o3tgIgUisgcEZmzevXqJpiBee5G6hNKgxTx25LUZbfddDXrMcfAeefB1Vf7miqZlwcff6zlgvygseL+ALAHMABYAdze0Ddwzo13zuU653K7NNVbMc/dSHXSvRpkrGjXTsslFxbCLbfAmWdqsoUP5Odrgym/pgEaJe7OuW+dc9udc5XAQ1SFXpYDPcN27eGNxRfz3I1UJlTq1+LtsaF5cy0X/I9/aFenY46BdesSbsbBB+u9X3H3Rom7iHQLe3oiEMqkmQqcJiItRaQX0AeY3TQT60FI3M1zN1KRlSvVMTFxjx0iGpZ59FHt8HTIIbB0ad2viyGdOukUQNKKu4g8AcwE9hKRZSIyBrhVREpEZAFwBHAJgHNuIfAU8CnwCnCBcy7+yecWljFSGUuDjB+//712d1q6VNskzZ+f0MPn5+s8rx+h//pky5zunOvmnMt0zvVwzk1wzp3hnOvvnNvfOfdr59yKsP3HOuf2cM7t5Zx7Ob7me1hYxkhlQuJuMff4cOSR6r1nZMChh8Lrryfs0Hl5sHYtfP55wg75M8FYoWqeu5HKhEr9Zmf7bUlw6d9fXehevWDECJg0KSGH9XMxUzDE3Tx3I5VZvBh69ND+oUb86NED3nlHFzqdfTbcdJNOZseRPn20mKUfcfdgiLt57kYqY5kyiWPnneGllzRF8vrr4dxzYevWuB3OzyJiwRB389yNVMZy3BNLixYalrnuOpgwQVv5/fBD3A6Xnw9ffQVNXavZUIIl7ua5G0lAcUkxOeNyyLgxg5xxORSXFEffuaJCS9aa555YRDQs89BDOsH6i1/o3yEOhJp3JDrubuJuGDGkuKSYwmmFlK0vw+EoW19G4bTC6AL/9dd6b+LuD+ecA9OmwZdfaqrkp5/G/BC5uXqxkGhxD0azDhGdjLKwjJEotmzRgt1ff6230lL4+muKek+hok31GG7F1gqKZhRR0L9gx/exHHf/GT4c3n4bjj1WYygvvKCefIxo1QoOOijxcfdgiDtYNyYjtlRW6mV6DfH++X758uqZFi1aQE4O5f0iT86Vry+PfBwT9+Rg4EBNlRwxAn71K3jkETjttJi9fV4e3Huv+gSJSooKjrhbk2yjoaxbF128lyyp3vxFRBtj9u6ti2J699ac6V699HG3bpCRQda4HMrWl+1wqKz2WZFtKC2FnXbSfDnDX3Jy1L0eORJOP11XtV5+eUwqdebnw+23w0cfVdWciTfBEXfz3I2ahIdOaor311/r0sFwOnRQod5vPzj++CoB791bFxjVw+UaO2wshdMKqdhaVee1TWYbxg4bG/kFVuo3uejYEV59FUaNgr/8RX/k775b+7Y2gdCk6nvvmbg3HPPc049Q6CRcsMMfRwqdhLztoUOre969eqm4N5FQXL1oykWUb/uOrHbdGXv0PyPH20EXMO29d5OPa8SQVq3giScgKwtuu03/jx5/XDs9NZKuXTXbNZGTqsESd/Pcg8e6ddHFO1LopHt3Fephw3YUby90Em8K+hdQUNpOL+/nToX+AyPvWFmpn2PEiMjbDf/IyIB//Uuv2P78Zw3FTZvWpE5Z+flaw8y5xFyoBUfcLSyTmoRCJ5HCJqWlO9bh7thRhbp/f118Ei7e9QydJIROnfT++1o6VK5cqf+ztoApebnwQi1bcPrpGk95+WWtKdAI8vJg8mS9WNtzzxjbGYHgiLuFZZKTykr45pvInndpqW4LD520bKkTW717V4VOwicvYxA6SQidO+v9d99F38cyZVKDkSPhzTd1HiYvD6ZObVTgPLyImIl7Q2jdGtav99uK9GTt2ujiXVYWOXTSuzccddSO4p2g0EncqY/nvnix3pu4Jz9Dh2qq5DHHaIjm8cfhxBMb9Bb77gvt2+uk6plnxsnOMIIj7ua5x4/Nm3fMOgl/HCl00rs37L+/dqIPzzrJykqe0Ek8CYl7XZ57RoaV+k0V9txTBf744+Gkk+Cuu+BPf6r3yzMy1OFP1KRqcMTdYu47UFxSTNGMIsrXl5PVPouxw8ZGztoID51EEu/lNdrgtmxZ5WkffHB18e7VS92TdKdFC23WXJvnXloKPXvqvkZq0KULvPEG/O53OtFaVga33lrvq838fPjrX9UfineEMTjibp57NUI1TkL51mXryyiccg7Mmk3B2h47LtgJ7xAvopNIvXpp6KTmgp3ddgtG6CTedOpUt+duIZnUo00bePZZuPhiXZm0dKmuaA3VuKqFvDydYpo5U6sexJPgiLt57tUomlFUbSENQMX2zRQtvJuCcajw9OqloZORI6t73ukSOok3nTvXHXM/7rjE2WPEjmbNdHFTdjZccYWut3jhhapwXBSGDNGXvv++iXv9sTz3akSrZVLeQWDdWgudJIJOnaKL+6ZN8O235rmnMiJanqBnT50hzc/XVMmcnKgvadsWBgxITBGxOq+tRWSiiKwSkU/CxjqJyOsi8pV339EbFxG5W0QWicgCEYmyeiMOhMIycW6blSpkSYfI4+2zTNgTRefO0cMyVuo3OPz2t1oTfuVKzaqZO7fW3fPyYNYs2LYtvmbVJ3A6CTimxthVwAznXB9ghvccYDjQx7sVAg/Exsx6EGq1Fx47TlemT2fsM2tpU1m9HkatNU6M2FOb5x7KcbcFTMHgsMPUHW/ZUssFT58eddf8fO3RMn9+fE2qU9ydc28DNf9DTwAe8R4/AowMG5/slA+ADiLSLUa21o612lM++QROO42CzIGMP2482e2zEYTs9tmMP3589BonRuwJxdwjXU3aAqbgse++8MEH0Levrp5+6KGIu4UXEYsnjY25d3XOhXpSrQS6eo+7A0vD9lvmjcWnf1U41iQbVq3SCbp27WDqVAq6d6dg0Gi/rUpfOnWC7dthw4YdQ2GLF2uz5jom4IwUo1s3bfxxyilQWAjl5drOL6yYTM+eenv/fc2mjBdNzmdzzjmgwYFuESkUkTkiMmd1LDrHprvnvnmzrphbtUqXR3fv7rdFRm0LmazUb3DxnCvGjIG//13LB9cIF+fnx99zb6y4fxsKt3j3q7zx5UDPsP16eGM74Jwb75zLdc7ldmlCpbWfSWfP3Tk491x1BSZP1qaNhv+E6stEiruXllq8PchkZmpY5qab4NFHtfJnWHmUvDxYtkxT5ONFY8V9KjDKezwKmBI2fqaXNTMUWB8Wvokv6ey533wzPPYY/O1vcPLJfltjhIjmuYdK/Vq8PdiIwHXXwaRJ8NZbcOihquhUFRGLp/den1TIJ4CZwF4iskxExgC3AL8Uka+Ao7znANOBUmAR8BBwflysjkRI3NPNc3/2WSgqgoICvTeSh2ie+zffaDE1E/f0YNQozZ5ZskRTJUtK2H9/zXmPp7jXOaHqnDs9yqZhEfZ1wAVNNapRpGNYZu5cOOMMre/yf/9n8dtkI5rnbpky6ccvfwnvvKPhmUMOoflzzzFkyLC4FhELToGQdAvLLF+u6VZdusDzz9erroWRYKKV/bUc9/TkgAM0VTIrC445hp2zbuajw3LIuDGDnHE5FJcUx/RwwSk/kE6e+6ZNKuwbNugkateudb/GSDzNm2u6YyRxz8jQL7mRXvTsCe+8Q/Ef8nh592uhRSUOr7DftEKAmK1FMc891ais1DoWH3+sTXz79/fbIqM2IpUgWLxYhT0z0x+bDH/p0IGioZvY0qKy2nDF1gqKZsRu3sw891Tjuuvguee01KhVFEx+IpUgsFK/aU/5hsg5kNEK/jUG89xTiUcfhX/8Q3PaL7nEb2uM+hCppruJe9qT1T5ySC7aeGMIjrgH3XN/7z045xw44gi47z7LjEkVatZ037hRVxHbZGpaM3bYWNpktqk2FuvCfsER91BziSCK+5IlWlogOxueecZitalETc/dSv0a6KTp+OPjW9gvODH35s31FrSwzIYNGlvfuhVefNEKTaUanTvD2rU6EZ6RoZOpYOJuUNC/IK5VWoPjuUPwWu1t3w6nnw6ff64ee9++fltkNJROnbT2z7p1+twWMBkJIljiHrQm2ZdfrsuW77sPhu2wINhIBWqWICgt1bb3dgVmxJlgiXuQPPfx42HcOLjoIvjDH/y2xmgsNUsQWKaMkSCCJe5N9NyLS4rJGRe/5cD15o034IILtD367bf7Y4MRG2p67osXm7gbCSFY4t4Ez724pJjCaYWUrS/D4X5eDpxwgf/ySzjpJNhrL/jPf6BZs7pfYyQv4fVltm/XzCcTdyMBBEbci0uKyTnmMzIGTGm41715M0WvXUnF1opqw7FeDlwn33+vmTHNm8O0aVqXxEhtwsMy33yjHXksx91IAIFIhQx53RWttZVV2foyCl8YA++9T8G2fTQV7fvv9T78ceh+82bKrwcirAuK5XLgWtm6VfsulpVpWKZXr8Qc14gvHTvq/fffW6aMkVACIe5FM4p29Lort1D0xf0UjPMGdtpJv2gdO6o3tddeVY87diRr222UVe7YDi3rB9FWWeedB7vuGp8P4BxceKGK+uTJVW1ajNSnWTPNjvnuO8txNxJKIMQ9mndd3kFg1bf65apjVefYkiz1/sN+JNpktGTst3vD7ddrTZeCArj44thXYrzrLs2Oufpqbb5hBItQCYLSUhX7nj3rfo1hNJFAxNxrLcLTpUu9lutHXA48cgIFj86Dzz6D0aO1xO7++8NRR+lq0crKOt+3Tl56CS67DH7zG+2UbgSPUAmC0lIr9WskDNHOeP6Sm5vr5syZ0+jX/xxzD/e6M9vEvFYD33+vHc3vuUc7IfXpo3noo0ZBu3YNf7+SEm2D3rcvvP22NlU0gsfw4SruGRkaHnz9db8tMgKCiMx1zuVG2hYIzz0RRXgA9cCuvFKLPz3xhMbsL7xQL7P/8hcob8Dk66pVcPzxmhEzdaoJe5AJ99wt3m4kiEDE3CH+RXiqkZkJp52mt5kzdSXpHXfo7aSTNC5/8MHRX795M4wcqQL/zjvQvXti7Db8oXNnvdLbssXE3UgYTfLcRWSJiJSIyDwRmeONdRKR10XkK+++Y2xMTVIOPhiefFK9sksvhVdf1VDL0KG6CGnrVqDG6te/d6Z440zNjDnoIJ8/gBF3OnVSYQcTdyNhxCIsc4RzbkBY3OcqYIZzrg8ww3sefLKy4NZbYdkyuPdejc+ffjr07k3x2NMonHpu1erXzAoKT8qkeK8tflttJILwImG2gMlIEPGIuZ8APOI9fgQYGYdjJC/t2mldmM8/14yavfaiaM2TVGyrXvOmgq2JXf1q+EeovgyY524kjKaKuwNeE5G5IlLojXV1zq3wHq8EukZ6oYgUisgcEZmzevXqJpqRhGRkwLHHwn//q/n2EUjY6lfDX0Kee8eOuubCMBJAU8X9EOfcQGA4cIGIHBa+0WmeZcRcS+fceOdcrnMut0uXLk00I7lJRDNcI3kp3vwhORdDxp/X+ltt1EgrmiTuzrnl3v0q4HlgMPCtiHQD8O5XNdXIVCcRzXCN5KS4pJjChTdT1gGc4F+1USPtaLS4i0hbEdkp9Bj4FfAJMBUY5e02CpjSVCNTnYTl4RtJR9GMIiq2Vy9DnfBqo0Za0pQ8967A8yISep/HnXOviMiHwFMiMgYoA05tupmpT0Lz8I2kIWrdI5tvMeJMo8XdOVcKHBBh/DvAGn4aBjqvUra+LOK4YcSTQJQfMIxkxeZbDL8wcTeMOGLzLYZfBKIqpGEYRjoS+KqQhmEYRnVM3A3DMAKIibthGEYAMXE3DMMIICbuhmEYASQpsmVEZDW6mjWV2AVY47cRjSAV7TabE0Mq2gypaXesbM52zkWsvJgU4p6KiMicaClIyUwq2m02J4ZUtBlS0+5E2GxhGcMwjABi4m4YhhFATNwbz3i/DWgkqWi32ZwYUtFmSE27426zxdwNwzACiHnuhmEYAcTE3TAMI4CYuNcTEVkiIiUiMk9E5nhjnUTkdRH5yrvv6LONE0VklYh8EjYW0UZR7haRRSKyQEQGJpndN4jIcu98zxOREWHbrvbs/kJEjvbJ5p4i8qaIfCoiC0XkIm88ac93LTYn7bkWkVYiMltE5ns23+iN9xKRWZ5tT4pIC2+8pfd8kbc9J4lsniQiX4ed5wHeeHz+N5xzdqvHDVgC7FJj7FbgKu/xVcA/fbbxMGAg8EldNgIjgJcBAYYCs5LM7huAyyPsuy8wH2gJ9AIWA818sLkbMNB7vBPwpWdb0p7vWmxO2nPtna923uNMYJZ3/p4CTvPGHwT+6D0+H3jQe3wa8KQP5zmazZOAkyPsH5f/DfPcm8YJwCPe40eAkf6ZAs65t4HvawxHs/EEYLJTPgA6iEi3hBhagyh2R+ME4D/OuS3Oua+BRcDguBkXBefcCufcR97jH4DPgO4k8fmuxeZo+H6uvfO10Xua6d0ccCTwjDde8zyHzv8zwDDxGj0nilpsjkZc/jdM3OuPA14TkbkiUuiNdXXOrfAer0Sbhicb0WzsDiwN228ZtX/R/eBC7zJ1YljIK+ns9i79D0Q9tJQ43zVshiQ+1yLSTETmAauA19EriHXOuW0R7PrZZm/7eqBzQg1mR5udc6HzPNY7z3eKSMuaNnvE5DybuNefQ5xzA4HhwAUiclj4RqfXV0mdV5oKNobxALAHMABYAdzuqzVREJF2wLPAxc65DeHbkvV8R7A5qc+1c267c24A0AO9ctjbX4vqpqbNIrIfcDVq+yCgE3BlPG0wca8nzrnl3v0q4Hn0n+zb0OWTd7/KPwujEs3G5UDPsP16eGNJgXPuW+8LUgk8RFU4IGnsFpFMVCSLnXPPecNJfb4j2ZwK5xrAObcOeBM4GA1dNI9g1882e9vbA98l1tIqwmw+xguLOefcFuBh4nyeTdzrgYi0FZGdQo+BXwGfAFOBUd5uo4Ap/lhYK9FsnAqc6c3UDwXWh4UTfKdGzPFE9HyD2n2alxXRC+gDzPbBPgEmAJ855+4I25S05zuazcl8rkWki4h08B63Bn6JzhW8CZzs7VbzPIfO/8nAG94VVMKIYvPnYT/6gs4RhJ/n2P9vJGL2ONVvQG80a2A+sBAo8sY7AzOAr4D/Ap18tvMJ9LJ6Kxq3GxPNRnRm/j40flkC5CaZ3Y96di3w/vm7he1f5Nn9BTDcJ5sPQUMuC4B53m1EMp/vWmxO2nMN7A987Nn2CfBXb7w3+kOzCHgaaOmNt/KeL/K2904im9/wzvMnwGNUZdTE5X/Dyg8YhmEEEAvLGIZhBBATd8MwjABi4m4YhhFATNwNwzACiIm7YRhGADFxNwzDCCAm7oZhGAHk/wEDwmrjYMJgDwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "'''\n",
    " Bitonic euclidean traveling-salesman problem\n",
    "\n",
    "Given a set of n points in the plane, we wish to find the shortest\n",
    "closed bitonic tour that connects all n points, that is, tours that\n",
    "start at the leftmost point, go strictly rightward to the rightmost\n",
    "point, and then go strictly leftward back to the starting point.\n",
    "Give a dynamic programming solution assuming no two points have the\n",
    "same x-coordinate.\n",
    "\n",
    "Note this is Problem 15-3 from \"Introduction to Algorithms\", Thomas H.\n",
    "Cormen, Charles E. Leiserson, Ronald L. Rivest, Clifford Stein,\n",
    "page 405, The MIT Press, 2009.\n",
    "'''\n",
    "\n",
    "import math\n",
    "from random import randint\n",
    "import matplotlib.pyplot as plt\n",
    "from functools import cache\n",
    "from time import time\n",
    "\n",
    "def generate_points(n):\n",
    "    '''Generate n points with different x-values, worted wrt x.'''\n",
    "    \n",
    "    while True:\n",
    "        points = [(randint(1, n ** 2), randint(1, n ** 2)) for _ in range(n)]\n",
    "\n",
    "        xs, ys = zip(*points)\n",
    "        if len(set(xs)) == len(xs):\n",
    "            break\n",
    "\n",
    "    return sorted(points)\n",
    "\n",
    "\n",
    "def dist(p, q):\n",
    "    '''Compute difference between two points in the plane.'''\n",
    "    \n",
    "    px, py = p\n",
    "    qx, qy = q\n",
    "    \n",
    "    return math.sqrt((qx - px) ** 2 + (qy - py) ** 2)\n",
    "\n",
    "\n",
    "def trace_time(f):\n",
    "    '''Decorator to print running time of function calls.'''\n",
    "    \n",
    "    def wrapper(*args):\n",
    "        start = time()\n",
    "        answer = f(*args)\n",
    "        finish = time()\n",
    "        print(f'Running time {f.__name__}(...) {finish - start:.3f} sec')\n",
    "        return answer\n",
    "\n",
    "    return wrapper\n",
    "\n",
    "\n",
    "def trace_calls(f):\n",
    "    '''Decorator to trace recursive calls.'''\n",
    "    \n",
    "    indent = 0\n",
    "    def wrapper(*args):\n",
    "        nonlocal indent\n",
    "        print(' ' * indent + f'{f.__name__}({\", \".join(map(repr, args))})')\n",
    "        indent += 1\n",
    "        answer = f(*args)\n",
    "        indent -= 1\n",
    "        return answer\n",
    "    return wrapper\n",
    "\n",
    "\n",
    "def memoize(f):\n",
    "    answers =  {}\n",
    "    def wrapper(*args):\n",
    "        if args not in answers:\n",
    "            answers[args] = f(*args)\n",
    "        return answers[args]\n",
    "\n",
    "    return wrapper\n",
    "\n",
    "calls = {}\n",
    "\n",
    "@trace_time\n",
    "def bitonic_tsp(points):\n",
    "    def d(i, j):\n",
    "        '''Distance between points[i] and points[j]'''\n",
    "        \n",
    "        return dist(points[i], points[j])\n",
    "\n",
    "    #@trace_calls\n",
    "    #@memoize\n",
    "    @cache  # Since Python 3.9\n",
    "    def solve(i, j):\n",
    "        '''Find two left-to-right paths from points[0] to points[i},\n",
    "           and points[0] to points[j], visiting all nodes 0..max(i,j).\n",
    "\n",
    "           Returns (length, path_i, path_j) where\n",
    "           length = sum of the lengths of the two paths\n",
    "           path_i = tuple of node ids on path 0..i\n",
    "           path_j = tuple of node ids on path 0..j\n",
    "        '''\n",
    "\n",
    "        calls[i, j] = calls.get((i, j), 0) + 1\n",
    "        if i == j == 0:\n",
    "            length, path_i, path_j = 0, (0,), (0,)\n",
    "        elif i < j:\n",
    "            length, path_j, path_i = solve(j, i)\n",
    "        elif i > j + 1:\n",
    "            length, path_i, path_j = solve(i - 1, j)\n",
    "            length += d(i - 1, i)\n",
    "            path_i = (*path_i, i)\n",
    "        else:  # i == j or i == j + 1\n",
    "            length = math.inf\n",
    "            for k in range(i):\n",
    "                length_, path_k_, path_j_ = solve(k, j)\n",
    "                length_ = length_ + d(k, i)\n",
    "                if length_ < length:\n",
    "                    length, path_i, path_j = length_, (*path_k_, i), path_j_\n",
    "\n",
    "        return length, path_i, path_j\n",
    "\n",
    "    return solve(len(points) - 1, len(points) - 1)\n",
    "\n",
    "n = 20\n",
    "\n",
    "points = generate_points(n)\n",
    "length, A, B = bitonic_tsp(points)\n",
    "\n",
    "print('Most frequent calls to solve:\\n',\n",
    "      sorted(calls.items(), key=lambda e: e[1], reverse=True)[:10])\n",
    "\n",
    "for path, color in [(A, 'r'), (B, 'b')]:\n",
    "#    for i, j in zip(path, path[1:]):\n",
    "#        plt.plot(*zip(points[i], points[j]), '-' + color)\n",
    "    plt.plot(*zip(*[points[i] for i in path]), '-' + color)\n",
    "        \n",
    "plt.title('Total length %.3f' % length)\n",
    "plt.plot(*zip(*points), 'og')\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}