{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Estimation of time delay and linear dispersion"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The signal of interest is\n",
"\\begin{align}\n",
" z(t)=A e^{-\\frac{(t-t_c)^2}{2b_w^2}}\\sin(2\\pi f t)\n",
"\\end{align}\n",
"of width $b_w$ and frequency $f$. In this example we will take $A=b_w=f=1$ and set $t_c=0$.
\n",
"The probability density function (PDF),\n",
"\n",
"\\begin{align} s(t) = B(z)(t) := \\frac{z^2(t)}{\\int_{\\Omega_s} z^2(t)dt} \\end{align}\n",
"\n",
"In this example, we are interested in estimating the time delay ($\\tau$) and linear dispersion ($\\omega$) parameters, i.e. $g_p(t) = \\omega t - \\tau$. Therefore $z_g(t)=z(\\omega t - \\tau)$ and $s_g(t) = B(z_g)(t)$. \n",
"\n",
"Normalized measured signal with noise $\\eta \\sim \\mathcal{N}(0,\\sigma^2)$ is given by,\n",
"\\begin{align} r(t) = B(z_g + \\eta)(t) \\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Parameter estimation in the CDT domain\n",
"Let, $\\alpha = \\omega$, $\\beta = -\\tau$, and $\\hat{s}$ and $\\hat{r}$ be the CDTs of $s(t)$ and $r(t)$, respectively. The Wasserstein distance between $(g_p^{-1})'r\\circ g_p^{-1}$ and $s$ is given by,
\n",
"\n",
"\\begin{align}\n",
"W^2((g_p^{-1})'r\\circ g_p^{-1}, s) = \\|g_p \\circ \\hat{r} - \\hat{s}\\|^2_{\\ell_2} = \\|\\omega\\hat{r}-\\tau-\\hat{s}\\|^2_{\\ell_2} = \\|\\alpha\\hat{r}+\\beta-\\hat{s}\\|^2_{\\ell_2}\n",
"\\end{align}\n",
"\n",
"Let, $\\vec{p} = \\left[\\begin{array}& \\alpha\\\\ \\beta\\end{array}\\right]$. The estimates are calculated as,\n",
"\\begin{align}\n",
" \\widetilde{\\vec{p}}= \\left[\\begin{array}& \\widetilde{\\alpha}\\\\ \\widetilde{\\beta}\\end{array}\\right] =\\left({\\bf X}^T{\\bf X}\\right)^{-1}{\\bf X}^T\\hat{s}\n",
"\\end{align}\n",
"where ${\\bf X}\\equiv \\left[\\vec{\\hat{r}}, \\vec{1}\\right]$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example\n",
"#### Create $s(t)$ and $r(t)$"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"SNR: 9.47544940682685 dB\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAl8AAAFBCAYAAABJi8prAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3debgcVZk/8O+bm+Rmg4QlLhAwERCJgAKZRERBWTS4ISPRqCxugzi4sijKT0ZQHwUdkRlwGJQRRCAgKkQnGBdQZBEIhEACBMIeEiEhYUnIdu99f3+cPlOn656qOtVd3VXd9/t5nvv0Vt1dN7n33G+/561ToqogIiIiovYYVvYOEBEREQ0lDF9EREREbcTwRURERNRGDF9EREREbcTwRURERNRGDF9EREREbcTwRURERNRGDF9EREREbcTwRZUlItuIyDMiskvGdteIyEnt2i8iojQcuyiLcIV7qioR+T6A7VX1E859PwCwp6rOdO7bC8BfAUxR1Rfav6dERBHf2FW7v2784tg1dLHyRZUkImMAfBrAxbGH/gnAHe4dqnofgEcBHN2evSMi8hORreEfu4DY+MWxa+hi+KLSiPEVEVkqIhtE5FkR+VXt4XcDGABwS23bESKyGcCBAL4hIioiS5yXmwvgI239BohoSBORSbWxaLaI3CAiGwGcAGfsqm2XNn5x7BqCGL6oTKcC+ASAfwXwegDvB/DH2mNvA3CXRvPi/QD2r12fAeDVAN7qvNYdAKaLyOhW7zQRUc2bapdfBfADAG8AsBPqxy4gffzi2DUEDS97B2hImwlgnqr+uXb7CQB/r11/DYCVdkNVHRCRVwN4CcCdOrhZcQWAEQB2APBIS/eaiMh4I4CNAGap6jIAEJGd4IxdQOb4xbFrCGLli8o0F8CXRORPIvIZEdneeWw0zKDm2gfAIk/wAoANzvOIiNrhTTAfIJc59/nGLiB5/OLYNQQxfFFpVPVHAHYH8HuYqcdHRGSP2sOrAWwTe8qbACxMeLlta5erit5PIqIEb4Q5WtHlG7uA5PGLY9cQxPBFpVLVZar6AwDTAAiAvWsPLQQwNbb5GwHcm/BSewJYoarPtGRHiYgcIjIWwC4A7o495Bu7gOTxi2PXEMTwRaUQka+KyMdFZKqIvA7ANwFsBvCX2ibzAewhIts5TxsO4PUisoOITIi95NtgKmhERO1gPyjeE7vfN3YByeMXx64hiOGLytILc4TQAgC3wnwqPMR++qutf3MHgNnOc06v3V4O4Lv2ThEZBeBIAD9py54TEZkx62FVXefemTB2AZ7xi2PX0MUV7qmyRGQmgPMATFXV/pTtTgRwhKq+s207R0SUgGMXZWHliypLVX8P4AIAkzI23QLg863fIyKibBy7KAsrX0RERERtxMoXERERURsxfBERERG1UUedXmj77bfXyZMnl70bRNQmd91112pVnVj2fhSB4xfR0JM0hnVU+Jo8eTIWLFhQ9m4QUZuIyBNl70NROH4RDT1JYxinHYmIiIjaiOGLiIiIqI0YvoiIiIjaiOGLiIiIqI0YvoiIiIjaiOGLiIiIqI0YvoiIiIjaiOGLiIiIqI0YvoiIiIjaiOGLiIioqhYtAr7yFeCll8reEypQR51eiIiIaEh505vM5cAA8IMflLsvVBhWvoiIiKruoYfK3gMqEMMXERFR1YmUvQdUIIYvIiKiqmP46ioMX0RERFU3jH+uuwn/N4mIiKqOla+uwvBFlXHHHcAtt4Rtu3kzcMUVwJo1rd0nIiKionGpCaoEVWDGDHP9xReBrbZK3/70081R10ccAVx7bev3j4iorX73O2D9+ug2px27CsMXVcLGjdH1xYuB/fdP3/4//9NcXndd6/aJiKg073tf/W1OO3YVRmmqhLVro+v33JO9/aZN5nLcuNbsDxFRpTB8dRWGL6oEN3wtXpy+rWp0feedW7M/RESVwvDVVRi+qBKefz66/sIL6dvaqhcA9Pe3Zn+IiCqFPV9dhf+bVAlu5WvDhvRt3cfd0EZE1LVY+eoqDF9UCW6Ievnl9G3dx9eurZ+GJCLqSgxfXYXhiyrBrXzlCV+bN2dXyoiIOh7DV1dh+KJKaLTyFX8uEVFXYvjqKgxfVAlukz3DFxF1hcsvBy65pJjXYvjqKlxklSrBnTrMG77cBVqJiCrj6KPN5cc+BowYEf48XyMrw1dXCap8ichMEVkqIstE5DTP470iclXt8dtFZHLt/sNE5C4Rua92ebDznL/UXvOe2tcrivqmqPO4y0fkDV/uc4mIKmdgIN/2vvDFpSa6SmblS0R6AFwA4DAAywHcKSJzVfV+Z7NPAVirqruKyGwAZwP4MIDVAN6nqitEZE8A8wHs6DzvY6q6oKDvhTqYW71i5YuIOp4buPIeks3KV9cLidLTASxT1UdVdTOAOQCOiG1zBIBLa9evAXCIiIiqLlTVFbX7lwAYJSK9Rew4dZd45SttrGLli4gqzw1fN90EfPrTwIsvhj2X4avrhfR87QjgKef2cgAzkrZR1T4ReQHAdjCVL+uDABaqqvun8mci0g/gVwC+rTr4J05EjgdwPADszHPJdC23ejUwYJaQ6E2I6Qxf1Ck4fg1hbvh617vM5TbbAN//fvZzGb66Xkjly/c/Hv/JSN1GRN4AMxX5Gefxj6nqXgDeVvs6xvfmqnqRqk5T1WkTJ04M2F3qRPEAlTb1yGlH6hQcv4YwX5/XihWD7/Nh+Op6IeFrOYCdnNuTAMR/gv5vGxEZDmA8gDW125MA/AbAsar6iH2Cqj5du3wJwBUw05s0RMUDVJ7wxcoXEVWO78SzoQGK4avrhYSvOwHsJiJTRGQkgNkA5sa2mQvguNr1owDcoKoqIhMA/C+Ar6nqLXZjERkuItvXro8A8F4Ai5v7VqiTsfJFRF3FV/li+KKazPClqn0APgdzpOIDAK5W1SUicpaIvL+22cUAthORZQBOAmCXo/gcgF0BfCO2pEQvgPkici+AewA8DeAnRX5j1FnyVL7ipxNi5YuIKscXvkKXi2jmudQRghZZVdV5AObF7jvDub4RwCzP874N4NsJL7tf+G5St7MBauxYYP16YMuW5G3tY6NHmyDGyhcRVY5v2jE0QLHy1fUYpakSbIDaaitzuXlz8rZ9feZy3DhzycoXEVUOpx0pBcMXVYINUFtvbS7TwpetfNnwxcoXEVVO0eEr7yr5VGkMX1QJrHwRUVcpetqR4aurMHxRJdgAZQNVSPgaO7b+uURElcHKF6Vg+KLSqearfHHakYgqj+GLUjB8Uen6+sxY09MDjBlj7uO0IxF1tGaWi2D46noMX1Q6W7kaNQoYOdJcD1lqgpUvIqos9nxRCoYvKp2tXPX2RuGLlS8i6micdqQUDF9UOl/lK0/4YuWLiCqH4YtSMHxR6WzlauRIYMQIcz2k4Z5HOxJRZRU97eh7PepYDF9UOhumRo7MV/myzflp/WFERKVopvLley4rX12F4YtKZ8PTiBH5wtfo0fXPJyKqDB7tSCkYvqh0vvAVcrQjK19EVFkMX5SC4YtKx8oXEXUdX48WG+6phuGLSmeDVqM9X2nbEhGVgkc7UgqGLypd3soXpx2JqPI47UgpGL6odG74CllqgtOORFR5nHakFAxfVLpGK18MX0RUWb6w5AtVPlznq+sxfFHpGm2457QjEVVWM2t1sfLV9Ri+qHR5l5rgtCMRVZ6vUsXwRTUMX1Q6W+VqdNqRRzsSUeWw8kUpGL6odDy9EBF1HV9YCu3b4umFuh7DF5Wu0Z6vUaPM5cAAxyUiqhhWvigFwxeVLs9SE6pR+HK3Z/WLiCqFPV+UguGLSucLXzZgxdnxrKfHLJnD8EVElVR05YtLTXQVhi8qnRu+hg+vvy/OhjK7HcMXEVVSMz1frHx1PYYvKp17tKMNVUmVLxuy4uGLRzwSUaVw2pFSMHxR6dyjHbOmHd1+L/sc9zWIiCqBDfeUguGLSuebdsxb+WL4IqJK8YWlpIEtjuGr6zF8UenyhC/2fBFRR/CFpcsuA66/Pvu5DF9dj+GLStdI+LKhi+GLiCopqbl+1qzs5zJ8dT2GLypdnqMdOe1IRB0hKSyJZD+X4avrMXxR6fKs88VpRyLqCElhaVjAn12u89X1GL6odHaZiJEj8087hpyOiIio7ZLCUkj44rkdux7DF5WORzsSUdcpuvLF8NVVGL6odGy4J6Kuw54vShEUvkRkpogsFZFlInKa5/FeEbmq9vjtIjK5dv9hInKXiNxXuzzYec5+tfuXich/iIT8RFI34umFiKjrNDPtyPDV9TJ/CkSkB8AFAA4HMBXAR0RkamyzTwFYq6q7AjgXwNm1+1cDeJ+q7gXgOACXOc/5LwDHA9it9jWzie+DOhinHYmo63DakVKEVL6mA1imqo+q6mYAcwAcEdvmCACX1q5fA+AQERFVXaiqK2r3LwEwqlYlezWArVX1NlVVAD8H8IGmvxvqSG6gsqGqv98//nDakYg6AsMXpQgJXzsCeMq5vbx2n3cbVe0D8AKA7WLbfBDAQlXdVNt+ecZrAgBE5HgRWSAiC1atWhWwu9Rp3EAlAvT0mNu+qn1S5Sv0rB1E7cTxawhjzxelCAlfvp+U+E9G6jYi8gaYqcjP5HhNc6fqRao6TVWnTZw4MWB3qdPEq1lpU4/xnq+saUqiMnH8GsKK7vniOl9dJSR8LQewk3N7EoAVSduIyHAA4wGsqd2eBOA3AI5V1Uec7SdlvCYNEXkCVVJQ47QjEVUKpx0pRUj4uhPAbiIyRURGApgNYG5sm7kwDfUAcBSAG1RVRWQCgP8F8DVVvcVurKorAbwkIm+uHeV4LIDrmvxeqEMlhS9foIpPO7LyRUSVVPS0IytfXSUzfNV6uD4HYD6ABwBcrapLROQsEXl/bbOLAWwnIssAnATALkfxOQC7AviGiNxT+3pF7bHPAvgpgGUAHgEQcKp36kaNVL7Y80VElVb0tKPvPupYw0M2UtV5AObF7jvDub4RwKBTtavqtwF8O+E1FwDYM8/OUnfKE6jy9IcREZWmmWlH33NDw5dqWHWNSsUV7ql0eSpfSdOO7PkiokopuucrJHz19QH77Qccc0z2tlQqhi8qXZ4+Lla+iKgjlDHtuGQJsHAh8ItfZG9LpWL4otI1s9QEe76IqJK4zhelYPii0vFoRyLqOkWHr5DKF3u9OgbDF5WO63wRUdcpo+eLOgbDF5Uuz1QiTy9ERB2BS01QCoYvKl0z63xx2pGIKokr3FMKhi8qHacdiajrcNqRUjB8UemKWOeLlS8iqhROO1IKhi8qnQ1UIdUsLjVBRB2hjMoXj3bsGAxfVDqeXoiIuk5S+AoJUb7nhvR8ua/NSlmlMXxRqQYGolOR2Q+EPL0QEXW8pPAT8kmx0cqXG9DYoF9pDF9UqnjVy73OyhcRdayk8JPUC+ZqtGrlvmfI+1BpGL6oVI2GL/Z8EVGltSJ8ZYUyVr46BsMXlSotfPH0QkTUsZKCUjPhKytQua/NylelMXxRqYqadmTPFxFVCitflILhi0oVD1PudZ5eiIg6VjwoTZ9uLpctyw5gRYQvVr4qjeGLShUPU+51nl6IiDpWvPI0dWp0/Zhj0p/baPhyAxcrX5XG8EWl4rQjEXWleFAaOTK6fuWV+Z5rZQUqVr46BsMXlSpv+GLDPRF1hHhQcnsrsrDnq+sxfFGp8h7tyKUmiKgjxMOPW/nKwp6vrsfwRaXiIqtE1JXiQakdlS/2fHUMhi8qlS985TnakT1fRFRJzUw7JgUn9nx1DYYvKpVvqYk8lS9OOxJRJZVR+WLPV8dg+KJSNXt6IU47ElElldHzxRXuOwbDF5Uq7zpfPNqRiDpCI+Hr3HOBf/7n5AGNla+uMTx7E6LWafRoR67zRUSVljbtOHq0/zknnWQut97a/zh7vroGK1/UUt/8JnDyycmPN7vOV0jP1xNPADNnArffHrTLRETNiwcld5AbNSr9uRs3+u/n0Y5dg+GLWubJJ4EzzwR++EPg2Wf92+Q92rGRnq//9/+A+fOBN785fN+JiJoSD0o9PdH13t705yYFJ67z1TUYvqhlrr02ur5ypX+bdpxeyB2Dnn8+eTsiosLEA5QbvrL6v4oIX6x8VRrDF7XMsmXR9RUr/Nu04/RCmzZF1598Mnk7IqLCxMPPMOfPbVb4SqpasfLVNRi+qGXWro2uZ4WvvOt85en5cqtu7j4REbVMPCjlCV+NLrLKpSY6BsMXtYw7xdfItGOeox3Twpcb/Bi+iKgtfJWvHXc013feOf25c+f67+e0Y9dg+KKWcYPO00/7t/Gt89Xo6YV845JqffBjzxcRtYWv4f6ii4p9zThOO3YMhi9qGTd8JYWeZhvuhw2Lqvm+D3obNgCbN/v3iYioZXyVLztYZYWoJFxqomswfFHLuIFr/Xr/Ns2eXihr+/j7svJFRG3h6/kSMdcbDUZcZLVrBIUvEZkpIktFZJmInOZ5vFdErqo9fruITK7dv52I3Cgi60Tk/Nhz/lJ7zXtqX68o4hui6nCrTOvW+bdp9mhH97qvRyz+vqx8EVFb+JaasOGrVZUv9nx1jMzTC4lID4ALABwGYDmAO0Vkrqre72z2KQBrVXVXEZkN4GwAHwawEcA3AOxZ+4r7mKouaPJ7oAratMlM+VmNVL7iYWpgIBpP3CVz8lS+GL6IqC18047tDF+sfFVaSOVrOoBlqvqoqm4GMAfAEbFtjgBwae36NQAOERFR1fWqejNMCKMhJD69l1X5Cllqwg1qdgxzn+sLX6x8EVEpfNOO7PmimpDwtSOAp5zby2v3ebdR1T4ALwDYLuC1f1abcvyGiPvnlDrdyy/X385T+UoKU75t3dtplS/705W0H0REhUqrfLHna8gLCV++UBSP3yHbxH1MVfcC8Lba1zHeNxc5XkQWiMiCVatWZe4sVYM9L+y225rLpNCT1sOVFL7cKpm7fVrP18SJ5jIeColaiePXEOZbaoI9X1QTEr6WA9jJuT0JQHy98v/bRkSGAxgPYE3ai6rq07XLlwBcATO96dvuIlWdpqrTJtq/oFR59pQ+29Xqn0U03PuCWtr2QBT6GL6oDBy/hrCyl5pg5avSQsLXnQB2E5EpIjISwGwA8eV35wI4rnb9KAA3qCb/lIjIcBHZvnZ9BID3Alicd+epumzla8IE82Fv40b/WJAnfCVNO4b0fL2idiwtwxcRtUUrph1Z+eoamUc7qmqfiHwOwHwAPQD+R1WXiMhZABao6lwAFwO4TESWwVS8Ztvni8jjALYGMFJEPgDgnQCeADC/Frx6APwJwE8K/c6oVDZ8jR4NjBljKlAvvwxstVX9dnmOdsyadmTli4gqI22dr0YrX+z56hqZ4QsAVHUegHmx+85wrm8EMCvhuZMTXna/sF2kTmSnHUeNAsaONSFo3bp84SvvtCN7voioMuJBSSR/+JoxAzjlFOD004GHHuLRjl2EK9xTS9jK16hRwLhx5rqv6d5Xzco62jFe+UqbdrTvyWlHImqaKjBrFnDoodlBKP64av6er912A446Kjy0uYFr7lyurVNhQZUvorxs5au311S+AH/TfRE9X2nTjvY9J0ww2/X1mQpZPMAREaXavBnYZx/g/tr64s8/D2yzTfL28crTwED+ni9byQoNbe7rXn45cN99wKJFYe9FbcXKF7WEW/kaM8Zcd1e8t1o97RjvPQNY/SKiBtx8cxS8gPqVnn3Swldo5cu+Rmhoiz9+771h70Ntx/BFLWFDT2+vCWBAVA1zFbnOl6/ytXmzuRw5kuGLiJoQ/3SXFb7iAWtgIP+0o618hYY2Ntl3DIYvagm34b6311zf6DnJVJ6jHZMqX2k9X+70J8MXETXMV1pPU+S0YyM9X1RpDF/UEu60o6185Q1fRfR8sfJFRIWIDzBZQScelPr78087NtPzRZXG8EUt4Vacig5feU4vxMoXERUiPsDkDULt6PnitGPHYPiilvBVvnw9X3mWmmjk9EK28sXwRURNabby1chSE5x27FoMX9QSvob7Vk07hvR8jRxpjngEkk/yTUSUKF75ynvkYVbPly9YMXx1LYYvaokiGu6LONrRN/3pq8AREaVqNnztvHN6iEoLX6EVs2eeSX+cKoOLrFJLNNNwb0+Bplp/dHYj63y5DfdpIZCIKFXeni/7+OWXm4Fuxgyz6GnSc31hLk/P1wMPAHPmpO8TVQbDF7VEaMN9WqDassWMWSNHmvsaOb0QK19EVIhGK1/TpgGve525nnfaMR6+0gLftdem7w9VCqcdqSXyNtyHNNE3u9SErXwxfBFRbo023A9z/symhSjf6+Xp+bJNrdQRGL6oJdyG+7w9X+5td7xrZNrRrXwxfBFRw+wnOSu08hUavprt+WL46igMX9QSbsN9SM9XSBN9s6cXYs8XETUsPnCEhi/3NERpISpk2jHtPe1aOtQRGL6oJZppuHdvh1S+2PNFRC0XH8BCG+59la+05npXnmnHnp70/aFKYfiilmhmhXv3Nnu+iKgS4gNHI5WvRqcdQ8JXfFqUKo3hi1qilQ33oacXUmXPFxEVJO+0o6/ylTbt2Gzli+GrozB8UUu0s+E+adqxvz86o0dPT3oFjogoVaM9X6HTjmmBzL5G2nsyfHUUhi9qidCG+zzna8w77ehWvdxLVr6IKLf4wBF6qp/QacdmK1++w72pshi+qCXyntuxiKMd42OP2+9l9wVg+CKiBhQx7cieL6ph+KKW8J3b0Rd68lS+8mzrvh8rX0TUtCIa7vMuNdFI+Bo7Nn2/qBIYvqgl3IZ7W3nyfTArouE+qecrXvlizxcRNSzv6YWKWGoi3vNlzw2Ztn9f/nL6flElMHxR4VT9DfdpRzvmmXZk5YuI2i4+wBTd8xWyyOqXvwysWeN/P/tp011s1X1vqhSGLypcX58ZR4YPN0cZplW+ipx2ZM8XEbVMoyfWLnKpCQBYscL/fu6Ad/HF5vqECen7SKVh+KLCuVUv97Koylfo6YXilS+ucE9EDSvyxNqhS034wlcSGw5HjgQ++EFz3bfyNFUCwxcVzm22B1pf+Qrt+eK5HYmoYY1WvopYasINcEnTnXbAGzEiOtWQfT5VDsMXFc5ttgeSK1+q0dgQPy0Ze76IqFKKrHyF9nwdfHD985K2A+o/bTJ8VR7DFxUuPu2YVPmy48Lw4YOr6kWs88XwRUSFKaLhPs9SE9/6FnDhhYNfI4k77cjwVXnDszchyidt2lE1GkeSphHd+4qcdmTPFxE1rN1LTZxySjRo5al8jRiRPB1AlcHKFxUuXvkaNsxfnUqqZAGtnXZkzxcRBduwwQwaeacd087tGFL5sgOW+7yk5wL1nzbd98zaTyoFwxcVLl75AqJxxJ16DKl8hYQ1LjVBRC2hCmyzjVmyoYiG+zxLTSRNNSZNJbrTjgCnHiuO4YsKF2+4B6LxwA0+SZUs974iTy/EaUciymXzZjNgbNoErF9f/1hWz1cRS01Y7vZJYSr+aZPhq9IYvqhw8WlH97qv8pV32rHR0wvZ7bZsYSWeiAK4wWXt2vrH2rHUhO+xpD4ut+cLYPiqOIYvKpxv2rGIylezPV8inHokohzcQSUejtqx1ITvvZLCV3zakU33lcbwRYULrXzlbbhv9vRC7n4wfBFRprTgkha+3CDV6FITSY9x2rErMHxR4UIrX2kN976pxGZPL+TuE8MXEWVqNHz5phzd2yFLTSQ9xmnHrhAUvkRkpogsFZFlInKa5/FeEbmq9vjtIjK5dv92InKjiKwTkfNjz9lPRO6rPec/RHj69W7hq3z5FlotquE+tOfL3SeGLyLKFC+nu0IqVcNif2IbnXZ0A1RS+Ip/OmX4qrTM8CUiPQAuAHA4gKkAPiIiU2ObfQrAWlXdFcC5AM6u3b8RwDcAnOJ56f8CcDyA3WpfMxv5Bqh6fEc7+kJPow33jfZ8ude51hcRZSq68pVnqYmkx5LCVHyAZPiqtJDK13QAy1T1UVXdDGAOgCNi2xwB4NLa9WsAHCIioqrrVfVmmBD2f0Tk1QC2VtXbVFUB/BzAB5r5Rqg60qYdm6l8NbvOl7tPrHwRUaZme76SKl95l5pwH8uqfNkBkQ33lRYSvnYE8JRze3ntPu82qtoH4AUA22W85vKM1wQAiMjxIrJARBasWrUqYHepbGkN976jHZttuE+adkyrfDF8UTtw/Opwrer5yqp8XXZZ8mOh4YuVr0oLCV++Xqz4T07INg1tr6oXqeo0VZ02ceLElJekqgitfOU9t2Peace0ni9OO1I7cPzqcGnhK2RpiEZ6vvbeGzj6aP/rAZx27BIh4Ws5gJ2c25MArEjaRkSGAxgPYE3Ga07KeE3qUHkrX62admTli4iaUvS0Y8hSE/HnxN+Lla+uEBK+7gSwm4hMEZGRAGYDmBvbZi6A42rXjwJwQ62Xy0tVVwJ4SUTeXDvK8VgA1+Xee6qktNMLNbPCfd7TC7Hni4iaUsZSE74D/xvp+WL4qjRPzaGeqvaJyOcAzAfQA+B/VHWJiJwFYIGqzgVwMYDLRGQZTMVrtn2+iDwOYGsAI0XkAwDeqar3A/gsgEsAjAZwfe2LukDaibVbUflizxcRtUQj4UsVmDPHXG9k2jGr8hU67ciG+0rLDF8AoKrzAMyL3XeGc30jgFkJz52ccP8CAHuG7ih1jtB1vlq91AR7voioKY30fH3ve8DXv26uJ1W+7PPd22mVL047dh2ucE+FS6t8tWKR1Tw9X5x2JKJgjVS+vvGN6HpS5QsYHN7sbYavIYHhiwqXVvlq57SjDXqcdiSihjQSvtyw45tCtM480/96vueEnNuR044dheGLCpe2wn0rGu7tB7y+vvoxygYsnl6IiBrSSPhyw1PaWfPOOqv+djOVr4GBaBrTvv+YMeby5ZeT94FKw/BFhQs9sXajla/49sOGReONO0bx9EJE1JRGwtf48dH1PKcsbmapCd/guNVW5vKll8L3gdqG4YsKl7bOVxEN92nbu31fPLE2ETXlj39Mfiyp4d4doDZsCH+v0IZ737Sjb1qA4avSGL6ocK2ofPX3D66qu3x9X1xqgogadv/9wHnnJT8e0vO1bl34+6VVvrLW+TrgADBYUQoAACAASURBVHNpezAAhq+KY/iiwuVdaiIkfKUFNd/27nux8kVEuT38cPrjSeErbfHVkNfzVb7cQOcLX4sWmUt3gLXh65ZbgIceamyfqGUYvqhwaQ33jZ5YO21bd3t32pFLTRBRw9xPbT5Fh6/Qhvv4tKNbFXOv2/B1xRXA7rs3tk/UMgxfVLjQE2vnmXZMq5L5tnffi5UvIsrNF74OPhjYeWdzPannq9F1tdKWmkhruHffz71uw5eVdiJwajuGLypc6Im18zTcZ1W+2PNFRME+/nHg6KPTt/GFr113Bd79bnM9qZeq1ZWvtPDlbhcPX3ma/6nlGL6ocK2ofLHni4gKMTAAXHopcPnl6UtJ+MLXpElRZerUU4GVK/2v34jQhvt4ZS2p0hYPX3ma/6nlGL6oUKrhi5vmabgPnXa02/X3my+R+ucwfBENcW4FKO86XpMm1Vem5s0bvE2z047NVL5c8fC1fn1j+0UtwfBFhXKDl/sBLq3y1YppR7fq5Y5lDF9EQ5y7wnLe8DVzZv3A5uujss/bYw/gxBPD96vRRVZZ+epICXUEosb4+qzc23nX+bKVrLwN93n2g4iGkNDKlxtqPvAB4L//G3jFK7LDl33e4sXp53aMa3SR1aTwNW5c/W1WviqFlS8qlG+ZCSB9na88la/QaUdfvxfA8EU05Lnhy12bJs4NPD09JngB6YHKDWN5Ti3kPtf3vLRFVpPC1+jR9bdZ+aoUhi8qlK/ZHmi88pV3nS9WvogoVSM9X25IS6t8uctF5A1fjS41kfQ92BNrWwxflcLwRYXyLTMBNL/Cfda0Y1rPl4vhi2iIa2Ta0d3ODVXxvjD7nDzTjVaji6yGVr447VgpDF9UqKRpR9+JtRtpuC+q58vtuSWiIaSRhvtGKl95hTbcx6dKOe3YkRi+qFBJocd3Yu1Gzu0YenohVr6IyCu058sNNaHhyz7HPcF1qLSG+yJ6vlj5qhSGLypUUZWv+DQij3YkokK0q+crhO+8jFmVL3cQBZLDV/x1WPmqFIYvKlRSw72v8lVkw31SzxfDFxHVaTZ8pTW/28dCK1/ua4UuNREavuJY+aoUhi8qVCMN963s+eK0IxHVaXba0XcC2fhzQitf7nuENtwnvWee90ry+OPAM8+EvR41hYusUqGyph0brXzlPb0QK19E5NVIw/2HPxxdd0NMvAqVd9rRF76KmnZMey+fdeuAKVPq94VahpUvKlRSr5UNR/aci0CxDffxacekytfw4WZsc/eDiIaQvNOO22wDnHyy/zlJVajQaUd3EGr1tGPWdqtXh70OFYLhiwqVVPkSGdx0n2epiaIa7t39YPWLaAjKu87XYYfVDzydWvlK+16B+u+Rla+WY/iiQiU13AOD+74aabhv9vRCAMMX0ZCW9/RC8SpWWuUrb8N9aOXLPU9jPHxlhSrfe/m47xv6mtQwhi8qVFLDvXufHa8aabhv9vRCvv0goiEkb+UrXolyn5NUhSq68jV/PjB+vLneqoZ793EOji3H8EW5pf1eFlX5sh8c+/vNmFTU6YWA7PDV388PfkRdZ8sW4Oyzgbvvju4L6fmKV7HSQkoz045pla/p04HFi831Vk07uq/DU4C0HMMX5XLlleZ8rf/93/7HG6l8+QKVSH0AK+r0Qr79cG3eDLz97cDOOwMvvOB/LyLqQBdcAJx2GnDTTdF9IdOO8SAV0vPVyLRj2lITgH+l6vhrnH562HtlPc7w1XIMX5TL2Web8eWEE/xr9iU13APJla+QqcSiTi8EpIevn/8cuPlmYOVKU+knoi5x332D72t22rEV63wlPde3WKL7GtOmAd/+dv1jf/6z/7183O+L4avlGL4o2IYNwAMPRLcffHDwNnkqTnmOYCzqaEfffrjcGYnf/97/XkTUgXxH8DUy7eg+55e/BB56aPBzip52BPynCXFfY9ttBz/n4IOBOXMG73fWvrDnq+UYvijYI4/Uf+jyha+yKl9F9Xy539OSJf73IqIO5C7XYDVS+YpXkC66aPBjzUw7hlS+3CCZ9Z5u/0bovrDy1XIMXxTsscfqby9dOnibtIb7eOjJs3xEuypf7vcU/36JqIP5Kl/HHpscNJKqWPHANnly9nOS5Kl89fSY11Wtf15W+HJXuA7dF4avlmP4omA2jIwZYy4ffnjwNmkN9/HKV+gRjG74KrLnKz6+vPwysGKFea3hw4FVq3guWqKukbRw6B/+4L8/5GhHIBoQ056TxA1yWZUvwN90H1r54tGOlcLwRcFs+DrgAHPpO/9q2rRjnhXugfqwlham3NdopvJlv58ddohOcfb44/73I6IOkxS+kipNIQ338dvNNNxnVb4Af99XUdOOaQcSUOEYvijYk0+ay+nTzaUvfKWFnvi40YrKVzM9X/b7eeUro5mEJ57wvx8RdRhfzxeQPKiENNzHbxdxeqGQ8JWn8hUfHEP2hZWvlgv6CRGRmSKyVESWichpnsd7ReSq2uO3i8hk57Gv1e5fKiLvcu5/XETuE5F7RGRBEd8MtdY//mEu3/Qmc5m38pW34b6ZacdmKl+vfKX5AszUIxF1saSgFLLOF+CvfLWi4R7wh6+s5lk23FdSwv9WRER6AFwA4DAAywHcKSJzVfV+Z7NPAVirqruKyGwAZwP4sIhMBTAbwBsA7ADgTyLyOlW1/8vvUFWeSr1D2HAydaoZH557zoQdNxAV2XBvxxk3fCVNOxZd+dp6a3Od4YuoSyRNO6ad6gLIF75a2XAPNNfzxfBVKSE/IdMBLFPVR1V1M4A5AI6IbXMEgEtr168BcIiISO3+Oaq6SVUfA7Cs9nrUgdyeqO23N9fj4aQVDfduz1foUhPNVr4mTjTXGb6IukTStGNS+Eqadnzzm+tv+8JXK1a4Bxrr+bID7K23AuvWhe3LddcBn/3s4AVdqTAh4WtHAE85t5fX7vNuo6p9AF4AsF3GcxXAH0TkLhE5PunNReR4EVkgIgtW8S9haV5+2fze9vaa87vaabn41GNIw/2mTfVHS7ey58sXvuy+hYSvZ5/1vx9RCI5fFZJU+Uqq8iRVvs48EzjvPODoo83tohru7fnMbNndp5GeL/doxwMPTH5t9/uYOxe48ELgkkuSt6emhPyE+GJ4/Kc4aZu05x6gqvsCOBzAiSLi/alQ1YtUdZqqTpto/yJS27nBRCSqfK1ZU79dSMP95s3R73lPT/IHvUbClx2T7H7kmXa038v227PyRcXg+FUhecNXUhVrzBjgC1+IDokuatrRDjZ2cPVpJnwBwMKFYftiPfdc8vbUlJCfkOUAdnJuTwKwImkbERkOYDyANWnPVVV7+SyA34DTkZXmhi8AmDDBXK5dW79daOUra+kI97HNm7N7vtz+MPsc9z2T9sNlv5cJExi+iLpO0tF+eRdZtXxHETbTcG8Hm7SQ3sy0Y559sUK/D8otJHzdCWA3EZkiIiNhGujnxraZC+C42vWjANygqlq7f3btaMgpAHYDcIeIjBWRrQBARMYCeCeAxc1/O9Qqq2uHRdgPZTZ8Pf98/XZpDfe+dbt84chyK19ZPV/xPtRGKl/2e2H4IupCSb1deacdLV/4Cql8nXji4PcAwsLXVluZy5deGvwaIZWvNL7wFRrcKLfMf1lV7RORzwGYD6AHwP+o6hIROQvAAlWdC+BiAJeJyDKYitfs2nOXiMjVAO4H0AfgRFXtF5FXAviN6cnHcABXqCpPY1xhtvq83XbmcpttzGU8fKU13OetfOWZdoxX4xtpuLffyzbbMHwRdZ284SureT4tfKUFnvPPNyfj/uMf84cv+6nX9ocBrQ1frHy1TFCsVdV5AObF7jvDub4RwKyE534HwHdi9z0K4I15d5bKEw9fSZWv0HW+Wh2+Gllqwq18jR9v3mvdOmDDBmD06OT9JKIOkBS+8i41YaVNO2b1fNmBzA08dnohLXyNH28u3YG3qGlH37QsK18twxXuKUhS+HJ7vvr70xdObUfPV0jlywbD+Adet+dLhNUvoq4S/4W3R/q0ovKVFb7ia2+phjXct7vyxfDVMgxfFCSk8rVhg7kcPdp/BKMbkNJ6sqw8PV95Kl82fNn9Bcz4u2mT2d4+zvBF1EXiFa7zzjOXRTbch67zFQ9fGzaYQau3t/5E3XG+ylfoCvdZfOEraW00ahrDFwUJ6fmyY1jSFF0ZPV++6U+7f+6Y6/Z72eDI8EXUReLhK6kEbrVy2jEevmwly4arJI1Uvpo52pGLrLYMwxcFyVv58snb89XItKMdX9Ma/+3+uZUvt9/LYvgi6iLxkGUHh3Y33LuP28Dz4ovmMit82cfbNe3I8NUyDF8UpIjw5S4HkXepiazKl/vaquk9XwxfREPQ+vX1t5NOdWFVufKVp+E+NHz5Gu4ZvlqG4YuChDTc2zDjm+oD6qtTeacd8/R89fWZD6A9Pf6Ku6/ny222t2zf62qe+p2o87lrYwHp045r1gAXXGCuF73OF1B/yh8g7NRCQGOVL047VhLDFwUpsuer1UtNpFW93P1L6vmybPjiGTaIOpw76ADAttumh69vfSu63oppx3jjfOi0ow1ndnuA044diuGLMr38shmf3ANxxowxY8+GDVHYCe352rQp7GjHRpeaSFtrzN2/rGlHGzRZ+SLqcPGq10c/mh6+3F/6Vkw7vuIV5vLZZ81l6LSjb6qUla+OxPBFmdyqlz0SUGTwgTdFN9y3uvKVFb5Y+SLqEuvWRddPPBE455z0hns3PLWi8mVPkmtPmmsrWVnTjr4VoputfK1eDaxcyfDVZgxflCk+5WjF+76yer6aWWoiT89XVuWLPV9EQ4ytfO2xhzm9z+jR6Q33bviqUuXLFxizBrys8DV5MrDDDvVTmRbDV8swfFGmpPAV7/tqpPKVdrRjo9OOzVS+3J4vTjsSdQkbvuyJqYH0acdGw1dow3288tXMtGPWwBufdlStv22PAr377sHPtUcwxadtqWkMX5Qpq/Jlg0uZi6z29JivgYFoLMmqfG3cGI1DaT1fzz03eLwiog7STPhq5bSjrXwtXGgud945/Xm+aces8BUPgnYwdfcXAO64Y/BzN28G9tnHTIe6h7ZT0xi+KFNo+Cqz58t9LVs9T6p8DRs2eAzzha9Ro4Bx48zY6qvIE1GHsD1f48ZF97Wi8pV32vGZZ0wwvPlm85xDD01/njvtaD8RZg288XO9ueHLnVb0LWh4xRXA4sXm+p13pu8b5cLwRZnstFv8fK/N9HzlObfjpk3ZPaXua9mglLQf7mN2n309X0B99YuIOpSv8tXKhvus8DVunAlLGzaYcNPXB7zudfV9Dz7DhtV/KgWyw1ecG7iSFpj18Z2wlxrG8EWZbPiyK75bRfR8hSw18fLL0e2033+7vR1n0/rJ4n1fvp4vgE33RF3h0UfNZei0oxu4WnFibZGo+vXUU+YyNDzFQ2Pe8JVU+aK2YviiTLYanVT5alXPVzxMpU05utuHVL7iC636ph0BNt0Tdby+PuDCC831XXaJ7o+fkyxJ0ie+ZqYdgajvy4avtE+LrnjTfTPhK0/liwrF8EWZsqYdW3W0ox1j7OtnjS32tfJWvlSj94gfbMS1vog63HPPRY3t//qv0f0i/gb2+O2k6lAzlS9gcOUrNHzF95nhqyMxfFGmpGnHePiy04NFrXBvXyc0fOWpfNmV+tevN199fea++P5w2pGow61ZYy53333wIJI09eiGkqSAUnTlK20wdBU57cjwVRqGL8oU2nBvDyhy2ypcdjkIoL6PK4kdF+3rp4Up97VCKl/2oKf165OnHAE23BNVmrtYXxIbvrbddvBjSU33brUrT+XLrtuV1TgPRJWvhx+u35csRU475un5cr9PahrDF6VSTe75ijfc+47mjotPDYaEL/v6oeHLrleYtr3dx3XrkpvtAVa+iCrrt7815erzz0/fLi18NVP5sp8k3VDy4IPm8vWvT98nANh/f3N5333mMm/lq93Tjr6g9vWvA0cfzYUQG8DwRansSbVHjYqm6qz4tGNI+LLji31O/DVd8Q94WeHLvpZ97bQPkmPHmks3fPkqXwxfRBVl+7c+//n07ULCV1E9Xw88YC732CN9nwDgve8Fdtwxup235yvPtOOvfx1dDw1fJ51Uf9v37/Dd7wKXXw48/XTy65AXwxelcvu94gf9NBK+7Lhhx8OQ8JV0O86+ln3tkGnHdeuS1/gCOO1IVFnu6uxpzj7bXKaFr9/+tv7+PD1fL75oBp3Nm4HHHjP9Xrvumr1fIv51x7K4gbG/34QpkfTK2ZFHAjNmmOu+db58R3QecUT97Xj4sqcSAcL/L+j/MHxRqqQpR6A+fKlGv4tp4ctWnGyoKzJ82Q9+dp/T9iO054uVL6KKCpnq+vvfo2qU77yJ//iHuTz11OjTI1AfuHbayf/aNny98IL5lPbwwyaE7LBDeJByB8BGph3dqlfWIqjxxVmBKFBNnRrd9+Mfm4HRHhAQ39ZyB8W03rukYOZbX20IYfiiVEnN9oAJQ6NGmd/JDRvCKl82fNmAlBa+4lX0rJYG+1r2te17+YT2fLHyRVQhfX3Ak08Ct9wSVm2xPVhA1AzqsktQAMB110XXbfh6//uBj33M/9rxgHX77eYyKaz5uANg3srXxo35+r3ctX4s+33utlt039NPm/2KHzmVFr6SgtSKFWba5Iwz6u+/+mqzz5demr3fXYrhi1KlhS+gvvplw1da6LGP5en5SrodF+/5Cg1fIdOOq1ezp5SodIcdBrzmNcBb3xodWZjmoYei65/9bPq2blCzoeTMM5PX7IqHk+XLzWWe8OWGpmYqX2kDqeUOepb9Pnt7zXkc99gDOOEEc1/8+3MrZkBY5es//9NMyX7rW/X3H3usufz4x7P3u0sxfBEefxy46qr6KXzLVpHia3xZvvAVUvmyWhG+kt7L5at8+cLX6NHmdbZsiY7QdN1/P3DttVwuh6jl+vuBv/wl33OWLjWXl19u1vlKc//9wN/+Zq67oSRJPJTZJSMmTQrfv0YqX27DfcgnXsuGqaTw9ZGPmH8Du//x12xk2jHpEyvPE8nwNdTdey/whjcAs2cDhx8erRFohVa+Vq0yv8fDhqWHpDzha/jw+rUKywhfQPIphu66C/infzK9rL5/O6LKWLWq88/jt2JF/ufYBUynTPE/7k57/frXwIEHml/skPAVZ8NXqytf7uk5sgYvlx303E+RaatdxxeKbWTasYqefhpYsKDsvWD4Guq+9KVowdO//Q34zW/qH7f9qHY9wDjbJ2WPNB43Lv1DTbwqltaqIFIfuIoMX76lJpLWRbR9p/FZjtNPj/7tbrwRuPLK9P0jKsVTT5lf4Le+tew9ac4TT+R/jj3dha/ZHjDTX3Pn1t93881Rk2dIRcmy045J0wQ+jVS+tt7aXL74Yr7wZStfNnz190d9cCHvHQ9fbiOsW/kaGAA++UngoouS/xiUWfmaNMl8an7ssfL2AQxfQ9p995nQsNVWwHe+Y+77yU/qt7Ghyl2OxmV/5+0HzLQpRyBf5QuoD2ehDfdJ7+UK7fkCzMFLQP1SNitWAPPnm0Boj2Q/77z0/SMqxe9/by7vvLPc/Qi1ahWwePHg+xsJXzZo2MDiE//Udd11ZmCYMSP5U6eP/aSaFPR8Gglf9vVfeKG58HXqqcBXvxr+3vHwZd8bqK983XIL8LOfAZ/5TPJrVWHa0e3xKwHD1xA2Z465/OhHgeOPN1XmG2+sPyio7PDVqsqXW4HPGr/s9+6Grz//2Vy+4x1mjcettjKV7GXL0veRqBBz5pj57pDT67h/GC+5BPjyl6t99MhrXgPstRdwzjn1c/lLlmQ/95vfBPbbL2pgtZWvpHOeAYN/8e+911wecUS+kGD3NU/4amTasdnwdcEFwL/9G3DuudFjIeEr3nDv/qGwP4f/8R+m6mUl9WKUFb7cn/uQk59bCxcCH/pQYx8AEjB8DVGqwC9/aa5/6EOmp+uAA8zv1w03RNuFhq9HHzWXtj8qiRuIhg+Plp5J0qrwZfdzzZrmwtehh5rx065HePXV6ftIVIiPfMQc6XHVVdnbukeDfOITwI9+FC2LUDWq0R/yr34V+PnPzfWBAeCyy8z1+fOT/3CeeSZw993m32ZgIOwooHjly06n5al6udpZ+bIBKOQ9bfh6/nngrLPqHwsJfvHKVzx8qQJf/GL9J1BbDQRMELb/tmWFL/fDSjxMpnnrW80fzOOPL2xXGL6GqEWLTH/oxImmxxQADjnEXNoDijZsMFNyI0YkN9zbRaPtQUVJ21luIAo5OjpP+IpPS6aNt+7iqXY8SOr5suHL9vuqAn/6k7l+6KHm8kMfMpchfwuJmmLLzEC00Oczz5jeJV9Fy3cYs60IVU18Qb2bbjKXS5eanqoddjC/dPGm9viaX1u21B8JmFblSPrFz9O75Wpn5ctW6fI03PtkDa6ACV/uz5cbvjZujPo3XI88El3fdlsz8BZ14Mc11wDve1++n2V3qjTP82xzb4F9YgxfQ5Sten3wg9H4/fa3m0sbvmz/6A47DD7wxXrVq8xlaOXL/f0vOnzlqXzZ0Lh6tZl67O1NHr/skdf2b95DD5kq2MSJwJ57mvve+U4zJt57bxREiZr24ouDpxbd3i0bMN7yFlN+9aV/u16My/4x8Vm/Pt8fSFtZ27Il+xx/quYXKGna8/HH62/b0GSPTnvzm81gFA84GzcODmAh/V5AcjNpSOXLt3xFuypft94a9Y7kmXb0sY2taa64wgx6v/+9+SQar3y5QcuyR4AC0RTks88WU/maNQv43e/MVKfPE08A//7v5me5r89U4dzw5Vs7KEvWVE0ODF9DkDvlOGtWdP+MGeb3/957zQdQ+3uTdpqyV7+6/nbRlS83rGX1k7mvJ5LeoD98eP149epXJ48Hu+xiLu2/h616HXJIFEp7e4EPfMBct/+2RE1ZvdqsPG5L0pZ72K39Y2I//fzxj4Nfxxe+kk7Z8I9/mE9USau6x518svmjvnCh+QWYNCm9N+sHPzCBxe03csXDl/1kaMPXtGnm8i1vqd/u5ZfrKxkvvRTW7wWYX/wTTxx8f0jl69prB9+X9X6uRk4vZMOX+3+Y9Ok4dL9e8xr//e4+Pf+8ec/DDzfTAQsXRo9t2FAftKykn70ipx2Tfpbf8hbglFPMz9qHP2wG+Ztvjh7/znfyn5Py/vujA1iaxPA1BN1+u/k9edWroilHwFSW9t/fXL/ppuh36XWvS36tvOHLDVAhR3Hbyhow+FRjaa89Zkz277e7r/Hvw7XTTiZcrVxpCg223yv+N9EGWYYvKsSll5oqwW231f9hc0+J436SB/wVJbfvxko6Wem8eeaH/Jprkv8wXXst8J73AEcfDfzwh1Gj6Lx55vErrkj+nr7yFXN58snRfVdeaVZAVx0cvuwvse1Rs+Hr9NNN75r9pf/zn+v/CH/hC8BJJ5nrWZUvADj/fPM9u0LC1+tfb/bFGjkyXyO3+wkxb+XLlbWALJAeviZP9t8fuk+XXFL/75Cm6PDl+5lfsybqE/nLX8wabkB9lezJJ5OrZmkOPzz/czwYvoagn/3MXB5zTPTB0nrHO8zljTdGZ+ZwT/sVlzd8uWEqZCFoN3C5z/VxDwpIC1NWaPjq6YmqXw8+aP5tgKjfyzrssGjq0XekPFEu8+dH191P2241IR6+XLfeak4XY/umXEnhy53iTKpkHHmkCVqXXx7d//3vR9dDG5lt2PnoR825/8491xyF53rpJdPkbMPXfvuZyx13NM3de+1lbn/ta+boSJf9NwutRLmhprc3/HlugMrbz+QOgiEhERgcvg47DDjooOzn7bJLcqBMOqIqNHytXDk4OCeJn6vtpz8F3v1ufz/VTTeZn4n4NLn7wcD3IcH9mXf7QOKVuC9/2QTHJBs2+Kt3BQgKXyIyU0SWisgyETnN83iviFxVe/x2EZnsPPa12v1LReRdoa9JjXvqKfPhdO7cwS0Y69ZFbQLHHTf4ubbv68YbzTpgQPqHqgkT6n8/s3q+3A9YO++cvi1QH7iywleeKhlQv69ZYc3+G1x6qfl799rXDv6wOHKk+TsCAD/+cf1jqmbh7KuvNh/SO2lBaCrBwoX1U4iLFkXX3T8GK1fWn7TYPbLx299OXssoKXy5h9L/7W8mAJ54oukDO+GE5D4odyo0fgLrlSujBlLXrFnAb38b3T755MF/ZNeurV980DZrWnaB1MceM4t6+oSGGnfQ2GOP8OpMyEmtkxx4oPk3+OEPTT9biAkT6nsmPvWpsOdtvbX5f/B9mo5/CrfinzCzZJ0/EzA/e+7P6b/8C3D99dGRrdaiRSZUnnWWCdfvfCdw8cVmWtseaACYP3Lf/a7phdy0yXzoOPLI6HE31LlVY+sTn4gqEnGHH974Ua9ZVDX1C0APgEcAvBbASACLAEyNbfOvAC6sXZ8N4Kra9am17XsBTKm9Tk/Ia/q+9ttvPyW/DRtUf/lL1cMOUzV/6s1XT4/qccepPvGE2e7f/s3c/5a3+F9n40bVceOi54uorlmT/t677x5tf8896dtu2RJte8YZ2d/Xf/1XtH1/f/b2dtsDD8ze9vOfj7Y/55z0bX/wg/p/1xNO8G+3eLH5Nxs5UvXhh819f/iD6hveUP/8ceNUP/1p1b//XXVgIHtfhyoACzRjXOiUr+Dx64wz6n9YANVXvlL16qvNL9Db3z74cft1wAGq556retttqsOGJW8HqL7vfaoHHaT6hS+YH8Tbb1edNcu/7Xvek/5a8a8pU8zlV7+quu225voeewze7uCD019n8uTo+mWX+f+9jjwy/TU++9mwf/e+vug5hx4a9hxV1W9+M3reNdeEP68Zv/pV9J5//Wu+5x5wQPTcpUtVn3kmedu1a83PScj/+aJFqj/+cfZ29o+Q7+uRR1Tf8Q7VL35R9TOfyfczZ3+mp07N/zxAdfRo875H/fmNBgAACeNJREFUHaX6yU+qzpyZvG1fX/A/d9IYJuaxZCKyP4Bvquq7are/Vgtt33W2mV/b5jYRGQ7gHwAmAjjN3dZuV3ta6mv6TJs2TRcEnJPpvPNMBdR+a/FL330hj5W1zcCAqeT39ZkDRnp6zEEXdp2stWtNH6GdgRg1ynxgGBgwrRj9/aY6NXOm+ZA1MGAqW7bKFXfccdGHkH32McvmpHn3u80HF8DsY1bLg/1AecYZZlmeNBddFC2UnPGjWvfaBx2UfQ7eq682fZgAcMcd5owTSRYuBPbdN7r917/W98u5jj3WLEk0ebKp7tkK+KteZXpAH30UuOeeaPvddzetLFttZfpmN20yXxs3mi8R83/t++rpic6B2Yqlc1rxmjvuaPpgw95f7lLVacXvRfuFjl+49Vaz6B5gSrIrVzb+pv/8z+Yf+5xz/M3hVdXTU79A54EHml86H3eQsIYPN4MRYHrKZs8Oe9/ddzf9FuecY1aAD3H88VF1LmSQKkJ/f1Stevjh9KOi4h58EPj4x03Debxx1WfLluSDAfbe2/z7iphq4ZVXRuX/UOPG1Z/su0iHHBI16RbpqaeCT6CeOIb5Epn7BeAoAD91bh8D4PzYNosBTHJuPwJgewDnAzjauf/i2utlvqbz2PEAFgBYsPPOOwclzf33byz4dvrXPvuo/uhH9ZWqRx5RnT27fruzzkr/91u40HwIAJI/bLp+8Quz7b77Bv336NixZvvrr8/e9vbbo/0Osd9+Ztsf/jB72xUrwj/IDAyYD2SAqaqlVeHWrlXde+/otcePV/3e90xV0XrgAdVTT1WdOLH8n5t2f73pTWH/l6qq6PDKVyPjlw4MmF8mEVPO3nHHxv6he3pUb7jBvObKlaoXXmheb/Fi1U99Krky9ra3qa5fr3r44aqjRtU/ZsviX/mK6uWXq65bV18qt+8LmPKv7/V/9CPV5ctVZ8yov//661WPPVb14x83v0S9vdFj3/te8r/Xc8+pvupVZrtddjHf8+rV0XOffjrs313V7NePf2wqjKHsALjXXuHPKcJf/qJ6ySXtea+TT1YdM6b+/2vYMLMPruXLzf/BGWeo/vzn9T8T9v/Ifn3iE6r332+mbdxpCPfroIMG/wxmfb32teYP16pV5tLeP2mSuZwwQfWQQ8wlMPgPZNLXKadE12+7LfifLmkMC6l8zQLwLlX9dO32MQCmq+rnnW2W1LZZXrv9CIDpAM4CcJuq/qJ2/8UA5sH0mqW+pk/oJ8c5c6IPi/aTe/yy0cfK2EbEVLhGjDAfCPv6zIcR+zV2rKnK2KZwn3vuMQdN7buvWVIiywMPmPaAww7L3lbVtF7MmJHdlwWYquTdd5tp+ZDKym9/az5UhXy4e/ZZU+076qjkNgbXrbeaipPt202zbp2p8B1+ePayFxs3mkLDpk1mHcB4q4q1ebM5in7pUtPbOTBgqpSjRpnL3l7zb9TXl/zV3598Fo9mZAwNDZs40RwoF2JIVr4Ac0TYqlXmaLqlS80h7gcfbKoyo0aZH5yNG00Z8YUXzJFbw4eb62PGmOfvuy/wxjcmv8eaNeZ1xo0zfTMvvmh+8Y86qn4JhDVrTJl4n33MwnZLlgDTp0ePP/CA2f7FF00z5OLFZrDad1/g7383g8KCBeawYXfwWbfOHBo8dqwZDOJrKC1aZPp1dt01WlAvycaNZnBUjao0ixeb9wjtpWrUwIDpz5s+PXnB1k6nan5WVq40Bz8ccogZ3JIa9a3nnjNTMrvsYv6dbr7ZTAeMHz/432rJEjMF8da3mp/5Z581C1EuX25uz5hhqnZr1pipjb/9Lfq5GDHC/L/7qlF3321+N/be20w7jB9vGn63bDH7c9BBZiro2WfN79vEieZneuxY87zHHjP/t1OmmCOHJ0wIW5S2JmkM68ppRyLqDkM2fBFRV0gaw0KOdrwTwG4iMkVERsI01M+NbTMXwHG160cBuKFWbpsLYHbtaMgpAHYDcEfgaxIRERF1ncyJGVXtE5HPAZgPc5Ti/6jqEhE5C2Yucy5ML9dlIrIMwBqYMIXadlcDuB9AH4ATVbUfAHyvWfy3R0RERFQtAV0xgKrOg+nVcu87w7m+EcCs+PNqj30HwHdCXpOIiIio23GFeyIiIqI2YvgiIiIiaiOGLyIiIqI2YvgiIiIiaiOGLyIiIqI2YvgiIiIiaiOGLyIiIqI2yjy9UJWIyCoAT5S9HzAnDV9d9k40qFP3nfvdXlXZ79eo6sSyd6IIFRq/gOr8/+bF/W6vTt1voDr77h3DOip8VYWILOjU88116r5zv9urU/ebwnTq/y/3u706db+B6u87px2JiIiI2ojhi4iIiKiNGL4ac1HZO9CETt137nd7dep+U5hO/f/lfrdXp+43UPF9Z88XERERURux8kVERETURgxfTRKRU0RERWT7svclhIh8X0QeFJF7ReQ3IjKh7H1KIyIzRWSpiCwTkdPK3p9QIrKTiNwoIg+IyBIR+WLZ+5SHiPSIyEIR+V3Z+0KtxTGstTpxDOP41XoMX00QkZ0AHAbgybL3JYc/AthTVfcG8BCAr5W8P4lEpAfABQAOBzAVwEdEZGq5exWsD8DJqroHgDcDOLGD9h0AvgjggbJ3glqLY1hrdfAYxvGrxRi+mnMugK8A6JjGOVX9g6r21W7+HcCkMvcnw3QAy1T1UVXdDGAOgCNK3qcgqrpSVe+uXX8JZiDYsdy9CiMikwC8B8BPy94XajmOYa3VkWMYx6/WY/hqkIi8H8DTqrqo7H1pwicBXF/2TqTYEcBTzu3l6JABwCUikwHsA+D2cvck2I9g/iAPlL0j1Docw9qi48cwjl+tMbzsHagyEfkTgFd5HjodwNcBvLO9exQmbb9V9braNqfDlJYvb+e+5SSe+zrmEzoAiMg4AL8C8CVVfbHs/ckiIu8F8Kyq3iUiby97f6g5HMNK19FjGMev1mH4SqGqh/ruF5G9AEwBsEhEAFP2vltEpqvqP9q4i15J+22JyHEA3gvgEK32WiPLAezk3J4EYEVJ+5KbiIyAGbguV9Vfl70/gQ4A8H4ReTeAUQC2FpFfqOrRJe8XNYBjWOk6dgzj+NVaXOerACLyOIBpqlqFk3imEpGZAH4I4CBVXVX2/qQRkeEwDbWHAHgawJ0APqqqS0rdsQBi/qJdCmCNqn6p7P1pRO2T4ymq+t6y94Vai2NYa3TqGMbxq/XY8zX0nA9gKwB/FJF7ROTCsncoSa2p9nMA5sM0fF5d9UHLcQCAYwAcXPt3vqf2aYyImsMxrPU4frUYK19EREREbcTKFxEREVEbMXwRERERtRHDFxEREVEbMXwRERERtRHDFxEREVEbMXwRERERtRHDFxEREVEbMXwRERERtdH/B3ZjheUKyI0WAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"# To use the CDT first install the PyTransKit package using:\n",
"# pip install pytranskit\n",
"from pytranskit.optrans.continuous.cdt import CDT\n",
"\n",
"N = 400\n",
"dt = 0.025\n",
"t = np.linspace(-N/2*dt, (N/2-1)*dt, N) \n",
"f = 1\n",
"epsilon = 1e-8\n",
"\n",
"# Original signal\n",
"gwin = np.exp(-t**2/2)\n",
"z = gwin*np.sin(2*np.pi*f*t) + epsilon\n",
"s = z**2/np.linalg.norm(z)**2\n",
"\n",
"# zero mean additive Gaussian noise\n",
"sigma = 0.1 # standard deviation\n",
"SNRdb = 10*np.log10(np.mean(z**2)/sigma**2)\n",
"print('SNR: {} dB'.format(SNRdb))\n",
"noise = np.random.normal(0,sigma,N)\n",
"\n",
"# Signal after delay and dispersion\n",
"omega = 0.8\n",
"tau = 10.3*dt\n",
"gwin = np.exp(-(omega*t - tau)**2/2)\n",
"zg = gwin*np.sin(2*np.pi*f*(omega*t - tau)) + noise + epsilon\n",
"r = zg**2/np.linalg.norm(zg)**2\n",
"\n",
"# Plot s and sg\n",
"fontSize=14\n",
"fig, ax = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(10,5))\n",
"ax[0].plot(t, s, 'b-',linewidth=2)\n",
"ax[0].set_title('$s(t)$',fontsize=fontSize)\n",
"\n",
"ax[1].plot(t, r, 'r-',linewidth=2)\n",
"ax[1].set_title('$r(t)$',fontsize=fontSize)\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Noise correction\n",
"Expression of noise-corrected CDF is given by,\n",
"\\begin{align}\n",
" \\tilde{S}_g(t) = \\frac{E[R(t)]\\{\\mathcal{E}_z+\\sigma^2(t_N-t_1)\\}-\\sigma^2(t-t_1)}{\\mathcal{E}_z}, \\quad t_1\\leq t \\leq t_N\n",
"\\end{align}\n",
"where $R(t)$ is the CDF associated with $r(t)$, $\\sigma$ is the standard deviation of noise, and $\\mathcal{E}_z$ is the total energy of the noise free signal."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAAFDCAYAAAAef4vuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxV1b3+8c83A4PIGOYACcoYQUEjIoIiOAEKWisOVNt7sdirtnoraqvVn6UDrdNttV5bqEOLFy1V69CiqMWKCg4og8zzEGaRSQIhw/r9sU8gpIGck5xknX3yvF+v89pnWDl5aOPOk7XX2ducc4iIiIhI1aT4DiAiIiISZipTIiIiItWgMiUiIiJSDSpTIiIiItWgMiUiIiJSDSpTIiIiItWQ5juAiIiIL2aWBtwJtAXqA4uAmcBS55wzsyZAZ+fcAo8xJcGZzjMlIiJ1lZmdAaxwzu0zsxSgD3AR0B34EtgIPO70y1KOQ2VKREREpBq0ZkpERKQSZtbczLaZ2cmVjHvRzH5YW7kkMahMSa3RzkhEEpWZZZrZE2a20swOmtl2M5tpZmdFhtwDTHfOrS73dQ+b2Ztlnvop8BMza1pb2cU/lSmJi8hOx0VuhZEd0nfLDdPOSEQSjpllAfOATOA7QA/gCmAuUGhmJwA3Ak9V8OVnAp+UPnDOfQGsAb5Vs6klkahMSbz0JShL7YAuwAvAH8ysL4B2RiKSwH4AFANXOuc+dM6ti2zvcs59DgwHSoAPS7/AzNLN7BBwLnBf5A/JxZGXXwOureV/g3ikMiXVFjls1wx40zm31Tm3HvgDYMApkWHaGYlIomoO1AOyj/H6IOCzcp/oKwbOjtw/i+APyYGRx58A/cysYfyjSiJSmZJ4OAPYCywEMLO2wEME5enzyBjtjEQkUT1GsA9baWafR5Ye9CnzehawpewXOOdKCPZZ+4BPI39I7oq8vBlIB9rXfHRJBCpTEg9nACcCe8wsn2Cn8w3gDufcksgY7YxEJCE55+YTLE84D/g7wUz6Z2b27ciQhsDBCr60L7CggnNQHSjzdVIHqExJPJwB/JHgZHcDgRnAZOfcb8qM0c5IRBKWc67YOfe+c+5+oDdHr9v8kuBQYHl9CBaul9cist0R96CSkFSmJB76ArOdc6siizVvAm42s95lxmhnJCJhYUADjux/5gE5FYw7jcjyhnJ6AZudc9tqJp4kGpUpqRYz60xQfL4ofS6yAH0ecH2ZodoZiUjCMbPnzOxeM+tvZllmdh7wKtAU+FVk2Aygp5lllPvyNKCHmbU3s2Zlnh8EvInUGSpTUl1nECw0X1ru+bcJztNSSjsjEUlEnxGskXodWE6wZGET0Mc5txAOn67lE+Cacl97b+S5PGAigJk1INj3Ta6N8JIYdG0+qRYzm0hwbpZu5Z4fCrwD9HLOLY48Nwd4zjn3RJlxY4BfEyw0/4Nz7r8iO6NtwMXOuY9q6Z8iInJMZnYJ8FsgxzlXfJxxtwCjnHMX1Vo48U5lSmqNdkYiEmZm9gPg1chShmONGQe855xbXnvJxDeVKalV2hmJiEiyUZkSERERqQYtQBcRERGpBpUpERERkWpI8/WNW7Zs6bKzs319exHx4LPPPvvSOdfKd47q0v5LpO453v7LW5nKzs5m7ty5vr69iHhgZsf84EGYaP8lUvccb/+lw3wiIiIi1aAyJSIiIlINKlMiIiIi1aAyJSIiIlINKlMiIiIi1aAyJSIiIlINKlMiIiIi1VBpmTKzp81su5ktOsbrZmaPmdkqM1toZqfHP6aIiIhIYopmZupZ4JLjvD4M6Bq5jQOerH4sERERkXCo9AzozrlZZpZ9nCGjgD875xzwkZk1M7N2zrktccooIj6VlEB+PuzbB3v3wo4dMGcOLF8OX38NDz0EHTv6TikiNcA5R2FhIfXq1QNg/fr1TJ06laKiIgoLCyksLCQ/P5+dO3finOOhhx6iffv2ADzxxBPMnj27wvft0aMH9913HwAHDx5k7Nixx8xw8803c8455wDwj3/8g6lTp1Y4rn79+jz99NOHH48fP54tWyquIsOHD2fMmDGV/OujF4/LyWQCG8s8zos892//AjMbRzB7RadOneLwrUWkSvbvD4rRokXw0UewahVs3gy7dsHu3cFrBw9CQQEcOnT89/rxj+tEmdL+S5LZl19+ydSpU9m1axf79+9n48aNfPzxx+zbt4///M//5Ne//jUA69at45577jnm+9x///2H73/44Yc8//zzFY4799xzD5epoqKiYxYkgBEjRhwuU8uWLTvm2EaNGh1Vpl5//XVWrFhR4di2bdsmXJmyCp5zFQ10zk0CJgHk5uZWOEZEquDQoWCmKD8fDhw4ctu3LyhImzfD4sVBadq7F7Zuje39TzgBGjc+cjvtNOjbF1q1gjpSLLT/krBbuXIlP/3pT9m1axclJSXce++9DBw4EIDvfe97vPTSSxV+3bx58w7fz8rK4u677yYtLY309HTS09Np0KABGRkZpKWl0a5du8Njb775ZkaMGFHhe7Zu3frw/fr16/Pcc88dM/eAAQMO3x8+fDht27atcFxa2tGV5qGHHmLfvn0Vju3Ro8cxv19VxKNM5QFl/yztAGyOw/uK1B1r1wYFZ8+eI7fdu4Pt118HM0n79wfFqPS50lt+fjCLVFIS/ferXx+aNoXmzWHo0KAcdeoELVpAs2bQpAk0bAgNGkC9emAV/c0kImFQWFjIHXfcweOPP37U82UPrV166aUsW7aMK664gkaNGtG2bVtOPfVU2rRpc1RBys7O5le/+lVU33fgwIGHy9rxpKenRz1L1LNnT3r27BnV2JEjR0Y1Lh7iUaZeA241sxeAs4A9Wi8lUgHnYOZMWLEimCFauRIKC2HLFliwoPrvf9JJkJERlKDSW+PGQTlq3RpycqB7dzjxROjQAVJTq/89RSShFRcXc9NNN/HMM88AcNppp/GTn/yERo0a0adPn8PjxowZww033EBKis6YVBWVlikzex4YDLQ0szzg/wHpAM653wPTgeHAKiAf+I+aCisSKsuWwaefwqxZ8MEHsG1bMLNUkbS04LBZs2bBrWnTI7fGjaFRo+BQW7NmwWzSiSceuZ1wQjB7FFkgKiJSasqUKTzzzDM0aNCA119/naFDh2IVzDSnp6d7SJc8ovk037WVvO6AW+KWSCSs/vlP+NvfghmnJUsgL+/fxzRtCldcESzY7tMnOIyWlganngrHWAcgIlJVY8aMoX79+jRu3JgLLrjAd5ykFY/DfCJ1h3PBYblNm4JF3aW3FSvgr38NXi/VtCkMGQKdO8N11wVlqWXLYL2SiEgtSE9P59prjzsnInGgMiUSrY0b4frr4b33Kn7dDG65BS6+GLp2hW7dQOsPRMSDOXPm0Lx587h/ak0qpjIlEo1p0+CGG4LzLkGwvql9++CWmRlszzknWOQtIuLZ3XffzQcffMAbb7zBxRdf7DtO0lOZEqnM228fKVJDhsCjjwanEhARSUCLFi3igw8+ID09nbPPPtt3nDpBxyBEjqewEMaODYrU978fLDJXkRKRBPXee+/Rv39/nHNcccUVNGnSxHekOkFlSuR4Hn44WCvVowf8z//4TiMickxbtmzh6quvZv/+/YwePZpJkyb5jlRn6DCfyLFMmwal16B6+GGd5FJEEpZzjssuu4xt27Zx/vnnM3XqVFK1z6o1mpkSqUhx8ZEi9ctfwjGuLyUiUtt2797NxIkTuemmmyiJXEbKzA5fBkZFqvZpZkqkIu+/D6tXQ1YW3Hmn7zQiIkBweZjzzz+f+fPnA3Dddddx3nnnAfDzn/+c7OzsY14IWGqOypRIeXl5MH58cP+664IzlIuIJIBnn32W+fPn07p1a372s5/RuXPnw68NGjTIY7K6Tb8lRMrasyc4/cHKlcElX26+2XciEREAtm/fzi23BFdvmzBhAuPGjfOcSEppzZRIKeeC0yCsXBlcK+/zz6FDB9+pREQAmDFjBgUFBQwaNEhFKsFoZkpk+3b4wx/gL3+BxYuhSRN48cXgOnoiIgnizTffBODKK6/EzDynkbJUpqRuWrYMnn0W3nwTFiw48nyjRvDcc8G19UREEsioUaMoKipi+PDhvqNIOSpTUvcsXgz9+kF+fvA4JQWGDw/WRw0aBCee6DefiEgFRo8ezejRo33HkAqoTEndUlwMd98dFKmhQ+H22+GCC6BBA9/JREQkpLQAXeqWp56Cf/wDTjghOMx36aUqUiKS0Jxz3HPPPbz66qsUFxf7jiMV0MyU1C3TpgXbRx7RJ/VEJBRWrFjBxIkTadWqFdu2bfMdRyqgmSmpO/bsgffeC9ZIXXWV7zQiIlGZOXMmAEOGDNGn+BKUypTUHW+9BUVFcM45kJHhO42ISFTKlilJTCpTUne89FKwvewyvzlERKLknOPdd98FYOjQoZ7TyLGoTEndsHs3vPIKmMHVV/tOIyISlU2bNrFz504yMjI46aSTfMeRY1CZkrrhrbegoADOOw86dfKdRkQkKosXLwbglFNO0XqpBKYyJXXDjBnBdtgwvzlERGJQUlJCnz59OOOMM3xHkePQqREk+TkXzEwBXHSR3ywiIjEYNmwYw/RHYMLTzJQkv2XLIC8P2rSBU0/1nUZERJKMypQkv/ffD7ZDhgTnmBIRCQHnHOvXr8c55zuKVEK/WST5LVwYbE8/3W8OEZEYbNy4kezsbLp06eI7ilRCZUqSX2mZOu00vzlERGJQ+km+TvoEcsJTmZLk5tyRMqX1UiISImVPiyCJTWVKktuGDcE1+Vq3Dhagi4iEhMpUeKhMSXLTrJSIhNSqVasA6N69u+ckUhmVKUluWi8lIiG1ceNGQGumwkBlSpLbggXBVjNTIhIiJSUlbNq0CYAOHTp4TiOV0RnQJbnpMJ+IhJBzjpkzZ7J161YaNGjgO45UQmVKkteBA7ByJaSlQc+evtOIiEQtNTWVQYMG+Y4hUdJhPkle69ZBSQlkZ0P9+r7TiIhIklKZkuS1bl2wzcryGkNEJFZTpkxh/PjxfP75576jSBRUpiR5rV8fbFWmRCRkXn75ZR555BGWL1/uO4pEQWVKkpfKlIiEVOkJO3v16uU5iURDZUqS19q1wVZlSkRC5MCBA6xatYq0tDSdsDMkoipTZnaJmS03s1Vm9qMKXu9kZu+a2TwzW2hmw+MfVSRGX3wRbHUpBhEJkaVLl+Kco1u3btSrV893HIlCpWXKzFKBJ4BhQA5wrZnllBv2E2Cac64vcA3wv/EOKhKTAwdg+XJISVGZEpFQWbRoEaBr8oVJNDNT/YBVzrk1zrlDwAvAqHJjHNAkcr8psDl+EUWqYPFiKC6G7t2hYUPfaUREolb6Cb5TdbLh0IimTGUCG8s8zos8V9YDwLfMLA+YDny/ojcys3FmNtfM5u7YsaMKcUWiVHoZGV2TT+JE+y+pLR06dKBv376cffbZvqNIlKIpU1bBc67c42uBZ51zHYDhwBQz+7f3ds5Ncs7lOudyW7VqFXtakWipTEmcaf8ltaX0/FJDhw71HUWiFE2ZygM6lnncgX8/jDcWmAbgnJsDNABaxiOgSJXMnx9s+/Txm0NERJJeNGXqU6CrmXU2s3oEC8xfKzdmAzAUwMx6EpQpzYOLH0VF8Nlnwf3TT/ebRUQkBu+88w4rVqzAufIHgCSRVVqmnHNFwK3ADGApwaf2FpvZBDMbGRl2B/BdM1sAPA98x+knQXxZuBDy86FLF2jd2ncaEZGoOOcYM2YM3bt3Z+XKlb7jSAzSohnknJtOsLC87HP3l7m/BDgnvtFEqmjGjGA7YIDfHCIiUSouLuaVV15h+/bttGvXjq5du/qOJDHQGdAluTgHf/pTcH/0aL9ZRESOY9asWfTr148WLVqQnp7ON7/5TQBuvfVWzCr67JckqqhmpkRCY/Xq4GSdzZvDxRf7TiMickxpaWnk5+eza9cuzIyMjAyuv/567rrrLt/RJEYqU5Jc3ngj2F5wAaTpx1tEEteAAQP4+OOP2b9/PxkZGaSmpvqOJFWk3zYSfmvXBueVmjcPHnsseG5U+ZP0i4gknkaNGtGoUSPfMaSaVKYkvA4ehKuugr///ejnR4yAa6/1k0lEpBLvvfceU6dO5Vvf+haDBg3yHUfiQGVKwuuXvwyK1AknwMCB0Ls39OsHV14ZXOBYRCQBvfrqq0yaNInWrVurTCUJlSkJpx074JFHgvtvvgnaIYlISMyPXKGhf//+npNIvOjPdwmnqVODE3NecomKlIiEhnOOBZFrh56ma4cmDZUpCaeXXw62N9zgN4eISAxWr17NV199RYsWLcjMzPQdR+JEZUrC59Ah+OgjMINhw3ynERGJ2osvvgjAsGHDdGLOJKIyJeGzaFFQqLp1g2bNfKcREYnav/71LwCuuOIKv0EkrrQAXcJn7txge/rpfnOIiMSoX79+OOfo1auX7ygSRypTEj6vvBJsBw70m0NEJEYTJkzwHUFqgA7zSbh8/TW89RakpgYn7BQREfFMZUrCZf16KC6Gk0+GVq18pxERidq+fftYvHgxe/bs8R1F4kxlSsJl/fpgm5XlN4eISIw++ugjevXqxeWXX+47isSZypSES2mZ6tTJbw4RkRjt2LEDgDZt2nhOIvGmMiXhsmFDsNXMlIiEzJdffglAy5YtPSeReFOZknDZuDHYduzoN4eISIxKZ6ZUppKPypSEy7ZtwbZtW785RERiVDoz1Uofnkk6KlMSLtu3B9vWrf3mEBGJkQ7zJS+VKQkXlSkRCSmVqeSlM6BLeJSUQGTNgc4xJSJh8+STT5KXl8fpuhRW0lGZkvDYvTs4YWfTplC/vu80IiIx6dGjBz169PAdQ2qADvNJeOgQn4iIJCCVKQmPrVuDrcqUiITMjh07GDt2LA8++KDvKFIDVKYkPDZtCraZmX5ziIjEaN26dTz99NM8//zzvqNIDVCZkvBQmRKRkNoamVlv166d5yRSE1SmJDxUpkQkpLZs2QKoTCUrlSkJD5UpEQmpNWvWAJCl64omJZUpCQ+VKREJqRUrVgDQvXt3z0mkJqhMSTg4B8uXB/dPPtlvFhGRGJWWqW7dunlOIjVBJ+2UcNi+HXbtgiZNQGsORCRkcnJyKCwspEuXLr6jSA1QmZJwWLo02ObkgJnfLCIiMZo2bZrvCFKDdJhPwmH16mCrKXIREUkwKlMSDqVnP2/f3m8OEZEYHThwgN27d+Oc8x1FaojKlIRDaZlq08ZvDhGRGE2fPp3mzZvzzW9+03cUqSEqUxIO27YFW5UpEQmZ7ZGLtGdkZHhOIjVFZUrCobRMtW3rN4eISIx27NgBQGtdpD1pqUxJOOgwn4iEVGmZatWqleckUlOiKlNmdomZLTezVWb2o2OMGW1mS8xssZlNjW9MqfN0mE9EQqr0MJ9mppJXpeeZMrNU4AngQiAP+NTMXnPOLSkzpivwY+Ac59wuM9NPjMTPwYOwZw+kp0Pz5r7TiIjERDNTyS+amal+wCrn3Brn3CHgBWBUuTHfBZ5wzu0CcM5tj29MqdMif9XRujWk6Mi0iISLZqaSXzRnQM8ENpZ5nAecVW5MNwAz+xBIBR5wzr1Z/o3MbBwwDqBTp05VySt1kdZLSQLQ/kuq6vHHHycvL4/OnTv7jiI1JJoyVdG1O8qfeSwN6AoMBjoA75tZL+fc7qO+yLlJwCSA3Nxcnb1MoqNP8kkC0P5Lqur888/3HUFqWDTHTPKAjmUedwA2VzDmVedcoXNuLbCcoFyJVJ8Wn4uISAKLpkx9CnQ1s85mVg+4Bnit3JhXgPMBzKwlwWG/NfEMKnWYDvOJSEitX7+e++67j6lT9SH3ZFZpmXLOFQG3AjOApcA059xiM5tgZiMjw2YAO81sCfAucKdzbmdNhZY6ZunSYKv1BiISMrNmzeLnP/85L730ku8oUoOiWTOFc246ML3cc/eXue+AH0ZuIvH16afB9swz/eYQEYnRvHnzAOjbt6/nJFKT9DlzSWx79sDKlVC/PvTq5TuNiEhMFi5cCECfPn08J5GapDIliW3ZsmDbs2dw0k4RkRBZuXIlAD169PCcRGqSypQkthUrgm23bn5ziIjEqKCggI0bN5KamkpWVpbvOFKDVKYksalMiUhIrV27FuccWVlZpGtmPalFtQBdxJvly4NtV522TETCpaCggNzcXDIzM31HkRqmMiWJbf78YNu7t98cIiIxOu200/i09NPIktR0mE8S1759wSf50tPhlFN8pxEREamQypQkrtJZqV69oF49v1lERGK0ceNGSkpKfMeQWqAyJYlr9uxge9ZZfnOIiMTo66+/Jjs7m4yMDA4dOuQ7jtQwlSlJXB98EGwHDvSbQ0QkRnPnzqWkpISTTjqJeppZT3oqU5K4Ipdh0MyUiIRN6cLzs7T/qhNUpiQx5efDpk3B4vPsbN9pRERismrVKgBycnI8J5HaoDIliWn16mDbuTOk6QweIhIua9euBeCkk07ynERqg8qUJKbIX3U6WaeIhNGaNWsA6Ny5s+ckUhtUpiQxvfhisO3Vy28OEZEYlZSUsHXrVgCytUyhTtDxE0k8S5fC888H66Vuvtl3GhGRmKSkpLBs2TKmTJlCw4YNfceRWqCZKUkshYVw003gHIwdC506+U4kIhKzDh068OMf/9h3DKklKlOSGJyDxx+H00+H99+Hdu3ggQd8pxIRiZpzjtdff52CggLfUaSWqUxJYrjrLvjBD2DRImjSBP76V2jTxncqEZGo/P73v6dFixaMHDmSwYMHU1RU5DuS1CKtmRL//vIXePjhYI3UpEnwjW8EhUpEJAR++9vfcvvttwPQuHFjbrjhBtJ0Spc6Rf9vi19ffQXf+15w/7HH4Dvf8RpHRCQW8+bNY/z48QBMnjyZG2+80XMi8UGH+cSvGTNg924YMCBYeC4iEiK33XYbRUVF3HLLLSpSdZjKlPg1c2awHTUKzPxmERGJwb59+/jwww9JS0vjl7/8pe844pEO84k/u3cHC80BLrzQbxYRkRg1atSIBQsWsHTpUpponWedpjIl/vz977BnDwwaBH37+k4jIhKTlJQUevXqRS9dqaHO02E+8af0+nsDB/rNISIiUg0qU+JP5Krq6KrqIhIyGzZs4JxzzuG2227zHUUSgA7ziT+lZUpXVReRkFmzZg2zZ8/G9MEZQTNT4tOaNcFWZUpEQmb9+vUAZGVleU4iiUBlSvwoKIDNmyElBTp29J1GRCQmGzZsAKCTLsYuqEyJL+vXBxc37tQpuIyMiEiIaGZKylKZEj+0XkpEQkxlSspSmRI/tF5KREJMZUrK0qf5xA/NTIlIiF1++eWsWLFCa6YEUJkSX3SOKREJsQcffNB3BEkgOswnfmhmSkREkoTKlPihNVMiElJ5eXl8/PHH7Ny503cUSRAqU1L79uyBXbugYUNo08Z3GhGRmLzwwgv079+fCRMm+I4iCUJlSmpf2UN8uhSDiISMPskn5alMSe3TeikRCTGd/VzKU5mS2qf1UiISYpqZkvKiKlNmdomZLTezVWb2o+OM+6aZOTPLjV9ESTo6LYKIhJjKlJRXaZkys1TgCWAYkANca2Y5FYxrDPwA+DjeISXJbNwYbLUjEpGQ2bt3L7t376ZBgwa0atXKdxxJENHMTPUDVjnn1jjnDgEvAKMqGPcz4EHgYBzzSTLaujXYtmvnN4eISIzKrpcyfYBGIqI5A3omsLHM4zzgrLIDzKwv0NE593czG3+sNzKzccA40MK9Oq20TLVt6zeHSAy0/xKAnj17sm7dOvbu3es7iiSQaGamKqre7vCLZinA/wB3VPZGzrlJzrlc51yupkfrKOeOlCmdY0pCRPsvAUhNTSUrK4vevXv7jiIJJJoylQd0LPO4A7C5zOPGQC/gX2a2DugPvKZF6FKh3bvh0CFo0gROOMF3GhERkWqLpkx9CnQ1s85mVg+4Bnit9EXn3B7nXEvnXLZzLhv4CBjpnJtbI4kl3HSIT0RC7Mknn+Tqq6/mnXfe8R1FEkilZco5VwTcCswAlgLTnHOLzWyCmY2s6YCSZDZHJjVVpkQkhGbPns20adPYtGmT7yiSQKJZgI5zbjowvdxz9x9j7ODqx5KkVXpahI4djz9ORCQBbd++HYDWrVt7TiKJRGdAl9qlMiUiIaYyJRVRmZLapTIlIiG2NbLuU2VKylKZktqlMiUiIVVUVMS2bdswM9pq3aeUoTIltWvdumCrkx6KSMhs3boV5xxt2rQhPT3ddxxJIFEtQBeJi4ICWLkSUlKgWzffaUREYlJcXMwVV1xBkyZNfEeRBKMyJbVn6VIoLg6KVMOGvtOIiMQkKyuLl19+2XcMSUA6zCe1Z9GiYKvLMIiISBJRmZLas3JlsO3e3W8OEZEqWLFiBZs3b8Y5V/lgqVNUpqT2rF4dbE8+2W8OEZEquPnmm8nMzOSNN97wHUUSjMqU1B6VKREJsSVLlgDQs2dPz0kk0ahMSe1RmRKRkNq7dy9btmyhQYMGZGVl+Y4jCUZlSmpHQQHs2AGpqdC+ve80IiIxWRc5R17nzp1JSdGvTjmafiKkdmzbFmzbtAnOMyUiEiKlZSo7O9trDklM+q0mtaNsmRIRCRmVKTkelSmpHZGLg6LrWYlICJU9zCdSns6ALrWjdGZKZUpEQujOO+/kqquuoqMu0i4VUJmS2lE6M6XDfCISQu3ataNdu3a+Y0iC0mE+qR2bNwdb7YxERCTJqExJ7Vi7Nthq8aaIhNC4ceP4wQ9+QH5+vu8okoBUpqR2qEyJSEgVFxfz1FNP8fjjj5Oenu47jiQglSmpec7B+vXBfX0SRkRC5quvvqKkpITmzZurTEmFVKak5m3aBAcPQosW0Lix7zQiIjHZvn07AK1atfKcRBKVypTUvJkzg21urt8cIiJVkJeXB0CHDh08J5FEpTIlNW/69GA7YoTfHCIiVbBhwwYAOnXq5DmJJCqVKalZJSVHZqYuvthvFhGRKlCZksqoTEnNWrQIduyAzEzo1s13GhGRmLVt25b+/fuTk5PjO4okKJ0BXWpW6azUkCFg5jeLiEgV3HLLLdxyyy2+Y0gC08yU1KzSMjV0qN8cIiIiNXLKRF4AABHnSURBVERlSmpOURG8915w//zz/WYREamCdevWsXTpUpxzvqNIAlOZkprz2Wewdy906QJauCkiIfS///u/5OTkMGHCBN9RJIGpTEnN2LABxo0L7l94od8sIiJVNGvWLAD69+/vOYkkMpUpiZ9Vq2DMGGjbFrKyYOFC6N4dfvIT38lERGL29ddfM3fuXFJSUhgwYIDvOJLA9Gk+iY9XXoHRo6GwMHicng4jR8KTT4IuwSAiITRnzhyKi4vJzc2lsS6FJcehmSmpPufgrruCIjVyJHzyCeTnw4svqkiJSGiVHuI777zzPCeRRKeZKam+d9+FlSuDE3O+9BKk6cdKRMKvtEyde+65npNIotPMlFTfpEnB9rvfVZESkaTgnKNt27YADBw40HMaSXQqU1I9BQXwj38E97/zHa9RRESqqqSkhOnTpzN69Gh2796NmTFlyhQeeughWrRo4TueJDiVKame996Dr7+G3r2DT/CJiITMvn37GDRoECNGjOCvf/0rhw4dAqBevXqMHz/eczoJAx2Tkep54olg+41v+M0hIlJFTz31FLNnz6Zp06Z873vfo2nTpr4jScioTEnVHTwYHOJLSYGbb/adRkQkZs45nnrqKSAoVVdeeaXnRBJGUR3mM7NLzGy5ma0ysx9V8PoPzWyJmS00s3+amY731AVffAHFxdCjB7Ru7TuNiEjMPvnkExYtWkSrVq247LLLfMeRkKq0TJlZKvAEMAzIAa41s5xyw+YBuc65U4EXgQfjHVQS0Lx5wbZvX785RESq6OWXXwZgzJgx1KtXz3MaCatoZqb6Aaucc2ucc4eAF4BRZQc45951zuVHHn4EdIhvTElI8+cH2z59/OYQEamiVq1a0a1bN0aMGOE7ioRYNGUqE9hY5nFe5LljGQu8UdELZjbOzOaa2dwdO3ZEn1ISk2ampA7R/is5jR8/nuXLlzN06FDfUSTEoilTVsFzrsKBZt8CcoGHKnrdOTfJOZfrnMttpcuMhFtxcXAhY1CZkjpB+6/kZlbRrzqR6ERTpvKAjmUedwA2lx9kZhcA9wIjnXMF8YknCWvFiuD6e506gU5oJyIhNHPmTJYuXUpxcbHvKBJy0ZSpT4GuZtbZzOoB1wCvlR1gZn2BPxAUqe3xjykJR+ulRCTEnHNcddVV5OTksGXLFt9xJOQqLVPOuSLgVmAGsBSY5pxbbGYTzGxkZNhDwInAX81svpm9doy3k2Sh9VIiEmJbtmzhq6++onnz5mRmHm8ZsEjlojppp3NuOjC93HP3l7l/QZxzSaJTmRKREFsYWfPZu3dvrZeSatO1+SR2zqlMiUioffHFFwCceuqpnpNIMlCZktjl5cHOndC8OXTsWPl4EZEEU3ZmSqS6VKYkdh99FGzPOAM0PS4iIaSZKYknlSmJ3RuRc7JeoKVyIhI+xcXFrFq1CoBTTjnFcxpJBlEtQBc5zLkjZWrYML9ZRESqIDU1lZ07d7J69WoaN27sO44kAc1MSWzmz4etWyEzE7TWQERCqn79+uTk5PiOIUlCZUpi8847wfaSS7ReSkREBJUpidXs2cH23HP95hARqaIbbriBwYMHM6/0FC8i1aQ1UxI9546UqQED/GYREakC5xwzZ85k06ZNNGzY0HccSRKamZLorVkD27dD69Zw8sm+04iIxGzVqlVs2rSJjIwMunXr5juOJAmVKYnehx8G2wEDtF5KRELp7bffBmDo0KGkpOhXoMSHfpIkejNnBlsd4hORkCotUxdeeKHnJJJMVKYkOgUF8Morwf3LLvObRUSkCoqKinj33XcBlSmJL5Upic7778OePcG5pXr08J1GRCRmS5YsYc+ePWRnZ5OVleU7jiQRfZpPolM6K3XxxX5ziIhUUfv27Zk8eTIlJSW+o0iSUZmSyr39NjzxRHB/+HC/WUREqqhly5bceOONvmNIElKZkootWwavvw4rV8L//V/w3B13wODBXmOJiIgkGpUpOdrs2XDvvTBrFpSdCr/6apg4UadEEJFQ2r9/P/fffz8DBgzgyiuv9B1HkozKlAScg1/8Ah54AIqLIT0drrkGzjwTBg6E3FzfCUVEquzDDz/k0Ucf5d1331WZkrhTmZLAH/8I990X3L/zTrjnHmjWzG8mEZE4efnllwG49NJLPSeRZKQyJXDo0JEi9ac/wQ03+M0jIhJHzjmmT58OwKhRozynkWSk80wJTJsG27bBKafA9df7TiMiElcrVqxg48aNtGjRgr59+/qOI0lIZaqucw4mTAju//CHWmAuIkmlqKiIb3/72wAMHjxY1+OTGqGfqrruiy+C0x+0bavDeyKSFPLz8w/fX716NR9//DFt2rThF7/4hcdUksxUpuq6yDoChg2DNC2hE5HwmjNnDqeddhpdu3Y9XKi6d+/OmDFjeO655+ihS2FJDVGZqsucg6lTg/sjR/rNIiJSDRMnTmTgwIEsXLiQvXv38tZbbx1+bcqUKVxwwQUe00myU5mqy15/PTjM16KFLhMjIqE1Z84c7rnnHgB+9KMfsWPHDi6//PLDr5vWgkoN03GdumrvXviv/wru33cf1KvnN4+ISBU98sgjANx5551MnDjRcxqpizQzVVc9+yxs3gz9+sH3v+87jYhIlaxfv55XX32VlJQUvq99mXiimam6qLgY/vCH4P748ZCa6jePiEgVffDBB7Rs2ZKhQ4eSmZnpO47UUSpTddHPfw5LlkCnTlp4LiKhNmbMGK677joKCgp8R5E6TGWqLikpgYcfDi5mnJICv/891K/vO5WISLWYGQ0aNPAdQ+owlalktncvzJ0L69YFM1GffALvvx+89rvfBeeWEhEJqYULF5KZmUlGRobvKFLHqUwli5074dFHYd48WLMG9uyBL7+EoqKjx2VkwNNP6/CeiITa/v37GTlyJJs2bWLWrFmcffbZviNJHaYyFXZr1gSzTH/+c1CoyjvjDOjePbiIcdu2MGIEtGlT+zlFROLogQceYP369fTp04czzzzTdxyp41SmwurgQfjtb4P1TwcPBs/16QP33AM5OdC8OTRrBiec4DWmiEi87d69m8cffxwzY/LkyaTpUljimX4CE11+frDu6aOP4Ouvg8ebNsE77wSH8QCuvRZuvx3OPBN0pl8RSXIvvPACBQUFDBkyhNzcXN9xRFSmvHEOCgpg2bKgFH31VXDbtSvYbtsGq1cHl3vZt6/i9zj9dJg4ES66qHazi4h4snnzZh544AEAvvvd7/oNIxKhMhUPRUWwYQMsXx6Unz17gvVLK1cGM0kHDlR8KymJ7v27dYPBg6F9++CwXUYG5OZC796aiRKROmXs2LFs27aNc889l9GjR/uOIwKoTAUzRAcOBIfQ9u8/UnQOHgxK0vbtwezR3r3BmLK3fftgxw7YuBEKC2P/3unpkJ0NHTsGa5xatDiybdkSTjoJunSBDh1UmkQkqS1evJjPPvuM3bt3U1BQwKFDhygsLKSwsJD69etz//33A/Dkk09y22238dRTT5GSoiuiSWKIqkyZ2SXAb4FU4I/OuV+Ve70+8GfgDGAncLVzbl18o8ZBfn5QjP72N1ixIrgtWBAUqurKzISuXeHUU6FVK2jcOJhRatIkmE1q2PDIrfSxLuMiIoJzjttuu41//vOfFb7esmXLw2UqOzubV199tTbjiVSq0jJlZqnAE8CFQB7wqZm95pxbUmbYWGCXc66LmV0D/Bq4uiYCA0H5KSoKboWF8OGHMHt2MKOUnx/MFm3bFqxFOngwWJtUUFDxqQMgKDYnnnh06WnQIJg5atIkOK1AixbBmPK3Vq2CIqVPzYmIVOjQoUMcOHCAwsJCioqKKCwsZN68eaSnpzNs2DDMjBdffJFbb72VZs2a0bBhQ9LT0w/fmjRp4vufIHJc0cxM9QNWOefWAJjZC8AooGyZGgU8ELn/IvA7MzPn4jDlM2sWPPdcUIbWroWlS498ii1WaWnB7FHfvsHZv9u0gQEDoFGjascUEanIoUOHuPvuu4FgBqbsFoJry/Xr1w+Af/3rX7z00ksVjk1PT+c3v/nN4a+bMGECW7ZsOeq9Su8PGTKEq68O/p5duXIlDz744L+NKfXTn/708AWCJ0+ezOzZsysc27VrV+69997D/6axY8ce8z1vuukmBg0adPj+M888Q2EFSyE6derE0qVLOeGEE2jWrBnPPffcv40RCYNoylQmsLHM4zzgrGONcc4VmdkeIAM4qvWY2ThgHAT/EUVl2TKYPPnfn09NDWaO0tKCNUWXXx4szG7QIJgtat062DZsGFx/rn794LxL6enRfV8RkTKqtP8CiouLjypB5Z1xxhmHy9SCBQv43e9+V+G4Bg0aHPU+f/nLX1iyZMkxx5aWqW3btvHHP/7xmN//v//7vw+Xqffff58pU6ZUOO6cc845XKaKi4uPW3wuuuiiw2Vq6NChTJ48mSZNmpCWlkZ6ejppaWk0a9aM//iP/9A5oiQpRPNTXNHK5/IzTtGMwTk3CZgEkJubG92s1aBBwQV569ULSlOPHsFhNS08FJFaVKX9F8GM0qOPPnr4sUU+TFK6LXv27vPOO4/HHnuswrGp5dZY3nfffezateuocaX3e/Xqdfhxly5dmDRp0lFfW3Z8+/btD9+/8cYbOf/88ysc27p166P+TX/+85+P+Z4DBgw4fH/UqFHs2rWLpk2bIpKsrLIjcWZ2NvCAc+7iyOMfAzjnJpYZMyMyZo6ZpQFbgVbHO8yXm5vr5s6dG4d/goiEhZl95pwL/VkWtf8SqXuOt/+KZnrnU6CrmXU2s3rANcBr5ca8Bnw7cv+bwMy4rJcSERERSXCVHuaLrIG6FZhBcGqEp51zi81sAjDXOfca8BQwxcxWAV8RFC4RERGRpBfVyj/n3HRgernn7i9z/yBwVXyjiYiIiCQ+reIWERERqQaVKREREZFqUJkSERERqQaVKREREZFqUJkSERERqQaVKREREZFqUJkSERERqYZKLydTY9/YbAew3ss3/3ctKXdR5pBQ7toV1tyQONmznHOtfIeoLu2/4iKsuSG82ZW7eo65//JWphKJmc0N4/XClLt2hTU3hDu7HF9Y/78Na24Ib3blrjk6zCciIiJSDSpTIiIiItWgMhWY5DtAFSl37Qprbgh3djm+sP5/G9bcEN7syl1DtGZKREREpBo0MyUiIiJSDSpTIiIiItWgMlWOmY03M2dmLX1niYaZPWRmy8xsoZn9zcya+c50PGZ2iZktN7NVZvYj33miYWYdzexdM1tqZovN7DbfmWJhZqlmNs/M/u47i9Qs7b9qlvZftS8s+y+VqTLMrCNwIbDBd5YYvA30cs6dCqwAfuw5zzGZWSrwBDAMyAGuNbMcv6miUgTc4ZzrCfQHbglJ7lK3AUt9h5Capf1XzdL+y5tQ7L9Upo72P8BdQGhW5Tvn3nLOFUUefgR08JmnEv2AVc65Nc65Q8ALwCjPmSrlnNvinPs8cn8fwX/YmX5TRcfMOgAjgD/6ziI1TvuvmqX9Vy0L0/5LZSrCzEYCm5xzC3xnqYb/BN7wHeI4MoGNZR7nEZL/qEuZWTbQF/jYb5Ko/YbgF2yJ7yBSc7T/qhXaf9W+0Oy/0nwHqE1m9g7QtoKX7gXuAS6q3UTROV5u59yrkTH3Ekzn/l9tZouRVfBcaP6KNrMTgZeA251ze33nqYyZXQpsd859ZmaDfeeR6tH+yzvtv2pR2PZfdapMOecuqOh5M+sNdAYWmBkEU82fm1k/59zWWoxYoWPlLmVm3wYuBYa6xD5xWB7QsczjDsBmT1liYmbpBDui/3POvew7T5TOAUaa2XCgAdDEzJ5zzn3Lcy6pAu2/vNP+q3aFav+lk3ZWwMzWAbnOuUS4SvVxmdklwKPAec65Hb7zHI+ZpREsMh0KbAI+Ba5zzi32GqwSFvyG+hPwlXPudt95qiLyl91459ylvrNIzdL+q2Zo/+VPGPZfWjMVfr8DGgNvm9l8M/u970DHElloeiswg2AR5LRE3xFFnANcDwyJ/G88P/LXkohUj/ZfNU/7r1qgmSkRERGRatDMlIiIiEg1qEyJiIiIVIPKlIiIiEg1qEyJiIiIVIPKlIiIiEg1qEyJiIiIVIPKlIiIiEg1/H/lFlCLQ2HB1wAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAl8AAAFDCAYAAAAEQ2tgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeZwUxf3/8Vct7HIuIIfIJRA1KHJF8ET8qaCi0aARFDyiAURF4u3XK5Ko8Y5HEsUkeCPeaERRvPBWEOQIoqiIqIhRllOWY9mlfn/M9NDT2z3dszs7M7u8n48Hj52p7uqpjaT5zKeqP2WstYiIiIhIdhTkegAiIiIiOxIFXyIiIiJZpOBLREREJIsUfImIiIhkkYIvERERkSxS8CUiIiKSRQq+RERERLJIwZeIiEgIY0w9Y8y/jTGzjTFPG2NOM8Y0j7efl+vxSe2i4EtERCTcAOBf1tp9gb8AewOvAW8Ai8M6G2N2Msb8aIzZLeS8Z4wxF2diwJK/jCrci4iI1CxjzG1Aa2vt711tfwV6WGsHu9p6Am8DXa2167I/UskGZb4kL+lboojUFcaYxsBo4H7PoX2Bj9wN1tqFwFLgtOyMTnJBwZfkq6uAl6y1X7kbjTF/NcZMdzVdC/zRGNM8q6MTkR2Oifk/Y8znxphNxpifjDFTInQ9BtgGvB+/TqExpgw4BLjGGGONMYtc508FRmT8F5C8oeBL8o4xphn+3xLB801R3xJFJIsuA34PjAX2BH5DbN1XmAHAx3b7Op8K4MD46/2BdsDBrvM/AvYzxjTKxKAl/yj4kpwyxnSMf+sbboyZYYzZDJyD61ti/LxU3xT1LVFEsmEwsYz8G9bab6y1M621/4zQrzPwg/PGWruNWMD1MzDbWvs/a+0a1/krgEKgfQbHLnlEwZfkWp/4z8uBvxJ7gqgTyd8SIfU3RX1LFJFsmApcaIx53RhztjGmdcR+jYDNnrZfAQus/1Nvm1z9pA5S8CW51pvYTWmYtdZZ49UJ17dECP2mqG+JIlLjrLV3Ad2A6cSmHr8yxuwFYIzpYYx53xizwBhzpTHmLVfXEmAnz+X6APMCPqpl/OfKjA1e8oqCL8m1PsTS+EtcbX7fEiH4m6K+JYpIVlhrl1hr/wr0AwzQyxhTH3gEGGOt7U1sbeoCV7d5QHfPpXoD/w34mB7ACmvtjxkdvOQNBV+Sa72J1bRx8/uWCMHfFPUtUURqlDHmcmPMmcaY7saYXwJ/BsqAt4DfArOstc461M9IDqxeAfYyxrRytdUH9jTGtDfGtPB83ABi2TWpoxR8Sc4YY5oAuwFzPYf8viVC8DdFfUsUkZrWgNja1DnAB8TuRwPj951ewHzXuXvjynzFn8r+CBjuOufq+PvlwE1OozGmIXACMLFGfgvJC6pwLzljjDkQeA9obq3d4GrvSexGtrO1dpWrfRnwDHAHsNFauzbe/hBQYa0dlb3Ri4jEGGMuAjpaay8xxgwgVn5iJ2vtJtc5g4G/Ad2ttRUprnUeMMRae2RNj1tyR5kvyaXewJfuwAsCvyWCzzdFfUsUkTwwCRhgjPmI2BTkx+7AC8BaOx24B+gYcq2twB9qZJSSN5T5krykb4kiUlsYY5pYa0uNMQa4EfjeWnt3rscl+UuZL8lL+pYoIrXIZcaYT4itVy0AJuR4PJLnlPkSERERySJlvkRERESySMGXiIiISBbVz/UA0tG6dWvbpUuXXA9DRLLk448/LrHWtsn1ODJB9y+RHU/QPaxWBV9dunRhzpw5uR6GiGSJMeabXI8hU3T/EtnxBN3DNO0oIiIikkUKvkRERESySMGXiIiISBYp+BIRERHJokjBlzFmsDHmc2PMEmPMFT7HGxhjnowfn2WM6RJv388YMz/+Z4Ex5oSo1xQRERGpi0KDL2NMPWLbvBwNdAdGGGO6e04bBayx1u4O3AncEm//BOhnre0DDAb+ZYypH/GaIiIiInVOlMzXfsASa+1Sa20Z8AQwxHPOEODh+OtngIHGGGOt3WitLY+3NwScvYyiXFNERESkzokSfHUAvnO9Xx5v8z0nHmytA1oBGGP2N8YsAhYC58SPR7mmiIiISJ0TJfgyPm3e3bgDz7HWzrLW7g3sC1xpjGkY8ZqxCxszxhgzxxgzZ+XKlRGGKyKSH3T/EhE/UYKv5UAn1/uOwIqgc4wx9YHmwGr3Cdbaz4BSoEfEazr9/m2t7Wet7demTZ3YZUREdhC6f4mInyjB12xgD2NMV2NMETAcmOo5ZypwRvz1UGCGtdbG+9QHMMZ0BroByyJeU0RERKTOCQ2+4mu0xgGvAJ8BT1lrFxljrjPG/CZ+2v1AK2PMEuBiwCkdcTCwwBgzH3gOGGutLQm6ZiZ/MRERkdpu4sSJnHnmmbz//vu5HopkUKSNta21LwEvedrGu15vBob59JsETIp6TREREdnuggsuYNOmTXTs2JH+/fvnejiSIapwLyIikqc2bdoEwCeffJLjkUgmKfgSERHJc8b4FQmQ2krBl4iISJ5T8FW3KPgSERHJcwq+6hYFXyIiInlOwVfdEulpRxEREckea5M3fWnXrl2ORiI1QZkvERGRPNOoUSMKCgo4/PDD2WmnnTjllFNyPSTJIGW+RERE8oyT+Zo2bRoNGzbM8Wgk05T5EhERyTNlZWUAfPnll5SVlVFeXp7jEUkmKfgSERHJU2PGjKFBgwbcdddduR6KZJCCLxERkTw1c+ZMABYt0vbHdYmCLxERkTynace6RcGXiIhInlOdr7pFwZeIiEieU/BVtyj4EhERyTNjx45Neq/gq25R8CUiIlIDxo8fz5VXXlml9VoHH3wwZ511VuJ9QYH+ua5L9F9T6r7rr4ebb871KERkB2Kt5frrr+fmm29OO/jatm0b69ato1+/fom2E088MdNDlBxShXup28rKYPz42OsrrsjtWERkh7Ft27bE66+++oq3336bs846i8LCwkh9zz33XAoKCrjxxhtZt24dAwYMqMnhSpYp+JK6raJi+2trQesmRCQLKlz3nh49egCwceNGLr300tC+ztZC27ZtY+jQoeyxxx41M0jJGU07St3mDb5ERLLAnflyzJ07N1Jf67pXDRs2jLvvvpv3338/Y2OT3FPwJXWbO/jyuRmKiNQEv+Ar6qJ5d/C1YMEC/vCHP/DMM89kbGySewq+pG5zL3R1B2IiIjXIL/iqV69epL7WJ0u/adOmao9J8ofWfEnd5g6+lPkSkSxxr/l69tlnWbt2LYccckikvn7Bl+p81S0KvqRushZmzYKddtrepuBLRLLEyXy1aNGCE044Ia2+Cr7qPk07St00ZQoceCD077+9TcGXiGRJcXEx8+fP56233kq7b+PGjZk8eXJSm4KvukXBl9RN06bFfq5atb1Na75EJEvq169P79696d27N23btsUYw2OPPRa5v4Ktuk3Bl9RNTZpUblPmS0SyqGfPnrRs2ZKffvoJgDfffDNy36OOOooPP/ww8V7BWN2iNV9SNyn4EpEcWrNmDZ988klSW9RSE6WlpQwaNCjxdGSDBg24++67Mz5GyR0FX1I3+QVfmnYUkSwpLS2t1BY1e1VeXs68efNo2rQpy5Yt8y1bIbWbph2lbmrcuHKbbmAikiUVPl/2ogZfztOOGzZsYNq0aXTt2jWjY5PcU/AldZOmHUUkh6pT4d7dd9y4cfTq1Ys//vGPGRub5J6CL8lvmzdXrZ9f5kvTjiKSJX7BV7qZL+f1woUL+e677zI2Nsm9SMGXMWawMeZzY8wSY8wVPscbGGOejB+fZYzpEm8/whjzsTFmYfzn4a4+b8WvOT/+Z+dM/VJSR8ycCY0awVVXpd/X7xumMl8ikiV+wddO7qLPKfgVWfWbxpTaK3TBvTGmHnAPcASwHJhtjJlqrf3UddooYI21dndjzHDgFuBkoAQ4zlq7whjTA3gF6ODqd6q1dk6Gfhepa264IfbzppvgxhvT6+tz81LwJSLZ4g6WSkpKaNWqVeS+fsGXFt3XLVEyX/sBS6y1S621ZcATwBDPOUOAh+OvnwEGGmOMtXaetXZFvH0R0NAY0yATA5cdgF8AFZXfjUrfHEUkSxo0aMB+++3H0KFD0wq8ABo2bEjHjh2T2hR81S1Rgq8OgHuyeTnJ2aukc6y15cA6wPu37URgnrV2i6vtwfiU4zVGFeQklRdfhBkzop+vzJeI5FDXrl2ZNWsWTz/9NB988AHvvfceW7dujdS3uLiYf//730ltCr7qlijBl19Q5P2XLeU5xpi9iU1Fnu06fqq1ticwIP7ndN8PN2aMMWaOMWbOypUrIwxX6gx3AHXccTBwYNX6OnTzkizT/WvH9sYbb/DMM8/Qv39/BgwYwO233x65b4sWLSguLk68V/BVt0QJvpYDnVzvOwIrgs4xxtQHmgOr4+87As8Bv7PWfuV0sNZ+H//5M/AYsenNSqy1/7bW9rPW9mvTpk2U30nEP/jStKNkme5fOy5rLRdccAHDhg1LtK1evTpS37KyMurXr5+U/Ro8eHDGxyi5E6XC/WxgD2NMV+B7YDhwiuecqcAZwIfAUGCGtdYaY1oA04ArrbXvOyfHA7QW1toSY0whcCzwerV/GxGHMl8ikkNz585l0aJFSW1RV9eUlJSw336xfMRRRx1Fx44dGT16dMbHKLkTGnxZa8uNMeOIPalYD3jAWrvIGHMdMMdaOxW4H5hkjFlCLOM1PN59HLA7cI0x5pp425FAKfBKPPCqRyzwmpjB30t2dAq+RCSHMlHnq3379rz00kuRi7NK7RFpb0dr7UvAS5628a7Xm4FhPv3+Avwl4LJ9ow9Tdkh+AZS1EOUGpqcdRSSHqlPh3gm+VqxYwT777MO9995L8+bN6d69e0bHKLmjcFpql6jlJ5T5EpEcqk5RVHedrwULFnDQQQdxySWXZGJYkicUfEn+qk4ApeBLRHLIL/MV9YlFFVmt+yJNO4rkjYoKqB/hr62edhSRHHIHS8OGDeOnn37igAMOiNRXwVfdp+BLahdlvkSkFnCmHXv27MnkyZMpLCyM3FfBV92naUepXaJmr/xuVLp5iUiWdOvWjfvuu4+bbroprcALoF27dtx9991JbQq+6hYFX5K/tOZLRGqp9u3bM2rUKAYMGMDAgQPp27cvX3/9daS+DRo0oF+/fkltCr7qFgVfUrtEzXxpzZeI5NgVV1xB//79mTFjBnPnzuWOO+6I3Ldr165J1fEVfNUtWvMl+UuZLxGppZYtW8Ytt9yS1Ba1/MSPP/7IH//4R5YtWwbALrvswkMPPZThEUouKfiS2kXBl4jUAgsXLqzUFjV7tW7dOu677z6aNWvGlClTaNOmDbvttlumhyg5pGlHqV007SgitUAm6nytX7+eiooKBgwYkNGxSe4p+JLaJWr2Sk87ikgO+QVaUacd3aUmRo8eze9//3v+9Kc/ZWxsknsKviR/VSd7pWlHEcmhTFW4X79+PQ899BAvvPBCxsYmuafgS2qX6qz50rSjiGSJX6D10EMPsXbt2tC+KrJa9yn4kvylzJeI1FLuYOmll16iTZs2QGwaMZ2+qdqk9lLwJbWLnnYUkVrAWku9evUYMWIERx99NCtXrgRg+vTpoX0bNWpE48aNk9qirheT2kHBl9Qu1dleSDcvEcmS4cOHU15ezqRJkygtLaV79+4ADBo0KLTvbrvtxscff5zUpsxX3aLgS2oXZb5EpJZYuXIlkydPpmnTpnz66adAtOALwBiT9F7BV92i4Evylyrci0gtdsEFF3DGGWcktdWvH17b3FrLL3/5SyZNmpRo69OnT8bHJ7mjCvdSu6jIqojUAs8++yyPP/54pfbly5eH9p0/fz777LMPAK1atWLo0KH885//zPgYJXcUfEn+UuZLRGqpkpIS3/aNGzeG9nVKTfTp04d58+ZldFySHzTtKLWLSk2ISC0QtEarvLw8tK8TfM2fP59jjjmGrVu3smXLloyOT3JLwZfULtXZXkjTjiKSJUHBV5SSEe4iqy+//DJFRUXsvvvuGRub5J6CL6ldlPkSkVogKMh68cUXQ/uqwn3dp+BL8pfWfIlILeUOlgoLCxOv05l2DLqe1H4KvqR2UfAlIrWAEyztvffeLF26NNFeUBD+z66Cr7pPwZfkr0zv7ag1XyKSJb169eLcc8/lT3/6Ex07dky0Rwm+unbtyrXXXpvUpuCrblHwJbVLdRbc6+YlIlkycOBAJkyYQO/evbnjjjsS7RMmTAjtu/POO3POOecktWlvx7pFdb6kdrn8cvjoo/DzNO0oIjn2wgsvcP7557Ns2bJEW6tWrSL1bdSoUdJ7v6lIqb2U+ZLaZfZsWLAg/DxNO4pIDi1fvpwbb7wxKfCC5MX3QX744QcmT57MIYccAkCLFi245557In3up59+yvfff5/2eCW7lPmS/BX0TS9ChWhlvkQklx566CFmzpzp2963b9+UfZcsWcK5554LwFVXXcVBBx3Er3/969DPXLduHXvvvTegTFm+U/AltU+EBasKvkQkl4LWaH377behfZ3AqWHDhvz2t78NDdYcThX81q1bRxyl5IqmHSV/BX1zq2rwpWlHEcmSTFS437x5MwcffDD/+Mc/Im2sXa9evZSfLflDmS+pfYwJP0dPO4pIDgUFQBs2bAjt654y3Lx5M+effz6NGzeu9ASkl1PGQk9G5r9ImS9jzGBjzOfGmCXGmCt8jjcwxjwZPz7LGNMl3n6EMeZjY8zC+M/DXX36xtuXGGP+bkyUf1Flh5LpzJeCLxHJkqDgq6SkJLSv33qtKGu4li9fDsTWfkl+C/1XzBhTD7gHOBroDowwxnT3nDYKWGOt3R24E7gl3l4CHGet7QmcAUxy9bkXGAPsEf8zuBq/h9RFQcFSlIWkmnYUkRxyZ5/cBVOrWuE+SvCl6cbaI0rmaz9gibV2qbW2DHgCGOI5ZwjwcPz1M8BAY4yx1s6z1q6Ity8CGsazZO2AZtbaD23sb9QjwPHV/m2kbgkKliLsjabMl4jkkhMI3XLLLYwfPz7RXpPbC2kCqfaIsuarA/Cd6/1yYP+gc6y15caYdUArYpkvx4nAPGvtFmNMh/h13NfskObYpa4LutlEyWAp+BKRHLrwwgs55ZRTaNeuXVK7U7srlYEDB7J582YaNmyYaFPpiLolSubLL5T2/i1IeY4xZm9iU5Fnp3FNp+8YY8wcY8yclStXRhiu1BlBQVaU4Msv0NK0o2SZ7l87rvbt29OnTx9effVVdtlll0T7rbfeGqm/N4sVJfiqXz+WT9lzzz3TGKnkQpTgaznQyfW+I7Ai6BxjTH2gObA6/r4j8BzwO2vtV67zO7r6+10TAGvtv621/ay1/dq0aRNhuFJnBGWqNO0otYTuXzu222+/nd/97nf8+OOPibYoFe6d84444oi0+jlTk07JCclfUYKv2cAexpiuxpgiYDgw1XPOVGIL6gGGAjOstdYY0wKYBlxprX3fOdla+wPwszHmgPhTjr8Dnq/m7yJ1TXUyXwq+RCSHJk+ezKWXXlqpPcrTjh9++CF9+/bltddeA+Dmm29mY4SdPZzgK8q6Msmt0P9C1tpyYBzwCvAZ8JS1dpEx5jpjzG/ip90PtDLGLAEuBpxyFOOA3YFrjDHz4392jh87F7gPWAJ8BbycqV9K6ohMZ7407SgiWfLhhx/6th911FGhfdetW8e8efPYb7/9+Prrrzn77LND+wC0bdsWgIULF1Ie5T4pOROpyKq19iXgJU/beNfrzcAwn35/Af4ScM05QI90Bis7mExnvrRgVUSyJGiNVjoV7j/66CP+8pe/MHHixEif2bZtWwoLC9m6davKTuQ55SYlf2U68xUl+ProI+jaFV5WIlZEqq462wu5+95///307NmTHj16RAqoVOW+dtD2QpK/Mv20Y5Tg67e/he+/h2OOUaZMRKosE5kvx6JFi1Je07Fq1arE5trKfOU3Zb4kf+WiztfWreHniIiEyMTG2lHbHZ9//nlanyO5o+BL8lemK9xHyWSpQrSIZIATfHXt2pXXXnuNTz/9FICvvvqKNWvWpOxb1eDLHfAp85XfFHxJ/sp05kvTiCKSJV27dmW//fbjzjvvZNCgQTRp0iRx7Mgjj0zZt3PnzowZM6ZSu4KvukPBl+Qv7e0oIrXU1VdfzaxZs+jVqxfvvvsumzdvThybM2dOyr69e/fmX//6V6X2dIKvRo0apTliySYFX5K/qpP5quqCexGRDFm6dClHHXUUhxxyCM888wx9+/at1vXCslnOOq/DDjtMwVeeU/Al+SsXdb605ktEMqC8vJxXXnmFL7/8EoiVgNi0aVOkvqtXr2b27NmJTbibNWvGeeedF7ptkLYXqj1UakLyl/Z2FJFa6swzz2Ty5MmJ91u2bKFTp058+umn3H///Sn7zpgxg2HDYnXLjzzySC6++OJIlfGd4Ov111+ntLQ0aZ2Z5BcFX5K/tOBeRGop7/qsiooKXnnlFYqLixk5cmSkvieeeCJPPfUUJmJG/rDDDqNx48Zs3LiRkpISBV95TNOOkr+U+RKRWsq7PquwsBBIr87XlClTqFevHq+88goffPBB6H6NRUVF7Lzzzr6fL/lFwZfkL2W+RKSW8ma+nOBr48aN3HbbbWn1Pfroo+nfvz+lpaWhn6vthWoHBV+Sv4KCpSiZLz3tKCI5FJT5Arj88st9+4wfP54BAwZQVlbmezys1MT777/P0qVLfT9f8ouCL8lfudheSE87ikgGuIOfX/ziF5xxxhmJ90FlIK6//nree+89XnnlFd/jYcHXjz/+mHitzFd+U/Al+UvTjiJSSzmB0r/+9S8WLFhA69atE8caNGiQsm9Q1ios+HIHXMp85TcFX5K/vDePQYNiP++8MzwAU/AlIjk0duxY7rvvPo444giaNm2adCwo+PrFL34BwEUXXeRbBT8soNL2QrWHSk1I/vIGS7vuGvu5di1Mngy/+130vqCnHUUkawYOHAjAxIkTefvttxk9enTiWHFxsW8fJ2Bq1apVIhBzi7q90K9//Wu6d+9epXFLdij4kvzlDZaKira//uqr9PqCKtyLSFZ9+umniQ2yDz/88ER70JovJ7hynlgMOh7EmXZs3ry5qtznOQVfkr9SBV+NG6fuq2lHEcmhadOm8dRTTyXeFxQU0KRJE0pLSytNQzq++eYbAB5//PHEU4sAt956K4MHD6Zly5YpP1PbC9UeCr4kf3mDpeoGX5p2FJEsmTBhAi+99FLifb169VizZg2bNm2ifn3/f3r33XdfZs+ezVVXXQXAbrvtxi233MJBBx1Eu3btQj9z1/jSjEmTJnHFFVdo6jGPacG95Ce/4MlVJ0eZLxHJZ94F7wUFBRQWFtKsWTMaB9y/vNON3377LT/99FPg+V6HHnooBxxwAABr166twqglWxR8SX7yy1JlY9pRa75EJAO867OC1nG5LViwIOn91q1bGTt2LKeddhpnnHEGa9asCb2GM+WoOl/5TcGX5Ce/4Mud+WrYMP3+mnYUkSzxy3zNnj2bQw89lEsvvdS3z+bNm33bX3zxRR555JHQ7YXWr1/Pd9995/v5kl+05kvyk1+Wyv1Nzi9DtWULDBsGxx+vaUcRySlv5mvnnXdm7dq1vP3220lbDVXnml6TJ0/m22+/BZT5ynfKfEl+8vvWtnXr9td+N6FHH4UXXoBRoxR8iUhOOZmnZs2ace6553LYYYdh4l8aq5qVUpHVukOZL8lPfjcO92azfoHUxo2pj2tvRxHJkoKCAgoKCnjuuecSNb6c4CssgxUkne2FlPnKbwq+JD+FZb78jrsDJ2W+RCSHXnvtNQA2bdrExo0badCgQWjw1bJlS1avXh14zagV7gHatm2b7pAlizTtKPnJ7yYTlvlyP02kOl8ikmPbtm2jX79+NGnShOnTp4cGX+7gqU2bNpWORw2+LrroIvr06VPVYUsWKPMl+am6ma+qbi8kIpIhpaWlfPrpp8D2aUgIDqLctbl69erFunXrmDNnDp06daJDhw4Uucvt+HCmGqOUtZDc0n8hyU9+wVO/fttf+9283MHX669XPq46XyKSJWeccQb9XPesgoICdt55Z04++eSkfR7dunbtmnj9xhtvMGfOHFq1asXzzz/Phx9+SMeOHVN+ppP5+uqrr/j5558z8FtITVHwJfnJL/gaOXL767DgK+o1RURqwJdffskXX3yReF9QUMBee+3FE088wZ/+9CffPt6pwlGjRlFSUsKvfvWrSJ950kkn0bRpU/7zn//w6quvVn3wUuMUfEl+8guu6tWDESNir8OmHaNeU0SkBlSlwr33nK3xpRZbt25ly5YtoeUjunbtyuDBgwE97ZjvFHxJfgq6yTg3p7AF934UfIlIlngDJWstmzZt4osvvkgUQvWaMmWK7zV69uxJw4YN+fzzz0M/19leaN68eWzZsqUqQ5csiBR8GWMGG2M+N8YsMcZc4XO8gTHmyfjxWcaYLvH2VsaYN40xG4wxd3v6vBW/5vz4n50z8QtJHREUfDnZrapkvlTnS0SqqKKigqOPPpqrr76a8vLytMo+OO/nz59Pt27dOOmkkyJ/JkSvD/b666/z5JNPAnDzzTdz/PHHR/ocyb7Q4MsYUw+4Bzga6A6MMMZ095w2Clhjrd0duBO4Jd6+GbgG8N/ICk611vaJ//mpKr+A1FFBNxknOKrKmi9lvkSkir744gumT5/OjTfeSGFhIRvdRZ19eAOl3XffPe0iq96nF8P6vfPOO0nvp0+fHulzJPuiZL72A5ZYa5daa8uAJ4AhnnOGAA/HXz8DDDTGGGttqbX2PWJBmEh0VZl2VPAlIjWkuLg46X3UrX5uvfVWXn/9dXbdddfIQZTDm/lKZ3shILH+S/JPlDpfHYDvXO+XA/sHnWOtLTfGrANaASUh137QGFMBTAH+Yqu654LUPbmadhQR8eGspXKE/XM1bNgw9t9/f4YPH06nTp2A6EGU48QTT0zql872QgB77LFHpM+R7IsSfPn9i+b9GxDlHK9TrbXfG2OKiQVfpwOPVPpwY8YAYwB23XXX8NFK3aDMl9QBun/VHatWrUp6HxZAXXnllQA8+uijrFy5klNOOSU0iMED/rwAACAASURBVCooKEi67iGHHJJoT9UvaEwNGjRIeb7kTpRpx+VAJ9f7jsCKoHOMMfWB5kDwBlWAtfb7+M+fgceITW/6nfdva20/a20/v+0WpI4KW/Pld+PT046SZ3T/qjsWL16c9D5q9ur000/n4osvZuHChaFBlLv9oosuokOHDkDVpx3vvffeSGOU7IsSfM0G9jDGdDXGFAHDgamec6YCZ8RfDwVmpJpCNMbUN8a0jr8uBI4FPkl38FKHhU071lSRVT3tKCI+trq3NyM8C7VgwQJmzZqVeL9+/frQIMp9zfvvv5+ZM2cCcO211/Lggw+GZk+9046lpaUpz5fcCZ12jK/hGge8AtQDHrDWLjLGXAfMsdZOBe4HJhljlhDLeA13+htjlgHNgCJjzPHAkcA3wCvxwKse8DowMaO/mdRumZx2bN0aSkqU+RKRKvMGX2FZqBEjRvDZZ58lnb/77rszY8YMmjZt6tvn3HPPTWSr1q9fz9NPP80BBxzAkCHeZ9z87bTTTknvBw0aFKmfZF+kjbWttS8BL3naxrtebwaGBfTtEnDZvtGGKHXKli3wt7/BkCHQrVvweZlacP/OO9CgAey/f/rBl7XKhIkIAPXrJ/9zGRRAObyZsWbNmlFcXMxhhx0W2Me7qP+///1vWmO85pprWLlyJf/4xz8AVbnPZ5GCL5GM+O47cNLml1+eOhgKOpZu5qtBg+190n3aceedYd48CNnMVkTqPnfmq1+/fjRq1Cjl+U5m7IYbbmDVqlUcccQRKc+31lYK2Jzg6amnnuLHH3/kpJNOom3btimv496iyPuQgOQPbS8k2XPCCdHPzVTmq6Ag9TqxVEpK4K670usjInWSs1XPWWedxezZs0PPd4KvoUOHcvvtt2OM4X//+x+XXnopt912W6Xzt27dyj333JPU5gRft912G+effz7ffPNN4OetWLGCq666ir/97W+JtnQzZ5I9Cr4kez7+OPq5VVlwH3R+Opkv73U90wAismPatGkTABMnTuSee+7h559/Tnm+k8XyZqJuv/12HnzwwcDz3ZwALkqpibVr13LTTTeF/BaSLxR8SX6qyrRj0FRkOgGbd41Efc3MiwiMHTs28XrcuHGsWbMm5fnewMn92i+I8mtLp8J9WVlZpbYmTZqkHKPkjoIvyU9VmXYMCr5SBWxhn6vgS0SAwsLCpPdRa24Z13KIVEGUX1s6G2svWLCgUtujjz6acoySO/qXRfJTVUpNBK0DSxWweSnzJSIRhAVfL7/8Mlu2bKF9+/aJtlRBlF/btddeC0SbdvQGh6AK9/lM/7JIfspk5iudaUdlvkTEh3cxfFjwtddee1VqSyf4OvTQQ+nZs2dSv1SfWV5eXqmtIGzXD8kZ/ZeR/BS2vVDUzFdBQXrTjt7MlxbciwiVtxcKq3DvJ1Xw5Q2sRo0aldheqKioyDez5ea35mvw4MFpj1GyQ8GX5CfnRtS8eexn/CZUrQX3YdOOy5aBdxGtMl8iAmzevDnpfVjm66KLLmLkyJFJC/OLiorYfffd6dKlS6XzGzZsmNR++umn869//QuAGTNmUFZWxoABAwI/zx189erVCyC0JpjkjoIvyU/Oje2Xv4SlS+GLL2LvUwVSYWu+wr6p/r//V7lNwZeIkH7w9cQTT/Dggw8m9evcuTNffvklr7/+eqXzi4qKOOuss5La3nrrrcjjcxeBdep76WnH/KXgS/KTc2MrKICuXaFx4+3vIb3MV9Rpx2+/rdymaUcRYXuRVYAePXr4ruly83vaMYx3jda8efMi93UyX0cffXSiTdsL5S8FX5KfnEDJu2C0OpmvdLcXco9DRHZoTgbr0EMP5bLLLgs936/IaiqbNm1KTDM6nOBp9OjR9OjRg5kzZwb2P+uss/jiiy+48cYbE20rV66M9NmSfQq+JD+5M19uVVlwX9XthQB8niASkR2PM6135plncuqpp4ae75f5Wr58OcXFxb5Zsw0bNrBs2bKkNif4WrZsGYsWLWL9+vWBn9eiRQv22GMPunXrlmjbuHFj6DglNxR8SX5yAilvyr4mpx39KG0vIsAhhxwCxIKv1q1bV3r60csv82WtZcOGDb5bE6V6AtIJ4I466ijfkhJuYRt+S35Q8CX5KV+mHZX5EhHgyiuvTAQ2a9euDc0qpapwH6XUBFSucA/w+eef+37eww8/zPDhw5k2bRq///3vAVW4z2d6lEvyU1WmHatbZLVevcqZLmW+RCTO2Vwbwp927Nu3L+vXr6e+64npdPd29BZZDToPYO7cuTz55JMceOCBPPDAAzzwwAMpxye5peBL8lNQ8OW8TyfzFXXasWFDKC1NblPmS0SAb775Jul9WJFVv3ISqSrV+13vmGOOAZKnLps1a+b7ec7TjmHFWCU/KPiS7Eg3gxS05qs6ma+waccGDSoHX8p8iQhwwgknJL0Py3z5SZX58l7v+++/T+wL6c58tWrVyvfaTvBVVFREeXk5EydOxBjDOeeck/Y4peZpzZdkh6dAYaigNV/V2Vg7LPPVokXlNmW+RITKeyeGBV9bt26lvLw8KdBKJ/PVsGHDxOvjjjsOiBVNrR9Q+Nl5GrOoqIiKigrGjh3L+eefn3KMkjvKfEl2uAoURhK25ivqxtrWRp92dM5r1gycR7qV+RIRkivIQ3jw1axZMzZv3kxpaSmN40WimzZtym233ZZ479apUyceeOABRo4cCUCDBg0Sx84991z2339/fv75Z8rLy5OOOdzTjvXixaFVZDV/KfMl2eG5cYWqyrSj382womJ7nx9+gOnTgz/T2Rutd+/tbcp8iQiVM1+77LJLyvP9Sk00btyYSy+9lLFjx/r2cZ/rDbDOOOMMDj30UJYuXerb1z3t6ARf27Ztq9IG4FLzlPmS7Khq8JXOtKNfW0VF8jWOPjo4A+YEX+efD+++u72/iOywrr32WurXr58UfB177LHstttuKftVZXshJ2g69dRTk6YXP//8cz755BOgchDo6NGjB2vWrKFt27YYYzDGYK1l27ZtietK/lDwJdmRbvCVqTpf7sxXGCf4OuwwuPxyuOUWBV8iO7Dy8nL+/Oc/A9CuXTsAvvjiCzp16hTa1y/zVVZWxpQpUygsLGTo0KFJ53/33XecfvrpQOXpwksuuSRpTH6uu+66pPf16tWjvLyciooKBV95SNOOkh25ynxt2xY9+HLWpRUVwa67xl5r2lFkh+UOdNauXQvE1l+98MILrFu3LmVfJ/PlDr42btzIKaecwqhRoyqdv9n1UJI3+HKvL4u6jkvrvvKbgi/Jjkyv+Yqa+apXr3IAF8TJfBUVxfqBMl8iO7CGDRsmyj08+OCDALzxxhucdNJJzJkzJ2XfdCvcu9uefvpp32tBcOZr/fr1rF+/PhFsOQvvFXzlJwVfkh2ZnnaMsuB+9GjYa69omS9rt4+xqAic9RbKfIns0IqLiwHo1atXUnuqpx39yktA+hXu/T4rKPgaNGgQzZs3TwSFzpORTZs2Dbyu5I7WfEl2ZHvasWVLmDgx9jpK8OWMr7Awdr4yXyLC9ory653yM3FhTxHef//9WGt9M19+gZu7ba+99gr8rKBMlrvOl+Q/BV9S8779dvuUXlRVqfPlbnNvsRFl2tE95QjKfIkIK1euZPbs2UDyondInfkyxiTqdXnbITzz5S0z4f6sAw44wPczneArqAir5BdNO0rN+s9/oHNniD/FE1nQmq+omS938BUl8+VebA/KfIkIW+L3hQ4dOvDBBx8kHavK9kKZCL6aNGnie21nOtLZ2/HII4+kV69e/PTTT2mPU2qeQmSpWRMmxH5+/XV6/apbaiLd4EuZLxHxcAIap2aWW6rgq6Kign/+85/Ur1+fs88+O9GeKvjaaaedEq9nzZqVdOzee+9l4sSJ7L777qFjdTJfn376Kd9//30igJT8osyX1KyAb2mhwqYdwzJf7nUPVZl2VOZLZIcXtLgdUq/52rp1K+PGjeOCCy5Iam/YsCHr1q1j5cqVlfq0b9+e119/HYDOnTsnHdtzzz3p3Lkz7733Hh999FHKsTrBl0pN5DcFX1KzfPYwi6QqC+4zkfly0v3KfIns8JyAZvny5UntBx54IIccckhgv6Dq9sYYmjVrlniC0mvgwIHMmjWL+fPnVzr2zjvv8Nhjj/HNN9+kHKuCr9ohUvBljBlsjPncGLPEGHOFz/EGxpgn48dnGWO6xNtbGWPeNMZsMMbc7enT1xizMN7n7yadPRik9qhu8JVOna+gzFdVph2V+RLZ4fllvk477TTOOussmjdvHtjPr7p9mNLSUmbNmkXjxo1p0aJF0rHbbruNKVOmBI4JYlOTjz32GC1btgQUfOW70DVfxph6wD3AEcByYLYxZqq19lPXaaOANdba3Y0xw4FbgJOBzcA1QI/4H7d7gTHATOAlYDDwcvV+Hck7VQ2+gtZ8VSXzlc60o9NPmS+RHZ430GnVqhWTJk0K7edX3d5pP/zwwzHG8OabbyYdW7x4MQcccAC/+tWvmDt3btKxZ555JvE6KJg67rjjkt47GTAFX/kpSli+H7DEWrvUWlsGPAEM8ZwzBHg4/voZYKAxxlhrS6217xELwhKMMe2AZtbaD23sK8IjwPHV+UUkT2V62rGmMl/uOl+gzJeI0Lp166QM16pVqzDG0LJly8C1V5B6U+23336bt956q1K7ky3z6xOlyKqXk/mKer5kV5TgqwPwnev98nib7znW2nJgHdAq5JruSXS/a0pd0LBh1fqFTTtmes2Xc4NyMl7KfIns8Dp27MiE+BPbI0aM4KijjgJgzZo1fB3wBPfatWsT+zR6M1+pVtekmqqMUmT1jjvu4K9//Stl8Sz+8OHDGTt2bGIaUvJLlFITfn9bvP/yRTmnSucbY8YQm55kV2ezY6k9qhq8VHfa8corK/dJxRmnMl+SQbp/1X6NGjUCYptiv/LKK4n2oKcdr732Wv7xj38AlQueuoMvb/X7VNmyKJmvq6++ms2bN3PeeecB8Mc//jH4l5Kci5L5Wg50cr3vCKwIOscYUx9oDqwOuWbHkGsCYK39t7W2n7W2X5s2bSIMV/JKVYOv6kw7XnMNDBpUuU8qzrSjN/Ol4EuqQfev2m316tWJmlubNm1KOhZU5+vFF1+koqKCDz74gJKSkkrHneBq6dKlSe2pMl/uz2rbtq3v53qfdpT8FiX4mg3sYYzpaowpAoYDUz3nTAXOiL8eCsywKYqgWGt/AH42xhwQf8rxd8DzaY9e8l+6ezo6qlNqwlWsEKjatKOT+dK0o8gO6+OPP+aWW24B4NVXX006FhR8OVOOHTr4r6Rx/mn0FkyNmvk6/vjKy6OttYngy1nrtWTJEt5//33WrFnjOw7JrdDgK76GaxzwCvAZ8JS1dpEx5jpjzG/ip90PtDLGLAEuBhLlKIwxy4A7gDONMcuNMd3jh84F7gOWAF+hJx3rpupmvqpSaiKoj/fabt5pR2W+RHZ4TkCzyy67VDoWFHytW7cOgO7du/seD5Jqwf1uu+2Wsq/76UonczZu3DgOPvhgZs6cmdY4JDsi5SettS8RKwfhbhvver0ZGBbQt0tA+xwql5+QuqYqwdfChXBFPH5Pp8J9ULbMq6Ki8jnKfImIR1Uq3Dds2JCff/6Z0tJSDj300EpPNu69994sWrSoUnDWp08f5s2bR2OfJ8Sfe+45tm7dSllZGVu2bKm096N3X0eApk2bAvDzzz8H/4KSM6pwLzUr3WlHa6FXL9iwIfY+nWnHoMyXl98NVWu+RMSloqKCE044AYD//e9/lY63apXqgf6YDz/8sFLb448/DlTOcDVp0oQ+ffrwy1/+0vda48ePp2nTptx+++2Vjvmt93Kq6G9w7qWSV7QyT2pWupmj//43+X1Q5uutt2KBkZOhguiZL78xBT3tqMyXyA7p1VdfDSzrcOONN3Lsscf6HnP3cco+uFWl8nx5eXnifL9+FRUVtGjRIilr5gRfynzlJ2W+pGalm/nyLg71ZrGcwGr5crjppuRj1cl8eacdnSDM5+YpInWfXwaqS5cuvPPOO5x66qmB/ZwF+kEefPBBIFbR3m3x4sWMHDmS2267rVKffffdN9HuNxX61ltvMWnSpKQ9KJ3g68ILL+TWW29NOSbJPgVfUrNSZY78pg693+qCMl8ADz+cfKw6mS/vtGOzZrGf+tYoskNy1ky57b777vTv35/27dsHZq5Gjx4dOHUI8Ne//tW3fcWKFTz44IO8/HLlZ8/c2Stv8LVt2zaGDBnCcccdl1QOw7159+WXXx44HskNBV9Ss7ZsCT4WVigVUgdf3gxXUGFWryjTjk7wtX596muJSK2zfPlyvvvuu5TnOMVV3U466STGjh1LYWEhEydODOzbuXPn0DEceuihSe+jlprwBn3u9+td9yt38OW9huSegi+pOU89Bc8+G3zc72bg/TYZNO3odyyoPIWXN/hauxbGjIm9djJfDRrE9ocsK0sdQIpIrdOpUyd23XXXwCcWAb755pvEayeQOemkkxKlHCZPnuzb75FHHmHlypWhY6jnXq9K1fd2dAdfBx54YOL1sGHDmDZtGgDt27f3Ld4quaMF91JzTj459XG/4CudzJdXVTNf48dvf+2uDt2sGZSUwBdfQI8e0Yq1ikhe2+pah1pRURFYEb60tDTx2inhsGnTpkRw9N577/Hee+9x8MEHJ/UbM2YMW7ZsYeTIkbRv3z5wHN41X1Er3KcKvpYtW5Z4vfPOO9OtWzeASqUpJPcUCkvuRMl8BZWagMxlvla4drbyBl8QK31x2WWprykitYKzLqpp06Ypt+JxBzyrV8d2y/v555+TgqPPP/88sN+9997L9ddfH3j977//3refX+bLnaE7/fTTk46lempySzxrv3Xr1pQ1yyT7FHxJzfELgs48c/vrKMGXd1og1Zqvqi64d6f/XUUKad58+2uf2joiUvs4wZffmi43b1Bzwgkn0LVr16Tgy2+zdKdfutN8Uacd+/btm3Kcjm+++YbTTjsNiK1xW7hwYVrjkZql4EvS9/XXsWm4J59MfZ5fqvu886BJk9hrz0a1QOWAzLXuAoi24D7dzJf7269f5ktE6gwn+CotLU25CN197JJLLuHZZ5+lfv36ScGRd39Gdz/vmq4wzZs3Z999901MFbr98Y9/pEmTJlx88cWVjgUFX+vXr2fevHmJ9yq2ml8UfEn6LrwQFi2C4cNTn1dUVLltzz3BWUux996Vj3tvJN60fqpvk1XNfCn4EtlhOBtfb9y4kZKSksDzgtZZuTNaRZ57nLuPXwYL4IUXXgBi2xC5HXTQQXz00Uf87W9/q9Tn3HPPZdasWXTo0IHp06cnHWvZsqXv53gze+41bJJ7Cr4kfVFrX/llvty1c378sfLxTJSaqE7myz3tqOBLpM5x18LamqIItDujtGTJksTrk10PEq3xFIWOkvVySkykOy05Z84cLrnkksT2RI6CgoLE5tn77rtvot0bfCnzlV8UfEn6olat9wZf//d/4X2cG163bjBoENx3X/LxKAvulfkSkQC9e/emeXw9Z6pF6G3atEm8dp+3//77c9JJJwGwaNGipD5O8JUqsPr222+BWObN27e8vNx3KvS5557jzPh6Wb8xt2vXjssvvzxpMb53g25lvvKLgi9JX1WCr/nzIWTbDWB78NW3L7z2GvTunXw8SqmJdDNf7m+p7uDLp8K1iNRuBQUFiU2xUwVfvV33Hm+GzJlu9O7dWFhYyNatW1Pup7jffvv5tj///PMUFhby29/+ttKxUaNGBY7lxx9/ZOzYsaxatYo//OEPiXZlvvKb6nxJ+qoSfLmn81JxvvUFpe1r4mnHoGlHzzdHEakbnBITUcsvtG3bNvF61qxZPProo8D2Ug4OYwz169dPWcLCyUB16dIlqT3q047e4Ovnn39OFFM97rjj+M1vfgNUru2lzFd+UeZL0leV4CvFzSiJk/kKCr5STTtWtciqO+Byj1PBl0id8+abb/LFF18Aqdd8bdy4kX322QeAE088MdH+/PPPJ157g6909O/fP+l91CKr3mybe23am2++mXjtDeKOPfbYKo9VMk/Bl6QvarG+6mS+ggKoVFOKVS2yGjTt6JTEcHhueiJS+yxfvjzxOlXm67XXXmPu3LnA9ickIXUgtGHDBg444ACOOOKI0HF4g6yoezumCr7uuuuupGP33nsvEKu6371799AxSfYo+JL0ZTLz5b35ZSPz5Q2iok47as2ESK3nBFI77bRTpak/N3dQ417D5W73BkJlZWXMmjWLOXPmhI5j0qRJSZm3qNOO3vIWUSrcR91e6Mknn+S1116LdK5Uj4IvSV8m13x50/bOjSRK5quq2wt5gy93oOf+XAVfInWOU2ri1FNPpUWLFoHnuQMedzFVd7bMO+2YboFVv+ArbNrRWd/lSBV8OU9U3n333cyYMSPlWEpKShg+fDhHHnlk+MCl2hR8SfoymfnyBl/ZWHDvDb7cj3a7tzNS8CVS5zjBl7fIqZcT1PzmN79h0KBBiXZ38HXhhRf69olaw8sdOKWadgwq2Oq9htdVV10FxAK7Z599NuVYtCA/uxR8SfqiBl9Rvv0FZb6iTDt6RS014f1M99Sn+0bmXfMVttatogI+/LDy9UUkbzjTjnfeeSdff/114HlOMOQt2eAOdpp5agFGyXy592Z0X2u//fbjvvvuY8yYMZX6fPvtt/z3v/9lxYoVlY55x+B2wAEH+H6WH/cTmta7p65knIIvSV/U4MudUfIGMo6gzFd1ph3DvnUGBXze197MV1jw9de/wkEHwRlnpD5PRHLGyXxVVFQkVa73CgqkUgUxUTJfc+bMYaeddqp0rd12241Ro0Zx+OGHV+rTpk0biouLOfzwwxk4cGDSsT322INTTjkFIFG/zOEutBoWfLmLynrXsknmKfiS9EUNvpz/s994I3i+PSZUJ/NV1e2Foma+0g2+7r8/9jNsw3ERyZkePXokXqd62jEokLrlllsS66ImTJiQdCzqmi/neFhA5GatZfHixXz11VeVjjkBmVPjy+F+UCCspllRURHFxcVA9UpoSDQqsirpSzfz5a1S71adzFfQ56W75itq5ivsRhm1nIaIZN369esZOXIkI0eO5Nhjj+XFF19MWedr0KBBTJ8+nXbt2iW1t2jRgpNPPplXX3210lONjRs35ve//30is+Xn3nvvTWzo7Q6+Fi9ezIwZM9hzzz19s1/OE4verNSWLVvo378/s2fPplOnTknHZs+enXgdJdD7+uuvadCgAU2CZiokY5T5kvSlm/lKFQylm/mKsrF2rqYdFXyJ5K2bbrqJKVOm8Otf/zpShfv27dtz1FFH0atXr0rHgrYXat26NQ888AC333574HXHjh2beO3Oqs2cOZPzzjuPhx9+2Ldf0Ge+//777Lnnnlx22WVJlfgh+aGCsODrp59+Yvz48dx5550pF/hLZij4kvSUlyc/EZhKWCAFmZ12TFVq4pVXoGVL/8/M1IJ7BV8ieeunn35KvHY2t466vZDbXXfdldjAujrTc5dccklSsJSq1AQEB19OUOU31Tl//ny6desGkHLLI4DVq1czYcKExNZJUrMUfEl0N94I7gJ/YU8zhpWNgMxOO6bKfB15JFx+uf9nugMu16JTTTuK1B3upwydyvWpph1nzZrFlVdembSdEMQyTQ5vILR582YWLFiQ2L4olXQq3EPVgq9u3bqxePFirLU89NBDKcfjXOeLL77wXVcmmaXgS6K7+urkrFdY8JXtacewIqtO3bGgNV877QSnnrq93RtMKfMlUmv16dOnUluqbNDcuXO5+eabmT59elK7O1s2derURCAHsTVTffr0YciQIaHj8QZLqSrcAxTG7y9lZWVJpSBSBV/pcP9ea9asqda1JJyCL6k6dymJVMdT3RTuuMO/T1DAVp3thZysXdC049/+VrkY7PLl0LNn7LUyXyK1lndbnkcffZQBAwYEZr+Cnnb0TlUuXLgwtI+fm2++mcWLFyfeh0071qtXjwsvvJBLL720RoIv95ow916WUjMUfEnVha3/ipL5euMN/z41mflKJ9vWoQN07Rp7rcyXSK21ZcsWjj/++MT70047jQ4dOvDZZ5/5nh+1zpc7GEt3e6GoFe4dd955J7feemtSgBYWfL399tvstddevsVbg8ai4KvmKfiS6rnvvuBjURbcQ3IGLayPO9jzZt7CsmZhwVfQFIQzFgVfIrXWm2++yX/+859K7UGBRlAWyxt8ud9XZ3shYwwFBQWR+3qvETSFunHjRhYvXpx4yCDKWBR81TwFX1I9qb5NRZl2hOTpvLAAKqgsBIQXWQ1a8+UEVUHjdG5qmnYUqbWCnkwMCjScTFTYtGO6ma8HH3ww8dod8IwZM4aKigr++c9/BvadNWsWb7zxRtLvctBBBzF16lQudx4o8ohSVgOSy1KoyGrNU/AlNSfKtCP4l3oIunmlCr7CAregNV9hn6nMl0itt2rVKt92Z7shr6BA6pe//GXSe3dQEyXzdeaZZ7LPPvsknR/V0KFDGTRoUFLZjPbt23Pcccex7777+vaJWk2/T58+jBgxAlDmKxsiBV/GmMHGmM+NMUuMMVf4HG9gjHkyfnyWMaaL69iV8fbPjTFHudqXGWMWGmPmG2PmeK8pdUDUzJd7wWtYAOWeavQGQ2EL7oOmHZ3rBE07Ou0KvkRqrdWrVyde/+IXv0i8DsryFBcX07lzZ1o69QHj/v73v7N27VpGjx4NpJ/5KikpYeXKlUD6wZfz0EA6mal0tjI66KCDGD58eKVK+ZJ5odsLGWPqAfcARwDLgdnGmKnW2k9dp40C1lhrdzfGDAduAU42xnQHhgN7A+2B140xv7TWOn8LDrPWlmTw95F8kqvMVyYX3EPVph2tDd9jUkSyxv2EYP/+/Vm6dCkQnOU555xzOOec7JNPKQAAIABJREFUcyq1FxYW0rx5c7p27UrXrl1p2rRp4tjee+/Nhx9+mLShtZd7A2t3QPTYY49x6623MmLEiMApRL9aX3PnzuX555+nb9++lfZ2hO3B17vvvsuPP/5YqQq+27hx4xg3blzgccmcKJmv/YAl1tql1toy4AnAW8RkCODsifAMMNDEHtkYAjxhrd1irf0aWBK/nuwIoi64Tyf48luc74ia+Qqq81XdBffuhwE2boy+E4CI1Dh38OUOuKq6vumqq65i6dKlSU8RFhcXc8ABB/huSeTVpUsXujpPUgNLlixhwYIFSZthe/nt7zh37lyuu+66SsVgHe6F+M50p59XX32VwsJCjjvuOKZNm8aECRP47rvvQn8PqZoowVcHwP1fYHm8zfcca205sA5oFdLXAq8aYz42xqR+BlZqp6ApxP088Xc6047ugMsbDDnXCQrcwqYdwzJfY8bA+vX+53jHVlwMv/518LkiklXu4Ovpp59OvD7ssMNyMRzOP//8pE27nanI1q1bB/bxy3yFlZpwf8aKFSsCr71161bKy8tZvHgx11xzDeedd15gGQ6pvijBl9/cifcrfdA5qfr2t9buAxwNnGeMOcT3w40ZY4yZY4yZ4/zllByIUs8rqN17U3j3XfjuO9h119j7TE07lsRnsFu18u9b1QX37ozYDTf4nwPJv4e18PLLwefKDkH3r/zxyCOP+LZ37NjRt/3mm2+mdevW3HbbbZE/4/PPP+e8887j73//e+i53tIQJfH7l3ta0ssJvtyZu7Dgq3Pnzpx22mmh43Gus2TJEubNmwfoqceaFCX4Wg64V991BLzhc+IcY0x9oDmwOlVfa63z8yfgOQKmI621/7bW9rPW9kv1l1JqmHeqzi3oyZigBfdFRdCxo/9C9nQW3LuDL2vhxx9jr3fe2b9vdacdIVbxPkiai2el7tP9K380atQoqZwCbN+yx09paSmrVq0KXBN211130bJlS6699tpE2/Lly5kwYQJTp04NvO51110HxDJx33//faI9SuaruLgYIGlqsiYq3DuWLFlSrWtKsCjB12xgD2NMV2NMEbEF9N6/WVOBM+KvhwIzbCzHOxUYHn8asiuwB/CRMaaJMaYYwBjTBDgS+KT6v47UmFSPHgc8qh264N4JdtzTjmFZKNcaiaRgZ/36WFDVtGnlDbEdTntpaXJ71GlHSL2IPmxNmIjkVPPmzROvjznmGLp3786MGTN8zw0rG1FWVsaaNWsodd1PopSauOaaaxg8eDAXXHAB//3vfxPtTvCVKkifMGECS5YsYeDAgYk252nLVEVW+/Xrx8EHH8ykSZMCr+0XfFU3oJNgocFXfA3XOOAV4DPgKWvtImPMdcYY59GK+4FWxpglwMXAFfG+i4CngE+B6cB58Scd2wLvGWMWAB8B06y1ybuXSn5JlX4OCszCAinnW2c6047du2+vqu/uF5b1AmjRIvbTu2ls1DpfkDr4UuZLJG/94Q9/SFoIv//++7NgwYLA4CusbIRf8dKgwqxefuUfokw7du3ald122y0pgxeW+fr222+58MIL+emnn1JOP/oVYU218bhUT6T/Za21LwEvedrGu15vBoYF9L0BuMHTthTone5gJYeqkvkKm0KsyrQjwPHHw+jRycGOE3yleIyaxo1jAd/mzbE/zg0sap0vSD/zpZITInlh5syZzJkTKylZWFhIo0aNgOB1TWGBlF/wFXUK0C/4Gj16NCtWrEg57eincePG7LLLLklZvbBx+lHmK7tU4V6iSRV8VTXz5Rd8RSlP4Rxz3yycis+pMl/GwE47xV6vXZv+OJ1rBPHLfLmnVEUkZ9xBVlFRUSJ7lO7ejo7qZL6c4+6A59prr2XixImJoNDPyy+/zMknn8wDDzyQaBs7diw//PAD48eP9+3jBFBLly7l//2//xd47X79+lUqRWFVLqfGKPiSaLw3qLPPhn79Yq/DMl9h047plJoA/6DN2Tok7FujE3y5px4zNe3o980y1YMKIpI17uDLGBMafFVl2jHqxtrpVJ13+/rrr3nqqaeYPXt25D7u8b/zzjuB53Xr1o1LLrkk8b558+a+RWYlMxR8STTeG1SDBtun7cIyX+lMO1Y18+VkspzgKohf8JWpaUe/G6ke1RbJC+4g64MPPggNvo455hhuuOEGDjroIN/jfsFXcXExPXv2pEuXLinH4g2+1q9fzxtvvMHChQtT9nOmFte6M/ch0pk6PPzwwxMlOVJV6Zfq02o6icab3WrQAJz0eNjTjtmYdnRuRs6i+iBVyXxVZ82XMl8iecH9VGJZWVmiWnxQ8HXEEUdwxBFHBF6vX79+3HTTTUmL+AcOHJj0BGMQb/C1ePFiBg0aRL9+/VJmtVrE72/r1q1LtN16663ceeed/N///R8XXXRRpT5RF81/9tlnzJw5k86dO/P4449rvVcNU+ZLovFueVFUFJ752rAh9jMo85XJacdMBF9R6nwp8yVSK7lrY5WVlVFcXEzbtm0TtbPc1qxZw9ChQ3n11VcDr9erVy+uuOIKjjnmmLTH8tBDD7F582ZGjBgBbA+mghbNO5zj7uBr3bp1/O9//2Pjxo2+faIGUTNmzGDkyJE888wzbNiwgRtvvJG77747Ul9JnzJfEo13W52mTVNnvi66qOYyXwUFsT/btsUCt8LC6MGXX7mJdOp8paLMl0he2rZtG8OGDWPy5MkAXH755bzzzjv873//8z3/+uuvZ8qUKUyZMoWlS5cm7cGYCU6lesf6+P21WbNmKfs5x9e77sfOtGdQkNWiRQv23ntvFi1alPLaThbu3nvvpX379qxYsSKpCKxkljJfEo0383XYYakzX3fdtf11OsFX2CJ9h7Ow3tmyJWrw5XzLdRdarcmnHZX5Esm5goICHn300cT7Hj16pDy/e/fuidevvfaa7zkrVqzg6aef5t133020TZ48mfr163PGGWf49gkSNfPlTJW6Hx4IK29Rv3599tprr9AxuBf/O3tAlunLY41R8CXReDNf++4bvubLUZUF9yFPC+FsFvvDD7GfUYMvv4AxU9OOynyJ5LVXX32Vs846K3S/RuP6/3lQUDNnzhxOOumkpGtVVFRQUVERWqJh+vTpHHjggVx99dVAepmvww47jAMPPDDpM1ONE1I/fbl+/XpWr17t++Slgq+ao+BLonGCrzFjYtmmgoLwNV+OdEpNRJl2BNhll9hPZ9qgOsGXc4MJCr7cAZd7b0mIFXc95BB48kllvkTy1KZNm1i8eDE9e/bk3//+N02aNOGzzz6jY8eOvrWv3JmlmqjztWnTJmbOnMknn8R21XOCr7DMV7t27ZgxYwYPP/xwpbE6WTGv8vJyioqKGDFihG+h1T333JNWrVolPZDgUPBVcxR8STRO8LXXXtun/Goi8xVlwT1kNvPljL9JE/8+7qDKWzR14kR4910YPlyZL5E8tWDBAvbaay+GDBmSaCsoKOD777/3Xff15ZdfJl4HZZSc9qrU+do5Xgz6p3hxaCeACst8+dkUv38FFWc1xvDoo4/yxBNP+P4uP8TvoXPnzq10rKysDGutgrAaoOBLonGCL/fNIWrmK9N1viA587V2bazCfUEBtGyZup93zOXlsQCpoCD2BKcfd/DlvQm5v6n6ZblKSqBvX7jjjtTjEpEa4zzp2LRp00RbqjpfG5wntala5ivsCUMn+Poxvi3aX/7yFxYuXMhxxx2Xsp+1lrVr1yaCNoDjjz+e8ePH86tf/cq3jzN+ay0lJSVJU6LbXJl8vxpjW7du5aijjqK4uDjpaVGpPj3tKNH4BV9RM19BUpWaCAu+3Jmvl16KBUiHHrp9TEG8wZcz9kaNgtdzuacavZkvdyHCr7+u3PfMM2OL++fOhYsvTj02EakRTjDlLivht3jdscceeyRe18S0Y9v4HrROEGWMCX0IAGIB1E7xcjnbtm3DGMOQIUOSMnpe7vVrbdq0YcOGDTSJZ/ndGa2pU6dSVFREt27dEm2HH344o0aNAmJ7Y6aqeybpUeZLoqlO5itIdRbcuzNfTlHDww4L/0zvmJ3aOKmqOafKfIX97j7rKEQku9LNfI0ePTrxOmx7IfdC9ajTjsXFxTRo0IDS0lLftVZBCgoKEp9b1anAra4vkE7gWVxcTI8ePRJFXCGWjRs5cmSVPkPCKfiScNZCfGFoRjNf1Zl2dGe+nJunT7HESqobfHkzX1pQL5L3nNIJ7sxXquDLOQYEZpX8Ml8DBgzg7rvv5uSTT045HmNMYurxvffe4/DDD2fcuHFRfpXE2JzA6c033+TFF1+MvOWQO/gqKiritttu489//jMQmw696qqraNKkCRs2bOD9999PnNumTZtI15doNO0o4aZNi2WY6tcHd7HB6ma+qlrhHpIzX87nBzztkyTbmS8RyamKigomTJgAQMeOHRPtTqHTLVu2YK1Nmp5zr4XyFkR17LPPPqxatSrpeM+ePenZs2ekcZ1++umUlpayYcMG3nzzzcAK9V4NGjRgw4YNieDrsssu4+OPP+ajjz5i3333De3vDr4aNWrEpZdeytVXX534/R999FF++OEHmjVrxs0335w4d7fddos0PolGwZeEcyojDxsG8bUKQP5kvpzsUy6CL2W+RPLaunXrEqUc3NmlgoICxo8fT2FhIdu2bUuaXrzYtT5z8+bNvptMFxYW0jLsAZ8UbrjhBgDeeustILhUhJd3rVrY045eW73Ze5I36v7hhx8qPfl47rnnRh6fRKPgS8KVlMR+9u6d3B6U+QopMJjgF3w5N4aw4Ktp01hpiNLS2JOO7vGkomlHkR1Ky5YtWbt2LRUVFZXWb1177bW+fdyL8D/88EMGDhwY6bMWLlzI7Nmz6dmzZ6QsFGxfuxWUYfPyTpdGCb723HNPFi9enPR5EAu6Xn75ZT766KNE22WXXZbUd+DAgZx99tls3rw58hglnNZ87cjmzo3VqQoLllativ106ns5gjJfUQMSZ9rRCb62bdv+xGCnTuH9O3SI/Xz99djPqmS+nLGnCr7cTztmetpx7tztG5CLSI2JusE0JAdfmwIy+9999x1HHHEEv/vd7xJtjz/+OKNGjWLatGmhn7Fy5Uo++OCDRJapJjNf7oDLnflatmwZp5xyCnPmzAns+8Ybb9CnT5+k9V9SfQq+dmR9+8Yq1r/5ZurznMyXN/gKynxFDb6czJdzM1i2LJaJatcuvF4X8P/bO/PoKKrsj38fSfAHRHaUJYGEzSQssoQoRgQBAYFBNlHUzIALKMioCAjOMIrjoOggHGUXUBSQJcABHJYAggs7URTDGpU9QBAIAcIScn9/vK7qV9VV3dWBdCdwP+f06Vreq7rdXf3qW/fddx+eflq+aw2Lk8ZLa6DMni9vLvsuXdzLN9PztXKl/A14+DbDBIUNGzZgyZIlHiMOVfFlJ9quXLmCtWvXGkSJ5l1S0zXYsWbNGiQmJmLEiBEAnIuvjz76CIsXL9bj15yIr169eqFy5cqYMmUKqmkPrfBvxKRVd+XUqVP1YH3GP1h8McDvv9vvy8tzj3SsUMG4T/uz791r9J459QaZux212LJ69ZzVN4uWgor5atoUWL1aLvsT8/Xii95tmTtXvm/Z4r0cwzD5ZuHChahTpw5Gjhzpsa9fv37o3r07jh07ZtiuihKrPGCA9WhHTXzFxMT4tOuJJ57QRzwCzrsdO3TogG7duulTETkRX++99x4yMjLQv39/wxRGdp8NAHr06GFYtxJqL774IkaNGuXx/TG+YfHFeJ8sevJkd1egnefr7Fngv/91b/e32/G996QQ2b9frjtouAB4TgdUUOILALSki+anP+046nczciSwfbsxLQfgOS8kdzcyTL6ZP38+WrVqhU8++cRnufT0dMtUDHbpJlRRsnv3bsvjhrnar2vXroGIkJubi/T0dABA3bp1fdofEhKiJ0xt0KCBYbJsp1y/fh1Xr16FECJfAfHexNcLL7xgWDeLLzVTvtWckYx3WHzdjly+bBQR3tI6DBniXrbzfAHAsGHuwHS1IRs0yP7YmvgCZBfioUNyOSrKvo6KWTAVpPjSbDU3Vtq6Oi3IqFFAfLyn+DI/OTqZrmPDBuCxx9xzWGqMHQt07Wo9nyTD3Aa89tpr+Pbbb/H3v//dtsz+/fuxaNEiANCFjoomWDIyMiyTjwLG3GAqmucrIyMDjRo1wtGjR3Ht2jVUrlxZzyDvC20U5eeff+71c6gkJyfj3XffxYEDB1CsWDFkZmbi8OHDhlQZVsybNw/vvfceDmntLNyfs3bt2gBkLq/k5GSkp6ejRYsW+Prrr/WyZvElhNDTT/Dcj/7Dox1vJ06fBpo3B9LTgXbt3Nu9/WkrVQKOHJHL5jgs8+jCQ4eAmjXdwqZ+feDjj+2PbRYn/oovcwNXUKMd1XOZc/FojXSnTsDs2UCtWu7v00p8qTY68XxpWftLlADmzXNv10Tx2rVAhw6+j8Mwtwh5eXk4d+6cwfNix37Nmw7gtBa7qqB5vjp06IDly5ejc+fOAIC3334bU6ZMQVZWliHbvUqY8vD4yy+/4NtvvwVgzCXmC62r0C6o34ovvvgCy5cvR4MGDVCnTh1UNPdI2DBt2jSsX78ezZo1Q40aNQC4xVe9evUMk4lrqN2nVgIrP/YzEvZ83U5MmiSFFwCkpLi328Vo5eQAWl9+RoZn+gez2Fm71ph3y5cYMj+JHjwo310Ng0/y4/nS4iquXJFxapr3ydeTqvb0e/68dXxbmTJS3KamuveZxZfZa+ZPt6PZ82V3TIa5xUlKSkKFChVw4sQJn2UTExNx//336/XMqJns1RF/bdq0wcKFC5GSkmIbS2X2iJ0/fx6lS5f2S3xpnq9du3bh7Nmzjup4m5PSG5o3Tk3mqgkquy5LVXCZxZfa1cjiy39YfDnh+nUZDO40f1Vhxc41bDe32JEjMk4pOtqdUV6lTBnpSdPo31+WdZpx3iy+NM9XQYovIYwpMrQGz6JLwkBoqDwfkfH7UhO8li5ttMlXt+PNiPkKxjWZkuL2hjJMAElJScFcbaCKA8qVK4dNmzbhzJkzljFVa7U0NQD+/PNPxwIIkJ6v1q1b6+stW7ZEVlYW5s+f7/gYmrB76aWX8OGHHzqqo04vdPDgQbRt29ZRl6U2r2W2Eu7Qs2dPZGVlYfr06R7liQjPPPOMvv74448b9m/btg2/ugZjWU3RxHiHxZcThg2TXWg+AjsLPSdPWm+3m9YiK0u+2wmTYsWAjRuBXr3c265ccQfOK5PYWmI+blaWPKZDNzrCwoxxY04DTrXRPufOORdfgFtMaZOMHzniFoxW586v+Nq8GbjvPqMXTUUVXIEWXz/8ALRvb5xmimECxLhx4/Rlp3FVQgjLeC8AeOutt/TlCRMmoHz58jh8+DAmT56McePG6Znx7dDiyQC3KPInEWlycrI+CtNpwLx2npycHJw8eRLr1q3DFgcjpjVPnSq+QkNDUbp0acu4NiGEnoMsJiYGd6uzmwD6dxMVFYWEhARHtjNubl/xZR555o2PPpLvDp9MggoRMH068NNPnvu0UYtmVE/O6dNukaCJL2VosgdCeMaCTZ0q331MLmspeMqW9R6DZkZtgJ2Kr7Jl5fuNiK9z54Dq1d0eIKsuVl/djurToiai8vKAHj2AbduAjh0995uPE+gnTq2RV7P+O6Woe46ZoJKRkYFVq1bp605EzqRJkzB8+HDs27fPcv/bb7+NOXPmGLb973//wz//+U8MHjzYp/gqW7YsiAhE5GiEo5nixYvrgf5ORZsqorQRnGW1Ns1hPadodaZqbbpCluv+0KxZM9uu2ezsbLRv397jO2ZuV/F14QJwzz1Az57+1bt4UT7x20xJUShYvRp44QWgSRNjsPulS25Pimtki2EfIOO7KlUCWrSQ607EF+ApvrSh2ap4sMJOfPmD2s13I+LLSVJXTUzFxnra7q/n6/p1z+7L48fl3JlafJc2bZIZtV6g01XkR3QBcl7Q+Hgencn4hIgsE3pu3rzZsO6ki3D+/PkYM2aM1zxUVatWNayHh4fjzJkziI6O9thnJjMzE3v27EF2djb69OmDunXrYr2vpNUmtNgtp54vLU9XVlaWLr7K+Gqj4e52vOBqM7788ksIIVC/fn1LcQXI0Y8AMHr0aCxdutSwTxOmpc3tnMLEiRORkpJi6L5kJLen+EpOloHnisvYEWfPyqDwwpzRV50m4pVX3N6GRYuk2EhIAJSZ6gG4b+Za/IM2z1d+xdfFi866D626JZ14oFRUz5eT0Y7qOc6ezZ/nywp/xZf5iTonR6aPsBiRpdcdNQr45RdjN7FdvJ4TTp70tMMX+RFfeXnyP/fjj+5BFQxjwzPPPIPIyEgPr9Phw4cN69WrV8eff/5pKdQ0zpw5AwBeRwSqXZLh4eG6l+ypp55CMW9peAA8+uijiIuLw9ixY5GWloYDBw44GoWpMXHiRL0r1annq3z58qhQoQJCQkJ075MTz1e5cuVQrlw5PWO/Ni1SWloa/rDpFdE8X6tXr/YQv9rvM2PGDKxbtw6AnDFg1qxZehnNPsaT21N8KZOI4vp1KVD+9jfgjTeCZ9PNwnyxazdq1zBoPPGEO2GouYwZrfHzJjoAa69RpUq+J8e26l701/Olip5Qh5lTtHOsXQv89ptcvlHxZZULyFx+wwZg2TK5bH5qv3TJ3tMFAFu3StF///03x/OVnS0HUXhL62F1E1FvdE6FmOsGCIATyzI+mTt3Lk6ePKnHG2lo4kubziYjIwNlypTR821ZoXnH7GK+ACni+vbtC0B6wbQ65hgnK1JdvQmjRo3SR0t68wSZUUdsOvV8DRo0CKdPn8Zbb73ll+fr1VdfxalTp1CyZEnMnj3bsM/u3Gr6CfNoRzVprSZYH374YfTp0we/udpVf4To7catLb4OHQL++lc5/Y2Kms/k3DkZt/PFF8AHH3jecOxiw+wuqn37ZLdfZqZv+4iACROAzz6Ty4MGueOl7Jg0Sc7HaGeX+Qau3fi0p5YHHpC5qFS0m7n657p82bnny2p0ooOGC4Cnx8hf8aU2vE5jxbRzjB/v33m9NapWAtQ8GnPoUJkwtWtXz8EbOTnuCcy9kZNjFF/59XxpovPsWetrKTcXaNZMvvLygMGDgWnT5P9Fw85LZ0ZNCeDr+mZuazRPSYkSJdCyZUvDPk183XPPPahcuTKuXbvmc1obTUh58wyVK1cOvXv3BiDzgk2aNAmA9DDlBydCSKOk0kb4E6iv4Y/nC5CTgQ8dOhRvvPEGhigJtO3ObZWWQiNTucfl5OQY5sfUcqBpQvSNW8GxcZO5tcVXz57Al18C5v5m9WZw9qzxpmf2SNjFFVhMVQFAZjqfPh146il7u7ZuBVatAmbNkoLr2WelN27CBN/zAQ4cCHz6qXE+wE2b3KJLu6lqnDkjbd29W+a4atxYdtWpGdm3bJHiT/0ePvvMPeejr8akVSvZrVS/vnubU/Flbjz9FV++vGtWWD0Fq6Mm7fAmvqyEn50YXLrUKPwAOUDC9KRvy83wfKldOlbX+Lp1MkZwxw7gq6+AceNkKhHVi9WggbNzqf+3KVPksRnGAq37Kzw8HFu3bsXvyryzmviqXr26HotVo0YNxMbGWua8unr1Ki5duoSQkBA93smOaIvRu07El9Wk2/54vrRA9djYWDxinqvWATExMejatSvqOZwPV/t+o6OjERkZqW+365ZVvxdz9+7zzz+Pyq4URJcvX9aPXbduXZQqVQorVqxAyZIl0apVK8uJxvP8GfR2C+JIfAkhOggh9gkh0oUQwy323yGEmO/av1UIEaXsG+Havk8I0d7pMW+Yixfd8U+pqcZRYWrKhTNnjOvqjeLrr+WweivsPFuaV03JH2MgJUV2HT36KOBydQMAlCHPtsJOHaWydKnsqvr+eyAxUb5++83YpQpIUaaJoiZN3J6mZcvc309mJjBzJrBrl7vegAHSGwj47nYsVkyO0MuP+KpQwRj75Ss9hRmbqT+8YhZ4Awc6q1eliv/nckqvXvapQMyo3cTZ2VL4rlxpXXbPHutRrqqH1Oq8S5a4l9W8RarXODNTCjNfAtCcINYuhQZz26PdwDMzM9G8eXPDvI2q+FITme7du1f30Fy6dAknXdez2uXoa+qd2rVr480339TX4+LiHHU7rlixwmObP+JL83wlJib6DO7XSEtLQ1RUFFq3bo2kpCQsWbIE3bp181lv+/btaNOmDQBg48aN2L59u74vyib84G0lvtns+YqPj8eAAQMASM+X1tW4f/9+VKxYEZ06dUJERATWr1+vd+tqrFixAiEhIVi8eLFPu29VfIovIUQIgIkAHgUQB6C3ECLOVOw5AGeJqDaAcQDGuOrGAXgSQD0AHQBMEkKEODxm/iECFi40bnMlg0NurrG7JDMT6NPHva6JrwsXpPfK7kZhF5+jelCuXpU3rg0b5PGSkmSXkxWrV7uXzd6rrCzZ9aPeRD/4AHjkEWDNGrmeni5HMZqDT7/+2u1datjQuK9xYykCAeD55+Xk1lY4daOrT081azqrAxhv3qpnxQljxgAPPQQM90O/33WXcV3JHeSVQYOkt9LMQw85P7c3IiKcxa2pHspZs+RIwr/8xfOa/PRTIC5OCvOVK2VM4+XLspz6kJGaCnzzjTuGa8oUY/fg8uXuZaXBBiD/I82aecYaAlLU3XuvHPihwgkZGRvMU9xoXVvXrl3TvVhVqlTB0KFD9amAAOCU69p//fXXUaNGDRw4cACnTp1Cw4YNHXuFGrg8uY8//jjS0tLQpEkTn3XatWuH95UBTGFhYYas+b7QxNcRP5IWh4WF4dChQ/j111/9Sm5qFqBq3Jed+PKW4R4w5hxbZDGAbefOnejYsSMqVqyIJ554AvNcU6T16NEDADDcn3b7VkPLUWL3AtAcwGplfQSAEaYyqwE0dy2HAjgNQJgJ7gR4AAANDElEQVTLauWcHNPq1bRpU/JJbi7RAw8QSQnmfr3wAtHevUTHj3vuU1///jfRggVEU6d6LwcQDR1K9OSTRKtXE+3aRXTiBFGxYtZl+/f3fTztVa8eUfXqRJ06EX3wAVFoKFHTptZlu3f3fqy4OPfy1q2e39eRI77t+eYb3987EdHcue46q1Y5q0NkPFdqqvN6+eXiRaKyZd3n9Bet3jvvyN/nxAn7sp9/7uw3L1VKXrv16vku++671tv79CHasYOoTRuizZuJIiKcX3Paa/Zs/+tor379iMaMIXr5ZWljlSrW5Z591o+vGjvIR7tQVF6O2q/bnN69exMA/dWxY0d9X25uLu3evdtQvn79+gSAOnToQEREXbt21esKIWjevHmOz71y5UoCQO3atfPL5k8++UQ/5/Dhw/2qu2zZMgJAbdu2pfT0dEd1Tpw4oZ8vNDSUxo0b56jesWPHDN8tABoyZAg99thjdOnSJcs6H3/8sV42JiaGoqOjadOmTURE9P7771NcXBwBoH79+lFERITH8RMTE6l27dr6elxcHBERPfjggwSA1qxZ48j2ooxdGybkPnuEED0BdCCi513rSQDuI6KXlTK/usocda3/BuA+AG8D2EJEs13bZwDQ+ke8HtOK+Ph4UuffsiUpSU5yDMhuLSeBzHa8847sQuvfP//HKAxMnmwfT+YrUP3iRd8TTwPSo+jKC4OzZ53Hb2nn79zZ6GUpSMaPB157TS77+A948I9/yO621FRnoyTnzgWeftp634wZ0hP00EOy2zY21nOASLCIjwd+/tnTm3qjtG7tOO5LCJFKRPE314Dg4Lj9ArBp0yYkJydb7hNCYOzYsfr62LFjbQPPmzdvrk8Lc/jwYYw3xxoqDB48WO/Omz9/PrZu3WpZLjIyEq+5/jt5eXmGwG0zvXr10udWdPKZDh48iPbt2+sTYleuXBm9e/dGQkICnnzySY96a9euxSOPPIKwsDAMGzYM4eHhGDFihL4/ISHB9nOY2b59OxISElCrVi3s2bPHMHG2N2bOnInnnnsOffr0wWeffeaojkZOTg6GDBmCAwcOYMGCBY4C569cuYJSpUrhustTvWDBAo+pf+zYvHkzHnjgAX3d1/2/YcOG2OUKR4mIiMDRo0dxxx134KWXXtKvpcTERHTr1g27d+9GZmYmlntpw5s3b446depgyZIlyM7OxrJly7zmRRs9erTuXZs4caLetWmmcePG+tydp06dMngjzQwYMAC1Xbkuly5dqk+IbqZSpUqGa2nEiBG4cuUKhg4diip+hKDYtmFWikx9AXgcwHRlPQnAJ6YyaQAilPXfAFSA7Fp8Rtk+A0APJ8dU9vUDsAPAjurVqzuTmqmp0gNVtSrRl1/m/0k+OpooJ0cec/lyopEjiebPJ5o1i6hmTft648cTLVwovVfq9tBQ+X733URDhhClpBBNm2Ysc+ed3m1q2JBoxgzpbSte3L29USPpoYqJIZozx+i5K15cev3sUL0dvXoRLVkiPy8gP4M/JCfLlz/ce68811df+VfvRsjNJRo1iui77/JXPy/PednsbKLISM/fMiqK6OpVY9kJE+R1u2gRUWysu2x4OFFYmPE3nTJF/uaHD9t7XNu08e+av/NOosaNiVasIMrMJGrd2r1PCPleuzZRkybu7Y895vu4sbHyGgWIatVy/NWhiHu+8tV+EdHkyZM9vAjaSwhhKNuoUSPbsv3799fLbdu2zbYcAPrxxx/1ss8995xtufj4eL3c9evXvR7z008/9fszad4g9dW3b1/b76pt27YEgMqXL08bN2401EtJSaE8h//V3Nxc/busWbOmozpERH/88QctXryYdu7c6bjOjTJz5kz9ezt48KBfdd98800CQG3atPFZdvfu3dSiRQv6/vvv6fTp0x6/S3R0tKH8mjVrvF4P6uuuu+6i2bNney1z/vx5/dgtW7a0Lde7d2+93N69e70ec/369XrZwYMH25arW7eu4bOFh4cTANq1a5df37ddG+YkMdJRAJHKegSA4zZljgohQgGUAXDGR11fxwQAENE0ANMA+eTowF4ZWP7TT0C1ajIFgPaKjJQByocOAXXryjiwxEQZlHzsmJxw+cIF6bE5fVrGQ2n99507y5dGUpKcbPuee2SqijvukJ6QUqXkk70QMgidSKaf+O47oHt3GdNUooScmgaQMWi1asl4rGPHZIzO0qXyGCVKyFicmBg5IvHhhwE1KDM1VR43Pl4mTwXkZwFkvJlWtmlT78HiTz8ts9GHhsqAd80TtXy5PLY/uPry/WLNGjnSr107/+vml5AQ4F//yn99f6ZACg93B6r/8IP8zZs2lQMVzE/XAwfKwQ5CyOtl717peS1bVl47O3fKWRbuvtuYYDY3VwbJ5+UB9erJwRTly8sBIykpsn6pUvLav+sueV3//rscsXjhgrzWqlb1HFyxbJkcTVu/vvTy7d0rr9HixeWxsrPltbxqlUyD0bWrTARbvrz0mt17r7S5Rw95fR06JP+Htwn5ar8gPQSqd8sbgwcPNgz7V2moxHlGRkZ6PWa1atX05V69eiEuzjoMVw1EN3vhzNx33336stPP1LlzZyQnJ+PcuXN6KgVvcVszZ87E4sWLERYWhubNm2PBggW4ePEikpKSLEcj2hESEoIVK1Zg4cKFePDBBx3Xi4qKso2ZKij69u2L+Ph4HD9+HDWsUv14YeTIkahfvz7a2w0mU4iNjcV3332nr2/atMmQbFUL4Ndo3bo15syZg2LFiqF69erYsmULSpcujZycHISGhoKI9Di1li1bomTJkl6vCTX/2IABA9ClSxfLcjExMfpypUqVvB6zlpJqqUuXLobrXsWcG2706NG4du2aPsLzRnHS7RgKYD+ANgCOAdgO4CkiSlPKDATQgIheFEI8CaA7EfUSQtQDMBdAAoCqANYBqAMZD+b1mFb447ZnGKboc7t2OzIMc2tg14b59HwRUa4Q4mXIYPkQADOJKE0I8Q6kO20ZZHfil0KIdEiP15OuumlCiAUAdgPIBTCQiK67DPI45s34oAzDMAzDMIUZR/OxENEKACtM2/6lLF+GjOOyqvsfAP9xckyGYRiGYZhbnVs7wz3DMAzDMEwhg8UXwzAMwzBMAGHxxTAMwzAME0BYfDEMwzAMwwQQFl8MwzAMwzABhMUXwzAMwzBMAGHxxTAMwzAME0BYfDEMwzAMwwQQn9MLFSaEEJkADgXbDgAVAZwOthH5pKjaznYHlsJidw0iqhRsI24Ghaj9AgrP7+svbHdgKap2A4XHdss2rEiJr8KCEGJHUZ1vrqjaznYHlqJqN+OMovr7st2BpajaDRR+27nbkWEYhmEYJoCw+GIYhmEYhgkgLL7yx7RgG3ADFFXb2e7AUlTtZpxRVH9ftjuwFFW7gUJuO8d8MQzDMAzDBBD2fDEMwzAMwwQQFl83iBBiiBCChBAVg22LE4QQHwoh9gohfhFCLBFClA22Td4QQnQQQuwTQqQLIYYH2x6nCCEihRDrhRB7hBBpQohXgm2TPwghQoQQPwkhvg62LUzBwm1YwVIU2zBuvwoeFl83gBAiEsAjAA4H2xY/WAOgPhE1BLAfwIgg22OLECIEwEQAjwKIA9BbCBEXXKsckwvgdSKKBXA/gIFFyHYAeAXAnmAbwRQs3IYVLEW4DeP2q4Bh8XVjjAMwDECRCZwjohQiynWtbgEQEUx7fJAAIJ2IfieiqwDmAXgsyDY5gogyiOhH13I2ZENQLbhWOUMIEQGgE4DpwbaFKXC4DStYimQbxu1XwcPiK58IIboAOEZEPwfblhvgWQArg22EF6oBOKKsH0URaQBUhBBRABoD2BpcSxwzHvKGnBdsQ5iCg9uwgFDk2zBuvwqG0GAbUJgRQqwFUNli1z8AvAmgXWAtcoY3u4loqavMPyBdy3MCaZufCIttReYJHQCEEOEAFgF4lYjOB9seXwghOgM4RUSpQohWwbaHuTG4DQs6RboN4/ar4GDx5QUiamu1XQjRAEA0gJ+FEIB0e/8ohEggohMBNNESO7s1hBB/A9AZQBsq3LlGjgKIVNYjABwPki1+I4QIg2y45hDR4mDb45BEAF2EEB0B/B+A0kKI2UT0TJDtYvIBt2FBp8i2Ydx+FSyc5+smIIQ4CCCeiArDJJ5eEUJ0APARgJZElBlse7whhAiFDKhtA+AYgO0AniKitKAa5gAh72izAJwholeDbU9+cD05DiGizsG2hSlYuA0rGIpqG8btV8HDMV+3HxMA3AlgjRBipxBiSrANssMVVPsygNWQAZ8LCnujpZAIIAlAa9f3vNP1NMYwzI3BbVjBw+1XAcOeL4ZhGIZhmADCni+GYRiGYZgAwuKLYRiGYRgmgLD4YhiGYRiGCSAsvhiGYRiGYQIIiy+GYRiGYZgAwuKLYRiGYRgmgLD4YhiGYRiGCSAsvhiGYRiGYQLI/wOve2uTxIIVBQAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# Calculate the noise-corrected CDF\n",
"R = np.cumsum(r) # CDF of r(t)\n",
"Ez = np.mean(z**2)*(t[-1] - t[0]) # Energy of the signal\n",
"Stilde = (R*(Ez+sigma**2*(t[-1]-t[0])) - sigma**2*(t-t[0]))/Ez\n",
"\n",
"# Preserve the non-decreasing property of CDFs\n",
"mind = np.argmin(Stilde)\n",
"Mind = np.argmax(Stilde)\n",
"Stilde[0:mind] = Stilde[mind]\n",
"Stilde[Mind:] = Stilde[Mind]\n",
"for i in range(len(Stilde)-1):\n",
" if Stilde[i+1]<=Stilde[i]:\n",
" Stilde[i+1] = Stilde[i] + epsilon\n",
" \n",
"Stilde = (Stilde - np.min(Stilde))/(np.max(Stilde) - np.min(Stilde))\n",
"\n",
"# Plot R and Stilde\n",
"fontSize=14\n",
"fig, ax = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(10,5))\n",
"ax[0].plot(t, R, 'r-',linewidth=2)\n",
"ax[0].set_title('$R(t)$',fontsize=fontSize)\n",
"\n",
"ax[1].plot(t, Stilde, 'k--',linewidth=2)\n",
"ax[1].set_title('$\\widetilde{S}(t)$',fontsize=fontSize)\n",
"\n",
"plt.show()\n",
"\n",
"# Calculate the noise-corrected PDF\n",
"sg = Stilde\n",
"sg[1:] -= sg[:-1].copy()\n",
"sg += epsilon\n",
"\n",
"# Plot R and Stilde\n",
"fontSize=14\n",
"fig, ax = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(10,5))\n",
"ax[0].plot(t, r, 'r-',linewidth=2)\n",
"ax[0].set_title('$r(t)$',fontsize=fontSize)\n",
"\n",
"ax[1].plot(t, sg, 'k--',linewidth=2)\n",
"ax[1].set_title('$\\widetilde{s}_g(t)$',fontsize=fontSize)\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Calculate CDTs using PyTransKit package"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlsAAAFECAYAAAAQrfuiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd5hV1d238XuBNI0dNCqoUTCI3aBGo0ZREyQCIqjYNbbEaHwfU9CoKGoSTTTGiA1LiIKg2IIRLCG2xwJiwQJqEFEnoFRB6jCw3j/WzMMEKQNzzuxT7s91zbXP3ucw+3cucPnda6+9VogxIkmSpPxolHUBkiRJpcywJUmSlEeGLUmSpDwybEmSJOWRYUuSJCmPDFuSJEl5ZNiSJEnKI8OWJEmrEEJoFEK4LYQwLYSwb9b1qDgZtiRJWokQQnNgKPAFcAhwYwihSx3/7KYhhC9CCDuu4XMPhRAuqnexKmiGLTUorxIlFYMQwibAPUD/GOOVMcbxQCfgiBDC6XX4Fb8BRsQYP6r1O68PITy5wuf6AZeFEDbOUekqQIYtNZg1XSV6JSipUMQYv4wxnhhjfKHWscUxxv+JMQ5c3Z8NIawPnAXcvcJb+wBjVjjPO8Ak4OScFK6CZNhSg6jjVaJXgpIKRkh+HUL4IISwsLpH/uE6/NEuwDLgperf0ySEUAkcDFweQoghhPdqfX44cELOv4AKxnpZF6DyEGP8EjhxhWOLgf+B/7oS7LrCH90HeH6FP/dOCKHmSvCWfNUsqez9CjgDOA+YCGwF7FmHP3cQ8HqMMVbvLwX2B8YC+wGfAotrfX4M6QKyRYxxYY5qVwGxZ0sNZg1XiV4JSio0nUm97aNijJ/EGF+NMd5ehz+3HTC1ZifGuIwU1L4CXosxfh5jnF3r81OAJsDWOaxdBcSwpYZU+yqxPdANeKb6vVVdCUK6EtwKOLDW7xoD7BtCaJHvoiWVreHA/wsh/DOEcG4IoWUd/1wLYNEKx/YCxtVq42qr6c2yPStRhi01pNVdJXolKKmgxBj/DHwbeJJ0kfhRCGHnOvzRGcCmKxzbE3hzFZ/frHo7fV3qVOEzbKkhre4q0StBSQUnxjgxxng90BEIwO4AIYRdQwgvhRDGhRAuCSE8V+uPvQl0WOFX7QG8vYrT7ApMiTF+kdvqVSgMW2owa7hK9EpQUsEIIfQJIZweQugQQtgJuBKoBJ4LIawH3AucE2Pcg/Qgz7haf/wpYOcQwua1jq0HtA8hbF39dHZtB5HaRZUow5Ya1KquEvFKUFJhaQb0IT1B+DKpPTqsus05BhgdY6x5aGcCtdqq6rmzxgC9a/2+S6v3K4Df1xysnn+wB3Bn3r6JMufUD2oQIYQ+pMlMxwBVwGlUXyVWf+Qp4LoQwuYxxpnVx/7vShBYUD19RA2vBCXlTYzxKuCqVby9O/BWrf1dgEdX+Ew/4KYQwu0xxqUxxsHA4JX8rjNJwe3V+taswmXPlhrK6q4SvRKUVExmAjsBhBAOIj38U3tqGmKMT5LmAWy9ht+1BLggDzWqgISVjz2WGl4IoTNwE9Ahxrh0NZ/7GdA9xviDBitOkqpVP9wzonr3JWDfGOP3MixJBc6eLRUMrwQlFYmFMcZ9SXMALgKGZFyPCpw9W5IkrYUQwpVAL9L406eAS6rnBpRWyrAlSZKUR95GlCRJyiPDliRJUh4V7DxbLVu2jNtvv33WZUhqQK+//vqMGGOrrOvIBdswqbysrv0q2LC1/fbbM3bs2KzLkNSAQgifZF1DrtiGSeVlde1XTm4jhhDuCSFMCyG8u4r3DwkhzAkhvFX90zcX55UkSSp0uerZGgj0Jy3MuSovxhiPytH5JEmSikJOerZijC8As3LxuyRJkkpJQz6NuH8IYVwIYWQIYZcGPK8kSVJmGipsvQFsF2PcA7gZeGxlHwohnBNCGBtCGDt9+vQGKk2ScsM2TNLKNEjYijHOjTHOq349AmhSvZDnip8bEGPsGGPs2KpVSTz9LamM2IZJWpkGCVshhG+GEEL1632rzzuzIc4tSZKUpZw8jRhCGAIcArQMIVQAVwBNAGKMt5MW7PxpCKEKWAj0ji7KKEmSykBOwlaM8YQ1vN+fNDWEJElSWXFtREmSpNpihKOOgr59Yd68ev86w5YkSVJt48bBE0/A7bdD8+b1/nWGLUmSpNqGDk3bnj1hvfqPuDJsSZIk1aishHurVx/s3Tsnv9KwJUmSVOOWW2DqVOjQAQ4+OCe/0rAlSZIEMGwY/PKX6fVVV0GaIrTeDFuSJEkjRsBJJ8GyZdCvXxqvlSOGLUmSVN6efz6FqyVL4Be/gMsvz+mvN2xJkqTy9fbbaU6tRYvg7LPhj3/M2e3DGoYtSZJUnhYsgBNOSBOXHn883HZbzoMWGLYkSVI5ihF+9jMYPx7at4e774bGjfNyKsOWJEkqP1dfDQMHQosWMGQIbLBB3k5l2JIkSeXlrrvgiiugUaM0W/yee+b1dIYtSZJUPu66Kw2EB+jfH7p1y/spDVuSJKn0xQjXXbc8aP3xj/DTnzbIqeu/uqIkSVIhmzsXfv1ruOOO9LThTTfBBRc02OkNW5IkqTQtWwb33Qd9+sAXX0DTpmn/uOMatAzDliRJKj3vvQdnnQWvvpr2998/jdHae+8GL8UxW5IkqXQsXZrGZu29dwpaW22VerNeeimToAX2bEmSpFLx4Ydw2mnLe7POPhuuvx422ijTsuzZkiRJxW3ZsjTofY89UtDaZhsYORIGDMg8aIE9W5IkqZh9/DGccQY8/3zaP/XUFLw22STbumoxbEmSpOLz+edwyy1w440wfz5ssUXqyerePevKvsawJUmSise4cSlgDRkClZXp2LHHwq23QsuW2da2CoYtSZJU+J5/Pi0ePWpU2g8BevSAiy6CAw/MtrY1MGxJkqTCNX06/PKXcO+9aX+DDeDMM+HnP4cdd8y2tjoybEmSpMITIwwdmpbVmTkTmjWDSy6BCy8sqMHvdeHUD5LyorIytZNLl2ZdiaSiM2UKHH00nHhiClqHHQbvvgtXXFF0QQsMW5LyZMAAOOEEOP74rCuRVDRihLvvhg4dYPjwNEfWnXfCM89A27ZZV7fOvI0oKefmzYN+/dLrk0/OthZJRWLaNDjlFHj66bR/1FFw223QunW2deWAYUtSzt1zD8yYAfvtV5BT3kgqNGPGQM+eUFEBm20Gf/lLuoUYQtaV5YRhS1JOVVWlKXAA+vQpmbZSUr7cf3+aAb6yEvbfH4YNS8vtlBDHbEnKqcceg8mToV076NYt62okFbRnn00LR1dWws9+Bs89V3JBCwxbknLsgQfS9rzzoHHjbGuRVMD+/e9067CqKs2j1b8/NG2adVV5YdiSlDOLFsGIEel1z57Z1iKpgM2aBV27wuzZaXvttVlXlFeGLUk58847sGAB7LwztGmTdTWSCtKMGWnerA8+gN12g8GDS74b3AHyknLmrbfSdq+9sq1DUoGaNQsOPTRNUNquXeoK33DDrKvKO3u2JOVMTdjac89s65BUoH73uxS02rdPC0uXwBxadWHYkpQzEyak7W67ZVuHpAI0f36apBTSrcOttsq2ngaUk7AVQrgnhDAthPDuKt4PIYS/hBAmhhDeDiHsnYvzSioss2albatW2dYhqQA9+2wa1LnPPrB3ecWAXPVsDQQ6r+b9I4F21T/nALfl6LySCsicOWm78cbZ1iGpAD35ZNp26ZJtHRnISdiKMb4AzFrNR7oD98bkVWCTEEL59B9KZWLu3LQ1bEn6mkmT0rZjx2zryEBDjdnaBvis1n5F9bH/EkI4J4QwNoQwdvr06Q1UmqRciNGeLdswaTXmz0/bMnj6cEUNFbZWtjpa/NqBGAfEGDvGGDu2ctCHVFQWLIClS6F585KdBHqNbMOk1Zg3L2032CDbOjLQUGGrAqg9xWFrYEoDnVtSAyj3Xi1Ja1DTs/WNb2RbRwYaKmwNB06tfirxu8CcGOPUBjq3pAZg2JK0WmXcs5WTGeRDCEOAQ4CWIYQK4AqgCUCM8XZgBNAFmAgsAM7IxXklFQ7DlqTVqglbZdizlZOwFWM8YQ3vR+BnuTiXpMJk2JK0SjEuv41Yhj1bziAvKScMW5JWqbISqqqgSZOyfILGsCUpJ2pmjzdsSfqaMu7VAsOWpBx54YW0dV1ESV9TxuO1wLAlKQcqK5evxPGjH2Vbi6QCZM+WJNXPXXfB7Nmw++7Qrl3W1UgqOPZsSdK6W7AArrkmve7bN9taJBUoe7Ykad316wdTp8Jee8Exx2RdjaSCZM+WJK2bp5+GP/wBGjWCm2+GsLJVUCWp5nFlw5Yk1d3o0dCrV3p95ZXwve9lWo6kQvbPf6bt3ntnW0dGDFuS1trYsfDDH8JXX0Hv3vCb32RdkaSCtXgxjBiRXnfrlm0tGTFsSVorjz4K3/9+mjH+mGPg3nuhceOsq5JUsIYMgZkz0+PK7dtnXU0mDFuS6iRGuPbaFLAWLIBTTkltaJMmWVcmqWAtWpSeogG46KKyHdiZk4WoJZW2adPgxz+GJ55IbeXvfw+//nXZtpuS6uqGG2DyZNh1VzjppKyryYxhS9JqPfFEClrTpsEmm8DAgdC9e9ZVSSp4L70EV1+dXt90E6xXvpHD24iSVmrGjHSr8KijUtA65BB4+22DlqQ6mDABunZNg+N/+lPo1CnrijJl2JL0X5YtS4PeO3SAQYOgeXO4/noYNQratMm6OkkF7+mn4dBD0xpe3brBX/6SdUWZK98+PUlfM3o0XHhh2kLqzbrzTmjbNtOyJBWDJUvg4ovhT39K+506padoyvj2YQ17tiQxdSqcfjp897spaH3zm2ls1qhRBi1JdTBnTurN+tOfUrj63e9SD9f662ddWUEwbkplbPFiuPFG+O1v09JlTZump7N/8xvYcMOsq5NUFGKEM85IA+Jbt4YHH4T998+6qoJi2JLKUIzw0EPQpw98/HE61r17ekp7xx2zrU1SEVmyJK3X9eijsNFG8NxzNiIrYdiSyszo0an36uWX0/4uu6TerSOOyLYuSUXk88/hjjvSz9SpadK9gQMNWqtg2JLKxCefwCWXpPGqAFtskabA+fGPHb8qqQ6qquDFF9NTMw89lHq1IF2xXXttmidGK2UTK5W4uXNTO/inP6UxWs2apZ6tiy9Ovf6StEoLFsAzz6TbhI8/DrNmpeONGkGPHnDBBemxZZeTWC3DllSiaubL6tMnTUoKcOKJ6SGh7bbLtjZJBWzWLPjHP+Cxx+DJJ2HhwuXv7bQTHHssnHMObLttdjUWGcOWVILefTdN2vy//5v2Dzgg9Wztt1+2dUkqUIsWpYB1330wcuTyW4QAHTumXqwePaB9e3ux1oFhSyoh8+ZBv35pwPvSpWlc1g03pPVfbR8l/Zdly9IYrEGDYNiwNFcWpFuEnTqlcNW9u0tH5IBhSyoRL7yQ1jL89NMUrM47L82ftckmWVcmqaB88UVaQmfQoNRg1Nh7bzj5ZDjhhDSzsXLGsCUVuSVL4Kqr0lisZctSe3n77bDPPllXJqmgfPVVWuj0hhtg/vx0bNttU8A66aS0IKrywrAlFbFJk9Kg99GjU2/WpZfCFVdAkyZZVyapYCxdCn/9K1x2WerVAujaFX75SzjwwHTbUHll2JKK1IsvpvZyzpy0QsagQfD972ddlaSCECO8/TY8/DAMHQr//nc6vt9+qXfrwAOzra/MGLakIjRiBPTsmR4g6tYtTdy86aZZVyUpU8uWwZgxKWA98kjq+q6x3XZpwr3jj/dpmQwYtqQi8/DD0Lt3msz5rLPS+KzGjbOuSlImamZ1f+SRNPHof/6z/L0ttoCjj4ZevdLEo44vyIxhSyoib7+dxrJWVcGvfgXXXedFqlR2Kith1Kh05fX3v8OMGcvfa9MGjjkm/Xzve16JFQjDllQkFixIT2QvWgRnnGHQksrKkiVp0tGHH07L5sydu/y9du3SuIJjjkkTkNowFBzDllQkLrsMxo9PEzj37297KpWNSZPSldaYMcuP7b57Clc9e6aFoG0QCpphSyoC06bBbbel14MGwfrrZ1uPpAYydCice27qyWrTBs4/P4Wstm2zrkxrwbAlFYHBg9Ptw6OOgu98J+tqJOXd5Mnwm9/AkCFpv2dPuPNOHzsuUs5kJhWBxx9P2xNPzLYOSXn2xRfw85/DTjuloNW8eXrkeNgwg1YRy0nYCiF0DiF8EEKYGEK4eCXvnx5CmB5CeKv656xcnFcqB5WV6cnuRo2gc+esq5GUF0uWpDW3dtwRbr45PXJ88snw3nvpNqJjsopavW8jhhAaA7cARwAVwGshhOExxvErfPSBGOP59T2fVG5mzEjt7pZbemErlaS3306PGL/xRto/6qi0ivzuu2dbl3ImFz1b+wITY4yTYoyVwFCgew5+rySWT6Gz+ebZ1iEpx5YuhWuuSdM1vPFGmuX9mWfSuAGDVknJRdjaBvis1n5F9bEV9QwhvB1CeCiE0CYH55XKwsyZaWvYkkrInDnQvTtcfnm6hfiTn8A778Dhh2ddmfIgF2FrZTeS4wr7jwPbxxh3B/4J/G2lvyiEc0IIY0MIY6dPn56D0qTiVxO2WrbMtg6tmW2Y6uTDD+G734UnnkhjA558Ms3tsuGGWVemPMlF2KoAavdUtQam1P5AjHFmjHFx9e6dwEofXo8xDogxdowxdmzVqlUOSpOKn7cRi4dtmFZr8WK45RbYd194//00Gelrr8EPf5h1ZcqzXISt14B2IYRvhRCaAr2B4bU/EELYqtZuN2BCDs4rlQVvI0pFbskSuPvuNJ3D+eenW4hHHw2vvJKePlTJq/fTiDHGqhDC+cBTQGPgnhjjeyGEq4CxMcbhwM9DCN2AKmAWcHp9zyuVC8OWVMQeeQT69IGJE9P+LrvAVVdBjx5O51BGcjKDfIxxBDBihWN9a72+BLgkF+eSyk3NbUTHbElFZsSINPM7pMWi+/WD446Dxo2zrUsNzuV6pAL3wQdpu9VWq/+cpAIyf36ajBTgkktSb9Z6/i+3XPk3LxWwqVNhzJi0YsdBB2VdjaQ6u+kmqKiAvfeGq6+2N6vMuTaiVMDuuittDz8cNtgg21ok1dHcuXD99en1H/5g0JJhSypUn36a2mmAX/wi21okrYWbb4bZs1N3dKdOWVejAmDYkgrQvHlw7LFp26MHHHJI1hVJqpNly+DOO9Pryy/3iUMBhi2p4CxalKbgGTMGtt8+TSwtqUi88gp88gm0bg2HHZZ1NSoQhi2pgFRUwA9+AKNGwZZbpjVpt9wy66ok1dmQIWnbuzc08n+xSvyXIBWIxx+HPfaAF19M0zw8/TS0bZt1VZLqbNkyGDYsvT7xxGxrUUExbEkZ+/hj6NULunWDWbPgyCNh3DjYffesK5O0VqZNSz+bbw577pl1NSoghi0pI/PmwaWXws47w8MPQ4sW6Wnxf/wDXMNYKkJTpqTtNts4MF7/xUlNpQY2dy4MGAA33ACff56OnXQSXHttGlMrqUhNnZq2LvegFRi2pAbyxRdpUulbb4U5c9KxffdNx7773Wxrk5QDhi2tgmFLyrOxY+GOO+C++2Dx4nTs4IPh4ouhc2fvNkglo+Y24tZbZ1uHCo5hS8qDefPg/vtTyHrjjeXHu3eHPn1g//2zq01SntizpVUwbEk5EiO8/jrcfTcMGpQCF8Bmm8Fpp8G558K3v51tjZLy6NNP09aeLa3AsCXV0/vvp3kM778fJk5cfvzAA1PA6tULmjfPrj5JDWDZMnj11fR6r72yrUUFx7AlrYPPPoMHHkgB6803lx/fcks44QQ46yzYZZfs6pPUwN5/P02Ut/XWaZ0tqRbDllRHM2fCQw+lgPXCC8uPb7wx9OyZQtahh0LjxtnVKCkjzz6btgcd5FMv+hrDlrQa8+bB8OEpYD31FFRVpePNm0PXrilgHXmktwmlsjdyZNr+4AfZ1qGCZNiSVlBZmYLV/fenoLVgQTreuHGaquGEE+Doo2GjjbKtU1KBmDdvec9W587Z1qKCZNiSSGNbX3ghBayHHoLZs5e/973vpTVle/WCLbbIrkZJBWrYsHRVdsABPomolTJsqWzFmBZ8Hjw4PU34n/8sf2/33VPA6t0bttsuuxolFYG7707bs87Ktg4VLMOWys7kyakHa/BgGD9++fHtt08B64QTYNdds6pOUlEZOxZeegm+8Q049tisq1GBMmypLMyYkXr6Bw9O7WKNzTeH449PC0Hvv78PEUlaC+PGLR+jdfrpKXBJK2HYUslatiyNWR0wAB59FJYsScdbtEgD3E86KT041KRJtnVKKkJvvAFHHJHm1urSBf74x6wrUgEzbKnkTJsGAwfCnXcun9G9USP44Q9TwDr6aNhww0xLlFTM5sxJc77MmpXmgBk2DJo1y7oqFTDDlkpCjPCvf329F6t16zRm9cc/hjZtsq1RUom49tp0VXfAAenx5aZNs65IBc6wpaK2bBk88ghcc00aPgGpF+uoo9K6hJ07w3r+K5eUC1VVMHQo/PnPaf/GGw1aqhP/N6SiVFWV1ib87W9hwoR07JvfhJ/8xF4sSTm2ZEkam3DttTBpUjp2yimw776ZlqXiYdhSUYkxhazLL18+HqtNG7j44hSyXDZHUs5UVcG998LVV6c5YwDatUsNzimnZFqaiothS0Xjk09Sz9WTT6b9HXeESy5JbZ49+ZJypqoqzXTcrx989FE61r499O0Lxx3navNaa4YtFbylS+Hmm+Gyy2D+fNhkE/jDH+CMMxyPJSmHFi1KXee//z188EE61q4dXHllmpDPkKV15P+qVNDmzUtTNYwalfaPOw5uuimNz5KknPj4Y7j99rTszsyZ6dgOO6SerJNO8qpO9ea/IBWsmqlsXnkFttwyzZvVtWvWVUkqCTHC009D//7wxBNpH2DvveGCC1LIcsZj5YhhSwVp9uw0OfPrr8O226aerbZts65KUtFbtCit2/WnPy1fHLVp03Sb8LzzYL/9XLdLOWfYUkG68MIUtHbcMQWt7bbLuiJJRW3mTLjtttST9cUX6dg228D558OZZ0KrVtnWp5Jm2FLBeeopuO++NI3DyJEGLUnraNmytPL8oEGpUVm4MB3fc0/4xS/SIFAfZVYDMGyp4PTtm7ZXXpkeBJKkOosR3nwzTd3wwAPw2WfL3zvyyBSyOnXyVqEalGFLBeXTT2HMGFh//TRGVZLqZOLE1IM1ZAh8+OHy49tuC717w6mnwi67ZFefylpOwlYIoTNwE9AYuCvGeO0K7zcD7gW+A8wEjo8xTs7FuVVannkmbTt3ToFLklZp4cK0OOpdd8Fzzy0/vsUW6RZh796w//5pwVQpQ/UOWyGExsAtwBFABfBaCGF4jHF8rY+dCcyOMbYNIfQGrgOOr++5VXqmTk3bb3872zokFbBx41LAGjQIvvwyHVt/fTj22DRlw6GHOjeWCkou/jXuC0yMMU4CCCEMBboDtcNWd+DK6tcPAf1DCCHGmolNpGT69LRt2TLbOiQVoGefhT594LXXlh/bZ5/0NGHv3rDxxtnVJq1GLsLWNkCtEYhUAPut6jMxxqoQwhxgc2BG7Q+FEM4BzgHYdtttc1Caik1N2PIpbBUj27A8qaiAX/4yDXiHtGbXKaekkLXHHtnWJtVBLm5kr+yRjhV7rOryGWKMA2KMHWOMHVv5f9uyZNhSMbMNy7GFC9NCqO3bp6DVogVccw1MmQJ/+YtBS0UjFz1bFUCbWvutgSmr+ExFCGE9YGNgVg7OrRJj2JLEggVwxx0paH3+eTrWsyfccIMT76ko5aJn6zWgXQjhWyGEpkBvYPgKnxkOnFb9uhfwL8draWVmVN9YNmxJZWj+/BSodtgBLrooBa3vfCfNdPzQQwYtFa1692xVj8E6H3iKNPXDPTHG90IIVwFjY4zDgbuB+0IIE0k9Wr3re16Vnhjt2ZLKzuLFKUw98AAMHw7z5qXjHTummY27dHECUhW9nDwbG2McAYxY4VjfWq8XAcfm4lwqXdOnQ2UlbLBBGpohqUQtWZIWPX3gAXj0UZgzZ/l7BxwAl16aZns3ZKlEOBGJCsY//pG2BxyQbR2S8uSdd+CWW9ItwZkzlx/fYw84/vg0EemOO2ZXn5Qnhi0VhGXL4M470+tevbKtRVIOLV2arqRuuinNk1Vj553T3FjHH+8sxip5hi0VhDvugFdfhW9+M7W9korc4sVw771w3XXw0Ufp2AYbwOmnw7nnwq67eptQZcOwpcw99xz86lfp9S23OAm0VPQmTYIjjkhbgG99C84/H3784zQhqVRmDFvK1JNPQo8esGhRaoePOSbriiTVS0UFHHYYTJ6cbhX27ZvWLGzcOOvKpMwYtpSJGNMdhrPPTg8mnX023H571lVJqpcvvlgetPbbD555BjbcMOuqpMzlYlJTaa189hl07ZqGbixZAhdemMZsNfJfo1S8pkyBTp3gww/T04UjRxq0pGr+700NprIS+veHXXaBJ55IY7PuvhtuvNFxslJR++gjOPBAGD8eOnSAp5+GTTfNuiqpYHgbUXlXVQX33QdXXZXuLgB07w633gpbb51paZLq691302D4zz+HffZJPVqbb551VVJBsWdLeVNZmcZldeiQBr/XjJd95JE0abRBSypyo0fDwQenoHXooWlWeIOW9DX2bCnnpk1LY7BuvTW1wQBt28IVV8AJJ/hQklQSRo1KXdTz56ft0KHQvHnWVUkFybClnIgRXn8dbrsNBg9O8xkC7LYbXHQRnHwyrOe/Nqk0PPpomv29shJOOQXuucf/wKXV8L8O1cucOSlc3XknvPVWOhYCdOuWnjI89FAHv0sl5Zln0rxZS5emiUpvuslHiaU1MGxprS1Zkh42GjwYHnsMFi5MxzffHE49Fc47L902lFRiqqrSpHhLl6ZlH667zqspqQ4MW6qTGNPahYMHwwMPwIwZy9/r1Cm1vz16QLNm2dUoKc/eegs++QS23x5+9zuDllRHhi2t1vvvp4B1//3LlzmD9FThySfDiSemdldSGXjhhbTt1MkxWtJa8L8WfU1FBTz4YApYr7++/PjWW6enCU86Cfbc04taqey88kraHnRQtnVIRUQASnsAABHzSURBVMawJQCmT4eHHoIhQ+DFF5cf32gj6NkzBaxDDnHaBqms1czl8q1vZVuHVGQMW2Xsyy/TE9xDh6Ypc5YuTcebNYOjjkpPdv/oR9CiRbZ1SioQX36Ztptskm0dUpExbJWZ+fPh8cdTwBo5Mk2TA2n4RZcuKWB17556tCTpv8yZk7Ybb5xtHVKRMWyVgcpKGDEiBazHH4cFC9LxENI419694ZhjXGVD0hoYtqR1YtgqUTUzuv/tb2kc1syZy9/bf/8UsI49FrbaKrsaJRWRZcvgq6/Sa7u+pbVi2CoxU6fCoEEwcCCMH7/8+K67pqkajj/eqRokrYO5c9NV3IYb+qSMtJYMWyWgsjINdB84MM3svmxZOt6yZZoH67TTYK+9nKpBUj14C1FaZ4atIvbFF3DHHWnx55onsps0SQPcTzsNjjwSmjbNtkZJJcKwJa0zw1YReu01uPnmtGxOzdOEu+wC556bJh1t2TLb+iSVoJqw5bQP0lozbBWRUaOgb194+eW0HwIcfTRccAEceqi3CSXlUc0cW/ZsSWvNsFUEXnsNLrkkhS1IF5ZnnQXnnedEzpIaiLcRpXVm2CpgEybAZZfBI4+k/Y03hj59Uk/WN76RbW2SykxFRdo6X4y01gxbBWjJErjiCrjuuvRkYfPmcOGF8Otfw2abZV2dpLL08cdpa3e6tNYMWwXm44/TIPfRo6FRozTovW9f2HrrrCuTVNYmTUrbHXbItg6pCBm2Cshjj6UpG+bOhTZtYPBgOOigrKuSJOzZkurBsFUgHn8cevWCpUuhRw+46y5vGUoqEFVV8Omn6bVLUEhrrVHWBQiefz6tU7h0aRoA//DDBi1JBWT06DSYdKedoEWLrKuRio5hK2Nz56YxWosXw09/Cr//vfNlSSowI0akbefO2dYhFSnDVsauvjotHr3ffmlWeIOWpIJSUQH9+6fXXbtmW4tUpAxbGfr8c/jzn1PA6t8fGjfOuiJJqiVGOPvs1AXftSscdljWFUlFybCVoeHD07jTLl2gY8esq5GkFdxzDzz5JGy6aVr13q53aZ3UK2yFEDYLITwTQvh39XbTVXxuaQjhreqf4fU5Zyl54om07dEj2zok6Ws+/hguuii9vvlmZ46X6qG+PVsXA6NijO2AUdX7K7Mwxrhn9U+3ep6zZNTMEWivlqSCMnMmHHlkun149NFw4olZVyQVtfqGre7A36pf/w04up6/r6zMmJG2LVtmW4ck/Z8FC6BbN/jgA9h9dxg40NuHUj3VN2xtGWOcClC93WIVn2seQhgbQng1hGAgI407nTkzvd5882xrkSQgNUwnnwwvv5yWsRgxAjbeOOuqpKK3xhnkQwj/BL65krcuXYvzbBtjnBJC2AH4VwjhnRjjRys51znAOQDbbrvtWvz64jNvXpojcIMN0kLTkopf0bdhDz4Ijz4Km2wCTz0F22yTdUVSSVhj2IoxHr6q90IIX4QQtooxTg0hbAVMW8XvmFK9nRRCeA7YC/ha2IoxDgAGAHTs2DHW6RsUqZpbiPZqSaWjqNuwTz5ZPiD+D3+AnXfOth6phNT3NuJw4LTq16cBf1/xAyGETUMIzapftwS+B4yv53mLnrcQJRWEBQvgyiuhfXuYMgX23RfOPDPrqqSSUt+wdS1wRAjh38AR1fuEEDqGEO6q/szOwNgQwjjgWeDaGKNhqzpsOTheUmYefjj1YPXrB4sWQe/e8Pe/QyOnYJRyaY23EVcnxjgT+NqUwjHGscBZ1a9fBnarz3lKkbcRJWUmRrjqqtSjBbDHHvCXv8DBB2dallSq6hW2tO4qKtLWni1JDSpG+M1v4NprUw/WjTfCz37memFSHhm2MvLUU2l7wAHZ1iGpjMyYkQbB33cfrLceDBoExx+fdVVSyTNsZWDaNHjxxXQh2blz1tVIKnlVVWltw8svh9mz03wzDz0EP/pR1pVJZcGwlYE//zm1fd26pfVdJSlv3noLTj8dxo1L+4cfnsZnObWD1GB85KSBjRkDf/xjet2nT7a1SCphVVXw29/CPvukoLXddunpw6efNmhJDcyerQY0Z056srqqCn7+c8drScqTDz+EU0+F0aPT/vnnpwHxG2yQbV1SmTJsNZCFC1PQ+vhj2HvvNEGzJOXcCy9Aly4wfz60bg1//Wu6dSgpM95GbABffQVHHglPPpnm1XrgAWjWLOuqJJWcV15Jg97nz4djj4V33jFoSQXAsJVnEyfCgQfC88/D1luni862bbOuSlLJGTs2Pd48bx6cdBIMGZIWlJaUOcNWHg0fDh07wttvw047pekeOnTIuipJJWfcOPjBD2Du3NSjNXCgk5RKBcSwlQdffQXnnQfdu6dB8T16wGuvwQ47ZF2ZpJLzxhvpVuHs2anRGTw4TVgqqWAYtnJs1CjYbTe47TZo0iQNhH/4Ydhoo6wrk1RyXngBDjkkzQzfpUsaENqkSdZVSVqBYStHZs2Cc89NF5iffJKeOBw7Fn71Kwgh6+oklZz77ku3Dr/6Ki258+ijPnkjFSjDVj0tWQL9+0O7djBgQLqo/O1v4dVXYffds65OUsmpqoJf/CLNo7V4cRqzMHgwNG2adWWSVsEb+/Xw5JNpTdcJE9L+oYemVTB23TXbuiSVqNmzUy/WM8+kcVk33ww/+UnWVUlaA8PWOpgwIV1YjhyZ9nfcEW64Ia116C1DSXkxcWKaQ+vDD6FVq7SQ9MEHZ12VpDowbK2FmTOhXz+49VZYujQNeu/bN62E4VAJSXkzYwYccQRMnpzGJwwfntY6lFQUDFt1sGRJerrwyitTL36jRqnnvl8/2GKLrKuTVNIWL063DidPThP3/etfsOGGWVclaS0YttZg5Mg0Luv999P+YYfBjTem6R0kKa8WL4ZevVLA2mKL9MShQUsqOoatVXj//RSyasZltW2bxmV17eq4LEkNoLISjjsO/vEP2GwzePrptLC0pKLj1A8rmD0b/ud/Us/VyJFpXNb118N77zkAXlIDuvjiNDZr003TbMl77JF1RZLWkT1b1ZYtgzvvhMsuS2NRQ4BzzoGrr3ZclqQGFmOaDR7SrcM998y2Hkn1YtgCxo+Hs8+Gl19O+9//Pvz5z7ZvkjLy4YcwZUqa4sHpHaSiV9a3ERcvTk8U7rlnClpbbQUPPgjPPmvQkpShxx5L20MPdeyCVALKtmfrvfegd2949920f/bZadHoTTbJti5JZS5GuOOO9PrUU7OtRVJOlGXP1sCBsM8+KWi1bZt6sgYMMGhJKgCTJsHHH6dbiJ07Z12NpBwoq7C1YAGccUb6WbgwXTS+9RYcckjWlUlStTffTNuOHaFx42xrkZQTZXMbcfr0NEfW6NHQogXccgucfrrDISQVmJqwtdde2dYhKWfKImxNm5Z6ryZMSMuJPf64M8BLKlDvvZe2zqsllYySD1szZ8Lhh6egteuuaRLmrbbKuipJWoUvv0zbVq2yrUNSzpR02KqqgqOPhnfegfbt0yTMTlAqqaB99VXaugaiVDJKeoD81VfD//4vbL21QUtSkZg7N2032ijbOiTlTMmGrfffh2uuSQPgBw1KgUuSCp49W1LJKdmw1bdvWu/wnHPSJMySVBTs2ZJKTkmGrXHjYNgwaNYsLSwtSUWhqipNAtioEay/ftbVSMqRkgxbgwen7ZlnQuvW2dYiSXVWcwvxG99wEkCphJRk2HriibTt1SvbOiRprdSELW8hSiWl5MLWnDkwfnyaJf7AA7OuRpLWQs14LQfHSyWl5MLW7Nlp26oVNGmSbS2StFbs2ZJKUr3CVgjh2BDCeyGEZSGEjqv5XOcQwgchhIkhhIvrc841mTMnbTfeOJ9nkaQ8sGdLKkn17dl6FzgGeGFVHwghNAZuAY4EOgAnhBA61PO8q1QTtrwwlFR0nGNLKkn1Wq4nxjgBIKz+qZl9gYkxxknVnx0KdAfG1+fcq2LPlqSiNWNG2m62WbZ1SMqphhiztQ3wWa39iupjeWHYklS0/vOftN0mb02kpAyssWcrhPBP4JsreevSGOPf63COlXV7xVWc6xzgHIBtt922Dr/66wxbkrJS7zasJmw5QaBUUtYYtmKMh9fzHBVAm1r7rYEpqzjXAGAAQMeOHVcayNakZnypYUtSQ6t3G1ZRkbb2bEklpSFuI74GtAshfCuE0BToDQzP18ns2ZJUtLyNKJWk+k790COEUAHsDzwRQniq+vjWIYQRADHGKuB84ClgAvBgjPG9+pW9aoYtSUXLsCWVpPo+jfgo8OhKjk8ButTaHwGMqM+56sqwJakozZ6dGrAWLWDzzbOuRlIOlewM8oYtSUXl7bfTdtddXYRaKjElF7bGV8/e1bZttnVI0lqpCVt77JFtHZJyrqTC1qxZ8OmnqRe+Xbusq5GktTBuXNruvnu2dUjKuZIKWzVt1W67QePG2dYiSXUWI4walV7vu2+2tUjKuZIKW88/n7bf+U62dUjSWpkwASZPhlatYJ99sq5GUo6VVNh6+OG07dYt2zokaa3cdVfa/uhH0KikmmVJlFDYeu01ePdd2GQT6NQp62okqY4mT4Y77kivf/7zTEuRlB8lE7ZuuCFtzzoLmjbNthZJqpMYU6O1YAEcdxzstVfWFUnKg3pNaloo3nsPHnwwhSwvDCUVjRDgJz+BqVPh5puzrkZSnpREz1a/fukC8eyzoU2bNX9ekgpGr17wzjuwxRZZVyIpT4o+bI0fD8OGQbNmcMklWVcjSevAQfFSSSv624jt26enECdPdu1WSZJUeIo+bDVqBMcck3UVkiRJK2fftSRJUh4ZtiRJkvLIsCVJkpRHhi1JkqQ8MmxJkiTlkWFLkiQpjwxbkiRJeWTYkiRJyiPDliRJUh4ZtiRJkvIoxBizrmGlQgjTgU+yrqMeWgIzsi4ix/xOha/Yv892McZWWReRC7ZhBafUvg/4nQrNKtuvgg1bxS6EMDbG2DHrOnLJ71T4Su37KDul9m+p1L4P+J2KibcRJUmS8siwJUmSlEeGrfwZkHUBeeB3Knyl9n2UnVL7t1Rq3wf8TkXDMVuSJEl5ZM+WJElSHhm26imE0DmE8EEIYWII4eKVvH9RCGF8COHtEMKoEMJ2WdS5Ntb0nWp9rlcIIYYQCvrJkbp8nxDCcdV/T++FEO5v6BrXVh3+3W0bQng2hPBm9b+9LlnUqcJm+1X47ReUXhtWlu1XjNGfdfwBGgMfATsATYFxQIcVPnMosH71658CD2Rdd32/U/XnNgReAF4FOmZddz3/jtoBbwKbVu9vkXXdOfhOA4CfVr/uAEzOum5/CuvH9qvw26+1+HsqmjasXNsve7bqZ19gYoxxUoyxEhgKdK/9gRjjszHGBdW7rwKtG7jGtbXG71TtauAPwKKGLG4d1OX7nA3cEmOcDRBjnNbANa6tunynCGxU/XpjYEoD1qfiYPtV+O0XlF4bVpbtl2GrfrYBPqu1X1F9bFXOBEbmtaL6W+N3CiHsBbSJMf6jIQtbR3X5O9oJ2CmE8FII4dUQQucGq27d1OU7XQmcHEKoAEYAFzRMaSoitl/FodTasLJsv9bLuoAiF1ZybKWPd4YQTgY6At/Pa0X1t9rvFEJoBNwInN5QBdVTXf6O1iN1wx9CunJ/MYSwa4zxyzzXtq7q8p1OAAbGGG8IIewP3Ff9nZblvzwVCduv4lBqbVhZtl/2bNVPBdCm1n5rVtLdGUI4HLgU6BZjXNxAta2rNX2nDYFdgedCCJOB7wLDC3iQaV3+jiqAv8cYl8QYPwY+IDVchaou3+lM4EGAGOMrQHPSmmNSDduvwm+/oPTasPJsv7IeNFbMP6SriUnAt1g+0G+XFT6zF2kwYLus683Vd1rh889RwANM6/h31Bn4W/XrlqQu7s2zrr2e32kkcHr1651JjVnIunZ/CufH9qvw26+1+HsqmjasXNsve7bqIcZYBZwPPAVMAB6MMb4XQrgqhNCt+mN/BL4BDAshvBVCGJ5RuXVSx+9UNOr4fZ4CZoYQxgPPAr+KMc7MpuI1q+N3+gVwdghhHDCE1HA5g7H+j+1XcSi1Nqxc2y9nkJckScoje7YkSZLyyLAlSZKUR4YtSZKkPDJsSZIk5ZFhS5IkKY8MW5IkSXlk2JIkScojw5YkSVIe/X9zs7OJDqSgsAAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# Reference signal\n",
"t0 = np.linspace(0, 1, N)\n",
"z0= np.ones(t0.size)+epsilon\n",
"s0 = z0**2/np.linalg.norm(z0)**2\n",
"\n",
"# Create a CDT object\n",
"cdt = CDT()\n",
"\n",
"# Compute the forward transform\n",
"s_hat, s_hat_old, xtilde = cdt.forward(t0, s0, t, s, rm_edge=False)\n",
"sg_hat, sg_hat_old, xtilde = cdt.forward(t0, s0, t, sg, rm_edge=False)\n",
"\n",
"# remove the edges\n",
"s_hat = s_hat[25:N-25]\n",
"sg_hat = sg_hat[25:N-25]\n",
"xtilde = xtilde[25:N-25]\n",
"\n",
"# Plot s_hat and sg_hat\n",
"fontSize=14\n",
"fig, ax = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(10,5))\n",
"ax[0].plot(xtilde, s_hat, 'b-',linewidth=2)\n",
"ax[0].set_title('$\\widehat{s}(t)$',fontsize=fontSize)\n",
"\n",
"ax[1].plot(xtilde, sg_hat, 'r-',linewidth=2)\n",
"ax[1].set_title('$\\widehat{s}_g(t)$',fontsize=fontSize)\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Calculate time delay and linear dispersion estimates"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"True value of time delay: 0.2575 seconds\n",
"Estimated value of time delay: 0.2694287224582905 seconds\n",
"\n",
"True value of linear dispersion: 0.8\n",
"Estimated value of linear dispersion: 0.8408065181070804\n",
"\n"
]
}
],
"source": [
"X = -1*np.ones([len(s_hat),2])\n",
"X[:,0] = np.transpose(sg_hat)\n",
"estimates = np.linalg.solve(np.matmul(np.transpose(X),X), np.matmul(np.transpose(X),np.transpose(s_hat)))\n",
"\n",
"print('\\nTrue value of time delay: '+str(tau) + ' seconds')\n",
"print('Estimated value of time delay: '+str(estimates[1]) + ' seconds\\n')\n",
"\n",
"print('True value of linear dispersion: '+str(omega))\n",
"print('Estimated value of linear dispersion: '+str(estimates[0]) +'\\n')"
]
},
{
"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.7"
}
},
"nbformat": 4,
"nbformat_minor": 2
}