L63 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",
    "from scipy.integrate import odeint\n",
    "import itertools\n",
    "from SparseARD import*\n",
    "from Oracle_Training import *\n",
    "import json\n",
    "np.random.seed(0)\n",
    "\n",
    "retrain = False\n",
    "\n",
    "n_trials = 10\n",
    "p = 5 # degree of polynomial library\n",
    "noise_percent = 1\n",
    "\n",
    "tol = 1e-8 # tolerance for ARD algorithm\n",
    "verbose = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lorenz_ode(X, t, params):\n",
    "    s,r,b = params\n",
    "    x,y,z = X\n",
    "    return [s*(y-x), x*(r-z)-y, x*y-b*z]\n",
    "\n",
    "n = 3\n",
    "m = 251\n",
    "dt = 0.05\n",
    "\n",
    "t = np.linspace(0,(m-1)*dt,m)\n",
    "\n",
    "params = (10,28,8.0/3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_maps, feature_descriptions = polynomial_feature_maps(n, p)\n",
    "d = len(feature_maps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if retrain:\n",
    "    \n",
    "    Xi = [{'ARD':[], 'ARDr' : [], 'ARDvi' : [], 'mSTSBL' : [], 'lSTSBL' : [], 'dSTSBL' : []} for dim in range(n)]\n",
    "    L1_err = [{'ARD':[], 'ARDr' : [], 'ARDvi' : [], 'mSTSBL' : [], 'lSTSBL' : [], 'dSTSBL' : []} for dim in range(n)]\n",
    "    L2_err = [{'ARD':[], 'ARDr' : [], 'ARDvi' : [], 'mSTSBL' : [], 'lSTSBL' : [], 'dSTSBL' : []} for dim in range(n)]\n",
    "    Added = [{'ARD':[], 'ARDr' : [], 'ARDvi' : [], 'mSTSBL' : [], 'lSTSBL' : [], 'dSTSBL' : []} for dim in range(n)]\n",
    "    Missed = [{'ARD':[], 'ARDr' : [], 'ARDvi' : [], 'mSTSBL' : [], 'lSTSBL' : [], 'dSTSBL' : []} for dim in range(n)]\n",
    "\n",
    "    for trial in range(n_trials):\n",
    "    \n",
    "        print('\\nTrial', trial+1, 'of', str(n_trials)+'.')\n",
    "\n",
    "        # get dataset\n",
    "        x0 = np.array([0,0,15]) + 5*np.random.randn(3)\n",
    "        X = odeint(lorenz_ode, x0, t, args=(params,), rtol = 1e-12, atol = 1e-12)\n",
    "\n",
    "        eps = 0.01*noise_percent*np.std(X)\n",
    "        measurement_noise = eps*np.random.randn(m,n)\n",
    "        X_obs = X + measurement_noise\n",
    "\n",
    "        Theta = np.hstack([f(X_obs) for f in feature_maps])\n",
    "\n",
    "        for dim in range(n):\n",
    "\n",
    "            print('Dimension', dim+1, 'of', str(n)+'.')\n",
    "\n",
    "            # True predictor\n",
    "            xi_true = np.zeros((len(feature_maps),1))\n",
    "            if dim == 0:\n",
    "                xi_true[2,0] = 10\n",
    "                xi_true[3,0] = -10\n",
    "            elif dim == 1:\n",
    "                xi_true[3,0] = 28\n",
    "                xi_true[7,0] = -1\n",
    "                xi_true[2,0] = -1\n",
    "            else:\n",
    "                xi_true[8,0] = 1\n",
    "                xi_true[1,0] = -8/3\n",
    "            xi = xi_true\n",
    "\n",
    "            y = FiniteDiff(X_obs, dt, 1, 6)[:,dim].reshape(m,1)\n",
    "\n",
    "            # Test each method\n",
    "            print('Training ARD.')\n",
    "            gamma_ARD,sigma2_ARD,mu_xi_ARD,_ = SBL([Theta,y], \\\n",
    "                                               sigma2=1, \\\n",
    "                                               estimate_sigma=True, \\\n",
    "                                               maxit=500, \\\n",
    "                                               verbose=False, \\\n",
    "                                               tol=tol)\n",
    "            ARD_results = gamma_ARD,sigma2_ARD,mu_xi_ARD\n",
    "\n",
    "            print('Training ARDr.')\n",
    "            lams = [np.linalg.norm(y)**2*l for l in [0,1e-1,1e0,1e1,1e2,1e3,1e4]]\n",
    "            gamma_ARDr, sigma2_ARDr, mu_xi_ARDr = Train_ARDr(Theta,y,\\\n",
    "                                                             xi_true,\\\n",
    "                                                             ARD_results=ARD_results,\\\n",
    "                                                             lams = lams,\\\n",
    "                                                             verbose=verbose)\n",
    "\n",
    "            print('Training ARDvi.')\n",
    "            alphas = [1,2,4,8,16,32,64,128,256]\n",
    "            gamma_ARDvi, sigma2_ARDvi, mu_xi_ARDvi = Train_ARDvi(Theta,y,\\\n",
    "                                                                 xi_true,\\\n",
    "                                                                 ARD_results=ARD_results, \\\n",
    "                                                                 alphas=alphas,\\\n",
    "                                                                 verbose=verbose)\n",
    "\n",
    "            print('Training STSBL.')\n",
    "            taus = [0,1e-5,1e-4,1e-3,1e-2,1e-1,1]\n",
    "            gamma_STSBL, sigma2_STSBL, mu_xi_STSBL = Train_M_STSBL(Theta,y,\\\n",
    "                                                                   xi_true,\\\n",
    "                                                                   ARD_results=ARD_results,\\\n",
    "                                                                   taus=taus,\\\n",
    "                                                                   verbose=verbose)\n",
    "\n",
    "            print('Training lSTSBL.')\n",
    "            taus = [np.inf,1e-8,1e-6,1e-4,1e-3,1e-2,1e-1,1,2,5,10]\n",
    "            gamma_lSTSBL, sigma2_lSTSBL, mu_xi_lSTSBL = Train_L_STSBL(Theta,y,\\\n",
    "                                                                      xi_true,\\\n",
    "                                                                      ARD_results=ARD_results, \\\n",
    "                                                                      taus=taus,\\\n",
    "                                                                      verbose=verbose)\n",
    "\n",
    "            print('Training dSTSBL.')\n",
    "            taus = [0,1e-2,1e-1,1,1e1,1e2,1e3]\n",
    "            gamma_dSTSBL, sigma2_dSTSBL, mu_xi_dSTSBL = Train_MAP_STSBL(Theta,y,\\\n",
    "                                                                      xi_true,\\\n",
    "                                                                      ARD_results=ARD_results,\\\n",
    "                                                                      taus=taus,\\\n",
    "                                                                      verbose=verbose)\n",
    "\n",
    "            # Record predictors and errors\n",
    "            Xi[dim]['ARD'].append(list(mu_xi_ARD.flatten()))\n",
    "            Xi[dim]['ARDr'].append(list(mu_xi_ARDr.flatten()))\n",
    "            Xi[dim]['ARDvi'].append(list(mu_xi_ARDvi.flatten()))\n",
    "            Xi[dim]['mSTSBL'].append(list(mu_xi_STSBL.flatten()))\n",
    "            Xi[dim]['lSTSBL'].append(list(mu_xi_lSTSBL.flatten()))\n",
    "            Xi[dim]['dSTSBL'].append(list(mu_xi_dSTSBL.flatten()))\n",
    "            \n",
    "            L1_err[dim]['ARD'].append(np.linalg.norm(xi - mu_xi_ARD, 1))\n",
    "            L1_err[dim]['ARDr'].append(np.linalg.norm(xi - mu_xi_ARDr, 1))\n",
    "            L1_err[dim]['ARDvi'].append(np.linalg.norm(xi - mu_xi_ARDvi, 1))\n",
    "            L1_err[dim]['mSTSBL'].append(np.linalg.norm(xi - mu_xi_STSBL, 1))\n",
    "            L1_err[dim]['lSTSBL'].append(np.linalg.norm(xi - mu_xi_lSTSBL, 1))\n",
    "            L1_err[dim]['dSTSBL'].append(np.linalg.norm(xi - mu_xi_dSTSBL, 1))\n",
    "\n",
    "            L2_err[dim]['ARD'].append(np.linalg.norm(xi - mu_xi_ARD, 2))\n",
    "            L2_err[dim]['ARDr'].append(np.linalg.norm(xi - mu_xi_ARDr, 2))\n",
    "            L2_err[dim]['ARDvi'].append(np.linalg.norm(xi - mu_xi_ARDvi, 2))\n",
    "            L2_err[dim]['mSTSBL'].append(np.linalg.norm(xi - mu_xi_STSBL, 2))\n",
    "            L2_err[dim]['lSTSBL'].append(np.linalg.norm(xi - mu_xi_lSTSBL, 2))\n",
    "            L2_err[dim]['dSTSBL'].append(np.linalg.norm(xi - mu_xi_dSTSBL, 2))\n",
    "\n",
    "            _,added,missed = sparsity_err(mu_xi_ARD,xi)\n",
    "            Added[dim]['ARD'].append(added)\n",
    "            Missed[dim]['ARD'].append(missed)\n",
    "\n",
    "            _,added,missed = sparsity_err(mu_xi_ARDr,xi)\n",
    "            Added[dim]['ARDr'].append(added)\n",
    "            Missed[dim]['ARDr'].append(missed)\n",
    "\n",
    "            _,added,missed = sparsity_err(mu_xi_ARDvi,xi)\n",
    "            Added[dim]['ARDvi'].append(added)\n",
    "            Missed[dim]['ARDvi'].append(missed)\n",
    "\n",
    "            _,added,missed = sparsity_err(mu_xi_STSBL,xi)\n",
    "            Added[dim]['mSTSBL'].append(added)\n",
    "            Missed[dim]['mSTSBL'].append(missed)\n",
    "\n",
    "            _,added,missed = sparsity_err(mu_xi_lSTSBL,xi)\n",
    "            Added[dim]['lSTSBL'].append(added)\n",
    "            Missed[dim]['lSTSBL'].append(missed)\n",
    "\n",
    "            _,added,missed = sparsity_err(mu_xi_dSTSBL,xi)\n",
    "            Added[dim]['dSTSBL'].append(added)\n",
    "            Missed[dim]['dSTSBL'].append(missed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "if retrain:\n",
    "    with open('./saved_results/L63_results.json','w') as f:\n",
    "        json.dump({'L2_err': L2_err, \\\n",
    "                   'L1_err': L1_err, \\\n",
    "                   'Added' : Added, \\\n",
    "                   'Missed' : Missed, \\\n",
    "                   'Xi' : Xi}, f)\n",
    "else:\n",
    "    with open('./saved_results/L63_results.json','r') as f:\n",
    "        results = json.load(f)\n",
    "    L2_err = results['L2_err']    \n",
    "    L1_err = results['L1_err']    \n",
    "    Added = results['Added']\n",
    "    Missed = results['Missed']\n",
    "    Xi = results['Xi']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA24AAAEtCAYAAABwAjOMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzde7ylZV3//9cbBkEYEJCRihomDERBQB0zv2pMXymVNA8YYSgM1Q+VL5Wl5glk8Fh9U6uvx0nkIEh5QMnMQyaTaGSOKdoUUip41sEDMQwHxc/vj/vesmax9t5rz16nvdfr+Xisx+x13de6r+tea31m3dd9Xfd1paqQJEmSJE2uXcZdAUmSJEnS3Gy4SZIkSdKEs+EmSZIkSRPOhpskSZIkTTgbbpIkSZI04Wy4SZIkSdKEs+EmSZIkSRPOhtsESLJHkvck+c8kn0nywSSHjLte0nKT5IVJPp/kR0meMO76SJIk9cuG2+R4Q1Xdt6qOAd4LvHncFZKWoX8Ejgc+Ou6KSJIkLYQNtzFI8uIkn0tyIkBV3VpVH+zI8i/AvD1uSdYnqVke3x9W/aWloDvOAKrqE1X1hQXuxziTJEljZ8NtxJIcD5wOfBz4lVmy/S5w+QJ2++vAQ7sexy2imtKS1mecLZRxJvUpyZvbixuvXsBrNiXZ1Ee+DUlqURXccX8zF2fWDGqf0jB0XUg8rMf2dR3bj2vTLkhy3cgrO4e2fhvGXY+laMW4KzCF/j/gr4DDgLtc+U/ygnbbIxewz89U1X8vpBJJdq+q2xa6bbH7lkZkzjjbScaZ1Ickd6e50AFwcpI/qqofjrNO0jJzE/A04Oyu9FPabXt3pL0U+IsR1UtDZo/bCCXZHfhl4G+B/wX8c9f25wAnAI+pqu0DLHdDe3XjyHbik23A2/vY9ugkVyW5JcmN7QQq9+l339I4zBdnQyzXOJMaTwT2Af4euBfw6PFWR1p2LgOemiQzCe0FkxOAd3VmrKovVNWnR1w/DYkNt9H6BeBHwO3A/sA/zWxI8ofAU4BfrqqF3jeza5IVXY9en+3lbZm/Brxmrm1JHg28D9gG/AbwTOBI4GNJDlrgvqVRmjXOFsk4k/pzKvA9YD1wC00vwA6SnJTkmiS3JdmS5Im9dpTkAUmuTHJrkq8lORtIj3wrkrygY59fT/KqJHt05TskyfuSbE+yNclfALsP4JilUXorcDDw8I60JwK70tVw6x4q2cbKS5N8oY2rG5J8LMnDO/L8ZpJPJ9nWXlD8XJKnd+332CT/mOSmJDe3FxWP7Mqza5KXJflGG3ObkhwxuLdh+jhUcrQeBnwKeCpwaVXdApDkp4FXAV8ErmgvoPywqtb2ud9reqS9D3hsV9pfVtVs3eU7bEuyua3PY2aGuCS5CrgWeDbwhwvYtzRKPeMMIMlZwDOAVcCRSV4LrK2qb/axX+NMmkeSn6K593NjVW1N8h7gSUn2q6rvtXmOA95GEz/PponHvwB2Az7fsa8DgI8A36RpDN4GPBdY3aPoi4HHAX9C08t+X5ohYmtoeiFIcjfgH4C7A/8H+DbwdOBJA3sDpNG4nmZ25KcBV7ZppwDvprkQOJfnAX8AvAj4DE3v+FqaC520DbiLgb+kibddgMOBfWd2kORXaS4kvo/mt3Zmv1cmOaqqvtKmbQBeCLwa+FBbzt/uxPGqZcNttI4AvgT8FvCLM4lV9VV6XEFcgCcCX+1K69Vr9+459vHjbUn2Ah4IvKLzvoSq+lKSjwPHLnDf0ij1jDOAqnoZ8LKd3K9xJs3vaTQnehe1zy+kGU3yG8Ab27RzaS6EPL6qfgSQ5D9pZlT+fMe+/gDYC3hUVX25zfcPNCetP5bkEe3+T62qmXI/nOS7wMVJjqmqz9A0/g4BHlpV/9K+9v3A5wZ07NIoXQS8KsnvAfvRXDB5TB+veyjwoa6LgO/t+PsXgO9X1bM60j7UtY+/AP6pqh4/k5DkCpoLkc8GnpVkP5oY3lhVz5nZT5I7gD/uo57qwaGSo/VTwC/RBMzn58u8AP9eVZu7Hr0mUfjGHPvo3LYfTUOyV/5v0l6VWcC+pVEyzqTxOQX4r6q6qn3+YeDrbTpJdgUeDLxzptEGzVIdwHVd+3oo8C8zjbY2383seJIJzT10twPv6hzKzJ0nmzMXcB4KfGWm0dbu70d4v6iWpnfQDPN9HHAyze/GP/bxuk8Cxyd5eZKHtz3R3dv3S3Jxkscm2bdzY5JDgXsDl3TF23bgKu6Mt/vTXHjpjq+/7v8Q1c2G22jtSzMk5AVjKn+u6ZM7t32vff4TPfL9BPCdBe5bGiXjTBqDJA8G7gdclmTf9oRvb5qJFB6aZvryA2iGRH6rxy66036yz3z3Au5GM0TsBx2Pb7fb77nA/UkTr6puAt5D08t9CnBJ58WQObwCOIfmXukrge8kOb8dmkxV/RPNrLA/QzPKY2uSDyc5qn39vdp/z2PHePsBza0DnfEGd40v420RbLiN1i7AX7RDIydWe0XzU8Cvt1dHAUhyMM0sfYOa7EEaBuNMGo9T23+fR3NhYuZxZpt+CnADzQnegT1e3532jT7zfQe4laYnr9fjTQvcn7RUXAT8Kk3v1kXz5AWgqn5QVX9SVfenaVz9Ac19oK/ryPPOqjqWZmTIE9t8H2gn5Jq5qPgCesfb49rtMyNEuuPLeFsE73EbkSSn0gTWv7cnaa8BXltV1w5g98fMXCnpsnkRa+ecTXPT6d8leT2wkua+hBtpJlKRJo5xJo1HO9zqJOATwPN7ZHkNd6479UngyUk2dNzj9hCaiUQ671+7Cnhukp+ZmeygvTf0cezoAzSNxXtU1VxDxa4CTkvyCx33uO0CnLiQY5UmyD/QDEX8flVtWeiL24m53pzkeJoZjbu3b6P5fTqE5r62e9Lch3odcERVzXWv2meBm2ni6yMd6ScttJ66kw23EUiyJ02382NpZrn6d+CvBnQyCc04515W0VzdXLCq+kA7a9A5NP8p3A5sAv6oqr6+M/uUhsk4k8ZqZojUs6tqU/fGJG8C3gCso/m+fwh4T5u+iuaCRffsrq8BzqCZ0GADd84qeUtnpqralORS4J1JXg38K82SIGuA44Hntf8PXEjTqLwsyQtphlI+g2ZWPWnJqao7aCb/6VuSy4GrgX+j6RF/AM19om9qt7+EplfsCpr7U38a+D3gM1W1tc3zf4DL2ws2b6f5DTyQZrTIl6vq1VX1/SSvAV6U5CaamH8w8NuLOugplypvmZAkSTuvPRn8JeAnqmp7j+33oBk69faqWp/kKTRTha8B/hs4C/h9gKpa1/G6B9Jc6X8wzRCtN9JcdH5xVXUuPrwL8Ls0s8neh6aRdx3wQeDlVXVjm+8Q4LU0DcibaZYl+I92vz9bVdct+s2QhiTJeuB84NBZJsciyTqaRtcvV9WHk1wArKuqNe32Z9Nc5DwU2BP4MnApTZz8oL2Y+HvAUTSTZH2bptF1ducFxSQPpVlS4GE0S2x8k2Zm2L+YmZyoHflyLvA7wD1oeuTPALYA51bVhgG8LVPFhpskSZIkTTgnJ5EkSZKkCWfDTZIkSZImnA03SZIkSZpwNtwkSZIkacKNZTmAAw44oNasWTOOoqWh+dSnPnVDVa0adz06GWtajiYt1owzLUfGmTR8C42zsTTc1qxZw+bNm8dRtDQ0Sa6fP9doGWtajiYt1owzLUfGmTR8C40zh0pKkiRJ0oSz4SZJkiRJE86GmyRJkiRNOBtukiRJkjThbLhJkiRJ0oSz4SZJkqQ5Jdk9yXlJrk9yU5JPJ3nMHPn/IMk3k9yY5C1Jdu/YtibJFUm2J7kmyXGjOQppabPhJkmSpPmsAL4CHAvcAzgbeHuSNd0ZkzwKeD7wSGANcAhwbkeWS4FPA/cEXgS8M8nErBknTSobbpIkSZpTVd1cVRuq6rqq+lFV/R3wJeBBPbKfCpxXVVuq6nvAS4H1AEkOAx4InFNVt1TVu4DPASeM5ECkJWwsC3BrciWZc3tVjagm0vI1X5yBsbZc+dlruUhyIHAYsKXH5iOAyzueXw0cmOSe7bYvVtVNXduP6FHG6cDpAKtXr15o/ebcPqw4G1eMj/P/Ft/r0ZQLNtzUpfOLlsQTCGkIuuPKWJsefvZaDpLsBlwCXFhV1/TIshK4seP5zN9799g2s/2g7p1U1UZgI8DatWsXFCjjirVpK3ecZU9bueBQSUmSJPUpyS7AW4HbgTNnybYN2Kfj+czfN/XYNrP9JiTNyYabJEmS5pVmjNh5wIHACVX1g1mybgGO7nh+NPCtqvpOu+2QJHt3be815FJSBxtukiRJ6scbgPsCj6uqW+bIdxHw20nul2Q/4CzgAoCquhb4DHBOkj2SPBE4CnjXUGsuLQM23CRJkjSnJAcDTweOAb6ZZFv7ODnJ6vbv1QBV9QHgT4ErgOvbxzkduzsJWAt8D/hj4MlVtXWEhyMtSU5OIkmSpDlV1fXAXNPprezK/2rg1bPs6zpg3aDqJk0Le9ykCZJkU5JbO65kfr5j228muT7JzUnek2T/cdZVkiRJo7OghluSQ9uTyos70jyZlAbrzKpa2T7uA5DkCOBNwNNobgrfDrx+jHWUJEnSCC20x+11wCdnnngyKY3MycB7q+qjVbUNOBt4UtesXJIkSVqm+m64JTkJ+D7wjx3JnkxKg/fKJDck+XiSdW3aEcDVMxmq6gs0a+gcNob6SZIkacT6argl2Qd4CfDsrk19n0wmOT3J5iSbt2514iBpFs8DDgEOAjYC701yb5qbvm/synsjcJeLJMaaJEnS8tNvj9tLgfOq6itd6X2fTFbVxqpaW1VrV61atfCaSlOgqj5RVTdV1W1VdSHwceB4YBuwT1f2fYCbeuzDWNNUSXJme7HitiQXdKSf3DHRz7Yk25NUkgfNsp9ZJweSJGnc5m24JTkGOA54TY/NfZ9MStopRTP98hbg6JnEJIcAuwPXjqle0iT5OvAy4C2diVV1ScdEPyuBM4AvAv82x77uMjmQJEmToJ913NYBa4AvJ4Gml23XJPcDPoAnk9JAJNkXeAjwT8APgd8AfhF4Fk2sXpXkETQnnS8BLqsqL5Jo6lXVZQBJ1gI/PUfWU4GLqqpGUjFJkgaon6GSG4F7A8e0jzcC7wMeBVwCPC7JI5LshSeT0mLsRtNrsBW4Afhd4AlV9fmq2gI8gybmvk0zHPmMcVVUWmqSHExzIeSiebL2mhxIkqSxm7fHraq200zzD0CSbcCtVbUV2Jpk5mTynsCHgdOGVFdpWWtj6sFzbH8b8LbR1UhaVk4BrqyqL82R53nAf9BMsnUSzeRAx7QTb+0gyenA6QCrV68eQnUlSdrRQtdxo6o2VNVTO56/rapWV9VeVfX4qvruYKsoSdKinQJcOFeGOSYH6pXXSYAkSSO14IabJElLSZKHAT8FvHOBL52ZHEiSpLGz4SZJWtKSrEiyB7ArzeRZeyTpvBXgVOBdc91/nWTfJI+aeW2Sk2nuifvgcGsvSVJ/bLhJkpa6s4BbgOcDT23/PgugbdCdSI9hkklemOT97dNZJwcaeu0lSepDP8sBSJI0sapqA7Bhlm23AvvOsu0VHX/POTnQctAu6TMrV0nQXJKcCawH7g9cWlXrZ8n3RpoLKDN2A26vqr3b7ZuAX6BZ9gbga66ZKPXHhpskSVOgu2GWxMaaFmJmoftHAXefLVNVPYNm+RoAklwA/Kgr25lV9eYh1FFa1my4SZIkaU4LWOj+x9o1fk8AHjvEqklTw3vcJEmSNAwn0Nw3+tGudBe6l3aCDTdJkiQNw6nARbXjmNznAYcABwEbaRa6v3evFyc5PcnmJJu3bt06/NpKE86GmyRJkgYqyc8AxwIXdaa70L2082y4SZIkadBOAf65qr44Tz4Xupf6ZMNNkqZEknkfktRLHwvddzsFuKBrHy50Ly2CDTdJmhJVtcNjtjRJ6qHnQvdJVifZlmT1TMYkD6WZefIdXftwoXtpEVwOQJIkSXOaa6F7YGVX3quAvXrsY9kvdC8Nkz1ukiRJkjThbLhJkiRJ0oSz4SZJkiRJE86GmyRJkiRNOBtukiRJkjThnFVSEkBfa3gNY7r4cZWr0Zrvc/YzliRpbjbcJAF3PXFOMpKT6XGVq9Hq/Ez9jCVJWjiHSkqSlrQkZybZnOS2JBd0pK9JUu3iwDOPs+fYz5okVyTZnuSaJMeN5AAkSeqDPW6SpKXu68DLgEcBd++xfd+q+mEf+7kUuAo4vn28M8mh7aLBkiSNlT1ukqQlraouq6r3AN/Z2X0kOQx4IHBOVd1SVe8CPgecMKBqSpK0KDbcJEnL3fVJvprk/CQHzJLnCOCLVXVTR9rVbbokSWNnw02StFzdADwYOBh4ELA3cMkseVcCN3al3di+5i6SnN7eV7d561ZHUkqShs+GmyRpWaqqbVW1uap+WFXfAs4EfiXJPj2ybwO60/cBbuqRl6raWFVrq2rtqlWrBltxSZJ6sOEmSZoWM2sQ9FpUbgtwSJLOHraj23RJksbOhps0YZIcmuTWJBd3pP1mkuuT3JzkPUn2H2cdpUmSZEWSPYBdgV2T7NGmPSTJfZLskuSewF8Cm6qqe0gkVXUt8BngnPb1TwSOAt41ymORJGk2NtykyfM64JMzT5IcAbwJeBpwILAdeP14qra8JJnzoSXjLOAW4PnAU9u/zwIOAT5AM9zx34HbgKfMvCjJG5O8sWM/JwFrge8Bfww82aUApMZs6yX2yLc+yR1d6yeu69jueonSTnIdN2mCJDkJ+D7wz8DPtcknA++tqo+2ec4G/jPJ3l0z4GmBqurHfyfZ4bmWjqraAGyYZfOlc7zuGV3PrwPWDaha0nIz33qJna6qqofPss31EqWdZI+bNCHaCRNeAjy7a9MRNNOSA1BVXwBuBw4bXe0kSdPM9RKl8bPhJk2OlwLnVdVXutKdplyStJQ8IMkNSa5NcnaSmRFerpcoLYINN02E+e41Wu73GyU5BjgOeE2PzU5TLklaKj4KHAnci6Yn7SnAc9ttXoiUFsF73DQRuu8tmsL7jdYBa4Avt43UlTSz492PZnKFo2cyJjkE2B24duS1lCRpDlX1xY6nn0vyEpqG2yvZiQuRwEaAtWvXTtVJgdSLDTdpMmwE/rrj+XNoGnLPpLlqeVWSRwD/RnMf3GVOTCJJWgKKO9dO/PF6iR2/YUcDbxtLzaQlxqGS0gSoqu1V9c2ZB81VyVuramtVbQGeAVwCfJtmSMkZY6yuJGnKzLZeYo98j0lyYPv34cDZwOXgeonSYi2JHrf57m8a1pC6fu6rmrLhfMvOpH7G7fTmnc/fhlckJUnjcxZwTsfzpwLnJnkL8B/A/arqy8AjgQuSrAS+BVwMvKLjdScBF9Csl/hlXC9R6ltfDbckF9ME4l7AN4E/rao3t9seSbNg8GrgE8D6qrp+kJUc11pL3ne1/PkZS5I0v3nWS1zZke85NMP9Z9vPdbheorRT+h0q+UpgTVXtA/wa8LIkD0pyAHAZTTf4/sBm4G+GUlNJkiRJmlJ99bi199j8+Gn7uDfwIGBLVb0DIMkG4IYkh1fVNQOuqyRJkiRNpb4nJ0ny+iTbgWuAbwB/T7Ng4tUzearqZuALuJCiJEmSJA1M3w23qjqDZja7R9AMj7yNBSyk6CKKkiRJkrRzFrQcQFXdUVUfA36aZn2pvhdSrKqNVbW2qtauWrVqZ+s7FZLM+5AkSZI0PXZ2HbcVNPe4baFZOBGAJHt1pGsnVdUOj9nSJEmSJE2HeRtuSe6V5KQkK5PsmuRRwFOAjwDvBo5MckK7KOOLgc86MYkkSZIkDU4/PW5FMyzyqzSLJf4Z8KyqurxdMPEE4OXttofQLKwoSZIkSRqQeRtuVbW1qo6tqn2rap+qun9V/VXH9g9X1eFVdfeqWtcurChJ0kgkObOd/Oq2JBd0pP9Ckn9I8t0kW5O8I8lPzrGfTUluTbKtfXx+JAcgSVIfdvYeN0mSJsXXgZcBb+lK3w/YCKwBDqaZOOv8efZ1ZlWtbB/3GXRFJUnaWX0twC1J0qSqqssAkqylmfV4Jv39nfmSvBb4p9HWTpKkwbDHTZI0LX6R+Wc9fmWSG5J8PMm6EdRJkqS+2HCTJC17SY6imfn4uXNkex5wCHAQzRDL9ya59yz7O729r27z1q1bB15fSZK62XCTJC1rSX4OeD/w+1V15Wz5quoTVXVTVd1WVRcCHweOnyXvxqpaW1VrV61aNZyKS5LUwYabJGnZSnIw8GHgpVX11gW+vIAMvlaSJC2cDTdJ0pKWZEWSPYBdgV2T7NGmHQR8BHhdVb1xnn3sm+RRHa89meaeuA8O/wikyTfbshs98p2a5FNJ/ifJV5P8aZIVHdtddkPaSTbcJElL3VnALcDzgae2f58F/A7NPWvndJwkbpt5UZIXJpmZeXI3miUFtgI3AL8LPKGqPKmUGrMtu9FtT+BZwAHAQ4BHAs/pyuOyG9JOcDkASdKSVlUbgA2zbD53jte9ouPvrcCDB1oxaRmZbdmNHvne0PH0a0kuAX5pyNWTpoI9bpIkSRqWXstw9LXshrO3Sjuy4SZJkqSBS3IasBb4s47kvpfdcPZWaUc23CRJkjRQSZ4A/DHwmKq6YSZ9IctuSNqR97hJkiRpYJI8Gvgr4Fer6nPzZHfZDalP9rhJkiRpTrMtu9Ej3/8GLgFOqKp/7drmshvSIthwkyRJ0nx6LruRZHW71MbqNt/ZwD2Av+9YhsNlN6QBcKikJEmS5jTPshsrO/LNOvW/y25Ii2OPmyRJkiRNOBtukiRJkjThbLhJEyLJxUm+keR/klyb5Hc6tj0yyTVJtie5IsnB46yrJEmSRsuGmzQ5Xgmsqap9gF8DXpbkQUkOAC6jueF7f2Az8Dfjq6YkSZJGzclJpAlRVVs6n7aPewMPArZU1TsAkmwAbkhyeFVdM/KKSpIkaeTscZMmSJLXJ9kOXAN8A/h74Ajg6pk8VXUz8IU2XZIkSVPAhps0QarqDGBv4BE0wyNvo5lm+caurDe2+e4iyelJNifZvHXr1mFWV5IkSSNiw02aMFV1R1V9DPhp4JnANmCfrmz7ADfN8vqNVbW2qtauWrVquJWVJEnSSNhwkybXCpp73LYAR88kJtmrI12SJElTwIabNAGS3CvJSUlWJtk1yaOApwAfAd4NHJnkhCR7AC8GPuvEJFIjyZnt8ODbklzQta3vpTSSrGnzbG9fc9zQKy9JUp9suEmToWiGRX4V+B7wZ8CzquryqtoKnAC8vN32EOCkcVVUmkBfB14GvKUzcSeW0rgU+DRwT+BFwDuTON5YkjQRXA5AmgBt4+zYObZ/GDh8dDWSlo6qugwgyVqae0NnPIk+l9JIchjwQOBXquoW4F1JnkVz0eSNwz8KSZLmZo+bJGm5WshSGkcAX6yqzkl/rp4lryRJI2fDTZK0XC1kKQ2X3ZAkTTQbbpKk5WohS2m47IYkaaJ5j5skabnaApw682SepTS2AIck2btjuOTRwNsWU4Gcm7kzbJg/T51Ti6mCNBBJzgTWA/cHLq2q9XPk/QPgecDdgXcBz6yq29pta4DzaSba+jJwZnsft6R52OMmSVrSkqxol8rYFdg1yR5JVrCApTSq6lrgM8A57eufCBxFc9IpaZbZW7u1y9k8H3gksAY4BDi3I4uzt0o7yR43SdJSdxZwTsfzpwLnVtWGJCcArwUuBj5Bx1IaSd4IUFXPaJNOAi6gWXbjy8CT2xlfF21nes3m7a2TRmiO2Vu7nQqcV1Vb2vwvBS4Bnj+s2VuzadP8ma64Ys58tW7dkim3r7LnKXdny86mTfCsZ911w7p18IQnwK23wtFHk2OO2XH7ox/dPJZYuX2V/f7337XcAZXdzYabJGlJq6oNwIZZts26lEZHg23m+XXAuoFWTpo+RwCXdzy/GjgwyT1Z4OytSU4HTgdYvXr1cGorLSE23CRJkjQo3TO0zvy9d49tM9sP6rWjqtoIbARYu3btvN3WQ+s1m9Byx1b2n//53OW2vUzLptxxl93Be9wkSZI0KN0ztM78fVOPbTPbe87eKmlHNtwkSZI0KFtoZmSdcTTwrar6Dh2zt3Zt7zXTq6Qu8zbckuye5Lwk1ye5KcmnkzymY/sjk1yTZHuSK5IcPNwqS5IkaZTmmL2120XAbye5X5L9aCYPugCcvVVarH563FYAXwGOBe4BnA28PcmaJAcAl7Vp+wObgb8ZUl0lSZI0HmcBt9BM9f/U9u+zkqxOsi3JaoCq+gDwp8AVwPXto3PW15OAtTSzt/4xA5y9VVru5p2cpKpuZsfZuv4uyZeAB9GswbGlqt4BkGQDcEOSw3utkyNJkqSlZ67ZW2kmHenM+2rg1bPs5zqcvVXaKQueVTLJgcBhNOORn0kzjSvQNPKSfIFmWtdrul7nlK6SJI1AX2vAbZg7386sPSdJGp4FTU6SZDeaRRQvbHvUZpvWde/u11bVxqpaW1VrV61atbP1lSRJkqSp03ePW5JdgLcCtwNntslO6ypJ0oTamV6zvnrrJEkj11ePW5IA5wEHAidU1Q/aTTtM+ZpkL+DeOK2rJEmSJA1Mv0Ml3wDcF3hcVd3Skf5u4MgkJ7RTxL4Y+KwTk0iSJEnS4PSzjtvBwNOBY4BvtlO+bktycjt96wnAy2mmdX0IzTSvkiRJkqQB6Wc5gOuBWQe8V9WHgcMHWSlJkiRJ0p0WNKukJEmSJGn0bLhJkiRJ0oSz4SZJkiRJE86GmyRp2eqYUGvmcUeS/zdL3vXt9s7860ZcZUmSeup7AW5Jkpaaqlo583e71ui3gHfM8ZKrqurhQ6+YJEkLZMNNmkaZdaLYH6v58lUNrDrSiDwZ+DZw5bgrIknSQtlwkyRNi1OBi6rmvOrwgCQ3AN8F3gq8sqp+uKhSz4d1V6zbIenEE0/kjDPOYPv27Rx//PF3ecn69esXVaQkafnxHjdpmlUt/LFYyZyP6iOPJtwEfsZJVgPHAhfOke2jwJHAvYATgKcAz51lf6cn2Zxk89atWwddXWkiJdk/ybuT3Jzk+iS/OUu+93fdK3p7ks91bL8uyS0d2z80uqOQli573KQJkGR34PXAccD+wH8DL6yq97fbH5lJUH0AACAASURBVAm8DlgNfAJYX1XXj6m60lJ0CvCxqvrSbBmq6osdTz+X5CU0DbdX9si7EdgIsHbt2rmvaJwGm87Z1HPTnnvuyaZNvbeddu5pc+5WGoPXAbcDBwLHAO9LcnVVbenMVFWP6XyeZBPwka59Pa6qPjzEukrLjj1u0mRYAXyFpkfgHsDZwNuTrElyAHBZm7Y/sBn4m3FVdGDG0dun0Zqsz/gU5u5t66UAu3glfjy5zwnA2VW1rao+Bvwt8LR5XrcGeATN0GNJi2DDTZoAVXVzVW2oquuq6kdV9XfAl4AHAU8CtlTVO6rqVmADcHSSw8dYZWnJSPK/gIOYezZJkjwmyYHt34fTXCy5fPg1lJaEw4A7qurajrSrgSPmed0pwJU9ersvSbI1yYeSHD3IikrLlQ03aQK1J4+HAVtofhSvntlWVTcDX2D+H0tJjVOBy6rqps7EJKvb+2tWt0mPBD6b5Gbg72l6ul8x2qpKE2slcGNX2o3A3vO87hTggq60k4E1wMHAFcAHk+zb/ULvJZV25D1u0oRJshtwCXBhVV2TZCXQ/Ys1649lktOB0wFWr17dK4s0Varq6bOkf5nmZHTm+XOA54yqXtKgpZ+lXnZ+SPI2YJ+utH2Am3rknanPw4GfAN7ZVYePdzx9ZZJTaYZTvrcrX//3kkpTwB43aYIk2YXmPoDbgTPb5AX9WFbVxqpaW1VrV61aNbS6SpImS1Xt8JgtbSddC6xIcmhH2tE0I0NmM9PbvW2+quP9pNK8Jq7hNv8M0eVM4VqW0lwqPY9mtq4TquoH7aYtND+OM/n2Au7N3D+WkiQNTDtM/zLgJUn2SvIw4PHMMulIkrsDv07XMMl2iPLDktwtyR5JngscAHy8x24kdZi4htu4zNcY7KfRKC3SG4D70kyRfEtH+ruBI5OckGQP4MXAZ6vqmnFUUpI0tc4A7g58G7gUeGZVbUnyiCTdvWpPoBnWf0VX+t40v3ffA74GPBp4TFV9Z6g1l5aBib3HbWd68208aT7zf0dq3jzDmLE8ycHA04HbgG923Kfw9Kq6JMkJwGuBi2nWcTtp8LWQtKycD+uuWLdD0oknnsgZZ5zB9u3bOf744+/ykvXr14+mblqSquq7NA2y7vQr6bhftE27lKZx1513C3DUsOooLWcT23AbFxuMGod2Me1Zv0ntIqVO/y9JkjSlbLhpKtlAl7TsnQabztnUc9Oee+7Jpk29t5127mnDq5Mkaad5j5skSZIkTTgbbpIkSZI04Wy4SZIkSdKEs+EmSZIkSRPOhpskSZIkTTgbbpIkSZI04Wy4SZIkSdKEs+EmSZIkSRPOhpskSZIkTbgV466AJC1rybxZar58VQOrjiRJWprscZMkLWtJNiW5Ncm29vH5WfIlyZ8k+U77+NOkj5a3JEkjYI+bJI3CzvSa2WYYpDOr6s3z5DkdeAJwNE1H6D8AXwTeOOS6SZI0L3vcJElqnAq8qqq+WlVfA14FrB9vlSRJatjjJknL1Tw9dvPeWwfL6f66Vyb5Y+DzwIuqalOPPEcAV3c8v7pNW7ScO57e03GVq+Upyf7AecCvADcAL6iqt/XItwF4EXBbR/JRVfXFdvsx7X7uC/wn8NtV9Znh1l5a+uxxkyQtd88DDgEOAjYC701y7x75VgI3djy/EVjZ6z63JKcn2Zxk89atW4dRZ2kSvQ64HTgQOBl4Q5LZLm78TVWt7HjMNNruBlwOXAzsB1wIXN6mS5qDPW6StNxN+f11VfWJjqcXJnkKcDzw/7qybgP26Xi+D7Ct6q5vYFVtpGkEsnbt2lnf4Dpn7vc+CT12v2jzlTvMsrU8JdkLOAE4sqq2AR9L8rfA04DnL2BX62jOP/+8ja2/TPIc4H8DHxhsraXlxYbbFJv/vKymaBSVpClSQK//3bbQTEzyr+3zo9s0SXAYcEdVXduRdjVw7Cz5H5fku8A3gNdW1Rva9COAz3ZdEPlsm77zDbdnPYt1++67Q9KJJ57IGWecwfbt2zn++OPv8pL169fDmjU7XeRYyx132RoLh0pKkpatJPsmeVSSPZKsSHIy8IvAB3tkvwj4wyQHJfkp4NnABSOsrjTJuocS0z7fu0fet9Pcv7YK+P+AF7c93Qvaj0OSpR311eOW5EyambXuD1xaVes7tj2SZszzauATwPqqun7gNdXQTPkoKknL227Ay4DDgTuAa4AnVNXnkzwCeH9VrWzzvonmXrjPtc/f3KZJuutQYtrnN3VnrKr/6Hj6z0n+AngycOkC99PXkGQA/vzP2bRuXc9Ne+65J5s2beq57bRZ0vs2rnLHXbbGot8et6/T/PC9pTMxyQHAZcDZwP7AZuBvBllBLT/J/I+ZYZqzb5ek+VXV1qp6cFXtXVX7VtUvVNU/tNuu7Gi0UY0/qqr928cf9bq/TZpS1wIrkhzakdbvcOLO4clbgKO6Jv05qs/9SFOtr4ZbVV1WVe8BvtO16UnAlqp6R1XdCmwAjk5y+GCrKUmSpHGpqptpLta/JMleSR4GPB54a3feJI9Psl8aPw/8Hs1MkgCbaHq/fy/J7u2oLoCPDP0gpCVusfe47bDmTRvUX6DHujeOU1a3qoU/JEnS2JwB3B34Ns2wx2dW1ZYkj0iyrSPfScB/0wx/vAj4k6q6EKCqbgeeAJwCfB/4LZrhy7eP7jCkpWmxs0quBLpbYT1vMF3QOGVJkiRNlKr6Lk2jqzv9SppzwpnnT+nO05X/08CDBl5BaZlbbI9b3zeYSpIkSZJ2zmIbbjNr3gA/Xpzx3niDqSRJd5Fkh8dsaVrauj9TP2NJg9BXw61d+2YPYFdg15n1cIB3A0cmOaHd/mKaRRWvGV6VJUlamqpq3oeWvu7P089Y0iD02+N2FnAL8Hzgqe3fZ1XVVuAE4OXA94CH0NyQKkmSJEkakH6XA9hQVel6bGi3fbiqDq+qu1fVuqq6bpgVlparJGe2M6/eluSCrm2PTHJNku1Jrkhy8JiqKUmSpDFY7D1ukgbHhe4lSZLU02KXA5A0IFV1GUCStcBPd2z68UL37fYNwA1JDvd+UkmSpOlgj5s0+fpe6F6SJEnLkw03afKtpFnYvlPPhe4Bkpze3iu3eevWrUOvnCRJkobPhps0+Ra00H1VbayqtVW1dtWqVUOvnCRJkobPhps0+VzoXpIkaco5OYk0IdpF7VfQsdA98EOahe7/b5ITgPfhQveSNPVybubPtGHufHWOi4FLS4k9btLkcKF7SUOTZIdHd5okabLZ4yZNiHZR+w2zbPswcPgo6yMtB0l2B14PHEezDuJ/Ay+sqvf3yLseOI/mosmMx1bVpuHXdPiq7F1Zjnam16yv3jpJE8eGmyRpOVsBfAU4FvgycDzw9iT3r6rreuS/qqoePsL6SZLUF4dKSpKWraq6uao2VNV1VfWjqvo74EvAg8ZdN2mpSbJ/kncnuTnJ9Ul+c5Z8z03y70luSvKlJM/t2n5dkluSbGsfHxrNEUhLmz1ukqSpkeRA4DBmn5X1AUluAL4LvBV4ZVX9cFT1kybc64DbgQOBY4D3Jbm6qrrjKcApwGdpZkH+UJKvVNVfd+R5XHsbgKQ+2XCTJE2FJLsBlwAXzjIr60eBI4HrgSOAv6GZ2fWVPfZ1OnA6wOrVq4dVZWlitEvRnAAcWVXbgI8l+VvgaTSTav1YVf1px9PPJ7kceBjQ2XCTFiybNk1l2TMcKilJWvaS7ELTg3Y7cGavPFX1xar6Ujuk8nPAS4Anz5LXhe41bQ4D7qiqazvSrqa5yDGrNFOWPoK79nJfkmRrkg8lObrHS0lyepLNSTZv3bp1MXWXlgV73CRJy1p74ngezfCu46vqB32+tGiGfEmClcCNXWk3AnvP87oNNB0F53eknQz8G018/T7wwSSHV9X3O19YVRuBjQBr1651WtQpVuvWzZsnyVBmz52v7GGV24s9bpKk5e4NwH1p7qm5ZbZMSR7T3gNHksOBs4HLR1NFaeJtA/bpStsHuGm2FyQ5k+Zet1+tqttm0qvq41V1S1Vtr6pXAt+n6ZWTNAd73CRJy1aSg4GnA7cB3+xYaPrpwJXAfwD3q6ovA48ELkiyEvgWcDHwipFXWkvKvGuibZg/z86sxTYG1wIrkhxaVf/Vph3NLBP9JPktmnvffrGqvjrPvu3dlvpgw02StGxV1fXMfUK4siPvc4DnDL1SU6ajsTwrFweffFV1c5LLgJck+R2aWSUfD/yv7rxJTqa56PFLVfXFrm2rgZ8BPkkz8ut3gQOAjw/3CKSlz4abJEkamu5G2SjvBxmlnek1m7e3bvKcAbwF+DbwHeCZVbUlySOA91fVzIWQlwH3BD7Z0XC/uKqeQXNP3Btolgm4FfgM8Jiq+s7oDkNammy4SZKkgeirIbJh7nxLZNjgVKqq7wJP6JF+JTv2Xv/sHPvYAhw1lApKy5yTk0iSJEnShLPHTZI0cOsAuqZQPvHEEznjjDPYvn07xx9//F1es379etaPoG4avikZNihJI2XDTVNoXff5ZF8nlHhKOTDrwJN6SZKkBbDhJkkauE0Amzb13LbnnnuyaZZtnHbacCokSdISZ8NtB/bETIdNs51PznlC6fnk4GwCT+olSZIWwIbb1BtXY9VG8ritA4crjsg6GNt7Pc6y6WP9LkmS1J8JbbiN66Tenphp4fnkFJrnQ68+8rAM156SJElLw4Q23KbL/I2IGuL55LgaqzaSx20TOFxxRDbB2N7rsZQ9z39Iy3UBZrXOh5zf9aN1BPDzwO3AJT1ecwzwgKHXbPk5H9ZdsW6HpP4udEtaiia04eZJvYZjvnNFTyiXt3UAxx67Q1pfwwb9z0VactLH0Ar/v5e0lExow20ares+n+zrhPK009aPpnqSJM2jzik4Z+48ecVoLpB1l7EsL8ydBpvO2dRz05wXus/1YpS0FNlwkzQ1NoHDQyVJ0pJkw21iODxUYzDOWVqcIUaSJKlvE9twG9c5neeSo+N7rbGYxi/eNB6zJEnLzMQ23KaR51YamT7u8xja/SDOOChJkrRgu4y7At2q5n5A+sgz+HL7KXupShb+WIxpfq81Rj2+SIE5HwP94o060HrVv8fxD+U/1AmTZP8k705yc5Lrk/zmLPmS5E+SfKd9/Gn6mZpQmhKDiqUkxyT5VJLt7b/HjO4opKVr4hpu06hXo2W+xzI8t5JGrqrmfGjZeB3NCmIHAicDb0hyRI98pwNPAI4GjgIeCzx9VJWUloBFx1KSuwGXAxcD+wEXApe36ZLmYMNtAs13MjmoE8px9W5KU6WPLuZ5e76005LsBZwAnF1V26rqY8DfAk/rkf1U4FVV9dWq+hrwKmD9yCorTbABxtI6mlt1/ryqbquqv6S5Kv2/h3wI0pI3kHvckuwPnAf8CnAD8IKqetsg9i2pYZxJO+Uw4I6qurYj7Wrg2B55j2i3debr1ZugCZU1PUa2HgH8PPDCWbYfAzygXYNuZ50POb9r3zPl3g5c0uM1bbmLlXNHNpp3ULF0BPDZ2vEq9Gfb9A8spoKZbXruIRtXueMse5zHPM0GNTlJZ9f5McD7klxdVVsGsfPuWwy6nw9rSFOvWxtGVfa4TNt7vcQ+45HGWa80P4PBmMb3elz/twArgRu70m4E9u4j743AyiTpOskkyek0w8FYvXr14Gq7DE1jjC9TA4mlhezHONNcpvG3dNENt46u8yOrahvwsSQzXefPX+z+YXz/oU/jD8m0vddL5TM2zpaXaXyvx1j2NmCfrrR9gJv6yLsPsK270QZQVRuBjQBr166dvi/xAozys6/r5inr5UMo85yCc+bJNKxyR2sgsZSk7/30G2e1bt18dR+KcZU7zrLHeczdpvG3dBD3uM3Wde7wEmlwjDNp51wLrEhyaEfa0UCvnuot7bb58knTaFCxtAU4qmvG1qNm2Y+kDoNouPXV5Z3k9CSbk2zeunXrAIqVpsqChpYYa1Kjqm4GLgNekmSvJA8DHg+8tUf2i4A/THJQkp8Cng1cMLLKShNsgLG0CbgD+L0kuyc5s03/yDDrLy0Hg2i49dXlXVUbq2ptVa1dtWrVAIqVpsqChpYYa9IOzgDuDnwbuBR4ZlVtSfKIdtjWjDcB7wU+B/w78L42TVJj0bFUVbfTLBVwCvB94LeAJ7TpkuYwiMlJftx1XlX/1aY5vEQaLONM2klV9V2aE8Xu9CtperNnnhfwR+1DUpdBxVJVfRp40JCqKS1bi+5xW2DXuaSdYJxJkiRNt0EtwN2z63xA+5bUMM4kSZKm1EDWcZut61zS4BhnkiRJ02tQPW6SJEmSpCHJOBaRS7IVuH4nX34AcMMAq7MUyp62csdZ9mLKPbiqJmoaxyUaa9NW7jjLXqrlTlSsLdE4G2fZlrs0yl5OcQbT9/kv1e/dtJW7oDgbS8NtMZJsrqq101T2tJU7zrLHecyTZto+A7/vy7/cSeT3znKXa9mTZto+/2n83k1DuQ6VlCRJkqQJZ8NNkiRJkibcUmy4bZzCsqet3HGWPc5jnjTT9hn4fV/+5U4iv3eWu1zLnjTT9vlP4/du2Ze75O5xkyRJkqRpsxR73CRJkiRpqthwkyRJkqQJZ8NtCUuScddBWu6MM2k0jDVp+Iyzpc2G2yIl2bXr+UgCIsl9gMuSrBhFeV1l/1qSlaMud9IkeeC46zAtjLPpZZyNzrjirC1rLLFmnN3JWBsdf9Om12LjzIbbIiRJVd2RZLckL0+ye41utpcvAn9XVT8cUXkAJPkAcEpVbRthmX+d5PhRldePJO8FXjruekwD42yk5U5UrBlnozPmOIMxxJpxdidjbXT8TRtZmcsyzpbsrJJti/VQ4Gbg61X1b216RhEASXatqjvavy8Edq2qpw673La8FTNBl2QX4M+BFw47IJK8HfiZqnpo+3zo73WSy4GfqKqHDLOchUhyGbBfVf3SuOsybMbZdMRZW85Exdo0xRmMN9bGGWdtmSOPNePsTtMUa/6mTcdv2nKOs5F3lQ5CkjcDhwA/AdwN+FGSjVX1Z1VVo/hStFdLdgHeCHyLtgU9ouDv3H+Ak4H7JHlSVd08jAKTXAycAKwbxv5nKfMi4LCqum/7fLeq+sGoyp+lTpcATwAObp/vXlW3jbNOw2KcTUecteVOVKxNU5zB+GNtzHEGI44142yHOk1NrI07zmDssTYVv2nLPc6W3FDJJO8BjgBOA44CHkWz8N1Lk7wAYIRdzvdoy38WzX8EMIL3tA38JHlRe+XmJ4E1wOVJ9hp0ee2Vi/sC7wJ+O8lxbT0qGc647CTvpAn2byU5rE0eadd+t/Z9eABwPvDaJIdW1W3tf8LLinE2HXHWljtRsTZNcQYTFWtjiTMYbawZZ3eapliboDgDf9M8d1xMnFXVknkA7wQ+3SP97sAzga8Bxw+x/AcCvwE8FnhQm3Yv4LPAhzvy7TqC9+IngRuB09vnuwHXAB8G9hpgOX8F/Fv79wqaALwUWNeRJwM+tvcDVwH3AV4JXAk8eBhlLaBOf9nxPtwHuAD4ALCmTdtlHPUa0rEaZ3eWsWzjrN3nRMXaNMVZezxji7VJirO2nKHHmnG2Q52mJtb8TduhLsv2N21a4mzkB7GIgz+0/RDeMfMBACs6tv8UcAXwh0Mq/83AR4D/AP4b+C/gue22VcDnhxmAvT5c4FTgb4GfbZ/vBnwO+CSw5wDKXAUc1JV2QPsZDCUAgWOBN3U8P5pmHPZHxxWAwM8AT+78bIEHA29Zbj90xtl0xFm7r4mKtWmKs/Y4xhZr446z2T7HYcaacbZDnaYm1sYZZ+3+/U0rzx1nPttBxNlS6g7/AnA2TZfnu5LsV1U/TDulaVV9HbiF5sMaqFm62N8EvKTtct4KPAI4MMkn2/rcMcg6VNWP2i7un+tI/hywL82Xg2rG8D6IZhzzAYspL8lbgEuA/drnu6a5qfYG4Bk0V1CekeTYQY3NTvI04FbgzPZ5qupqmuEM/wa8KsmDq4Y7pKWrTm+k+dw/0z7fraruqKpPAm8AvgG8Mcma9jNaSjHVi3G2zOOsLWeiYm0K4wzGFGuTEGftPkcWa8bZDnWatljzN22Z/6ZNXZyNsvW5ky3WlcABHa3VXwLeTtP1vV+bfrf237fSTDc6yPL76WL/1TbtJ2m6aQ9eZJk9rwoArwb+Afh/Hcf8Apru9j0GeMyX01yhOhS4Z9e2Xdp/70lz5eTvgYcP6H3+JHDczLF0vg/A/WiunnwEeOiIvnvvAT4O/Dywd/d70P79YJohAVcBq0dRryEdq3F2Z/qyjbOO93piYm2a4qw9lrHF2jjirPv71ZU+9Fgzznao09TE2jjjrOPz9zet6zs26Fibxjgb+gEs8uDfCHwQ+BRw8syB03SJvp1mzOz+bfqzaLqc1wyw/H672J/dkbaobu6uL1znuOifpbmh9bj2S3ol8BKamXouYEDjs4HfAa7sStsLOBDYrSt9VVv2QYss82LgYzSzHM36/tHc5Ppm4H3AHrP9JzWg9+GFwCe60nZv/71bV/pDgdcuNPgm5WGcTUectfuaqFibpjhrj2FssTaOOGv3MbZYM852KGtqYm2ccdbu09+08tyxI22gcTaUig/o4C8H/olm+syXA7cBvzjzBW2/dG+nGSv6fOA7wAMGXIdd2nIuBWbWX+gOwL8HLhxQeZ2B1zku+gvt4/c7tv8O8Dqam0x/BGwcUB2eA7xt5ksG/ALwifY/tncDv9xu23XmPVpkeT8FfIj2ylibdhBNF/PJwM915T8cOHAE37//y5037+5OM4ziQ+3n/Taa9UE68w/sqtUoH8bZdMRZu4+Ji7VpibO27mONtVHH2cxxdfw98lgzznYoZypibdxxNvM9GmWsjTvO2v167lijibOhHsAiDvxNwFc7nv8k8M/AozrSdgEeTnPV4lbggQMsf+Rd7F2B9x7abnOa8cD3boPiFuDFXa97dBuoRw6ifODE9j++V9GsL/ItmisCTwL+keYK0sBungUOA/6H5orIvsCvA98D/p1mHPZm4N5j+A6+iuZG4uOBF9P8534+8LI2AC9q/3NasjdvG2fTE2dtmRMXa9MQZ+1xji3WxhFn7b7GFmvGWc96LftYG2ectfv2N81zx6HH2UgPqM+DvidNN/aFwC+1aQcBN7UfxMuBpwL7tNuOYoDd+Yy/i72fcdGP7dq2YoDl/2T7Zfso8HrgpI5tx9GMib7noMpr9/va9j+WK4GrgT9q0x8C/CvwlDF8D/eluTryNZqrZus7tj0fuGzUdRrw8RlnUxZn7b4nKtaWe5y1xzG2WBt3nLX7HVusGWc71GlZx9o446zdn79pnjuOJM5WMGGq6jtJXkRzleDkJIcA59B8QF+imanmFcBzknyqqn57UGW3i+TtC7yGdsrOJF+pqo8m+ShNN/sZwJ8luRZ4LnBcVV03oPIPpRkP/N8zs+0kWVFVP6yqW9r6nUhzpWFm5pyqqoEsLtju7xtJ/pjmPd69qm7uyPJA4Os0V6kG6fdprhL9D/CtqvpXgKr6RJKZ7vyRqqrvA7+Z5EDg5qra1rF5N+AHSfYAbqs2IpcS42wq4wwmLNaWe5zB+GJt3HHW1mFssWac7Wi5x5q/aVP5mzadcTaqVuh8D5rxyLtwZzfyEcB5wDeBS7ry/gRNF/TPDbD8kXex03WDJKMfF31253vILFdfaG5sfRZNl+/RA/7cZ+0uprlKdD0DvgI83+dBx420dNxM2n4+/wf4LnD/UdRpCMdonE1hnM0c9xzbRhZr0xBn7XGMLdbGEWedn21XGSOJNeNs9s9jOcfaOOOs3ae/aZ47jjTOhn5AfR70s2laxm+n6c5e3aYf1AbgeQx3VfuRd7Gz47jkcYyLPrN9z/+l+8vdBv5vt3+v+//bO79QzaoyDj/rTJPjqIQKhho2ExlaNzM2CDrhjDoqBBajXgiijhhZ5sXgRZKNhggduymZIIfMMsEwahgKapIjkV7IQBd2URSVOIUaNUKBkSXp28X7Hs8+35xzOjOz159v798DG873h/2utfb37HXWXv/wMbu/Bjb1XO5rOn9/jNj4EbgSHw/8txO9wa0iDTfhT222LJGmB4Ar4++tUfZ/oOeJzKUOeTZOz+bz2/m7uGtj8izyUc21Gp7Feaq5Js8WpWE0rtX0LOKoTlv8mf53tPyeZcvUMRbAFfFD+Bm+cd/rwBeBzfhmgI/hG+ldnzENF+ATNb8F3A78GZgFPg3cE69/BTzWQ6yueFXGReNLxD4bP/A/srCD+wHgNx3ZTw8ZTnTZ1mXHUsdN5if4cIJTgZvxiawfyfy7O4BPYn0OeAu4tvPZd4HD8+kGzsKX192YM02Z8yvPBu5ZnKsp18bmWeSjqmslPYt4VV2TZ+/EHZVrtT2LNKhOs3yuybMl4ufM3DEWxN4Q7HR8/O/38G7Fz4cM+4GvAKf2GLN4F/uEeFWWrcWfzJwUZXttnP8VYEvcDJbcc+IE4s3iTydOWuKzL8dNbe1E+k7J/Hv7MTDXef0N4AV8j4/LgIPzaWJiWMI0H/JsuJ7FuZpybayeRX6KulbDs8nrVsM1efZOjFG6VtqziKk6Tf87zr8u7lm2zB1HYXwOfyLw3ni9K36Mv8d3IP8tcE6P8Wp3sdcYF30Ji7tzb8NX+lmHPyV4EfhAfLb2RGJNxL06ru3OSQHxJ2bdcdjZKxTgE3Htd06UxZMs7DFy1KaZQzjk2XA9i/M149qYPYs8FXOttmcRq6hr8mxRzNG6VtKzOL/qNP3vWNWzrJk8jkJ5AR8TuxV4Ge9eXIdv5Nfbsq2di16li50646Ifj/x+E99jYybef4jYnBHfA+RFepy424l/WdxIr5sUMD7fApxc6He2AdiDLyG7I35jh4G7SsSvfciz4XoW52/CtbF7FmVQxLWankX8oq7Js6Nijdq1Up5FLNVppv8da3o23zKsSkppxszeTintwlereQ+wx8z2BITLtwAABaZJREFUZY67Fy/4e/Dd5DfjmxI+hHctfxwfw/ugLV7S81jjzJjZ2xPvXYAvWwu+nOmX8JvAS3iX//wKNL0sW5tSuhTfvf0JfIxywsdHn413a++O783hy8peZD0tFdtJwzbgUbzr/qCZvRHv34tvnni1mR3pM+YKaTkPuAV/knMhcK+ZPVwidi3k2Tg8i/M34doYPYM6rpXyLGJVdU2eLZmW0bmmOm0cdZo8m6BkK3EVrdlz8dbrU4Xile5ir71s7VZ8QuVV+Njv3cBr+NOUqzrfe1/GMt8W5Xt9vL47yvzizNd6AxNPZeK9+/AnRfPjw2cY0Pj/ZcpCni18d5CexfmLuybPjiqPYq6V9ixi1NzyYLSeRRy5tpBv1WkL3x1knSbPOvFLBDnGQrotLs75heKVGspSfQ5CxLs8BNwZrzcBNwCnsMK+GD2nYRs+Pvr7+K73H80cbzbK/tvAJyc+24h3fR8CriuR/xYOeZY9v9U9i7jFXJNny5ZLMddKeRaxqrs2Rs8inlw7ukxUp+XNb3XX5FnELhlslQX1fmCO2JsiY5z5cbq78Fbza8BnMsarvmxtJy3b8T0lbqh4nbcDfyXz/jH4KkP3R9l/Dfgn8BTw2c53zsKfnvyOzrKuQz7k2Tg866Qjq2vybMWyye5aac8iVhOujcmziCPXli4X1Wn5y7i6a/KskTluk6SUTrYYw1og1rl41/YhM7sxc6xicxBWkZZtwCP4ZNP9ZvZmznjLpGGdmf27QJwz8E0gZ/Gb3x58bPJa4GF8H5D/EksHm9nh3GlqAXk2Ds8iHdldk2fLU8q1kp5FvCZcG5NnEUeuLYHqtHHUaaP3rFaruaWDQl3sVJiD8H/SswP4JXBa7WuQIW/30RnjjU9q/SmwPl7fij9JOYQ/wbqRgS2T3Nohz+SZPCtyTUoOzWzGtSF7FvmTa40dqtOG59o0eFa9kFo4KNTFHrGKzUFYZXrW1y7/DHm6qyPWhnhvO77S0mb8icnLcdM9D7gT+HDtdA/9kGf1r0HPeZJnDR4lPYt4zbg2RM8iX3KtwUN1Wv1r0HOepsKzJodK1iB3F3utZWvHSEppI773yIX4E5EdZvZSSunrwE141/b9KvvyyLPhIM/apcSQMblWDrnWLqrThsO0eDZTM3hL5K7kbGEvjjl80uMztS/+EEkprQFexcce3w48CzwfQj6JDy/4qpntSym9q15Kx4k8GwbyrG1yexYx5FoB5FrbqE4bBtPkmRpuhTGzV4AHgItSSufXTs9QSCldklJaY2Zvmdl/gCP40rl34ptHPgP8BR+XfQ2AZdiQVbSBPMuDPBOTyLU8yDXRRZ7lYRo9U8OtDj8H/gT8vXZChkBK6XF8gu4jKaUHY2jBd/AJpXeY2a345N0fAT8A/pFSOjOllKolWpRAnvWIPBMrINd6RK6JZZBnPTKtnmmOWyVKLls7dFJKl+JPRp7AN0VMwNPA2cC7zWx3fO9p4EPAxWZ2pFJyRUHkWX/IM7EScq0/5JpYDnnWH9PqmRpuYhCklLYC+4C78dV+TsP33DgDuMbM5uJ755jZq9USKsQUI8+EKINcEyI/0+iZGm5iMKSULgf24qv+HEgpbQI+CBwE3uhM8hVCHCfyTIgyyDUh8jNtnqnhJgZFSmk78CjwBTP7YeXkCDFI5JkQZZBrQuRnmjzT0rFiUJjZL1JKn8Inm64F9pvZm7XTJcSQkGdClEGuCZGfafJMPW5ikKSUdgCzwBVm9nrt9AgxROSZEGWQa0LkZxo8U8NNDJaU0noz+1ftdAgxZOSZEGWQa0Lkp3XP1HATQgghhBBCiMbRBtxCCCGEEEII0ThquAkhhBBCCCFE46jhJoQQQgghhBCNo4abEEIIIYQQQjSOGm5CCCGEEEII0ThquAkhhBBCCCFE46jhJoQQQgghhBCN8z8Sw0mYZ6AEVAAAAABJRU5ErkJggg==\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([np.sum([L2_err[dim][key] for dim in range(3)], axis = 0) \\\n",
    "               for key in L2_err[0].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([np.sum([L1_err[dim][key] for dim in range(3)], axis = 0) \\\n",
    "               for key in L1_err[0].keys()], positions = np.arange(6), \\\n",
    "              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([np.sum([Added[dim][key] for dim in range(3)], axis = 0) \\\n",
    "               for key in Added[0].keys()], positions = np.arange(6), \\\n",
    "              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([np.sum([Missed[dim][key] for dim in range(3)], axis = 0) \\\n",
    "               for key in Missed[0].keys()], positions = np.arange(6), \\\n",
    "              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
}