summaryrefslogtreecommitdiff
path: root/notebooks
diff options
context:
space:
mode:
authorMokhtar Z. Alaya <mzalaya@Mokhtars-MacBook-Pro.local>2020-01-08 18:54:07 +0100
committerMokhtar Z. Alaya <mzalaya@Mokhtars-MacBook-Pro.local>2020-01-08 18:54:07 +0100
commit3e77515b4f19cf1c37b2f971a54b2fe5efe9daef (patch)
tree80ed034012902c4e83c00aa4f0469f1ac121c438 /notebooks
parente00f46aa2ea11f0e88a5b2005caa7518ca109357 (diff)
add illustration for screenkhorn
Diffstat (limited to 'notebooks')
-rw-r--r--notebooks/plot_screenkhorn_1D.ipynb207
1 files changed, 207 insertions, 0 deletions
diff --git a/notebooks/plot_screenkhorn_1D.ipynb b/notebooks/plot_screenkhorn_1D.ipynb
new file mode 100644
index 0000000..6126346
--- /dev/null
+++ b/notebooks/plot_screenkhorn_1D.ipynb
@@ -0,0 +1,207 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "%matplotlib inline"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "# 1D Screened optimal transport\n",
+ "\n",
+ "\n",
+ "This example illustrates the computation of Screenkhorn: Screening Sinkhorn Algorithm for Optimal transport.\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Author: Mokhtar Z. Alaya <mokhtarzahdi.alaya@gmail.com>\n",
+ "#\n",
+ "# License: MIT License\n",
+ "\n",
+ "import numpy as np\n",
+ "import matplotlib.pylab as pl\n",
+ "import ot\n",
+ "import ot.plot\n",
+ "from ot.datasets import make_1D_gauss as gauss\n",
+ "from ot.bregman import screenkhorn"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Generate data\n",
+ "-------------\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "#%% parameters\n",
+ "\n",
+ "n = 100 # nb bins\n",
+ "\n",
+ "# bin positions\n",
+ "x = np.arange(n, dtype=np.float64)\n",
+ "\n",
+ "# Gaussian distributions\n",
+ "a = gauss(n, m=20, s=5) # m= mean, s= std\n",
+ "b = gauss(n, m=60, s=10)\n",
+ "\n",
+ "# loss matrix\n",
+ "M = ot.dist(x.reshape((n, 1)), x.reshape((n, 1)))\n",
+ "M /= M.max()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Plot distributions and loss matrix\n",
+ "----------------------------------\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAADCCAYAAABnjpSEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3iUZfbw8e8hoQRBkGKhgxQFUoBIExRYpCwuiivSFLAs4q6K+pMVyyKLZXV1Zd1F8WXtFRQbrgVEYUFQMAFCR6oSbNQAQiDlfv84MzGElEkyyTMzOZ/rmmvaM/Ocqee5uzjnMMYYY4qrktcBGGOMCU+WQIwxxpSIJRBjjDElYgnEGGNMiVgCMcYYUyKWQIwxxpRItNcB5FWvXj3XrFkzr8MwxhgDJCcn73XO1c/vvpBLIM2aNSMpKcnrMIwxxgAi8m1B91kVljHGmBIJKIGIyAAR2SwiW0VkUj73VxWR2b77l4tIM9/tlUXkJRFZKyIbReTu4IZvjDHGK0UmEBGJAp4CBgJtgREi0jbPZtcDB5xzLYFpwKO+24cCVZ1zsUAn4EZ/cjHGGBPeAmkD6Qxsdc5tBxCRWcBlwIZc21wGTPFdngNMFxEBHHCaiEQDMcAJ4FBwQjfGBEtGRgapqamkp6d7HYrxSLVq1WjUqBGVK1cO+DGBJJCGwK5c11OBLgVt45zLFJE0oC6aTC4DfgCqA7c75/YHHJ0pkXnzYMMGmDABKlkrlwlAamoqNWvWpFmzZuixn6lInHPs27eP1NRUmjdvHvDjAvl7ye/blHcK34K26QxkAQ2A5sD/iUiLU3YgMk5EkkQkac+ePQGEZAryzDPw29/CHXfAyJFw/LjXEZlwkJ6eTt26dS15VFAiQt26dYtdAg0kgaQCjXNdbwR8X9A2vuqqWsB+YCTwiXMuwzn3M7AUSMy7A+fcTOdconMusX79fLsbmyI4B/fdBzfdBAMHwoMPwuzZejktzevoTDiw5FGxleTzDySBfA20EpHmIlIFGA7MzbPNXGCM7/KVwOdOFxr5Dugj6jSgK7Cp2FGaIt16Kzz0EFx3Hbz3Htx7L7z8MixZAhddBEeOeB2hMYV76KGHaNeuHXFxcSQkJLB8+XKvQzpFjRo1APj++++58sorC9zu4MGDPP3004U+V/fu3QFYtGgRl156abHieO+999iw4ddm6MmTJ7NgwYJiPUdQOOeKPAG/Bb4BtgH3+m6bCgz2Xa4GvAVsBVYALXy31/Ddvh5tdJ9Y1L46derkTPHs3OmciHPjxjmXnX3yfe++6xw49//+nzexmfCwYcMGT/e/bNky17VrV5eenu6cc27Pnj1u9+7dpX7ejIyMUj9HbqeddlpA2+3YscO1a9cu3/syMzNPur5w4UI3aNCgYsUxZswY99ZbbxXrMYHI73sAJLkC/q8DamJ1zn3knGvtnDvXOfeQ77bJzrm5vsvpzrmhzrmWzrnOztdjyzl3xHd7O+dcW+fcY0HJeuYkL7yg53ffDXlLoZddBu3awXPPlX9cxgTqhx9+oF69elStWhWAevXq0aBBAwA+++wzOnToQGxsLNdddx3HfQ17zZo1Y+/evQAkJSXRq1cvAKZMmcK4cePo168fo0ePJisrizvvvJPY2Fji4uL497//DUBycjIXX3wxnTp1on///vzwww+nxLVjxw66devGBRdcwF/+8pec23fu3En79u0BWL9+PZ07dyYhIYG4uDi2bNnCpEmT2LZtGwkJCUycOJFFixbRu3dvRo4cSWxsLPBraQbg0KFDDBkyhLZt2zJ+/Hiys7NP2WbOnDmMHTuWZcuWMXfuXCZOnEhCQgLbtm1j7NixzJkzp8j36/7776djx47ExsayaVPpK4NCbioTUzxZWZpA+vaF/KYQE4Hrr9dG9XXrwPedN6ZAt90Gq1cH9zkTEuCf/yz4/n79+jF16lRat25N3759GTZsGBdffDHp6emMHTuWzz77jNatWzN69GhmzJjBbbfdVuj+kpOT+eKLL4iJiWHGjBns2LGDVatWER0dzf79+8nIyOCWW27h/fffp379+syePZt7772X559//qTnmTBhAjfddBOjR4/mqaeeyndfzzzzDBMmTGDUqFGcOHGCrKwsHnnkEdatW8dq3xu5aNEiVqxYwbp16/Lt5bRixQo2bNhA06ZNGTBgAO+8806BVWTdu3dn8ODBXHrppadsU9T7Va9ePVauXMnTTz/N448/zrPPPlvo+1gU6+QZ5j77DL77TpNEQa65BipXtlKICV01atQgOTmZmTNnUr9+fYYNG8aLL77I5s2bad68Oa1btwZgzJgxLF68uMjnGzx4MDExMQAsWLCA8ePHEx2tx8t16tRh8+bNrFu3jksuuYSEhAQefPBBUlNTT3mepUuXMmLECACuueaafPfVrVs3Hn74YR599FG+/fbbnP3m1blz5wK7yHbu3JkWLVoQFRXFiBEj+OKLL4p8jfkp6v264oorAOjUqRM7d+4s0T5ysxJImHvuOahTBy6/vOBt6tXTqqxXXoFHHgFfLYEx+SqspFCWoqKi6NWrF7169SI2NpaXXnqJhISEArePjo7OqerJ2/30tNNOy7nsnDulh5Fzjnbt2vHll18WGVdRvZNGjhxJly5d+PDDD+nfvz/PPvssLVqcMlrhpJiK2of/eu7bA+liq00WBfNXEUZFRZGZmVnk8xXFSiBhbN8+7XF1zTVFJ4UbbtDt5+btP2dMCNi8eTNbtmzJub569WqaNm3Keeedx86dO9m6dSsAr7zyChdffDGgdfrJyckAvP322wU+d79+/XjmmWdy/jD3799PmzZt2LNnT04CycjIYP369ac89sILL2TWrFkAvPbaa/k+//bt22nRogW33norgwcPZs2aNdSsWZPDhw8H/PpXrFjBjh07yM7OZvbs2fTo0QOAs846i40bN5Kdnc27776bs31Bz1/Y+1UWLIGEsVdfhRMnCq++8uvbFxo3tmosE5qOHDnCmDFjaNu2LXFxcWzYsIEpU6ZQrVo1XnjhBYYOHUpsbCyVKlVi/PjxANx///1MmDCBnj17EhUVVeBz33DDDTRp0oS4uDji4+N5/fXXqVKlCnPmzOGuu+4iPj6ehIQEli1bdspjn3zySZ566ikuuOAC0goYUDV79mzat29PQkICmzZtYvTo0dStW5cLL7yQ9u3bM3HixCJff7du3Zg0aRLt27enefPmDBkyBIBHHnmESy+9lD59+nDOOefkbD98+HAee+wxOnTowLZt23JuL+z9KgtSVJGnvCUmJjpbD6RozkFcHMTEwIoVgT3m/vvhgQdg505o0qRMwzNhZuPGjZx//vleh2E8lt/3QESSnXOnDAAHK4GErXXr9HTddYE/5tprNfG88UbZxWWMqTgsgYSpzz/X80GDAn9Ms2Y6JmThwjIJyRhTwVgCCVMLF8K552q7RnH07g1ffAEZGWUTlzGm4rAEEoays2HxYvANvC2WXr3gl1/AmpmMMaVlCSQMpaTAgQMlSyD+Hn1WjWWMKS1LIGFo0SI9L0kCqVcPYmN/fQ5jjCkpSyBhaOFCaNkSGjUq2eN79YKlS3UMiTGhYN++fSQkJJCQkMDZZ59Nw4YNc66fKKMv6sqVK/nkk08C2rZHjx4581r179+/0EGCTzzxRKGjxq+99lo2b95MZmYmtWvXLlXM7777Lo895t0ctZZAwkxWlrZ/9O5d8ufo3RuOHoWvvw5eXMaURt26dVm9ejWrV69m/Pjx3H777TnXq1SpUuTjs7Kyir3P4iSQ3ObNm0fNmjULvL+wBJKVlcULL7xAmzZtir1fODXmIUOGBDRQsaxYAgkzq1frCoMlqb7yu+giPbd2EBMOfve739GpUyfatWuXM3us/+j9vvvuo3PnzqxYsYK5c+fSpk0bevbsyS233MLlvgnijhw5wtixY+ncuTMdOnTggw8+4NixY0ydOpXXXnuNhISEnKnQ/Y4ePcrQoUOJi4tj+PDhJyWERo0acfDgQQ4fPszAgQOJj4+nffv2zJkzh2nTpvHzzz/Ts2dP+vbtm2+cuUszALfffjsdO3bkkksuYd++fcDJJZ4ff/yRli1b5hvzs88+mzPT7o4dO+jduzdxcXFccsklOZNDXn311UyYMIHu3bvTokWLk6ZEKS2bTDHMlKb9w69uXR3FvmiRLoNrzEm8mM+9EC+99BJ16tTh6NGjJCYm8vvf/56aNWuSlpZGx44defDBBzl69CitW7dm6dKlNGnShKuuuirn8VOnTmXAgAG8+OKLHDhwgC5durBmzRomT57MunXr+Gc+cU2fPp0zzjiDNWvWsGrVKhITTx2I/dFHH9GsWTM+/vhjANLS0qhVqxb/+Mc/WLJkCbVr1yYzM/OkOPNKS0uja9euTJs2jcmTJ/PAAw/kGw9ATEzMKTHnno79j3/8IzfccAOjRo1i5syZ3HbbbTmJ8eeff2bp0qWsXbuWq666KmeqlNKyEkiYWbgQWrcG31o7Jda7t7aD+NaaMSZkTZs2jfj4eLp160ZqamrO3E9VqlTJ+SPcsGEDbdq0oWnTpohIzhTsAPPnz+ehhx4iISGB3r17k56eznfffVfoPhcvXszVV18NQIcOHWjXrt0p28TFxfHJJ58wadIkli5dSq1atfJ9rtxx5hUdHc3QoUMBLSmUdBp3gOXLlzN8+HAARo8ezZIlS3Luu/zyyxER4uLi2L17d4n3kZeVQMJIZqauce77jpRKr17w5JM6j1bPnqV/PhNBvJrPPR8LFixg8eLFfPXVV8TExNCjR4+c6qSYmJic6c4Lm9PPOcd7773Hueeee9LtRa0rUtQ07ueffz5JSUl89NFHTJw4kUsvvZR77rnnlO1yx1nUPvzXC5uqviSq5pquO5jzH1oJJIysWgWHDpWuAd3voot0tUJrBzGhLC0tjTp16hATE8P69ev5uoCeH+3atWPz5s3s2rUL5xyzZ8/Oua9///7861//yrm+atUqoOAp0QEuuuiinOnbU1JS8p3qfffu3dSoUYNrrrmGO+64g5UrVxb5vHllZGTwzjvvAPD666/nTOOee6r63O0zhT13165defPNNwF49dVXucjf2FmGLIGEEf8BUzCm969TR9tBAljczRjPDBo0iKNHjxIfH8/UqVPp0qVLvttVr16d6dOn07dvX3r27EmDBg1yqpTuv/9+jh49SmxsLO3atWPKlCkA9OnTh5SUFDp06HBKI/rNN9/Mvn37iIuLY9q0afm2gaSkpHDBBReQkJDA3//+95zSx7hx4+jbty99+/Yt8vXVqlWLlStX0rFjR7744gvu8zVKTpw4kSeffJLu3btz4MCBnO0Li3n69OnMnDmTuLg4Zs+ezbRp04rcf2nZdO5hZMQIWLYMvv02OM83fjzMmqWj2osorZsIFwnTuR85coQaNWrgnOPGG28kNjaWW265xeuwwopN5x7BkpIgnwOhEktM1C7BudajMSZszZgxg4SEBNq2bcuxY8f4wx/+4HVIEc8a0cPEgQOwdWvx1v8oij8ZJSXpyHZjwtnEiRM9HVRXEVkJJEz42ueCWgJp107XUrcaQ2NMSVgCCRP+P/lOnYL3nJUr6/guSyAGgtu904Sfknz+lkDCRFIStGihvaeCKTERkpN1jRFTcVWrVo19+/ZZEqmgnHPs27ePatWqFetx1gYSJpKSoHPn4D9vYiI89RR88w2cd17wn9+Eh0aNGpGamsqePXu8DsV4pFq1ajQq5hTflkDCwN69sHMn3HRT8J/bXyWWlGQJpCKrXLkyzZs39zoME2YCqsISkQEisllEtorIpHzuryois333LxeRZrnuixORL0VkvYisFZHilZEMvgGpQW1A9zv/fIiJsXYQY0zxFZlARCQKeAoYCLQFRohI2zybXQ8ccM61BKYBj/oeGw28Cox3zrUDegEZQYu+gvAnkI4dg//c0dHQocOv+zDGmEAFUgLpDGx1zm13zp0AZgGX5dnmMuAl3+U5wG9EZwXrB6xxzqUAOOf2OeeKv/JLBZeUBK1aQTEXLwtYYqJ2Ey7BmjzGmAoskATSENiV63qq77Z8t3HOZQJpQF2gNeBEZJ6IrBSRP5c+5Ion2CPQ80pM1BUKN20qu30YYyJPIAkkv1mS8vb1K2ibaKAHMMp3PkREfnPKDkTGiUiSiCRZL5CT/fQT7NpV9gkErB3EGFM8gSSQVKBxruuNgO8L2sbX7lEL2O+7/X/Oub3OuaPAR8ApNfnOuZnOuUTnXGL9+vWL/yoiWFk2oPu1bg01algCMcYUTyAJ5GuglYg0F5EqwHBgbp5t5gJjfJevBD53OiJpHhAnItV9ieViYENwQq8YkpJ0ptwOHcpuH1FR2kBvCcQYUxxFJhBfm8bNaDLYCLzpnFsvIlNFZLBvs+eAuiKyFbgDmOR77AHgCTQJrQZWOuc+DP7LiFxJSdCmDdSsWbb7SUzUZbAzM8t2P8aYyBHQQELn3Edo9VPu2ybnupwODC3gsa+iXXlNCaxeXT5LznboAOnpsHmzTrJojDFFsbmwQtj+/dqAHh9f9vvy7yMlpez3ZYyJDJZAQtiaNXpeHgnkvPOgShVLIMaYwFkCCWGrV+t5QkLZ76tyZa268u/TGGOKYgkkhKWkwFln6ak8xMdbCcQYEzhLICEsJaV8qq/84uN14OKPP5bfPo0x4csSSIjKyID168un+srPvy8rhRhjAmEJJERt2gQnTpR/CQQsgRhjAmMJJET5/8TLM4GccQY0bmwJxBgTGEsgIWr1aqhaVUehl6eEBOuJZYwJjCWQEJWSAu3b64JP5Sk+Xkejp6eX736NMeHHEkgIcq78e2D5xcfrwlLr15f/vo0x4cUSSAj68UfYs6d8e2D5WU8sY0ygLIGEIH8bhBclkBYtdG0QawcxxhTFEkgI8h/9x8WV/74rVYLYWCuBGGOKZgkkBKWkQLNmULu2N/tPSNAYXN6Fi40xJhdLICFo9Wpvqq/84uMhLQ2+/da7GIwxoc8SSIg5dgy++cb7BAJWjWWMKZwlkBCzbh1kZ3vTA8svNlbXYbcEYowpjCWQEONlDyy/006DVq2sJ5YxpnCWQEJMSgrUrKmN6F6ytUGMMUWxBBJiUlK0+24ljz+Z+HjYvh0OHfI2DmNM6LIEEkKyszWBeNn+4eePwb8uuzHG5GUJJITs3AmHD3vb/uFnPbGMMUWxBBJCvFgDpCANG0KdOpZAjDEFswQSQlJStO2jfXuvI9FuvNaQbowpjCWQELJ6tXafrV7d60hUfDysXavTuxtjTF6WQEJIqDSg+yUk6Mj4LVu8jsQYE4osgYSIgwe1ET0U2j/8rCHdGFMYSyAhwt9dNpQSyPnn65K6NiLdGJOfgBKIiAwQkc0islVEJuVzf1URme27f7mINMtzfxMROSIidwYn7MgTSj2w/KpW1SRiJRBjTH6KTCAiEgU8BQwE2gIjRKRtns2uBw4451oC04BH89w/Dfi49OFGrpQUqFcPGjTwOpKT+dcGMcaYvKID2KYzsNU5tx1ARGYBlwEbcm1zGTDFd3kOMF1ExDnnRORyYDvwS9CijkD+NUBEvI7kZPHx8MorukZ7/fpeR2MCkp4OmzfDrl16SkvTI5PGjaFpU2jePPS+aCYsBZJAGgK7cl1PBboUtI1zLlNE0oC6InIMuAu4BCiw+kpExgHjAJo0aRJw8JEiM1Oncf/Tn7yO5FS5G9L79vU2FlOIY8fgk0/grbfggw/gyJGCt23RAq66Sk8JCZZMTIkF0gaS37cr72KnBW3zV2Cac66QbzM452Y65xKdc4n1K+Bh7jffwPHjodX+4Wc9sULc8eMwbRo0agRXXAHz58PIkTB7Nnz1FaSmajL55hv47DOYMUMHGz32GHTsCF27wuLFXr8KE6YCKYGkAo1zXW8EfF/ANqkiEg3UAvajJZUrReTvQG0gW0TSnXPTSx15BFm1Ss9DaQyIX/36Wvvhj9GECOdg1iy45x7t/92vH/zf/0GfPtp1Lq9WrfTUpw+MHw9798Kbb8LDD8PFF8PvfgePPqq9JowJUCAlkK+BViLSXESqAMOBuXm2mQuM8V2+EvjcqZ7OuWbOuWbAP4GHLXmcKjkZqlWDtnm7JoSITp00RhMi9u+Hyy/XkkatWjBvnp769cs/eeSnXj344x+1ZPK3v8H//qdHME8/rcnJmAAUmUCcc5nAzcA8YCPwpnNuvYhMFZHBvs2eQ9s8tgJ3AKd09TUFS07WqqJAf/vlrVMnbZM9fNjrSAxffgkdOsDHH8MTT8DKlZo4Sqp6dZg0Sacb6NtXG+KGDtWRrcYUIaC/LOfcR8BHeW6bnOtyOjC0iOeYUoL4Il52tlYPXXON15EUrFMnPShNSYEePbyOpgKbMQNuvVV7Uy1dChdcELznPvNMbXx/4gm4+25NTB9+aFVaplA2Et1jW7fqkX3Hjl5HUjB/bFaN5RHnYMoUrXIaOFCPOIKZPPwqVYI774QlS7RXV8+esGJF8PdjIoYlEI/5/5Q7dfI2jsI0aABnn20JxBNZWXDzzfDXv8J118E772i7R1nq2lVLOLVqaaP7p5+W7f5M2LIE4rHkZJ0ypF07ryMpnDWkeyArC0aP1obtP/8Znn22/BrKWrTQJNKyJQwaBO+9Vz77NWHFEojHVq6EuDioXNnrSArXqRNs2gS/2HwC5cM5uOkmeP117Wr76KPlP+Dv7LNh0SL98IcNgwULynf/JuRZAvGQc5pAQrn6yq9TJ23wtwGF5cA5LXH85z9w773aqO2V2rXho4+gTRvtOvzll97FYkKOJRAPbdum0xSFSwIBq8YqF3/7Gzz+uHapfeABr6OBM87QEe7nnAO//a0dRZgclkA8FA4N6H4NGsBZZ1kCKXMvv6yljquvhn/9K3TmqTr7bK3CqlFDk8ju3V5HZEKAJRAPJSdDlSqh34AO+j/WsaMlkDL1xRdwww3a8+n557VbbShp2lSrsw4dgsGDrUHMWALx0sqVEBurSSQcdOoEGzbA0aNeRxKBtm+HIUN0qvU5c0K3V0VsrM7BtXq19hDLzvY6IuMhSyAeCacGdD9/Q7p/+V0TJGlpOplhVhb897/a5hDKBg2Cf/xDx6Tcd5/X0RgPWQLxyI4dcOBA+CUQsGqsoMrO1nlsvvkG3n5bZ8wNBxMmwLhx2uD/1lteR2M8YgnEI+HUgO7XqJFO724JJIgeeeTXOah69/Y6msCJwL//Dd26wbXXwsaNXkdkPGAJxCPJyVrN3b6915EETsRGpAfV/PlaBTRypE5XEm6qVNHSx2mnafvNoUNeR2TKmSUQj/hn5a5a1etIiqdLF11+16Z2L6WdO2HECO2CN3Nm6HTXLa6GDXX1w61btSRia4lUKJZAPJCRoZOcdu/udSTF1727VtsvX+51JGHsxAldjzwzUxuiTzvN64hKp1cvnWrlnXd0eV1TYVgC8cDq1ZCeDhde6HUkxde1qx4sL1vmdSRh7K674Ouv4YUXwqfRvCh33KHVWHfdZUcXFYglEA/4/3zDsQRy+uk6FMASSAm99x7885+6MNQVV3gdTfCIwHPPaU+LYcO0i6GJeJZAPLBsmQ7qbdDA60hKpnt3bcOxMWTFtHOnthN06gR//7vX0QTfGWdoe8j331t7SAVhCaScOafLLIRj6cOve3ftcLN+vdeRhJGMDBg+XLPum2+GX++JQHXurMnx/fd1Li8T0SyBlLNdu3QeunBPIGDVWMVy333aNvDcc7pYUySbMEHnyvrzn3X5XROxLIGUM/+fbjg2oPu1aKEz81oCCdD8+XpUfuONcOWVXkdT9kR0Msgzz9T2EOvzHbEsgZSzZcu012ZsrNeRlJyIlkIsgQTgxx91qpL27StWF9e6deG113TRm3AcJGkCYgmknC1dqoPxymtp67LSvbuOHfvpJ68jCWHZ2Tpj7eHDOoNtTIzXEZWviy6CyZN1jZNXXvE6GlMGLIGUoyNHdDG3cG7/8PO/BlvhtBCPPw6ffqrddsNh0ZeycN99cPHF8Mc/wpYtXkdjgswSSDn6+mudsTsSEkjHjjoVklVjFWD5cl1ZcOhQ+MMfvI7GO1FR8Oqr+mUZPhyOH/c6IhNElkDKkf/Ptls3b+MIhmrVIDHREki+Dh7UP8uGDcN7nqtgadRIR92vXAl33+11NCaILIGUo2XLtCajdm2vIwmO7t0hKckOKk/inPa22rUL3ngjcj7s0ho8GG65RTsSfPih19GYILEEUk4yMmDJEujRw+tIgqdHD00eX33ldSQh5LnndKDggw9GRlEzmP7+d0hIgDFjdDCUCXsBJRARGSAim0Vkq4hMyuf+qiIy23f/chFp5rv9EhFJFpG1vvM+wQ0/fHz1lXbG6dfP60iCp3dv7U02b57XkYSItWv1KLtvXx1EZ05WrZr2RktP1zVQMjO9jsiUUpEJRESigKeAgUBbYISItM2z2fXAAedcS2Aa8Kjv9r3A75xzscAYoML25Zs3T9sTf/MbryMJntNP14NsSyBoF7urrtIqq1dfhUpWuM9XmzbwzDOweDH89a9eR2NKKZBveWdgq3Nuu3PuBDALuCzPNpcBL/kuzwF+IyLinFvlnPved/t6oJqIROgkQIWbP1+nQq9Vy+tIgqt/f20b3bPH60g89qc/webNOnjurLO8jia0XX01XHcdPPQQLFjgdTSmFAJJIA2BXbmup/puy3cb51wmkAbUzbPN74FVzrlTmlxFZJyIJIlI0p4I/Cfau1cbmyOp+srP/5o+/dTbODz14os6WG7yZOhTYWtpi+ff/4a2bWHUKPjhB6+jMSUUSALJrw9i3nmaC91GRNqh1Vo35rcD59xM51yicy6xfv36AYQUXhYs0M45/ft7HUnwdeyos1ZU2GqsNWt0kFyvXvCXv3gdTfioXl07Gxw5ol2erT0kLAWSQFKBxrmuNwK+L2gbEYkGagH7fdcbAe8Co51z20obcDiaP1+XSkhM9DqS4IuK0jbj+fMr4PIPaWnw+99ru8cbb+ibYQLXtq2Ok1m8GO65x+toTAkEkkC+BlqJSHMRqQIMB+bm2WYu2kgOcCXwuXPOiUht4EPgbufc0mAFHU6c06Pzvn0j9/+lf3+dM9G/+UsAABCRSURBVHDtWq8jKUfOwdixsGOHHkmffbbXEYWnUaO0BPfYY7qmugkrRSYQX5vGzcA8YCPwpnNuvYhMFZHBvs2eA+qKyFbgDsDf1fdmoCXwFxFZ7TudGfRXEcLWr9cF2iKx+srP3w5SoaqxHn9cl6d97LHIGtzjhSee0IWoxo6Fb77xOhpTDOJCrN4hMTHRJSUleR1G0PzjH3DnnfDdd9C4cdHbh6v27fUgvEJ0qlmwAAYMgCFDtPRR0acqCYbvvtMGtbPO0kFTNWt6HZHxEZFk51y+FfDWWb2MzZ8P558f2ckDtIS1ZAkcPep1JGVs2zYd73HeebpokiWP4GjSRJPx5s26fkp2ttcRmQBYAilDx45p+2AkV1/59e8PJ07AokVeR1KGDh+Gyy7TpDF3rh0lB1ufPjpX1vvvw/33ex2NCYAlkDL04Yc6a8OgQV5HUvZ69tT/07ff9jqSMpKdrUfGmzbpkXKkr2vulZtv1kGGDz4Ib73ldTSmCJZAytBrr2m7QO/eXkdS9mJitElgzhxNmhHnnnv0yPiJJyJrPppQIwJPP61TPY8ZAytWeB2RKYQlkDJy4AB89JGOkYrU7rt5jRoFhw7p644ozzwDjz4K48frZImmbFWtCu++q0dfv/uddpU2IckSSBl5+21tExg50utIyk+fPnDmmVryihj//a/OczVokE6/YY3m5ePMM+Hjj3UdhIEDYf9+ryMy+bAEUkZefx1atYrM0ecFiY7WEteHH+qifGEvORmGDdM1LGbN0hdoyk+bNlptuGOHdl6IyLrR8GYJpAzs3q29kUaOrHgHrCNH6iJTYT+oeMMGHetRv76WQmrU8DqiiqlnT3jpJfjiC+0+nZHhdUQmF0sgZWDWLJ3poiJVX/l17gznnqslsLC1fTtccomWOD79FM45x+uIKrbhw+Gpp+CDD2D0aMjK8joi42MJpAy89ppWXbVu7XUk5U9EE+fnn+sULmEnNVV7WaWna/Jo1crriAzofFmPPqpHZzfeaAMNQ4QlkCDbuBFWrdIeSRXVqFFaAps92+tIismfPPbt04m92rf3OiKT25//DPfdp+vO33yzJZEQYAkkyJ55Rms+hg3zOhLvtGkDF1yg70XY1DZs36717T/8oL1/KlLvh3Aydaomkhkz4Prrw+gLFpksgQTRzz/Df/6jA5YrerX5xIk6seq773odSQA2bYKLLtJBLJ9/Dhde6HVEpiAi8MgjMGWKrgQ5apQ1rHvIEkgQPfmkVp3fdZfXkXjviiu0Dejhh0N8oamkJLj4Yv0TWrTISh7hQETnynrsMa0nHTJEVzY05c4SSJCkpcH06XDllVqFU9FFRcGkSdoeFLLrhLz7rpY8qlfXWS9jY72OyBTHnXdqPenHH+vnuHu31xFVOJZAguTpp7UG5O67vY4kdIwaBY0aaSkkpDinC0L9/vcQF6frT1jWD0833qjjdLZsgS5dYPVqryOqUCyBBMHRozoL9YAB0KGD19GEjipVtC1kyRI9hYRfftGV7yZO1OLiwoW6iJEJXwMH6kBDEV0dMqwHIYUXSyBB8OyzsGePTthqTnbDDVCvHjz0kNeRoKPLO3eGV17RRthZs3QaYRP+4uN15t4OHbToe+ONNvVJObAEUkq7d2t7Xq9e2gvUnKx6dT3YnzfPw7VCnIMXXtC+xXv36jKR998PlezrH1HOOUdLlHfdBTNnQteuOjDLlBn7BZWCc/CHP+jcTzNneh1N6Lr9dl3u+qabtKRWrnbt0pl0r7tOE8iqVdC3bzkHYcpNdLR28/3vf3VgaEIC/O1vkJnpdWQRyRJIKTz/vHYAeeQRm/GiMJUr63x4aWmaRMqlW292tmb1du3gf//TPtaffw4NGpTDzo3nBg2C9et1PZF77tHSiDWwB50lkBL67js9su7VS2dVMIVr3x7++letxirzKU7+9z8dz3HjjXq+di3ceqtVWVU0Z52lS2S+9Zb+YDt21CqDn37yOrKIYb+oEjhxQjvyOKelEPtfCsydd2pPyz/9CXbuLIMdrF+vXXN79dK2jtdfh88+s/XLK7orr4TNm2HCBB293rKl9uo4dMjryMKe/fUV04kTuizBwoW6QF3z5l5HFD6io7UqKztb14kPWhJJTtah7+3bawP5gw/qH8aIERVvQRaTvzPO0L7269fr0pn33QdNm2pnClvtsMQsgRSDP3m8/74mj7FjvY4o/LRpo7OkHzxYyiRy4oTWhfXurdVUCxfC5Mn6hPfea91zTf5at9Yf8Ndf63dn6lRo0kTXu1+1yuvowo4lkAAdPXpy8rB2j5JLTIQFCzSJ9OoFW7cG+EDndO6qiROhcWNdaGjnTl0n4ttvtZGlbt0yjNxEjMREXTZz7Vr9Yb/8sraRdOmi00pYO0lAxIXYTHeJiYkuKSnJ6zBOsmABjBunSzNPn651+Kb0kpN14b/jx+GBB7Sd+5Rlx0+cgGXLtLvbW2/phxAdrb1sbrpJn8AaoUxpHTigA0xnztRqrkqVdJLNIUOgf3/tZllBq0NFJNk5l+8so5ZACvHjjzq31Ysvasn3P//ROdtM8KSm6mJzH3ygB4UzpqWTKMmaNJYs0aqpI0c0afzmN3q0ePnlUKeO16GbSOScJpC33tIq0s2b9fZmzfRgpXt3PVWghFLqBCIiA4AngSjgWefcI3nurwq8DHQC9gHDnHM7fffdDVwPZAG3OucKnZvV6wSSlaXtsP/5j/6pOafr10yeDNWqeRZW5MnK0q6VW7bg1q1n59w1HFq6hvMz11IFXd8hu0VLKvW/BPr104bP00/3OGhT4Wzbpn8I8+bpdP9paXp73bo6SDEuTmdxPv98TSp16kRcYilVAhGRKOAb4BIgFfgaGOGc25Brmz8Ccc658SIyHBjinBsmIm2BN4DOQANgAdDaOVfgMmLlmUCys3UBuq1bdULWpUv1wHffPqhfH8aM0aorGyRYDOnpWh1w8KB2pf35Zz398IMWN3bt0sSxY8fJCwGdfTYZbeNIkQ68vKUbs7/rSlrVs7jgAl3f6cIL9TfapIlO0mhMucvO1qlRvvwSli+HlBRYtw6OHft1m9q1tWtm48Z6atRIx6Oceaae6tTRbWrX1jUPwkBpE0g3YIpzrr/v+t0Azrm/5dpmnm+bL0UkGvgRqA9Myr1t7u0K2l9pEsiezftJefhDsrLIOWVkwIkMOHEc0o/DL0d0QtZDh/T/LSPXDAfnnK3JIi4eOnXMpz6+LBX2OeS+r6jLuc/znrKzfz3Pfcr9hmVl6bQPmZn65vlPx4/rKT1dT8eOac+Co0fh8OFfT8eP5/8aKlXSuYoaNdIf1rnn6pt97rnQtq3+uHK9lOXLtRZh6VJYufLXXFOpkj7FOefob7FOHahVSztd+U+VK+spOlp/o1FR+jiRU09+BV02pjCSnUWNn7Zx+o/fUPOnrZz+0xZq7NlB9QOpVN+3i6pHDxb42IxqNcioVpOMajXJrFaDzCrVyaoSQ1blGLIqVyOrclWyK1cjK7oKLqoy2b6TqxRFdlQ0rlIUTqL0vFIUTipBpUo4qYQTAf955cp0nXl9yV9jIQkkkL/IhsCuXNdTgS4FbeOcyxSRNKCu7/av8jy2YT4BjgPGATRp0iSAkPK3N/lb+r48usSP50ffKVSmHveKiP77Vq6sh/uVK0PVqiefqleHmjX1j79mzV9PtWtrn/szztBpeP1HXvXqBZyRRXTmia5d9fqxY5pEtmzRgsuOHdpJZs8eraJOS9Ntjh0L8dUPTQSKAlr7Tqeqzi/UZw9n8jNn8RNncCDndHr6IWqmH+Z0DlGDI8RwjBgOU52fqMpxqnKcGI5ThRNUJiPnPJrirQN/lBgoRQIpTCC/6PyOx/L+TAvaJpDH4pybCcwELYEEEFO+Wg5uy96vthIdTc6pSpUw6qRT2KFvIIfL/su5z/2n3IfguQ/J/Zf9h+r+6yEkJubXaqzCOKedtjIyfi1AZWWdXNjKXSDL/bj8LhtTeqf5Ts2K/chs4JjvdBJ/LUJWFpKd9et57tv5tdZBBJqW9mUUIJAEkgo0znW9EfB9Aduk+qqwagH7A3xs0FSuUZV6Xc4tq6c3IU7k1wKSMZFL0JKP920ogRxqfg20EpHmIlIFGA7MzbPNXGCM7/KVwOdOG1fmAsNFpKqINAdaASuCE7oxxhgvFVkC8bVp3AzMQ1Pe88659SIyFUhyzs0FngNeEZGtaMljuO+x60XkTWADkAn8qbAeWMYYY8KHDSQ0xhhToMJ6YYVWa6kxxpiwYQnEGGNMiYRcFZaI7AG+LeXT1AP2BiGccGbvgbL3wd4DP3sfSvYeNHXO1c/vjpBLIMEgIkkF1dlVFPYeKHsf7D3ws/ch+O+BVWEZY4wpEUsgxhhjSiRSE8hMrwMIAfYeKHsf7D3ws/chyO9BRLaBGGOMKXuRWgIxxhhTxiIqgYjIABHZLCJbRWSS1/GUFxFpLCILRWSjiKwXkQm+2+uIyKcissV3fobXsZY1EYkSkVUi8l/f9eYistz3Hsz2zecW0USktojMEZFNvu9Et4r2XRCR232/hXUi8oaIVKsI3wUReV5EfhaRdbluy/ezF/Uv3//lGhHpWNz9RUwC8a2c+BQwEGgLjPCtiFgRZAL/55w7H+gK/Mn32icBnznnWgGf+a5HugnAxlzXHwWm+d6DA+jyypHuSeAT59x5QDz6flSY74KINARuBRKdc+3ROfyGUzG+Cy8CA/LcVtBnPxCd4LYVuh7TjOLuLGISCLps7lbn3Hbn3AlgFnCZxzGVC+fcD865lb7Lh9E/jIbo63/Jt9lLwOXeRFg+RKQRMAh41nddgD7AHN8mFeE9OB24CJ3gFOfcCefcQSrYdwGdKDbGt7xEdeAHKsB3wTm3GJ3QNreCPvvLgJed+gqoLSLnFGd/kZRA8ls58ZTVDyOdiDQDOgDLgbOccz+AJhngzIIfGRH+CfwZXYsHdFXMg845/8LFFeE70QLYA7zgq8p7VkROowJ9F5xzu4HHge/QxJEGJFPxvgt+BX32pf7PjKQEEtDqh5FMRGoAbwO3OecOeR1PeRKRS4GfnXPJuW/OZ9NI/05EAx2BGc65DsAvRHB1VX58dfyXAc2BBuiSgAPz2TTSvwtFKfXvI5ISSLmufhhqRKQymjxec86947v5J3+R1Hf+s1fxlYMLgcEishOtvuyDlkhq+6oxoGJ8J1KBVOfcct/1OWhCqUjfhb7ADufcHudcBvAO0J2K913wK+izL/V/ZiQlkEBWToxIvrr+54CNzrknct2Ve6XIMcD75R1beXHO3e2ca+Sca4Z+9p8750YBC9FVMiHC3wMA59yPwC4RaeO76Tfogm4V5ruAVl11FZHqvt+G/z2oUN+FXAr67OcCo329sboCaf6qrkBF1EBCEfktetTpXznxIY9DKhci0gNYAqzl1/r/e9B2kDeBJuiPaqhzLm8DW8QRkV7Anc65S0WkBVoiqQOsAq52zh33Mr6yJiIJaEeCKsB24Fr0YLHCfBdE5K/AMLSH4irgBrR+P6K/CyLyBtALnXX3J+B+4D3y+ex9yXU62mvrKHCtc65Yq/lFVAIxxhhTfiKpCssYY0w5sgRijDGmRCyBGGOMKRFLIMYYY0rEEogxxpgSsQRijDGmRCyBGGOMKRFLIMYYY0rk/wNBXX1GSmneeAAAAABJRU5ErkJggg==\n",
+ "text/plain": [
+ "<Figure size 460.8x216 with 1 Axes>"
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2deZwcVbm/n7dnsicsCQHCEhJ2NBqWEQEVI4sKLqAXRXNlUwyiqIDIol4FrnpVvIIKsguIC1wQIXD9KQJyQUEwQTZFQFliIEBYRZYkM3N+f7zV090z3ZmurqquXr7P59Ofnj59qupNZ+adZ06d8x4LISCEEKL5FPIOQAghuhUlYCGEyAklYCGEyAklYCGEyAklYCGEyAklYCGEyAklYCG6BDN7i5ndn3ccooQSsOh6zGy+mS0ys3+Z2TIz+39m9uaE53zEzPZIK8Y6rhfMbPPV9Qkh3BxC2KrB8z9iZivNbJ1h7XdG157VyHm7HSVg0dWY2dHAacDXgfWAmcAPgH3yjCttzKw3hdM8DHy47JyvAyakcN6uRQlYdC1mtiZwMvCpEMIVIYSXQgirQghXhxA+H/UZZ2anmdnj0eM0MxsXvbeOmV1jZs+b2bNmdrOZFczsYjyRXx1Z9bFVrj3PzJaa2bFm9lRk3vua2d5m9kB0vi+U9d/RzG6NrrXMzE43s7HRezdF3e6Krrd/2fmPM7MngAuKbdExm0XX2D56vYGZPW1m81bzkV0MHFj2+iDgRw19+AJQAhbdzc7AeOAXq+nzRWAnYFtgLrAj8KXovc8BS4HpuD1/AQghhAOAJcB7QgiTQwjfqnHu9aPrbwh8GTgX+AiwA/AW4MtmtmnUdwA4Clgnint34JP4BXeN+syNrndp2fmnApsAC8ovHEL4O3Ac8BMzmwhcAFwYQrhxNZ/FH4A1zGwbM+sB9gd+vJr+YhSUgEU3Mw14OoTQv5o+/w6cHEJ4KoSwHDgJOCB6bxUwA9gkMuebQ7ziKquAr4UQVgGX4Mn1uyGEF0MIfwb+DLweIISwOITwhxBCfwjhEeBs4K2jnH8Q+EoIYUUI4ZXhb4YQzgUeBG6L/h1frCPmogXvCfwVeKyOY0QNlIBFN/MMsM4o46MbAI+WvX40agM4BfgbcK2ZPWRmx8e9fghhIPq6mCCfLHv/FWAygJltGQ13PGFm/8THrCtuiFVheQjh1VH6nAvMAb4fQlhRR8wXA/OBg9HwQ2KUgEU3cyvwKrDvavo8jv8JX2Rm1EZkqp8LIWwKvAc42sx2j/qlXWbwTNw4twghrIEPd9gox6w2BjObjN+APB840cymjhZECOFR/Gbc3sAVdcQtVoMSsOhaQggv4GOvZ0Q3wCaa2Rgz28vMiuO2PwO+ZGbToylYXyYa9zSzd5vZ5mZmwD/xcdqi0T4JbEp6TImu8S8z2xo4fNj7jVzvu8DiEMKhwP8CZ9V53MeA3UIIL8W8nhiGErDoakII3wGOxm+sLQf+ARwBXBl1+SqwCLgbuAe4I2oD2AK4DvgXbtM/KLuJ9V944n7ezI5JIdRj8D/9X8SHDS4d9v6JwEXR9T442snMbB/gncAnoqajge3N7N9HOzaE8PcQwqIYsYsamAqyCyFEPsiAhRAiJ5SAhRAiJ5SAhRAiJ5SAhRAiJ9Io0CE6gHXWWSfMmjUr7zCEaCsWL178dAhheqPHKwELAGbNmsWiRZpZJEQczOzR0XvVRkMQQgiRE0rAQnQbL70Emv/fEigBC9ENPPggHHggbLYZTJ4M06bB298OV1yhZJwjSsBCdDL9/XDiiTBnDlx5JWy3HZx8Muy3Hzz8MPzbv8G73gVLluQdaVeim3BCdCr9/XDAAXDJJfDhD8N//zfMmFH5/umnw3/8B7zlLXDjjTB7dm7hdiMyYCE6kYEBH3K45BL4xjfgpz+tTL4Avb1w5JFw003w4oswb55bsWgaSsBCdCInnQQ/+5kn3+OOW33f7baD66/3JLzPPvDKiM0zREYoAQvRadxwA3z1q3DIIaMn3yLbbeeWfM89cPTR2cYnhlACFqKTWL4cPvIR2HJL+P734x37znfC5z8PZ50Fl1+eTXyiAiVgITqJY4+Fp5+GSy+FSZPiH//Vr0JfHxxxBLzwQvrxiQqUgIXoFG69FS68EI46CubObewcY8fCmWfCU0/5OLLIFCVgITqBgQG31g028GllSejrgwUL4Hvfg3vvTSc+URUlYCE6gQsvhDvu8Lm+kycnP9/XvgZrruk2LTJDCViIdmfFCl/dtuOOsP/+6Zxz2jT44hfhuuvg//4vnXOKESgBC9HunH++LyX+z/8Es/TOe/jhvnjjP/5D9SIyQglYiHbmlVd8uODNb4Y990z33BMmwBe+ADff7CYsUkcJWIh25rzz4PHH07ffIh//OGy8MXzlK7LgDFACFqJdWbUKvv1tt99587K5xrhxvpru1lvh97/P5hpdjBKwEO3KZZf52O+xx2Z7nUMO8Zty3/pWttfpQpSAhWhHQvCEuM02Xs83SyZO9DnGV18N992X7bW6DCVgIdqR3/wG7rrLazcUmvBj/KlP+U25U07J/lpdhBKwEO3IqafC+uvD/PnNud706T4U8ZOf+DJlkQpKwEK0Gw88AL/6lc/THTeuedf99Kdh5Uo455zmXbPDUQIWot04/XQYM8brNTSTrbf2jTzPPNNnYIjEKAEL0U68+KLXfdh/fx+CaDaf+YzPO77iiuZfuwNRAhainbjoIk/Cn/50Ptffay/f2v5738vn+h2GErAQ7UII8IMfwBve4IV38qBQgE9+Em65xWdhiEQoAQvRLtx0k8/DPfzwfOM4+GAYP963LhKJUAIWol046yxYa630Sk42ytSpHsOPf+zDIaJhlICFaAeefBJ+/nO3z4kT847GLfxf//J5waJhlICFaAd++EOf+nXYYXlH4uy4I2y7rU9JU5W0hlECFqLVGRyEc8/1imdbb513NI4ZfOITcPfdcNtteUfTtigBC9HqXHcdPPxw69hvkfnzYdIkrYxLgBKwEK3OOefAOuvA+96XdySVTJniSfiSS+D55/OOpi1RAhailXniCbjqKr/51sy6D/Vy2GG+LZJuxjWEErAQrcwFF0B/v28N1IrssANsvz2cfbZuxjWAErAQrUr5zbctt8w7mtosWAD33KObcQ2gBCxEq9KqN9+Go5txDaMELESrcvbZrXnzbTjlN+NeeCHvaNoKJWAhWpFly2Dhwta9+TacBQt0M64BlICFaEVa/ebbcPr6dDOuAZSAhWg1Bgb85tvb3tbaN9+Gs2CBr4z7wx/yjqRtUAIWotW49lp45BFf6ttOzJ8Pkye7BYu6UAIWotU46yxYbz3Yd9+8I4nHlClwwAFw6aXw7LN5R9MWKAEL0UosWQLXXAMf+xiMHZt3NPE57DB49VXfOkmMihKwEK3Eeef5Tax2ufk2nLlzYeed3eJ1M25UlICFaBVWrvSbb3vtBbNm5R1N4xx+ODzwAFx/fd6RtDxKwEK0Cldc4cV3jjgi70iS8YEP+AKS00/PO5KWRwlYiFbh9NNh883hHe/IO5JkjB/vU9Kuvtpnc4iaKAEL0Qr86U/w+9/Dpz7lW7+3O8UpdNo5ebV0wP+0EB3A6af7ZpsHH5x3JOmw8cY+je6883yJsqiKErAQebN8uddQOOAA33a+U/j0p+GZZ1QfYjUoAQuRN2eeCStWwJFH5h1Jurz1rb5z8qmnakpaDZSAhciTV1+FM86Ad72rdXY8TgszOPpo+Mtf4Ne/zjualkQJWIg8+clP4KmnPFF1IvvvDxtsAN/5Tt6RtCRKwELkxeCgJ6Ztt/XKZ53I2LE+Fvyb38Bdd+UdTcuhBCxEXlxzjf95/rnP+Z/rncphh3mVtG9+M+9IWg4lYCHyIAT42tdg9mz40IfyjiZb1l7blydfein87W95R9NSKAELkQc33AC33w7HHQe9vXlHkz1HHQVjxsC3vpV3JC2FErAQefD1r8OMGXDQQXlH0hxmzICPfhQuvBCWLs07mpZBCViIZnPTTW7AxxzjdRO6hWOP9aGX//qvvCNpGZSAhWgmIcCXvuRGePjheUfTXGbN8kLz556rIj0RSsBCNJNrr4Wbb/YkPGFC3tE0ny99yYsNnXRS3pG0BErAQjSLov1usgkcemje0eTDRhvBJz8JP/oR/PWveUeTO0rAQjSLSy6BRYvgxBPbc7+3tDjhBJg0yWeAdDlKwEI0g5df9oSz3XZw4IF5R5Mv06fDF78ICxfCddflHU2uKAEL0Qy+/W34xz/gtNM6o+B6Uj77WV+EctRR0N+fdzS5oe8EIbJmyRJfhrvffrDrrnlH0xqMHw+nnAL33gtnn513NLmhBCxEloRQmm52yin5xtJqvP/9sMcePib82GN5R5MLSsBCZMmll8Ivf+l1H9p5q/ksMPM94/r7fWZEFxZtVwIWIiuefho+8xl4wxu8JKMYyWab+ZzghQvhssvyjqbpKAELkQUheO2DF17wjSl7evKOqHU56ijo6/OdlP/xj7yjaSpKwEJkwRlnwNVXe/Wv178+72ham95e+OlPYeVK+MhHYGAg74iahhKwEGlzxx1eaGfvvX0IQozOFlv4L62bboKTT847mqahBCxEmixbBvvsA+uuCxdc0Nk7XaTNgQd6ec6TT4bLL887mqbQBZWghWgSr74K73sfPPss/P73noRF/RRnRTzwgCfjTTeF7bfPO6pMkQELkQYrV/pCi9tug4sv9o02RXzGj4df/MKXK++1V8cX7FECFiIp/f0wfz787//CmWf6AgPROOut52U7zWD33eHvf887osxQAhYiCS+95An35z+HU0/1qVQiOVtt5YV6VqyAt7wF7rwz74gyQQlYiEZ58knYbTc33zPOgCOPzDuizmLOHLjxRp9Dveuu8Otf5x1R6igBC9EI110Hc+fCPff4mOUnP5l3RJ3JnDlw661eOW2vveDLX+6o6mlKwELE4cUX4eij4e1vh6lT/abbe9+bd1SdzUYb+aySAw+E//xPmDcP/vznvKNKBSVgIephYAB+/GPYZhsf6z3sMPjjH+F1r8s7su5g8mTf0v7ii+G++3yWyTHHeL2NNkYJWIjV8fLLvqDita+FAw7wub233uqzHSZNyju67uMjH4H77/f/i+98x4cmjj8eHn0078gaQglYiOH098Nvfwuf+hRsuKEX1Rk71ldnLVoEO+2Ud4TdzTrrwA9/6MXc3/Uur7O86ab+9Y9+BM8/n3eEdWOhC2twipH09fWFRYsW5R1GPrz4ok9z+uMfvRbBjTd6FbOJE31Z8Sc+4VOhtKy4NVmyBM49Fy66yKup9fR4CdC3vQ3e+EavtLbBBpn8/5nZ4hBCX8PHKwEL6LAEPDAAr7ziwwcvvuiP557zJcJPPeXTx5Ys8T9bH3gAli4tHTt7tk/+f+c7/aFhhvYhBP8lunAh3HAD3H57qbLammv63OLZs2HmTJgxwxd8TJsGa6/t70+Z4v/fEyd6hbY6ErYSsEiFvokTw6LNN2/eBYd/3xVfl7eHMPIxOOiPgYHSo7/fH6tW+cT90aYpmfkP4MyZXoVr6629ZGRfH6y/frr/TpEfL78Md90Fixf7kub77/dfukuW+PfJ6jCDceN86GnMGE/IPT2lx6abwvXXJ07AKsYjnHHjoJkJGEYaRvF1ebtZ5aOnp/Rc/hgzxh/jx/u/ZcIEf0yZ4o+11/ZpY+uu63UGevWt3/FMnAg77+yPckLwceInn4RnnvG/jP75T/9L6eWX/a+nV1/1JL1ypf9i7+8v/cIfHEyt0JK+C4Wz2WZwxRV5RyFE9pj5L+S11847Es2CEEKIvFACFkKInNBNOAGAmb0I3J93HKOwDtDqS5/aIUZojzjbIcatQghTGj1YY8CiyP1J7uY2AzNbpBjToR3ibJcYkxyvIQghhMgJJWAhhMgJJWBR5Jy8A6gDxZge7RBnx8eom3BCCJETMmAhhMgJJWAhhMgJJeAux8zeaWb3m9nfzOz4vOMpYmYbm9lvzew+M/uzmX02ap9qZr8xswej59zXk5pZj5n9ycyuiV7PNrPbohgvNbOxOce3lpldbmZ/jT7PnVvtczSzo6L/53vN7GdmNr4VPkcz+6GZPWVm95a1Vf3szPle9LN0t5ltP9r5lYC7GDPrAc4A9gJeA3zYzF6Tb1RD9AOfCyFsA+wEfCqK7Xjg+hDCFsD10eu8+SxwX9nrbwKnRjE+B3wsl6hKfBf4VQhha2AuHmvLfI5mtiHwGaAvhDAH6AE+RGt8jhcC7xzWVuuz2wvYInosAM4c9ewhBD269AHsDPy67PUJwAl5x1Uj1quAPfHVejOithn4ApI849oo+iHcDbgGMHz1Vm+1zziH+NYAHia64V7W3jKfI7Ah8A9gKr447BrgHa3yOQKzgHtH++yAs4EPV+tX6yED7m6K3/hFlkZtLYWZzQK2A24D1gshLAOIntOpC9g4pwHHAoPR62nA8yGEYlHivD/TTYHlwAXRMMl5ZjaJFvocQwiPAd8GlgDLgBeAxbTW51hOrc8u9s+TEnB3U63kf0vNSzSzycDPgSNDCP/MO55yzOzdwFMhhMXlzVW65vmZ9gLbA2eGELYDXqI1hm2GiMZQ9wFmAxsAk/A/54fTUt+bVYj9f68E3N0sBTYue70R8HhOsYzAzMbgyfcnIYRiseInzWxG9P4M4Km84gPeBLzXzB4BLsGHIU4D1jKzYp2VvD/TpcDSEMJt0evL8YTcSp/jHsDDIYTlIYRVwBXALrTW51hOrc8u9s+TEnB380dgi+hu81j8xsfCnGMC/I4ycD5wXwjhO2VvLQQOir4+CB8bzoUQwgkhhI1CCLPwz+6GEMK/A78F9ou65R3jE8A/zGyrqGl34C+00OeIDz3sZGYTo//3Yowt8zkOo9ZntxA4MJoNsRPwQnGooiZ5Dbzr0RoPYG/gAeDvwBfzjqcsrjfjf77dDdwZPfbGx1ivBx6MnqfmHWsU7zzgmujrTYHbgb8BlwHjco5tW2BR9FleCazdap8jcBLwV+Be4GJgXCt8jsDP8HHpVbjhfqzWZ4cPQZwR/Szdg8/qWO35tRRZCCFyItEQRKtO4hdCiHagYQOOJvE/gM/NXIqPJ344hPCX9MITQojOJcmOGDsCfwshPARgZpfgU0lqJuB11lknzJo1K8ElRVKWL4fnnoMtt6xsX/uO2fFPNnxbeQArDHs5fOv56P1h7VY8V6FQee7oden98i3rh52r0FPxeuiYnp7Kcw61F4/351AYFkNPZSyhUPZv6xnWNvQ6eo6ODT1W/f3oebB3+HHFZyqPB0J0qcFhfYa/LvYrvV95rsHeyvdHPJf9M0f07Q01zh0q34+eidqJXlv02np92nRPT/TcW5xGDWPG+NTfMT0DAIzt9edxPcVnf39C7yoAxkfPk3pWentP9Lp3BQATC94+pedVACZHz2sUXhm65sTCiqjN35sy9OznmmIheu0fyOTCeAAK6z9YbepZ3SRJwNUmHb9xeCczW4Avy2PmzJksWpRoBw+RkM9/Hs44A4b/N+xZ+ED8k5X/9VRMaiH6QYqSYxiMfuAKw94frEyexb/EbDB6v5jYotfFRGeln1MoDDsXA9FzT3RMFMpA1F5MxEUGBiteWjQiF6hsLybiodiA4homi/qWXhcp9i2ek2HvR+9GywwGR/wkFnuGEW2FqG2wxuvhFD+dwahfIeo3WLV3jfhqxFU6d/VrD//7uvTajxygGinvlJbG6aJEzGBxXUiUxJOeNsGxdU06DiGcE0LoCyH0TZ8+PcHlRBr098OYMXlHIYSAZL8bWnoSv6jOihUwNouaUkUbbpIJe5/oiyabcHl8zTPhkUfLhGOSgQnnacAtO4lf1Obll2HixLyjEEJAgt8JIYR+MzsC+DWuGj8MIfw5tchEJvzrXzB5coYXaJIJQ5Vx4WaZMIwYF87ehMt7d64JQzUbbmETTkiiUEIIvwR+mUokoik88wxMm5Z3FEIISP1Xi2h1nngCXvvaJlwoaxOG2jMkZMJVXw+nNU24dHRbmHBCVIynixgchIcfhtkNTPkVQqRPC/wOEM3ioYd8FsTwRRiZkpUJw+hzhTMyYRh9rnDaJgz1zBVO14TLo6xF2iZc2dYkE87olPUgA+4ifvc7f95ll3zjEEI4MuAuYuFCWHdd2GabHC6esglDjFVzKZsw1L9qLi0Thjir5tIxYW+rb1xYJtwYMuAuYdkyT8AHHVR5H0sIkR8y4C7hpJNcQj/+8RodzCprO2RFSibsp4pZP0ImXPX1iPOXfR13hkRSEy6PfuTrzjNhuVAX8H//B2efDUceCVtskXc0QogiSsAdzp//DB/4AGy6KZx8ct7RCCHK0RBEB3PvvbDbbtDbC7/8JUyaNMoBQ8MCbTAUAY2Xskw6FAENl7JsdCiisk/UI+OhiFIUGorIChlwBzI4CD/4Aey0kyffG2+ErbYa9TAhRJORAXcYDz8Mhx4KN9wA73gHnHcebLRRHQdaocxC28CEIXlR9zYyYe/DsD5RD5nwMNrHhGXAHcKDD/oMh622gj/+Ec49F/7f/6sz+QohckEG3ObceSd84xtw2WW+08XHPw7HHw8bbzz6sSMo7rHWDiZc3qfZJgypbW9UvwlDetsbta4Jr+76nWjCSsBtyNKlcMkl8NOfwp/+BFOm+F5vRx4J66+fd3RCiHpRAm4Tnn0WLr/ck+5NN7kwvuENcOqpcPDBsNZayc5vBRuyzHYwYW9KZ3ujuCYMKZSyjG3CkP5Gn6s34fK24WRlwpXn7nwTVgJuUfr7fefiX/8arr0WbrsNBgZ8jPfEE2H+fNh887yjFEIkQQm4hXj0UU+2114L110Hzz/vItjX5+O6738/bLddSQ7TpmiU7WDC3pTuRp/1mrAfE4XTJBP2tiLNMeFabRXXGIooLRMuxdUWJpyQ1ouoS+jvh3vugVtu8cett/oUMvCZC+9/v08j2313bSEkRKeiBNwknnsO/vCHUsK97TZ46SV/b8YMeNOb4LOfhbe/HbbeOjvLrUnZPOC2MOGyuJpvwpDVlve1TLiyrUi2JuxXyHZ7o5EmPDKuTjbh1omkg3j+eZ+dsHgx3HGHPz/wgL/X0wNz58Ihh3hh9F12gZkzc0i4QojcUQJOyLPPlpJs8fnvfy+9v/HGsMMOcOCBnmzf8IaMt4VvlIKVTC+uCUP2NjzchCviaLYJQ2YbfdYwYT8mWSnLuCbsZ2zORp8VWy/ViKsTTTj/CNqE4oaWd91V+XjkkVKfWbM82X70o/68/fYwfXpeEQshWh0l4Cq89JLfICtPtPfcAy++6O8XCr6x5RvfCIcf7ol2++1h6tR8405M0RzjmjA0b1y4/PxZbXk/mgmXn6tJJux9Gquk1rgJl9o6Y8v7WiYMeaXCrk7AIfiqsuFW++CDpZ/zNdaA17/ehxDmzoVtt4XXvhYmTsw3diFE+9M1CXjFCvjLX0Ym22efLfWZPduT7Pz5/jx3rg8rdMMNMjMrbXgZ14TL+rTEDImsTRhS395oNBMuj695JjzyaJlwunRkAn75ZS9Ss3ixrya74w7461997i3AhAnwutfBv/1bKdG+/vVuu0II0SzaPgG/8oqbbDHZLlrkpluUmfXW8xti731vKdluvnn10q5dTaFQMqy4JuyNFX062YS9T/RFs0wYMtvos7YJl/fuXBOG/GZItF0Cfvxx32Tyxht9McO993qNBPAZB3198L73+fMOO8AGG3THEIIQov1o+QS8bJkn2+KjuKBhjTV8y513v9sTbV+fL+FVsm0QsyHji2vC3tTCq+ZSNmGIXz9CJlweQauZcOnoZptwyyXgEHxHhx/9yAvS3H+/t6+xBuy6KyxYAPPm+WwEDSMIIdqZlknATz8NF18MP/yhDytMmOA7+h56aCnh9rZMtB1IoVAyvJgmDG1SPyItE4bs9pmrYcIQv35EUhOG7HbXqGXC5VHWIm0TrmxrrgmPelYz2xj4EbA+/m8+J4TwXTObClwKzAIeAT4YQngubgBLlsAxx8CVV8KqVbDjjnD22bD//rDmmnHPJoQQ7UM9ab0f+FwI4Q4zmwIsNrPfAAcD14cQvmFmxwPHA8fFufjTT8Oee/o47xFH+BLeOXPi/hNEGphZ6S5/XBOG+PUj2tmEIbXdNeo1YYhfPyKpCUN6u2vUa8Le1mgltfYz4VHPFkJYBiyLvn7RzO4DNgT2AeZF3S4CbiRGAl6xwm+gLVniY71velPMyIUQos2Jlc7NbBawHXAbsF6UnAkhLDOzdWscswBYADBz5syh9iee8Dm7fX1eIUzkTMGG7C22CUPjldTa0IShgfoRMuGK1yPOX/Z1Vjsur26357wqqRVG7+KY2WTg58CRIYR/1ntcCOGcEEJfCKFvellpsE02gXPO8bm8Bx/swxFCCNFN1JWAzWwMnnx/EkK4Imp+0sxmRO/PAJ6Ke/GPfhS+/nX42c98wcQHPwi/+lVpYYUQQnQy9cyCMOB84L4QwnfK3loIHAR8I3q+qpEATjgB3vMeOP98n4Z22WW+oOLgg2HffTXft2lYYejP5dhDEdB4Kcs2HIrwUyUt6h5zKAIy3PK++lBEZZ+oR8ZDEaUo2mUoIhn1GPCbgAOA3czszuixN5549zSzB4E9o9cNMWcOnHoqPPaYJ+A5c+BrX/Px4WnTPEH/9397UR3ZsRCiU6hnFsTvqD4uD7B7msGMGwf77eeP4hLk3/7Wn6+5xvusuaaviJs3z2dOzJ0L48enGUWXUr4lUUwThvqXLXeECUPyou5tYMLeh2F9oh4y4VRo2bVlM2bAhz/sD3A7LhbhufFGuPpqb+/tdWMu1oPo6/NSk+PG5RW5EELUR8sm4OFsuKEXSp8/318/9hjcfnupBOUvfuHjyABjxngSLibk7bf3XSxkyquh0MPQ7/WYJux94hXwaWsThvS2N6rXhCGzLe9rmzAkL+re+ia8uuvXU8oyCW2TgIez4YZedvJ97/PXIcCjj5YS8uLF8D//41PdwL93t9qqVBO4uL3Q+uvn928QQnQ3bZuAh2Pm2wfNmuVjyOBJ+aGH/OZdcQui3/3Op70VWXfdyqQ8dy5svduwYmoAABYhSURBVLVbdFdRMIr2FduEoeFSlm1pwuV9mmTCEL+AT3IThvS2N6rPhMvbhpOVCVeeu5FSlo3TMQm4Gmaw2Wb++MAHSu3PPgt33125N9z3v+/LowHGjoXXvGZkYp42LZ9/hxCiM+noBFyLqVN9FsW8eaW2Vau82Ht5Uv7Vr+Cii0p9NtywNHRRTMpbbFFZpbBd8WI8xVcxTRgaLmXZjibsTQ2WsmzQhP2Y6NpNMmFvK9IcE67VVnGNoYjSMuFSXI0U8ElCVybgaowZ4zfqXvva0o0+gCefHLmT8rXXljb4nDzZE/IOO/jNvh128LFm1S4WQoyG0sQorLcevP3t/ihS3OL+zjvhT3/yG37nnuu7MYMXk99221JC3mEH2GabFh9X7ukZMqu4JuzHxCzg08Ym7E0Ji7rHNmHIasv7WiZc2VYkWxP2KzRWyrJxEx4ZV7NMWAm4AcaNg+2288chh3jbwIBvn7R4sd/0W7zYhy/OOMPfnzDBi83vsos/dt5ZY8pCdDtKwCnR0+M37l7zGjjgAG8bHIQHH/RkfPvtcMstcMoppeGLrbYqJeRddvHZF7mNJ5sNmVRcE4YG6kckNWE/Wf3/vkaoZcJlcTXPhCGzjT5rmLAfE69+RFIT9jMmK+oe14Qr+8ZfNZcEJeAMKRQ8yW61VWlc+eWXfZ7yLbf4Y+FCuOACf2/aNNhjj9KQx0Yb5Re7ECJ7lICbzMSJXsti1139dQhuybfe6nUvrr0WLr3U33vNa0rJeNddYdKkDAMrN8u4JgyNV1Jr1ITLY262CZfH0SwTLj9Xk0zY+8SrH5HchEtt7bTlfaN0wASq9sYMttwSDjoILrzQl1jffTd8+9tuwGedBXvv7VPn9tnHk3PxZp8Qor2RAbcYZl7H4nWvg899Dl55BW6+GX75Sy/VuXChT33bd18f1thjj5RmV/QURm5/U68JQ8OV1Bo24bI+TTfhims2yYQheSW1mCZcHl/zTHjk0Z1swjLgFmfCBB+COO0038D0hhvgQx/y8px77+2LQ048EZ55Ju9IhRBxUQJuI3p64G1v8znHTzwBV17p09lOOsn32Dv6aB/CaIhCwc2np+B2V/7o6fF5wmaYmVtcwbyC2tAjarOCP6LXpWMK/iieM3o99P5QHMPOE2EFq6zF4I2VRlw8dzMIodKIw2DF+HQYDBUr54beHwz+GDpNcBseHPRH8bzR6+L73id6DD/X4ED08NdD/QcG/FE8Z/ExMOiP6Bo2GLDBshgGyh7RMTY46DY8EGCg2uvoMTAYPQIWvVfxfvQo9PujdFz5g+gRvR70vwQKA4FC2fvDXxf7ld73R/E8hX433NL5qzyK1xret9/8MezcSVECblPGjfMx4auugnvu8apw3/sezJ7tiVjjxEK0PkrAHcCcOb6f3oMP+s28U0/1lXi33FL/OUKhzFJjm7A13YQrbLgLTLjchptmwoM5mPBge5lwUpSAO4jZs3144oYbvLjQm9/s48PNWL0rhIiPZkF0IG97m09lO+IIHx9++WX45jdHkcJCYehuuA3/vTzK7AiIXz8i6ewIqGOucAvVj0g8OwIS1xSOOzsC4tePSDo7ApLXFI47O6I8ylqsvpJa4ygBdyhTpvi84kmTfPnz2LHw1a/mHZUQohwl4A7GzIsBvfoqfP3rsNdevpN0VXpGWk+9Jgx1zBVO24Qhfv2IdjZhSFBJrTEThtHnCqdtwtB4JbVGTdjbklRSaxyNAXc4Zj47YuZMOPRQWLky74iEEEVkwF3A5Mnw3e/66rkrr4QPfrBKJ7MKC4YYJgwNV1Jr2ITLz9UFJgxVxoVlwhWvGzVh79N4JbUkyIC7hPe8xxdrnH123pEIIYooAXcJhQIceKBXXHvhhbyjEUKAEnBX8da3+l/Gt9468r1QvjiiuBAjWiQRClbfQo26li2ntFAD6l+23AELNfxU9S1bTm2hRpxlyykt1Ii3bDmdhRoNLdYYLBvySYAScBex447+fMcd+cYhhHB0E66LmDLFNxl9+OEqb/ZY6YZJ8WZPvTfloPFSlg3elIPRF2t01E05GH2xRlY35SD1Le9r3ZTzPgzrE/XI+KZceRSNlLJsBBlwl7HJJl7WUgiRPzLgLmPaNFi+fGR7KBRGWkoLm3B5fF1hwpC8qHtcE4bMtryvbcJQ77LlTjBhGXCXsfba8NxzeUchhIAYBmxmPcAi4LEQwrvNbDZwCTAVuAM4IISgdVYtzsSJNWoFl48BxzRhP6a+ZcupmTDELuDT1iZc3qdJJgz1L9ZIz4QhbgGfpCZc3jacrE04jgF/Friv7PU3gVNDCFsAzwEfSykmkSETJvg+c0KI/KnLgM1sI+BdwNeAo80nZO4GzI+6XAScCJyZQYwiRcaMgf4qyyhDwcq8IZ4Je586ly2nZcLQcCnLdjRhb2qwlGWDJuzHRNdukgl7W5HmmHCttoprDEU0spRlEuo14NOAYyl9ItOA50MIxR/lpcCGqUQkMqW314u1CyHyZ1QDNrN3A0+FEBab2bxic5WuVXXBzBYACwBmzpzZYJgiLXp6KodRi4SeAgxZbNTWwibsx6x+rnAnmbA3JSzqHtuEIast72uZcGVbkWxN2K/QeCnLJNRjwG8C3mtmj+A33XbDjXgtMyv+UzcCHq92cAjhnBBCXwihb/r06SmELJJQKPiqUyFE/oxqwCGEE4ATACIDPiaE8O9mdhmwH56UDwKuyjBOkRJm1YUu9Bjlv98hhglD46UsGzRhv350rW4w4bK4mmfCkNb2RvWasB9T36q5tEzYz9h4KcskJDnPcfgNub/hY8LnpxOSyJJaCVgI0XxirYQLIdwI3Bh9/RCwY/ohiSypWQSsxyo2aHHqM2FoYNVcUhOG5EXdGzVhPxmZMtyEK+JokgmXn6tJJux96l01l5YJl9oaWTWXBK2EE0KInFAtiC6jlgFXzgMuUp8JV7Q1y4Qh/Y0+6zVhaN64cPn5s9ryvpYJQ/JKajFNuDy+5pnwyKObZcIyYCGEyAkZsADcgIvENWFvq2+ucHomDJlteT+aCZf1aYkZEhmZsPeJvmiWCUNmG33WNuHy3s01YRmwEELkhAxYADDYayO22q7XhL1PvFVzyU0YMtvyfjQT9saKPp1owhBn1Vx3mnBSZMBCCJETMmAB+Bhw0QjimnBln+aYsMc89GbxoCisbE3Ym1p41VxaJgzZ7TNXw4Sh/lVzaZkwJKuklgQZsBBC5IQMWABEtSCc+CYMSSupxTVhiLFqLmUThhir5trZhCH9feZGMWGof9VcWiYMySqpJUEGLIQQOSEDFgCEHhjuBvWbMDRcSa1RE4bMdlwe1YSh/lVzbWzCEGPVXJeacFJkwEIIkRMyYAEUx4Cr/3YfzYSBmjMkOtKEy8/VwSbsp8pmx+WaJgwZ7rhc3YQr+0Q9YlVSaxwZsBBC5IQSsBBC5ISGIARQ/NNw9Tcaag1FVDsy86EIyG7L+9GGIqC9tzeqdygCMtjos/WGIrwPw/pEPeoqZdk4MmAhhMgJGbAAYLDHyiaXy4T9/eomDPUvW25rE4YMt7yvYcKQ2Zb3tU0YkhXwaRwZsBBC5IQMWAC+EGNkqb36TBjiF/BJasJ+TH3LltM24fL4OtqEy/s0yYSh/sUa6ZkwpFPKMj4yYCGEyAkZsAAqx4DjmjA0UsAnmQl7nzoXa6RtwpD+Rp8taMLelM72RvWasB8TXbtJJuxtReKZcFJkwEIIkRMyYAFUHwOWCdcwYchuy/sWMmFvymjL+5omDFlteV/LhCvbitRnwkmRAQshRE7IgAVQeRc6vgmXt0XnyNqEIcMt71dvwn5Mnavm2tmEy+JqnglDZht91jBhP6a+VXNpm7AMWAghckIGLACiguyV1G/CEHfVXFIThgZWzaVkwn7t6FrNMmE/GZky3IQr4miSCZefq0km7H3qXTVXacJJkQELIUROyIAFAIM9tX8bj27CUO+qubRMuKKt2SYMqW/0OaoJQ/PGhcvPn9WW97VMGFLf3mg0Ey6Pr5FKakmQAQshRE7IgAUAoccYjDw0rgmXtzXLhL2tvrnCqZswZLblfU0TLuvTEjMkMjJh7xN90SwThoSV1BqnLgM2s7XM7HIz+6uZ3WdmO5vZVDP7jZk9GD2vnUpEQgjRJdRrwN8FfhVC2M/MxgITgS8A14cQvmFmxwPHA8dlFKfImMFeKES/3+OacPW2bE3Y+8RbNZeeCUNmW97XMmFvrOjTiSYMcVbNtYIJJ2NUAzazNYBdgfMBQggrQwjPA/sAF0XdLgL2TSkmIYToCuox4E2B5cAFZjYXWAx8FlgvhLAMIISwzMzWrXawmS0AFgDMnDkzlaBF+ngtCCeuCXufxupHNGrClX2abcKQ2Zb3NUzYm1p41VxaJgzZ7TNXw4Sh/lVz1eajJ6Ges/QC2wNnhhC2A17ChxvqIoRwTgihL4TQN3369AbDFEKIzqMeA14KLA0h3Ba9vhxPwE+a2YzIfmcAT2UVpMie8pVwcU3Y+ySrpBbfhEuRNtuEof5Vc2mZMMRYNdfOJgzp7zM3iglD/avmqs1HT8KoBhxCeAL4h5ltFTXtDvwFWAgcFLUdBFyVSkRCCNEl1DsL4tPAT6IZEA8Bh+DJ+3/M7GPAEuAD2YQomkH1WhBOa5pweSTNNWGIsWouLROG+lfNtbEJQ4xVcy1gwkmpKwGHEO4E+qq8tXsqUQghRBeilXACGH6Hv5JWNOHyI5tuwpDZjss1Tbj8XB1swn6qbHZcrmnCkKiSWhLS8WghhBCxUQIWQoic0BCEAKKlyKNstV1rKMLbah2TzVBEtSM7eigC2nt7o3qHIiDF7Y2yH4pIigxYCCFyQgYsgGFLkWOasLc1VsqyLU0YstvyvoYJQ/3LltvahCGDjT5HMWFIVsAnATJgIYTICRmwACD0Bhga23XqNWFovJRloyZcLb5ONuHy+DrahMv7NMmEIX4Bn7RMWAYshBA5IQMWQHEpcqV11mvCFX2bZMKQ/kaf9ZqwH1PfsuXUTBjS296ohU3Ym9LZ3qheE/Zjoms3UMoyCTJgIYTICRmwABi2Lb1MGGqbsPeJV8AnsQlD+ht9tqAJe1NGW97XNGFIUsAnCTJgIYTICRmwACD0hDL7HGqNnlvRhMvbonN0sAn7MTEL+LSjCZfF1TwThjRKWTaCDFgIIXJCBiyA4jxgpz1MeGRcQ+fI2oQhwy3vq5uwXz+6VrNM2E9Gpgw34Yo4mmTC5edqpH5EAmTAQgiREzJgAVQacJF6TRgar6TWuAmX4mi2CUMDq+aSmjAkL+oe14SheePC5efPasv7WiYM6VRSawAZsBBC5IQMWDg9gVqOM5oJex+nWSZc3tZsE65oa5YJQ3rbG9VrwmV9WmKGREYm7H2iLxqppJYAGbAQQuSEDFg4ZWPA8U0Y4s6QSGrC1duaY8LeFq9+RHIThsy2vK9lwt5Y0acTTRji149Y3Sa2cZABCyFETsiABQBWZQy4fhMu790cE/Y+6e6uUa8Je594q+aSmzBktuV9DRP2phZeNZeWCUOiSmpJkAELIUROyIAFANY7SPH3cVwThmrjwtmacPn1m23ClX2aY8IQv35EUhOGNqkfkdSEIWEltcaRAQshRE7IgAUAPT2DZb/T45kwpFc/oj1MuBRps0wY4tePSGzCEL9+RBuaMDRQP2IwHXeVAQshRE7IgAUAPb0l85IJV563WiXitGoK123CkNmOyzVNuPxcHWzCfqoEldQSIAMWQoicUAIWQoicqGsIwsyOAg7F/wq7BzgEmAFcAkwF7gAOCCGszChOkTFjxvQz/NuhlYciKs9d69rZDEWUH9nRQxGQvKh7OwxFQLJSlgkY1YDNbEPgM0BfCGEO/r/8IeCbwKkhhC2A54CPpROSEEJ0B/XehOsFJpjZKmAisAzYDZgfvX8RcCJwZtoBiuYwpqd8Ynk8E67VBtmZsLeltb1RPBOudmTmJgzZbXlfw4ShgQI+7WjCkKiUZRJGPUsI4THg28ASPPG+ACwGng8hFL81lwIbVjvezBaY2SIzW7R8+fJUghZCiE5gVAM2s7WBfYDZwPPAZcBeVbpW/RUXQjgHOAegr6+vCb8GRSOM7a22tLJ1TdjPlayou0w4ulYNEy6Pr6NNuLxPA6Usk1CPR+8BPBxCWB5CWAVcAewCrGVmxZ/QjYDHU4lICCG6hHrGgJcAO5nZROAVYHdgEfBbYD98JsRBwFVZBSmyZ1zP6oqLrN6Eof4ZEu2x5f3qTbhafFmbsB8Tr4BPYhOG5EXd28CEvSlBKcsE1DMGfBtwOT7V7J7omHOA44CjzexvwDTg/FQiEkKILqGuWRAhhK8AXxnW/BCwY+oRiVwY11OHctY0YWi0lGU7mjCkt71RvSbsfRorZdmwCUN62xu1sAl7U4JSlgnQSjghhMgJFeMRAEzoXRWj98hvm0ZXzcmEW9eE/ZiYBXza0YTL4mqklGUSZMBCCJETMmABwPhYBlykW024vC06R9YmDBlueV/dhP360bWaZcJ+MjJluAlXxBHPhJMiAxZCiJyQAQsAJvUkLWSXTiW1ek0Y0i/qXr8Jj4xr6BwZmTCkUEktrglD8qLucU0YmjcuXH7+JJXUEiADFkKInJABCwAm9DQyBlyN5piw94nO3HQTLo+jOSZc0dYsE4b0N/oczYTL+rTEDIl6KqklQAYshBA5IQMWAEzqXZHyGbM2YWi0klpSEy5va5YJe1vCSmqxTRgy2/K+lgl7Y0WfljbhhMiAhRAiJ2TAAoCJhZUZfTdkZcLlvZtrwtXbsjVh79NYJbXGTRgy2/K+hgl7UwuvmqtSSS0JMmAhhMgJGbAAYErPq6UXbWDCpSPLezfHhL1PurtrjGbClX2aY8IQv35EUhOGNqkfYem4qwxYCCFyQgYsAJhcbsBFWtiEIf195uo14cprN8uES5E2y4Sh8UpqDZswxK8fkXcltQTIgIUQIidkwAKANQqv1H5TJlzHtbM24fJImmTCkNmOyzVNuPxc7WDCCZEBCyFETigBCyFETmgIQgAwsVDHUmQNRdRx7WyGIsqP7OihCGiv7Y0SIgMWQoickAELANYoVJmGVosWMOFabZC9CVeeu9a10zXhakdmbsKQ3Zb3NUwYEpSyzMOEEyIDFkKInJABCwCmxDHgIl1qwt6W9kafMuGhf187bnnfIDJgIYTICRmwAGBKIeGWRKl/J63ehCFJKctkJuznSmt7o/pMuFp8WZuwH5O0qHtME4b0tjdqAxOWAQshRE7IgAUAUyxAUguGJpowpLe9UTwTrujbJBOG7La8r2XC3iet7Y3qNGHIbsv7FjRhGbAQQuSEDFgAMKXQC4ORWrWBCUN2W97LhPMzYT8mWVH3djJhGbAQQuSEDFgAMLkwHojmAsuEo/ejs7aECZe3RefI2oQhwy3vq5uwXz+6VrNM2E9GHsiAhRAiJyw0MfOb2XLg0aZdUMRhkxDC9LyDEKKbaGoCFkIIUUJDEEIIkRNKwEIIkRNKwEIIkRNKwEIIkRNKwEIIkRNKwEIIkRNKwEIIkRNKwEIIkRNKwEIIkRP/H6aeTnC7cT/+AAAAAElFTkSuQmCC\n",
+ "text/plain": [
+ "<Figure size 360x360 with 3 Axes>"
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "#%% plot the distributions\n",
+ "\n",
+ "pl.figure(1, figsize=(6.4, 3))\n",
+ "pl.plot(x, a, 'b', label='Source distribution')\n",
+ "pl.plot(x, b, 'r', label='Target distribution')\n",
+ "pl.legend()\n",
+ "\n",
+ "# plot distributions and loss matrix\n",
+ "\n",
+ "pl.figure(2, figsize=(5, 5))\n",
+ "ot.plot.plot1D_mat(a, b, M, 'Cost matrix M')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Solve Screened Sinkhorn\n",
+ "--------------\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 21,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Epsilon = 0.014767606916367452\n",
+ "\n",
+ "Kappa = 3.3854408965782907\n",
+ "\n",
+ "|I_active| = 30 \t |J_active| = 30 \n"
+ ]
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZhcZZXH8e/p6iXpJEASAgQCJGwBiYalQUDFyKKAC6ggwsimEERRWRRRZ0ZAcURRUEH2TUVgQJTIMOwyqEAwYd93QiBAICSELL2+88d5b1V1dXe609tbXf37PE89N/fWrVunK+mTU+99FwshICIig68qdQAiIsOVErCISCJKwCIiiSgBi4gkogQsIpKIErCISCJKwCL9zMw+YmZPp45joJjZKWb2h9V9TjpSApayY2aHm9mjZrbczF43s/PMbK343Plm9l58NJlZc9H+/w5CbMHMNlvVOSGEv4cQpvby+h82s3vMbImZLTKzf5rZDr2LVsqdErCUFTM7ETgD+A6wJrATsDFwm5nVhhC+GkIYHUIYDfwEuCbbDyHsnS5yZ2bVfXjtGsCNwG+AccAGwKlA42DFUE7MLJc6hoGmBCxlIyagU4FvhBBuDiE0hxBeAr6AJ+Ev9eKaM8xsvpmdZGZvmtkCM9vPzPYxs2dilfn9ovN3NLN7zWxxPPccM6uNz90dT3s4VtwHFl3/u2b2OnBZdiy+ZtP4HtvF/fXN7C0zm9FJuFsAhBCuCiG0hhBWhBBuDSE8UhTfUWb2pJktNbMniq77UozhEWCZmVXH9/qTmS00sxfN7JtF16kys5PN7Hkze9vM/tvMxsXnJsdK/zAzmxfj/UEXn2+NmV0V36c2Hq41s9/FGB83s4ai87cys7vi5/u4mX2m6LnL47edm8xsGfCxeOxcM/ufeL3ZZrZpD//6y18IQQ89yuIB7AW0ANWdPHcFcFXJsVOAP3RzzRnxmv8J1ABHAQuBPwJjgK2BlcAm8fzt8aq7GpgMPAkcV3S9AGzWyfXPAOqAkfHY/KJzjorXqQduAc7sItY1gLfjz7o3MLbk+QOAV4EdAAM2AzaOz70EPARsGGOoAubGn7sW2AR4AfhEPP844D5gUoz7guzzjT93AC6K15qOV+FbFX/u8bn/AS4HckXPrQT2AXLAfwH3xedqgOeA78eYdgOWAlPj85cDS4APxfhHxGOLgB3j38mVwNWp/63210MVsJSTtYG3QggtnTy3ID7fG83A6SGEZuDqeJ1fhRCWhhAeBx4HPgAQQpgbQrgvhNASvPq+APhoN9dvA34YQmgMIawofTKEcBHwLDAbmAh0Wk2GEN4FPkwh+S00s1lmtm485UjgZyGEfwX3XAjh5aJL/DqE8EqMYQdgQgjhtBBCUwjhhXjNL8ZzjwZ+EEKYH0JoxBPn/iXNF6cGr8IfBh7GE3FmDeBm4HngiBBCa9Fz/wgh3BSP/b7odTsBo4GfxpjuxJtcDip67Q0hhH+GENpCCCvjsetDCPfHfxdXAtt09vkNRRXRViQV4y1gbTOr7iQJT4zP98bbRQkiS5BvFD2/Ak8MmNkWwC+BBrxircYryVVZWJQsunIRMAuYGRNep0IITwKHx1i2xCvNs/EktSGe8LryStGfNwbWN7PFRcdywN+Lnv+zmbUVPd8KrFu0/3rRn5cTP6NoJ7yiPSjE8nUVrxsRE/v6wCshhOL3fBlv6+7sZ+hJHEOaKmApJ/fiX3U/V3zQzEbhX8nvGIQYzgOeAjYPIayBf122bl6zyikFzWw0nkQvAU7J2lq7E0J4Cv8KPi0eegVYVftncRyvAC+GENYqeowJIexT9PzeJc+PCCG82pPYgFvx5oU7iir07rwGbGhmxXlnI7xZpbOfoeIpAUvZCCEswW/C/cbM9oo3eCYD1wLz8a+zA20M8C7wXqxAjyl5/g28PXV1/AqYG0I4Em8zPb+zk8xsSzM70cwmxf0N8cr3vnjKxcC3zWx7c5uZ2cZdvOf9wLvxxtxIM8uZ2bSiLm3nA6dnrzezCWa27+r8UCGEn+Ft6XeYWU+ah2YDy4CT4t/tDODTeLPQsKQELGUl/lJ/HzgTT4Sz8Wpt91V9de9H3wYOxm8OXQRcU/L8KcAV8S7+F7q7WExqewFfjYdOALYzs3/r5PSlwAeB2bEXwH3AY8CJACGEa4HT8aS3FPgL3l2tg9jk8mm8vfRFvPnmYrxrH/h/CrOAW81saXyvD3b383TyPj+KcdzeXWUfQmgCPoN/m3kL+C1waKz0hyXr2HwjIiKDQRWwiEgiSsAiIokoAYuIJKIELCKSiAZiCABrr712mDx5cuowRIaUuXPnvhVCmNDb1ysBCwCTJ09mzpw5qcMQGVLM7OXuz+qamiBERBJRAhYZbpYtA/X/LwtKwCLDwbPPwqGHwqabwujRMH48fPzjcP31SsYJKQGLVLKWFjjlFJg2Df7yF9h2WzjtNNh/f3jxRfj85+GTn4R581JHOizpJpxIpWppgUMOgauvhoMOgl/8AiZObP/8OefAf/wHfOQjcNddMGVKsnCHI1XAIpWotdWbHK6+Gn76U/jjH9snX4DqajjuOLj7bli6FGbM8KpYBo0SsEglOvVUuOoqT77f/e6qz912W7jjDk/C++4LKzos6iEDRAlYpNLceSf8+MdwxBHdJ9/Mttt6lfzoo3DCCQMbn+QpAYtUkoUL4Utfgi22gN/8ZvVeu9de8J3vwPnnw3XXDUx80o4SsEglOekkeOstuOYaGDVq9V//4x9DQwMceywsWdL/8Uk7SsAileLee+Hyy+H442H69G5P71RtLZx3Hrz5prcjy4BSAhapBK2tXrWuv753K+uLhgaYORN+/Wt47LH+iU86pQQsUgkuvxweeMD7+o7uh1XbTz8d1lzTq2kZMErAIkNdY6OPbttxRzjwwP655vjx8IMfwO23w//9X/9cUzpQAhYZ6i65xIcS/+hHYNZ/1z3mGB+88R//ofkiBogSsMhQtmKFNxd8+MOw5579e+2RI+H734e//90rYel3SsAiQ9nFF8Nrr/V/9Zs56ijYcEP44Q9VBQ8AJWCRoaq5Gc4806vfGTMG5j3q6nw03b33wj//OTDvMYwpAYsMVdde622/J500sO9zxBF+U+5nPxvY9xmGlIBFhqIQPCFutZXP5zuQ6uu9j/Ff/wpPPjmw7zXMKAGLDEW33QYPP+xzN1QNwq/x17/uN+V+/vOBf69hRAlYZCg66yxYbz04+ODBeb8JE7wp4sorfZiy9AslYJGh5pln4OabvZ9uXd3gve83vgFNTXDhhYP3nhVOCVhkqDnnHKip8fkaBtOWW/pCnued5z0wpM+UgEWGkqVLfd6HAw/0JojB9s1ver/j668f/PeuQErAIkPJFVd4Ev7GN9K8/957+9L2v/51mvevMErAIkNFCPDb38IOO/jEOylUVcHXvgb33OO9MKRPlIBFhoq77/Z+uMcckzaOww+HESN86SLpEyVgkaHi/PNhrbX6b8rJ3ho3zmP4wx+8OUR6TQlYZCh44w3405+8+qyvTx2NV+Hvvef9gqXXlIBFhoJLL/WuX0cfnToSt+OOsM023iVNs6T1mhKwSLlra4OLLvIZz7bcMnU0zgy++lV45BGYPTt1NEOWErBIubv9dnjxxfKpfjMHHwyjRmlkXB8oAYuUuwsvhLXXhs9+NnUk7Y0Z40n46qth8eLU0QxJSsAi5ez11+GGG/zm22DO+9BTRx/tyyLpZlyvKAGLlLPLLoOWFl8aqBxtvz1stx1ccIFuxvWCErBIuSq++bbFFqmj6drMmfDoo7oZ1wtKwCLlqlxvvpXSzbheUwIWKVcXXFCeN99KFd+MW7IkdTRDihKwSDlasABmzSrfm2+lZs7UzbheUAIWKUflfvOtVEODbsb1ghKwSLlpbfWbbx/7WHnffCs1c6aPjLvvvtSRDBlKwCLl5tZb4aWXfKjvUHLwwTB6tFfB0iNKwCLl5vzzYd11Yb/9UkeyesaMgUMOgWuugUWLUkczJCgBi5STefPgxhvhK1+B2trU0ay+o4+GlSt96STplhKwSDm5+GK/iTVUbr6Vmj4ddt7Zq3jdjOuWErBIuWhq8ptve+8Nkyenjqb3jjkGnnkG7rgjdSRlTwlYpFxcf71PvnPssakj6ZsDDvABJOeckzqSsqcELFIuzjkHNtsMPvGJ1JH0zYgR3iXtr3/13hzSJSVgkXLw4IPwz3/C17/uS78PdVkXOq2cvEoV8DctUgHOOccX2zz88NSR9I8NN/RudBdf7EOUpVNKwCKpLVzocygccogvO18pvvENePttzQ+xCkrAIqmddx40NsJxx6WOpH999KO+cvJZZ6lLWheUgEVSWrkSzj0XPvnJ8lnxuL+YwQknwBNPwC23pI6mLCkBi6R05ZXw5pueqCrRgQfC+uvDL3+ZOpKypAQskkpbmyembbbxmc8qUW2ttwXfdhs8/HDqaMqOErBIKjfe6F/PTzzRv65XqqOP9lnSzjgjdSRlRwlYJIUQ4PTTYcoU+OIXU0czsMaO9eHJ11wDzz2XOpqyogQsksKdd8L998N3vwvV1amjGXjHHw81NfCzn6WOpKwoAYuk8JOfwMSJcNhhqSMZHBMnwpe/DJdfDvPnp46mbCgBiwy2u+/2Cvjb3/Z5E4aLk07yppf/+q/UkZQNJWCRwRQC/Pu/e0V4zDGpoxlckyf7RPMXXaRJeiIlYJHBdOut8Pe/exIeOTJ1NIPv3//dJxs69dTUkZQFJWCRwZJVvxtvDEcemTqaNCZNgq99DX73O3jqqdTRJKcELDJYrr4a5syBU04Zmuu99ZfvfQ9GjfIeIMOcErDIYFi+3BPOttvCoYemjiatCRPgBz+AWbPg9ttTR5OUErDIYDjzTHjlFTj77MqYcL2vvvUtH4Ry/PHQ0pI6mmT0L0FkoM2b58Nw998fdt01dTTlYcQI+PnP4bHH4IILUkeTjBKwyEAKodDd7Oc/TxtLufnc52CPPbxN+NVXU0eThBKwyEC65hq46Saf92EoLzU/EMx8zbiWFu8ZMQwnbVcCFhkob70F3/wm7LCDT8koHW26qfcJnjULrr02dTSDTglYZCCE4HMfLFniC1PmcqkjKl/HHw8NDb6S8iuvpI5mUCkBiwyEc8+Fv/7VZ//6wAdSR1Peqqvhj3+Epib40pegtTV1RINGCVikvz3wgE+0s88+3gQh3dt8c/9P6+674bTTUkczaJSARfrTggWw776wzjpw2WWVvdJFfzv0UJ+e87TT4LrrUkczKIbBTNAig2TlSvjsZ2HRIvjnPz0JS89lvSKeecaT8SabwHbbpY5qQKkCFukPTU0+0GL2bPj9732hTVl9I0bAn//sw5X33rviJ+xRAhbpq5YWOPhg+J//gfPO8wEG0nvrruvTdprB7rvD88+njmjAKAGL9MWyZZ5w//QnOOss70olfTd1qk/U09gIH/kIPPRQ6ogGhBKwSG+98QbstptXvueeC8cdlzqiyjJtGtx1l/eh3nVXuOWW1BH1OyVgkd64/XaYPh0efdTbLL/2tdQRVaZp0+Dee33mtL33hv/8z4qaPU0JWGR1LF0KJ5wAH/84jBvnN90+85nUUVW2SZO8V8mhh8KPfgQzZsDjj6eOql8oAYv0RGsr/OEPsNVW3tZ79NHwr3/B+9+fOrLhYfRoX9L+97+HJ5/0Xibf/rbPtzGEKQGLrMry5T6gYuut4ZBDvG/vvfd6b4dRo1JHN/x86Uvw9NP+d/HLX3rTxMknw8svp46sV5SARUq1tMDf/gZf/zpssIFPqlNb66Oz5syBnXZKHeHwtvbacOmlPpn7Jz/p8yxvson/+Xe/g8WLU0fYYxaG4Ryc0lFDQ0OYM2dO6jDSWLrUuzn9618+F8Fdd/ksZvX1Pqz4q1/1rlAaVlye5s2Diy6CK67w2dRyOZ8C9GMfgw9+0GdaW3/9Afn7M7O5IYSGXr9eCVigwhJwayusWOHNB0uX+uOdd3yI8JtvevexefP8a+szz8D8+YXXTpninf/32ssfamYYOkLw/0RnzYI774T77y/MrLbmmt63eMoU2GgjmDjRB3yMHw9jx/rzY8b433d9vc/Q1oOErQQs/aKhvj7M2WyzwXvD0n932X7x8RA6Ptra/NHaWni0tPijudk77nfXTcnMfwE32shn4dpyS58ysqEB1luvf39OSWf5cnj4YZg714c0P/20/6c7b57/O1kVM6ir86anmhpPyLlc4bHJJnDHHX1OwJqMR1xdHQxmAoaOFUa2X3zcrP0jlytsix81Nf4YMcJ/lpEj/TFmjD/GjvVuY+us4/MMVOuffsWrr4edd/ZHsRC8nfiNN+Dtt/2b0bvv+jel5cv929PKlZ6km5r8P/aWlsJ/+G1t/TbRkv4Vitt0U7j++tRRiAw8M/8PeezY1JGoF4SISCpKwCIiiegmnABgZkuBp1PH0Y21gXIf+jQUYoShEedQiHFqCGFMb1+sNmDJPN2Xu7mDwczmKMb+MRTiHCox9uX1aoIQEUlECVhEJBElYMlcmDqAHlCM/WcoxFnxMeomnIhIIqqARUQSUQIWEUlECXiYM7O9zOxpM3vOzE5OHU/GzDY0s7+Z2ZNm9riZfSseH2dmt5nZs3GbfDypmeXM7EEzuzHuTzGz2THGa8ysNnF8a5nZdWb2VPw8dy63z9HMjo9/z4+Z2VVmNqIcPkczu9TM3jSzx4qOdfrZmft1/F16xMy26+76SsDDmJnlgHOBvYH3AQeZ2fvSRpXXApwYQtgK2An4eoztZOCOEMLmwB1xP7VvAU8W7Z8BnBVjfAf4SpKoCn4F3BxC2BKYjsdaNp+jmW0AfBNoCCFMA3LAFymPz/FyYK+SY119dnsDm8fHTOC8bq8eQtBjmD6AnYFbiva/B3wvdVxdxHoDsCc+Wm9iPDYRH0CSMq5J8ZdwN+BGwPDRW9WdfcYJ4lsDeJF4w73oeNl8jsAGwCvAOHxw2I3AJ8rlcwQmA49199kBFwAHdXZeVw9VwMNb9g8/Mz8eKytmNhnYFpgNrBtCWAAQt/0zL2DvnQ2cBLTF/fHA4hBCNilx6s90E2AhcFlsJrnYzEZRRp9jCOFV4ExgHrAAWALMpbw+x2JdfXar/fukBDy8dTblf1n1SzSz0cCfgONCCO+mjqeYmX0KeDOEMLf4cCenpvxMq4HtgPNCCNsCyyiPZpu82Ia6LzAFWB8YhX+dL1VW/zY7sdp/90rAw9t8YMOi/UnAa4li6cDMavDke2UIIZus+A0zmxifnwi8mSo+4EPAZ8zsJeBqvBnibGAtM8vmWUn9mc4H5ocQZsf96/CEXE6f4x7AiyGEhSGEZuB6YBfK63Ms1tVnt9q/T0rAw9u/gM3j3eZa/MbHrMQxAX5HGbgEeDKE8Muip2YBh8U/H4a3DScRQvheCGFSCGEy/tndGUL4N+BvwP7xtNQxvg68YmZT46HdgScoo88Rb3rYyczq4997FmPZfI4luvrsZgGHxt4QOwFLsqaKLqVqeNejPB7APsAzwPPAD1LHUxTXh/Gvb48AD8XHPngb6x3As3E7LnWsMd4ZwI3xz5sA9wPPAdcCdYlj2waYEz/LvwBjy+1zBE4FngIeA34P1JXD5whchbdLN+MV7le6+uzwJohz4+/So3ivjlVeX0ORRUQS6VMTRLl24hcRGQp6XQHHTvzP4H0z5+PtiQeFEJ7ov/BERCpXX1bE2BF4LoTwAoCZXY13JekyAa+99tph8uTJfXhL6auFC+Gdd2CLLdofH/vAlDQBiQxht7Vd21nXsx7rSwLurNPxB0tPMrOZ+LA8NtpoI+bM6dMKHtJH3/kOnHsulP417Fl1QJqARIaxvrQB96jTcQjhwhBCQwihYcKECX14O+kPLS1QU5M6ChGBviXgsu7EL51rbITapHNziUimLwm4bDvxS9eWL4f6+tRRiAj0oQ04hNBiZscCt+DTx10aQni83yKTAfHeezB6dOooRAT6dhOOEMJNwE39FIsMgrffhvHjU0chItDHBCxDz+uvw9Zbp46izFmfehb1jUamDiuajGcYaWuDF1+EKeryK1IWVAEPIy+84L0gSgdhDFs9rXRtMOuUtu5PyahaHvJUAQ8j//iHb3fZJW0cIuJUAQ8js2bBOuvAVluljiSxrirfLipdq1pFpdxddVz62rZuqtbQ9XuF0tda/okuXqAKudypAh4mFizwBHzYYVClv3WRsqAKeJg49VQviI46KnUkCfWw8u1Q8WbPd1IJW3bN0v/VetuTIqta2zqpaksr2qwiLqma85WyKuSyp1poGPi//4MLLoDjjoPNN08djYhklIAr3OOPwwEHwCabwGmnpY5GRIqpCaKCPfYY7LYbVFfDTTfBqFGpI0pkdZseSpocOm1mKDlWOKfza3QbS5fNC4XmA4vHQmkzRZf77ZsoOtzE66rLm5omBo0q4ArU1ga//S3stJMn37vugqlTu32ZiAwyVcAV5sUX4cgj4c474ROfgIsvhkmTUkdV3rqtfHO5eLxQr1iuqv1rsv14rmXnZtvsWqWVclfy1WtRNdra6peIx0JrrGDbWtu9xuJ5IW6za1iskPOv6+rmXWeVsariAaEKuEI8+6z3cJg6Ff71L7joIvjf/1XyFSlnqoCHuIcegp/+FK691le6OOooOPlk2HDD7l87bHUzeKJDm2/Wzpsrel1WFVf7r5DFbf411f58KKmM85VyfI+Qf68uYmktqjyzKjRf+frWWmKlm22ztuCWlnh+a8l+W/vj2fmdVMbddmlTZdwnSsBD0Pz5cPXV8Mc/woMPwpgxvtbbccfBeuuljk5EekoJeIhYtAiuu86T7t13e+Gxww5w1llw+OGw1lqpI6wAXfRY6NAWDB0r3xrfhlj55vdrcp1vYyUc4nuGXNZG3HV4WS8IWrO2Xt9WtcRKuDmrdLP9lnZbsv1YCYfm5ni9Lirm4j932aNCFXFfKAGXqZYWX7n4llvg1lth9mz/XZg6FU45BQ4+GDbbLHWUItIXSsBl5OWXPdneeivcfjssXuxFWEODt+t+7nOw7bZp5wsfkrKqLPvgsnbM0rbgrKrLlXzAVSU9GOjYTpyvfGt9yelQ679abSPiNla+bbV+fmvcttX6ddqqY2UZQwpVnVTE2Y+R9XaIP0ZVS6yEm+O2yZ/INbbF417FVq2MFXBTrHSbvALOV8qNTf42WSVMcbVcUh3n26NjjwtVxL2iBJxISws8+ijcc48/7r3Xu5CB91z43Oe8G9nuu2sJIZFKpQQ8SN55B+67r5BwZ8+GZcv8uYkT4UMfgm99Cz7+cdhyS1W5Za+0P29WCWfbrOKN29aRvm0ZGSvgOt+21PnrW2v9Mm3ZtqQi9veI23wlHLexKM1XwLFpt7rR93NxW70iVsYrs61XtVUr4rYxVsQrmwrvGatksvbibL+0vbhk22XvCVXE7SgBD4DFi713wty58MADvn3mGX8ul4Pp0+GII3xi9F12gY02UsIVGY6UgPto0aJCks22zz9feH7DDWH77eHQQz3Z7rCDloVPpou24NAW+/lWxf14p986zKnQSZ/c0v2sDTe2I4dYybbVZJWub5tH+ral3s9vifutdfEt47a1tvA+ITYzh5Km69JKONfk16pq9G2u0Y9Xr/AXVi/3/ZoV3l5ds6w1HvdtbkWhDbhquVfDtjJWx7GdON9e3Ny+HTmUVsbZ/BX5jhVqIy6mBNxD2YKWDz/c/vHSS4VzJk/2ZPvlL/t2u+1gwoRUEYtIuVMC7sSyZX6DrDjRPvooLF3qz1dV+cKWH/wgHHOMJ9rttoNx49LGLf2krX0jazanQvFE7aEqzvnQ2r5t00r66ObPz0a+5bJeD348q4hbR/h+82h/XetI37aNLKoY62KFWht7N1RlpW+sMmMl39Qa36Mxxhi3Vcv9+Zpl/nz18rh9z4/XvlcdtzX5t6xe5o3SNe/FduLlzXHbGK8d24tXxjK7q4o4++y66zUBw6oqHtYJOAQfVVZa1T77bOHfwBprwAc+4E0I06fDNtvA1ltDfX3a2EVk6Bs2CbixEZ54omOyXbSocM6UKZ5kDz7Yt9One7OCbpBVmF62BednEQMs135+BWuOvR+y0XL50WhxW1IZWzbfQhZCfFlbLD6zyjc3pjn/niPrvcpcY+RKAEbXePU5sjpWpbFib4uNxCta/GLvNXsVu3SlNywvX+bl9sql/nxuqb95zVIPpvbdwoi/mqV+rbql1XG/pt02tyyriP09LKuEs8q4tPdEd70m/CegnQquiCsyAS9f7pPUzJ3ro8keeACeeqrQc2bkSHj/++Hzny8k2g98wKtdEZHBMuQT8IoVXslmyXbOHK90swme1l3Xb4h95jOFZLvZZu2H9csw1V0lTKzO4untvgg1W/vnsnkjYj/gqmzWs+qsX3D7kW+5pqzPbnxdLAitpPiryhXmZRhV51XlhJHegXz9+iW+X+s3J9bMrQCgJl6sNUa8PHateCd2uViwck3fLveKY+F7vlTKe0tGAtC4uDb/nrVL/BpNS6ra7dfFKrn23dKKOFbCWRvxiqwLRqyIs0q4Ke5nn33Rt4su55uowEp4yCXg117zRSbvussHMzz2WOHbzIQJPmz3s5/17fbbw/rrqwlBRMpT2SfgBQs82WaPbEDDGmv4kjuf+pQn2oYGH8KrZCurbTUr4WL5gV4lFXBpRZyLvR+qs/7BcdtaE3s7xH7CbfleEXHEXHPhq1pza/sOwCOrvIpcu/o9ADasfRuAcTnfH2VN7c5fGfzXfWGrV76vNPkY95dW+vb5sd5n8pV318y/ZvFir46bFnkVXbvYY2h+J7YXj/LtiPrYwyJWxjWxR0VV7FFRtdzbrfOVcPYNoXRkHUW9TrLPvS37uSuvbbjsEnAIvqLD737nE9I8/bQfX2MN2HVXmDkTZszw3ghqRhCRoaxsEvBbb8Hvfw+XXurNCiNH+oq+Rx5ZSLjVZROtVKQeVsLQsRq2kuPZfla7WRfVWrD4j9razw/cFivk5upCe+y7Ob/GwhqvFteo9Tbf8TXeJtwaez9kle+6Od+OiX2Wc7GCXBneAGBRrW9fGzEGgGfqfTb/p0cXZvV/emrvfKwAABKISURBVMy6AMwb4xNOLx3tFXFLfU3cxhF9ca6LESNi/HGui5o4F0YubquWxb7J+RVDYltxU6G3R5ftwxXYNtxtSjOzDYHfAevhP/GFIYRfmdk44BpgMvAS8IUQwjurG8C8efDtb8Nf/uI9VnbcES64AA48ENZcs/vXi4gMVT2pKVuAE0MID5jZGGCumd0GHA7cEUL4qZmdDJwMfHd13vytt2DPPb2d99hjfQjvtGmr+yOI9LNuKmEobp/M5osoqcZKViLOKuBc6Yi5bDBbHBpnsZ3X2uJ8Dq2FdrbG2JvhjeCVSVusCNtKJofIxYvW2FsAjIi9IuqrvP/v6Lhds8qPj6vyXhTjcj5JxDrV7+avlVXXa9VOBOD5Om8vXjjC25FX1HmF3lZbMt9xbNOui/u1sRdITewVkq+ES1ePhvyouZDNV4xvs/kk8v20K6BtuNsEHEJYACyIf15qZk8CGwD7AjPiaVcAd7EaCbix0W+gzZvnbb0f+tBqRi4iMsStVquqmU0GtgVmA+vG5EwIYYGZrdPFa2YCMwE22mij/PHXX/c+uw0NPkOYSNnpUEkVKq0Oo+ay6iyr1rKKuKRCziriqriCcU1LtmqFV5K5Jv+VzMXeD1VNhcqwqsmPNTZ6Bfv6St9f1uivXdzk/XjfHuPttG/We5W6Re3rAKxfHSvdeMkRFnssxPbncVXettxcXRgemvXGyMcQq+uqOP/Em+bv0Zjz6rytw1p32X5cIy9eria2cedK192jqCrOjmXbbHWOrnpJ5Fc6yb65lH8lvOr1uYuY2WjgT8BxIYR3uzs/E0K4MITQEEJomFA0NdjGG8OFF3pf3sMP9+YIEZHhpEcJ2Mxq8OR7ZQjh+nj4DTObGJ+fCLy5um/+5S/DT34CV13lAya+8AW4+eb2i7KKiFSqnvSCMOAS4MkQwi+LnpoFHAb8NG5v6E0A3/sefPrTcMkl3g3t2mt9QMXhh8N++6m/r5SRdl9p2zc9ZFNCdlh6J3/TLX5dzr4+xyaIXNxmS8pnTRFVjd7Nq3pl4Vd0ZX5Cdf+K3RgnwHkvDnp4cqk3TSxYy7uVPb+Gf+PMupVNHuEDNbIBGxNy/kV2hGWDIGKTRih0fctu6I3J+UCKCbU+yOPdkd7csXKMv/fiOAVmc4s3RVTF/exmomVNNqHzlFNV1ASR77oXt50OBW93fOg2RfSkAv4QcAiwm5k9FB/74Il3TzN7Ftgz7vfKtGlw1lnw6quegKdNg9NP9/bh8eM9Qf/iFz6pjqpjEakUPekF8Q86/ueT2b0/g6mrg/3390c2BPlvf/PtjTf6OWuu6SPiZszwnhPTp8OIEf0ZhUgPdHGDrsPw5dJuavlpGON+vLFUlS0N3xQryDh1ZPWKQjVavdyr4qZl2XDfOFHOu3Ei9iX+i/DOmnHynbV87asX1vSVAiaM9i5lE+u98p04wifzGRvXKKqPgyJyRcNMmuM8mUtaveJd0VaIB6AmF6fNrPGfq7nO97Mllqqy5ZHipENVLdnkQ/EmYycT2IeSrnvZjcuQ7dPeUK6Ey3Zs2cSJcNBB/gCvjrNJeO66C/76Vz9eXe0VczYfREODTzVZV5cqchGRninbBFxqgw18ovSDD/b9V1+F++8vTEH55z97OzJATY0n4Swhb7edr2KhSlkGTL66Wr224fwkNPmJ3UsWv1xRqCSqlseq+L32U0A2xwlwGhfHdtgxsSJew//BLx/jr3thjLcNvzzGB1PUj/J23TEjvPItneAdup/kPesC19ba/kZN1t2sdKL5/BJMcdsyIg7IaC4aTBHfo6q1fde9/MT4+TdpXxEPxUp4yCTgUhts4NNOfvazvh8CvPxyISHPnQv//d/e1Q38Rt7UqYU5gbPlhdZbr+v3EBEZSEM2AZcy8+WDJk/2NmTwpPzCC37zLluC6B//8G5vmXXWaZ+Up0+HLbf0KlpktXXVNpxVxFmbcNYWnHXxyfZjWzAllTBALv45W/4ntyxOERknRa8d7b/OLaP8mk2jYyU8Olab2fP1vl0+yivk9+pjjHHRz6q6wp3uqjgBUDbgJJsQpy1Wl63ZQJFs4c/YU6MqFtHZRPNZeZqtxBRKpuXMJqwH8pPYk03YU5MN087a0WNFnMsq3bifLXxqJYNkyrgSrpgE3Bkz2HRTfxxwQOH4okXwyCPt14b7zW98eDRAbS28730dE/P48Wl+DhGpTBWdgLsybpz3opgxo3Csudkney9OyjffDFdcUThngw0KTRdZUt588/bziIi009OKuC1b+DPrJ5xVwoWJyvNV8cq4/E9cALMqLgOUi23CbfW+Xzcyq3i9kmyOlXFz7KEQVyjKTyXZWhcn1Cm6gd1aG+OKhXrJvD9UZ0VlVtBnvR5iMZPLttkSTC3Z+dnipB2r0NLqOD9MOQ5rtmx4c1tsP84+u9L29SHQZXVYJuDO1NT4jbqtty7c6AN4442OKynfemvh3sno0Z6Qt9/eb/Ztv723NWvuYhHpjtJEN9ZdFz7+cX9ksiXuH3oIHnzQb/hddJGvxgw+mfw22xQS8vbbw1ZbqV1Z6L4izm/bj5yDQlWcn6Yxax+uzZb9idvYSyI3wivh6rrYRhwr4tYR2dYryZaRsT23LlbGdYWetlnvhWzJpA6VcHZqflpN31a1xIo3W4Eobqsb47Sc+UVJYyXcUvhc8n2Cs44RpeuMlU7SU9V+f7XbgiFZe7AScC/U1cG22/rjiCP8WGurL580d67f9Js715svzj3Xnx850ieb32UXf+y8s9qURYY7JeB+ksv5jbv3vQ8OOcSPtbXBs896Mr7/frjnHvj5zwvNF1OnFhLyLrt47wu1Jw8z+XbLbD5La7ebbyOGjj0nmkuWeI+9BbL2L1sep7aMk6ZXxeerY8Ucan2/LVbEbTVZW3DhH2E2wXp+0dDqGF9WTOarz+znibtt7SecL1TEWQUcp7XMts1FlX78s+Wn9CyZ2rMMei/0FyXgAVRV5Ul26tRCu/Ly5d5P+Z57/DFrFlx2mT83fjzssUehyWPSpHSxi8jAUwIeZPX1PpfFrrv6fgheJd97r897ceutcM01/tz73ldIxrvuCqNGpYtbBsmqJoEvrYqz6jP2iw3ZaLpctuBl7CXQ2H4peIs3IyxWylllnPW3ra4pSguxT26I1XFbtp/LeiZkFXEX08VkFXF+CaZs7ofYdzerdjupgLMeINnMcZTMCZH/rNq6qIjzn0/JcStpC4ZkfYP1hTcxM9hiCzjsMLj8ch9i/cgjcOaZXgGffz7ss493ndt3X0/O2c0+ERnaVAGXGTOfx+L974cTT4QVK+Dvf4ebbvKpOmfN8q5v++3nzRp77KHeFRWt04qsi/kmSnpQ5Cu9XGwrrsqq1zjPRC4bcRbTQKyYq4r7UFa1H5VWFbf549k1rH3f3a7kezjk++qWzBJH0Yi30sq3ZN6MfE+R0oo401VlXEZUAZe5kSO9CeLss30B0zvvhC9+0afn3GcfHxxyyinw9tupIxWR1WVhENs8Ghoawpw5cwbt/SpZY6OP1Lv0Uq+KR42CmTO9at5gg9W/3p5VB3R/kpS3Dv1ls36wXRzPqteSRTCteAma7M/ZNapy7V+bH6WWa7/fVZtwprRHQ3EeyireruZQDtk3gJJeEqUVcf49suuVzEzXmdXMh7e1XdvND7pqqoCHqLo6bxO+4QZ49FGfFe7Xv4YpU+CEE9ROLDIUqA24Akyb5uvpnXaaL3J61lneRHH55d6/WIaJ7kbZZaxk3ol8X95sCfnC/BP5eRhKquQQ963DaLSq9vuFC/UwZjpUrpTO9ZBVvCX73Va+ZUgVcAWZMsWHRN95p/fR//CHvX24gvqti1QUVcAV6GMf865sxx4Lp57qzRFnnNF1ESIVqsv/edu3gXbZmwI6qWTbV7j5VShK/3GVVMy90qFXQ0mlmz/eedtuh8q3H9t++4sScIUaM8abIEaN8uHPtbXw4x+njkpEiikBVzAznwxo5UpvG957b19JWoa5rqq9kvkoig/ldVYlQ9e9HqyLVs5V9ZLors22i0q2y7berirfMmibUxtwhTPz3hEbbQRHHlmYt0VE0lMFPAyMHg2/+pWPnvvLX+ALX0gdkZS1VVWGnVTJxYc7sC6uNRCrVayqjbfdeekr34wq4GHi05+GjTeGCy5IHYmIZJSAh4mqKjj0UJ9xbcmS1NGICCgBDysf/ah/+7r33tSRyJAXQg8fbYP36GmMZUQJeBjZcUffPvBA2jhExOkm3DAyZowvMvrii6kjkWGjzCrOcqMKeJjZeGOf1lJE0lMCHmbGj4dFi1JHISKgBDzsjB0L77yTOgoRgdVIwGaWM7MHzezGuD/FzGab2bNmdo2Z1Q5cmNJf6us1V7BIuVidCvhbwJNF+2cAZ4UQNgfeAb7Sn4HJwBg50teZE5H0epSAzWwS8Eng4rhvwG7AdfGUK4D9BiJA6V81NYU1DUUkrZ5WwGcDJ1GYSHQ8sDiEkP0qzwd6sRKZDLbqap+sXUTS6zYBm9mngDdDCHOLD3dyaqcd/sxsppnNMbM5Cxcu7GWY0l9yucKKLiKSVk8q4A8BnzGzl4Cr8aaHs4G1zCwbyDEJeK2zF4cQLgwhNIQQGiZMmNAPIUtfVFUVFpkVkbS6TcAhhO+FECaFECYDXwTuDCH8G/A3YP942mHADQMWpfQbMw1OEikXfekH/F3gBDN7Dm8TvqR/QpKBpAQsUj5Way6IEMJdwF3xzy8AO/Z/SDKQtDCnSPnQSDgRkUSUgIcZVcAi5UMJWEQkESVgEZFElIBFRBJRAhYRSUQJWEQkESVgEZFElIBFRBJRAhYRSUQJWEQkESVgEZFElIBFRBJRAhYRSUQJWEQkESVgEZFElIBFRBJRAhYRSUQJWEQkESVgEZFElIBFRBJRAhYRSUQJWEQkESVgEZFElIBFRBJRAhYRSUQJWEQkESVgEZFElIBFRBJRAhYRSUQJWEQkESVgEZFElIBFRBLpUQI2s7XM7Doze8rMnjSznc1snJndZmbPxu3YgQ5WRKSS9LQC/hVwcwhhS2A68CRwMnBHCGFz4I64LyIiPdRtAjazNYBdgUsAQghNIYTFwL7AFfG0K4D9BipIEZFK1JMKeBNgIXCZmT1oZheb2Shg3RDCAoC4XaezF5vZTDObY2ZzFi5c2G+Bi4gMdT1JwNXAdsB5IYRtgWWsRnNDCOHCEEJDCKFhwoQJvQxTRKTy9CQBzwfmhxBmx/3r8IT8hplNBIjbNwcmRBGRytRtAg4hvA68YmZT46HdgSeAWcBh8dhhwA0DEqGISIWq7uF53wCuNLNa4AXgCDx5/7eZfQWYBxwwMCGKiFSmHiXgEMJDQEMnT+3ev+GIiAwfGgknIpKIErCISCJKwCIiiSgBi4gkogQsIpKIErCISCJKwCIiiSgBi4gkogQsIpKIErCISCJKwCIiiSgBi4gkogQsIpKIErCISCJKwCIiiSgBi4gkogQsIpKIErCISCJKwCIiiSgBi4gkogQsIpKIErCISCJKwCIiiSgBi4gkogQsIpKIErCISCJKwCIiiSgBi4gkogQsIpKIErCISCJKwCIiiSgBi4gk0qMEbGbHm9njZvaYmV1lZiPMbIqZzTazZ83sGjOrHehgRUQqSbcJ2Mw2AL4JNIQQpgE54IvAGcBZIYTNgXeArwxkoCIilaanTRDVwEgzqwbqgQXAbsB18fkrgP36PzwRkcrVbQIOIbwKnAnMwxPvEmAusDiE0BJPmw9s0NnrzWymmc0xszkLFy7sn6hFRCpAT5ogxgL7AlOA9YFRwN6dnBo6e30I4cIQQkMIoWHChAl9iVVEpKL0pAliD+DFEMLCEEIzcD2wC7BWbJIAmAS8NkAxiohUpJ4k4HnATmZWb2YG7A48AfwN2D+ecxhww8CEKCJSmXrSBjwbv9n2APBofM2FwHeBE8zsOWA8cMkAxikiUnGquz8FQgg/BH5YcvgFYMd+j0hEZJjQSDgRkUSUgEVEElECFhFJRAlYRCQRJWARkUSUgEVEElECFhFJRAlYRCQRJWARkUSUgEVEElECFhFJRAlYRCQRJWARkUSUgEVEElECFhFJRAlYRCQRJWARkUSUgEVEElECFhFJRAlYRCQRJWARkUSUgEVEElECFhFJRAlYRCQRJWARkUSUgEVEElECFhFJRAlYRCQRJWARkUSUgEVEElECFhFJRAlYRCQRJWARkUSUgEVEErEQwuC9mdlC4OVBe0NZHRuHECakDkJkOBnUBCwiIgVqghARSUQJWEQkESVgEZFElIBFRBJRAhYRSUQJWEQkESVgEZFElIBFRBJRAhYRSeT/AQ6KGRM5MOlWAAAAAElFTkSuQmCC\n",
+ "text/plain": [
+ "<Figure size 360x360 with 3 Axes>"
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "# Screenkhorn\n",
+ "\n",
+ "lambd = 1e-2 # entropy parameter\n",
+ "ns_budget = 30 # budget number of points to be keeped in the source distribution\n",
+ "nt_budget = 30 # budget number of points to be keeped in the target distribution\n",
+ "\n",
+ "Gsc = screenkhorn(a, b, M, lambd, ns_budget, nt_budget, uniform=False, restricted=True, verbose=True)\n",
+ "pl.figure(4, figsize=(5, 5))\n",
+ "ot.plot.plot1D_mat(a, b, Gs, 'OT matrix Screenkhorn')\n",
+ "\n",
+ "pl.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "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.7.4"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 1
+}