naive_sinkhorn.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "#### Este notebook contiene una implementación ingenua del algoritmo de Sinkhorn\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.spatial.distance import cdist\n",
    "\n",
    "\n",
    "def Sinkhorn_Knopp(mu, nu, C, reg = 1e-1, iteraciones_maximas = 1e10, error_minimo=1e-1):\n",
    "    '''\n",
    "    Esta función realiza el algoritmo de sinkhorn-Knopp para calcular\n",
    "    la solución al problema de transporte óptimo (Kantorovich) regularizado.\n",
    "\n",
    "    input:\n",
    "            mu:  Distribución source\n",
    "            nu:  Distribución target\n",
    "            C:   Matriz de costo\n",
    "            reg: Factor de regularización\n",
    "            iteraciones_maximas: Cantidad máxima de iteraciones\n",
    "            error_minimo: Error minimimo a partir del cual se considera \"buena\" la solución\n",
    "    \n",
    "    output: \n",
    "            T: Matriz de transporte\n",
    "            error: Error obtenido con T\n",
    "    '''\n",
    "   \n",
    "    v = np.ones_like(mu)\n",
    "    xi = np.exp(- C / reg)\n",
    "    \n",
    "    error = 1\n",
    "\n",
    "    iteraciones = 0\n",
    "    # A continuación cominenza el proceso iterativo\n",
    "    while error > error_minimo:\n",
    "        u = np.divide(mu, np.dot(xi.T, v))\n",
    "        v = np.divide(nu, np.dot(xi, u))\n",
    "\n",
    "        # Medimos que tan lejos está T de mu y T de nu\n",
    "        diag_u = np.diag(u)\n",
    "        diag_v = np.diag(v)\n",
    "        T = np.dot(np.dot(diag_u, xi), diag_v)\n",
    "\n",
    "        error = np.linalg.norm(T.sum(axis=1) - mu, 1) + np.linalg.norm(T.sum(axis=0) - nu, 1)\n",
    "\n",
    "        iteraciones += 1\n",
    "\n",
    "        if iteraciones % 1000 == 0: \n",
    "            print(f\"Iteración {iteraciones}, Error {error}\")\n",
    "\n",
    "        if iteraciones > iteraciones_maximas:\n",
    "            print(f\"Límite de iteraciones, error alcanzado: {error}\")\n",
    "            break\n",
    "        \n",
    "    return T, error\n",
    "\n",
    "\n",
    "def plot_transporte(mu, nu, T, threshold=1e-8):\n",
    "    '''  \n",
    "    Función para plotear el transporte óptimo en el caso de distribuciones bidimensionales\n",
    "\n",
    "    input:\n",
    "            mu:  Distribución source\n",
    "            nu:  Distribución target\n",
    "            T:   Matriz de transporte\n",
    "            tthreshold: límite al partir del cual no graficamos el segmento\n",
    "    \n",
    "    output: \n",
    "            T: Matriz de transporte\n",
    "            error: Error obtenido con T\n",
    "    '''\n",
    "\n",
    "    plt.figure(figsize=(12, 8))\n",
    "\n",
    "    # Puntos\n",
    "    plt.scatter(mu[:, 0], mu[:, 1], label='mu', color='blue')\n",
    "    plt.scatter(nu[:, 0], nu[:, 1], label='nu', color='red')\n",
    "\n",
    "    # Segmentos\n",
    "    mx = T.max()\n",
    "    for i in range(nu.shape[0]):\n",
    "        for j in range(mu.shape[0]):\n",
    "            if T[i, j] / mx > threshold:\n",
    "                plt.plot([nu[i, 0], mu[j, 0]], [nu[i, 1], mu[j, 1]],\n",
    "                        alpha = T[i, j] / mx, color='black')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prueba simple"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0., 0., 1., 1.])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu = np.array([1, 1, 0, 0, 0])\n",
    "nu = np.array([0, 0, 0, 1, 1])\n",
    "\n",
    "C = (mu[:, np.newaxis] - nu)**2\n",
    "\n",
    "\n",
    "T, err = Sinkhorn_Knopp(mu, nu, C, 1e-1, 2)\n",
    "\n",
    "np.dot(T.T, mu)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAp8AAAH5CAYAAADHrVXSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAApvElEQVR4nO3dfXBV9Z348U+IcqNbkuq6JIFEsWJ9AqRileCPKl1axjouDOOU0Y5gV+3awo7KTruNtcXVtnG0dnEqlT6spV3L0mpBdyjFUiw+lDhWhBl83KFSeTCJdUZzMW0DJuf3B0NqNIHckHyT6Os1c0bvud+T+7lzJvrOyb03RVmWZQEAAAkMG+gBAAB4/xCfAAAkIz4BAEhGfAIAkIz4BAAgGfEJAEAy4hMAgGSOGOgBeqK9vT1eeeWVGDFiRBQVFQ30OAAAvEOWZbFnz54YNWpUDBvW/fXNIRGfr7zySlRXVw/0GAAAHMLOnTujqqqq2/uHRHyOGDEiIvY/mdLS0gGeBgCAd8rn81FdXd3Rbd0ZEvF54FftpaWl4hMAYBA71EskveEIAIBkxCcAAMmITwAAkhGfAAAkIz4BAEhGfAIAkIz4BAAgGfEJAEAy4hMAgGTEJwAAyYhPAACSEZ8AACQjPgEASOaIgR5gsGlri3jssYiGhojKyoipUyOKiwd6KgCAAgzioDmsK5+33nprFBUVxXXXXXfQdffdd1+ceuqpUVJSEuPHj481a9YczsP2m5UrI8aMiZg2LeKyy/b/c8yY/fsBAIaEQR40vY7P3//+9/G9730vJkyYcNB1GzdujEsvvTSuvPLK2Lx5c8yaNStmzZoVzzzzTG8ful+sXBlxySURu3Z13r979/79g+R8AQB0bwgETVGWZVmhB7355ptx1llnxXe/+934+te/HhMnTozFixd3uXbOnDnR0tISq1ev7tg3efLkmDhxYixdurRHj5fP56OsrCyam5ujtLS00HEPqa1t/w8E7zxPBxQVRVRVRWzfPmiuWAMAdDbAQdPTXuvVlc/58+fHRRddFNOnTz/k2vr6+netmzFjRtTX13d7TGtra+Tz+U5bf3rsse7PU0RElkXs3Ll/HQDAoDREgqbgNxytWLEinn766fj973/fo/WNjY1RXl7eaV95eXk0NjZ2e0xdXV38x3/8R6Gj9VpDQ9+uAwBIbogETUFXPnfu3BnXXntt/PSnP42SkpL+milqa2ujubm5Y9u5c2e/PVbE/jeB9eU6AIDkhkjQFHTlc9OmTfHqq6/GWWed1bGvra0tHn300bjrrruitbU1it/xGoKKiopoamrqtK+pqSkqKiq6fZxcLhe5XK6Q0Q7L1Kn7XwKxe/f+K9LvdOAlElOnJhsJAKAwQyRoCrry+Y//+I+xdevW2LJlS8d29tlnx2c+85nYsmXLu8IzIqKmpibWr1/fad+6deuipqbm8CbvQ8XFEXfeuf/fi4o633fg9uLF3mwEAAxiQyRoCorPESNGxLhx4zptf/d3fxd///d/H+PGjYuIiLlz50ZtbW3HMddee22sXbs27rjjjnjhhRfipptuiqeeeioWLFjQt8/kMM2eHXH//RGjR3feX1W1f//s2QMzFwBAjw2BoOnzv3C0Y8eOGDbsb007ZcqUWL58edx4441xww03xMknnxwPPPBAR6wOJrNnR8ycOWj/IAAAwKEN8qDp1ed8ptbfn/MJAMDh6dfP+QQAgN4QnwAAJCM+AQBIRnwCAJCM+AQAIBnxCQBAMuITAIBkxCcAAMmITwAAkhGfAAAkIz4BAEhGfAIAkIz4BAAgGfEJAEAy4hMAgGTEJwAAyYhPAACSEZ8AACQjPgEASEZ8AgCQjPgEACAZ8QkAQDLiEwCAZMQnAADJiE8AAJIRnwAAJCM+AQBIRnwCAJCM+AQAIBnxCQBAMuITAIBkxCcAAMmITwAAkhGfAAAkIz4BAEhGfAIAkIz4BAAgGfEJAEAy4hMAgGTEJwAAyYhPAACSEZ8AACRTUHzefffdMWHChCgtLY3S0tKoqamJX/3qV92uX7ZsWRQVFXXaSkpKDntoAACGpiMKWVxVVRW33nprnHzyyZFlWfz4xz+OmTNnxubNm+OMM87o8pjS0tJ48cUXO24XFRUd3sQAAAxZBcXnxRdf3On2N77xjbj77rvjiSee6DY+i4qKoqKiovcTAgDwntHr13y2tbXFihUroqWlJWpqarpd9+abb8YJJ5wQ1dXVMXPmzHj22WcP+bVbW1sjn8932gAAGPoKjs+tW7fGBz7wgcjlcnHNNdfEqlWr4vTTT+9y7SmnnBL33HNPPPjgg3HvvfdGe3t7TJkyJXbt2nXQx6irq4uysrKOrbq6utAxAQAYhIqyLMsKOWDv3r2xY8eOaG5ujvvvvz9++MMfxiOPPNJtgL7dvn374rTTTotLL700brnllm7Xtba2Rmtra8ftfD4f1dXV0dzcHKWlpYWMCwBAAvl8PsrKyg7ZawW95jMiYvjw4TF27NiIiJg0aVL8/ve/jzvvvDO+973vHfLYI488Mj7ykY/Etm3bDroul8tFLpcrdDQAAAa5w/6cz/b29k5XKQ+mra0ttm7dGpWVlYf7sAAADEEFXfmsra2NCy+8MI4//vjYs2dPLF++PDZs2BAPPfRQRETMnTs3Ro8eHXV1dRERcfPNN8fkyZNj7Nix8cYbb8Ttt98eL7/8clx11VV9/0wAABj0CorPV199NebOnRsNDQ1RVlYWEyZMiIceeig+8YlPRETEjh07Ytiwv11Mff311+Pqq6+OxsbGOOaYY2LSpEmxcePGHr0+FACA956C33A0EHr6AlYAAAZGT3vN33YHACAZ8QkAQDLiEwCAZMQnAADJiE8AAJIRnwAAJCM+AQBIRnwCAJCM+AQAIBnxCQBAMuITAIBkxCcAAMmITwAAkhGfAAAkIz4BAEhGfAIAkIz4BAAgGfEJAEAy4hMAgGTEJwAAyYhPAACSEZ8AACQjPgEASEZ8AgCQjPgEACAZ8QkAQDLiEwCAZMQnAADJiE8AAJIRnwAAJCM+AQBIRnwCAJCM+AQAIBnxCQBAMuITAIBkxCcAAMmITwAAkhGfAAAkIz4BAEhGfAIAkMwRAz3AYNPWFvHYYxENDRGVlRFTp0YUFw/0VAAABRjEQSM+32blyohrr43Ytetv+6qqIu68M2L27IGbCwCgxwZ50BT0a/e77747JkyYEKWlpVFaWho1NTXxq1/96qDH3HfffXHqqadGSUlJjB8/PtasWXNYA/eXlSsjLrmk83mKiNi9e//+lSsHZi4AgB4bAkFTUHxWVVXFrbfeGps2bYqnnnoqPv7xj8fMmTPj2Wef7XL9xo0b49JLL40rr7wyNm/eHLNmzYpZs2bFM8880yfD95W2tv0/IGTZu+87sO+66/avAwAYlIZI0BRlWVcT9tyxxx4bt99+e1x55ZXvum/OnDnR0tISq1ev7tg3efLkmDhxYixdurTbr9na2hqtra0dt/P5fFRXV0dzc3OUlpYezrhd2rAhYtq0Q6/77W8jLrigzx8eAODwDXDQ5PP5KCsrO2Sv9frd7m1tbbFixYpoaWmJmpqaLtfU19fH9OnTO+2bMWNG1NfXH/Rr19XVRVlZWcdWXV3d2zF7pKGhb9cBACQ3RIKm4PjcunVrfOADH4hcLhfXXHNNrFq1Kk4//fQu1zY2NkZ5eXmnfeXl5dHY2HjQx6itrY3m5uaObefOnYWOWZDKyr5dBwCQ3BAJmoLf7X7KKafEli1borm5Oe6///6YN29ePPLII90GaG/kcrnI5XJ99vUOZerU/W8C272765dJFBXtv3/q1GQjAQAUZogETcFXPocPHx5jx46NSZMmRV1dXZx55plx5513drm2oqIimpqaOu1ramqKioqK3k3bT4qL93/6QMT+8/J2B24vXjxoPh4LAODdhkjQHPZfOGpvb+/05qC3q6mpifXr13fat27dum5fIzqQZs+OuP/+iNGjO++vqtq/fxB8LBYAwMENgaAp6NfutbW1ceGFF8bxxx8fe/bsieXLl8eGDRvioYceioiIuXPnxujRo6Ouri4iIq699to4//zz44477oiLLrooVqxYEU899VR8//vf7/tn0gdmz46YOXPQ/kEAAIBDG+RBU1B8vvrqqzF37txoaGiIsrKymDBhQjz00EPxiU98IiIiduzYEcOG/e1i6pQpU2L58uVx4403xg033BAnn3xyPPDAAzFu3Li+fRZ9qLjYxykBAEPcIA6aw/6czxR6+rlRAAAMjH7/nE8AACiU+AQAIBnxCQBAMuITAIBkxCcAAMmITwAAkhGfAAAkIz4BAEhGfAIAkIz4BAAgGfEJAEAy4hMAgGTEJwAAyYhPAACSEZ8AACQjPgEASEZ8AgCQjPgEACAZ8QkAQDLiEwCAZMQnAADJiE8AAJIRnwAAJCM+AQBIRnwCAJCM+AQAIBnxCQBAMuITAIBkxCcAAMmITwAAkhGfAAAkIz4BAEhGfAIAkIz4BAAgGfEJAEAy4hMAgGTEJwAAyYhPAACSEZ8AACQjPgEASEZ8AgCQTEHxWVdXFx/96EdjxIgRMXLkyJg1a1a8+OKLBz1m2bJlUVRU1GkrKSk5rKEBABiaCorPRx55JObPnx9PPPFErFu3Lvbt2xef/OQno6Wl5aDHlZaWRkNDQ8f28ssvH9bQAAAMTUcUsnjt2rWdbi9btixGjhwZmzZtio997GPdHldUVBQVFRW9mxAAgPeMw3rNZ3Nzc0REHHvssQdd9+abb8YJJ5wQ1dXVMXPmzHj22WcPur61tTXy+XynDQCAoa/X8dne3h7XXXddnHfeeTFu3Lhu151yyilxzz33xIMPPhj33ntvtLe3x5QpU2LXrl3dHlNXVxdlZWUdW3V1dW/HBABgECnKsizrzYGf//zn41e/+lU8/vjjUVVV1ePj9u3bF6eddlpceumlccstt3S5prW1NVpbWztu5/P5qK6ujubm5igtLe3NuAAA9KN8Ph9lZWWH7LWCXvN5wIIFC2L16tXx6KOPFhSeERFHHnlkfOQjH4lt27Z1uyaXy0Uul+vNaAAADGIF/do9y7JYsGBBrFq1Kh5++OE48cQTC37Atra22Lp1a1RWVhZ8LAAAQ1tBVz7nz58fy5cvjwcffDBGjBgRjY2NERFRVlYWRx11VEREzJ07N0aPHh11dXUREXHzzTfH5MmTY+zYsfHGG2/E7bffHi+//HJcddVVffxUAAAY7AqKz7vvvjsiIi644IJO+3/0ox/FFVdcERERO3bsiGHD/nZB9fXXX4+rr746Ghsb45hjjolJkybFxo0b4/TTTz+8yQEAGHJ6/YajlHr6AlYAAAZGT3vN33YHACAZ8QkAQDLiEwCAZMQnAADJiE8AAJIRnwAAJCM+AQBIRnwCAJCM+AQAIBnxCQBAMuITAIBkxCcAAMmITwAAkhGfAAAkIz4BAEhGfAIAkIz4BAAgGfEJAEAy4hMAgGTEJwAAyYhPAACSEZ8AACQjPgEASEZ8AgCQjPgEACAZ8QkAQDLiEwCAZMQnAADJiE8AAJIRnwAAJCM+AQBIRnwCAJCM+AQAIBnxCQBAMuITAIBkxCcAAMmITwAAkhGfAAAkIz4BAEhGfAIAkMwRAz3AYNO2ty22fvex+PMfGuLokypj/BemRvHw4oEeCwDgPaGgK591dXXx0Y9+NEaMGBEjR46MWbNmxYsvvnjI4+6777449dRTo6SkJMaPHx9r1qzp9cD96YkvrYymo8fExOunxZS7LouJ10+LpqPHxBNfWjnQowEAvCcUFJ+PPPJIzJ8/P5544olYt25d7Nu3Lz75yU9GS0tLt8ds3LgxLr300rjyyitj8+bNMWvWrJg1a1Y888wzhz18X3riSyvjnNsviYq2XZ32V7TtjnNuv0SAAgD0gaIsy7LeHvynP/0pRo4cGY888kh87GMf63LNnDlzoqWlJVavXt2xb/LkyTFx4sRYunRpjx4nn89HWVlZNDc3R2lpaW/H7Vbb3rZoOnpMVLTt6rLG26MoGoqrouLP2/0KHgCgCz3ttcN6w1Fzc3NERBx77LHdrqmvr4/p06d32jdjxoyor6/v9pjW1tbI5/Odtv609buPxahuwjMiYlhkMbptZ2z97mP9OgcAwHtdr+Ozvb09rrvuujjvvPNi3Lhx3a5rbGyM8vLyTvvKy8ujsbGx22Pq6uqirKysY6uuru7tmD3y5z809Ok6AAC61uv4nD9/fjzzzDOxYsWKvpwnIiJqa2ujubm5Y9u5c2efP8bbHX1SZZ+uAwCga72KzwULFsTq1avjt7/9bVRVVR10bUVFRTQ1NXXa19TUFBUVFd0ek8vlorS0tNPWn8Z/YWq8UlwV7VHU5f3tURS7i6tj/Bem9uscAADvdQXFZ5ZlsWDBgli1alU8/PDDceKJJx7ymJqamli/fn2nfevWrYuamprCJu1HxcOLY8fCOyMi3hWgB27vXLjYm40AAA5TQfE5f/78uPfee2P58uUxYsSIaGxsjMbGxvjLX/7SsWbu3LlRW1vbcfvaa6+NtWvXxh133BEvvPBC3HTTTfHUU0/FggUL+u5Z9IHJt82OJ794fzQWj+60v6G4Kp784v0x+bbZAzQZAMB7R0EftVRU1PWvpX/0ox/FFVdcERERF1xwQYwZMyaWLVvWcf99990XN954Y/zxj3+Mk08+OW677bb41Kc+1eMh+/ujlt7OXzgCAChcT3vtsD7nM5WU8QkAQOGSfM4nAAAUQnwCAJCM+AQAIBnxCQBAMuITAIBkxCcAAMmITwAAkhGfAAAkIz4BAEhGfAIAkIz4BAAgGfEJAEAy4hMAgGTEJwAAyYhPAACSEZ8AACQjPgEASEZ8AgCQjPgEACAZ8QkAQDLiEwCAZMQnAADJiE8AAJIRnwAAJCM+AQBIRnwCAJCM+AQAIBnxCQBAMuITAIBkxCcAAMmITwAAkhGfAAAkIz4BAEhGfAIAkIz4BAAgGfEJAEAy4hMAgGTEJwAAyYhPAACSEZ8AACQjPgEASEZ8AgCQTMHx+eijj8bFF18co0aNiqKionjggQcOun7Dhg1RVFT0rq2xsbG3MwMAMEQVHJ8tLS1x5plnxpIlSwo67sUXX4yGhoaObeTIkYU+NAAAQ9wRhR5w4YUXxoUXXljwA40cOTI++MEPFnwcAADvHcle8zlx4sSorKyMT3ziE/G73/3uoGtbW1sjn8932gAAGPr6PT4rKytj6dKl8Ytf/CJ+8YtfRHV1dVxwwQXx9NNPd3tMXV1dlJWVdWzV1dX9PSYAAAkUZVmW9frgoqJYtWpVzJo1q6Djzj///Dj++OPjv//7v7u8v7W1NVpbWztu5/P5qK6ujubm5igtLe3tuAAA9JN8Ph9lZWWH7LWCX/PZF84555x4/PHHu70/l8tFLpdLOBEAACkMyOd8btmyJSorKwfioQEAGEAFX/l88803Y9u2bR23t2/fHlu2bIljjz02jj/++KitrY3du3fHT37yk4iIWLx4cZx44olxxhlnxF//+tf44Q9/GA8//HD8+te/7rtnAQDAkFBwfD711FMxbdq0jtsLFy6MiIh58+bFsmXLoqGhIXbs2NFx/969e+Pf/u3fYvfu3XH00UfHhAkT4je/+U2nrwEAwPvDYb3hKJWevoAVAICB0dNe87fdAQBIRnwCAJCM+AQAIBnxCQBAMuITAIBkxCcAAMmITwAAkhGfAAAkIz4BAEhGfAIAkIz4BAAgGfEJAEAy4hMAgGTEJwAAyYhPAACSEZ8AACQjPgEASEZ8AgCQjPgEACAZ8QkAQDLiEwCAZMQnAADJiE8AAJIRnwAAJCM+AQBIRnwCAJCM+AQAIBnxCQBAMuITAIBkxCcAAMmITwAAkhGfAAAkIz4BAEhGfAIAkIz4BAAgGfEJAEAy4hMAgGTEJwAAyYhPAACSEZ8AACRzxEAPMOi0tUU89lhEQ0NEZWXE1KkRxcUDPRUAQI8N5pwRn2+3cmXEtddG7Nr1t31VVRF33hkxe/bAzQUA0EODPWcK/rX7o48+GhdffHGMGjUqioqK4oEHHjjkMRs2bIizzjorcrlcjB07NpYtW9aLUfvZypURl1zS+UxFROzevX//ypUDMxcAQA8NhZwpOD5bWlrizDPPjCVLlvRo/fbt2+Oiiy6KadOmxZYtW+K6666Lq666Kh566KGCh+03bW37f0TIsnffd2DfddftXwcAMAgNlZwp+NfuF154YVx44YU9Xr906dI48cQT44477oiIiNNOOy0ef/zx+M///M+YMWNGl8e0trZGa2trx+18Pl/omIV57LF3/4jwdlkWsXPn/nUXXNC/swAA9MJQyZl+f7d7fX19TJ8+vdO+GTNmRH19fbfH1NXVRVlZWcdWXV3dv0M2NPTtOgCAxIZKzvR7fDY2NkZ5eXmnfeXl5ZHP5+Mvf/lLl8fU1tZGc3Nzx7Zz587+HbKysm/XAQAkNlRyZlC+2z2Xy0Uul0v3gFOn7n8b2O7dXb9Qoqho//1Tp6abCQCgAEMlZ/r9ymdFRUU0NTV12tfU1BSlpaVx1FFH9ffD90xx8f7PH4jYf2be7sDtxYsHzwdkAQC8w1DJmX6Pz5qamli/fn2nfevWrYuampr+fujCzJ4dcf/9EaNHd95fVbV//2D4YCwAgIMYCjlT8K/d33zzzdi2bVvH7e3bt8eWLVvi2GOPjeOPPz5qa2tj9+7d8ZOf/CQiIq655pq466674ktf+lL88z//czz88MPx85//PH75y1/23bPoK7NnR8ycOXj/JAAAwCEM9pwpyrKuXhXQvQ0bNsS0adPetX/evHmxbNmyuOKKK+KPf/xjbNiwodMx119/fTz33HNRVVUVX/3qV+OKK67o8WPm8/koKyuL5ubmKC0tLWRcAAAS6GmvFRyfA0F8AgAMbj3ttX5/zScAABwgPgEASEZ8AgCQjPgEACAZ8QkAQDLiEwCAZMQnAADJiE8AAJIRnwAAJCM+AQBIRnwCAJCM+AQAIBnxCQBAMuITAIBkxCcAAMmITwAAkhGfAAAkIz4BAEhGfAIAkIz4BAAgGfEJAEAy4hMAgGTEJwAAyYhPAACSEZ8AACQjPgEASEZ8AgCQjPgEACAZ8QkAQDLiEwCAZMQnAADJiE8AAJIRnwAAJCM+AQBIRnwCAJCM+AQAIBnxCQBAMuITAIBkxCcAAMmITwAAkhGfAAAk06v4XLJkSYwZMyZKSkri3HPPjSeffLLbtcuWLYuioqJOW0lJSa8HBgBg6Co4Pn/2s5/FwoULY9GiRfH000/HmWeeGTNmzIhXX32122NKS0ujoaGhY3v55ZcPa2gAAIamguPz29/+dlx99dXx2c9+Nk4//fRYunRpHH300XHPPfd0e0xRUVFUVFR0bOXl5Yc1NAAAQ1NB8bl3797YtGlTTJ8+/W9fYNiwmD59etTX13d73JtvvhknnHBCVFdXx8yZM+PZZ5896OO0trZGPp/vtAEAMPQVFJ+vvfZatLW1vevKZXl5eTQ2NnZ5zCmnnBL33HNPPPjgg3HvvfdGe3t7TJkyJXbt2tXt49TV1UVZWVnHVl1dXciYAAAMUv3+bveampqYO3duTJw4Mc4///xYuXJl/MM//EN873vf6/aY2traaG5u7th27tzZ32MCAJDAEYUsPu6446K4uDiampo67W9qaoqKiooefY0jjzwyPvKRj8S2bdu6XZPL5SKXyxUyGgAAQ0BBVz6HDx8ekyZNivXr13fsa29vj/Xr10dNTU2PvkZbW1ts3bo1KisrC5sUAIAhr6ArnxERCxcujHnz5sXZZ58d55xzTixevDhaWlris5/9bEREzJ07N0aPHh11dXUREXHzzTfH5MmTY+zYsfHGG2/E7bffHi+//HJcddVVfftMAAAY9AqOzzlz5sSf/vSn+NrXvhaNjY0xceLEWLt2bcebkHbs2BHDhv3tgurrr78eV199dTQ2NsYxxxwTkyZNio0bN8bpp5/ed88CAIAhoSjLsmyghziUfD4fZWVl0dzcHKWlpQM9DgAA79DTXvO33QEASEZ8AgCQjPgEACAZ8QkAQDLiEwCAZMQnAADJiE8AAJIRnwAAJCM+AQBIRnwCAJCM+AQAIBnxCQBAMuITAIBkxCcAAMmITwAAkhGfAAAkIz4BAEhGfAIAkIz4BAAgGfEJAEAy4hMAgGTEJwAAyYhPAACSEZ8AACQjPgEASEZ8AgCQjPgEACAZ8QkAQDLiEwCAZMQnAADJiE8AAJIRnwAAJCM+AQBIRnwCAJCM+AQAIBnxCQBAMuITAIBkxCcAAMmITwAAkhGfAAAkc8RADzDotLVFPPZYRENDRGVlxNSpEcXFAz0VAECPDeac6dWVzyVLlsSYMWOipKQkzj333HjyyScPuv6+++6LU089NUpKSmL8+PGxZs2aXg3b71aujBgzJmLatIjLLtv/zzFj9u8HABgCBnvOFByfP/vZz2LhwoWxaNGiePrpp+PMM8+MGTNmxKuvvtrl+o0bN8all14aV155ZWzevDlmzZoVs2bNimeeeeawh+9TK1dGXHJJxK5dnffv3r1//2A5YwAA3RgKOVOUZVlWyAHnnntufPSjH4277rorIiLa29ujuro6/vVf/zW+/OUvv2v9nDlzoqWlJVavXt2xb/LkyTFx4sRYunRpjx4zn89HWVlZNDc3R2lpaSHj9kxb2/4fCd55pg4oKoqoqorYvn3wXLMGAHibgc6ZnvZaQVc+9+7dG5s2bYrp06f/7QsMGxbTp0+P+vr6Lo+pr6/vtD4iYsaMGd2uj4hobW2NfD7faetXjz3W/ZmKiMiyiJ07968DABiEhkrOFBSfr732WrS1tUV5eXmn/eXl5dHY2NjlMY2NjQWtj4ioq6uLsrKyjq26urqQMQvX0NC36wAAEhsqOTMoP2qptrY2mpubO7adO3f27wNWVvbtOgCAxIZKzhT0UUvHHXdcFBcXR1NTU6f9TU1NUVFR0eUxFRUVBa2PiMjlcpHL5QoZ7fBMnbr/RRC7d++/Jv1OB14kMXVqupkAAAowVHKmoCufw4cPj0mTJsX69es79rW3t8f69eujpqamy2Nqamo6rY+IWLduXbfrB0RxccSdd+7/96KizvcduL14sTcbAQCD1lDJmYJ/7b5w4cL4wQ9+ED/+8Y/j+eefj89//vPR0tISn/3sZyMiYu7cuVFbW9ux/tprr421a9fGHXfcES+88ELcdNNN8dRTT8WCBQv67ln0hdmzI+6/P2L06M77q6r27589e2DmAgDooaGQMwX/haM5c+bEn/70p/ja174WjY2NMXHixFi7dm3Hm4p27NgRw4b9rWmnTJkSy5cvjxtvvDFuuOGGOPnkk+OBBx6IcePG9d2z6CuzZ0fMnDl4/yQAAMAhDPacKfhzPgdCv3/OJwAAh6VfPucTAAAOh/gEACAZ8QkAQDLiEwCAZMQnAADJiE8AAJIRnwAAJCM+AQBIRnwCAJCM+AQAIBnxCQBAMuITAIBkxCcAAMkcMdAD9ESWZRERkc/nB3gSAAC6cqDTDnRbd4ZEfO7ZsyciIqqrqwd4EgAADmbPnj1RVlbW7f1F2aHydBBob2+PV155JUaMGBFFRUX9/nj5fD6qq6tj586dUVpa2u+PR99zDoc252/ocw6HPudw6Et9DrMsiz179sSoUaNi2LDuX9k5JK58Dhs2LKqqqpI/bmlpqW+4Ic45HNqcv6HPORz6nMOhL+U5PNgVzwO84QgAgGTEJwAAyYjPLuRyuVi0aFHkcrmBHoVecg6HNudv6HMOhz7ncOgbrOdwSLzhCACA9wZXPgEASEZ8AgCQjPgEACAZ8QkAQDLiEwCAZN638blkyZIYM2ZMlJSUxLnnnhtPPvnkQdffd999ceqpp0ZJSUmMHz8+1qxZk2hSulPIOfzBD34QU6dOjWOOOSaOOeaYmD59+iHPOf2r0O/BA1asWBFFRUUxa9as/h2QQyr0HL7xxhsxf/78qKysjFwuFx/+8If9t3SAFXoOFy9eHKecckocddRRUV1dHddff3389a9/TTQtb/foo4/GxRdfHKNGjYqioqJ44IEHDnnMhg0b4qyzzopcLhdjx46NZcuW9fucXcreh1asWJENHz48u+eee7Jnn302u/rqq7MPfvCDWVNTU5frf/e732XFxcXZbbfdlj333HPZjTfemB155JHZ1q1bE0/OAYWew8suuyxbsmRJtnnz5uz555/PrrjiiqysrCzbtWtX4snJssLP3wHbt2/PRo8enU2dOjWbOXNmmmHpUqHnsLW1NTv77LOzT33qU9njjz+ebd++PduwYUO2ZcuWxJNzQKHn8Kc//WmWy+Wyn/70p9n27duzhx56KKusrMyuv/76xJOTZVm2Zs2a7Ctf+Uq2cuXKLCKyVatWHXT9Sy+9lB199NHZwoULs+eeey77zne+kxUXF2dr165NM/DbvC/j85xzzsnmz5/fcbutrS0bNWpUVldX1+X6T3/609lFF13Uad+5556b/cu//Eu/zkn3Cj2H7/TWW29lI0aMyH784x/314gcRG/O31tvvZVNmTIl++EPf5jNmzdPfA6wQs/h3XffnX3oQx/K9u7dm2pEDqHQczh//vzs4x//eKd9CxcuzM4777x+nZND60l8fulLX8rOOOOMTvvmzJmTzZgxox8n69r77tfue/fujU2bNsX06dM79g0bNiymT58e9fX1XR5TX1/faX1ExIwZM7pdT//qzTl8pz//+c+xb9++OPbYY/trTLrR2/N38803x8iRI+PKK69MMSYH0Ztz+L//+79RU1MT8+fPj/Ly8hg3blx885vfjLa2tlRj8za9OYdTpkyJTZs2dfxq/qWXXoo1a9bEpz71qSQzc3gGU8sckfwRB9hrr70WbW1tUV5e3ml/eXl5vPDCC10e09jY2OX6xsbGfpuT7vXmHL7Tv//7v8eoUaPe9Y1I/+vN+Xv88cfjv/7rv2LLli0JJuRQenMOX3rppXj44YfjM5/5TKxZsya2bdsWX/jCF2Lfvn2xaNGiFGPzNr05h5dddlm89tpr8f/+3/+LLMvirbfeimuuuSZuuOGGFCNzmLprmXw+H3/5y1/iqKOOSjbL++7KJ9x6662xYsWKWLVqVZSUlAz0OBzCnj174vLLL48f/OAHcdxxxw30OPRSe3t7jBw5Mr7//e/HpEmTYs6cOfGVr3wlli5dOtCj0UMbNmyIb37zm/Hd7343nn766Vi5cmX88pe/jFtuuWWgR2OIed9d+TzuuOOiuLg4mpqaOu1vamqKioqKLo+pqKgoaD39qzfn8IBvfetbceutt8ZvfvObmDBhQn+OSTcKPX9/+MMf4o9//GNcfPHFHfva29sjIuKII46IF198MU466aT+HZpOevM9WFlZGUceeWQUFxd37DvttNOisbEx9u7dG8OHD+/XmemsN+fwq1/9alx++eVx1VVXRUTE+PHjo6WlJT73uc/FV77ylRg2zPWsway7liktLU161TPifXjlc/jw4TFp0qRYv359x7729vZYv3591NTUdHlMTU1Np/UREevWret2Pf2rN+cwIuK2226LW265JdauXRtnn312ilHpQqHn79RTT42tW7fGli1bOrZ/+qd/imnTpsWWLVuiuro65fhE774HzzvvvNi2bVvHDw4REf/3f/8XlZWVwnMA9OYc/vnPf35XYB74YSLLsv4blj4xqFom+VucBoEVK1ZkuVwuW7ZsWfbcc89ln/vc57IPfvCDWWNjY5ZlWXb55ZdnX/7ylzvW/+53v8uOOOKI7Fvf+lb2/PPPZ4sWLfJRSwOs0HN46623ZsOHD8/uv//+rKGhoWPbs2fPQD2F97VCz987ebf7wCv0HO7YsSMbMWJEtmDBguzFF1/MVq9enY0cOTL7+te/PlBP4X2v0HO4aNGibMSIEdn//M//ZC+99FL261//OjvppJOyT3/60wP1FN7X9uzZk23evDnbvHlzFhHZt7/97Wzz5s3Zyy+/nGVZln35y1/OLr/88o71Bz5q6Ytf/GL2/PPPZ0uWLPFRS6l95zvfyY4//vhs+PDh2TnnnJM98cQTHfedf/752bx58zqt//nPf559+MMfzoYPH56dccYZ2S9/+cvEE/NOhZzDE044IYuId22LFi1KPzhZlhX+Pfh24nNwKPQcbty4MTv33HOzXC6XfehDH8q+8Y1vZG+99VbiqXm7Qs7hvn37sptuuik76aSTspKSkqy6ujr7whe+kL3++uvpByf77W9/2+X/1w6cs3nz5mXnn3/+u46ZOHFiNnz48OxDH/pQ9qMf/Sj53FmWZUVZ5lo5AABpvO9e8wkAwMARnwAAJCM+AQBIRnwCAJCM+AQAIBnxCQBAMuITAIBkxCcAAMmITwAAkhGfAAAkIz4BAEjm/wPXDjvJ686VOwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mu = np.zeros((5, 2))\n",
    "mu[:, 0] = np.array([1, 1, 0, 0, 0])\n",
    "mu[:, 1] = np.array(([0, 1, 2, 3, 4]))\n",
    "\n",
    "nu = np.zeros((5, 2))\n",
    "nu[:, 0] = np.array([0, 0, 0, 1, 1])\n",
    "nu[:, 1] = np.array(([0, 1, 2, 3, 4]))\n",
    "\n",
    "# Plot the points\n",
    "plt.figure(figsize=(8, 6))\n",
    "plt.scatter(mu[:, 0], mu[:, 1], label='mu', color='blue')\n",
    "plt.scatter(nu[:, 0], nu[:, 1], label='nu', color='red')\n",
    "\n",
    "\n",
    "C = cdist(mu, nu, 'euclidean')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Límite de iteraciones, error alcanzado: 0.18485459975258176\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[1.98896019e-01, 8.69976288e-03, 3.55406655e-07, 3.96745319e-10,\n",
       "        1.74772864e-14],\n",
       "       [1.10398041e-03, 1.91274019e-01, 4.60506594e-04, 3.05280141e-06,\n",
       "        1.34480943e-10],\n",
       "       [8.38357602e-10, 5.08349143e-05, 1.69665456e-01, 2.83950632e-07,\n",
       "        7.42816387e-11],\n",
       "       [4.19042365e-10, 2.54091842e-05, 8.48051164e-02, 1.96753710e-01,\n",
       "        3.03335539e-03],\n",
       "       [1.96279914e-14, 1.19016904e-09, 3.97228116e-06, 3.22536619e-03,\n",
       "        1.96966644e-01]])"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu_d = np.ones(mu.shape[0]) / mu.shape[0]\n",
    "nu_d = np.ones(nu.shape[0]) / nu.shape[0]\n",
    "T, err = Sinkhorn_Knopp(mu_d, nu_d, C, 1e-1, 2)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA9UAAAKTCAYAAAAXPCvaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABPaklEQVR4nO39e5xkdX3g/7+ru6cvc+me7vHLzMA03iCoCKKoMLjcsqNEkMDy8BGWzUMwEbNGNBg2FyEaNprYKCaL2RDxkojZ3VkUMmACiCEjAxjGh8rl8cMbG1REcWbUoS8zPTN9qTq/P4Yq+lLV3XW6urouz+fjUY/uOn1O1ad4HLrrNedT52SSJEkCAAAAKFvLcg8AAAAA6pWoBgAAgJRENQAAAKQkqgEAACAlUQ0AAAApiWoAAABISVQDAABASm3LPYCFyOVy8bOf/SzWrFkTmUxmuYcDAABAg0uSJPbt2xdHHnlktLSUPh5dF1H9s5/9LPr7+5d7GAAAADSZn/zkJ7Fp06aSP6+LqF6zZk1EHH4x3d3dyzwaAAAAGt3IyEj09/cXerSUuojq/JTv7u5uUQ0AAEDVzPcRZCcqAwAAgJRENQAAAKQkqgEAACAlUQ0AAAApiWoAAABISVQDAABASqIaAAAAUhLVAAAAkJKoBgAAgJRENQAAAKQkqgEAACAlUQ0AAAApiWoAAABISVQDAABASqIaAAAAUhLVAAAAkJKoBgAAgJRENQAAAKQkqgEAACAlUQ0AAAAptS33ABpFNhvx4IMRu3ZFbNwYcfrpEa2tyz0qAACAGtGg0SSqK2Dbtogrr4z46U+fX7ZpU8QnPhFx0UXLNy4AAICa0MDRtKjp39ddd11kMpl43/veN+d6t956a7zsZS+Lzs7OOOGEE+Luu+9ezNPWlG3bIt761un7RkTEM88cXr5t2/KMCwAAoCY0eDSlPlL9zW9+Mz71qU/FiSeeOOd6Dz30UFxyySUxMDAQb3nLW2Lr1q1x4YUXxiOPPBKvfOUr0z59TchmD/9jS5JERCQRcaDws8PLIn7v9yK2bGmIWQ0AAADlyWYj3vve5wPpOSsjIpMkEZlMxPveF3HBBXUbTZkkmfHqFmD//v3xmte8Jv72b/82/vzP/zxOOumkuOGGG4que/HFF8fo6GjceeedhWWnnnpqnHTSSXHTTTcV3WZsbCzGxsYK90dGRqK/vz+Gh4eju7u73OEumR07Is4+O39vNCKKj63F6eAAAIBmlcvNWjQcEaunLrjvvoizzqrSgBZmZGQkenp65u3QVLl3xRVXxHnnnRdbtmyZd92dO3fOWu+cc86JnTt3ltxmYGAgenp6Crf+/v40w1xyu3YtbL3y/9kCAACgAZSIoVlLFxpXNajs6d+33HJLPPLII/HNb35zQevv3r071q9fP23Z+vXrY/fu3SW3ufrqq+Oqq64q3M8fqa41GzdOvbcyDv97S373SArf/9M/HT6xXTGZTCYymcys7wEAAGpRfrJzkiTTvi/qwQcjzj8/8pWTee62cuZ60+OqrpQV1T/5yU/iyiuvjHvvvTc6OzuXakzR0dERHR0dS/b4lXL66YdPWPfMMxFJkokZExgikzn881/7tSRaWmbveHMR2wAAwHIqK56fM6thfu3XIvN8NBXb4HA0lToKWQfKmv798MMPx89//vN4zWteE21tbdHW1hb3339//PVf/3W0tbVFNpudtc2GDRtiz54905bt2bMnNmzYsLiR14DW1sNngI84vC9Mlb9/ww0RbW2ZaGlpiZaWlmhtbS38t2ttbS0snxnNSZJELpeLXC4X2Ww2JicnY3JyMrLZbGF5io/DAwAATJOP5nx75Psj/32+Pab2Rz6a842T75z89/nGybS1LSya6vQkZRFlRvV//I//MR5//PF47LHHCrfXvva18Zu/+Zvx2GOPRWuR/xCbN2+O7du3T1t27733xubNmxc38hpx0UURt90WcdRR05dv2nR4+VyXXMvvhGIbAABYaksaz3PNrF1MNNWBsqZ/r1mzZtZlsFatWhXr1q0rLL/00kvjqKOOioGBgYiIuPLKK+PMM8+Mv/zLv4zzzjsvbrnllvjWt74Vn/70pyv0EpbfRRcdPgP8gw8e/nz9xo2HZy+k/ceWUjvlzGkXU3f2UtPKTSMHAIDmUpFp21OWVUSlo6mGpL5OdSlPP/10tEy5htRpp50WW7dujQ984ANxzTXXxLHHHht33HFH3V+jeqbW1qU/A7zYBgAA8moynudSjWhaBqmuU11tC70+GNPNFduliG0AAKgtaePZ+/rFWWiHVvxINbVjriPb+a8z/6d0ZBsAAJaHeK5PoroJlZrmIbYBAGDpiefGIqopENsAAFA54rk5iGrmJbYBAKA08dzcRDWpiW0AAJqJEwFTjKim4pYytos9LgAAVJJ4phyimqqpdGznv4ptAADSEM9Ugqhm2aWN7alfiz2W2AYAIEI8s7RENTVLbAMAUA7xzHIQ1dQdsQ0A0NzEM7VEVNMwxDYAQGMRz9QDUU3DE9sAALWt3HgudnUY78lYLqKapiW2AQCqazHx7D0WtUpUwwxz/cIu9YdAbAMAPE8800xENZSh1NQisQ0ANCPxDKIaKkJsAwCNTDxDaaIalpDYBgDqRbH3I+IZ5ieqYRmIbQBgucz1XqMU7zmgNFENNWQpYnuuxwUAGpd4huoQ1VAHKhXbMx9LbANA/RPPsLxENdSxcmO72P1ijyW2AaD2iGeoTaIaGpDYBoD6JZ6hvohqaCJiGwBqh3iGxiCqAbENAEtIPENjE9VASWIbABYubTwXu2oHUD9ENVA2sQ1AM6tEPPt7B41DVAMVM1ds578Wu+yX2AagFolnYCFENbDkSk1nE9sA1ALxDCyGqAaWjdgGoJrEM7AURDVQc8Q2AIshnoFqEtVA3RDbAEw11wkyS/H7H6g0UQ3UvUrHdv6rS5sA1AbxDNQyUQ00rLSxPfVrsccS2wBLQzwD9UhUA01HbAMsL/EMNBJRDfAcsQ1QWeIZaAaiGmAeYhtgbuIZaGaiGiAlsQ00G/EMMJuoBqgwsQ3Uu3Lj2e8qoJmJaoAqEdtArRHPAIsnqgGW2VxvTEu94RXbQDnEM8DSEdUANazU5w/FNlCMeAaoPlENUIfENjQ38QxQO0Q1QAMR29A4iv1/Kp4Bao+oBmgCYhtq11z/D5bi/0WA2iGqAZpYpWJ75mO5Fi3MJp4BGpOoBmCWcmO72P1ijyW2aQbiGaC5iGoAFkxsw/PEMwARohqAChDbNDLxDMBcRDUAS0ZsU0/EMwBpiGoAqk5ss5zEMwCVJKoBqBlim0pKG8/2GwDKIaoBqHkLje2ZX8V2cxDPACwnUQ1A3ZortvNfxXbjEM8A1CJRDUDDKfWZV7FdH8QzAPVEVAPQNMR2bRHPADQCUQ1A0xPbS2uuz76X4r8jAPVCVANACWK7PHOdpb2UZvjvAkBjE9UAUKZKx3b+a71c/1g8A8DzRDUAVEja2J76tdhjLVdsi2cAmF9LOSt/8pOfjBNPPDG6u7uju7s7Nm/eHF/+8pdLrn/zzTcX/qDmb52dnYseNADUk/zfwJaWlmhtbY3W1tZoa2srfN/S0jIrQPMRm8vlIpvNRjabjcnJycL3uVxuwaE7n/zz5J9rcnKy8Fz55cWiv6WlpfCa2traCq8pv1xQA9AMyjpSvWnTprjuuuvi2GOPjSRJ4vOf/3xccMEF8eijj8bxxx9fdJvu7u544oknCvf9gQWAw6p9ZNuRZwCovLKi+vzzz592/y/+4i/ik5/8ZHz9618vGdWZTCY2bNhQ1qDGxsZibGyscH9kZKSs7QGgni02tmcGc5Ik804jF88AkE5Z07+nymazccstt8To6Ghs3ry55Hr79++PF77whdHf3x8XXHBBfOc735n3sQcGBqKnp6dw6+/vTztMAGgYM6eRz5xmnZ++nZ+yPXXq9syf5aN75mOZtg0A5ckkZX4Y6/HHH4/NmzfHoUOHYvXq1bF169Y499xzi667c+fO+Pd///c48cQTY3h4OD7+8Y/HAw88EN/5zndi06ZNJZ+j2JHq/v7+GB4eju7u7nKGCwB1r9xp21OjOJPJzDpqvZBt6+ls5ACwFEZGRqKnp2feDi07qsfHx+Ppp5+O4eHhuO222+Kzn/1s3H///fGKV7xi3m0nJibi5S9/eVxyySXx4Q9/eMHPudAXAwD1LG08lxvAC73cVyWeCwDq1UI7tOxLarW3t8cxxxwTEREnn3xyfPOb34xPfOIT8alPfWrebVesWBGvfvWr48knnyz3aQGgYVQrnhfyeAsdWy1f+gsAltOir1Ody+WmTdWeSzabjccff7zkdHEAaCTLHc9plDpJmdgGgOLKiuqrr7463vzmN8fRRx8d+/bti61bt8aOHTviK1/5SkREXHrppXHUUUfFwMBARER86EMfilNPPTWOOeaYGBoaiuuvvz5+/OMfx+WXX175VwIAy6Qe47lcYhsAiisrqn/+85/HpZdeGrt27Yqenp448cQT4ytf+Uq88Y1vjIiIp59+Olpanj+h+ODgYLzzne+M3bt3R29vb5x88snx0EMPLejz1wBQS+YKx1KaISDFNgDNruwTlS0HJyoDoFrE89JqhqP6ADSGJTtRGQA0AvG8PCp1ZHvmY5V6XABYaqIagIYmnutDubFd7H6xxxLbACw1UQ1AQxDPjUlsA1DrRDUAdUU8EyG2AagdohqAmiSeSUNsA1BtohqAZSWeqQaxDcBSEdUAVEXaeBYvLCWxDcBiiWoAKko80wjmiu3812KX/RLbAM1HVAOQinimGZX6uIHYBmheohqAOYlnmJ/YBmheohqAiBDPsBTENkDjE9UATabUm/i5eDMPlSW2ARqHqAZoUHOdvbgUb85heS1lbBd7XAAWT1QD1DnxDI2v0rGd/yq2ARZPVAPUCfEMzJQ2tqd+LfZYYhtg4UQ1QI0Rz8BiiW2A6hHVAMtEPAPVJrYBKk9UAywx8QzUOrENkJ6oBqiQcuO52Fl5vfEEaonYBpifqAYo02Li2RtJoBGIbYDniWqAEsQzQHnm+t1X6neq2AbqnagGmp54Blh6pT7iIraBeieqgaYhngFqj9gG6p2oBhqOeAaof2IbqBeiGqhLc72hKsUbK4D6J7aBWiOqgZomngFYiKWI7bkeFyBPVAM1QTwDsBQqFdszH0tsA3miGqgq8QxALSg3tovdL/ZYYhuaj6gGloR4BqAeiW2gXKIaWBTxDEAzENtAKaIaWBDxDACziW1AVAPTpI3nYmdLBYBmJbaheYhqaFKViGd/1AGgPGIbGo+ohgYnngGg9s0V2/mvxS77JbZh+YlqaBDiGQAaT6mPVYltqB2iGuqMeAYAxDbUDlENNUo8AwDlEttQfaIaltlcJyYpxR85AKAcYhuWjqiGKhHPAECtWcrYLva40IhENVSYeAYA6l2lYzv/VWzTiEQ1pCSeAYBmkza2p34t9lhim3omqmEe4hkAYG5im2YmquE54hkAoLLENs1AVNN0xDMAwPIS2zQSUU3DKjee/UIGAFheYpt6JKqpe+IZAKCxiW1qmaimbohnAACmmus9Xqn3jmKbShPV1BzxDADAYpU6D47YptJENctGPAMAUG1im0oT1SypYr+MxDMAALVGbJOWqKYi5vpFU4pfOAAA1DqxzXxENWURzwAAsDSxPdfjUrtENUWJZwAAKF+lYnvmY4nt2iWqm5x4BgCApVdubBe7X+yxxPbyE9VNQjwDAEDtEdv1T1Q3GPEMAAD1T2zXD1Fdp8QzAAA0H7Fde0R1hWSzEQ8+GLFrV8TGjRGnnx7R2rr4x00bz/7nAACA5lEXsb1U0bTMRHUFbNsWceWVET/96fPLNm2K+MQnIi66aGGPIZ4BAIBKW2hsz/xa8diuRDTVqJZyVv7kJz8ZJ554YnR3d0d3d3ds3rw5vvzlL8+5za233hove9nLorOzM0444YS4++67FzXgWrNtW8Rb3zp934iIeOaZw8u3bZu+PL+D5nK5yGazkc1mY3JysvB9LpebtRNnMploaWmJlpaWaG1tjba2tmhra4vW1tbCckENAAAs1MzGyHdG/vt8Y0ztjHzH5FtmcnKy0DL55UUPDJYbTXUmk8x3OHSKf/7nf47W1tY49thjI0mS+PznPx/XX399PProo3H88cfPWv+hhx6KM844IwYGBuItb3lLbN26NT760Y/GI488Eq985SsXPMiRkZHo6emJ4eHh6O7uXvB2Sy2bjXjRi/L7RhIRB6b89PB/1iOPTOLb306itdWRZwAAoD6lnlmby0W84hWR2bUrptbNyojD9zOZw0esf/SjmpsKvtAOLSuqi+nr64vrr78+3vGOd8z62cUXXxyjo6Nx5513FpadeuqpcdJJJ8VNN91U8jHHxsZibGyscH9kZCT6+/trLqp37Ig4++z8vdGIKD62TObwrRTxDAAA1KN5czJJDt/yd5/7+q2IeM3U9e67L+Kssyo7uEVaaFSXNf17qmw2G7fcckuMjo7G5s2bi66zc+fO2LJly7Rl55xzTuzcuXPOxx4YGIienp7Crb+/P+0wl9SuXeWtnz/6PPMGAABQj0o1Tv6IdpIkkUQUbnm/nPlA5cZVDSn7RGWPP/54bN68OQ4dOhSrV6+O22+/PV7xilcUXXf37t2xfv36acvWr18fu3fvnvM5rr766rjqqqsK9/NHqmvNxo1T762MiOHnvp8eynfdFXHGGVUaFAAAwBJJkiRGRkbi2WefjcHBwRgcHIxnn302hoaGYnJycvYGP/hBtN50U/RERF9ErH3u69Ez15seV3Wl7Kg+7rjj4rHHHovh4eG47bbb4rLLLov777+/ZFin0dHRER0dHRV7vKVy+umHp/8/80xEkmQiYvW0n+c/HvCmN9XcxwMAAABKSpIkhoeH49lnn429e/cWIvrZZ5+NbDY7a/3W1tZob2+P3t7e6Ovre/7W0xM9//zP0fKzn02bBl6Qj6bTT6/Cq1oaZUd1e3t7HHPMMRERcfLJJ8c3v/nN+MQnPhGf+tSnZq27YcOG2LNnz7Rle/bsiQ0bNqQcbm1pbT18Bvi3vvXwvjB1H8nP6r7hBkENAADUplwuV4jnqbfBwcGi8RxxOKD7+vqit7c31q1bVwjo7u7uaGkp8gnjv/7rho6mRV+nOpfLTTup2FSbN2+O7du3x/ve977CsnvvvbfkZ7Dr0UUXRdx2W/FLrt1wQ91fcg0AAGgAuVwuhoaGYnBwsHDkOT9te754zt/yAd3d3V3euaEaPJrKiuqrr7463vzmN8fRRx8d+/bti61bt8aOHTviK1/5SkREXHrppXHUUUfFwMBARERceeWVceaZZ8Zf/uVfxnnnnRe33HJLfOtb34pPf/rTlX8ly+iiiyIuuCDiwQcPf75+48bDsxfq+B9bAACAOpTNZmdN287Hcy6XK7pNW1vb9CnbU448V+zEyg0cTWVF9c9//vO49NJLY9euXdHT0xMnnnhifOUrX4k3vvGNERHx9NNPTzvcf9ppp8XWrVvjAx/4QFxzzTVx7LHHxh133FHWNarrRWtrzZ0BHgAAaFDZbDaGhoYK0bx3794YHBycM55XrFhRNJ7XrFlTnasSNWg0Lfo61dWw0OuDAQAANJJsNls4QdjU29DQUMlrRM+M5/y07dWrV7ukbxkW2qGL/kw1AAAAi5PNZgsnCJs6bXt4eLhkPLe3t8866rxu3bpYvXp10fVZGqIaAACgSiYnJ6cdec4H9MjIyJzxPPUs2/mbeK4NohoAAKDCJiYmik7bHh4eLrlNR0fHtOna+duqVauqOHLKJaoBAABSmpiYmBXO+SPPpXR2dhadtr1y5coqjpxKEdUAAADzGB8fLzpte9++fSW36erqKnq2bfHcWEQ1AADAc8bHx2cddd67d2/s37+/5Db5eJ45bburq6uKI2e5iGoAAKDpjI2NFZ22PVc8r1y5ctp07d7eXvGMqAYAABpXPp6nXqbq2WefjdHR0ZLbrFq1qui07c7OziqOnHohqgEAgLp36NChotO2Dxw4UHKb1atXR29v76xp2x0dHVUcOfVOVAMAAHXj4MGDRadtzxfPU6dt579vb2+v4shpVKIaAACoOfl4njlt++DBgyW3WbNmTdFp2+KZpSSqAQCAZXPgwIFp07Xzl62aK567u7tnhXNvb694ZlmIagAAYMmNjo4WnbZ96NChktv09PQUzrCdn7bd29sbK1asqOLIYW6iGgAAqJj9+/cXjjZPnbo9NjZWcpuenp5pR53XrVsXa9euFc/UBVENAACUbf/+/dOmbee/Hx8fL7p+JpOZFc/5I89tbbKE+mXvBQAAStq3b1/RadtzxfPatWuLxnNra2uVRw9LT1QDAECTS5KkZDxPTEwU3aZYPOenbYtnmomoBgCAJlEsnvNn3C4Vzy0tLUWPPItnOExUAwBAg0mSJEZGRorG8+TkZNFtWlpaCmfannrr6ekRzzAHUQ0AAHUqSZIYHh6eNWV7rnhubW2NtWvXFi5RNTWeW1paqvwKoP6JagAAqHG5XC5GRkamnWU7H8/ZbLboNq2trSWPPItnqBxRDQAANSKXy8068rx3794YGhqaM55nhnNfX190d3eLZ6gCUQ0AAFWWy+ViaGio6LTtXC5XdJu2trbo7e2dNW27u7s7MplMlV8BkCeqAQBgiWSz2RgeHp41bXtoaGjOeJ56iar8NZ7FM9QmUQ0AAIuUzWYLR57zAT04ODhnPK9YsaLotO01a9aIZ6gjohoAABYom83G4ODgrGnbQ0NDkSRJ0W3a29sLR5unTt1evXq1eIYGIKoBAGCGycnJovE8PDw8bzxPnbadj2egcYlqAACa1tR4nvq555GRkTnjeWo0r1u3Lnp7e8UzNClRDQBAw5uYmJh25Hnv3r0xODgYw8PDJbfp6OiYdabtvr6+WLVqVRVHDtQ6UQ0AQMOYmJiYNWU7f+S5lM7OzsLR5qkRvXLlyiqOHKhXohoAgLozPj5eOMP21Knb+/btK7lNV1fXtCPO+YDu6uqq4siBRiOqAQCoWfl4njpt+9lnn439+/eX3Karq6votG3xDCwFUQ0AwLIbGxsrOm17rnheuXLlrHju7e0Vz0BViWoAAKpmbGxs2lm287fR0dGS26xatWrWUee+vr7o7Oys4sgBihPVAABU3KFDh2aF8969e+PAgQMlt1m9enXReO7o6KjiyAHKI6oBAEjt4MGDReP54MGDJbdZvXp10c88t7e3V3HkAJUhqgEAmNfBgweLTtueK57XrFkz7Szbvb294hloOKIaAICCAwcOTDvL9uDgYOzduzcOHTpUcpvu7u5ZR517e3vFM9AURDUAQBMaHR0terbtueK5p6cnent7p03d7u3tjRUrVlRx5AC1RVQDADSw/fv3F43nsbGxktv09PRMm7adj+e2Nm8dAWbymxEAoAHk43nm557Hx8eLrp/JZArxnA/o3t5e8QxQJr8xAQDqyL59+6adZXtwcHDeeF67dm3Rzzy3trZWefQAjUdUAwDUmCRJpsXz1NvExETRbVpaWmLt2rWFM2znp26vXbtWPAMsIVENALBMkiSJkZGRwhm28+E8ODg4bzxPnbbd19cXPT094hlgGYhqAIAllo/nqdO28/E8OTlZdJuWlpbCUeept7Vr10ZLS0uVXwEApYhqAIAKSZIkhoeHi07bzmazRbdpbW0tGs89PT3iGaAOiGoAgDLlcrmi8Tw4OFhWPK9bty66u7vFM0AdE9UAACUUi+e9e/fG0NDQnPE886hzPp4zmUyVXwEAS01UAwBNL5fLxdDQUNEjz7lcrug2bW1ts+K5r69PPAM0GVENADSNbDZbNJ6HhoZKxvOKFSuit7e3cJbt/G3NmjXiGQBRDQA0nnw8T71MVT6ekyQpus2KFSuKTttevXq1eAagJFENANStbDYbg4ODRY88l4rn9vb2otO216xZU+XRA9AIRDUAUPMmJyeLxvPw8PC88Txz2vbq1aurPHoAGpmoBgBqRj6eZ07bHh4eLrlNR0fHtOnafX190dvbK54BqApRDQBU3cTEROHIcz6gBwcH54znzs7OotO2V61aVcWRA8B0ohoAWDITExOzpmw/++yzMTIyUnKbzs7OWVO2+/r6YuXKlVUcOQAsjKgGABZtfHy8aDzv27ev5DZdXV2zpm339fVFV1dXFUcOAIsjqgGABZsaz1M/97x///6S26xcubLotG3xDEAjENUAwCxjY2PTjjjv3bs3BgcH54znVatWFY3nzs7OKo4cAKqrrKgeGBiIbdu2xfe///3o6uqK0047LT760Y/GcccdV3Kbm2++OX7rt35r2rKOjo44dOhQuhEDABVz6NChotO2R0dHS24zNZ6nnm1bPAPQjMqK6vvvvz+uuOKKeN3rXheTk5NxzTXXxJve9Kb47ne/O+eZN7u7u+OJJ54o3M9kMulHDACU7dChQ4WjzVOnbR84cKDkNqtXr552xHndunXR29sbHR0dVRw5ANS2sqL6nnvumXb/5ptvjiOOOCIefvjhOOOMM0pul8lkYsOGDQt+nrGxsRgbGyvcn+sMoQDA8w4ePDhr2vazzz4bBw8eLLnNmjVrik7bbm9vr+LIAaA+Leoz1flrSfb19c253v79++OFL3xh5HK5eM1rXhMf+chH4vjjjy+5/sDAQPzZn/3ZYoYGAA3twIEDRadtzxfPMy9V1dvbK54BYBEySZIkaTbM5XLx67/+6zE0NBRf+9rXSq63c+fO+Pd///c48cQTY3h4OD7+8Y/HAw88EN/5zndi06ZNRbcpdqS6v78/hoeHo7u7O81wAaAuHThwYNp07fxtrnOTdHd3F522vWLFiiqOHADq28jISPT09Mzboamj+nd/93fjy1/+cnzta18rGcfFTExMxMtf/vK45JJL4sMf/vCCtlnoiwGAejU6OjornPfu3TvtH5ln6unpmTVlWzwDQGUstENTTf9+z3veE3feeWc88MADZQV1RMSKFSvi1a9+dTz55JNpnhoA6tr+/fuLTtsuFc+ZTCa6u7uLTttua3NlTABYbmX9NU6SJN773vfG7bffHjt27IgXv/jFZT9hNpuNxx9/PM4999yytwWAerF///6i07bHx8eLrp/JZApHnqcG9Nq1a8UzANSwsv5KX3HFFbF169b40pe+FGvWrIndu3dHxOHpZ11dXRERcemll8ZRRx0VAwMDERHxoQ99KE499dQ45phjYmhoKK6//vr48Y9/HJdffnmFXwoAVN++ffumnWV7cHBw3nheu3Zt0Wnbra2tVR49ALBYZUX1Jz/5yYiIOOuss6Yt/9znPhdvf/vbIyLi6aefjpaWlsLPBgcH453vfGfs3r07ent74+STT46HHnooXvGKVyxu5ABQJUmSFOJ55m1iYqLoNi0tLUXjee3ateIZABpI6hOVVZMTlQFQDUmSxMjISNF4npycLLrN1HieOm27p6dHPANAHVvSE5UBQD2bGs9TP/c8ODhYMp5bW1unHXnOX6Zq7dq102ZoAQDNRVQD0LCSJInh4eFpl6jKf+Y5m80W3aa1tTV6e3tnTdvu6ekRzwDALKIagLqXy+WmxfPUI89zxXP+BGFTp213d3eLZwBgwUQ1AHUjl8vF0NBQDA4OTpu2PTQ0VDKe29raph15zgd0d3d3ZDKZKr8CAKDRiGoAak42m501bTsfz7lcrug2bW1ts6Zsi2cAYKmJagCWTTabjaGhoVnTtueK5xUrVhSN5zVr1ohnAKDqRDUASy6bzRZOEDYznktd2bG9vb3wmeep07ZXr14tngGAmiGqAaiYmfGcn7Y9PDw8bzxPva1bty5Wr15d5dEDAJRPVANQtsnJyaJHnueL56ln2c7fxDMAUM9ENQAlTU5OzgrnfDyX0tHRUTSeV61aVcWRAwBUh6gGICYmJmZdpurZZ5+NkZGRktt0dnYWnba9cuXKKo4cAGB5iWqAJjI+Pl502vZc8dzV1VX0bNviGQBAVAM0pPHx8aLTtvft21dym66urqLTtru6uqo4cgCA+iKqAepYPp5nTtvev39/yW1Wrlw57RJV+UtWiWcAgPKJaoA6MDY2NusyVYODg3PG86pVq4pO2+7s7KziyAEAGpuoBqghhw4dKjpte3R0tOQ2q1evLhrPHR0dVRw5AEBzEtUAy+DgwYNF4/nAgQMlt8nH88zPPbe3t1dx5AAATCWqAZbQ1Hie+rnngwcPltxmzZo10y5Rlf/Ms3gGAKg9ohqgAg4cOFD0yPNc8dzd3T1rynZvb694BgCoI6IaoAyjo6NF4/nQoUMlt+np6SkcbZ56xu0VK1ZUceQAACwFUQ1QxNR4njpte2xsrOQ2PT090446r1u3LtauXSueAQAamKgGmtr+/funHXHOB/T4+HjR9TOZzKx4zh95bmvzKxUAoNl4Bwg0hX379hWdtj1XPK9du7ZoPLe2tlZ59AAA1CpRDTSMJEli//7906Zr528TExNFtykWz/lp2+IZAID5iGqg7iRJMuvI8969e2NwcLBkPLe0tBQ98iyeAQBYDFEN1KwkSWJkZKTotO3Jycmi27S0tBTOtD311tPTI54BAKg4UQ0su3w8z5y2PTg4WDKeW1tbY+3atYVLVE2N55aWliq/AgAAmpWoBqomSZIYHh6edpbtwcHBePbZZyObzRbdprW1teSRZ/EMAMByE9VAxeVyuUI8zzzyPFc8zwznvr6+6O7uFs8AANQsUQ2klsvlYmhoKAYHB6dN3R4aGioZz21tbdHb2ztr2nZ3d3dkMpkqvwIAAFgcUQ3MK5vNzpq2nY/nXC5XdJu2trZpl6jKX+NZPAMA0EhENVCQzWZjaGho1mWq5ornFStWFJ22vWbNGvEMAEDDE9XQhLLZbOEEYVNvQ0NDkSRJ0W3a29sLR5unTt1evXq1eAYAoGmJamhg2Wy2cIKwqdO2h4eH543nqdO28/EMAABMJ6qhAUxOTk478pwP6JGRkTnjeWo0r1u3Lnp7e8UzAACUQVRDHZmYmCg6bXt4eLjkNh0dHbPOtN3X1xerVq2q4sgBAKAxiWqoQRMTE7PCOX/kuZTOzs7C0eapEb1y5coqjhwAAJqLqIZlND4+XnTa9r59+0pu09XVNe2Icz6gu7q6qjhyAAAgQlRDVYyPj8866rx3797Yv39/yW1WrlxZ9FJV4hkAAGqHqIYKGhsbKzpte754LvaZ587OziqOHAAASENUQwr5eJ56mapnn302RkdHS26zatWqWdO2e3t7xTMAANQxUQ1zOHToUNFp2wcOHCi5zerVq4tO2+7o6KjiyAEAgGoQ1RARBw8eLDpte754LjZtu729vYojBwAAlpOopqnk43nmtO2DBw+W3GbNmjXTzrLd29srngEAgIgQ1TSoAwcOTJuunb9s1Vzx3N3dXXTa9ooVK6o4cgAAoJ6Iaura6Oho0Wnbhw4dKrlNT0/PrHDu7e0VzwAAQNlENXVh//79haPNU6duj42NldwmH89TP/fc29sbbW12ewAAoDLUBTVl//7906Zt578fHx8vun4mk5l25Dl/mSrxDAAAVIPqYFns27ev6LTtueJ57dq1Radtt7a2Vnn0AAAAh4lqlkySJCXjeWJioug2LS0tsXbt2ujt7Z02bXvt2rXiGQAAqDmimkUrFs/5M27PF89Tp2339fVFT0+PeAYAAOqGqGbBkiSJkZGRovE8OTlZdJuWlpbCdZ2n3tauXRstLS1VfgUAAACVJaqZJUmSGB4enjVle654bm1tLRrPPT094hkAAGhYorqJ5XK5GBkZmXaW7Xw8Z7PZotvMjOf8tO3u7m7xDAAANB1R3QRyudysI8979+6NoaGhOeN55lHndevWRXd3d2QymSq/AgAAgNokqhtILpeLoaGhotO2c7lc0W3a2tpmxXP+yLN4BgAAmJuorkPZbDaGh4dnTdseGhoqGc8rVqyYdZmqvr6+WLNmjXgGAABISVRXSHY8G4//7YNx4Ae7YuVLN8YJ7z49WtsXd2mobDZbOPKcD+jBwcF547nYtO3Vq1eLZwAAgAoT1RXw9T/aFkf/1ZVxUvanhWU/+4NN8fRVn4hTP3bRvNtns9kYHBycNW17aGgokiQpuk17e3vRadtr1qyp2OsCAABgbmVF9cDAQGzbti2+//3vR1dXV5x22mnx0Y9+NI477rg5t7v11lvjgx/8YDz11FNx7LHHxkc/+tE499xzFzXwWvH1P9oWr7/+rRExPX43ZJ+JDde/Nb4etxXCenJyMgYHB2NwcHDa1O3h4eF543nmtO3Vq1cv9UsDAABgHmVF9f333x9XXHFFvO51r4vJycm45ppr4k1velN897vfjVWrVhXd5qGHHopLLrkkBgYG4i1veUts3bo1LrzwwnjkkUfila98ZUVexHLJjmfj6L+6MiKSmHoxqcmIGIwkfhkRP//L340vbW6N4X3DMTIyUjKeOzo6Zl2mqq+vr+R/VwAAAJZfJilVeQvwi1/8Io444oi4//7744wzzii6zsUXXxyjo6Nx5513FpadeuqpcdJJJ8VNN91UdJuxsbEYGxsr3B8ZGYn+/v4YHh6O7u7utMOtuMdu2BEn/f7ZERExGhHbI2IwIkZmrPfT//Te2HT2sRER0dnZWXTatngGAACoHSMjI9HT0zNvhy7qM9XDw8MREdHX11dynZ07d8ZVV101bdk555wTd9xxR8ltBgYG4s/+7M8WM7SqOPCDXYXvOyLiqee+z0REZ0Ssi4i+iOjpOCp+9T/9p+jr64uVK1dWe5gAAAAskdRRncvl4n3ve1+84Q1vmHMa9+7du2P9+vXTlq1fvz52795dcpurr756Wojnj1TXmpUv3Vj4vi0itkREdxyO6a4p6z12yimxadOm6g4OAACAJZc6qq+44or49re/HV/72tcqOZ6IOPz54o6Ojoo/bqWd8O7T42d/sCk2ZJ+Jlkji5c8tz1+4KheZ2NW6KU549+nLNUQAAACWUMv8q8z2nve8J+68886477775j0Cu2HDhtizZ8+0ZXv27IkNGzakeeqa0treGk9f9YmIOBzQU+Xv/+SqGxZ9vWoAAABqU1lRnSRJvOc974nbb789vvrVr8aLX/ziebfZvHlzbN++fdqye++9NzZv3lzeSGvUqR+7KL7xh7fF7tajpi3f1bopvvGHty3oOtUAAADUp7LO/v3ud787tm7dGl/60pemXZu6p6cnuroOf4r40ksvjaOOOioGBgYi4vAltc4888y47rrr4rzzzotbbrklPvKRj5R1Sa2FnnVtOWXHs/H/u/GBOPCDXbHqmCPjhHef7gg1AABAnVpoh5YV1ZlMpujyz33uc/H2t789IiLOOuuseNGLXhQ333xz4ee33nprfOADH4innnoqjj322PjYxz4W55577kKfti6iOiIK16Au9d8JAACA+rAkUb1cRDUAAADVtNAOTXWiMgAAAEBUAwAAQGqiGgAAAFIS1QAAAJCSqAYAAICURDUAAACkJKoBAAAgJVENAAAAKYlqAAAASElUAwAAQEqiGgAAAFIS1QAAAJCSqAYAAICURDUAAACkJKoBAAAgJVENAAAAKYlqAAAASElUAwAAQEqiGgAAAFIS1QAAAJCSqAYAAICURDUAAACkJKoBAAAgJVENAAAAKYlqAAAASElUAwAAQEqiGgAAAFIS1QAAAJCSqAYAAICURDUAAACkJKoBAAAgJVENAAAAKYlqAAAASElUAwAAQEqiGgAAAFIS1QAAAJCSqAYAAICURDUAAACkJKoBAAAgJVENAAAAKYlqAAAASElUAwAAQEqiGgAAAFIS1QAAAJCSqAYAAICURDUAAACkJKoBAAAgJVENAAAAKYlqAAAASElUAwAAQEqiGgAAAFIS1QAAAJCSqAYAAICURDUAAACkJKoBAAAgJVENAAAAKYlqAAAASElUAwAAQEqiGgAAAFIqO6ofeOCBOP/88+PII4+MTCYTd9xxx5zr79ixIzKZzKzb7t27044ZAAAAakLZUT06OhqvetWr4sYbbyxruyeeeCJ27dpVuB1xxBHlPjUAAADUlLZyN3jzm98cb37zm8t+oiOOOCLWrl27oHXHxsZibGyscH9kZKTs5wMAAIClVrXPVJ900kmxcePGeOMb3xj/9m//Nue6AwMD0dPTU7j19/dXaZQAAACwcEse1Rs3boybbrop/vEf/zH+8R//Mfr7++Oss86KRx55pOQ2V199dQwPDxduP/nJT5Z6mAAAAFC2sqd/l+u4446L4447rnD/tNNOix/84AfxP/7H/4j/9b/+V9FtOjo6oqOjY6mHBgAAAIuyLJfUev3rXx9PPvnkcjw1AAAAVMyyRPVjjz0WGzduXI6nBgAAgIope/r3/v37px1l/tGPfhSPPfZY9PX1xdFHHx1XX311PPPMM/EP//APERFxww03xItf/OI4/vjj49ChQ/HZz342vvrVr8a//Mu/VO5VAAAAwDIoO6q/9a1vxdlnn124f9VVV0VExGWXXRY333xz7Nq1K55++unCz8fHx+O//bf/Fs8880ysXLkyTjzxxPjXf/3XaY8BAAAA9SiTJEmy3IOYz8jISPT09MTw8HB0d3cv93BKyv+nzGQyyzwSAAAAFmOhHbosn6kGAACARiCqAQAAICVRDQAAACmJagAAAEhJVAMAAEBKohoAAABSEtUAAACQkqgGAACAlEQ1AAAApCSqAQAAICVRDQAAACmJagAAAEhJVAMAAEBKohoAAABSEtUAAACQkqgGAACAlEQ1AAAApCSqAQAAICVRDQAAACmJagAAAEhJVAMAAEBKohoAAABSEtUAAACQkqgGAACAlEQ1AAAApCSqAQAAICVRDQAAACmJagAAAEhJVAMAAEBKohoAAABSEtUAAACQkqgGAACAlEQ1AAAApCSqAQAAICVRDQAAACmJagAAAEhJVAMAAEBKohoAAABSEtUAAACQkqgGAACAlEQ1AAAApCSqAQAAICVRDQAAACmJagAAAEhJVAMAAEBKohoAAABSEtUAAACQkqgGAACAlEQ1AAAApCSqAQAAICVRDQAAACmJagAAAEhJVAMAAEBKohoAAABSEtUAAACQkqgGAACAlEQ1AAAApNS23ANoGNlsxAMPROzaFXHkkRGnnx7R2rrcowIAAKgJ2WzEgw8eTqaNGxsnmUR1JWzbFnHllRE//enzyzZtivjEJyIuumj5xgUAAFADGjmZyp7+/cADD8T5558fRx55ZGQymbjjjjvm3WbHjh3xmte8Jjo6OuKYY46Jm2++OcVQa9S2bRFvfev0vSMi4plnDi/ftm15xgUAAFADGj2Zyj5SPTo6Gq961avit3/7t+OiBfyTwo9+9KM477zz4l3velf8n//zf2L79u1x+eWXx8aNG+Occ85JNeiakc0e/ueWJIkkIkafW5yJiEiSw3d+7/citmxpjHkNAAAAZchmI9773ufz6HkrI0kykclEvO99ERdcUL/JlEmS2S9vwRtnMnH77bfHhRdeWHKdP/7jP4677rorvv3tbxeW/ef//J9jaGgo7rnnnqLbjI2NxdjYWOH+yMhI9Pf3x/DwcHR3d6cdbuXt2BFx9tkRcTio/7/nFmdmrtfZGdHinHAAAEBzyeUiDh0q9pOfRURP4d5990WcdVaVBrVAIyMj0dPTM2+HLnnp7dy5M7Zs2TJt2TnnnBM7d+4suc3AwED09PQUbv39/Us9zHR27VrYeun/3QIAAKBu5XIlfzLt3kLTqhYt+YnKdu/eHevXr5+2bP369TEyMhIHDx6Mrq6uWdtcffXVcdVVVxXu549U15yNGwvfroyIXzz3/ayE3rYt4owz5nyoTGbW8W0AAICak8vlIkmSyOVy027F/Nu/5U9E1jLllomINdPWm5JWdacmz/7d0dERHR0dyz2M+Z1++uFT1j3zTGSSJFbN/Hkmc/jn55wT0doaaWbai20AAGA5lIrnluc+2tpa5EPQLS0thVsmk4m3vKUlNm1qiWeeKT6BN59Mp5++1K9m6Sx5VG/YsCH27NkzbdmePXuiu7u76FHqutLaevgc8G996+G9Yepeko/hG24ofOK+VCDPFdulfia2AQCASsjH8syAnsvUeM4HdEuJ80iVkUx1ack/U7158+bYvn37tGX33ntvbN68eamfujouuijittsijjpq+vJNmw4vX8AZ0jOZTNHbXJIkKXoDAAAoJpfLxeTkZExMTMTY2FgcPHgwRkdH4+DBgzE2Nhbj4+MxOTk5LahbWlqira0t2tvbo7OzM7q6umLVqlXR1dUVHR0dsWLFimhtbS0Z1BEVSaaaVvaR6v3798eTTz5ZuP+jH/0oHnvssejr64ujjz46rr766njmmWfiH/7hHyIi4l3velf8zd/8TfzRH/1R/PZv/3Z89atfjS9+8Ytx1113Ve5VLLeLLjp8DvgHHzz8CfuNGw/PX1jkP7c4sg0AAJRr6tHmqUef5zLzyHP+6HOlLFEy1YSyL6m1Y8eOOPu5y0hNddlll8XNN98cb3/72+Opp56KHTt2TNvm93//9+O73/1ubNq0KT74wQ/G29/+9gU/50JPZd7sfGYbAACaRy3GcyNZaIcu6jrV1SKqF0dsAwBA/So3nvMfJxXPi7PQDq3Js39TWaaRAwBA7Zt5lu18QJcinmuDqG5iYhsAAKovTTxPPcO2eK4toppZxDYAACzezHjOB3Qp4rk+iWoWTGwDAMB0+UvbiufmJapZNLENAECjW0w8zwxo73kbi6hmyYhtAADqjXimXKKaqhPbAAAst3w4zwzohcbz1ID2nrS5iWpqxly/jEr9chPbAADMRTyz1EQ1daHco9tiGwCguYhnlouopq6JbQCA5jI1mmd+X8rMeM6Hs/eAVIKopiGJbQCA+iaeqReimqYitgEAakux6drimXoiqiEqF9tzPRYAQDOrVDy3tLRUcdQwP1ENc3D5LwCA8qSJ55nXdxbP1BNRDSmIbQCg2c2M53w4i2eajaiGChLbAECjEc8wN1ENVSC2AYBalzaeiwU0NBNRDctIbAMA1VYsnnO53JzbiGcoTVRDDRLbAMBizby+cznxPDOggdJENdQRsQ0AzCSeYXmJamgAYhsAGt9i43lqQAOVI6qhgYltAKg/M6/vLJ6htolqaEJiGwCWn3iGxiCqgQKxDQCVN/MSVWniOR/QQO0R1cC8xDYAzE88Q3MS1UBqYhuAZjTzRGH5gC4lk8nMur6zeIbGIaqBihPbADQC8QwshKgGqkZsA1CL0sTzzOs7i2doXqIaWHZiG4BqmBnP+YAuRTwDCyGqgZoltgEoV5Iksy5RJZ6BpSSqgbojtgFYTDzPDGh/C4DFENVAwxDbAI1HPAO1TlQDDW+uN1Gl3pSJbYDqyofzzIBeaDxPDWi/q4FqEtVAUyv36LbYBlgc8Qw0GlENUITYBlgc8Qw0C1ENUAaxDTDd1Gie+X0pM+M5H85+NwL1SFQDVIDYBhqdeAYoTlQDLCGxDdSbYtO1xTNAaaIaYBlUKrbneiyAuVQqnltaWqo4aoDaI6oBaohrbQOVliaeZ17fWTwDlCaqAeqA2AbmMzOe8+EsngGWlqgGqGNiG5qPeAaoLaIaoAGJbah/aeO5WEADsHRENUATEdtQe4rFcy6Xm3Mb8QxQO0Q1AGIbqmDm9Z3LieeZAQ1A7RDVAJQktqF84hmguYhqAMomtmHx8Tw1oAGoX6IagIoR2zSimdd3Fs8ATCWqAVhyYpt6IJ4BSENUA7BsxDbLYeYlqtLEcz6gAUBUA1BzxDaVUG48ZzKZWZeoEs8AzEdUA1A3xDbFzDxRWD6gSxHPAFSSqAag7ont5iCeAahFohqAhiW261OaeJ55fWfxDEC1iGoAmo7Yrg0z4zkf0KWIZwBqkagGgOeI7cpLkmTWJaoWGs8zA9p/UwBqkagGgHmI7fmJZwCalagGgJSaMbbFMwBMJ6oBoMLmisVS8VlrsZ0P55kBvdB4nhrQ4hmARtaSZqMbb7wxXvSiF0VnZ2eccsop8Y1vfKPkujfffHPhkhb5W2dnZ+oBA0A9m/k3MX8rJX9keOatUpIkiWw2G5OTkzE+Ph6HDh2KAwcOxIEDB+LQoUMxNjYWExMTkc1mC8+byWSitbU1VqxYER0dHdHZ2RkrV66MlStXRmdnZ7S3t0dbW1u0trYKagAaXtlHqr/whS/EVVddFTfddFOccsopccMNN8Q555wTTzzxRBxxxBFFt+nu7o4nnniicN8fWACYrtyp5OUe2a7Ekef8UWd/xwHgeWVH9V/91V/FO9/5zvit3/qtiIi46aab4q677oq///u/j/e///1Ft8lkMrFhw4YFP8fY2FiMjY0V7o+MjJQ7TABoCOXGdrHrO891jef8c4hnAEinrOnf4+Pj8fDDD8eWLVuef4CWltiyZUvs3Lmz5Hb79++PF77whdHf3x8XXHBBfOc735nzeQYGBqKnp6dw6+/vL2eYANAUcrlcTE5OxsTERIyNjcXBgwcL07bHx8enTdvOR3hLS0u0tbVFe3t7yWnbPgcNAAtXVlT/8pe/jGw2G+vXr5+2fP369bF79+6i2xx33HHx93//9/GlL30p/vf//t+Ry+XitNNOi5/+9Kcln+fqq6+O4eHhwu0nP/lJOcMEgIaR/8xzPpxnfuZ5ZjznjzoXC+di8Zx/jqX83DYANLIlP/v35s2bY/PmzYX7p512Wrz85S+PT33qU/HhD3+46DYdHR3R0dGx1EMDgJpR7LPO833meeYlqvK3cp+33J85ig0Azysrql/wghdEa2tr7NmzZ9ryPXv2LPgz0ytWrIhXv/rV8eSTT5bz1ADQEJYrnktpxmttA0AllfUXub29PU4++eTYvn17YVkul4vt27dPOxo9l2w2G48//nhs3LixvJECQB2ZOW07/3nnqdO2Jycnp12qqqWlZdqlqrq6umLVqlXR1dVVdNr2Uir30l/512waOQDNpuzp31dddVVcdtll8drXvjZe//rXxw033BCjo6OFs4FfeumlcdRRR8XAwEBERHzoQx+KU089NY455pgYGhqK66+/Pn784x/H5ZdfXtlXAgDLYOaR5/xR5/mOPBc7+lwPHNkGgOnKjuqLL744fvGLX8Sf/umfxu7du+Okk06Ke+65p3DysqeffnraG4PBwcF45zvfGbt3747e3t44+eST46GHHopXvOIVlXsVALDEmi2eyyW2AWhWmaQO5mWNjIxET09PDA8PR3d393IPB4AGViye57vO88zrOzdyPFdKmrcfYhuAalpohy752b8BoBZNPUmYeK4+R7YBaBSiGoCGtth4nhrQLD2xDUC9EdUANISZl6gSz41FbANQq0Q1AHVFPDOV2AZguYlqAGpSJeI5H9A0H7ENQLWIagCW1cxLVIlnlpLYBqDSRDUAVVFuPGcymVnXdxbPLBWxDUBaohqAipp5lu18QJcinqllYhuA+YhqAFJJE88zr+8snqlXYhuAPFENwJzEMyyc2AZoPqIagEiSZNZZtvMBXYp4hoUT2wCNS1QDNJHFxPPMgPbGHhZPbAPUP1EN0IDEM9Q3sQ1QP0Q1QB1bbDxPDWhvvKH2iW2A2iOqAepAPpxnBrR4BiLmDuRSvyfENkBliGqAGiKegUor9+i22AYoj6gGWAZTo3nm96XMjOd8OHujC6QhtgEqQ1QDLCHxDNQbsQ1QHlENUAHiGWh0YhugOFENUIZin3VOE88tLS1VHDXA0qlUbM/1WAC1TFQDFJEmnmde31k8A83M5b+AZiGqgaZW6hrP4hlgaYhtoNGIaqApiGeA2ia2gXolqoGGkjaeiwU0AMtPbAO1TlQDdalYPOdyuTm3Ec8AjUNsA7VCVAM1bTHxPDOgAWh8YhuoNlEN1ISZ13cWzwBUktgGloqoBqpqsfE8NaABYLHENrBYohpYEjOv7yyeAagnYhtYKFENLIp4BqCZiG1gJlENLMjMS1Slied8QANAoxHb0LxENTCNeAaAyhHb0PhENTSpcuM5k8nMur6zeAaAdMQ2NA5RDQ1u5lm28wFdingGgOUjtqH+iGpoEGnieeb1ncUzANQmsQ21S1RDnZkZz/mALkU8A0DjEtuw/EQ11KAkSWZdoko8AwALJbahekQ1LKPFxPPMgPYHDwCYj9iGyhPVUAXiGQCoZWIb0hPVUEH5cJ4Z0AuN56kB7Q8SALDc5no/Uur9jdim2YhqSEE8AwDNrtyj22KbRiWqYQ7iGQCgPGKbZiOqIWJaNM/8vpSZ8ZwPZ38AAABmE9s0KlFNUxHPAAC1RWxT70Q1DanYdG3xDABQP8Q29UJUU9cqFc8tLS1VHDUAAGlVKrbneiwoh6imLqSJ55nXdxbPAACNy7W2WS6impoyM57z4SyeAQBIQ2yz1EQ1yyJtPBcLaAAAKJfYplJENUtKPAMAUE/ENuUS1VREsXjO5XJzbiOeAQCoF2KbUkQ1ZZl5fedy4nlmQAMAQL0T24hqihLPAACQnthuHqK6yS02nqcGNAAAMDex3XhEdZOYeX1n8QwAALVDbNcvUd1gxDMAADQOsV37RHWdmnmJqjTxnA9oAACgvojt2iGqKyWbjXjwwYhduyI2bow4/fSI1tZFP2y58ZzJZGZdoko8AwBAc6jl2F6iZFp2oroStm2LuPLKiJ/+9PllmzZFfOITERddtKCHmHmisHxAlyKeAQCAhVru2K5AMtWsVB+cvfHGG+NFL3pRdHZ2ximnnBLf+MY35lz/1ltvjZe97GXR2dkZJ5xwQtx9992pBluTtm2LeOtbp+8dERHPPHN4+bZt0xbncrmYnJyM8fHxOHToUBw8eDBGR0fj4MGDMTY2FhMTE5HNZgtBnQ/ntra2aG9vj87Ozli5cmWsXLkyurq6oqOjI1asWBGtra2CGgAAKEv+YN3M21ySJCl6K6XMZKo7mWSuV1/EF77whbj00kvjpptuilNOOSVuuOGGuPXWW+OJJ56II444Ytb6Dz30UJxxxhkxMDAQb3nLW2Lr1q3x0Y9+NB555JF45StfuaDnHBkZiZ6enhgeHo7u7u5yhru0stmIF70o4qc/jSQiDkz5US5/O/LIyD36aCQtLfMeeZ55fWdHngEAgFpSZj5GNhvxildk4mc/m/mTlRGRiUzm8BHrH/2o9qaCL7RDy47qU045JV73utfF3/zN30TE4SOv/f398d73vjfe//73z1r/4osvjtHR0bjzzjsLy0499dQ46aST4qabbir6HGNjYzE2NjbtxfT399deVO/YEXH22RERsT8i1izrYAAAAOrF/ohYVbh3330RZ521bIMpaqFRXdb07/Hx8Xj44Ydjy5Ytzz9AS0ts2bIldu7cWXSbnTt3Tls/IuKcc84puX5ExMDAQPT09BRu/f395QyzenbtKnzreDIAAEA6U9Kq7pR1orJf/vKXkc1mY/369dOWr1+/Pr7//e8X3Wb37t1F19+9e3fJ57n66qvjqquuKtzPH6muORs3Fr5dGRFDcTiuW2JGZN99d8QZZ1RzZAAAAMvugQcizj232E9WTrs3Ja3qTk2e/bujoyM6OjqWexjzO/30wx8AeOaZyCRJ9Mz8ef4DAm96U+19QAAAAGCJvelNhWSKYh88zifT6adXf2yVUtb07xe84AXR2toae/bsmbZ8z549sWHDhqLbbNiwoaz160pr6+FzwEcc3humyt+/4QZBDQAANKVmSKayorq9vT1OPvnk2L59e2FZLpeL7du3x+bNm4tus3nz5mnrR0Tce++9JdevOxddFHHbbRFHHTV9+aZNh5fX+0XXAAAAFqHRkynVJbUuu+yy+NSnPhWvf/3r44YbbogvfvGL8f3vfz/Wr18fl156aRx11FExMDAQEYcvqXXmmWfGddddF+edd17ccsst8ZGPfKQxLqk1VTYb8eCDhz9hv3Hj4fkL9fzPLQAAABVUb8m00A4t+zPVF198cfziF7+IP/3TP43du3fHSSedFPfcc0/hZGRPP/10tLQ8fwD8tNNOi61bt8YHPvCBuOaaa+LYY4+NO+64Y8FBXTdaW2vvHPAAAAA1olGTqewj1cuhLo5UAwAA0DCW5DrVAAAAwPNENQAAAKQkqgEAACAlUQ0AAAApiWoAAABISVQDAABASqIaAAAAUhLVAAAAkJKoBgAAgJRENQAAAKQkqgEAACAlUQ0AAAApiWoAAABISVQDAABASqIaAAAAUhLVAAAAkJKoBgAAgJRENQAAAKQkqgEAACCltuUewEIkSRIRESMjI8s8EgAAAJpBvj/zPVpKXUT1vn37IiKiv79/mUcCAABAM9m3b1/09PSU/HkmmS+7a0Aul4uf/exnsWbNmshkMss9nJJGRkaiv78/fvKTn0R3d/dyDwdmsY9S6+yj1Dr7KLXOPkqtq6d9NEmS2LdvXxx55JHR0lL6k9N1caS6paUlNm3atNzDWLDu7u6a30FobvZRap19lFpnH6XW2UepdfWyj851hDrPicoAAAAgJVENAAAAKYnqCuro6Ihrr702Ojo6lnsoUJR9lFpnH6XW2UepdfZRal0j7qN1caIyAAAAqEWOVAMAAEBKohoAAABSEtUAAACQkqgGAACAlEQ1AAAApCSqy3TjjTfGi170oujs7IxTTjklvvGNb8y5/q233hove9nLorOzM0444YS4++67qzRSmlU5++hnPvOZOP3006O3tzd6e3tjy5Yt8+7TsFjl/h7Nu+WWWyKTycSFF164tAOk6ZW7jw4NDcUVV1wRGzdujI6OjviVX/kVf+9ZUuXuozfccEMcd9xx0dXVFf39/fH7v//7cejQoSqNlmbzwAMPxPnnnx9HHnlkZDKZuOOOO+bdZseOHfGa17wmOjo64phjjombb755ycdZSaK6DF/4whfiqquuimuvvTYeeeSReNWrXhXnnHNO/PznPy+6/kMPPRSXXHJJvOMd74hHH300Lrzwwrjwwgvj29/+dpVHTrModx/dsWNHXHLJJXHffffFzp07o7+/P970pjfFM888U+WR0yzK3UfznnrqqfiDP/iDOP3006s0UppVufvo+Ph4vPGNb4ynnnoqbrvttnjiiSfiM5/5TBx11FFVHjnNotx9dOvWrfH+978/rr322vje974Xf/d3fxdf+MIX4pprrqnyyGkWo6Oj8apXvSpuvPHGBa3/ox/9KM4777w4++yz47HHHov3ve99cfnll8dXvvKVJR5pBSUs2Otf//rkiiuuKNzPZrPJkUcemQwMDBRd/zd+4zeS8847b9qyU045Jfmv//W/Luk4aV7l7qMzTU5OJmvWrEk+//nPL9UQaXJp9tHJycnktNNOSz772c8ml112WXLBBRdUYaQ0q3L30U9+8pPJS17ykmR8fLxaQ6TJlbuPXnHFFcmv/uqvTlt21VVXJW94wxuWdJyQJEkSEcntt98+5zp/9Ed/lBx//PHTll188cXJOeecs4QjqyxHqhdofHw8Hn744diyZUthWUtLS2zZsiV27txZdJudO3dOWz8i4pxzzim5PixGmn10pgMHDsTExET09fUt1TBpYmn30Q996ENxxBFHxDve8Y5qDJMmlmYf/ad/+qfYvHlzXHHFFbF+/fp45StfGR/5yEcim81Wa9g0kTT76GmnnRYPP/xwYYr4D3/4w7j77rvj3HPPrcqYYT6N0Extyz2AevHLX/4ystlsrF+/ftry9evXx/e///2i2+zevbvo+rt3716ycdK80uyjM/3xH/9xHHnkkbN+sUElpNlHv/a1r8Xf/d3fxWOPPVaFEdLs0uyjP/zhD+OrX/1q/OZv/mbcfffd8eSTT8a73/3umJiYiGuvvbYaw6aJpNlH/8t/+S/xy1/+Mv7Df/gPkSRJTE5Oxrve9S7Tv6kZpZppZGQkDh48GF1dXcs0soVzpBqIiIjrrrsubrnllrj99tujs7NzuYcDsW/fvnjb294Wn/nMZ+IFL3jBcg8HisrlcnHEEUfEpz/96Tj55JPj4osvjj/5kz+Jm266abmHBhFx+PwpH/nIR+Jv//Zv45FHHolt27bFXXfdFR/+8IeXe2jQMBypXqAXvOAF0draGnv27Jm2fM+ePbFhw4ai22zYsKGs9WEx0uyjeR//+Mfjuuuui3/913+NE088cSmHSRMrdx/9wQ9+EE899VScf/75hWW5XC4iItra2uKJJ56Il770pUs7aJpKmt+jGzdujBUrVkRra2th2ctf/vLYvXt3jI+PR3t7+5KOmeaSZh/94Ac/GG9729vi8ssvj4iIE044IUZHR+N3fud34k/+5E+ipcUxNpZXqWbq7u6ui6PUEY5UL1h7e3ucfPLJsX379sKyXC4X27dvj82bNxfdZvPmzdPWj4i49957S64Pi5FmH42I+NjHPhYf/vCH45577onXvva11RgqTarcffRlL3tZPP744/HYY48Vbr/+679eODtof39/NYdPE0jze/QNb3hDPPnkk4V/8ImI+H//7//Fxo0bBTUVl2YfPXDgwKxwzv8jUJIkSzdYWKCGaKblPlNaPbnllluSjo6O5Oabb06++93vJr/zO7+TrF27Ntm9e3eSJEnytre9LXn/+99fWP/f/u3fkra2tuTjH/948r3vfS+59tprkxUrViSPP/74cr0EGly5++h1112XtLe3J7fddluya9euwm3fvn3L9RJocOXuozM5+zdLrdx99Omnn07WrFmTvOc970meeOKJ5M4770yOOOKI5M///M+X6yXQ4MrdR6+99tpkzZo1yf/9v/83+eEPf5j8y7/8S/LSl740+Y3f+I3legk0uH379iWPPvpo8uijjyYRkfzVX/1V8uijjyY//vGPkyRJkve///3J2972tsL6P/zhD5OVK1cmf/iHf5h873vfS2688caktbU1ueeee5brJZRNVJfpf/7P/5kcffTRSXt7e/L6178++frXv1742Zlnnplcdtll09b/4he/mPzKr/xK0t7enhx//PHJXXfdVeUR02zK2Udf+MIXJhEx63bttddWf+A0jXJ/j04lqqmGcvfRhx56KDnllFOSjo6O5CUveUnyF3/xF8nk5GSVR00zKWcfnZiYSP77f//vyUtf+tKks7Mz6e/vT9797ncng4OD1R84TeG+++4r+v4yv19edtllyZlnnjlrm5NOOilpb29PXvKSlySf+9znqj7uxcgkiXkfAAAAkIbPVAMAAEBKohoAAABSEtUAAACQkqgGAACAlEQ1AAAApCSqAQAAICVRDQAAACmJagAAAEhJVAMAAEBKohoAAABSEtUAAACQ0v8f/zQwOIcLhMsAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1200x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_transporte(mu, nu, T, threshold=1e-8)"
   ]
  }
 ],
 "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}