demo_synthetic.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison of online NICA algorithms on synthetic sparse uniform sources"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from nica_algorithms import bio_nica_direct, bio_nica_indirect, two_layer_nsm, nonnegative_pca\n",
    "from util import permutation_error, add_fill_lines"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load desired dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = '3-dim_synthetic'\n",
    "# dataset = '10-dim_synthetic'\n",
    "\n",
    "S = np.load(f'datasets/{dataset}/sources.npy')\n",
    "X = np.load(f'datasets/{dataset}/mixtures.npy')\n",
    "\n",
    "s_dim = S.shape[0]\n",
    "x_dim = X.shape[0]\n",
    "samples = S.shape[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Online algorithms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "trials = 1 # set the number of trials\n",
    "\n",
    "bionica_indirect_err = np.zeros((trials,samples))\n",
    "bionica_direct_err = np.zeros((trials,samples))\n",
    "nsm_err = np.zeros((trials,samples))\n",
    "npca_err = np.zeros((trials,samples))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trial 1 of out 1:\n",
      "Running online algorithms...\n",
      "PP.T close to degenerate\n",
      "PP.T close to degenerate\n",
      "PP.T close to degenerate\n",
      "PP.T close to degenerate\n",
      "PP.T close to degenerate\n",
      "PP.T close to degenerate\n",
      "PP.T close to degenerate\n",
      "PP.T close to degenerate\n",
      "PP.T close to degenerate\n",
      "PP.T close to degenerate\n",
      "PP.T close to degenerate\n",
      "PP.T close to degenerate\n",
      "PP.T close to degenerate\n",
      "PP.T close to degenerate\n",
      "PP.T close to degenerate\n",
      "PP.T close to degenerate\n",
      "PP.T close to degenerate\n",
      "PP.T close to degenerate\n",
      "PP.T close to degenerate\n",
      "PP.T close to degenerate\n",
      "PP.T close to degenerate\n",
      "PP.T close to degenerate\n",
      "PP.T close to degenerate\n",
      "Computing errors...\n",
      "Final Bio-NICA (indirect) permutation error: 1.3906235143426384\n",
      "Final Bio-NICA (direct) permutation error: 0.027559837414870128\n",
      "Final 2-layer NSM permutation error: 0.038030891897274\n",
      "Final Nonnegative PCA permutation error: 0.07893823057132256\n"
     ]
    }
   ],
   "source": [
    "for i_trial in range(trials):\n",
    "    \n",
    "    print(f'Trial {i_trial+1} of out {trials}:')\n",
    "    \n",
    "    # Random permutation of the mixtures\n",
    "\n",
    "    idx = np.random.permutation(samples)\n",
    "    \n",
    "    # Noncentered whitening for Nonnegative PCA\n",
    "\n",
    "    sig, U = np.linalg.eig(np.cov(X))\n",
    "\n",
    "    X_white = U@np.diag(1./np.sqrt(sig))@U.T@X\n",
    "    \n",
    "    n_dim = s_dim\n",
    "    \n",
    "    eta = 1e-3\n",
    "    \n",
    "    bionica_indirect = bio_nica_indirect(s_dim, x_dim, n_dim, eta0=eta, decay=0)\n",
    "    bionica_direct = bio_nica_direct(s_dim, x_dim, eta0=eta, decay=0)\n",
    "    nsm = two_layer_nsm(s_dim, x_dim, eta0=eta, decay=0)\n",
    "    npca = nonnegative_pca(s_dim, x_dim, eta0=eta, decay=0)\n",
    "    \n",
    "    bionica_direct_Y = np.zeros((s_dim,samples))\n",
    "    bionica_indirect_Y = np.zeros((s_dim,samples))\n",
    "    nsm_Y = np.zeros((s_dim,samples))\n",
    "    npca_Y = np.zeros((s_dim,samples))\n",
    "    \n",
    "    print('Running online algorithms...')\n",
    "    \n",
    "    for i_sample in range(samples):\n",
    "\n",
    "        x = X[:,idx[i_sample]]\n",
    "        x_white = X_white[:,idx[i_sample]]\n",
    "        \n",
    "        bionica_indirect_Y[:,i_sample] = bionica_indirect.fit_next(x)\n",
    "        bionica_direct_Y[:,i_sample] = bionica_direct.fit_next(x)\n",
    "        nsm_Y[:,i_sample] = nsm.fit_next(x)\n",
    "        npca_Y[:,i_sample] = npca.fit_next(x_white)\n",
    "        \n",
    "        # if a neuron hasn't fired by the 100th iteration, we flip the feedforward weights:\n",
    "        \n",
    "        if i_sample==100:\n",
    "            for j in range(s_dim):\n",
    "                if bionica_indirect_Y[j,0:i_sample].sum()==0:\n",
    "                    bionica_indirect.flip_weights(j)\n",
    "                if bionica_direct_Y[j,0:i_sample].sum()==0:\n",
    "                    bionica_direct.flip_weights(j)\n",
    "                if nsm_Y[j,0:i_sample].sum()==0:\n",
    "                    nsm.flip_weights(j)\n",
    "                if npca_Y[j,0:i_sample].sum()==0:\n",
    "                    npca.flip_weights(j)\n",
    "        \n",
    "    print('Computing errors...')\n",
    "\n",
    "    bionica_indirect_err[i_trial,:] = permutation_error(S[:,idx],bionica_indirect_Y)\n",
    "    bionica_direct_err[i_trial,:] = permutation_error(S[:,idx],bionica_direct_Y)\n",
    "    nsm_err[i_trial,:] = permutation_error(S[:,idx],nsm_Y)\n",
    "    npca_err[i_trial,:] = permutation_error(S[:,idx],npca_Y)\n",
    "\n",
    "    print(f'Final Bio-NICA (indirect) permutation error: {bionica_indirect_err[i_trial,-1]}')\n",
    "    print(f'Final Bio-NICA (direct) permutation error: {bionica_direct_err[i_trial,-1]}')\n",
    "    print(f'Final 2-layer NSM permutation error: {nsm_err[i_trial,-1]}')\n",
    "    print(f'Final Nonnegative PCA permutation error: {npca_err[i_trial,-1]}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAAFQCAYAAADz1MZXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOydd3hUxfrHP7Ob3nujBQgdBOm9qAgoIIjYCyp4L6iIhZ8NRQXLVe+1Yi+goiBKUURAgdBROkiHEEgPISG97e78/jibZDd1E5JNYT7Psw97ZubMec8h+91333lnRkgpUSgUCkXdo6tvAxQKheJKQQmuQqFQ2AkluAqFQmEnlOAqFAqFnVCCq1AoFHZCCa5CoVDYCSW4CoQQC4UQ883vhwghTtS3TZYIIX4XQtxnh+tIIUREXV+nPhBCRAohpta3HVc6SnCbAEKI74QQCUKIDCHEycv5YEkpt0opO9SmfZeLlHKMlHJRfdtxpSCEiBZCXFffdjRFlOA2DV4HwqWUXsB4YL4Qolc926SoJkIIh/q2QVG3KMFtAkgpj0gp84sOza+2FbUXQlwthNgnhMgUQiwFXCzqhgshYi2Oo4UQs4UQh4QQ2UKIL4UQweaf+ZlCiD+FEL4W7fsLIXYIIS4JIQ4KIYZb1EUKIeYJIbabz10vhAgw17mYPfWL5nN3CyGCLc6ban6vE0LMEUKcE0IkCyG+EUJ4m+vCzWGB+4QQ54UQKUKI5y2u31cIsdPcf4IQ4kMhhJMtz1gIMUUIEWW2+6wQ4i4b7LF6lhbP8zrz+5eEED+Z7zsDmCKE8BNCfC2EiBdCpAkhVlqcO1YIccBs/w4hxFWV2DtSCHFcCJEuhPgQEBZ1bYUQG83POkUIsVgI4WOu+xZoCfwqhMgSQvyfuXyZECLR3N8WIUQXW56bohRSSvVqAi/gIyAHTWz3AR4VtHMCzgGPA47ALUAhMN9cPxyItWgfDewCgoFmQLK5/6sBZ2AjMNfcthlwEbgB7ct8pPk40FwfCZwB2gOu5uM3zHX/An4F3AA90Avwsjhvqvn9A8BpoA3gASwHvjXXhZvv/3Nz/92BfKCTub4X0B9wMLc9BsyyuFcJRJTzzNyBDKCD+TgU6GKDPVbP0uJ5Xmd+/5L52U8wPy9X4DdgKeBr/v8ZZm7b0/zs+5mfz33mvpzLsTfAbO8t5j4eBwwWzzDC/H/jDAQCW4B3y7PRouwBwNN8zrvAgfr+m2+Mr3o3QL1q8T9T+yAOBuYAjhW0GQrEA8KibAeVC+5dFsc/Ax9bHD8KrDS/f7pIbCzq1wH3md9HAnMs6mYAa83vHzDbcVU5NkdaiMUGYIZFXQezaBWJqASaW9T/DdxewbOYBaywOK5McC8BkwDXUnWV2WP1LC2ep6XgbrGoCwVMgG85NnwMzCtVdgKzIJcqvxfYZXEsgNiiZ1hO+wnA/vJsrKC9j/lZedf333xje6mQQhNCSmmUUm4DmgPToXiEP8v8ugsIA+Kk+ZNj5lwVXSdZvM8t59jD/L4VMNn8k/eSEOIS2hdAqEX7RIv3ORbnfosmzkvMP6ffFEI4lmNLWCl7z6GJW3BV1xBCtBdCrDb/NM4AXkPzBitFSpkN3Ab8G0gQQvwmhOhYDXsqI8bifQsgVUqZVk67VsCTpZ5tC/P1SxNm2a/5/7r4WAgRJIRYIoSIMz+H76jkOQgh9EKIN4QQZ8zto81VVT47hTVKcJsmDphjuFIb4fcwvxYDCUAzIYSwaN+ylq4bg+bh+li83KWUb1R1opSyUEr5spSyMzAQGIvmqZUmHk18imiJ9nM5qZy2pfkYOA60k9oA43NYxDarsG+dlHIk2pfHcbSwRVX2ZKOFSABNuNB+wlt1bfE+BvAriqeWIgZ4tdSzdZNS/lBO2wQ0MS66rrA8RhtklWi/JryAu7F+DqWXELwTuAm4DvBG+yUBNj47RQlKcBs5Zm/ldiGEh9kTGQXcgRZbLY+daIIwUwjhIIS4GehbS+Z8B4wTQowy2+JiHjhqbsN9jBBCdDOLUgbaz3JjOU1/AB4XQrQWQnigealLpZQGG+zzNPedZfZQp9tyU0IbJBwvhHBHiwlnWdhWmT0nARchxI1mb30OWgy0XKSUCcDvwEdCCF8hhKMQYqi5+nPg30KIfkLD3dyvZzld/QZ0EULcLLTMh5lASKnnkAVcEkI0A2aXOj8JLSZt2T4fLR7vZr5HRQ1Qgtv4kWjCEQukAW+jDQStKrexlAXAzcAUc/vb0AZ6Lt8QKWPQPKHngAtoXtlsbPs7CwF+QhPEY8BmNAEvzVdo4YctwFkgDy2ObAtPoXlrmWgCttTG83TAk2jebCowDC3+XKk9Usp0c7svgDg0j9cqa6Ec7kH7sjmONkg2y9zXHmAa8CHa/9tptP/DMkgpU4DJwBtoItkO2G7R5GW0Qbh0NHEu/f//OjDHHLp4CvgGLVQSBxxFG0RV1ABhHcpTKBQKRV2hPFyFQqGwE0pwFQqFwk4owVUoFAo7oQRXoVAo7ESTXCzDx8dHRkQ0yVX2GizZ2dm4u7vXtxlXFOqZ25+9e/emSClL51LbTJMU3ODgYPbs2VPfZlxRREZGMnz48Po244pCPXP7I4SoalZmpTSpkIIQYpwQ4rOsrKz6NkWhUCjK0KQEV0r5q5TyIQ8Pj6obKxQKhZ1pUoKrUCgUDRkluAqFQmEnmpTgqhiuQqFoyDQpwVUxXIVC0ZBpUoKrUCgUDRkluAqFQmEnlOAqFAqFnWhSgqsGzRQKRUOmSQmuGjRTKBQNmSYluAqFQtGQUYKrUCgUdkIJrkKhUNgJJbgKhUJhJ5qU4KosBYVC0ZBpUoKrshQUCkVDpkkJrkKhUDRklOAqFAqFnVCCq1AoFHZCCa5CoVDYCSW4CoVCYSealOCqtDCFQtGQaVKCq9LCFApFQ6ZJCa5CoVA0ZJTgKhQKhZ1QgqtQKBR2QgmuQqFQ2AkluAqFQmEnlOAqFAqFnVCCq1AoFHZCCa5CoVDYCSW4CoVCYSealOCqqb0KhaIh06QEV03tVSgUDZkmJbgKhULRkFGCq1AoFHZCCa5CoVDYCSW4CoVCYSeU4CoUCoWdUIKrUCgUdkIJrkKhUNgJJbgKhUJhJ5Tg1jJGk9HqvdFkxCRNHEg+QE5hTj1aplAo6huH+jagLkgzpLEldgt9Q/ri4uBil2tKKXln7zssOroIkzRV2M7PxY/O/p159OpH6ezf2aruROoJcg259AjqUdfmKhSKeqBJCm6mKZOHNzwMgLPemUFhg7ivy31cHXQ1Qogy7Y0mI7sSdhHoFkiETwQ6oTn+UkqWnVxGYnYiPYJ60D2wO97O3sV1nx76lKScJLIKslgbvdYm21LzUtkWtw1PJ0/eHPpmcfm09dPYlbCLT0d+WlxWYCzgYu5FQtxDyrVboVA0Lpqk4FqSb8xnY8xGNsZsBGBY82F4OnkS6h5Kj6AebDi/geWnlhe37+jXkdcHv06EbwSbYzczb9c8AHoG9eTlgS8XC25CdgILDiyosV1Dmg0pfl9gLCA2MxaAl3e8zENXPcSk9pN4a/dbLDmxBIC7O93NI1c/gruje42vqVAo6pcmKbiuOleaezQnNiu2TN3m2M2Vnts9sDvB7sEADG8xnLWT1vLZoc9Yfmo5uxJ2Ee4dDsChC4cq7OOBrg8gpaSTfyfyDHkkZCfQNaArXfy7sD1+Oz8c+4Hugd2L22+N3Vpsa3x2PKEeoQA83fdpJJKlJ5by3bHv+O7Yd/i5+HFNy2voF9qPgWED8XLyqtazUSgU9UeTFNxAh0B+n/Q7hy8cZsmJJayOWl1pXNWSpSeWsuH8Bj657hM6+HVgX9K+Yg944ZGFjGw1En9XfySScK9wojOiCXQN5FL+Jb4Z8w1dA7pW2v/4tuMZ33a8Vdm1ra5lyY1LiEqPYmSrkcVx5zxDHk56J7ydvUnPTwe0kMRPJ3/ip5M/AVpM+IleTzC69Wic9c7Vek4KhcK+CCllfdtQKUIId+AjoACIlFIuruqcDh06yBMnTliVpeSmsC9pH6ujVrMrYReuDq5kF2aTb8wn3CucRWMW4e3kzcELB8kx5DAobFCDiJvGZ8Uzf9d8tsZtrbRdsFswM3vOZEz4GBz1jnayroTIyEiGDx9u9+teyahnbn+EEHullL1rfH59CK4Q4itgLJAspexqUT4aeA/QA19IKd8QQtwDXJJS/iqEWCqlvK2q/ssT3PIwmAwIBHqdvsb3Yk+i0qPYHLOZs+lnWXF6Rblt/Fz8GN92PHd1uosQ9xC72aY+/PZHPXP7c7mCW18hhYXAh8A3RQVCCD2wABgJxAK7hRC/AM2Bw+ZmRmzAOT8Vsi6AR2Cl7Rx0jSui0sa7DW282wDwyqBXOJJyhPf3v8/fiX9jMBkALeSw8MhCFh5ZSM+gnoxtO5ZxbcbZLT1OoVBUTL2FFIQQ4cDqIg9XCDEAeElKOcp8/Ky5aSyQJqVcLYRYIqW8vYL+HgIeAugVquv19798SQoeSnzYKDK9OtTx3dQvuaZcIjMi2ZG1g0vGS2XqPXWezAmbg5verc5syMrKQu20YV/UM7c/I0aMaJQebnk0A2IsjmOBfsD7wIdCiBuBXys6WUr5GfAZQO8wvdTJQkITNxCauAFaDYbhT0ProbZbk5kErj6gM8dDdRaT8kxGMOSBU6kUrcI8uHgKMhKgIBM6jgUHZ5ASctPA1RfqKC48hjEYTUa2xW3ji8NfcODCgeK6Ea1HcMPgG+rkukWon7f2Rz3zxkdDEtzylEhKKbOB+6vTUT6lBo3ObYNF2zQBHPEcBHW2Fr6EQxD5hlbm0wp2lZNfO+RJuPZF7X3ORXi7nfbePwJSo6C8LIjnk7R/f34Q/vlZe99mOPR9CDrcUOviq9fpGdZiGMNaDCM5J5k1UWtYeXolM6+eWaZtSm4KAa4BtXp9hUJROQ1JcGOBFhbHzYH46nQghBgHjHMKiWBS/nRmOqxgmN4iX/b4aji+mriAgXzb+i1u7t2K9sGe5G3/CJcTv1XeuWdoyXvLMMzF0+W3H/MWZF8A7+bQ405wD4KjqyAqUnuVpt+/of908A237WarIMgtiCldp3Bfl/vKZFuk5aUxdsVYegf3Znr36XQJ6FIr11QoFJXTkBav2Q20E0K0FkI4AbcDv1Sng6Jde4PddLi0HciDxmcZlz+f3419rNo5X/iHddt2cf07Wwh/ZjXX7e1fZd/GzhNLDgps2Ib999maJywEtL1WE97MSr4//voEFvS3FvNaoLzUtq//+Zrswmw2x25m9pbZxQNuCoWibqkXD1cI8QMwHAgQQsQCc6WUXwohHgHWoaWFfSWlPFKT/l0dYPHU/py/mMPrvwcz/Z829DCc5t8OvzJav5sAkcFXjm8youAdQNBXHCs+N0H6ESpSOWwK53vjtfQQZ9hp6szxL44yZWA4Q9oHkpjlh/+Y72h1ZAEybj8CE8JYAJMXQpeJkHRUi+/6tiq6Yc3LdfWBuL2wfzEY88sa3muKdZhh/3ew6mG4aQF0vwNqKX0toyADgUAieaLXE40uW0OhaKw0+IkP1aEopBAWFjYtLi6uuDwzr5APN57m6x3RNDPG4Usm+2Q7SsLGRc9AO9ZhwlRD579LmBfdmnnTrbk3wzsE0czHtWyjvAzY8Ars/ty6fPDjMHR2yWDcnq9g9ePa+8BOMPwZ6DSuVoT3dNpp1p1bx4zuM8p4wQeSD9A9sHu1Jn6oARz7o565/WmUEx/qmoomPsRdyuXLrWeJv5TL9V2CCfF2YcvJFNJzCxnZOYhh7YMwSUlOgZECg4lPN5/hi21nL8uWiCAPBkcE4OXqyB19WxDqXY4AF2EoAAcn7b2U8OtjsG+RdRtHd2jeS/OkO40H99od+NoSu4WHNzxMr+BePNX7qSqnKhehPvz2Rz1z+6MEtxxsnWlmKwnpuSzbE8sHG09RaLy85+XiqOPBwa2ZPjwCD+cqfsqnRsH6F7RBtvLixkIHbgFauts1c8Cv9WXZlmvIZeKqicRllfw6GBM+hkd7PkoLzxaVnKk+/PWBeub2RwluOdS24FryT1w651NzaOnnxpLd5+kY4sXAtv6k5RSSnltAdEoOqw7GczCm7ASE0ky8uhmPXduO8IAqllzMSoZdH8PuL8G8iI0Vzl4w8wC4+9fwrsyXKcji/f3vs+zEMgyyZCDNUefIlC5TmNptKm6O5U+eUB9++6Oeuf1RgmtBRTHc+iI2LYcdZy7y895Y/jqbWm4bnYCh7QO5rXcLruscjKO+ktixIR9i/oK4fXByLZzfqZVf95IW/y1CStjxAeRdgh53gX/batkdnR7N+/vf549zf1iVB7sF81SfpxjValSZ+K768Nsf9cztjxLccqhLD/dyuJRTwG+HE/h0cxTnUyve36xdkAf/ueUqerb0rbzDtGgtk6FZb2g/qiTDwVgI73aDzATtuEU/uPoeuOq2khixDRxIPsCbu9/kcMphq/IBoQN4rt9zxWsDg/rw1wfqmdsfJbjl0FAFtwijSfLT3hgW7jjHsYSMCtt1DPGkfxt/Zo/qgHtV8V5LDi2D5VPLlnuGQv8Z0HcaOFYyeGeBSZpYdXoV7+57l9S8Ei/dUefIg90eZGq3qTjrndWHvx5Qz9z+KMEth4YuuEVIKTkSn8HafxJZtDOazLzyJyD4uzvx6DURTOzZHG9XG9a6LciBE2u06cSn1kPpiQ0uPlp6WY87oUV/63UiKiCjIIMF+xew5MQSq8XcW3m14oX+L5B7Ild9+O2MElz7owTXgoYWw60O2fkG9p+/RPylXFYfTmDbqQuYSv3XOOgEw9oHcnPP5lzbKQgXRxvycbNTYP+3Wkw352LZ+iFPwbUv2GznsYvHmL9rPodSrLcY6uPeh7fHvo2fi5/NfSkuDyW49kcJbjk0Fg+3MmJSc1i4I5ovK8gDdnHUMbpLCHf2a0Xf1jaIXEGOltO76yO4dL6kfPIi6DKh5Nho0FZCc6542T+jychPJ3/ivX3vkVmYWVzu6+zL8/2fZ1T4qKrtUVw2SnDtjxLccmgKgltEboGRpbvP8/O+OA7HlZMSBvQJ9+XBwW0Y2TkYva6K2WEmk5bpsP9bOP4bTN8B3s1K6hMPwyeDIbgbtBuphR0C2pXb1YWcC/xn939YF73OqnxGjxlM7z69WvepqD5KcO2PEtxyaEqCa8m5i9ms2B/HT3tjiU3LLVPfzMeV2/u04N6B4bbFeo0G0JcajDuyApZNKTkWOuh4ozbY1nJAuUtKbondwvORz3PJeAlnvTM/j/+ZVl6tqnl3iuqiBNf+KMEth6YquEVIKdkfc4nvdp7jl4PxGEoFez1dHJgyMJzx3cNoF+xZvc53fwG/PVl+XchV2hKSXSdpC6tb8PvG3/nL5S/aeLfh3i73Vu+aihqhBNf+KMG1oDEPmtWUxPQ8vtkZzeK/zpOeW1imvk+4L/cPas2oLiFVhxuKyE2Dk+vg0I9wZkPZelc/LdxwzQvgo035LfrwSynLTIpYenwpnk6ejGk9pkHshNxUUIJrf5TglkNT93DLI6/QyK8H41mw6TTRF8tOqugY4skzYzoyrH1g9UQv6ai2Vu+hpdpgmiWP7IWACKDiD39UehSTf5lMgamAa1tey7xB8/B0qqbXrSgXJbj253IFtyEtQK64DFwc9Uzu3YL1jw/j9Zu7MbxDIA4WHu3xxEymfL2bMe9t5dtd57iUU2Bbx8GdYfz78PhRbYshr+ZaubMX+LWxbmsogDMbrRZRX7B/AQUm7VoJ2Qm4Otg24UKhaIo0WQ/3+LFjCBsS+psyReGGr7dHk1tovcO8q6Oe+waG8+9hbfBxs326LyYjnNsOKSehT8lstsjISIb7JcPyadpU4yFPQLvryTYV8L89/2PF6RUsGbuE9r7ta+nuFMrDtT8qpFAOXV1c5c9du+L/r4cImDatyva5hw+T8uECsnfsoMPePQgnTYBkQQG5R46AyYRTq1bofX0R+trZdcGeJGXksWDTaZbtiS0jvB7ODkwd0poHB7fG08WGzIYKiNy0ieEnXoREiwkRPi1h+LPQ7VYSci8Q6hFqdY5Jmtget53BzQar2G4NUIJrf1RIoQJMWVm4dOxUfCwLCsg9cABTdrZVO2k0EvvIo2Rt3owsLLRKexJOTmRv38G5u+4mbtbjVueZ8vOJf/Y5Tg4eQvJ775W9fl4eDeXLLNjLhVdu6srOZ6/huRs60i6oZFJDVr6Bd/88xdA3N/Hp5jPkFhgr6alihDRAy/6gt8heuHQeVk6Hj/oTena7lgNswdITS5mxYQYzN83kQs6FGl1XoWhMNCnBFUKME0J8BiCcnXEfPKi4LvfwYc7dex8nevUm/dfVmAq0uKLQ6/G7fwoAzu0iMFy4gDRqomO8dInUhQu1848cwZRbMmh0ackSMn7/HWNKCqJ0LisQM306J3r24uwtk8k/dapMvSm/nD3N6hgfNyceGtqWdbOG8v4dVxNhIbxpOYW8/vtxhr61ic+3RJGZVzbjoTKkzhFueAseOwADHgFXi5XOLp7Stor/ZBAc/QWkJDYzlnf2vgNAZEwkN626iV/P/NpgvqQUirqgSQlu0a69huAg2v6x3upnqluvXjRf8CHC0ZH0FStIXbiIwqRkAHwnT9Ya6R3I2bu3OGyQf/o0zhERCEdHzUPev7+4P49hw/C85hoA8o4eLWOL4cIFZG4uef/8U2aygDQYONmvP8e6diP5nXcxZtmwC3AtotMJxncPY92sofx3cnda+JUMZF3IzOfVNcfo99oG5qw8TGxaxctIlotXGIx6FR47BMOf0wbXikg+Cj/eA19ej5+TNze1vam4KrMgk+e2PcejGx8lOSf5cm9RoWiQNCnBLUK6uOAYFGRdZjJRGBePcHUle8cOLn7xBUWbR+rc3el0/BhtVq7Ae9y44nPcevcmfMkPdNi3l47/HMZjyODiOmNWNhlr1gAQ8vxzZWwwZZWELvQ+PlZ1uYcPI/PywGAgfeXK4phxsa3Gmv2sry56nWBSr+ZseGI4r07sSoiXS3FdToGR73ad55q3NzN/9VHSc6rn8eLiBcOfhscOwpAntb3YimjWEzdnT57v/zxfjfqKZh4lU4s3x25mwqoJrDq9Snm7iiZHkxw0qygPVxqNZG3eQmFsLF7jxuLgW8UC31VQmJSM3scbnbNz+fUJCRQmJODavbvVYFvakqUkvvwySEnw88/jd8/dVufFTJ9B/okTuA8Zgu9dd+LS3j4j+3mFRpbvi+Pr7Wc5lWztdXs4O3DfwFZMG1J+VkOVAzjZKbD9XW2L+Bm7wDO4uCqnMId3/3qDH86ssDplSLMhzB0wl2D34NK9KVCDZvWBylIoh8Yw8aEwMZGLn31O4MxHrTxgU14eJ/sP0DxgnY7WK5bj0qGDXW2TUrLjzEXeXn+C/eet92bzdHbggcGteXBIa7wsshps/vAX5pZd/DwzCT7oxe6Q9rzoZiQ2v2Shc09HT/6v7/9xU9ubVCZDKZTg2h+VpdBIcQwJIeTFF8qEG/KOHkOaB9R8Jk2yEltTdjYX3n+f4z2u5linzhhSy+6TZkhJwZCWdlm2CSEYFBHA8ukD+eiunrQJLAkHZOYbeG/DKQa/sZEPN54iK7/8RdMrpLydJjbNh4JM+pzfy88nDnKXsaRNZmEmL2x/gRkbZpCYnVjTW1IoGgRKcBsYbj2vJmLDn/hNmYJzuwgKk5KK63Tu7sjCQs37lRJDYlkBSnz1VU4NGEjUxJutzq0JQghu6BbKn48P44M7rqathfBm5Bl4e/1JhvxnI59sPkO+oYa/lExGyM8CoYVc3KTkmfMn+DohiRamkj/PbXHbmLhqIitOrVCxXUWjRYUUGhFSSs5cex2F8fEA+E+bRtCTTxTXF8bHc3rk9WA0gqMjHS0mcYAWwzampeEQEFCj6xtNkl8OxvHen6fKrNfg5QSPXd+Zu/q1tG0nitJkJsG2d2DPV2DUPPxcIXjf15vF3l5Y/pUOChvESwNfIsQ9pEb30VRQIQX7o2K4FjT11cKklBScPYtwcsZ4MQWHkFAcg0uyMbK2bSf5P/8h/9QpXDp3pvXyn63Ozzt2jLO3TMa5TWtcOnch5JWXKxzwqwyD0cTy/XG8v+FUmXV5g72c+b9RHbm5Z7OaxVzT47TtgPYuBIPW9z5nZ14I9OO8Y0nM2N3Rndm9Z3Nzu5uv2NiuElz7owS3HJqqh2srhYmJGDMyymQ3pH7/PUmvzAPAY/hwWnzy8WVdp8Bg4qe9sXy48RTx6dYrifVu5cvLN3WhS5h3zTrPugAbX4F93wKSXCH4IKw13zkZkWZ/18/Fj18m/IK3cw2v0chRgmt/1KCZogyOISHlppIVWnj9AY8+UqY+c+NGq1l4VeHkoOPOfi3ZNHs4d3dyIsizxFvecy6NcR9s48VV/1Q/hxfAIxDGfwDTt0O7UbhKyf8NfZ1FYxYV7yYxp/+cK1ZsFY0T5eFeQUgpMcTHk382GvcB/a1yg6WURI0dR8GZM+h9fWnx6Se4XnWVzX1HRkbSZ8BgPth4mi+3RVFoLPm78nN34unRHZjcqwU6WxdBL835XdCiHwhBriGXP879wfg247TVyTqN17Z9F4KMggy8nLyq7q8JoDxc+6M8XIXNCCFwbNYMj8GDyqx6lr1tGwVnzgDaQj/ObdtWu393ZweeGdORtbOGMqRdycBcanYBT/98mJs/3sGh2EuV9FAJLfsXT5F2dXBlfNvxcHQlHF6mTRdeOJbdR5Yy6qdR/HjiR5XJoGiQKMFVAODSpQsBMx/FITQUr3Fj0bm7W9UbMzPJ2r7dJiFrG+jBNw/05ZO7e9LMpySn9kDMJW5asJ3nVhwmLdvGBdArY/v7xW9zzm/nxZ0vkVWYxbxd8/j44OXFpxWKuqBSwRUaLexljKL+cPDzI3DGDCL+/IOgJ8tuIpn2wxJiHpzK2QkTydxQzj5npRBCMLprKH8+MYxHRkTgpNf+1KSE7/86z4j/RrL4r3MYTZfhid79M/SbDjoHLuj1OJi/DDxNkluy8rRdiRWKBkSlgis1d2alnWxRNACEXo/e03rPMVNeHqmLFgGQf+IExvQMm/tzddLz1KgOrHt8KMM7BBaXX/13O+oAACAASURBVMop5PkV/zBhwXb2n6/hzDg3PxjzBkzfSavw4SyLT+T+Sxk8e/EiQX++DJ8OhXM7a9a3QlEH2BJS2CWE6FPnligaLDIvD69R1yNcXXEICcF77I3lNKrcU20d4M7XU/rw+b29rZaDPByXzsSPdjB72UHiL+VW0kMlBLaHu37C5fYfeAJfxmWZJ2UkH4GvR8PKh/n4wEcsOb4EkzRV3pdCUYdUmaUghDgKtAfOAdmAQHN+bR/CtjMqS6FuMKanU3A+BtduXa3KC5OSOH7rrTSfPgOfmyeWWW6yNHmFRj7ZfIaPI8+QbygRQCe9jvsGtuKRa9rh7VrD7X4K82Dnh7D1v1CoCe/hqyZwd9YBTNJEn5A+vDzwZVp4Nv5ImcpSsD91PvFBCNGqvHIp5bmaXrSuUYJrXxJfmUfa998D4DlyJM0/eL+KMzRiUnN4+dej/HnMes0HP3cnHr+uHbf3bYmjvobjupdiYN2zcHYLT/Qexx9xW4qrXB1cmdVzFrd3vB2daLzjxkpw7Y9dZpoJIboDQ8yHW6WUB2t6wbqkqU/tbYiYCgo4M2o0hoQEAFp8/hkeQ4ZUcZY120+n8Na6ExyIsU4ZaxPgzpyxnbim42Wsh5sRT767P58c/ISv//kaoyxZ3L23bydeGf5fWng1Tm9XCa79sYeH+xgwDVhuLpoIfCal/KCmF61rlIdrX0w5OeyZP5/QS+naNkaltxQqKKgyzCCl5LfDCby+5jhxpWK513QMYu64zrTyd6/gbNs4knKEOdvncPrS6eIyV3Q81n06d3R/qNF5u0pw7Y89Jj48CPSTUr4opXwR6I8mwAoFADo3N3Kuv54WHy0oI7b5Z85w+trruPTz8kpzeIUQjL0qjA1PDuPp0R3xdCnZmHPj8WRGvrOFd/44SV5hzbcf6hLQhaXDP2BaehZ6sy25mHjj4AIeWHo9MRfVl7SibrFFcAVg+VduNJcpFJUipSRx3nwMFy6Q8PzzJL36WpXnuDjqmT68LZFPDeeOvi2K998sMJh4b8Mprn9nCxuP13ydXyfv5sy8fQ2LdS2JsFgzYm9+EpN+ncTiP5/CZLLPnnKKKw9bBPcr4C8hxEtCiJeAXcCXdWqVoklguHCBwvPntQO9Hp9bJtl8rr+HM6/ffBUrZwziquYlC9ScT83hgYV7mLpoDzGp1dxRuIigTnS5dw1Lh73HQwWOJd6uELwRt44HFvUh5uSamvWtUFRCVTPNdMBfwP1AKpAG3C+lfNcOtikaOY5BQbT5bTX+06bi/8ADuHTsWO0+urfwYcWMQbw6satVqtifx5K47n+beX/DqRqHGZzaj+bRB/5mccQ9RFikp+3VFTJp+2wWr7hT5e0qahVbBs12SikH2MmeWkENmtmfmgzg5OzbR8aa3wl6fFaZtRvKIzW7gDfXHmfJ7hir8nB/N14a34XhHYIqOLNqCnIu8unv/+LLzOMYLeLQvYN788qgVxpk3q4aNLM/9hg0Wy+EmCSu1GX1FXWCKS+PhOeeJ+2774gafxO5h/+p8hw/dyfemHQVy2cMpGuzkiUYoy/mMOXr3fzr2z3EptUszODk5s+jk35i8dD/ESFcisv3JO1h0i+TWBO1RptNp9ZnUFwGtgjuE8AyIF8IkSGEyBRC2D6ZXqEoh4zffqMgOhoA46VLOPj72Xxuz5a+rHp4MPNu6oKXRTbDuiNamGHBptMUGGoWCujS5nqW3rWdaV3uR2/e2DLPkEeYRxicWg8f9Ydjq6ucyqxQlIctMdzRUkqdlNJJSuklpfSUUl4ZKzwr6gzvm28m7D9voPP2Jmj2bBzDwqp1vl4nuGdAOBufGs7kXs2Ly/MKTby17gQ3vr+Vv8+W3UbeFpz0Tszs/QSLb1hMhE8E93a+lx7+XeGPF+HiKVh6F3w9BmJ216h/xZWLiuEqaoWaxhMNqanofXwQOuvv/pz9+3Fu3Rq9j49N/ew9l8qclUc4lmD94+uWXs15ZkxHAjyqv1kmQIGxAInE+cIpTWTztf53ujjTwmCgefuxcO1c8K/+gu2Xi4rh2h8Vw1U0ahz8/MqIrSE1ldiHH+HM2HFkrF9vUz+9Wvnx6yODeHFsZ9ydSnaz+GlvLCPejmTh9rMYjNUPMzjpnXDWO0NIV5h5APrPIMXBmdlBAdzcLJQfYjZgWtAP1j0PeenV7l9xZVGdGG6BiuEq7EHSq69hTE3FmJJC0muvY8q1bdlGB72OBwa35s8nhzG6S0hxeWaegZd+Pcq4D7fzV9TFmhvm7g+jX+fVPhNI1+vJ1elY6O1FnjRoK5S91wN2LoDCGi4zqWjyVCm45pitTkrpqGK4CnvgdeONOARpKV6hr7yMztW1ijOsCfV25ZN7erHw/j6E+7sVlx9LyOC2z3bx+NIDJGfmVdJD5TzYayZtvbUQwiu6YNyKwnK5qbDuOdhU9Yw6xZVJlYJr3mbnbiHEC+bjFkKIvnVvmuJKxfOaEbRZ/Suhr87HY+jQMvXSYFtq1vAOQax7fCizR3XAxbHkT33F/jiufXszX2yNqlE2Q9eAriwdt5R3R7xLv/s3waQvwaelVql3JvXqO6rdp+LKwJaQwkfAAOBO83EWsKDOLCqFEKKNEOJLIcRP9rqmov7Re3nhM6nsVOCC6GjOjBpNxjrbYrvODnoeHhHBxieHc+NVocXlmfkG5v92jNHvbmHDsaRq7/LrrHfm2pbXajsJd7sFHtkDY95iQ+/bGL1+CkuPLy2ZpXYpBiLfgDwVibvSsUVw+0kpHwbyAKSUaUDla+2ZEUJ8JYRIFkL8U6p8tBDihBDitBDimcr6kFJGSSkftOV6iqaNNJmIf34OhXFxxD32GMnv2j7DPMzHlQV39uS7B/vRJqBkVltUSjYPLtrDvV/9zd5zNdxbDcDBmUvdb2VexmFyDbnM/2s+09ZPIy4rDna8D5Gvw3vdYedHYMiv+XUUjRpbBLdQCKEHJIAQIhCw9XfYQmC0ZYG5rwXAGKAzcIcQorMQopsQYnWpV83naiqaHIXx8RTGxmoHej2eI0dWu4/B7QJYO2soc27sZLUE5NZTKUz6eAf//nYv5y5m18i+1LxUvJ1LFtr5O/FvJq6cwNITy7QPTG6qtgvFB71h3zdgLKzRdRSNF1vycO8CbgN6AouAW4A5UsplNl1AiHBgtZSyq/l4APCSlHKU+fhZACnl61X085OU8pZK6h8CHgIIDAzs9eOPP9pinqKWyMrKwsPDo86vI3Jy8PzpJ4w+vmSPH3dZfWUUSFaeLmDTeQOWnwK9gGHNHbgpwglv5+plQxbKQtZcWsOGjA1Ii1575xuZn5xIM0PJQjt5zkHEtJhAQuh1mPTVzxO21zNXlDBixAi7bLHTEbgWbR3cDVLKYzZfoKzg3oI2e22q+fgetLDFIxWc7w+8CowEvqhKmEFNfKgP7J2EL6Uss9h57sGDFCYl4XX99dXq62RSJh9tOs3KA/FW5Z7ODjx2XTvuGxhe7b3VDl04xJztczibfra4zE048uSlTCanJlsvKO0eBAMeht4PgIvtCUBq4oP9scfEB6SUx6WUC6SUH1ZHbCugPJehQtWXUl6UUv5bStnWFrFVXBmUFltTXh7xTz9D3MzHiHviSQxptsdj2wd78u7tV7N8xkD6ti5Z06FoYG3Me1vZdCK5WgNrVwVexbJxy7i/6/3FW/fkyELmebswrXM/4t39SxpnJ8Ofc2HJnRX0pmgq1McmTrGA5Vp3zYH4CtpWCyHEOCHEZ1lZWbXRnaIRcfGzz4sXw8navBlp42QJS3q29GXpQ/35+v4+tAksGVg7nZzF/V/v5rZPd7H3nO3rMzjrnXmi1xN8M+YbWnu3Li7/KzeBiWGB/Nj/HqSnxRoSfdTYcFOnPgR3N9BOCNFaCOEE3A78UhsdSyl/lVI+pOJaVx5+992L94QJAAQ9/X/VXgynCCEEIzoEsfaxoTx/Qyc8nEsG1v6OTmXSxzuZsXgv/8TZPo23e2B3zdvtYuHtGnKYl7SZh7oOJH7UPGh3PXQab32ilLDlbUiNqtG9KBoedSq4QogfgJ1AByFErBDiQSmlAXgEWAccA36UUh6pSzsUTR+9tzdhb7xOq+8X4zN5cpl6Q0pKtfpzctAxbWgbNj41jLv7t8RBVxLCWHM4kbEfbGPqoj2cTrbt15Sz3pknemvebrhXeHH5rsS/ufns9yzrcyuy9K7BJ36HjfO0rIbFt8I/y8FQgKLxYkuWws3Af4AgtPirAGRDnN4rhBgHjAsLC5sWFxdX3+ZcUTTkARxDaipRY8fhPnAgwc8/h4Ovb7X7OH8xhzfXHWf1oQSrcr1OcGfflsy6rh3+Nq5IlmfI46MDH7Ho6CKrLXwGhA7g5YEvE+oRqnm3nw6FxEPWJ3u3hKFPwVW3Ebl9V4N95k2Vyx00s0VwTwPjamGwzG6oLAX705AFN/bxx8n8fS0ALl27Er7sxzKDbrZyODadT7ecKSO8Xi4OzLquPXf1b4mzg76Cs605kHyAF7a/QHRGdHGZv4s/ayetxUXvDGc2apMmoiLLnuzqx/mA4bSc/Bp4hZatV9QJ9shSSGpMYqtQWCILC9E5lUyMDJz5aI3FFqBbc28+vLMnqx8dzIA2JZkGGXkGXll9lGve3sySv8/btBRkj6AeLBu3jCldpiDMyTtTu03FxcFFmzIccS3cuwpm7odhT4ObRWZDbiotY5bD/zrC97dB1oUa35PCftji4b4HhAArgeI5iVLK5XVrWs1RHq79acgeLkDmxo3k7ttH0FNP1VqfUko2HEtm/m9Hib5ovZdamLcLd/ZryV39WuHrXvVM+APJB/j51M+8PPDl4oG1MuRnwp6v4O/PIT0GiUAgwdkbnjoBjtVbVU1RfewRUvi6nGIppXygphetK1QMt/5o6IJbEQXR0aR88QXBs2ej9/au+oRyyDcY+XbnOT6OPMPF7LKDWr1a+TKkXQDju4fRJrB6GTSn007z3r73eK7fc1psF7SNLE+uJf33+XhnmH98PnMeXCzsz0gAZw9w9qzRPSnKp84FtzGiPFz70xgFV5pMnLvnXnL37sUhMJCwt97CvX+/GveXnW9g4Y5ovtx2ltRyhBdgTNcQnhjZnnbBVQuhwWTgnjX38M/Ff3B3dOe1wa9xTctriusjIyMZ3sIErn4Q1sP65J+nwcl10PMe6HkfBLav8X0pSqjzGK4QorkQYoV51a8kIcTPQojmVZ2nUDR0srdtI3fvXgAMaWnovS7PG3R3duDhERHsevZa3rmtO92bl/WYf/8nkVHvbmHmD/vZeeZipbPX9ibt5chFLWOywFhAS8+WZRu1vaas2GYlw5EVkJ+u7USxoA98cxMc+1Vt817P2DJo9jXaxIQwoBnwq7lMoWjUeAwdSvMPP0AfEEDAQw/h0rlzrfTr5KBj4tXNWfXIYLY9PYJXJ3a12nnCJOGXg/Hc8fkuRrwdyceRZ7iQWXbJxn6h/Vg0ZhGtvFoxo8cMInwjbDMgPRZ8W1mXRUXC0rvh3W4Q+R/ITLyMO1TUFFtiuAeklD2qKmsIqBhu/dEYQwpFGNLS0Lu7I5ysB7fyz5zBISCgxrHd0hyOTeet9SfYcrJsRoGjXjDuqjAeviaCtqXivLmGXBx1jjjoHKzKP1jzAZOHTSbEPYQymIxwegPs/RpOrgVZKmtC6KBFPxjyFLS77rLv7UrBHmlhKeYtdvTm193AZezEV3eoqb2KmuDg61tGbE15ecQ+/AhR48aTtXlzrVynW3NvvnmgL6sfHcw9/VvhaTFtuNAoWb4/juv+t5mrX1nP4r/OFaeWuTq4lhHbmIwYFqYsZOKqiaw4taJsaEKnh/bXwx0/wGOHNGF1DyyplyY4vxPi9tTKvSlswxYPtyXwIdo2OxLYATwmpTxX9+bVDDVoZn8as4dbHklv/IfUhQsB0Lm7E7HhT/Q+PrV6jdwCI2sOJ/D93+cr3G1iUIQ/b97SnWY+JSlfUkoeXP8guxN3F5cNbjaYuQPmlu/tFmEogGO/aKll53YAEsZ/AD3vtW53ZiP4R5Ts06YoRmUplIMSXPvT1AQ3448/SHzpZYwXLxLyysv43nprnV5v3/k03vvzFJvLCTcA9An3ZVSXEG7p1RwfNyf2Ju1l9p+zuWAoae/p6MnsPrOZEDGh6skdmUlwdBW0GgAh3UrKC/Pgf5203SnCekLn8dD5JvBrUxu32eipM8EVQvyflPJNIcQHlLNerZRyZk0vWtcowbU/TU1wQYvtXlr2E/7TppYRsPIWQK8NDsVe4outZ/nlYPkrlro56bm7fyvuHdCKo/u3sd9jP4uPLbbaXWJIsyHMHTCXYPfg6htw4AdY+e+y5c16QZeJ2qpmgR2q328ToS4Fd5yU8lchxH3l1UspF9X0onWFGjSrP5qi4FaEITWV81OnEvT443gMGVIn1yg0mlhzOIHHlhwot14n4OogPU+M642TezRzd75ITGZMcX21vF1Lzm6B7e9D1CYwVZBCFtJNE98uE684z9ceM80ml96/rLyyhoTycO3PlSS4sbMeJ3OtthiO/7/+RdDjs+r0esmZeaz7J5Fvd53jZFLZ5SADPZ2Z3DuYbPdfWHFmae14u7lp2vKQR1dp2Q6mcja8bHsN3LOiurfTqLFHlsKzNpYpFE2ewqRkcv7+u/jYrVfPOr9mkKcL9wwIZ+1jQ/n83t4MivC3qr+Qmc9Hm86zeM3VdNM/Q6BLyeLrW+O2MnHVRFaeXlmtLYJw9YUed8KdS+HJEzDufeg4Fiw3u+xyc9nzVkyHt9vDoR+1OLHCispCCmOAG4BbgaUWVV5AZyll37o3r2YoD9f+XEkeriE1lcR589C5uRH26qv1YsOppEze+HkHB1N1pGSVmkYsCggJ30i2S6RV8dDmQ3mx/4s1i+0WkZeh5fUeWQkTFmjCXITJCPMCQZbsTIxPK80TbjsCWg0C94CaX7sBUJcx3O5AD+AV4EWLqkxgk5TS9l367IwSXPtzJQluEbKwEOHoaFVWEB1NQWwcHoMH1fn1IyMjGTRkKOuOJPLNznP8fdZ6vzW9axQuYT+hcyop93T05Om+TzO+7fjaH/SL2Q1fVjGJwqelNuGiRT+4+u5Gt8KZPWK4jlLKcgI4DRcluPbnShTc0kijUVsMZ98+fCbfQtDTT6Ovw0k4pZ/5sYQMvtl5jpX748gtNHuZogDnoLU4+e2wOrdP0CDeGPYKQW5BtWeQlBC9FfYu0tZy0DmAseyUZQAc3eCZGNBbTOjISQUHZ3ByL/+cBoA9BLcd8DrQGXApKpdSNrjhSZWlUH8owYXU7xaTNH++duDgQOtlP+LSqVOdXa+iZ56eU8iyvTF8u+sc58zr9OrdonAJtfZ2fTKncUObkfQJ92NQRACuTrbtVGEzhgKI/VubSHF2CyQcBKM5/BE+BKasLmlbmAu/PQUHvgMPc8jjmjlaGppnJZM57MzlCq5D1U34GpgLvAOMAO4Haj8BsRaQUv4K/NqhQ4dp9W2L4srDa8xocv7+m8z16wn417/qVGwrw9vNkalD2vDAoNasP5rI7J8O4Sg6kho1C+eg33Hy20lhendi4tvyaWwUn26JwsPZgXHdw7i1d3O6N/dBp6uFj7iDE4QP1l6gCXDSYTj/F3iWiiM7upas3ZtlHmz75VHzDbWE8EFaLnBAO2g12NozbkTY4uHulVL2EkIcllJ2M5dtlVLWTQJiLaBCCvZHebgaUkqyNkXiMXhQ2fUZsrPRudfez+XqPvPMvEI2HEtm0b4N7I9yAaO1LUKfiTR64O/uzLD2gYzoGMTQdoF4uzlW0GMt8/NUOGxDtqneCW76CK4quztzXWMPDzdPCKEDTgkhHgHi0HbwVSgUpRBC4HnNiDLlptxczk66Bbd+/QiaPRu9h/3jlJ4ujky4uhkTrr6XnAIDW06msO98Gn8eTSLqYipu4R9hyg8mNfFmlu8vYPl+LSzXNtCdQREB9G/jz4A2/jZtGVQjJn2hpZ/F7dFmvB3/TVvTtzTGAnC3To3jwglY8S8oyIaA9hDcFUKvgtDu4NVM2yOuAWCL4M4C3ICZwDy0sMK9lZ6hUCisuPDe+1oGQ3Q0uQcO0HrFcoTOljT4usHNyYHRXUMY3TWEZ8d0ZOYfLxCZkIbOKQ0350/JPvMEoMV0z1zI5syFbL7ZWbJeVYCHM9d2DOLOfi1pG+SBu5O+drIenNyg9VDtBVoYIvEQnNsOsbvh7FbIu6SJqtUNBWgxYmmClJNwfLV1vVdzCOqkxYe9m8GQJ7UBOjtji+CGSyl3A1lo8VuEEJOBv+rSMIWiqSBNJgrjS9ZG8Lvn7noV2/II9XYD887vM3pMp7B1RzYdT+ZgbDkeJpCSlc/SPTEs3VMynXhIuwAGtg1gRMdAOgR71o4AOzhB897aq4icVOv8X4CzkWXX/LUkI1Z7AegcYdgz1vU7P4KD32tTlX1aanFjZ09tXzivZuAbbr1nXA2xJYa7T0rZs6qyhoSK4dofFcOtHCklGb+tISsykrC33qwVMartZ/53wt/8ce4Pnuv3XLF9uQVGdkalsPPMRfacS+NwbDoGU9Uz1pz0OgZF+NOntR+juoTQJsC9Thb7sSI/E1JOwYXjkHgYEg5p3nF+hnU7n1Yw61DJscmoLVu5bErl/XuGIZ46XjcxXIuZZs2EEO9bVHkBamMkhaIaCCHwHnsj3mNvLFNnSE3l4udfEPDww/US2y2ib2hf+oZaTyB1ddLj759IcvxXfDZlDp6OfhxPzGTl/jgW7ojG2UGHlFBgtPYuC4wmNp24wKYTF3hzrbXz4+ygw0mvY0j7AJr7utE33I8OIZ4093W9PFF29oRmPbVXESYTpJ2F5GPaLheZCdri7Jbo9Nq2RFXh5ldz28xUFlKIB/YA44G9FuWZwOOXfeU6wCIPt75NUShsJvGVeWSuXUvm+vWEvfE6bn361LdJxeQZ8pizbQ7RGdHsSdrDs/2e5cbWN9KjhQ8vje8CaN771lMpHI5L52hCBr8fTqAyJzjfYCLfYGLNYW1ftc+2RAEQ4uVCiLcL/u5OCAEOOh13929Fz1Y+uDnVMA1MpwP/ttqrMrpOgqDO2gac6TGQEa8NwOVnaDHh9Dhw86+8Dxuo8C6klAeBg0KIxVLKRuHRqjxcRWMjZ8+e4pXHCuPiMOVVMDOrntget53ojGgAMgoyeHbrs6yLXsfcAXMJcNXWRRBCMLR9IEPbl2zhc+5iNpuOJ7PttBaOyC4wlte9FYkZeSRm5FmVrT2iiXKnUC+6N/emZ0tferbyoXWAB/rayBUuwitMe1WGlDDl8mLvtsRwz1L+AuQNbqZZESqGa39UDLdmSCnJWL2axPmv4jnyOsKKZqrZgL2e+V8Jf/Hi9heJzy4Z+PNy8ir2dqsKA0gp2R9zia0nU7i2UxDODjp+PRjPuiNJdAr1JDYtl2MJGTaJchEujjrCfFy5uoUv3q6OODoI2gV5EujpjIuDjrZBHgR41H4Wgj2m9lr60S7AZMBPSvliBafUO0pw7Y8S3MujMDkZnasrek9P6/KkJPQeHuVOmLDnM88uzOadve+w9MRSq/IRLUbw4oAXi73dmlJgMBGblkNSRj6xaTmk5xay/mhSmQV5qkOIlwsdQz3pEuaFq6OeQqPEz92JdkEeuDs7EOjpTJCnMw56273WetnTTAixTUo5uKYXrWuU4NofJbi1jzQaOXf3PRguXCD01Vdx72c9oFUfz3xXwi7mbp9r5e16O3vzbN9nuaH1DXWSiZCZV8jhuHQOxFxi37k0DsSkk5JVe6GXMG8XwnxcaRfsSadQTzqGeJGWU0CwlwvNfFwJ8HAqvq86n2kmhLBM/9IBvQHPCporFIpaIvWbb8ndvx+AmKlTafvnnzgG1+8kz/6h/Vl+03L+t+d//HjyRwDS89N5ZuszrI9ezwsDXrhsb7c0ni6ODGyr5fiCFqJIyshnw/EkjsRnEODhzJnkLNYfTcTVUU9GXvWGnOLT84hPz2NPBTsne7k4kG8wMaLD5T97W4b+/mvx3gBEoy1KrlAo6hCHgAB0Xl6YMjLwf+ihehfbItwd3XlhwAuMDB9p5e1ujNnI3uS9PNf3Oca0HlNnebdCCEK8XbirX6sK2xiMJs5cyOZATBrJGflkFRgwGCXJmfmcv5jN0YQMCo22/bovEvCiAbzLoUrBlVKWnRiuUCjqHO9xY3Hr25fUr74k4F8P1bc5ZSjydv+7578sO6ktOpOen87TW59m/bn1zOk/p9a9XVtx0OvoEOJJh5DKf4xn5BUSk5rD1lMpFBhMnEjM5LfDCXVmly2DZj5oayeEYyHQapt0hSUqhmtfTLm5HLr7Hto/83SDyNvdGb+TuTvmkpBdIlbezt517u3WNVJK4tPz2HcujSd/PMip126o800k16CJ7WG0CRBFL4VCUU9cePc9nI8c4dy995H87rv1bQ4Dwgaw4qYVTG5fsmRikbf71OanMFW2zkEDRghBMx9XxnUP4+SrYy67P1tiuC5Syicu+0oKhaJWMKank75ypXYgJU7Nm9evQWbcHd15ccCLjGw10srbbebZDJ1oWIv11Be2PIVvhRDThBChQgi/oledW1YDhBDjhBCfZWVl1bcpCkWdoff2pvUvq8jv2hX3IUPwnjSpvk2yYkDYAJaPX84t7W+hjXcbHu7xcH2b1GCwJYb7MPAqcImSGWdSzTRTWKJiuPYnctMmhvbtW2ZShOHiRQqio3Hr1aueLCshpzAHN0c3q7LknGT2J+9nVPioerKq5lxuHq4tHu4TQISUMlxK2dr8arBiq1BcMQhR7gy0xFfmce7ue0h6/XVMubn1YFgJpcVWSskrO1/hqc1P8UTkE1zMvVhPltUPtgjuESCnrg1RKBSXT8batWSuWwdSkrroG3IPHKhvk6xYG72WzbGbAfjj3B9EpUfVs0X2xRbBNQIHhBCfCiHeL3rV3IyqYwAAIABJREFUtWEKhaL6uHbvjvtgbda99y2TcB8woJ4tsmZws8FMaqfFnO/oeAd9Quo/pc2e2JKlsNL8UigUDRzH0FBafP4ZGb/8gsc115Spl4WFCEc77cJbDp5Onrw08CVGtx7NVQFXlalPyk4iyC2o0ebtVkWlgiuE0AP3SCmvs5M9CoXiMhFC4H3TTWXKpdHIufvvx7VrNwJnPYbOxaUerNPoH9q/TFlaXhq3rr6VnkE9eb7/8/U2S60uqTSkIKU0AjlCiMvfPU2hUNQrqd98S+6evaQuXEj05FuRBQX1bZIVr/31Gql5qfx5/k8mrJrAmqg11GQ1w4aMLSGFPOCwEOIPILuosCFP7VUoFNZIKcnZs6f42PP66xFOTvVokTUGkwF3x5KMi4ayJkNtY8ug2W/AC8AW1NRehaJRIoSg+YcfEPLKy7j17t3gFsNx0Dnw0sCX+PS6TwlxDyku33B+Q5PydqsUXCnlIuBHYJeUclHRq+5NUygUtYkQAt9bb6Xlt9+U8W5NubmkfPoZpvz63VNtYLOBrBi/glva31JcVuTtzto0i5TclHq07vKpUnDNO+EeANaaj3sIIX6pa8MUCkXdUF4GwIV33+XCO+9wduLN9Z676+HkwdwBc/l05KeEuocWl2+M2ciEVRP4Leq3Ruvt2hJSeAnoiza1FynlAaB1HdqkUCjsSN6xY6R+8y0ABVFR5J8+Xc8WaQwMG1i8JkMRRbtLNFZv1xbBNUgp00uVNc6vF4VCUQbnjh0JmTsXnZtbg1sMp8jb/WzkZ03C27VFcP8RQtwJ6IUQ7YQQHwA76tguhUJhJ4QQ+N5+G61/+YXQ+fPKhByMGRn1HtutaL3dZ7Y+w2ObHms03q4tgvso0AXIB74H0oFZdWmUQqGwP07Nm+EYHFymPOGFFzk7aRK5hw/Xg1UlFK23+9nIzwhzDysu3xSziZtW3sTqqNUN3tutUHCFEC5CiFnAm8B5YICUso+Uco6UMs9eBgohJgghPhdCrBJCXG+v6yoUCsj4/ff/b+/Ow6qq9oCPfxeDAgIOSF5n9KpMMouCI2iihbMWqTdN85qaenuse9HUMstrvdlrUmZZGqamds3XodTMAkdKIc0ZHEJFjckJB2Ra7x8HdsyI4DmHw/o8D09nD2fv31kdfyzW3vu3yPjhB7LOnSfxuZFkXbpk6JB09XYHb+LZDn/NZXs76zaz9s3SCuMYq/J6uKvQTYl+HHgKWFTZgwshVgohUoQQJ4qt7y+EiBdCnBNCzCzvGFLKzVLKfwIvAGGVjUFRlEeXl/kAYaMrsdhg2FDqtGpl4Ih0CmYO/jzkc6236/83f3q26GngyMpX3pNmblJKDwAhxArg0CMcPxL4GPiqYEV+fYalQF8gCTicf5uZObCw2PvHSylT8l/PyX+foih60mDoEGz8O5G6JIInwsMNHU4JBTMHR/wWwT/c/lFiKh8ppVEVwilzxgchxG9SSt+ylh/6BEI4Ad9JKTvmLwcC86SU/fKXZwFIKYsn24L3C+Bd4Ecp5e5yzjMRmAjg6Ojo980331Q2VKUK7ty5g62traHDqFUM3uZ5edivWcu9Xj3Jad3acHGUIVfmsixlGZ3rdca/nn+1JN7g4OAqzfhQXg/XSwhxO/+1AKzzlwW6KXbsH/GczYHLhZaTgC7l7D8NeBKoL4RoJ6X8tLSdpJTLgeWgm2JHTfeiX2qKHf0zdJunr/ySlIMHsf71VxpPmoTjtKkGi6U0y48tJ/5SPPGZ8aTYpfB+z/cN3tstM+FKKc0f0zlL+8RlXlqUUkYAD1XwPP+puIHNmjWrcF9FUR5d3r17pH/+uW4hNxfMjOfPdoDs3Gy2nd+mLbs5uBk82cLD3RZW3ZKAloWWWwBXq+PAUsptUsqJ6k9bRXm8zGxscPpmAzb+/tR1daXxROMqhmNpbsmGARt4zvk5PB09GeM2xtAhAQ9XnrG6HQbaCyHaAFeA54BRBohDUZQqqNOyJa1WRZJ740aJWSTy7t8n648/sHJzM1B0ugksZwfMJis3Cwuzoqnu4u2L/J76OwPbDtRrz/ex9nCFEOuAGMBZCJEkhHhRSpkDTAV+AE4D30gpTz7OOBRFeTyEmRkWDg4l1qcsXswfzzxLypIlBi90Xse8aGW03Lxc5h6Yy+z9s5n28zRS7qWU8c7q91gTrpRypJSyqZTSUkrZQkq5In/9dillBynl36WUC6rrfEKIgUKI5Xfu3KmuQyqKUkn3YmO5sXoN5OaSvuxTMn76ydAhFfHt2W85knIEgD1JexiyZQhbz2/Vy1NqhhjDfWzUGK6iGJ6FoyPWfro7SOv17IFd//4GjqioAW0HMNJlpLackZWht96uSSVcRVEMr07r1rT+6iuazJlD0/nzS4yRGrregY2lDa93eZ2V/VbSwraFtl4fvV2TSrhqSEFRjIMwM6PRP0Zj+be/FVkvpeTqq6+RGvGRwcd2/f/mz7eDvi21tzv156mPpbdrUglXDSkoinHL2LGD29u3k/bJJ/zxbBh59+8bNJ6yert7k/YyZMsQtpzbUq29XZNKuIqiGLdb27drr609PDCztjZgNH8p6O2OcvnrDtWMrAzmHJjD1J+nknw3uVrOoxKuoih602LJEprMmkmdNm14Ivw/hg6nCBtLG2Z1mVVqb3folqFsPre5yucwqYSrxnAVxbgJc3MajR1L261bMC829Cdzc7nxzTfI7GwDRadT0Nsd7TpaW5eRncHcA3OrfGyTSrhqDFdRaobiT6YBXF/1FX++8SZ/PBtG5pkzBojqLzaWNszsPJMv+31JSztdJYLmts2rfFyTSriKotRMWZcukbpkCQAPTp82moclOv2tExsHbmS062je7vZ2lY+nEq6iKAZn2bw5jv/6F6JuXaMrhlPQ2/X/m3+Vj2WI4jWKoihFCHNzHMaPwzYoCPJySxbDefAAYWZW6lBETWJSPVx10UxRara6bdtQt127EutTPviAxLDnyIxPMEBU1cekEq66aKYopqegGE7mqVP8MWIEmfHxhg7pkZlUwlUUxfRkXbykDSXUCwygbocOBo7o0amEqyiKUWswfBht/t8m6nXrVmoxnJpEJVxFUYxe3b//nVYrvii1GE7KkiU8OHvWQJFVjkq4iqLUWBk7dpC+7FP+GDactM+WG7z0Y0VMKuGquxQUpfaQWVkkv/d/dK+zs8lOSjL64QZh7L8RHoWzs7OML3YlMzs7m6SkJDIzMw0UlWnLzMzEysrK0GHUKuW1uZWVFS1atMCyht+3WpEH585xddbr5KSn0Xbr1hL1GaqbECJOStnpUd9fax58SEpKws7ODicnJ6P/LVgTZWRkYGdnZ+gwapWy2lxKSXp6OklJSbRp08YAkelP3XbtcFr3NdlXr5ZaDCf78mXqODkZJrhSmNSQQnkyMzNxcHBQyVYxeUIIHBwcas1fc8LCgjqtWpVYfz1yFRcGDiLt88+ROTkGiKykWpNwAZVslVqjtn/XH1y4QOqSJcjsbFI/+L/c/HaToUMCalnCVRSl9qjbvr3uv26uNBg21MDR6KiEq0fm5uZ4e3vj5eWFr68vBw8eBODq1auMGDGiUseaN28eNjY2pKT8NdFd4UeaC79OSEjg6aefpl27dri6uvLss8+SnPzXlCH/+te/aN68OXl5eWWe78iRI0yYMAGArVu38u6775YbX2JiIl9//XWlPlNN8dxzz3G2htz3WVvVbdsWpw3rcXzlFZotXGg8RW+klCb306FDB1ncqVOnSqzTt3r16mmvd+7cKXv27PnIx3rzzTdly5Yt5X/+859Sj1/w+v79+7Jdu3Zy69at2raff/5ZHj9+XEopZW5urmzZsqXs0qWLjIqKKvN8I0aMkEePHi1z++3bt4ssR0VFydDQ0Ep9ppycnErtXxnZ2dnVdqzo6Gg5YcKEajveoyre5sUZw3feWCUvWiTTvvhC5lXyOwfEyirkJpPq4dak+3Bv375Nw4YNAV1vsGPHjoDu4t64cePw8PDAx8eHqKioMo8xfvx4NmzYwPXr18vc5+uvvyYwMJCBAwdq64KDg7XzRUVF0bFjRyZPnsy6detKPUZGRgbHjh3Dy8sLgMjISKZOnQrACy+8wPTp03nyySdp27YtGzduBGDmzJns27cPb29vFi9eTG5uLv/+97/x9/fH09OTzz77DIDo6GiCg4MZNWoUHh4eREdHExQUxIgRI3BxcWH06NHazexxcXH06tULPz8/+vXrx7Vr1wAICgoiNjYWgLS0NJzyr0pHRkbyzDPPMHDgQEJCQrh+/TpDhgzB09OTgIAAjh07Buj+Whg/fjxBQUG0bduWiIgIAO7evUtoaCheXl507NiRDRs2ANCjRw92795NjpFciFEq597hw6R/sYKU9xdxcdRocm/f1tu5Teq2MCnlNmCbs7PzP8vbz2nm948thsR3Q8vcdv/+fby9vcnMzOTatWv8/PPPJfZZunQpAMePH+fMmTOEhISQkJBQ6v2Wtra2jB8/niVLlvDWW2+Ves4TJ07g5+dXZkzr1q1j5MiRDB48mNdff53s7OwS927GxsZqCbo0165dY9euXVy5coVBgwYxYsQI3n33XRYtWsR3330HwPLly6lfvz6HDx/mwYMHdOvWjZCQEAAOHTrEiRMnaNOmDdHR0Rw5coSTJ0/SrFkzunXrxoEDB+jSpQvTpk1jy5YtODo6smHDBmbPns3KlSvLjAsgJiaGY8eO0ahRI6ZNm4aPjw+bN2/m559/ZsyYMRw9ehSAM2fOEBUVRUZGBs7OzkyePJmdO3fSrFkzvv9e9325desWAGZmZrRr147ff/+93LZVjNP1NWsh/5e4WYP6mOnxdkaT6uEaO2tra44ePcqZM2fYuXMnY8aMKfEo4v79+3n++ecBcHFxoXXr1iQklF0DdPr06axatYrbj/BbOisri+3btzNkyBDs7e3p0qULu3btKrHftWvXcHR0LPM4Q4YMwczMDDc3tyJjw4Xt2rWLr776Cm9vb7p06UJ6ero2Dtq5c+ci94t27tyZFi1aYGZmhre3N4mJicTHx3PixAn69u2Lt7c377zzDklJSRV+xr59+9KoUSOgaNv27t2b9PR0LYmGhoZSt25dGjduzBNPPEFycjIeHh7s3r2b8PBw9u3bR/369bXjPvHEE1y9erXC8yvGp/mi93H813TMGzXSezEck+rh1iSBgYGkpaWRmppaZH3xBFxg9uzZWk+roFcG0KBBA0aNGsUnn3xS6vvc3d3Zs2dPqdt27tzJrVu38PDwAODevXvY2NgQGlq0l25tbV3uPZ1169atMH4pJR999BH9+vUrsj46Opp69eqVeTxzc3NycnKQUuLu7k5MTEyJY1tYWGgX/IrHWfjYpcVW8I+ttHN26NCBuLg4tm/fzqxZswgJCeGNN97QzmNtbV3qZ1WMm7C0pPHkyTQaMwazYt89KSUZu37Erk9vhEX1p8damXDL+7NfX86cOUNubi4ODg7cu3dPW9+zZ0/Wrl1L7969SUhI4NKlSzg7O7NgwQIWLFhQ6rFmzJiBv79/qWOKo0aNYuHChXz//fdaIt25cyfNmzdn3bp1fPHFF4wcORLQjVm2adNGS7wFXF1d+eCDDyr1+ezs7MjIyNCW+/Xrx7Jly+jduzeWlpYkJCTQvPnDz4Lq7OxMamoqMTExBAYGkp2dTUJCAu7u7jg5OREXF0fnzp21MeTSFLTt3LlziY6OpnHjxtjb25e5/9WrV2nUqBH/+Mc/sLW1JTIyUttWcG6l5iqebAFub9/O1Vdfw8rTk2b/XVDq7BNVOme1Hk0pV8EYrre3N2FhYaxatQpzc/Mi+0yZMoXc3Fw8PDwICwsjMjKySO+rNI0bN2bo0KE8ePCgxDZra2u+++47PvroI9q3b4+bmxuRkZHY29vzww8/FOnN1qtXj+7du7Nt27Yix3BxceHWrVtFEmhFPD09sbCwwMvLi8WLFzNhwgTc3Nzw9fWlY8eOvPTSS5W66FSnTh02btxIeHg4Xl5eeHt7a7fVvfbaayxbtoyuXbuSlpZW5jHmzZtHbGwsnp6ezJw5k1WrVpV7zuPHj9O5c2e8vb1ZsGABc+bMASA5ORlra2uaNm360PErxi8nPZ3kt98BIPPYMa6vXVvt56g1xWtOnz6Nq6urgSKq+RYvXoydnZ12L25xtamWwuLFi7G3t+fFF180aBwVtbn6zleOzM4mfcVKUpcuxcKxcanFcFTxGkUvJk+ezP/+9z9Dh2EUGjRooF18U0yHsLSk8aSXsO0dTN6dO6UWw6kqlXCVh2JlZaWSTL5x48YZOgTlMbIqY86064XG8B+VGsNVFEWpwIMLF7j9Q8lbJitLJVxFUZQKZJ48RdN5b1b5OCY1pCCEGAgMbNasmaFDURTFhNQfOKBajmNSPVwp5TYp5UTbxzzNhqIoyqMwqYRr7EylPOPatWvx9PTE09OTrl278vvvv1cqdkP78MMPizxsUlOlpaXRv39/Q4ehVEZVSo0Z648qz/h4yzMeOHBAXr9+XUop5fbt22Xnzp0rLBVoLHJycmTr1q1lamqqwWKorlKRt2/fli+88ILcv39/qduN4TtvalDlGWummlyesWvXrlrsAQEBZRaRuXPnjvZZPD09+fbbbwFdhTIPDw86duxIeHi4tr+trS3h4eH4+fnx5JNPcujQIa1k4tatWwFdycXBgwfTv39/nJ2di1RJGzJkCH5+fri7u7N8+fIix33jjTfo0qULCxYs4OrVqwQHBxMcHFxmu8FfpSe7du1apPQkwPvvv6+VmnzzTd3FlML/HwEWLVrEvHnzAF0Jyddff51evXqxZMkSfvrpJ3x8fPDw8GD8+PHaU4JOTk68+eab+Pr64uHhwZkzZwDYs2eP9pSij4+P9tTfkCFDWPsYnohSHg+TumhWKVELYU/5sxZofMfCoIii67ZOh98KPRraayYEzyr3MKZYnnHFihU89dRTpW57++23qV+/PsePHwfgxo0bXL16lfDwcOLi4mjYsCEhISFs3ryZIUOGcPfuXYKCgnjvvfcYOnQoc+bM4ccff+TUqVOMHTuWQYMGAX+Vc7SxscHf35/Q0FA6derEypUradSoEffv38ff35/hw4fj4ODA3bt36dixI/Pnzwdg5cqVREVF0bhx4zLbpcC1a9fYv38/Z86c0UpP7tq1i7Nnz3Lo0CGklAwaNIi9e/fSqpSJDAu7efMme/bsITMzk/bt2/PTTz/RoUMHxowZw7Jly3jllVcA3aPav/32G5988gmLFi3iiy++YNGiRSxdupRu3bpx584drKysuH//Pp06ddIeOVaMn+rh6pGplWeMiopixYoVvPfee6Uef/fu3bz88svacsOGDTl8+DBBQUE4OjpiYWHB6NGj2bt3L6Crl1AwJunh4UGvXr2wtLTEw8ODxMRE7Th9+/bFwcEBa2trhg0bxv79+wGIiIjAy8uLgIAALl++rJV/NDc3Z/jw4ZVuHyi99OSuXbvYtWsXPj4++Pr6cubMmYeacicsLAyA+Ph42rRpQ4f8G+zHjh2rtQHAsGHDAPDz89M+d7du3ZgxYwYRERHcvHkTi/xKVqpMZM1Se3u4BlbTyzMeO3aMCRMmsGPHDhwcHMjIyGDp0qV8/vnnAGzfvh0pZYlao2V9PgBLS0ttfzMzM61oj5mZWZFCN8WPKYQgOjqa3bt3ExMTg42NDUFBQVrMVlZWJYoEPazSSk9KKZk1axYvvfRSkX2TkpKKXHgsq1RkeW1Q+JwFZSJBN4NGaGgo27dvJyAggN27d9O8eXNVJrKGqb093OBZMO/Ww/0UH04A3brC+1QwnFBc4fKMhRWUEARKlGc8evRokWRbYMaMGXz22Wdllmc8ePCglqxBl2iPHz+ulWdMTEwkMTGRP/74g127dpW4gu/q6sq5c+e05UuXLjFs2DBWr16t9dIAXn75ZS3GZs2aERISwscff6xtv3HjBl26dGHPnj2kpaWRm5vLunXr6NWrV6Xa7scff+T69evcv3+fzZs3061bN27dukXDhg2xsbHhzJkz/PLLL2W+v3jpyMrq168fK1eupGAqpytXrpCSkkKTJk1ISUkhPT2dBw8eaLNdFOfi4kJiYqLWpqtXr66wDc6fP4+Hhwfh4eF06tRJG9tNSEgodzYOxbjU3oRrAKZSnnH+/Pmkp6czZcoUvL296dSp9OJJc+bM4caNG3Ts2BEvLy+ioqJo2rQpCxcuJDg4WLs9bvDgwQ/VfgW6d+/O888/j7e3N8OHD6dTp07079+fnJwcPD09mTt3LgEBAWW+f+LEiTz11FPaRbMJEyZoc6I9jJCQEEaNGkVgYCAeHh6MGDGCjIwMLC0ttYtzAwYMwMXFpdT3W1lZ8eWXX/LMM8/g4eGBmZkZkyZNKvecH374odaO1tbW2rh5VFRUib9IFOOlyjMqD8VYyjNGRkYSGxtbpOdcW2VkZBAaGsqWLVu0u0YKU9/56lfV8oyqh6s8lMmTJ1fY01b0Ky0tjRkzZpSabBXjpHq4SrWoTQXIjYUqQK5/qoerKIpSQ6iEqyiKoicq4SqKouiJ0SdcIYSrEOJTIcRGIcRkQ8ejKIryqB5rwhVCrBRCpAghThRb318IES+EOCeEmFneMaSUp6WUk4BngUcerDa0y5cvExwcjKurK+7u7ixZsqTMffVdzzcoKKjIvbSxsbEEBQUBuqfPRo8erRWb6d69u3bDvxCiyDxnOTk5ODo6MmBA9RRrVhRT87gf7Y0EPga+KlghhDAHlgJ9gSTgsBBiK2AOLCz2/vFSyhQhxCBgZv6xaiQLCws++OADfH19ycjIwM/Pj759++Lm5qbXOLQycWZFf9empKSwY8eOEoVolixZQpMmTbQCNPHx8Vpxm3r16nHixAnu378P6J4Aa968uR4+haLUTI814Uop9wohnIqt7gyck1JeABBCrAcGSykXAqV2jaSUW4GtQojvga9L20cIMRGYCODo6Eh0dHSR7fXr16/S45xVZWtrS/v27bUY2rdvz9mzZ2nZsmWp+2dkZHDnzh1GjhzJzZs3yc7OZu7cuYSGhvL222/j4ODAlClTAN2TX46OjkyePJklS5awadMmsrKyGDBgALNnz+bixYsMHz6cHj16cPjwYb7++usila1yc3OZNm0ab731Ft27d+fu3bvk5uaSkZHBxYsXadWqlRZ3s2bNyMrKIisrC4A+ffqwceNGBg4cyFdffcWwYcM4ePCgQdu6tij4f1SWzMzMEv8OFMMyRPGa5sDlQstJQJeydhZCBAHDgLrA9rL2k1IuB5aD7j7cgj+JC5w+fbrIPYufHP2EZb8vq3TwpenVohcf93n4zndiYiLHjx8nODi4zPso7ezssLa2ZuvWrdjb25OWlkZAQABhYWFMmTKFYcOGER4eTl5eHps2beLQoUPExMRw6dIl4uLitLKBR44coVWrVpw9e5ZVq1bxxRdflDiXubk5wcHB7Nixg9jYWOzs7DA3N8fOzo5JkyYREhLCtm3b6NOnD2PHjqV9+/bae8eMGcP8+fPp378/p0+f5qWXXuLQoUPqnlw9qOg+XCsrK3x8fPQYkVIRQyRcUcq6Mp++kFJGA9GPKxh9u3PnDsOHD+fDDz/E3t6+3H2llLz++uvs3bsXMzMzrly5QnJyMk5OTjg4OHDkyBGSk5Px8fHBwcGhSNnAgnOdPXuWVq1a0bp163LrC4Cu9sE777xTpNyit7c3Fy5cYNeuXezevRt/f39iYmK0G+o9PT1JTExk48aNPP3001VsHUUxbYZIuElA4b+jWwDVUtDT2Gftzc7OZvjw4YwePVqreXr58mVtNoZJkyYVKWKydu1aUlNTiYuLw9LSEicnJ63k34QJE4iMjOTPP/9k/PjxQNllAxMTE7XSgOXp3bs3c+fOLVFpy9bWlmHDhjFs2DDMzMzYvn17kSeYBg0axOzZs9mzZw/p6emP0DKKUjsYIuEeBtoLIdoAV4DngFHVcWAp5TZgm7Oz8z8r2neK9xSmeE+pjtM+FCklL774Iq6ursyYMUNb37Jly1JLLgLcunWLJ554AktLS6Kiorh48aK2bejQobzxxhtkZ2fz9de6Ye1+/foxd+5cRo8eja2tLVeuXCkxe0NFZs+ezaRJk2jbti0ABw4cwM3NjYYNG5KVlcWpU6coPlwzfvx46tati4eHhxozVJRyPNaEK4RYBwQBjYUQScCbUsoVQoipwA/o7kxYKaU8+TjjMAYHDhxg9erVeHh44O3tDcB///vfcv8MHz16NAMHDqRTp054e3sXKfdXp04dgoODadCggVbiMSQkhNOnTxMYGAjoeqZr1qypVPHtp59+usjsDufPn2fy5MlIKcnLyyM0NLTE7AktWrTQLuApilI2kypeU2hI4Z9Xrlwpss3UCnnk5eXh6+vL//73vyIXsQxFFa/RP1W8Rv9U8ZpCpJTbpJQT9f3ggL6dOnWKdu3a0adPH6NItoqiPBw1p1kN5ObmxoULFwwdhqIolWRSPVxFURRjZlIJVwgxUAixvOBZf0VRFGNiUgm3tozhKopSM5lUwlUURTFmKuHqkRCCV199VVtetGgR8+bNM1xAFUhMTNQeqgBd2cbp06dX+bjR0dHUr18fHx8fXF1deeutt7Rthw4domfPnjg7O+Pi4sKECRO4d++etn3w4MHafcaKUtOYVMI19jHcunXrsmnTJtLS0gwdykMpnnA7depEREREtRy7R48eHDlyhNjYWNasWUNcXBzJyck888wzvPfee8THx3P69Gn69++vVcS6efMmv/32Gzdv3uSPP/6oljgURZ9MKuEa+xiuhYUFEydOZPHixSW2Xbx4kT59+uDp6UmfPn24dOkSAC+88ALTp0+na9eutG3blo0bNwK6XmJQUBAjRozAxcWF0aNHU/AQS1xcHL169cLPz49+/fpx7do1AA7WS1HgAAALCklEQVQfPoynpyeBgYH8+9//pmPHjoAusfbo0QNfX198fX05ePAgADNnzmTfvn14e3uzePFioqOjGTBgAHl5eTg5OXHz5k0tfi8vL5KTk0lNTWX48OH4+/vj7+/PgQMHym2TevXq4efnx/nz51m6dCljx47VerBCCEaMGEGTJk0A+Pbbbxk4cCDPPfcc69evf+T/D4piMAUFqU3pp0OHDrK4U6dOFVlOifhInnJ2kaecXWRKxEcl9v9z4bva9rQVK0tsvzpnrrb9+voNJbaXpl69evLWrVuydevW8ubNm/L999+Xb775ppRSygEDBsjIyEgppZQrVqyQgwcPllJKOXbsWDlixAiZm5srT548Kf/+979LKaWMioqS9vb28vLlyzI3N1cGBATIffv2yaysLBkYGChTUlKklFKuX79ejhs3Tkoppbu7uzxw4ICUUsrw8HDp7u4upZTy7t278v79+1JKKRMSEqSfn592jtDQUC3+wsvTp0+XK1fq2uWXX36RQUFBUkopR44cKfft2yellPLixYvSxcWlRDsUPk5aWpps3bq1PHHihBw6dKjcvHlzme3Xp08fuXfvXhkfHy89PDweqs1N2e3bt8vdXvw7r1QdECurkJvUgw96Zm9vz5gxY4iIiMDa2lpbHxMTw6ZNmwB4/vnn+c9//qNtGzJkCGZmZri5uZGcnKyt79y5My1atAB0ZRQTExNp0KABJ06coG/fvoCuSHXTpk25efMmGRkZdO3aFYBRo0bx3XffAboqZlOnTuXo0aOYm5uTkJBQ4ecICwtj/vz5jBs3jvXr12vVz3bv3s2pU6e0/W7fvl3qI6j79u3Dx8cHMzMzZs6cibu7e7nnS05O5ty5c3Tv3h0hBBYWFpw4cULrpStKTaASrgG88sor+Pr6Mm7cuDL3EeKvssF169bVXstCtS8Krzc3NycnJwcpJe7u7sTExBQ53o0bN8o81+LFi2nSpAm///47eXl5WFlZVfgZAgMDOXfuHKmpqWzevJmoqChAV+MhJiamyC+T0vTo0UNL+AXc3d2Ji4tj8ODBJfbfsGEDN27coE2bNoAuka9fv5533nmnwlgVxViY1BhuZS6aOU6biuuZ07ieOY3jtKkltjeZGa5tdxhfMjE2fXu+tr1h2LOVirNRo0Y8++yzrFixQlvXtWtXbVxy7dq1dO/evVLHLODs7ExqaqqWcLOzszl58iQNGzbEzs5Oq3VbeAz01q1bNG3aFDMzM1avXk1ubi6gm3WirClchBAMHTqUGTNm4OrqioODA6CrWPbxx3/NflFW6cnSTJ06lVWrVvHrr79q69asWcOff/7JunXr2LlzJ4mJiSQmJhIXF6fGcZUax6QSrjTyi2aFvfrqq0XuVoiIiODLL7/E09OT1atXlzurb3nq1KnDxo0bCQ8Px8vLC29vb+0i2IoVK5g4cSKBgYFIKalfvz4AU6ZMYdWqVQQEBJCQkKAVK/f09MTCwgIvL69SL/SFhYWxZs0awsLCinyO2NhYPD09cXNz49NPP33o2Js0acL69et57bXXcHZ2xtXVlX379nH9+nUuXbpUZMaKNm3aYG9vXyQ5K4qxM6nyjAWcnZ1lfHx8kXWqVJ1uyp2CX0bvvvsu165de+TEXpwqz6h/qjyj/lW1PKMaw61Fvv/+exYuXEhOTg6tW7cmMjLS0CEpSq2iEm4tEhYWVuTPf0VR9MukxnArYorDJ4pSGvVdN04mlXDLu0vBysqK9PR09UVUTJ6UkvT09Ie6vU/Rr1pz0Sw7O5ukpCRtmnGlemVmZqp/4HpWXptbWVnRokWLSs/arJRPXTR7SJaWltpN80r1i46OxsfHx9Bh1CqqzWsekxpSUBRFMWYq4SqKouiJSriKoih6YpIXzYQQGUB8hTtWn/rALT2+/2H2L2+fsrZVZn3xdY0BfVZWV22u2ryifR5HmztLKR/9kcqq1HY01h+qWLPyEc63XJ/vf5j9y9unrG2VWV98nWpz1eaqzSv+UUMK1WObnt//MPuXt09Z2yqzvqqfuapUm+ufavMqMtUhhVhZhXvllMpTba5/qs31r6ptbqo93OWGDqAWUm2uf6rN9a9KbW6SPVxFURRjZKo9XEVRFKOjEq6iKIqeqISrKIqiJyrhKoqi6EmtSLhCiCFCiM+FEFuEECGGjsfUCSFchRCfCiE2CiEmGzqe2kIIUU8IESeEGGDoWGoLIUSQEGJf/vc9qKL9a2zCFUKsFEKkCCFOFFvfXwgRL4Q4J4SYCSCl3Cyl/CfwAqDmmHkElWzv01LKScCzgLpP9BFVps3zhQPf6DdK01PJdpfAHcAKSKrw4Pp8NLCaHzPsCfgCJwqtMwfOA22BOsDvgFuh7R8AvoaOvSb+VLa9gUHAQWCUoWOvqT+VaXPgSeA5dJ2KAYaOvSb/VLLdzfK3NwHWVnTsGtvDlVLuBa4XW90ZOCelvCClzALWA4OFznvADinlb/qO1RRUpr3z998qpewKjNZvpKajkm0eDAQAo4B/CiFq7L9tQ6tMu0sp8/K33wDqVnRsU5vxoTlwudByEtAFmIauB1BfCNFOSvmpIYIzQaW2d/5Y1jB0X8DtBojLlJXa5lLKqQBCiBeAtEKJQKkeZX3XhwH9gAbAxxUdxNQSrihlnZRSRgAR+g6mFiirvaOBaP2GUmuU2ubaCykj9RdKrVLWd30TsOlhD2Jqf3YkAS0LLbcArhooltpAtbf+qTY3jGppd1NLuIeB9kKINkKIOuguImw1cEymTLW3/qk2N4xqafcam3CFEOuAGMBZCJEkhHhRSpkDTAV+AE4D30gpTxoyTlOh2lv/VJsbxuNsd1UtTFEURU9qbA9XURSlplEJV1EURU9UwlUURdETlXAVRVH0RCVcRVEUPVEJV1EURU9UwlWMnhBithDipBDimBDiqBCiy2M+X7QQotJlJYUQzkKIyPxiSQcfR2xKzWZqtRQUEyOECAQGoCur+UAI0RhdeTxj1APYB3gC6mEEpQSVcBVj1xRd9asHAFLKtIINQog3gIGANbrauy9JKaUQIho4AvgBjsAYYBbgAWyQUs4RQjgBO4FfAR8gARgjpbxX+OT5M4S8ha7y2XlgnJTyTrF9egAfAa2AZMAOyBNCxEopVQF2RaOGFBRjtwtoKYRIEEJ8IoToVWjbx1JKfyllR3RJt/DUMllSyp7Ap8AW4GWgI/CCEMIhfx9nYLmU0hO4DUwpfOL83vQc4EkppS8QC8woHqCUcp+U0htd0nYDdgNPqWSrFKcSrmLU8nuTfsBEIBXYkF/zFSBYCPGrEOI40BtwL/TWgsIix4GTUspr+b3kC/xV9emylPJA/us1QPdipw9Al0APCCGOAmOB1qXFKYSwATKl7ln59kD8o3xexbSpIQXF6Ekpc9HV143OT65jhRDrgU+ATlLKy0KIeejmlSrwIP+/eYVeFywXfO+LFxIpviyAH6WUI8uLTwixFXABGgghjgFOQKwQYqGUckPFn1CpLVQPVzFq+Vf+2xda5Q1c5K/kmiaEsAVGPMLhW+VflAMYCewvtv0XoJsQol1+LDZCiA7FDyKlHAR8DkwGpgOfSim9VbJVilMJVzF2tsAqIcSp/N6jGzBPSnkTXZI7DmxGV6+0sk6j6y0fAxoBywpvlFKmopuUcV3+Pr+g68mWpie6hN0D2PMIsSi1gCrPqNRK+XcpfJd/wU1R9EL1cBVFUfRE9XAVRVH0RPVwFUVR9EQlXEVRFD1RCVdRFEVPVMJVFEXRE5VwFUVR9OT/AyJyhHuPvr3CAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "linewidth = 3\n",
    "\n",
    "t = list(range(samples))\n",
    "\n",
    "fig = plt.figure(figsize=(5,5))\n",
    "\n",
    "ax = plt.subplot(1, 1, 1)\n",
    "\n",
    "add_fill_lines(ax, t, bionica_indirect_err, plot_kwargs={'ls': '-', 'lw': linewidth, 'label': 'Bio-NICA (interneurons)'})\n",
    "add_fill_lines(ax, t, bionica_direct_err, plot_kwargs={'ls': '--', 'lw': linewidth, 'label': 'Bio-NICA (2-compart. neurons)'})\n",
    "add_fill_lines(ax, t, nsm_err, plot_kwargs={'ls': '-.', 'lw': linewidth, 'label': '2-layer NSM'})\n",
    "add_fill_lines(ax, t, npca_err, plot_kwargs={'ls': ':', 'lw': linewidth, 'label': 'Nonnegative PCA'})\n",
    "\n",
    "ax.loglog()\n",
    "\n",
    "plt.title(f\"{s_dim}-dimensional source data\")\n",
    "plt.grid()\n",
    "plt.legend(loc = 'lower left')\n",
    "plt.xlabel('Sample #')\n",
    "plt.ylabel('Permutation error')\n",
    "plt.xlim((1e2,samples))\n",
    "plt.ylim(bottom=10**-3.5,)\n",
    "# plt.ylim(10**-2.5,10**1)\n",
    "\n",
    "plt.savefig(f'plots/{dataset}.png', dpi=300, transparent='true', bbox_inches='tight')\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}