Linear example statistics.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "import numpy as np\n",
    "import itertools\n",
    "from SparseARD import *\n",
    "from AICc_Training import *\n",
    "import json\n",
    "np.random.seed(0)\n",
    "\n",
    "retrain = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 250\n",
    "n = 250\n",
    "s = 25\n",
    "noise_percent = 10\n",
    "kappa = 10**2\n",
    "n_trials = 100\n",
    "\n",
    "tol=1e-8\n",
    "verbose=False\n",
    "\n",
    "S = np.diag(np.logspace(-np.log10(kappa), 0, np.min([n,m])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if retrain:\n",
    "    \n",
    "    L1_err = {'ARD' : [], 'ARDr' : [], 'ARDvi' : [], 'mSTSBL' : [], 'lSTSBL' : [], 'dSTSBL' : []}\n",
    "    L2_err = {'ARD' : [], 'ARDr' : [], 'ARDvi' : [], 'mSTSBL' : [], 'lSTSBL' : [], 'dSTSBL' : []}\n",
    "    Added = {'ARD' : [], 'ARDr' : [], 'ARDvi' : [], 'mSTSBL' : [], 'lSTSBL' : [], 'dSTSBL' : []}\n",
    "    Missed = {'ARD' : [], 'ARDr' : [], 'ARDvi' : [], 'mSTSBL' : [], 'lSTSBL' : [], 'dSTSBL' : []}\n",
    "    \n",
    "    for trial in range(n_trials):\n",
    "    \n",
    "        print('Trial', trial+1, 'of', str(n_trials)+'.')\n",
    "\n",
    "        # Get random matrix\n",
    "        A = np.random.randn(m,n)\n",
    "        U,_,Vt = np.linalg.svd(A, full_matrices = False)\n",
    "        Theta = U.dot(S).dot(Vt)\n",
    "\n",
    "        # Get random sparse xi\n",
    "        nnz = np.random.choice(n,s,replace=False)\n",
    "        xi = np.zeros((n,1))\n",
    "        xi[nnz,0] = np.random.randn(s)\n",
    "\n",
    "        # Get output\n",
    "        y = Theta.dot(xi)\n",
    "        sigma = 0.01*noise_percent*np.std(y)\n",
    "        y = y + sigma*np.random.randn(m,1)\n",
    "\n",
    "        # Test each method\n",
    "        gamma_ARD,sigma2_ARD,mu_xi_ARD,_ = SBL([Theta,y], \\\n",
    "                                           sigma2=1, \\\n",
    "                                           estimate_sigma=True, \\\n",
    "                                           maxit=500, \\\n",
    "                                           verbose=verbose, \\\n",
    "                                           tol=tol)\n",
    "        AICc_ARD = AICc(Theta,y,gamma_ARD,sigma2_ARD)\n",
    "        ARD_results = gamma_ARD,sigma2_ARD,mu_xi_ARD,AICc_ARD\n",
    "\n",
    "        gamma_ARDr, sigma2_ARDr, mu_xi_ARDr = Train_ARDr(Theta, y, ARD_results=ARD_results)\n",
    "        gamma_ARDvi, sigma2_ARDvi, mu_xi_ARDvi = Train_ARDvi(Theta, y, ARD_results=ARD_results)\n",
    "        gamma_STSBL, sigma2_STSBL, mu_xi_STSBL = Train_M_STSBL(Theta, y, ARD_results=ARD_results)\n",
    "        gamma_lSTSBL, sigma2_lSTSBL, mu_xi_lSTSBL = Train_L_STSBL(Theta, y, ARD_results=ARD_results)\n",
    "        gamma_dSTSBL, sigma2_dSTSBL, mu_xi_dSTSBL = Train_MAP_STSBL(Theta, y, ARD_results=ARD_results)\n",
    "\n",
    "        # Record errors\n",
    "        L1_err['ARD'].append(np.linalg.norm(xi - mu_xi_ARD, 1))\n",
    "        L1_err['ARDr'].append(np.linalg.norm(xi - mu_xi_ARDr, 1))\n",
    "        L1_err['ARDvi'].append(np.linalg.norm(xi - mu_xi_ARDvi, 1))\n",
    "        L1_err['mSTSBL'].append(np.linalg.norm(xi - mu_xi_STSBL, 1))\n",
    "        L1_err['lSTSBL'].append(np.linalg.norm(xi - mu_xi_lSTSBL, 1))\n",
    "        L1_err['dSTSBL'].append(np.linalg.norm(xi - mu_xi_dSTSBL, 1))\n",
    "\n",
    "        L2_err['ARD'].append(np.linalg.norm(xi - mu_xi_ARD, 2))\n",
    "        L2_err['ARDr'].append(np.linalg.norm(xi - mu_xi_ARDr, 2))\n",
    "        L2_err['ARDvi'].append(np.linalg.norm(xi - mu_xi_ARDvi, 2))\n",
    "        L2_err['mSTSBL'].append(np.linalg.norm(xi - mu_xi_STSBL, 2))\n",
    "        L2_err['lSTSBL'].append(np.linalg.norm(xi - mu_xi_lSTSBL, 2))\n",
    "        L2_err['dSTSBL'].append(np.linalg.norm(xi - mu_xi_dSTSBL, 2))\n",
    "\n",
    "        _,added,missed = sparsity_err(mu_xi_ARD,xi)\n",
    "        Added['ARD'].append(added)\n",
    "        Missed['ARD'].append(missed)\n",
    "\n",
    "        _,added,missed = sparsity_err(mu_xi_ARDr,xi)\n",
    "        Added['ARDr'].append(added)\n",
    "        Missed['ARDr'].append(missed)\n",
    "\n",
    "        _,added,missed = sparsity_err(mu_xi_ARDvi,xi)\n",
    "        Added['ARDvi'].append(added)\n",
    "        Missed['ARDvi'].append(missed)\n",
    "\n",
    "        _,added,missed = sparsity_err(mu_xi_STSBL,xi)\n",
    "        Added['mSTSBL'].append(added)\n",
    "        Missed['mSTSBL'].append(missed)\n",
    "\n",
    "        _,added,missed = sparsity_err(mu_xi_lSTSBL,xi)\n",
    "        Added['lSTSBL'].append(added)\n",
    "        Missed['lSTSBL'].append(missed)\n",
    "\n",
    "        _,added,missed = sparsity_err(mu_xi_dSTSBL,xi)\n",
    "        Added['dSTSBL'].append(added)\n",
    "        Missed['dSTSBL'].append(missed)\n",
    "        \n",
    "    if True:\n",
    "        with open('./saved_results/Linear_results.json','w') as f:\n",
    "            json.dump({'L2_err': L2_err, \\\n",
    "                       'L1_err': L1_err, \\\n",
    "                       'Added' : Added, \\\n",
    "                       'Missed' : Missed}, f)\n",
    "        \n",
    "else:\n",
    "    with open('./saved_results/Linear_results.json','r') as f:\n",
    "        results = json.load(f)\n",
    "    L1_err = results['L1_err']\n",
    "    L2_err = results['L2_err']\n",
    "    Added = results['Added']\n",
    "    Missed = results['Missed']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3IAAAEtCAYAAABNkFLqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeZhcVZ3/8fcni2EJEYQQEQ0RBMGogDQ/ZVxoBxyGjI4LDuKAEEYNGtFBwXVYEkDRccRlUDQadnREDaCijDCSiA46RBE0I6IMBBDQDmIkCYTt+/vjnJbqSnV3dXfVvXWrPq/nuU933Xv6nlPLt+ueezZFBGZmZmZmZlYdk8ougJmZmZmZmY2NK3JmZmZmZmYV44qcmZmZmZlZxbgiZ2ZmZmZmVjGuyJmZmZmZmVWMK3JmZmZmZmYV44qcmZmZmZlZxbgi14EkbSbpMkm/kvRzSf8paeeyy2XWbSR9SNKvJT0u6TVll8fMzMysWa7Ida6zI2KPiNgL+BbwpbILZNaF/guYB/yg7IKYmZmZjYUrch1A0smSfiHpUICIeCgi/rMmyY+BUVvkJM2XFMNsf2pX+c2qoD7OACLiJxFx6xjP4zgzMzOz0rkiVzJJ84AFwI+Avxkm2TuBy8dw2n8A9qvbDpxAMc0qrck4GyvHmVmTJH0p3+w4cwx/s1zS8ibSLZIUEyrg0PMN3qyZ06pzmrVD3Y3F3Roc7685fmDed56k2wsv7Ahy+RaVXY4qmlJ2AYy3Al8EdgM2aRmQ9MF87IAxnPPnEfHbsRRC0rSI2DjWYxM9t1lBRoyzcXKcmTVB0uakGx8Ah0t6X0Q8WmaZzLrMA8CbgJPq9h+Zj21Vs+804NMFlcvazC1yJZI0DXgF8E3gr4D/rjt+AnAIcHBEbGhhvovy3Y/n5olU1gGXNHHsbyVdJ+lBSWvzhCzPbvbcZmUYLc7amK/jzCx5LTAD+A6wPfC35RbHrOssA46QpMEd+QbKIcA3ahNGxK0RcUPB5bM2cUWuXC8CHgceBp4CrBg8IOk9wBuBV0TEWMfdTJY0pW5r9F5fnvP8e+CTIx2T9LfAFcA64A3A24HnAj+UtOMYz21WpGHjbIIcZ2bNOQq4H5gPPEhqJRhC0mGSbpa0UdIqSa9tdCJJe0u6VtJDkn4n6SRADdJNkfTBmnPeLekTkjarS7ezpCskbZA0IOnTwLQWPGezIl0I7AS8pGbfa4HJ1FXk6rtW5lg5TdKtOa7WSPqhpJfUpPlHSTdIWpdvMP5C0jF1591f0n9JekDS+nyT8bl1aSZLOl3SPTnmlkua27qXofe4a2W5Xgz8FDgC+EpEPAgg6enAJ4D/A67JN1gejYi+Js97c4N9VwCvrNv3mYgYrnl9yDFJK3N5Dh7sEiPpOuAW4HjgPWM4t1mRGsYZgKQTgbcBM4HnSjoL6IuIe5s4r+PMbBSSnkYaO7okIgYkXQa8TtI2EXF/TnMg8GVS/BxPisdPA1OBX9ecazvg+8C9pMrhRuC9wOwGWV8EvAr4GKkVfg9Sl7I5pFYKJD0JuArYHHgH8AfgGOB1LXsBzIqxmjT78puAa/O+I4FLSTcGR/J+4N3AvwA/J7We95FufJIrdBcBnyHF2yRgd2DrwRNI+jvSjcUrSN+1g+e9VtLzI+LOvG8R8CHgTOB7OZ9vjuP5WuaKXLnmArcB/wS8bHBnRNxFgzuMY/Ba4K66fY1a9S4d4Rx/OSZpS+AFwEdqxzVExG2SfgTsP8ZzmxWpYZwBRMTpwOnjPK/jzGx0byJd+F2QH59P6m3yBuDzed9i0o2RV0fE4wCSfkWasfnXNed6N7AlcFBE3JHTXUW6iP0LSS/N5z8qIgbzvVrSH4GLJO0VET8nVQZ3BvaLiB/nv/0u8IsWPXezIl0AfELSu4BtSDdQDm7i7/YDvld3U/BbNb+/CPhTRBxXs+97def4NLAiIl49uEPSNaQbk8cDx0nahhTDSyLihMHzSHoM+GgT5bQG3LWyXE8DXk4KoF+PlngMfhkRK+u2RpMy3DPCOWqPbUOqWDZKfy/5rs0Yzm1WJMeZWXmOBH4TEdflx1cDd+f9SJoM7At8fbASB2lpEOD2unPtB/x4sBKX061n6EUnpDF4DwPfqO36zBMXn4M3dPYD7hysxOXzPY7Hm1o1fY3ULfhVwOGk743/auLvrgfmSfqwpJfklur649tIukjSKyVtXXtQ0q7ALsDFdfG2AbiOJ+LteaQbMfXx9R/NP0Wr54pcubYmdSH5YEn5jzRdc+2x+/PjpzZI91TgvjGe26xIjjOzEkjaF3gOsEzS1vkCcCvSxAz7KU2Xvh2pC+XvG5yift8OTabbHngSqUvZIzXbH/Lxbcd4PrOOFxEPAJeRWsGPBC6uvTkygo8Ap5DGWl8L3Cfp3NyVmYhYQZp19hmkXiADkq6W9Pz899vnn0sZGm+PkIYa1MYbbBpfjrcJcEWuXJOAT+eulB0r3/H8KfAP+e4pAJJ2Is0C2KrJI8zawXFmVo6j8s/3k25UDG7H5v1HAmtIF3yzGvx9/b57mkx3H/AQqaWv0faFMZ7PrCouAP6O1Pp1wShpAYiIRyLiYxHxPFJl692kcaSfrUnz9YjYn9Rz5LU53ZV5gq/Bm4wfpHG8vSofH+xBUh9fjrcJ8Bi5kkg6ihRov8wXbZ8EzoqIW1pw+r0G76TUWTmBtXtOIg1i/bakzwHTSeMa1pImZjHrOI4zs3Lk7lmHAT8BPtAgySd5Yt2r64HXS1pUM0buhaSJSWrHv10HvFfSMwYnT8hjS1/FUFeSKo9PjoiRupZdBxwt6UU1Y+QmAYeO5bmadZCrSF0X/xQRq8b6x3miry9JmkeaMbn++DrS99POpHFx25LGsd4OzI2Ikca63QSsJ8XX92v2HzbWctoTXJErgaQtSM3UryTNovVL4IsturiE1E+6kZmku59jFhFX5lmJTiH9k3gYWA68LyLuHs85zdrJcWZWqsEuVcdHxPL6g5K+AJwN9JM+798DLsv7Z5JuYNTPHvtJYCFpgoRFPDFr5YO1iSJiuaSvAF+XdCbwP6QlSOYA84D35/8D55MqmcskfYjU9fJtpFn7zConIh4jTSbUNEmXAzcCPyO1mO9NGmf6hXz8VFKr2TWk8a1PB94F/DwiBnKadwCX5xs4l5C+A2eRepPcERFnRsSfJH0S+BdJD5Bifl/gzRN60j1OER5iYWZmZq2TLw5fDjw1IjY0OP5kUlerSyJivqQ3kqYmnwP8FjgR+GeAiOiv+bsXkFoC9iV16fo86ab0yRFRuxjyJOCdpNlqn02q9N0O/Cfw4YhYm9PtDJxFqlCuJy2D8L/5vM+MiNsn/GKYtYmk+cC5wK7DTLaFpH5SJewVEXG1pPOA/oiYk48fT7rpuSuwBXAH8BVSnDySby6+C3g+adKtP5AqYSfV3mCUtB9pCYMXk5b0uJc08+ynByc7yj1jFgNvAZ5MarFfCKwCFkfEoha8LD3FFTkzMzMzM7OK8WQnZmZmZmZmFeOKnJmZmZmZWcW4ImdmZmZmZlYxrsiZmZmZmZlVTEcsP7DddtvFnDlzyi6GWUv99Kc/XRMRM8suxyDHmXWjToszcKxZd+q0WHOcWTcaa5x1REVuzpw5rFy5suximLWUpNWjpyqO48y6UafFGTjWrDt1Wqw5zqwbjTXO3LXSzMzMzMysYlyRMzMzMzMzqxhX5MzMzMzMzCrGFTkzMzMzM7OKcUXOzMzMzMysYlyRM+tQki6SdI+kP0u6RdJbao4dIOlmSRskXSNppzLLamZmNh6S1tVtj0n697LLZVYFrsiZda4zgDkRMQP4e+B0SftI2g5YBpwEPAVYCXy1vGKamZmNT0RMH9yAWcCDwNdKLpZZJXTEOnJmtqmIWFX7MG+7APsAqyLiawCSFgFrJO0eETcXXlAzM7PWeD3wB+DasgtiVgWVa5GTNOpm1i0kfU7SBuBm4B7gO8Bc4MbBNBGxHrg176//+wWSVkpaOTAwMNa8HWdmbebvNLMhjgIuiIhodHAi32lWjDL/p3Vqvu3Mu3IVuYgYsg23z6wbRMRCYCvgpaTulBuB6cDauqRrc7r6v18SEX0R0Tdz5syx5j1irJlVgaRj84XfRknnDZPmFEkh6cCafdMknZPHqN4r6T3tKJ+/08wSSbOB/YHzh0szke80K0aZ/9PKum4p8zm7a6VZh4uIx4AfSjoCeDuwDphRl2wG8EDRZTOrgLuB04GDgM3rD0rahdSd6566Q4uAXYGdgKcC10j634i4sq2lNetdRwI/jIjbyi6IWVVUrkXOrIdNIY2RWwXsObhT0pY1+82sRkQsi4jLgPuGSXIW8H7g4br9RwKnRcT9EfEr4IvA/LYV1MyOZITWODPblCtyZh1I0vaSDpM0XdJkSQcBbwS+D1wKPFfSIZI2A04GbvJEJ2ZjI+kfgIcj4jt1+7cBnkbNWNT8+ybjUGv+xmN3zMZJ0l8BO+LZKs3GxBU5s84UpG6UdwH3A/8GHBcRl0fEAHAI8OF87IXAYWUV1KyKJE0HPgIc1+Dw9Pyzdixqw3Gogzx2x2xCjgKWRYSHCJiNgcfImXWgXFnbf4TjVwO7F1cis66zGLhwmPE46/LPGcBDNb/7ItOsDSLimLLLYFZFbpEzM7NedADwrjwj5b3AM4BLJL0/Iu4nTX6yZ036PfE4VDMz6yBukTMzs64laQrpu24yMDmPK32UVJGbWpP0euA9wHfz4wuAEyWtBGYBbwWOLqrcZmZmo3GLnJmZdbMTgQeBDwBH5N9PjIj7IuLewQ14DLg/Iga7VZ4C3AqsBlYAH/fSA2Zm1kncImdmZl0rIhaR1oQbLd2cuscbgX/Km5mZWcdxi5yZmZmZmVnFjFqRkzRN0lJJqyU9IOkGSQcPk3a+pMckravZ+lteajMzMzMzsx7WTNfKKcCdpKnQ7wDmkWb2el5E3N4g/XUR8ZLWFdHMzMzMzMxqjVqRi4j1DB1f8G1JtwH7ALe3p1hmZmZmZmY2nDGPkZM0C9iN4dfT2VvSGkm3SDopT/1sZmZmZmZmLTKmSpakqcDFwPkRcXODJD8Ankuarnku8FXSej1nNDjXAmABwOzZs8dWajMzMzMzsx7WdIucpEnAhcDDwLGN0kTE/0XEbRHxeET8AjgVeP0waZdERF9E9M2cOXMcRTczMzMzM+tNTbXISRKwFJgFzIuIR5o8fwAaZ9nMzMzMzMysgWZb5M4G9gBeFREPDpdI0sF5DB2SdgdOAi6fcCnNzMzMzMzsL5pZR24n4BhgL+DemvXhDpc0O/8+OMjtAOAmSeuB7wDLgI+0q/BmZmZmZma9qJnlB1YzcvfI6TVpTwBOaEG5zMzMzMzMbBhjXn7AzMzMzMzMyuWKnJmZmZmZWcW4ImdmZmZmZlYxrsiZmZmZWakkHSbpV5LWS7pV0kvLLpNZp2tqHTkzMzMzs3aQ9ArgY8AbgP8Bdii3RGbV4IqcmZlZj5FGmow6iYgCSmIGwGLg1Ij4cX78uzILY1YV7lppZmbWYyJiyDbcPrN2kzQZ6ANmSvqtpLsknSVp87LLZtbpXJEzM7OuJelYSSslbZR0Xs3+F0m6StIfJQ1I+pqkHWqOS9LHJN2Xt39VM81YZjZWs4CpwOuBlwJ7AXsDJ9YnlLQgx/PKgYGBYktp1oFckTPrQJKmSVoqabWkByTdIOngfGyOpJC0rmY7qewym3Wou4HTgXPq9m8DLAHmADsBDwDn1hxfALwG2BN4PvBK4Jg2l9WsFz2Yf/57RNwTEWuAM4F59QkjYklE9EVE38yZMwstpFkn8hg5s840BbgT2B+4g/SFdomk59Wk2ToiHi2jcGZVERHLACT1AU+v2f/d2nSSzgJW1Ow6CvhERNyVj38CeCvw+XaX2ayXRMT9ku4C3J/XbIzcImfWgSJifUQsiojbI+LxiPg2cBuwT9llM+tSLwNW1TyeC9xY8/jGvK8hd/kym5BzgXdK2l7SNsBxwLdLLpNZx3NFzqwCJM0CdmPohebqPCj8XEnblVQ0s8qT9HzgZOC9NbunA2trHq8Fpg83Ts5dvswm5DTgeuAW4FfADcCHSy2RWQW4ImfW4SRNBS4Gzo+Im4E1wL6kcT37AFvl443+1q0EZiOQ9Czgu8A/R8S1NYfWATNqHs8A1oWnczRruYh4JCIWRsTWEfHUiHhXRDxUdrnMOp0rcmYdTNIk4ELgYeBYgIhYFxErI+LRiPh93v83kmbU/71bCcyGJ2kn4GrgtIi4sO7wKtJEJ4P2ZGiLuJmZWalckTPrULkL11LS1MyHRMQjwyQdbCHw1OhmdSRNkbQZMBmYLGmzvG9H4PvAZyOi0QQmFwDvkbSjpKcBxwPnFVZwMzOzUXjWSrPOdTawB3BgRAxOz4ykFwJ/An5DmkL9M8DyiFjb8Cxmve1E4JSax0cAi0k3QHYGTpH0l+MRMT3/+oV8/Bf58ZfyPjMzs47gipxZB8pdvo4BNgL31syvcAzwOPARYHvgz8BVwBtLKKZZx4uIRcCiYQ4vHuHvAnhf3szMzDqOK3JmHSgiVjNyV8mvFFUWMzMzM+s8HiNnZmZmZmZWMa7ImZmZmZmZVYwrcmZmZmZmZhXjipyZmZmZmVnFuCJnZmZmZmZWMa7ImZmZmZmZVYwrcmZmZmZmZhXjipyZmZmZmVnFuCJnZmZmZmZWMa7ImZmZmZmZVYwrcmZmZmZmZhUzakVO0jRJSyWtlvSApBskHTxC+ndLulfSWknnSJrW2iKbmZmZmZn1tmZa5KYAdwL7A08GTgIukTSnPqGkg4APAAcAc4CdgcWtKaqZmZmZmZlBExW5iFgfEYsi4vaIeDwivg3cBuzTIPlRwNKIWBUR9wOnAfNbWmIzMzMz6xqSlkt6SNK6vP267DKZVcGYx8hJmgXsBqxqcHgucGPN4xuBWZK2bXCeBZJWSlo5MDAw1mKYmZmZWfc4NiKm5+3ZZRfGrArGVJGTNBW4GDg/Im5ukGQ6sLbm8eDvW9UnjIglEdEXEX0zZ84cSzHMzMzMzMx6WtMVOUmTgAuBh4Fjh0m2DphR83jw9wfGVTozM7MJkHRs7v2xUdJ5dccOkHSzpA2SrpG0U82xaXnCrj/nCbzeU3jhzXrLGZLWSPqRpP6yC2NWBU1V5CQJWArMAg6JiEeGSboK2LPm8Z7A7yPivgmV0szMbHzuBk4HzqndKWk7YBlpAq+nACuBr9YkWQTsCuwEvBx4n6S/LaC8Zr3o/aQJ8nYElgDfkrRLfaLxDsuRNOrWLmXla72h2Ra5s4E9gFdFxIMjpLsAeLOk50jaBjgROG9iRTQzMxufiFgWEZcB9TcUXwesioivRcRDpIrbnpJ2z8ePBE6LiPsj4lfAF/HkXWZtERE/iYgHImJjRJwP/AiY1yDduIblRMSQbbh97TBa3mYT0cw6cjsBxwB7AffWzCh0uKTZ+ffZABFxJfCvwDXA6ryd0r7iF8t3VczMusaQybkiYj1wKzA334h8GptO3jW30BKa9a4AfGFlNoopoyWIiNWMHEzT69KfCZw5wXJ1pPo7J5J8N8XMrJqmA/V9s9aSJueaXvO4/lhDkhYACwBmz57dulKadTlJWwMvBFYAjwJvAF4GHFdmucyqYMzLD5iZmXWB+sm5yI8fyMdg08m7hp24yzMxm43bVNI41gFgDfBO4DUR4bXkzEYxaoucmZlZF1oFHDX4QNKWwC6kcXP3S7qHNGHXVTnJnjReP9XMJiAiBoB9yy6HWRW5Rc7MzLqWpCmSNgMmA5MlbSZpCnAp8FxJh+TjJwM31ayRegFwoqRt8gQob8WTd5mZWQdxRc6sA+U1rJZKWi3pAUk3SDq45viw61+Z2RAnAg8CHwCOyL+fmFsBDgE+DNxPGqNzWM3fnUKa/GQ1aezOx/OEXmZmZh3BXSvNOtMU4E5gf+AO0jTMl0h6Hmn8zjLgLcC3gNNI61+9qJyimnWuiFhEWlqg0bGrgd2HObYR+Ke8mZmZdRxX5Mw6UJ4KfVHNrm9Lug3YB9iWvP4VgKRFwBpJu9d0CzMzMzOzLuaulWYVIGkWsBtpsoVh179q8HcLJK2UtHJgoH6mdTMzMzOrKlfkzDqcpKnAxcD5ucVtOkPXt4Jh1rjylOhmZmZm3ckVObMOJmkScCHwMHBs3j3S+ldmZmZm1gNckTPrUJIELAVmAYdExCP50CrSmlaD6f6y/lXhhTQzMzOzUrgiZ9a5zgb2AF4VEQ/W7B9t/SszMzMz63KuyJl1oLwu3DHAXsC9ktbl7fAm1r8yMzMzsy7n5QfMOlBErAY0wvFh178yMzMzs+7nFjkzMzMzM7OKcYucmZmZFSbN4zS8iCioJGZm1eaKnJmZmRWmvqImyZU3M7NxcNdKMzMzMzOzinFFzszMzMzMrGLctdLMOspo42fAY2jMzMzMXJEzs47i8TNmZmZmo3PXSjMzMzMzs4pxRc7MzMzMSidpV0kPSbqo7LKYVYErcmZmZmbWCT4LXF92IcyqwmPkzMysZ0maA3wO2A/YCHwdOC4iHpW0F7AU2AP4FfDmiPh5SUU162qSDgP+BPw38KySi2NN0vLlIye45ppR00R/f+vzbSLv8eTbadwiZ2ZmvexzwB+AHYC9gP2BhZKeBFwOXARsA5wPXJ73m1kLSZoBnAocP0q6BZJWSlo5MDBQTOEqStKom1WfW+TMzKyXPRM4KyIeAu6VdCUwF+gnfUd+KtK0qZ+RdALw18CVZRXWrEudBiyNiDtHqmBExBJgCUBfX5+nMx5B0TNAt61VrUPz7RRukTMzs172aeAwSVtI2hE4mFRRmwvcFEOvfG7K+zfhlgKz8cldmA8EPll2WcyqxhU5MzPrZStIlbM/A3cBK4HLgOnA2rq0a4GtGp0kIpZERF9E9M2cObONxTXrOv3AHOAOSfcCJwCHSPpZmYUyq4KmKnKSjs13GjdKOm+EdPMlPSZpXc3W36rCmpmZtYqkScB/AsuALYHtSOPhPgasA2bU/ckM4IEiy2jWA5YAu5DGqO4FfB64AjiozEKZVUGzLXJ3A6cD5zSR9rqImF6zLR936czMzNrnKcAzSGPkNkbEfcC5wDxgFfB8DR2w8/y838xaJCI2RMS9gxvpJspDEeE+ymajaKoiFxHLIuIy4L42l8fMzKwQEbEGuA14u6QpkrYGjgJuBJYDjwHvkjRN0rH5z75fSmHNekRELIqII8ouh1kVtGPWyr0lrQH+CFwInBERj9YnkrQAWAAwe/bsNhTDzMxsVK8DPgW8n1RxuwZ4d0Q8LOk1wJeAj5LWkXtNRDw83oy0uInpvheNni5O8WR9ZmbW+orcD4DnAqtJg8e/CjwKnFGf0FPImplZ2fIC3/3DHLsB2KfQApmZmTWppRW5iPi/moe/kHQq8F4aVOTMzMx60Xhb1Jpq0TMzs57R7uUHAvA3j5mZmZmZWQs1u/zAFEmbAZOByZI2k7RJa56kgyXNyr/vDpwEXN7KApuZmZmZmfW6ZlvkTgQeBD4AHJF/P1HS7LxW3OBsJQcAN0laD3yHtDbPR1pcZrOeMNz6jZLmSIq69RpPKrGoZmZmZlawpsbIRcQiYNEwh6fXpDsBOGHCpTIzeGL9xoOAzRsc37rRjLBmZmZm1v3asfyAmbVARCwDkNQHPL3k4piZmZlZB3FFzqy6VksK4CrgvXlxYzOzTYw64+Uir19nZlY17Z61csKkkTeIJtKYdZU1wL7ATqQ1rrYCLm6UUNKCPM5u5cDAQIFFNDMzM7N2coucWcVExDpgZX74e0nHAvdImhERf65LuwRYAtDX1+fb6WY9bjytal6/zsysM1WmIhfjuAR1a5z1iMHo8CfezMzMrEdUpiJn1mvyWo1TqFm/EXiU1J3yT8BvgG2AzwDLI2JtWWU1MzMzs2K5ImfWuU4ETql5fASwGPg1aX3G7YE/kyY7eWPhpTMzMzPrIVq+fPRE11wzYrro729VcVyRqwI10Uc0xtP31DraKOs3fqW4kpiZmZlZp3FFrgLqK2mSXHEzMzMzMyvBeFrVmmrNG6OOX37AzMzMzMzMhnJFzszMzMzMrGJckTMzMzOz0ki6SNI9kv4s6RZJbym7TGZV4DFyZmZmRTkX+q/pH7Lr0EMPZeHChWzYsIF58+Zt8ifz589n/vz5sB7YspBSmhXtDODNEbFR0u7Ackk3RMRPyy6YWSdzRc7MzHqepMNIy33MBu4F5kfEtZIOAD6b9/8k719dXknNuk9ErKp9mLddgHFV5EadVGKU6eGhjZNZtGFq+lbkO968ATjuOPq33nrIrmZuUDFnzvjyKzvfDuKKnJmZ9TRJrwA+BrwB+B9gh7x/O2AZ8BbgW8BpwFeBF407s6Nh+SnLGx7aYostWD7ShZZb46yLSfocMB/YHLgB+E6DNAuABQCzZ88usnhmHckVOTMz63WLgVMj4sf58e/gLxeNqyLia/nxImCNpN0j4uZSSmrWpSJioaR3AvsB/cDGBmmWAEsA+vr6Rl2Hqawp4quUb0vy/tSnWD5M3iPdoDq6qvl2EE92YmZmPUvSZKAPmCnpt5LuknSWpM2BucCNg2kjYj1wa95ff54FklZKWjkwMFBU8c26SkQ8FhE/BJ4OvL3s8ph1OlfkzMysl80CpgKvB14K7AXsDZwITAfW1qVfC2xVf5KIWBIRfRHRN3PmzPaW2Kz7TSGNkTOzEbgiZ2ZmvezB/PPfI+KeiFgDnAnMA9YBM+rSzwAeKLB8Zl1N0vaSDpM0XdJkSQcBbwS+X3bZzDqdK3JmZtazIuJ+4C7SLHn1VgF7Dj6QtCWplWBVg7RmNj5B6kZ5F3A/8G/AcRFxeamlMqsAT3ZiZma97lzgnZKuBB4BjgO+DVwKfFzSIcAVwMnATZWd6GSca9iZtVNEDAD7l10Osypyi5yZGSBp1M261mnA9cAtwK9IU59/OF9gHgJ8mNRS8ELgsLIKaWZmVsstcmZmQMTQnnWSNtln3SkiHgEW5q3+2NXA7oUXqh3GuYbd0YuPbl+ZzMxs3CpSkeunfpmI5rqDzB93js3dfI8R0/ka0MzMzMzM2sFdK83MzMzMzCqmIi1yyxluEfYRu4O0oDfIeFrVPJTGzApcqT8AACAASURBVMzMzMzayS1yZmZmZmZmFVORFjkza6smmpFjtHQeFGpmZmZWmKZa5CQdK2mlpI2Szhsl7bsl3StpraRzJE1rSUnNzMzMzMwMaL5r5d3A6cA5IyWSdBDwAeAAYA6wM7B4AuUzsyJFjH0zMzMzs8I1VZGLiGURcRlw3yhJjwKWRsSqiLiftMjq/IkV0czMzMzMzGq1eozcXODymsc3ArMkbRsRo1UCzczMup4We2pjMzObuFbPWjkdWFvzePD3reoTSlqQx92tHBgYaHExzKpvpLGpkg6QdLOkDZKukbRTScU0MzMzsxK0ukVuHTCj5vHg7w/UJ4yIJcASgL6+Pg+0sY6gZmZvLG5c2ODY1IOAzQd3StoOWAa8BfgWqQvzV4EXFVUwMxu7OGX0/x2SivwfY2ZmFdbqFrlVwJ41j/cEfu9ulVYVETFkG25fQWUZbmzq64BVEfG1iHgIWATsKWn3wgpnZmZmZqVqdvmBKZI2AyYDkyVtJqlRa94FwJslPUfSNsCJwHktK62ZQRqLeuPgg4hYD9ya9w/hLsxmZmZm3anZrpUnAqfUPD4CWCzpHOB/gedExB0RcaWkfwWuIXUF+0bd35nZxE0H6mtla2kwFtVdmM2sTE1N7LJo5HTNdEk1M+tFTVXkImIRqftWI9Pr0p4JnDmhUpnZSOrHopIfbzIW1cyslmfMNDPrHq2e7MTM2m8Vac1GACRtCeyS95uZdZzxtKq50mlmNrLKVOSamEzQrKvkcahTqBmbCjwKXAp8XNIhwBXAycBNEXFzaYU1s442WkXKs2VaWSRNAz4HHAg8Bfgt8KGI+G6pBTOrgMpU5MrRT3//0D2HHnooCxcuZMOGDcybN2+Tv5g/fz4wv4CyWQ9oODY1IhblStxZwEXAT4DDSijfxDWz3MNo6XzxaS0gaVfgF8DXI+KIvO8fgTOA7YCrgH+KiD+WV0qzrjQFuBPYH7gDmAdcIul5EXF7mQUz63StXn6g5SJG3kBNpDGrnohYFBGq2xblY1dHxO4RsXlE9PvLzmzCPgtcP/hA0lzgC8CbgFnABlKrgZm1UESsz993t0fE4xHxbeA2YJ+yy2bW6dwiN6LlLF/e+MgWW2zB8mEOHn102wpUqA5bHNvarB+ob4JupgV6fisyH8/nyP2trUUkHQb8Cfhv4Fl59+HAtyLiBznNScCvJG0VEZ5YyKxNJM0CdmOi476PO47+rbcesqupXlVz5kwo26rlO3/+fFi7Fp785Allr+EumNuslHzLeo8b6PgWOStPmYtjSxpxMzNrBUkzgFOB4+sO1a/XeCvwMOkC08zaQNJU4GLg/Ebjvr02qtlQbpGzjlRfSfRA/PZbDgzXBD1SC3TXNEFbrzoNWBoRd9bdJJpOWp+xVsP1GiFdYAILAGbPnt2GYpp1N0mTgAtJN0yObZRmTGujfupTLK+f6CAbsVfVRFt4KpYvMKHWuBgmz0Htun4bLd925l3ae9yAW+TMzKwnSdqLNFPeJxscHtN6jRGxJCL6IqJv5syZrS2oWZdTuouylDQe9ZCIeKTkIplVglvkzMysV/UDc4A7cmvcdNJSH88BrgT2HEwoaWdgGnBL4aU0635nA3sAB0bEg2UXxqwqXJEzM7NetQT4j5rHJ5Aqdm8Htgeuk/RS4GekcXTLPNGJWWtJ2gk4BtgI3FvTxfmYiLi4tIKZVYArcmZm1pMiYgNpWQEAJK0DHoqIAWBA0ttIEy9sC1wNeECoVcpok4N1wtjziFgNeBYzs3FwRc7MzIy0dmPd4y8DXy6nNGYT54nDzNqjrOUW6rkiNwrPdG9mZmZmZp3GFbkO01zFMUZM55ttZmZmQzWzBqlbq8xsJKUue9CAK3LDaOb1dxcFMzOzanA3QzPrNq7IdajxfLe4G6iZmZmZWW9wRc7MzMza51zov6Z/yK5DDz2UhQsXsmHDBubNm7fJn8yfP7+YspmZVdiksgtgZmZmZmZmY+MWOTMrVT9A3eDhZu7Wzy+gbGbWAkfD8lOWNzy0xRZbsHyYabyPXuxl+8zMRuIWOTMzMzMzs4pxi5yZlWo5wDB35Ee6W8/RvltvZmZmvcsVOTN7gqc+tQJ4Pa8nNHot6vf1ymthZmZjU7mKnL/0rFW8+LpZObye1xN69XmbmdnEVa4i5y89szZoIq56+WLbzMzMrNNUriLXG/rrJ/Frcs2d+QWUrft48XUzMzMzqxrPWmlmZmZmZlYxbpHrSMuHm8Rv5DV3PImfmZmZmVlPcIucWUVJWi7pIUnr8vbrsstk4yNpxM3MzMysnlvkzKrt2Ij4UtmFmLAer6x4FkczMzMbq6Za5CQ9RdKlktZLWi3pH4dJt0jSIzUtBOsk7dzaIpuZtYA04hZNpDEzs4mTdKyklZI2Sjqv7PKYVUWzXSs/CzwMzAIOB86WNHeYtF+NiOk12/+1oqBm1tAZktZI+pGk/rILM2YRo24aLZ3ZBEiaJmlpvkn5gKQbJB1cc/wASTdL2iDpGkk7lVlesy51N3A6cE7ZBTGrklErcpK2BA4BToqIdRHxQ+CbwJvaXTgzG9H7gZ2BHYElwLck7VKbQNKCfJdz5cDAQBll7HxNVCZdeexqU4A7gf2BJwMnAZdImiNpO2BZ3vcUYCXw1bIKatatImJZRFwG3Fd2WcyqpJkxcrsBj0XELTX7biR96TXyKkl/BO4BzoqIsxslkrQAWAAwe/bs5ktsZgBExE9qHp4v6Y3APODfa9IsIVXy6Ovrcw3ErE5ErAcW1ez6tqTbgH2AbYFVEfE1SMMHgDWSdo+Im4sua5Vpsbsit1tTr/GikdPFKd31NaHhpgB3vtYlmulaOR1YW7dvLbBVg7SXAHsAM4G3Aifni8tNRMSSiOiLiL6ZM2eOochmNowAfLVkNgGSZpFuYK4C5pJuXAJ/qfTdmvfX/51bv62h0Wal9cy0zXOcmQ3VTIvcOmBG3b4ZwAP1CSPif2se/rekTwOvB74y7hKatVU//f1D9xx66KEsXLiQDRs2MG/evE3+Yv78+cD8Aso2PElbAy8EVgCPAm8AXgYcV2a5rFqauYDspdkzJU0FLgbOj4ibJU0H6q8WG97IdOt3Y8208HT7LK1Fz0o7nla1qrSYNhtnUf/FXqdd78Fo+bYr77LytfI1U5G7BZgiadeI+E3etyfpbuVo3EJgo2ruZmSMmK4H/zdNJQ0M3x14DLgZeE1EeC05a5qXPXiCpEnAhaSJvY7Nu5u+kWlmZla0UStyEbFe0jLgVElvAfYCXg38VX1aSa8GfgD8CdgXeBfwoZaW2KylljNcl/ItttiC5cMcPProthWoKRExQIoxM5sgpabJpaSZmedFxCP50CrgqJp0WwK70NyNTDNrkqQppGvSycBkSZsBj0bEo+WWzKyzNbsg+ELSlLB/IM0o9PaIWCXppcB3I2J6TndYTjcNuAv4WESc3+IyW5uM3jI2cqsYTKxlbDx/66EFNhH9QH3f2ma61s4voGxWqLNJ47sPjIgHa/ZfCnxc0iHAFcDJwE2e6MSs5U4ETql5fASwmKETEZlZnaYqchHxR+A1DfZfS5oMZfBxw4lNbOxcQTEza7+8LtwxwEbg3ppxg8dExMW5EncWcBHwE9INSzNroYhYhCttZmPWbIuc9RC3jBk0ngijdl83jKVaDgzXt3akrrWl9621lomI1YwwljsiriaNRTUzM+sorsh1mGaujXt5QgIrjj9j1jKj3OmJJtL04oxGZmZmI2lmHTkzMzMzMzPrIG6RMzOzYrjftpmZWcu4ImcdopoLc5uZmZmZlcEVOTMzM7NucC70X9M/ZFdzN0XNrIpckbM6ZbWMVXNhbjMzMzOzMrgiZz3PQ3DMzKwrHA3LT1ne8NCIN0UX+66oWRW5Imd13DJmZmZmZtbpXJGrgNEWZgav+TUeXrPPADfJmpmZWSW5IlcBvVKR8PW0WRs1EWCjLsw9gf9F/UD9ANxmxt/OH3eOZmZm3c0VOTPrTaNUStwaa2bjocWj3DRZNHqaOMX/e8xsdK7I2SaKbhlzF0frKSW3jJW1KPdyYLgBuCONv/UAXDMzs8ZckRuD0caquaJhZmZmML5WtVFb88zMargiNwbdXlEb7em5VcyshUpoGesHj1MzMzPrEq7ImZmZWVdpqmVr0cjpPE7NzDqdK3JmZhS7zEc/lNIythzKHafmqWmtV5wL/df0D9nVTIybmY2FK3JmZnR/1+m/cGXKekgvjlOrevnNrHmuyJlZRymyZawsywFWrBi6c8UKeMc72GLweL369FVS4aUeJD0FWAr8DbAG+GBEfLncUlnHOxpWMDRmVwys4B2L35EevHzTP1mxusIxPkGOM7PxcUXOzDpKp17QV14Tr2snV6hK9FngYWAWsBdwhaQbI2JVucWqLs8A3R7NtD52cIw7zszGwRU560j+oreu5QpVZUjaEjgEeG5ErAN+KOmbwJuAD5RauAor7LNd0ji10SpUju+hHGdm4+eKnA2rzC5u/pIza79e6MY6QbsBj0XELTX7bgT2L6k8NkYrbq/r3njFCt4x8I7U9nN7g/SXrYC9CymaPcFxZjZOrsjZsHrxAq5KF7YeU9A9ymqBLuuzXKE4mw6srdu3FtiqPqGkBcACgNmzZ7e/ZBNUofdg/EaabPVJoxxvoTJf64r0bml7nPXie1DWc/Zr3Xhfu/J2Rc6sRod8qTXLYwq6RMU+dxNWoee7DphRt28G8EB9wohYAiwB6Ovr6/gnWKH3YFw6aQ24Ml/rirzPbY+zXnwPei3fMvMu8zlPKi1nMxu3mjEFJ0XEuoj4ITA4psDMWuMWYIqkXWv27Qn4ZolZ6zjOzMbJFTmzahpuTMHc2kSSFkhaKWnlwMBAoQU0q7qIWA8sA06VtKWkFwOvBi4st2Rm3cNxZjZ+rsiZVVNTYwoiYklE9EVE38yZMwsrnFkXWQhsDvwB+ArwdndfNms5x5nZOHiMnFk1NT2mwMzGLyL+CLym7HKYdTPHmdn4NNUiJ+kpki6VtF7Sakn/OEw6SfqYpPvy9q9qNJWLmU2UxxSYmZmZ9bBmu1bWzo53OHC2pLkN0i0g3VHZE3g+8ErgmBaU08xqeEyBmZmZWW8btSI3xtnxjgI+ERF3RcTvgE8A81tYXjN7gscUmJmZmfWoZsbIDTc73v4N0s7Nx2rTNWq5q9ziqWadxmMKzMzMzHqXRlvETtJLga9FxFNr9r0VODwi+uvSPgbMjYib8+NdSWN5JsUIGUkaAFaP8zlsB6wZ599OVFl591q+ZeY9kXx3ioiOmSpygnEG1XwPqphvmXlXMd+OijOo7HeaP+/dn+9E8+6oWKtonJWZd6/lW2behcVZMy1yY5kdrz7tDGDdSJU4gIn8Y5C0MiL6xvv3E1FW3r2Wb5l5l/mcW22iX8C99h748979+bZLFb/T/Hnv/nzLzrvVqhhnZebda/mWmXeR+TYz2clYZsdblY+Nls7MzMzMzMzGadSK3Bhnx7sAeI+kHSU9DTgeOK+F5TUzMzMzM+t5zS4/0HB2PEkvlbSuJt0XgG8BvwB+CVyR97XTkjafvxPz7rV8y8y7zOfcaXrtPfDnvfvz7US9+B702nPuxde60/Tie9Br+ZaZd2H5jjrZiZmZmZmZmXWWZlvkzMzMzMzMrEO4ImdmZmZmZlYxrsh1EUkquwxm3c5xZtZ+jjOzYjjWqs0VuRaSNLnucWHBIenZwDJJzawN2Mp8/17S9CLz7FSSXlB2GXpFWbFWVpzlvB1rOM6K5DjrXY6z4vjasbdNNNZckWsRSYqIxyRNlfRhSdNGWwi9xf4P+HZEPFpUhpKuBI6MiHWjJm5tvv8haV6ReY5G0reA08ouRy8oOdYKjzMoJ9YcZ73NcVZovh0Va46z4vjasTidFmfQmljrmlkrc412V2A9cHdE/CzvV7uDQtLkiHgs/34+MDkijmhnnnX5TxkMQkmTgE8BH2pnkEi6BHhGROyXH7f9dc75XA48NSJe2O68miVpGbBNRLy87LK0W5lxlvMpLdbKiLOcV+Gx5jgrl+OsN+Is59NRsdZLcQa+dvS1Y3laFWuFd1toB0lfAnYGngo8CXhc0pKI+LeIiHZ/UPLdlEnA54Hfk2vXRX1Agdo8BBwOPFvS6/KC7i0l6SLgEKC/1eceJd8LgN0iYo/8eGpEPFJkGRqU6WLgNcBO+fG0iNhYZpnapew4g9JjrdA4g3JizXFWLsdZb8RZzrejYq2X4gzKjzVfOxaj0+Isl6FlsVb5rpWSLgPmAkcDzwcOIi3Ed5qkDwIUFBBPznkfR/qnAAW9vvmfgST9S767swMwB7hc0patzCvf1dgD+AbwZkkH5jKE1L5+3ZK+TvoH8HtJu+XdhXa7qZdfi72Bc4GzJO0aERvzP+au0kFxBiXFWpFxBuXEmuOsXI6z3oiznG9HxVovxRl0VKz52rGH4gzaEGsRUdkN+DpwQ4P9mwNvB34HzGtT3i8A3gC8Etgn79seuAm4uibd5IJeix2AtcCC/HgqcDNwNbBli/L4IvCz/PsUUkB+BeivSaM2PLfvAtcBzwbOAK4F9m1Xfk2W6TM1r8WzgfOAK4E5ed+kMsrVpudaWpzlfDom1oqIs3zewmPNcVbu5jgbUpaujbN8zo6KtV6Ks/x8fO34RD6+diz2s9fyWCv8SbTwxdg1vylfG3xDgCk1x58GXAO8pw15fwn4PvC/wG+B3wDvzcdmAr9ud0A2erOBo4BvAs/Mj6cCvwCuB7aYYH4zgR3r9m2XX/+2BSSwP/CFmsd7kvpx/6CsgASeAby+9v0F9gXO6bYvvzLjLJ+/1FgrOs5qnlehseY4K3dznPVGnOVzdVSs9VKc5efha8dN9/nasZjPXltircpN5rcCJ5GaSL8haZuIeFR5CtWIuBt4kPTmtcwwzfFfAE7NzdMDwEuBWZKuz2V5rJVlyOd8PDeJP6tm9y+ArUkfFiL1Ad6H1A96u/HmJekc4GJgm/x4stIg3TXA20h3WN4maf9W9u2W9CbgIeDY/FgRcSOp+8PPgE9I2jeivU3zdWX6POm9/3l+PDUiHouI64GzgXuAz0uak9+jKscYlBRn0BmxVmScQTmx5jjrCI6zLo+znE9HxVoPxhn42tHXjt32nVZkbbRFNdrpwHY1tdmXA5eQmsq3yfuflH9eSJritFV5N9Mc/3d53w6kJt2dWpBvw7sGwJnAVcC/1zznD5Ka6Ddr0XO+nHT3aldg27pjk/LPbUl3Vr4DvKSFr/X1wIGDz6X2dQCeQ7q78n1gv4I+e5cBPwL+H7BV/euQf9+X1I3gOmB2EeVq03MtLc5q3v9CY63MOMvnLDzWHGflbo6zIfu7Ns5qXuuOibVeirP8XHzt+MR+Xzt20Xda259Ai1+MzwP/CfwUOHzwhSA1oV5C6nf7lLz/OFIz9ZwW5d1sc/zxNfsm3Cxe9wGs7Vv9TNIg2QPzh/Za4FTSbEDn0YL+3cBbgGvr9m0JzAKm1u2fmfPdsQX5XgT8kDSL0rCvIWng7JeAK4DNhvun1aL3/0PAT+r2Tcs/n1S3fz/grLEGY6dsZcZZPmfhsVZmnOU8C481x1m5m+OsN+Isn6ujYq2X4iw/B187+tqxNl1Xfae1peBtejEuB1aQpuv8MLAReFk+pvwhvITU1/QDwH3A3i3Mf1LO4yvA4NoP9QH5HeD8FuZZG4i1fatvzds/1xx/C/BZ0qDVx4ElLcj/BODLgx844EXAT/I/uUuBV+Rjkwdfoxbk+TTge+Q7Z3nfjqQm6cOBZ9Wl3x2YVcDn7+M8MRh4Gqnbxffye/5l0voktelbdge5yK3sOBv8HBUZa2XHWT5vobHmOCt3c5z1Rpzlc3RcrPVKnOWy+9rR145d/Z3W1ifQwhfiC8BdNY93AP4bOKhm3yTgJaQ7Gw8BL2hR3qU0x9cF4mXkpnZSn+JdcqA8CJxc93d/mwP3uRPNGzg0/xP8BGl9k9+T7ha8Dvgv0h2mVg963w34M+mOydbAPwD3A78k9eNeCexSwmfwE6TByfOAk0n/7M8FTs8BeUH+h1XZAeFlxlk+d+GxVmac1eZfdKw5zsrbHGe9E2c5z46LtV6Is/w8fe3oa8eu/04r9AmN80XYltTsfT7w8rxvR+CB/MZ8GDgCmJGPPZ8WdQGg5K4v+bzN9K1+Zd2xKS3Ke4f8wfsB8DngsJpjB5L6U2/birzq8j0r/6O5FrgReF/e/0Lgf4A3lvA53Jp09+R3pDtr82uOfQBYVnSZWvz8SouzfL6yu5mVFmf5XIXHmuOs+M1x1ntxls/dUbHW7XGWn4evHX3t2BPfaVPocBFxn6R/Id1FOFzSzsAppDfsNtJsOB8BTpD004h4cyvyzQv2bQ18kjw9qKQ7I+IHkn5AapJfCPybpFuA9wIHRsTtrcg/l2FXUp/i3w7O6CNpSkQ8GhEP5jIeSroTMTgzT0TEhBc7zOe6R9JHSa/vtIhYX5PkBcDdpDtYrfbPpLtIfwZ+HxH/AxARP5E02PxfqIj4E/CPkmYB6yNiXc3hqcAjkjYDNkaO0CopK86g/FgrM85qzldGrDnOCuY468k4gw6LtW6PM/C1I752LD3Oct7tj7WiaqVj3Uj9mSfxRNPzXGApcC9wcV3ap5KarZ/VorxLaY6nbsAlBfatJk3H+6yaxw3vzJAGyR5Hah7esw3v+7DNy6S7SKtp8Z2r0d4PagbnUjM4Nb8/7wD+CDyviDK14TmWFmf5nIXHWplxls9Veqw5zordHGe9GWeDz3uEY4XFWi/EWX4evnb0tWP9sa7+Tmv7Exrni3A8qeZ8Can5e3bev2MOyKW0aBarBnmX0hzP0H7NhfatJq218Tjw4/oPev4n8Ob8ez+pv+8vgb3a8NpPrvn9JeSFKIEDSP2J/9CKf3qjlOFw0l2dvgZlWgwckH9/cX7tf0OLJyEoaiszznI+hcdamXGWz1N6rDnOit0cZ70ZZ4PPt+b3wmOtl+IsPw9fO/rasee+09r2pCb4gvx1/nBcSVpM8AHgX4C9SYsTLiUt7HdIm/LfnTTo80vAm4E7gDOABcD78+OfA0tblF9tIBbet5o0He2K/GH/LU+sLn8psKom8LfJgdGKaWKH7Yud/+lcQeqCMB14E2lw7Nw2f+4uJQ2K/QHwGPCqmmPnA7cPlhvYnjSd7zPbWaY2P99S4yyXobBYKzvO8rkKjTXHWfmb46z74yyfq6NirdfiLD8PXzv62rHnvtPa9sRa8MJ8JgfcNqQ+xF8mNUO+LwfIN4CPAdNblF8pzfF1gVj4NLmkuzbT8uv6qnzu3wF9+R9Dw/UuJpjnGaS7F9MaHPtI/kc3ta6MW7b58/ZN4Kqax58DbiCtMfIy4LuDZaKuG0OVt6LjLOdZeKyVHWf53IXGmuOsczbHWffGWT5XR8Var8ZZfj6+dvS1Y20Zu/47rW1PrgUvzjtIdwxm5cfz84fzFtIK6b8CntaivErt+pLzKrRvNWnhwdqm36NJMwltRrqDcCuwcz42dbz5DJP33+T39rX1AUm6o1bbj7vtXzLA3+f3/7V1r8fFPLHOySYLeXbDVmSc5fOX3c2sjLFCpcSa46xzNsdZ98ZZPl/HxFovx1l+Tr529LVj2+Ms59ERsdbWJ9mCF+kGUr/aFwN3kZojNyMtLtjKKZnLbo4vtG81cF5+vktIa3xMyvs/Sl4okrT+yK20cMB9XRlelv+xvq4+IPPxPmDzgj5nc4ATSVPWHpg/Y7cDxxaRf9lbUXGW8yot1oqOs3yOUmPNcdY5m+Ose+Msn78jYq3X4yy/Br529LVjEZ+zjoi1wZpiR5E0KSIelzSfNCPOk4ETI+LzbczzM6Q34f2kle73Ji2Q+FFSM/Q8Uh/g02Lo9KHjyWtSRDxet2930jS5kKZPPYX0j+E2UheBwRluJjxNrqS/Iq0sfwGpj7NIfat3IDWDH5fTXUWawvYF0aIpoOvKsT/wRVJT/3cj4sG8/0OkxRz/JiIGWp3vMGWZDRxJutOzB/ChiPhUEXmXpYw4y/kWEmtlx1nOr/RYc5yVy3HWG3GWz98RsdaLcQa+dsTXjr35nVZkrXEctd0dSbXb/yggr0K7vuQ8ypwm98WkwZmvIPUbPw5YQ7rb8oqadE9v8+u+f36ND8mP35Nf9//X5nznUHfXJu87iXQ3abB/+SS6bAxBg9eisDjL+RXdzazsqd9LjzXHWfmb42xI2q6Ms3z+wmPNcbbJ6+FrxyfS+tqxtfl2XKy1PYMWvGhH5zdr1wLyKrLrSyf0rX55DsjX5sd7Aa8HtmSENTnaUI79SX2svwo8COzT5vzOyK/9OcCr6449k9RU/mPgdUW9BmVvRcZZzq+ori+lx1nOr/RYc5yVvznO2v58S4+znG9hseY4G/Z18bVj+55vz8VZzq8jY62wjCbwwu0EXEVeG6NNeQz2851PqlGvAd7W5udV+pTUuRz9pPUsXl/y+9z//9u7exc7qjAA489rjOhCmggBFURBBa3UQogWRlnIvxCwSKEgBguxEzQgFulUkiKBFIogCBrERgkJoo0E/AfsVNCIpFBQ/AjKa3FGuAR2I3LPx8x9fjBFloU59948y5y5Z2aAH6n8DBvKXYyOT+/9G8CvwHvAcyu/c4ByduUrVm4ju+StRWfTfpq2Nkpn01i6t2ZnfTc724zOVsZRtTU72/W98dix7jg2prNpP8O2NuQ1cteKiFtyWgNbeT93UL4Kv5SZRxrsr9na6uuM43HgNOXi1XOZebXWvq4zjpsz848G+9lPeTDlCcofw5cpa5v3Am9SnkPyF9OtijPzm9pjGkGrzqZ9NWttlM6msXRvzc76srPN6GwaR/XW7GxnHjt67Ljm/YzZWs+Z9Igbbb+Ob762epexbANfAvt6fwaVXt8rrKwRp1wk+zGwNf37KOVMyyXKGa4jLPDWzCNtrVobqbNp/4ttq9xGDwAAAsJJREFUzc7G2+zMzuysyWfiseMAn0OF1zd8a93fpNE2Gi19Wdlfs7XV/2EsW73f/0qv6/mV0O6afnaIcienhyhnVL6b/hDfCRwDHug97qVvLVsbqbNpPItrzc7G3Oys/2ew5tdkZwNuHjv2/wwqvK5ZtDaLpZWttfg6vtctqTdRRNxNef7J/ZQzJtuZ+XVEnAKeonwVftz3vr3ardlZO3Y2LjtbDjsbl8eOyzKX1m7oufNR1Q5x2se/zwK5QLmI8mLv/wxLFBF7gMuUtctPA58DX0yBvktZjvB6Zp6JiBv7jXQz1W7Nztqws7HZ2TLY2dg8dlyOObXmRK6zzPweeBV4OCLu7T2epYiIgxGxJzP/zsw/gSuUW/UeozzQ8iLwA2Vt92GArPDgSo3BzuqwM62yszrsTNeytTrm2JoTuTF8CnwL/NR7IEsQEW9TLvg9HRGvTUsR3qJcoPpsZh6lXAz8EfA+8HNE3BoR0W3QasHO1sjOtAM7WyM70y5sbY3m2prXyA2i5S2ply4iHqWcOXmH8pDGAM4DtwE3ZeYL0++dB+4DHsnMK52Gq4bsbH3sTDuxs/WxM+3G1tZnrq05kdMiRcRjwBngRcrdhPZRnvmxHzicmRem37s9My93G6g0Y3Ym1WdnUhtzbM2JnBYrIp4ATlLuKvRhRDwI3AN8Avy+ctGwpP/JzqT67ExqY26tOZHTokXEIeAs8FJmftB5ONIi2ZlUn51JbcypNW9Pq0XLzM8i4hnKxat7gXOZebX3uKQlsTOpPjuT2phTa34jp40QEdvACeDJzPyl93ikJbIzqT47k9qYQ2tO5LQxImIrM3/rPQ5pyexMqs/OpDZGb82JnCRJkiTNjA8ElyRJkqSZcSInSZIkSTPjRE6SJEmSZsaJnCRJkiTNjBM5SZIkSZoZJ3KSJEmSNDNO5CRJkiRpZv4BOlQp/RBnuioAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "figure(figsize = (15,4))\n",
    "\n",
    "lw = 2\n",
    "label_fs = 16\n",
    "\n",
    "whis = [0,100]\n",
    "showfliers = False\n",
    "\n",
    "subplot(1,4,1)\n",
    "title(r'$\\ell^2$ Error', fontsize = 16)\n",
    "bp1 = boxplot([L2_err[key] for key in L2_err.keys()], positions = np.arange(6), \\\n",
    "              patch_artist=True, whis=whis, showfliers=showfliers)\n",
    "for box, whisker,cap in zip(bp1['boxes'], bp1['whiskers'], bp1['caps']):\n",
    "    box.set(color='blue', linewidth=lw)\n",
    "    box.set(facecolor = 'w' )\n",
    "yticks(fontsize = 12)\n",
    "xticks(np.arange(6), ['ARD','ARDr','ARDvi','M-ST','L-ST','MAP-ST'], fontsize = 12, rotation=45)\n",
    "\n",
    "for median in bp1['medians']:\n",
    "    median.set(color='k', linestyle = '--', linewidth=1.5)\n",
    "\n",
    "##########################################################################################\n",
    "##########################################################################################\n",
    "##########################################################################################\n",
    "\n",
    "subplot(1,4,2)\n",
    "title(r'$\\ell^1$ Error', fontsize = 16)\n",
    "bp2 = boxplot([L1_err[key] for key in L1_err.keys()], positions = np.arange(6), patch_artist=True, whis=whis, showfliers=showfliers)\n",
    "for box, whisker, cap in zip(bp2['boxes'], bp2['whiskers'], bp2['caps']):\n",
    "    box.set(color='red', linewidth=lw)\n",
    "    box.set(facecolor = 'w' )\n",
    "yticks(fontsize = 12)\n",
    "xticks(np.arange(6), ['ARD','ARDr','ARDvi','M-ST','L-ST','MAP-ST'], fontsize = 12, rotation=45)\n",
    "\n",
    "for median in bp2['medians']:\n",
    "    median.set(color='k', linestyle = '--', linewidth=1.5)\n",
    "\n",
    "##########################################################################################\n",
    "##########################################################################################\n",
    "##########################################################################################\n",
    "\n",
    "subplot(1,4,3)\n",
    "title(r'Added', fontsize = 16)\n",
    "bp3 = boxplot([Added[key] for key in Added.keys()], positions = np.arange(6), patch_artist=True, whis=whis, showfliers=showfliers)\n",
    "for box, whisker, cap, median in zip(bp3['boxes'], bp3['whiskers'], bp3['caps'], bp3['medians']):\n",
    "    box.set(color='green', linewidth=lw)\n",
    "    box.set(facecolor = 'w' )\n",
    "yticks(fontsize = 12)\n",
    "xticks(np.arange(6), ['ARD','ARDr','ARDvi','M-ST','L-ST','MAP-ST'], fontsize = 12, rotation=45)\n",
    "\n",
    "for median in bp3['medians']:\n",
    "    median.set(color='k', linestyle = '--', linewidth=1.5)\n",
    "\n",
    "##########################################################################################\n",
    "##########################################################################################\n",
    "##########################################################################################\n",
    "\n",
    "subplot(1,4,4)\n",
    "title(r'Missed', fontsize = 16)\n",
    "bp4 = boxplot([Missed[key] for key in Missed.keys()], positions = np.arange(6), patch_artist=True, whis=whis, showfliers=showfliers)\n",
    "for box, whisker, cap, median in zip(bp4['boxes'], bp4['whiskers'], bp4['caps'], bp4['medians']):\n",
    "    box.set(color='c', linewidth=lw)\n",
    "    box.set(facecolor = 'w' )\n",
    "yticks(fontsize = 12)\n",
    "xticks(np.arange(6), ['ARD','ARDr','ARDvi','M-ST','L-ST','MAP-ST'], fontsize = 12, rotation=45)\n",
    "\n",
    "for median in bp4['medians']:\n",
    "    median.set(color='k', linestyle = '--', linewidth=1.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}