point_estimation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import os, random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Estimator():\n",
    "    def __init__(self, observations, lam):\n",
    "        self.w_star = np.mean(observations)  # w_star is the optimal global model\n",
    "        self.w_local = observations  \n",
    "        self.K, self.lam = len(observations), lam\n",
    "        \n",
    "        \n",
    "    def solve_mtl(self, lam):\n",
    "        # empirical solutions with the l2 regularizer\n",
    "        w_local_mtl = np.zeros(self.K)\n",
    "        for i in range(self.K):\n",
    "            w_local_mtl[i] = (lam / (1+lam)) * self.w_star + (1/(1+lam)) * self.w_local[i]\n",
    "        return w_local_mtl\n",
    "    \n",
    "    \n",
    "    def get_local_test_error(self, w, w_true):\n",
    "        return (w-w_true)**2\n",
    "    \n",
    "    \n",
    "    def get_overall_loss(self, w_true):\n",
    "        # compute (w_k - w_\\lambda_empirical)^2 given the observations \n",
    "        w_local_mtl = self.solve_mtl(self.lam)\n",
    "        test_loss = []\n",
    "        for i in range(self.K):\n",
    "            test_loss.append(self.get_local_test_error(w_local_mtl[i], w_true[i]))\n",
    "        return test_loss\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Simulator():\n",
    "    def __init__(self, num_mali=0, tau_a = 0.5, theta=0, tau=0.1, sigma=1, K=20, n=20, num_trial=1000):\n",
    "        \n",
    "        self.optimal_lambda = sigma**2 /n *(K /(K * tau**2 + num_mali * (tau_a**2-tau**2)/(K-1)))\n",
    "        #print('optimal lambda: ', self.optimal_lambda)\n",
    "\n",
    "        sigma_n_square = (sigma**2) / n\n",
    "        self.sigma_w_square = 1/(1/sigma_n_square + (K-1)/(K * (tau**2) + sigma_n_square + (num_mali * 1.0 /(K-1)) * (tau_a**2-tau**2)))\n",
    "        #print(\"sigma_w_square (min obj value):\", self.sigma_w_square)\n",
    "        #print(\"min variance:\", 2 * (self.sigma_w_square**2))\n",
    "        \n",
    "        self.w_underlying = np.zeros((num_trial, K))\n",
    "        for i in range(num_mali):\n",
    "            self.w_underlying[:,i] = np.random.normal(theta, tau_a, num_trial)\n",
    "        for i in range(num_mali, K):\n",
    "            self.w_underlying[:,i] = np.random.normal(theta, tau, num_trial)\n",
    "            \n",
    "            \n",
    "        self.w_empirical = np.zeros((num_trial, K))\n",
    "        for j in range(num_trial):\n",
    "            for i in range(K):\n",
    "                self.w_empirical[j][i] = np.random.normal(self.w_underlying[j][i], sigma / np.sqrt(n))\n",
    "                \n",
    "        self.num_mali = num_mali\n",
    "    \n",
    "    \n",
    "    def run(self, lambdas, num_trial=1000):\n",
    "\n",
    "        w_underlying, w_empirical = self.w_underlying, self.w_empirical\n",
    "        \n",
    "        sols = {}\n",
    "        variance={}\n",
    "        for lam in lambdas:\n",
    "            loss = 0 \n",
    "            var_ = 0\n",
    "            for trials in range(num_trial):\n",
    "                e = Estimator(w_empirical[trials], lam)\n",
    "                l_vec = e.get_overall_loss(w_underlying[trials])  \n",
    "                loss += np.mean(l_vec[self.num_mali:])\n",
    "                var_ += np.var(l_vec[self.num_mali:])\n",
    "            \n",
    "            sols[lam] = loss / num_trial\n",
    "            variance[lam] = var_ / num_trial\n",
    "            #print(lam, loss/num_trial, var_/num_trial)\n",
    "        return sols, variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  0.    0.2   0.4   0.6   0.8   1.    1.2   1.4   1.6   1.8   2.    2.2\n",
      "   2.4   2.6   2.8   3.    3.5   4.    4.5   5.    5.5   6.    6.5   7.\n",
      "   7.5   8.    8.5   9.    9.5  10.   20.   30.   40.   50.   60.   70.\n",
      "  80.   90.  100.  110.  120.  130.  140.  150.  160.  170.  180.  190. ]\n"
     ]
    }
   ],
   "source": [
    "lambdas = np.concatenate((np.arange(0, 3, 0.2), np.arange(3, 10, 0.5), np.arange(10, 200, 10)))\n",
    "print(lambdas)\n",
    "\n",
    "K=8\n",
    "plotting = Simulator(sigma=0.8, n=30,  tau=0.1, tau_a=0.1, K=K) \n",
    "mse, variance = plotting.run(lambdas = lambdas)\n",
    "\n",
    "K=16\n",
    "plotting2 = Simulator(sigma=0.8, n=30, tau=0.1, tau_a=0.5, K=K, num_mali=int(0.5*K)) \n",
    "mse2, variance2 = plotting2.run(lambdas = lambdas)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, sys\n",
    "import re\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "import numpy as np\n",
    "import matplotlib.ticker as mticker\n",
    "from matplotlib import rcParams\n",
    "matplotlib.rcParams['font.family'] = \"sans-serif\"\n",
    "matplotlib.rcParams['font.sans-serif'] = \"Arial\"\n",
    "from mpl_toolkits.axisartist.axislines import Subplot\n",
    "import matplotlib.patches as mpatches\n",
    "from matplotlib.lines import Line2D\n",
    "from matplotlib import rc\n",
    "rc('text', usetex=True)\n",
    "\n",
    "\n",
    "plt.rc('font', size=14)          # controls default text sizes\n",
    "plt.rc('axes', titlesize=16)     # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=13)     # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=12)    # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=12)    # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=14)    # legend fontsize\n",
    "plt.rc('figure', titlesize=16)   # fontsize of the figure title\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def errorfill(x, y, yerr, color=None, alpha_fill=0.15, ax=None, label=None):\n",
    "    ax = ax if ax is not None else plt.gca()\n",
    "    if color is None:\n",
    "        color = next(ax._get_lines.prop_cycler)['color']\n",
    "    if np.isscalar(yerr) or len(yerr) == len(y):\n",
    "        ymin = y - yerr\n",
    "        ymax = y + yerr\n",
    "    elif len(yerr) == 2:\n",
    "        ymin, ymax = yerr\n",
    "    ax.plot(x, y, color=color, label=label)\n",
    "    ax.fill_between(x, ymax, ymin, color=color, alpha=alpha_fill, linewidth=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/ipykernel_launcher.py:20: UserWarning: Attempted to set non-positive left xlim on a log-scaled axis.\n",
      "Invalid limit will be ignored.\n",
      "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/ipykernel_launcher.py:35: UserWarning: Attempted to set non-positive left xlim on a log-scaled axis.\n",
      "Invalid limit will be ignored.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgYAAADICAYAAACTQilNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3da5RcZZkv8P+769pV3V1V3ekkdAhJuhNBQ7gkAQYQBBNEBxS5K+pRRw06nkFh1pLhHD+cWefMcWBmIc44DkS8IYKSgI7jcEaTjCOiEkiCgsAg6Y5A0p100ve678tzPuxdlUqnL1XdVbXr8v+t1au6dr1V/fbu3s9+9nvbSkRAREREBACa2xUgIiKi2sHEgIiIiPKYGBAREVEeEwMiIiLKY2JAREREeUwMiIiIKI+JARVFKbVXKbXF7XoQUWUopb6glNrhdj3mopTqU0rdXWRZxq158LpdASIiohLcCqDf7Uo0MiYGRERUN0Rkp9t1aHTsSiAiIqI8JgZ0AqXUA0qpUedrh1KqZ4ZyUaXUNqdcX2E/nlJqi7Mt99rmgtf25voyndf3zvQziGhh5jgWe5zjb9QZW3BewWsPKKW2Tfms/BiEOY7/vUqpGwpiSdT5vq/gmF8/V/1m+ay9qsh4Q/PDxIDylFJ7nW9XiUgMwGwDfHYB+IFTbgOAWwtO8CMANhR8xrYp770LwI3O65jj5xDR/M12LO4AsEdEYiJyBYBowWsPALhhymfd6mwHZj/+AeDrzuMqAJsBbBSRXqf8jTg+RmCuWHHCZ4nIWIm/I80DEwMCADhZdo+I3Jo7+ERkp4icNMjHKbsewN1KqT4Ae2EHlc3O+7YXfMZWANHcFYLjsYID/Ac4MSARUZnMdCwWHu8FxXcUvG8fgH6l1A0A4By/HSKyfa7j37GnIJaMAehxrvyjItJfUKe5YsXUzyr6d5zP/iIbBx9STg+KH+nbA2Cnc5VxEqVUFHarwExNen1TnncU+XOJqASzHIs9APbN8fYHYLcSbHceHyt474zHv6MwydiplPoS7Kv5HqXUTtgthmNFxIoTPms6RX4GlYAtBpTTD/uAL7bsxulecA7SAwCeE5ENIrKhTPUjohLMcSwWc7xvBbDZ+ZybcLwbYcbjv8AJV/cico+I9AKIwb4Q2FJCrJi2pQBgvKkUJgYEID8FqF8plTv44TT9nZSFF5TdNqXsetgHfRTO1UiuKZKIqm7GYzE35S93vDsn2MJuBTjN8zthX+mPON0Lcx3/J1FKbc7FEeczR+aqXzl+R5o/JgZUaBMA5GYlALgZM3cvTFvWGZOwHUCf0/943gzvJ6IKKuJY3ABgo3P8fh12EjDVAwC24HhrQU4psQIA7iwoOwZgazliBeNNZSgRcbsOREREVCPYYkBERER5VZuV4PT9jMHuD+rP9VcVU8bpnxqBPVimw5mSQkQNgLGBqLZUpSvBWfjiztycWaXUNhG5sZgyzqCYbbmpMUopERFV8UoTUcUxNhDVnmp1JWzGlCkn04xgnbaMiIwVHPjrYU+hIaLGwNhAVGOqlRhEAQwXPM81/RVdxjnwb56yUhcR1TfGBqIaUzcrH4rIPqVUv1Kqz1ko4wTOTTW2AEA4HN5wxhlnVL2ORG7Yu3fvMRHpcrsebmFsIDrZQuJCtRKD3KChnA6cPOd12jIFa3TvdJbQhFJq89R7cjuDjrYCwMaNG2XPnj1l/yWIapFS6nW367AAjA1EFbCQuFCtroTHABRm8tFpRh7PVGYjTr7JTrFr+hNRbWNsIKoxVVvgaOrSugVLcu4FsMnJ+GcqswXH+xX7RWT7bD+LVwXUTJRSe0VkrrXraxZjA1H5LSQuVG2MwdTmvYLtG4oow9HGRA2KsYGotnDlQyIiIspjYkBERER5TAyIiIgoj4kBERER5TExICIiojwmBkRERJTHxICIiIjymBgQERFRHhMDIiIiymNiQERERHlMDIiIiCiPiQERERHlMTEgIiKiPCYGRERElFe12y4TERFR5YhpAqYJsawFfQ4TAyIiohomlgUxDMAwIIZxPAEwTXtb7vsyYWJARETkMtF1+yt38tcNwHCeL7AFoFRMDIiIiKpEstl8EmBls0AuIRBxu2p5TAyIiIgqQHQdViYLyWYgmYydFFT56n8+mBgQEREtkIhA0mlY6bSdAGQyZe33ryYmBkRERPMg2SysVMpOBlKpmuoOWAgmBkREREUQy4KkUseTAV13u0oVwcSAiIhoBlY2C0km7UQgnW6YVoHZMDEgIiIqICKwEglY4+P2zIEmw8SAiIgIgBgGzMlJWJOTdTtwsByYGBARUVOz0mmYExOQZLIpugrmwsSAiIiaTr67YGICVibjdnVqChMDIiJqGuwumFtDJgZZM4s/jv9xznKxYAyRQKTyFSKi2mCaMMfHZy2ifD5ooVCVKkTVYmUydutAIsHugjk0ZGJQrNH0KAAwOSBqEmIYMEZG5iznXbQInra2KtSIKsnuLkjCmhhnd0EJmjoxAOzkQCmFdn+721UhohphHDsGiMDTzrhQj8Q0YU5MwIrH7dsVU0maPjEAgJGUfQXB5ICIcozhYTs5iLBFsV6IrsMcG2N3wQKVlBgopdpFZKJSlXHTSGoECgptfjYfEpWqUWODMTJiJwfRqNtVoVmIM3bEmphgQlAGWonlDyilGvbMOZwaxmR20u1qENWjho0NxugojNFRt6tB0xARmBMT0A8dgjk+zqSgTEpNDLYDuLsSFakVw6lhxLNxt6tBVG8aOjaYY2NFDVqk6rHSaRgDAzCGhzntsMxKHWPwMwBfV0r1AtgBYCz3gog8WM6KuelY6hiUUgj7wm5XhaheNHxsMMfHARF4OzvdrkpTE9OEOTICM84LuEopNTG4C0A/gE4AHyjYLgAa4uDPOZo8CoTA5ICoOE0RG8yJCTs5WLTI7ao0JXNyEubICMSy3K5KQyspMRCRjZWqSC06mjwKFVII+bjYCdFsmik2mJOT9oDERYuglHK7Ok1BslkYw8Ow0mm3q9IUSh1jAMAegayUemejDjYqdDR1FEk96XY1iOpCs8QGMx6HcfQoB7tVmIjAHBuDPjDApKCKSl7HQCl1P4AtzlNRSm0TkQ/M9h7nfTfA7neMAugXkX3FlnG2A8B5AJ4Tke2l1nu+RARHU0fRhS62HBDNotlig5VIwBDAu7iLLQcVINksjGPHuGKhC0pqMVBKfQlAB4CYiGiw+xM7lFL/d4739QC4QkR2OgfuXcWWUUqtBzAmIttF5E7YA5yqOqk4lxyw5YBoes0aG6xkAsbQEFsOyswcH7dbCZgUuKLUroQbROQmERkHABEZA3ATgFvneN9mFIxSBvIHdTFlegBcUbB5xNlWVbnkIGWkqv2jiepB08YGK5mEceQIk4MyEF2HPjgIY2SE+9NFpSYGJ83TcQLAXO1oUQDDBc+nO4CnLVNwNQDnaqBjuqbGahARDCWHmBwQnaypY4OVSsE4fJij5RfAjMc5lqBGlJoYPDa1aVAp9c8AdpavSrO6G8CG6V5QSm1RSu1RSu0ZHh6erkhZ5JKDtMF/XqICdREbjlZwkSIrnWZyMA9imtCPDNmDObnvakKpicGdAN6llHpNKfVTpdQw7Ga+T87xvrEpzztgz3kuuowzyOgBEZn6PgCAiGwVkY0isrGzwguQMDkgOkldxIaujo65fo8FsTIZOzngSnxFsVIpu5UgmXC7KlSg1MRAnPnKn4F9JbBFRNYUcfOUxwD0FjyPTtPkN2MZpdRmAPtEZJ9SKuoMRnKVJRaGkkPImBwcQwTGhjwmB3MTERijo9APH+ZtkWtQqdMVDyilVorITpTQRCgiY0qpbc5BDBSsqa6U2gtg00xlnEFGDwAYc6YE9YhIrMR6V4QlFo4kjmBJeAkCnoDb1SFyE2NDASubhT54GL6lS6C8vLt9IU5DrH2l/sduB3AP7KuCkjgBY7rtG2Yr41wZ9E7dXityycHS8FL4PX63q0PkFsaGKUTPQj98hMlBATMehzk8zLEENW6+N1HqQYPeKGU+ClsOmBxQk2JsmIboWeiDg/AuXgwt0LytimJZditBgmMJ6gFvolQmpphMDqiZMTbMQAwDxuAgPIu64GltvpuyWZmMPeNA192uChWp1MTgnUUMJmpaTA6oiTE2zEJEYBwdguhReGM1MQyiKozRUZhjUyeVUKnEsmBNTsIcH4c1MQlzcsJ+PjEJKx6HlYjDiidgJhKwEglIcmGr9M538OHkgn5qA2NyQE2KsaEI5tgYRNfh7Wrs+yuIrttdB1ysaFZiGDCGh2EMDcE8dgzGyCjMkRGYo6MwnEdzbAxWPA7MNC5D06C1tkILh6G1hqGFwvCccsqC6lW1wYfNhMkBNSHGhiJZiQQMw4B38eKGHJTIAYbHiQjM0VHoBw9BHzgE/dAAjKEjMIaOwjh6FObICDBl6Wfl98PT0QFPLAbfqacieOZaeCKR/JcWicDT1gatrQ2e1jaoUMv0Sea3vzXvenPwYYXkkoOl4aXweXxuV4eo0hgbSmBlMvagxCVLoPkb4+JBTBPGseGmXKxITBP6wYPI9PVBf/NN6AMDdjIwOAhJFSyh7/XCt2QxvF2L0bJ+PbxdXfAu7rIfFy2Cp6MDWjjsemsSBx9WkCkmDicOMzmgZsDYUCIxDBgDA/B2dUEL1/egRCuZhDE83BSLFYmuI/vGG8js70O2bz8yr+1H9sABSG5dBk2Dd8kS+JYtQ3DdmfB1L4Nv2TL4lnXbXUgej7u/QBFKSgyclc2oBEwOqBkwNsyPiEAfGoI3FoMnWtU7RpeFmKbdJx6Pu12VihBdR/bAAWT277cTgf37kTlwAHASINXSgkBvL9re8x4EVvcisHo1fN3dUL76jvUld3AppVYC2JxrHlRKnQOgnyOSZ8bkgJoBY8P8GaOjEF2HZ9Ei15uRi2XGEzBHhhtu6WdzYgLJPXuQ3L0byT178yP8tdZWBFavRuT91yDQuxr+1b12EqCVemeB2ldSYqCU+iSATwOI4HjzYC/sZsSby1u1xmKKicPJw1gaYnJApbPEgojAgvMoFiypncFdjA0LZ8bj9oyFRYuganjcgeg6jJERWAucEldLsgcPIvnMbiR3P4P0y68AlgVPLIbWSy9By/r1CKxZA++SJXWTtC1UqS0GdwLYCGBvboOIPK6U2lrWWjUo02Jy0MwssWCKaT9aZv65iNiPEJiW8+iUyyUBdYCxoQysTAb6wAC09nZ4otGauhoVEXvq3Pg4ZMpI+nojpon0Sy/ZrQK7d0M/NAAA8Pf0IHrzTQhdcAECa9bU1P6vplITgw4RGZ8ma2qONKoMmBw0HkssGJYBwzJgignTMmGIccLJP5cANDDGhjIREXshm3gcnlgMnrY2t6tkT0EcHa3rwYVmPI7Unr1OF8Eee20ArxctZ5+N9ve/H+Hzz4d38WK3q1kTSk0MdimlroM90hgAoJT6AUq4mxo5yUHiMLpCXQh6g25Xh+YgIjAsA7ql578Kk4A6uaKvNMaGMrOn/x2DNTEBT2cntGD1Y4WVSNitBNls1X92OegDg0ju3o3E7meQ/v1LgGlCa29H6MI/QfiCC9By7rnQQiG3q1lzSk0MPgVgF4BepdRPYTcd9gPYVO6KLYQlFtLG3CtuBTwB1/qMTDFxJHkEncFOtPpbXakDHSciMMSAbur5x1wCYFj1e5VURXURG+qRlc3CGhyEFg7D29FR8UWRRARWIgFrfLzuEgKxLGRefgWJXBfBm28CAHwrViB6/XV2F8Hpp9fFlEE3lTpdcRzARqXUJgA9AO4RkV0VqdkC9I334caf3DhnuVXtq3D7htuxKrKqCrU6mYjgWOoYDMtANFh/U5XqmWEZyJgZZMwMsmYWGTPT6E39FVUvscGMxxH/5dOzltFCIbSce07N9S9biQT0ZDK//K0KBst6YSPZrL3W/uRk3c00MIaGMPGznyG+YweMoaOAx4OWdevQ/qfvQeiCC+Bb4BLBzUY1YjBc/rblcvtDt89aRrd0/KT/J4jrcXzkrR/BNauvgUe5l0WGfWEsaqmfqUr1xBLrhAQgY2ZgWvUV+GazKrpqL9cRKM6ZwRbZtnLlnOUCb3sbuj7/OfiXL698peZJeTzQQiFooRBUywzL4s7BymYhyaTdQlBvrQOmieSePZh88kkkn9sDAGg55xy0XbEZofPPr/tFoxYq2NMz77jQkInBWeeeJT/+zx/PWW48M46v/vareGbwGZzZeSZu33A7FofcG3wS9AbR1dIFj8ZmroUoTACyZhZZs74CXqmYGBRv/emnyy/uv3/WMulX/4CRBx+ElU4jdsstiN54Q83f00Bpmp0c+HxQHo/dVO712q0eIvZ9C0x7AKzoOiSTsb/q8H4G5uQkJv/93zHxb0/COHIEno4OtF35LrRdeSV8S5a4Xb2awcRgimITA8Buzt/1xi5sfdGeVfXpsz6Ny5df7tqVu1fzYkloCWcsFMm0zHwSkEsEmm0wIBOD4m1Yt05+/eO5Y4MxOorh++9H4qlfwr9qFbo+/zkE3vKWKtSQZpI9eBDjP/wR4rt2QTIZBNetQ/t734vwhX9S84mbGxaSGDT93lRKYfOKzThz0Zm4b999+PK+L2P34d3487P/HJFApOr1MSwDg4lBdIW60OJtqfrPr3WGZSBtpJE200gbaQ4MpIrwxmJYctddSFx2GY7909dw6PY7ELn2WsQ+/CFXZgc0s/RLL2Fs++NIPvMMlM+H1ssvR/v7r0FglTtjw5pBWRIDpVR7vS97ujS8FH/z9r/Bj/b/CA+//DBeGX4Ft62/DRuXVP9CzBILRxJH0B5oRywQa+pxByKCtJlGykghZaSgm7rbVaIS1HtsCF94IYLr1mHkG9/E+OOPI/HrX6Prc7eh5eyz3a5aQxMRJJ99FmOPPYbMy69Aa29H9JZb0H71VfDGYm5Xr+GVNOxWKfXaNNvOBbCtbDVykUd5cP2a63HvZfeiPdCOv/7NX+Nrv/taUVMfK2EiM4GBxAAyZsaVn+8WwzIwmZ3EkcQRvDH5Bo4kjmAiM8GkoIY1cmzwtLai63O34ZS//RIAYPCv7sLRr3ylYW8c5CYxTcSfegqHPvvfceR//TXMY8fQ+ZnP4LTvfBsdH/kwk4IqKbXFoHPqBhF5XinVUP2bqyKrcO877sV3X/ku/mX/v+B3Q7/DHRvuwOkdp1e9LrqpYzA+iEgggmgg2rCtBxkzg6SeRMpINfxgwQbV8LGh5eyzcerX/gmj3/sexp/4IZLPPodFn/1zhC+6yO2q1T0xTSSeegqjj34f+ptvwrd8Obr+8i/Retk7OH7ABUXtcaXUz2CvaBZxFi8p1ANgpNwVc5vf48cnzvwEzltyHr6878v4wi+/gJvfcjNuOv0meLXq/6OOZ8aRNJJY1LIIAU+g6j+/3HKLUCWNJFJ6CqY0zvTBZtJssUELBtH5iU+g9dJLcfS+r+DI//4/CF/ydnR+5jO8mp2HfELwyKPQDx6Eb+UKLL7rrxC++GIuQuSiYs9w22CveX4FgO1TXhtBAy97elbXWfjqO7+KB154AI+++ij2HNmDOzbcgVPbTq16XXKtB63+VkQCEfi0+pq5YFgGUkYKST2JtJnmgkKNoSljQ2DNGiz7yn0Y2/44xh55BKnnf4vOLVvQunlTw7bqlZNYFhJPP43Rh79ntxCsXIHF//N/IHzRRTW3sFQzKmm6olLqMRG5qYL1KYtSpiuW4ulDT+Nrv/0aMlYGf7b2z/Cnq/7UtSCglEKrz04Q3GjBKFbGzCClp5A0kuwiqJBamK5YL7Gh2OmKpci++SaO3vcVZF5+GS3r12PRbX/B+fQzEBEkf/MbjH73YWT/+Ef4li9H7MMfQvjtb2dCUGZVW8dAKXUOgH4RmVBKtQPY4ry0tZZGHlcqMQCA4dQw/uH5f8C+oX1Yv3g9bjv3NnS2nNS9WjW5BKE90F4TLQimZSJlpOyZBOwiqIoaSQzqIjZUIjEA7CvgiSefxMg3vwVYFjo+9lG0v/e9bA53iAhS+/Zh9DsPIfPaa/At60bsQx9C+NJLuY8qpJqJwXMAbhSRPzr9iZ2wmwtHReTm+VSgEiqZGAD2P/mTB57EN1/6JgJaAJ8957O4eNnFFft5xQp6gwj7wgj7wtBUdbLv3FiB3LoCbBWovhpJDOoiNlQqMcgxhoZw9B+/itSePQicfjq6bv88/CtWVOzn1YPUiy9i9KGHkP79S/AuXozYh25B66ZNTAgqrJqJwYiIdCilIrAP+phzhTAsIu5dNk9R6cQg5+DkQdy79168NvYaLl9+OW4961aEfe6vz62UQou3BS3eFgQ9wbKtoigi0C39hJUGOYXQfTWSGNRFbKh0YgDYx0n85/+J4QcegJVMInrTTYjeeEPTLYyUfuW/MPrdh5B6/rfwdHQg+sEPoP3KK6F87rdsNoNqrnyYG2G8GcCBgibCphxtc2rbqbjn0nvw2KuP4Qd/+AFePPYiPr/+8zi7y93FT0QEST2JpJ4EYK/P4Pf44fP44NN88GpeeJQHSilo0PLjJEQEFixYYsG0TJhi5m87rFv2bYg5YJBmwNjgUEqh7Z2XI7R+PY49cD/GHnkEkz/9KWIf+TDaNm9u+CvlzB/+gNHvfQ/JZ5+D1t6Ojk99Eu1XXQUtUP+zqZpFqYnBTqfJsAfAlwDAuc3qnnJXrF54NS9ueest2Lh0I+7dey+++Ksv4upVV+Ojaz+KoLc2rhBMMfMrBxJVCGPDFJ5oBEvuvBPpq67G8De+gWP3fQXjP/oROj/+Z2g5b2PDzV7I7N+P0Ye/h+Tu3dBaWxH76EcRed97oYVCbleNSlTyTZSUUtcDGMvda905+FFL916vVlfCVBkzg4deegg/7v8xusPduH3D7Tij44yq14OaSy10JQD1ERuq0ZUwHRFB4ulfYeRb34IxOIjg2Wej85OfQGD16qrXpdzSr/4BY9//PpLPPAOttRWR665F5H3XQAszIXBTtW+itBd2c2HOMID++fzwRhPwBPCpsz6FC065APftuw93PnUnrn/L9fjgGR+siRkDRBXG2DADpRRaL3k7wn9yASae/H8YfeQRHPqL29B6+eWIffS/1d30RhFB+oUXMLZtO1J799otBB/+MCLvvwZa2P1xVrQwJSUGSqlPAvg0gAiAB53NvQDuAlAzI4/dllsU6cEXH8S2P2zDc4efwx0b7sCqCO8GRo2JsaE4yudD5Jr3oW3zJoxt22bfRviXv0TkmmsQvfkmeNra3K7irMQ0kfjVrzD+xBPIvPoHeGIxdHz84/YYArYQNIxSZyW8BmAjgL0isrpge02NPHarK2E6zw4+i3/87T8ino3jlrfegutWXweP1tiDj6i6aqEroV5ig1tdCTMxjh7FyEPfRXzXLmjhsD1y/+qrofn9blftBOb4OCZ/+jNM/Nu/wRgagre7G9Frr0XrFZs5qLBGVbMroUNExqcZNNNYo2jK6PxTzsdXO76Kf/7dP+Ohlx/Cs4efxe3rb0d3a7fbVSOXmWIiY2TyA0PTZhoZw54GmjbTyJgZZM1s/jFrZpG1stBN3X609FqaLsrYMA/eri4s/ss7ELn2/Rj5xjcx8vUHMfbo9xG+9FK0bd6EwBlnuDZIUSwL6RdewOSOHUj88mmIriN41lno/PStCJ1/fsPPrmhmpSYGu5RS18G+aQoAQCn1AzToeujlEglEcOd5d+IXB3+B+1+4H7f9/DZ8bO3H8J5V74FH8eCqVyKChJ7ApD6JyewkJrITiGfj+ceEnkDCSCCpJ+3v9YR90ygnEZjP7bS9mhd+zW9PP9V8tTR2hbFhAQI9PTjlb/4PUi+8gMl//yniu3Zh8skn4VvWjdZNm9D6zndWZRyCiCDb14fE008j/p+/gHHkCLRwGG1XvgvtV1/d9Is1NYtSuxIiAHYBWA9gB+ymw34Am2pp2dNa6kqYqnBJ5d5ILz511qewtnOt29WiAmkjjZH0CEbToxjJOI/pgsfMKMbSY5jUJ2GJNePntHhb8itRhrwh+9EXQsgbQtAbtB89QQS9QbR4WxDwBBD0BhH0BBHwBBDwBhDwBOD3+PPJwHQrWtZIV0JdxIZa60qYiZVIIvGrpzG56z+QfuEFAEBw3Tq0btqE8EUXlnUsgpVMIvXii0g9/zySzz4HY3AQ0DS0nHMO2jZvQuiii9hdUIeqtvJh/k32NKQe2GujFzUVSSl1A4AxAFHnfftKKaOUWg9go4hsnetn1XJiANhZ+VOHnsK3X/o2jqWO4ZJll+Bjaz+GxaHFbletKWTNLAYTgxiID5zwOJwexkh6ZNr1HryaF7FADB3BDsSCMcQCMbT529Dmb0O7vx2t/la0+9vz28K+cNVag2ohMcip9dhQL4lBIf3IEcT/4+eI79oJ/dAAAMC7dCn8vT0I9PbC32M/ejo7Z+12EBFY8Tj0gQFk+/qR6etDtq8Pmf37AdOE8vsRPGsdwhdfjPCFF8ITiVTrV6QKqPZ0xdy85PxBr5Rqn+2qQCnVA+AKEbnVeb4NwI3FllFKbQZwK4Dn5lPfWqOUwjtOfQcuWHoBnnjtCTz+2uPYfXg3rl99Pa5bc13NLIxUz3RLx5HEEQwkBjAQH8BAYgCD8UEcih/CsdQxyPEWb7T729Hd2o1VkVXYsHgDYsHjCUDusc3X1nAL0lQCY0P5+ZYsQeyDH0D0Azcj81+vIvXC75Dt60O2vx/JX/06X05rbYUWCkEFAlB+v/3o88FKJGCOjsIcHwcM43j5cBj+nh5Er78OLeeei8Db3lZzgx7JHaVOV3xNRNZM2XYugL8FcOUsb90MO9svfN/6KVcGM5YRkZ1OcIiWUt9aF/QGcctbb8HmFZvx7Ze+jUdffRQ73tiBj6/9OC5ZdglPRHMwLRNDqSEcih+yr/rjg/lEYCg1dEIzf6uvFd2t3VjbuRbdrd04JXwKulu70R3uRqu/1cXfojEwNlSeUgrBt56B4FuPL5pmJZPI9h9Apr8f+ptvwEqlIdkMJJOBlclC9Cw8sRj8q1bCE43BE4vCt3gx/L298C5ZwhhTQ5SmAZoGKA1KUyd/79yW+ng5Zf/9lMo/P2HbApTaYnDStCMReV4pNaZwcmAAAA14SURBVFdzRRT2Yic5I7CbG/eVWKYhLQ4txhfO+wKuWnUVtr64FX+35+/wk/6f4ENv/RDOWnRW0x+8IoJjqWM4MH4AByYO4MD4Abw+8ToGE4Mn3Na5xduC7nA31sTW4B2nvsM+8Ttf7f52F3+DpsDY4AItFELwzLUInslxSrVCKQV4PPYJ3OMBNM2ewaFpUF6vfbL32Cf33Em+1mZ4FJUYKKV+Bnu0ccS5pWqhHhy/gQotwNpFa3HvZfdi1+u78PArD+OLv/oieiO9uHbNtbi4+2J4tXn1/NQVSywMxAfQN96HvjH768D4AUzqk/kyS0NLsaJ9BS7svhDd4W6c0noKusPdiAaiTZ9EVRtjAzUTpRTg9UJ5vfbJfOr3Hk/NneTno9gzzTbY85GvALB9ymsjmHtKUm7QUE4HTl4qtZgyM1JKbQGwBQC6l9fvGgEe5cG7Vr4Lly2/DD9/8+f40f4f4e/3/D2+0/IdvK/3fXjXinch5GuMFcYssXAofgh9Y314bew19I31oX+8Pz/4z6f5sKJ9BS7qvgirIqvQE+nBivYVDfP7N4i6ig3Lu+s3NlB1KKXsMRo+X/6xkU76xSgqMRCRrwP2QJ/c9yV6DMDdBc+j04w8LqbMbHXcCmArYM9KmEcda4rf48eVK6/EFSuuwJ7De/DE/ifwjd9/A9//r+9j02mbcOmpl+ItsbfUzRWyJRYGE4PYP7bf/hrdj77xvnwS4Pf40RPpwablm9Ab7UVvtBfL25Y3RStJPau32LBh3bq6jw1UPsrnh/L7oE1NBJrcvKYrzusH2aOH80Rkp7N9L+y5zmOzlMmNPI4CuDu3fSa1Pl1xvl4deRU/3P9D7D68G4ZlYHFoMS5ZdgkuWXYJeiI9NZMkmJaJg/GD6B/vR/9YP/aP7Uf/eD+SRhKA3RLQE+nB6uhqrI6uxprYGpzaeiqXip6nWpquOB/VjA31OF2RFk75fAUnfjsZUD5fzcTMSlBKVXcdg1rXqIlBTjwbx+7Du/HUwafwu6O/gykmloaWYt2idVi7aC3O7DwTi0OLK/5Pb4mFkfQI3ph8A6+Pv47XJ1/H6xP2l27ZS/X6Nb/dDRDtweqInQSwJaC86j0xqCYmBo1Nadrx6ZoFLQBKO3lhsEa3kMSA0bkOtfpbsem0Tdh02iZMZCfw64FfY8/hPfjN4G+w440dAIBFLYuwqn0VTms/DcvbluPU1lPR2dKJSCBS9DK6ppiIZ+MYSY9gJD2CoeQQhpJD+SmBg4nBE5b1jQViOK39NFzVcxV6Ij3oifSwJYCIKkb5fNACAahg0H7kOgxlwcSgzrX72/Hule/Gu1e+G5ZYeGPiDfx++Pd4efhlvDH5Bp4feh6GGCe8p81nr8wX8AYQ9AShlLJXRRMLGdO+iU9STyKux09YCAiwB0cuDS/FKeFTcHbX2ehu7cbytuU4re00RAJcKY2IKkf5fNCCQTsRCAbt6X9UdtyrDURTGlZGVmJlZCWu7rkagN3ffzh5GIcmD+XX/R/NjCKlp/J38RMReyQuFBZ5FuXX8c8t99sR7EBHsANdoS7EgjHe+ImIqkJ5PHYS0NJiJwIcGFgVDZkYeDUvOltmvwW8QJA1s0jqyVlvhFPvPJoHy1qXYVnrMrerQkQ0Jy0QgGppgRYK8eZNLmnIxEBTGtr8xd19TIKCpJFEPBuf9uY5RNQ4lM835+2LxRJIKgkrmYRYjXvRUCvyrQKhMLSWYNOsFVDLGjIxKIVSKn9rXMMykNATiOtx6KbudtWIqNw0DVqoiAWyWsMQEUgqBSvpJAmmOff7qCjK57NbBHI3fWrgaYP1qOkTg0JezYtIIIJIIIKMmUE8G0dST56wHj8RNQelFJRz8gIAK52GlUjYSYJhzPFumkoLBI4nA5w9UNOYGMwg4Akg0BJAR7ADKSOFyewkuxqImpjmjIRHZyesTOZ4kqCzdXE6SqnjYwVaWjiDoI7wLzUHpRRCvhBCvhBMy0RcjyOhJ5A1s25XjYhcogUC9sC4jg5Y2aydJCSSEL2544LSNCcZCEMLtTTlwkKNgIlBCTyaJ9/VkDWzdpKQTbCrgaiJaX4/NL8fiMUg2aw9JiGRgJVtjiRBeb3HuwiCQY4XaABMDObJ7/Gjw9OBWCCGlJFCXLdnNTTiEtNEVBzl98Pj98MTjUJ0/fjAxUymoWJDbryACoXspIgaChODBZra1ZA0kkjoCaSNtNtVIyIXKZ8PnkgEnkgEYlmQdBpWOm3PdKiz1gTl9doLDLWEOKWwCTAxKCOP5kGbvw1t/jbolo5E1p76aFgcwUzUzJSmnTDDQUwznyhYqXTNjU1QHg9UIAitxVl6mK0CTYWJQYX4NB+iwSiiwSinPhLRCZTHAxUOQwuHAQBiGPnWBMlmIbpeta4HpZR9F8JAAMoZVMmlh5sbE4MqmDr1MaEnkDSSDdXnSHNTSkGDBihAQeXvTwEcfz7b97nnuc/KfU/1T3m98LS2Aq2t+W2i6yd+GUb+seTPVwrweu2ExOu1b0XsfMHn44BBOgETgyoqHI9giYW0kUbSSCKlp9iSUAM0peVPxJrS8s81pR3fBvv73Ek+d4LOva6gTjjpFz4SlSJ/4p5CRAARwLIAy7KfO99DKftL0wDlpI5eL6cNUkmYGLhEU1o+SUALkDbSSBkpJI0kl2NeoNxJ2qM8+Uel1AnPcyf+wu958qZ6oApP/gDbjajsmBjUiKA3iKA3iBhi0C0dST2JjJlBxsg0fWuCpjR4NA88ypM/kXs1b/6EPl0SQERE88PEoAb5NB8igUj+uW7qSJtpZMwM0ka6YWY5FJ7wvZr3hJO+R3nyr2mKzaBERNXCxKAO+Dw++Dw+tMG+lbRpmchaWeiWjqyZhW7q0C0dltTGLWI9ygNN0/In/NyVfP6E75z0ecInIqo9TAzqkEfzoEVrQQtaTthuWiZMMWFYRj5RMMWEaZmwxIIFy04eBLBgJxHTzYyYOnBu6kC8fNO9Zg/GK7yyZ1M+EVF9Y2LQQDyaBx544PdwMRIiIpoftuUSERFRHhMDIiIiymNiQERERHlMDIiIiCiPiQERERHlMTEgIiKiPCYGRERElMfEgIiIiPKYGBAREVEeEwMiIiLKY2JAREREeUwMiIiIKE9Nd3e9eqeUOgrg9SKKRgCMl/DRc5Wf6fXpthezrfB54feLABwror7FKmU/FFO2nPthtn1Sy/thttdL/dvP9fx0EWmbq8JUdGyoVlyY6bVS/va1cjwUU76UuDDddu6HuZ+XJy6ISNN+AdhazvIzvT7d9mK2FT6f8v0et/ZDMWXLuR/m2Cc1ux9me73Uv32190Ozf1UrLpTjf6FWjoeF7Idit3M/VG8/NHtXwr+WufxMr0+3vZht/zrLa+VUymcXU7ac+2G2fVJu5dwPs71e6t++mOdUPtWKCzO9VsrfvlaOh2LKlxIXptvO/TD387Lsh4bsSmh0Sqk9IrLR7Xq4jfvBxv1AAP8PcrgfbAvZD95yV4aqYmvuG6VUFMAWAP0A+kVkn2u1qr6p+2EjgPUico97VXLFtPsBwD4R2elarajaGBdsjAu2eccFthjUOaXUF2D3K40ppe4WkTvdrpNblFI9AG5t8n2wBcBOEelXSu0QkSvcrhNVH+PCcYwLpccFthjUKKXUegAbRaQw67sBwBiAKI5fBZxXkAn3VL+mlVXCfmhoxe6H3OtO+X5XKksVw7hgY1ywVSouNPvgw5qklNoM4C7Yf9jcth4AV4jIThHZ7ryeey1XLooGUup+aFTz3A83A2jaK6RGxLhgY1ywVTIuMDGoQU7/z44pmzfDzgLznOzvOQAdzqYxNJAS90PDKnU/OFcMX8Lx/wtqAIwLNsYFWyXjArsS6kcUwHDB8xHYTYRbAWxRSo3B/qM3upn2wz7YB8V6pVSPiDR6M/q0+8G5YrgbdnPhPrDVoNExLtgYF2xliQtMDOqciIwBaLbRttNy+tG2zlmwgTnNh9vdrge5i3HhOMaF0uMCuxLqx9TmwA405+Ay7gcb9wMB/D/I4X6wlWU/MDGoH48B6C14Hm2GUbfT4H6wcT8QwP+DHO4HW1n2A9cxqEHOaNNbYfcX3Z1bjMLZntfoi9dwP9i4Hwjg/0EO94OtkvuBiQERERHlsSuBiIiI8pgYEBERUR4TAyIiIspjYkBERER5TAyIiIgoj4kBERER5TExICIiojwmBrRgSqkHlFJb3K4HEdUWxob6xAWOaMGc+76PAog5N28hImJsqFNsMaAFcw74nQBucrsuRFQ7GBvqExMDKpdtsNftJiIqxNhQZ9iVQGVR0GTYKyLNeLtTIpoGY0P9YYsBlUVBk+ENbteFiGoHY0P9YWJAZeHc6rMfbDIkogKMDfWHiQEtmHPgrxeRWwF0KKXWu10nInIfY0N9YmJAC+Ic+DeKyD3OpscA3OxilYioBjA21C8OPqR5cw78O0XkioJt6wHsEpGYezUjIjcxNtQ3thjQQtwN4MbCDSKyD0C/UooDjYiaF2NDHWOLAREREeWxxYCIiIjymBgQERFRHhMDIiIiymNiQERERHlMDIiIiCiPiQERERHlMTEgIiKiPCYGRERElMfEgIiIiPL+P019fx43oboeAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 540x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pylab as pl\n",
    "\n",
    "colors = pl.cm.Reds(np.linspace(0, 1, len(variance.values())))\n",
    "\n",
    "f = plt.figure(figsize=[7.5, 3])\n",
    "\n",
    "\n",
    "plt.title(\"\")\n",
    "\n",
    "ax = plt.subplot(1, 2, 1)\n",
    "#plt.errorbar(lambdas, list(mse.values()), yerr = [np.sqrt(a) for a in list(variance.values())], label='clean')\n",
    "#plt.errorbar(lambdas, list(mse2.values()), yerr = [np.sqrt(a) for a in list(variance2.values())], label='adversarial')\n",
    "errorfill(lambdas,np.array(list(mse.values())), \n",
    "          np.array([np.sqrt(a) / np.sqrt(8) for a in list(variance.values())]), label='clean', color='#2ca02c')\n",
    "\n",
    "plt.ylabel('test error', fontsize=16)\n",
    "plt.xlabel(r'$\\lambda$', fontsize=16)\n",
    "plt.xscale('log')\n",
    "plt.ylim(0.002, 0.03)\n",
    "plt.xlim(0, 100)\n",
    "plt.title('clean')\n",
    "#plt.vlines(2.1333, -0.01, 0.016, linestyles='dashed', color='#2ca02c')\n",
    "#plt.annotate(r'$\\lambda^*$', (1.95, 0.018))\n",
    "#plt.legend(frameon=False)\n",
    "\n",
    "ax = plt.subplot(1, 2, 2)\n",
    "errorfill(lambdas,np.array(list(mse2.values())), \n",
    "          np.array([np.sqrt(a) / np.sqrt(16 * 0.5)  for a in list(variance2.values())]), label='adversarial', color='#d62728')\n",
    "    \n",
    "\n",
    "plt.ylabel('test error', fontsize=16)\n",
    "plt.xlabel(r'$\\lambda$', fontsize=16)\n",
    "plt.xscale('log')\n",
    "plt.ylim(0.002, 0.03)\n",
    "plt.xlim(0, 100)\n",
    "plt.title('adversarial')\n",
    "#plt.vlines(1.1851, -0.01, 0.018, linestyles='dashed', color='#d62728')\n",
    "#plt.annotate(r'$\\lambda^*$', (1.05, 0.02))\n",
    "plt.tight_layout()\n",
    "#plt.legend(frameon=False)\n",
    "\n",
    "#f.savefig(\"pe_fig.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "mse = [[], [], []]\n",
    "var_ = [[], [], []]\n",
    "taus = np.arange(0.0001, 1, 0.05)\n",
    "for t in taus:\n",
    "    plotting = Simulator(tau = t) \n",
    "    mse_all, variance_all = plotting.run(lambdas = [0, plotting.optimal_lambda, 1000])\n",
    "    mse[0].append(mse_all[0])\n",
    "    mse[1].append(mse_all[plotting.optimal_lambda])\n",
    "    mse[2].append(mse_all[1000])\n",
    "    var_[0].append(variance_all[0])\n",
    "    var_[1].append(variance_all[plotting.optimal_lambda])\n",
    "    var_[2].append(variance_all[1000])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "mse2 = [[], [], []]\n",
    "var_2 = [[], [], []]\n",
    "\n",
    "for t in taus:\n",
    "    plotting = Simulator(tau = t, num_mali = 3, tau_a=0.1, theta=0, sigma=0.5, K=10, n=10, num_trial=2000) \n",
    "    mse_all2, variance_all2 = plotting.run(lambdas = [0, plotting.optimal_lambda, 1000])\n",
    "    mse2[0].append(mse_all2[0])\n",
    "    mse2[1].append(mse_all2[plotting.optimal_lambda])\n",
    "    mse2[2].append(mse_all2[1000])\n",
    "    var_2[0].append(variance_all2[0])\n",
    "    var_2[1].append(variance_all2[plotting.optimal_lambda])\n",
    "    var_2[2].append(variance_all2[1000])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.0047908706110883345, 0.004572787437747128, 0.004598569454165577, 0.004770618076253796, 0.004583165623642554, 0.004950057301709997, 0.004579454900567387, 0.004748955971692328, 0.004754136956152046, 0.004702054178998425, 0.004859857264197733, 0.00487910056548816, 0.00472463159759325, 0.004840676338118819, 0.004900521205953544, 0.004692821184486249, 0.004927749494060927, 0.004988820789018073, 0.004927956363230694, 0.0047553294234049725]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgwAAADICAYAAACEYLmEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO29e3wb13Xg/70zAN+kQNIPWbKjmHT8iK3EJuUkdpo0icg8um2epNXutt3utibbbffZ1oqTtEkfsUt1++t22/4SSn0ku203EpU2DzdpQ7pxk22cxCRtJ87Ljui3JVsSBfFNApi7f9wZYACCBEAABAY8388HnJk7d2YOhrhnzpx77rlKa40gCIIgCMJmWJUWQBAEQRCE6kcMBkEQBEEQciIGgyAIgiAIORGDQRAEQRCEnIjBIAiCIAhCTsRgEARBEAQhJ2IwCDlRSk0ppYYqLYcgCJVBKXWXUmq80nLkQil1Sik1kmdd0WsFEqq0AIIgCIJQIoaBmUoLUauIwSAIgiDUBFrriUrLUMtIl4QgCIIgCDkRg0FIopQaVUpdcD/jSqmuDepFlFJjbr1T/n5ApdSQW+bt6/Ptm/L6Qt39UxtdQxCE8pKjrXa57fOCG7twq2/fqFJqLONcyRiHHPphSik14NM1EXf9lE8n9OSSb5NzTak89ZFQOGIwCIBpfO7q1VrrdmCzwKH7geNuvV5g2PfgnwV6fecYyzj2bmDQ3U+O6wiCUD42a6vjwKTWul1r3Q9EfPtGgYGMcw275bC5fgA45i6vBvqAA1rrbrf+IKkYhFy6JO1cWutogd9RKBAxGARcq7tLaz3sNTqt9YTWel3wkFu3BxhRSp0CpjDKpM897qTvHEeBiPfG4HLC17CPk66IBEHYJjZqq3594Ks+7jtuGphRSg0AuO27Q2t9Mpd+cJn06Zoo0OV6CiJa6xmfTLl0Sea58v6OW7lfggQ9CoYu8o8s7gIm3LeOdSilIhgvwkauv1MZ2x15XlcQhBKySVvtAqZzHD6K8SqcdJcnfMduqB9c/MbHhFLqXszbf5dSagLjgYzmoUvSzpWNPM8h5Il4GAQwxkK+sQQzwIFsO9zG+STwkNa6V2vdWyL5BEEoITnaaj764CjQ557nDlLdERvqBx9p3gCt9RGtdTfQjnmBGCpAl2T1LIDoo3IgBoPgDUWaUUp5jR7XRbjOKvfVHcuo24Np7BHctxPPZSkIQtWxYVv1hiZ6+sB98Pq7J3Dd/BMYz8Cs202RSz+sQynV5+kZ95yzueQrxXcUtoYYDILHQQBvlARwiI27KbLWdWMeTgKn3P7LWzc4XhCECpJHW+0FDrjt+xjGOMhkFBgi5V3wKESXABz21Y0CR0uhS0QflR6lta60DIIgCIIgVDniYRAEQRAEISdiMAiCIAiCkBMxGARBEARByIkYDIIgCIIg5EQMBkEQBEEQchLITI+XdHbol1+1t9JibIxSYNdVWgphhzM1NXVOa31ppeWoFvLWG/OnYTkKHV0Qqi+/YOXEssAKV1oKIUBspjcCaTC8/Kq9TE78XaXF2JzGdmjurLQUwg5GKfV0pWWoJvLWG8sX4K8HYNdl8L4/B8suv3DlZNeVEG6otBRCQNhMb0iXRLlYiUJ8rdJSCIJQKI3t8IZfhRe/DY+drLQ0xbN4ttISCDWCGAzlQmtYfKnSUgiCsBWufQe87HZ48E9MF0WQia/CysVKSyHUAGIwlJPYijRUQQgiSsGb7gbtwAP3mheAILN0HpxEpaUQAo4YDOVGGqogBJO2PfC6X4an/wWe+FKlpSkOx4Gl2dz1BGETxGAoN44Di+cqLYUgCFvhVYfgshvhq79vRk4EmZWLpntCELaIGAzbweo8rC1VWgpBEArFsuEtvwGrc/Avf1hpaYpHAiCFIhCDYbtYPBv8flBB2Ilc8gro+Tn4/n3wzNcrLU1xxFZgZa7SUggBRQyG7SIRM+O7BUEIHgd+HiL74IF7ILZcaWmKY+m86SoVhAIRg2E7Wb4guRkqRDQapbu7u2rPJ1Q5oXp484dg7nn45milpSkOJwHLEgCZD6I30hGDYTvRWvoQBSGo7O2BG98Lj/w1vPS9SktTHCsX5eVFKBgxGLab2LL0IQpCULn9P0FjB/zT75huxqAiLy/CFhCDoRIsnZPcDBXm6NGj9Pb20tvby8TERFp5d3c33d3dHD58GIDBwUF6e3vp7+8nGg340LoqQyk1oJTqc5c9m9TrUUoNZZSNKKW6lFKRzH1lo74VfvQwnPuB8TQEmdiyGcEl5M1O1xuBnHwq8DiOCTxquazSkpSV33jiOR5bKG+A2E0tjfzOK64s6Jjp6WnGxsaYmpoCoLe3l/vvv5+ZmRlGR0c5deoUADMzMwCMjY0BcPLkSe69915GRkZK+A12LkqpLqBfaz3sbo8Bg1nq9QHDwEMZu3qAcWAauLO80vrofgt0vQW+edSsR162bZcuOYvnoK7FZLasIqpRd4jeCKrBEPQoZTDdEvWtEG6stCQ7juPHjzM8PJzcPnToECdOnODUqVMcOnQoWd7V1QXAxMQEo6OjzMzMJMuEktAHpL16KaV6tNbT/jKt9YRrXEQyjh/VWldmdqg33gV/MwBf/ii8++NV98DNGydhMkDKzLo5Eb0RVINh7jkz4qCxvdKSFMfCS+btJKjKJgeFvvlvF5nuwfPnzxOJZD6LDNPT04yMjDA2NsbMzAz33nvvdoi4U4gA533bs0AXxmOQD12u9yECsK3GQ8ulcPt/hgc+Ct/7HLzyXdt26ZKzEoWGNrDDlZYkSTXqDtEbQY1hcBLwpQ+CE6+0JMUhuRkqwuDgIKOjZmhcNBrl5MmT9PX1cejQIY4fP56sNz09zeTkJD09PUQiESYnJyslspAFrfURrfWEayjcrZTKrr3LxY3vhj09JgNkkNO/SwBkXojeCKqHoXU3PPsN+Pr/b6KWg8zyBdM1UUXWfa3T19fH9PQ03d3dRCIRRkZGki7D4eHh5Djp4eFhhoaGkgFOteJWrCKipHczdAAz+RyolBoAurTWR9yirN4JNxhyCOBlV+4pVt4MISx48wfhUz9l5pp4e4D7qNeWYG0R6porLUnVInoDlA5guuIDN+/Xkx+6FR77NLz99+Ca/kqLVBzhRti1t9JSCDWGUmpKa32g0nJshOsRGPEFPY5rrbM2ZvfBH/EMBG9EhRfv4H7X3s2ud+Dm/Xpy4u9K+RUMk38BX/9T+LE/gK43lf7824UdMtksa7SLVMiPzfRGMLskAN7w67D7VXD/b8H5U5WWpjhkeJOwA9FaR4Exd1hlH5B8RVdKTXldDO6+fqDfXfcMhS53OOZdwOHt/wYut/wMdF4D//x7wW7Hibh0kQqbElyDwQ7D249AuAm+8GvBbqhg+hAlv7uww3BjEJIfX3mva1B4dQa11v0ZdU66nyP+8m3HDpsZLZfOB39Gy+ULxnAQhCwE12AAE6n89hGYfx7GfwN0gB+4jmMSOgmCEDwuvwlu+Vn47mfhiS9VWpqto7XoIWFDgm0wAOy5BX7kV+Gpr8JDf1ZpaYpjZc5MPysIQvB47S/C5fvhy79rJqkKKqsLJghSEDIIvsEAsP8OuO5fmcxrT3210tIUx+JLxsoXBCFY2GF42z2Agn/8QLDnmlg8K3pIWEdtGAxKwZs/AJdeC1/6EESfqbREWye+JoFHghBU2vbAWz4ELz4G3/hYpaXZOomYSegkCD7KbjDkM7mMu2/AnUxmYEsXCjXAO/47WLYJggyyS235QrDfTgRhJ3NNv5kGe/qT8MyDlZZm6yzNSgCkkEZZDQbf5DLJbGxZ6vQAUTfa+TBwbMsZ29r2wNvuhQtPwj/9dnBdalrD3Asyo2UNMD09zZEjR3JXzCAajTI4uG4eJiEo/Mh/g45umPiwGT0RRLQOruwBp1r1Rrk9DFknl8mo04UZY+3hZWzbGle9Fm77FfjhODz8v7d8moqTiLlGQ4BHflQR7e3tyalme3t70yaRmZmZob29PPOS3Hnnndx1110FHxeJRDh06NCWlIZQBYQbTTzD6gKM/2ZwR3Ctzu/oLlLRG+mU22DYaHKZJD7Pgpf5rSNztrqCueVnofsgPPjH8Ow3izpVRYmvwvzp4HpKqoz777+f8fFxpqamknPag5ld7sknn0zWO3r0aNpxmdv5cvLkybRZ7Pzn2+ic09Opn/7AwEBajnohYHReA2/4VXj268F+eVk8DysXKy1FxRC9kaLagh5HgKzpXZVSQ0qpSaXU5Nnzs5ufRSk4+GFofzn84/th7nTpJd0uYsuw8GKlpag5hoaG6OjoYGLC5PvxzzrnTTCz0Xa+jI6O0tfXl/XaIyPp8w5Eo1GOHj3K8ePHkzKBUUp+ZSAEjBvfa15evv6ncObblZZm6yycNcO+dzg7XW+Ue/KpvCeXcYMdR7XWWfdrrY8CR8HkhM955bpmEwQ59jPwD4fhfX8Gdl2h8lcHqwugXoKWyyotSUGcueceVr/3/bJeo/6G69n9gQ9s6djBwUHGxsY4cOAAvb29nDp1iuHhYaanp+nt7aWvr49oNJq2PTIywtGjR5PKYGRkJGvjBpIz1mVjYGCAiYmJdcdmTqHb1dXFxMTEhucRqhyl4M0fgpe+a4Za/uTfmMnmgsjCS+b7bIP81aw7drLeKLfBcAJffnjM5DHrzB43P/y01nrG1y2R16x1m9K+D/p+y4ya+OofwJvWxVwGh5U5UDY0d1Zakpqhq6uLsbGxtLLR0VEmJiaYmppKlvm3p6enGRsbS2739vZy//33p71p5MPw8DCHDx9ONvxIJMLQ0FByylyP7u5uTp0K+FwpO52GNnjrPfC3vwAP3GPWgzrB08JLgIL6lkpLUjF2st4oq8GgtY4qpca8CWPImFwGOIiJaRgFoso0oi6tdekiSbrebGIaHv5fZrKq6/9VyU697SxfAMuCxvIE2pSarb75bxczMzMFTz17/PjxtMCnQ4cOceLECYaGhtbV7ejo2PA83htINBpNUxoDA+mjijs6OtKUkBBQrniVyQT59T81gdmvfHelJdoaWpsuUmVBXVPZLlPNumMn642yxzDkmlxGaz2tte52t3tLaix43PbLsKcHHvgonHui5KffVhbPS19iiRgfH6e/v7Cp0TNdf+fPbzzsbHY2e6zNxMQEs7OzDA8Pc+LEiU2vNzs7W/BbiFCl9P4cXPka+MoRmC3egVoxtDbB2LHlSktSEXay3qi2oMfyYIVMfoa6FvjiXbC2UGmJimPxLKwtVlqKQHP06FFmZmbWWeZg3HzRaDTZyP3bg4ODyX7IaDS6zhXoJ9ubwsTEBDMzM/T19TE0NJQzMCoajdLZKd1QNYGyoP+3zQy7/3g3xAM8b4yXK2aHzX2z0/XGzjAYAJovMTNbzj0PE78V7KGKWsP8mR1r4W+VgwcPJsdTT01Nbeiy6+vro7e3l8OHD6/b7uvro7+/n+7ubg4ePMjIyMiG7smuri5mZlJvkjMzM4yPjyfdkJFIZF2dTB566KENFYsQQJovhb7fhvM/hP8b8KmwtTb6NL5aaUnKiuiNFEoH8MF54Ob9enLi77Z28MN/Zeasf/1/gVt+prSCbTdKwa4rIVRfaUmELJw8eZKZmZktJWDx8JTUVlBKTWmtD2z54jVGUXqj1PzfP4RH/gre8fvQ/ZZKS1MclgVtV0IooKPQqoxq1hs7x8PgcfO/MeOiv/bH8HzAx7d7bkGZd6IqKTaBykYJXIQa4LZfgcteaVLYBzlPDJhstHPPm4nzhKKpZr2x8wwGpeDgb5o38398v4kHCDJOwjRWmSSmKjl27FjSRVkI0WiU8fHxot4yhCrGmwrbceBLHwQn4O03qYfk5aUUVKve2HldEh7nfwhj/xYuuwHe9THTgINMqM64Ba2dZwMK2ZEuiXSqqkvC4/F/MAbDgZ+H1/2HSktTPHbI6CG73Cl+hHIhXRLZ6LwG3vxBeOFhMzY66MTXYP6FYAdzCsJO49q3ww3vhMm/gGe/UWlpiicRN54GmWm3Jtm5BgPAdT8G+wfNxDCn/qnS0hRPbEUmqxIChVJqQCnV5y43zGOrlOpRSg1t5diq5413mXlv/v6/wg/vr7Q0xZOIidFQo+RtMCil2sopSMX4kf8Gl98IEx+BC09XWpriWVty07cKQmkoV9tXSnUB/W5St5NA1tztbqbYu/HNS5PvsYEg3AjvOQqd18I/3AVTnwi+0R9fc42GgE7rLWSlEA9DDZi+WbDrTH4GOwRf/PXayG2wOg+L5yotRSCIRqN0d3cXXaeY+gGgXG2/DzNBXZJsngI3Q+z4Vo4NDE0d8J6PwyveCg/+MXz5d4IfQFjDRsNO1RuFGAy/p5Q6rpQK6FRrm9B6Bbz1oyZd6wP3BN+6B1iOwlKOacAFIT/K1fYjgD9H7ixmbplyH1udhBqMHrr1F+C7n4XP/8fgp4GPr0psVQ1RiMFwDBgELiqlHPdTO51UL7sNXjMEP/gCfOfTlZamNCzNwsrFSkshBJ/abvvVhLLgtb9kZtl94WE4+XNw8dlKS1UcsRWTL0aMhsCTt8Ggte7QWlsZH7ucwm07t/4C7Hs9fOW/w+lHKy1NaVg4a7ooBI4cOUJ3dzf9/f0MDg5y5MiRdXWOHj1Kb28vvb29TEwk50pLTvrS3d2dNuvc4OAgvb299Pf3r5tgplYoY9vPvGEdQL6zMuV1rFJqSCk1qZSaPHs+QB6363/cDPdejprh3y88XGmJiiO2bNLZB9BoEL2RoqBREkqpX1NKPeF+7imXUBVDWSbPe8tl8JlfhO9+ptISlYaFl0ww5A7GS2hy6tQpIpEId99997rkJv4566empjh8+HBaYx4ZGeHUqVNMTk4mlYJXf3h4mHvvvXdbv9N2Uqa2fwLwd9pGtNb5pl/N61it9VGt9QGt9YFLOzeeNrgq2dsDg5+A+l3wmV+CH3yx0hIVx9qimRo7QIjeSCfv7BpKqXsBBRwAOoEhpdTHtNa/VC7hKkJjBAY/aZKp/NPvwOlvwY/eZfoXg4o3HW3LZVC/jSEoX3w/nPl2ea+xez+84/dyVvPPYd/f38/k5CQ9PekxchvNWX/HHXfQ0dGRnC52eHiY8fFx+vr6mJiYYHR0NO38tUa52r7WOqqUGnNHQQCM+K45BRx06/QB/UBEKTXtjozY8NiaIvIyYzR84ddh/EOme+LWO03G2iCyugC8CK2Xb16vSnSH6I10CvEw9Gmt36+1vqi1ntFavx+jQGqPxnb4iT822de+91k4+e+D34+oNcy/uGNHT/T09DA5OQnA1NQUBw6s/+kWMmd9Z2cn09PTjIyMcOzYMY4dO1ZagauLsrV99+Gf/PjKe7XWUV+dQa11f0adrMfWHA274F1/aropvjlqDIcgzxC5Oh+Yod+iN9IpJH9nmknrjs0OqJmbB5ZtUrXu3g/jvwnHf9oEInW9qdKSFcdy1Cib1t3mO5aTPN78twuvUff399Pf37/uLQFMv+LIyAgDAwPJOevHx81ovtnZVP/32NgYIyMjybeNSCSSVCo1ys5q+9WIHYaDH4HIPpOZdv40/NgfmJebILIyZ7qAmy/Jvr9KdIfojXQKMRh+Tyn1EOBNozUMFD47RtB4+Rvg0F/DF++CL/wq9PxbY0hYAc6VHluG6DNmOGk4wF0tBTI9bbq4JycnOXXqFKOjo2n7+/r6mJ6epru7m0gkkpyz3lMaw8PDTE5O0tfXR09PD11dXckgpyC5FbfAzmz71YZScODfm4nzJj5sgiF/4o+g/epKS7Y1lqOAgubOSkuyKaI3UuQ9+ZT7VnEJ4KVnPa61rkjobkUmkYmvwlf/wAy53NsLb7sXmqr7h54TpYyF37Cr0pKUncOHD9Pd3c3Q0BDRaJTDhw/T39/PwMBApUUrG6WafKqa2n4xVOXkU1vlzLfh7/8bODF4xxG48jWVlmjrNHWYTxUieiOdgjI9ev2X7idwCqMoQvXw5g+YbokXH4NP/evgD3XS2gy7XHgpkMOdCqWjwyilSCRCJBJJcxcKm7Kz2341snu/Cc5uvhQ+9ysm0VNQWZo18VVVOveE6I0UkumxUK7/cRj4pMn//nfD8MhfB/9huzIHF58zM83VKCMjIxw/fjw5VhpgaGgox1GCi7T9aqRtD7zvL2DvAfin34av/U/QAU3DvDoP0aerLtGc6I10CumSmCU1+YvGBD3pSiRvqgrX4uo83P9bMPNl6D4IB38T6loqK1OxWBa07Ia6pkpLIpSAEnZJVE3bL4aq0BvlIBGDr/y+6S7tegv0fTjYuijcAM2XQaiu0pLsSErVJXG1P8tbTWZ6LIT6VnjH78Pr/wvMPAAnfhbO/7DSUhWH45jo6+ULlZZEqC6k7VczdhjedLeZeffJfzbdpeXOYVBOYitmGPvi+eB7b2uMQgyG2h3nvFWUglt+Bt79cVhbMFHLP/hCpaUqDq1NQ50/U5OzzAlbQtp+taMU3Pxv4L3HTLfEp38eJv+iauMCcqK1eXGJPm0yRApVQSBjGBw0q47DmnaIaU0CjUMFLdG9PXDob+CyV8L4b5gMkdv4I9do4lqzph1WHIclJ8FCIsF8PM7FeJzZWJzZWIzZWIwLMVO2kEiwmEiw5CSS9zKuNdq7j6sLMPdc1U+x62hNQmvijmbNcVh1P2uOQ8wx5XHH1ElojaM1+XbDCUmqpu0LObji1fCT/weuOWjyNXz2PwQmSVJWEnGYO21eYGo4xiooBDKG4VWvulHf94UTWfcpFEoZ4RRguSlUN8oyk1nu39YZS289We7eu+R+J07r1DFaH/0rEi27if7oh4jv6Ume0yIllxFLJdeT+5RK1tfgPuSMUeBocHA/WuNgDHFdYmPJu4cWoCwL1XI5qq4FS4GFwlJgK2X2K7B9ZSRlMoadf53kurl3qXXvvuoNjk3V1xn1S/FdzTKVbVeRsfT9hlRaue9+ZSn37/PK/ef2nz+fLEiF1mkJhySGwUfNxjBkQ2v43ufhq0fAqjMxVl1vqrRUxaGUGcreGMldV9gym8Uw5J19SGtdnQNlM/AeOqmCbXybVDbRA7/I0pW30/mV36Hzvl9mfv9PEu0dMsMyA4J3Dx2ARAIuvkCisR2nMRA/gYLwjA6d/JOt0s72SASl7Qs+lIJXvtN4HL70AZN0bv+gibkK6rw4WpvU9qvzZl6cAOnUWqHQ2SrblFI3+9bb8jhmQCnV5y7X59VM1etRStXEeJW13a/izHs+ycIN76bt2/+H3Z/594TP/aDSYhWFvXwBe/50cPtEhaLYStsXqoD2fTDwl3DzT8O3x2ojODu+CtFnjfEgcVbbSt4Gg1LqvcA/Ad5sGd2+9Y2O6QL63clhTgJ3b1Cvz91XM74mHW7iwut/nZfe/v9hrc2z+7O/QNvDfwlOcPvhrNgSobnngj3xjVAwW2n7QhVh18GP/Fczod5KFE78DHzrRPA9Z8tRExS5ulBpSXYMhXgYRtx+jQsAbra3XImw+4C0qbyyeRncmebGC5AlMKxc+TpOv++vWOo6SGTqGJd/fphQ9KlKi7VllBMnPPc8anW+0qII28dW2r5Qbey73QREXnkrfGXEdFMEfQi1kzABkXMvSFDkNlCIweCPxUMptYvccVgRwD/X5yw7UNHo+jbOv/kjnH3L7xKae4Hdf/dztDx2IrhZ2dCEFl/CXjwb/LcUIR+20vaFaqSpE378f8CP/Co8/TX41E/Bcw9VWqriWVsy3oblC6KTykghBsP7lVKTQJdS6teASeCe8ohVmyx3vYXT7/srVvfcSsfX/weXfeE/mbiAgGKtzmHPvyBxDbWPtP1aQllw87+GwU9AuAk+80vw4J9U/RDqnHg5ZC4+a5I/CSUnb4PBjUG4EziKebt4q9b6b3McFs3Y7gBmCpLQRSk1pJSaVEpNzp49h37hBfTZs+hoFL2wgF5ZQcer3yXlNHVy9q1HOP+Gu6k7932u+Nufofnxvw+sVWzFVwjNPYuKSwOtVbbY9oVq59Lr4dBfwyvfBVN/CX/78+ZhG3Tia2ZunIWXJCiyxOQ9rBKSfZeFzFR3AhjxbUe01tOFXNN37aMYhcWrbrjWjMSPrZmPvx4KwmEIhdI/to2yChoUUj6UYvG6n2BlTy+d//y7dH7lozQ+9c/M/sj7cap0mtfNUE6C0NwLJJo6cXbAVNk7kS20fSEIhBvhLb8BL7sNvvy78Kl/A296P1z3Y5WWrHhW5kwG3uZLTSp/oWjK+gTVWkeBMXdYZR8+40EpNaWUirjrfUA/0O+uF3NVY0QsL8H8HFyYhbMvwZnTxivx0kvo2fPoi65nYnkZHYuhE9vvVk+07uGlf/UnXHjtf6Tx+W9yxd/+NI1PPrDtcpQGjb10jtDcc6jVhcB6TARhR3JNnwmIvOQVJlvt+G/WRkpmxzFTZ198PvhdLlVAQR6GreCOgMhW3ptRZxvy1WuIx8wn616V9EYQsiEUdtddD4UqQ5yXspjf/1MsX/k6Ov/5d7j0/g+weM3bmb3tv6IDaBWr+Cqh+IvoZRunfhdOfRtYgUoIKAhFkcCkIdfarDvaZCpNgJuaHLNOKmOr4xrYycyjvmy1ys2qCqkyL8uqd8xG2W3dhLJp2UiVf9t/ndYr4D2jZg6Kh47BmUfhrffA5TeW9X5tFcdNu2aSzLkJ2NyEc9pNK+t1SOj4PHplgURDO05DO1qptMyyXvrSje6Pud8que7dUwv/PTQZb/3/t1qjKINBKdWmtZ4rlTCVx2dQZEk1oC07rYvDb1QU290Rb7+aF995lF0Pf4K2Rz5Jw/MPMbf/J1m4/l3oAE5Vq5wE9vIs9vIFnPoWEg0RMx5cqAmC3PaTD/FkunU3s2ny4e3LAOp7oKz76FRy8lTK8uI8a/7rpgq3z1unUPCqn6Ph8lto//JHsE/+Oxb3H2LxpjtwWnanDY3JlW5/o7pp676N1MM++3327jEUcZ9jZ9FLURJNl6DDTVs7RwFYGcac3+jwGyZWRpp5v+FhKbXueL8xkw/53q1cP7W8DQal1ENa61sziseAt+V7jsDjJGAtAWvrrQmtLJ8hETJxFOEwKlSATWaFuNj7Cyy/7PVEHvoY7d/8U3Y9/JcsXPdO5m+6g0TL7hJ+mU9d+cQAACAASURBVO1CY63OY63O44SbcOrb0HXNlRZKKIBytn2l1AAmODoCzGSLcdqojlJqBBjFDNe+w41z2pA1R/PMiiQd2wzvQbx8+atYec8n6Hjwf9D87U/R/O1PsbzvjczfdAerl786/UkfMFQiRmj+NE5dK4mmzrJ6QL05dJIEvKs259PMTQHb7a7f7NvVCWQqkZ2LdtwgzIxiMF4Iz4hwl5sZEmuX3sBLP/Y/CZ/7AW3f/j+0fmeM1u+MsdR1kLn9P0XskuvK+lXKhRVbwootoe0wTkMEp6410Iqn1il32/dlgh12t8eAwQLq9GASvk1jRnHkINjKervR9W2cf9NvEj0wROt3P03zDz5H01MPsNZ5LfM3DrLY1Rfo+RystXlUbAmnqcN0nQo5yef1tx84hFEcI6S8IFHgcJnkqi28bo6V5WRR0pDwj+jIMCRil1zH+Td/hOitv0jrd8Zo+f5naT71JVau6GHuVf+alStfZ8ZUBwyViGEvnsVamsWpbzUjK6yyh9MIhVPutp81E2yGl2GzOqPukE+hjCRadhN9zS9zsefnafrhP9D6nTE6v/JRIt/8UxaufzcLN7yHRPOllRZzSyidMLpodQ6nrhWnrll00SbkvDNa608Dn1ZKfVxr/YvbINPOIUsA5npDwiZe38mF1/wKF2/5d7R8/3O0Pnacy/7x11iLXM38/p9i8Zq3BjI+QOkE9koUe+UiTl0zTsMudFBn0qtBtqHtb5QJdjrPOl3uqKqIK+/mxkMigZ497+sM1+kd447jBiK4+5RlPGCW5Yts85cB+Pcrt/PZMmVeHatMAdPbjA41sHj9u1m87l3UvzBF63dO0PbIJ2l79H+z1PUW5m+8g7XLqjNAMhcqvoodX8VeOocTakDXtYjxkIVC7sbHvRWl1Fswbx/3BjXwqarZYCSHY4eY2/s25q56K80vfI2274/R+dV7iEx+nPkbB1m4/j04DUF0rWmstQWstQWcUANO/S50ffACPWuYqmz7Wusj3ro7THvCHcqNr3wIGALYu/syWCkgwZh2jEFRgkymJqe2lTIuLMtnXKiMfSqLsWLWq8LwUIrVvQdY3XuA0NxztHz307T84D6aT42zeukrmb/xDpaufjPY4UpLuiWs+ArEV8R4yEIhd2BMKdWLcU+OAkeA+5E4hu0jEU9OsLLY3svi63poOP8YbU9+jsjkKG0Pf4KFq9/G/CsHiLfvg1CoOhRMAVjxFaz4Cnr5PE6DOywzgN0uNUa52r4XyOiRLRNs1jpuIGSXz2jI5p3IkvCtgmgHEsVlHlxveGQxLkIhs27bYFsou3wPunjblURf95+52PMLND/xRVq/e5JLHvgI8W/+CQs3vIeF69+F0xi8ZHQeYjykU8i3vqi1nnMt9hGt9Z8ppYbLJZiQB0qxcsl+Vi7ZT3juadqevI/WU39P6w8/z9Lu1zLX9U7WOm9Idm0QroO6usJGblQI5cSxl85jLV9A17WQaNgVyG6XGqFcbT+fTLBZ67iGsN+46NhqFtnAUaDhocGMBAjZZmlbYIfM0ldezMuFrmtm4cYBFl75Xhqe+wat3zlBZOoYux7+BIvd/czfOBjYYG0PMR4KMxhm3YlnDpHegIUqINa2j/Ov/mWi1/0UrU99kZZnxmk+83VW2q9nruudLF/Wm3xT1yioc+Mk6uogXL1GhNIOanXOBCWFm0ycwzaMnxbSKEvb11pHlVJjvuyuaZlggYMb1XGNhgF3FEUXEoC9Od6Q8E3Qlp30SmB7xkUq90xeuWaUxcpVt7Fy1W2Eok/R+p2TND/xRVqe+AIrl7+a+ZsGWd73xsA/ZHeq8aB0nuNClVJXA8PAuNb6fne7xw2M2lZedcO1+vOf/JPtvmygUPFlWp77Mm1P3kdo+SyxpiuYf1k/S1e8jkRjlohmZSVzR1BXB3Xhsroyi0Hbdaa7QoZlbso1zY1TWusDxZ6nmtp+MYjeKBIvcV04lBoqHgqh7M3zGKjVeVoev4/W736a0PwLxJsvZ/6V72XxunfW3NwztWA8bKY38jYYIDkuu0tr/YhSahegKxH4JA2/AJwETS9+g7Yn76M++gQAq5FXsLj7to2NBw9luR6IcHKZSzlsJ1rZOA1tbvrpYDbOclIqgwGqp+0Xg+iNMqGsdUYE4dD6Fw4nQeOz/0LrYydoOD2NY9ezdM3bTHdFR3dlZC8jQTUeNtMbhWR6fC/wAUyX2K0YN+D7MW5KoVqxbJauuJ2lK24ntHiapjNfp+n0g3R8/3/R8f3/ZYyHK25naffrSDRekn6sdmB1xXy8Ist2uzNShkSlZgFVOoG9fAF7OYrjxTkEOJFMtSJtX9gU7cDamvn4i7PMHLy093aW972R8OwpWr9zgqYf/gMtP/gcq5fdxPJVt7Gy9zWsXXJ9Tcw/U4vdFoV0STyhtX6FUupLWuu3umXZUsaWHXlTKJ7Q4mmaTn+d5jNfo27uKQBWItca42L360g0duZ/MjuU8kJs9HaxTehQPU6oER1qQIcbd/QIixJ2SVRN2y8G0RtVhOuNsBJLtDz9JZqe/Sp15x9HoUnUt7Gy51ZWrnwtK1e+hkTzZZWWtqSkjIeWqjSMSuJhIJXlzYzsMW7JinQgJ4A57c4UBljGlsVyP5mzudUq3qx3ccBBkcBsJ7QijjcjXqpcAWEgjCbUtIdw93sJd7+HpsXTtJ55kOYzD9LxvU/Q8b1PsNJ+HUu7b8vPePCGe/ozWebrpszre5qZ/uLa/bjrCd+690mwDCxjuTPHYYVQoToINaDsepQdTk7sYmVM/uKfbc6bNMb7HVkZ+1Cp31rQhq5ugapp+0KN4OaacVDM7Xkbc3vehrV6kYbzj9F4/lEaTj9M85P3A7C2ax8re25l+arXsbanFx1wL2J2z0N1Gg+ZFKK936+UmgR2uRHTw1QoMvlxJ0T/Un6pSC00IcB2lyHAVt66xnbLQspbN/scIOF72Ma1Sj6cE6QeyF65v8xc12e4sN6wSVsqva6MpAzGAMhmGDil0tnqMrji1YSvGOYVS8/y42cf4MfOPsD1rvHwrbYbeeCyH+VfLnsjF+o7jdGhNGH3fnnbIVf2uHtf4suk7hNxEjpBnFXiyiKuFHHcj1KbGgLFpc2JAytA+brb/QZqmqHhNzJ8s8ytM06USvtPZs72l7JHVPb9GQeW4UleNW1fqF2c+l0s7Xk9S3teD1oTnn+GxnOP0nD2EVp/8FnavncSbYVZ6Xgly5f3sHLFAWKRq4130xvNYQcrq2aa8RBuRIebq9p4KDTo8RZMbneAk1rrJ8siVQ6uvP56/R///JhJwgY4WiXX05cq+dCN+966497D3S3zGwFePUW6cWEnDQ9/ma/cNThs99hs8mht5MkmZ/K7eHWT508ZNlbGdbKVe7Imy337NBBzv18MRSzLdly75Sg6lp6j58WvcuClr/LyBTPk/Tu7buQrl/0oD1z6Rs7UX0LMrR93j3dcGTyDzDPGshlomYZcyFKEbIuQpbAti5BtEbYsQgps5R3nfgBbqeS6V+41s7T7q9PXM/cl7DDaCuFYYfNBmSmPk/V16jhvm9R0xt66Wfq309f9UyU7Wvtk075pe9OXJLfTpz1ev5+M/ZoTt7yilEGPVdH2i0G6JIKLSqxSP/tdGs8+SsO5R6hbeB6AeH07K5e8muVLb2blkv04dW6it1AoNTw0FDaJrNyEVtUUuL0RlTQeSjZKYt3BSrVVZJTEddfoz33yTyoWbLcTCS28QNOZB2k+/SB1808DsNJ+PUtXuN0WDWXM5hZKn6DLe6Mo15uEtuvQoQbjLgw1BDbFbSlHSWRSqbZfDGIw1A728jnX+/AoDee/hR1bRKNY29WVNCBWI6/YIMBQJb0RyY+X1CpkV91w8u02HkpiMGQLclJK/aPW+m0lkLEgbmpo1GMvfzk01KMaG6CxIX3Z0IBqrHeXDaiWJmhqQrU0oZqboKUZ1dxoIvwD5L6qFkILz9N85us0nf4adfPPoFGstl/H0hW3s3zpLcSbLt+eAJKM+Ahvu9T/U23ZaC+QMlQfmAmyShj0WDVtvxjEYKhRdIK6izM0nn2EhnOPUh99AqUdnFAjK503sXzJq1m59Gajl/LFzjQoQikvRQW7PbbDeCgq6NEdf93trt/s29VJheaRUB0RQu9+OyyvoFdWzHLZLJ0LF93yVVheSflosxGyobnZGBJNjaiWZmg2RoVqakgaHEnDxG+gNDRAY33VZkgsJ/GWvVy85n1cvOZ9xng4/SBNZx6k47t/AYATamStdR9rbfuItb7cXb4MbZc4WGmDSbq019i9t4iQnVQAW3FHKieBWluAtQVzfmUlDQfvU4ujMaqx7ReDA1zUKtmNFnO7ITfa9ndpakwXptdtmPqkl6W6OVPdnXFfWSJLmdc9Clm66pRO6+Kzc5R5XY/+ri0nYzv5ca+psxzj/yj3vBYm3iq57sqZWk/FY6XqrN+/0fHp3bUq1X3nkyWz+zm53nIjTstN6Kt/mnBsgd2z3+LK8w9z5bmH6XzxIQCijVfwTOctPNXZw7Md+1kLNaVdxzMBLEDFtBtz5GDhALGknACWpVC2hVIWyrKwbDO/h2XbZtsXs5SKaUrFMCWv4+7PjHfzjjH1VKpsdQnFEoqz6HA92A3Gi2KH0MpGZXhU/HbNZiZOvuZPTg+DUup9mPHWfcBDvnNHMZnfjuV5rZKR75uC1hpWVtGLS7C0jF5YQi8uwsISemkZvbBo9i0uoReX0AupdVZW8xMmHFrn0aCxwRggTY3Q1IhqdtcbNygLVX+fWj6EFp6nYfZ71M09RXj+Kermn8GKm5ETGkW8+QrW2vax1vpyYu4y0dCxzcNZVCpIKs2oCBVl/KW6MRrR4YaqGGtdrIehGtt+MYSve6Xu/PjfbNv1UnFPOiNuJ1t8j3lo+YOp0wyPrEHX4h3NidZ0Lz/Lmy48xJtmH+L10YdpclaIKZvJtpv4cvutPNDxGr7d8gozskvgxbfcUpIuiY9rrX+xpJJtke1wLWrHMcaG66nQyyuw4i1X05dJT8cqemUFvbRiDJSlZVha3tzLAVBflzIkfEvVYrwfSS9IS7PpTmltNl0rjQ3V3aWiHULLZwnPPUXd3NPUzT9FeO5pwssvJaskwq2uB2Ifa22uN6LlSrAqFDdge0lmUl6JpBuygJgZfzeGE2qoSEKpEnZJVE3bL4ZrL7lEf+wd78BKJLAdB8txsJ0Eyl23Egksx0E5DkprlGWhw6Fkv7YOu2mQw+7vIRRChd2l+5uxQjZWOIwVCmGFva4yGxUOpXWfqWQ3Wii1LxxKTgqVT7vWmqyB3Cidept16yrfJ31bZylL/6R5UHwelQTmjT+xyX6dY793vCZ9WLzKMnIs20izzNFm/rLUG3rqONuJ0XLh+0TOPcKuc4/SNG9id+PhVlZarmS1aTcrzbtZbbyClaYrWG66nES4OcOjodI8MN40YP7g9qzlaccrtG3jWBbYFo6y0bZl9IalzBKSgdKpc68/n//3sBFp9ZSFtlz/jmWhlQ3KBsvI8Z+u3lOeoMdKEaS+yKSXwzUetM+QyFq26JW5Ho+l5Y1PblumC6Wl2Xy8+IwWv4Hh29/aAq3NFY8SVrFF6uafcT0RT7vLZ7Ac072glU2sZS9rrfuItb3cGBKt+3DqK5x3Phl9nZqQx+vyyBUoVYlujHIGPQaRm1pb9ckDvcYn7BmArsImc10pMyNkIg6xOMQT6EQcYgmIxyEeR8dT68SLG/ybhlKuwWqnDArvNxd2f2shO2VgeAaJa8yo+nrTXVpfb7pRvWVDnfGC1tdDQx2qoR7q680Lyw4NILdWL9B47tvUn/8O4cXThJbOEFq9kFYnUddGrGk38ebdxJt2p607da3lE84/usMfT+Eaq+Vi31U3lSRxk7AFlFKpYMzO9oKP14kELLrdJ+6HxaXU+oK3voRzbhaeetYYGrH1fftJmptQba2othbzaW2BtpZUWWtqneamkisTHW5mteMGVjtu8H9RQotnqJt7ijrXiGg4/xgtL3w1WSVe3+52ZRhvRKx1H7HmPds37Eg7EFuDLLdWg2ngXsxEKGNsOA4qtgwxX3IrfzdGqM50Y4hbtGxYV15B48gHynJurTUkEq5xEUe7y+S2b51Y+nZqPZHVEFm/7TNS1pbR/rJYHL26BisrBU2BTV1dhpFR5xoZbldrU6OJ6/J7QRt9ZY2N0NRQ8ZeRQnHq21nc+0YW974xWabiy4SWXiK0dIbw4hlCS6cJL52h4fx3CD3/lbTjE+Fm4k2uIeEaEd66U7eruO5WLz5rbX33ePqU5ZbRG5YyZZbfALZBqZL9X8RgqHKUbbsP85aCjtOra+uMCz2/CHPz6LkF85lfwHn+DHpuARaXsvu0lDJeiWzGRKtrcLS1uuUtRnlspZEom3jLXuIte1ni9clia22OurmnTbeGa0i0nfs2Sps3OscKE2u5yhgSrici0dBJoqG99EGWufCCMLOEv2ivcdupN0HsNZS9TMj2j05Uxl2obHOMCiW3sWzjPvQtxcCoDpRSKa8T1ZEGU8fixru5upq+XPFvr5ntbPsWl3HOX/B1yeYR15XZvZoZz9XYAM2esZFRzwssr3A3qw41EmvbR6xtH5n+XZVYTRoTocUzhJfOEFo6Q130CZpOf83trDA4dgPxZtcj4XolPO9Eor69+Labx5Tlye8ErlGR+cliZGyCGAw1iqqvM28JeXo1dCIBC4tJQ8IYFca4YD617jz1HHp+YeOuEttygzrdoaveqJPmJnfd7KO5MaO8CVW3Pm7BqWtj5ZL9rFyy31cYI7zwfFpcRONLk7Q89+X0Y0ONJOojJOrbfUtvPUKiwaw74W2YJjvZuNfW7dLeuHBL+Rq1QinLdZ2r9Mae9sbgGhhWKNkX6W0LOxvlxlCo1uaSnE8nEsZ4WFqGJbPUy976ktm36JUtm1iuuQWcF88m47pI5HjA+T2yjRlB5Gll7naTO3S+sdF4SJoajUckiy4pyT2w64m1XkWs9ar1O50YoaWzxjPhMyjq5p6m6cWHki85AI5VR7zpcp8RcUXSM2ECwcvgqdGO28229VOIVhEA15Oxqw21qy2v+joWN4aDZ1h4RsXCoulC8UabzC/gnHkpGZuxaWROOOzmzMhiTDQ3pRkZiZYmVptuQF3Za+pbFvbqBcLzz2KvzmKvRrFXL2CvRLFXo9RdnMFevYCVyOLeU3bKiPAbFQ3tGeW7yhSMqbMOD81xBKCSLkcsK9knr5RbLgglRNl2MkZqK2itTVepa0ykxXK5Xgy9tGyCy5e84fLLRoe8eM4Eli+tbN7d6hGy1xkXfuODxnrXyGhI39fgGSWmK4ZwAXldrDDxlj3EW/awkrnPSRBaOUdo8UyaQRFaPE3j2UdQTuo7aStEvPHyZBeHP2Yi3nhpRdNGi8EgbAkVDqE6ItARyfsY7TjuW4gX0Llk3kgWXSNjYcm8qbjGhnP+Ajz9nKtEcrhDG+qNIvNGmDS6eTIaroTGa5JuT6vZIhSOEQqvEbJXsdUSIRaxE4vYsYuEll+iPvo49lr2JIaJcGuGMZHhtXDXdahxG4aLahOUJwgBQCllYiXq6lCRrQcw63g81U2ytOKOTFtOy8eTZoR4OXrOXzAGiruNk0eMh235hsvXG4MiMy9PRr4e1VCf1fiINV7mJo96dcYXcrBXZn0xEynvRMP5x9JecrSyiTdeSrzJGBSJxktxQk1o281Ma5uYKLPemCwrlZEhBoOwbSjLSnoMKHDGWh1PuAZGhrHhrS8spfJtLK/gzEZheTmpVPyu0A1ND9t9K2naB411hFtt7BYIN2lCDQ6h+jghZw07vkJo+Sz1+mlCzoKbjicdx67P7rWobyfRkDIskrnvBUHICxUKmfip1sLiuvxorWEtltQRemmTYfPJ5Qp6eRUWl3DOXUirnxd5Gx97oOFqY3xcWoddFyOs5glzkXBilnDsPKHlF2mJPp7Mc5MLxwqjQ404doNZJo0Ln2Hh7t8MMRiEQKBCNrS1otq2NoxJx+IpF2fGm8dGy/jyCrHTK6wspxQGMe+NXgHNQBNWnSbUkCDU6JhlkybUogg1rxBqfIFQ/TM0hGPY9npvgEaRsFpI2M1oqw5thc3oCbsObde7QzHNaAoddlNTW+7+dUtzrOOWCYKQHaUU1NeZIM0ivB3gek5X11LGhZd92Mvhk1y6uXr8RsjicrrxsbqWO28PuMbHHuzmMFaDjapX2PVg1SusOjP9jRXWqBBYIQfLdrDshPlYCSxrAUtFCRNHsYZFDEuvuVkyNkYMBmFHYALACh9tkknK8Egl7DJejGViyyvE/IphaRkupuqo1SVsZ5EQS6ZLpMEh1Jgg1JjArouibI1l6+TSCmmU5SuT1ioIVYeyrFS8RJFox4G1WB6JAl0jZHUVJ2aG08ZjMViKw8WY0VOxmDvUVqPX3KG7sZjPIPHmPfZGk2nX2bnxnHJlV0FKqQFMKtkIMKO1nt5KHUGoBkpleMRjceIrK+kKYN0Y/lj6GP21GMRXUbFVVGINlVhxl2sm6ZWzhqVjKCeG0nHgkdJ86TJSjH4QvSHUGsqyUvkvKH2iOpMvxEkaE0nDIhZzjYoY/Ox/2PD4shoMSqkuoF9rPexujwGDhdYRhFojaXgU0Q/rx8lYAvBX1T2ZZDH6QfSGIBSO8rKIuvMXFRqWXe5oqz7MG0ASpVTPFuoIglB7FKMfRG8IwjZTboMhApz3bc8CXVuog1JqSCk1qZSanI1eLLmggiBsO8Xoh7z0hiAIpSMwYVRa66PAUTCTT1VYHEHYIsrkZ/B/krtUxtKdY08p31SDvjkIvTSu/vqZ59jhKKWGgCGAvbsLHMsrCEIa5TYYvIAkjw5gZgt1BCEHOR7Eyl8no767uel5LAutvEle3CXKne9BJed10P79lmXmfUChvMQpSvkMASu5nbzuzqIY/dCRx7HyoiFUAJ/Rvu4lINsLQuYxvvOsOz7XdazkOVIvG55+8V40lO/lwn8cOXPClNtgOAGM+LYjWSKZ86kjVC0qNR1w8uOWZ3s4qw2Og/X1k2XZHu4qlQYZUg9z76Gt3Dfx5HEqVe497P0P6jI/yOVJlZUt6wel1EwexwqVxGtHVuohlvZgSivLaGO5tjetu65CfufK9nBPe8BnrPj2b2XCrDR95OktrOT9WK/TUvX9+9LkKTNlNRi01lGl1JhSqs8tSjZwpdQUcHCzOjuTLA/SDatu1lCyNBL/w9jKtEgzHupJa9RKL894UJcaMwuj10jMZEoo703e9r3le+t2doUjVD3F6Ict6Q2lzGyhhQtaWFmuB9JWyvIh18Nt0/Jsb7BuQdpDP4t+8HeNKWWGBtYoGz3kdaEP+YBmdy17DIPWemKD8t5cdYrH/aF7DVlrtv6ul/lQ3egNOqNRZbWw/Y0wJWelp3UtHSrpmtfJB7qdmqbZMwb8hoFM1bwjKUY/FKo3QuEwnVfuzXYes/SXZdZJq79xPdBobcqzfrT2rSePKORrFIRyLQDv/cBSKhkO4/rgzLpK1bZcVWV5b89kfv/0+7VuqdO309b0BsdkWabOp9f9b9L/B/ndv9TD3p06Hk9H+fXS1h/yyr1/aY5UvPubft9J1vOVq83PpbKca9133Oi7ZyvLYvTmupOBCXpMw7ahbVfqTTlt6l9AWXlbudpvTGj/zzX18K+dB3kpMG/72g6DCqHtENoKmzd+dykPfqEasZSipUpn8dSsfyhmrbfBDn9x8oFf8Cj7AJB8IfO9ZCjb9Me7LyVaWWhco8B96GuvuxG/AZeyXlyNn3K8ZDzIUw/t9Ae59wDfKQTSYLBsm9ZdbZDtn0fqH+htZ7NWvY+jSe71vxl4xyfPk8Xdv5EFCOk/omR5lh8jGef29nlS6cw3Fp/8AI5vO/O7OJD2Xf3fad33S8pioe0wyg6hrRBYYZQVcpOTh1CuS1dlHu+755tR1P4cB6ucZy+BDP56eVbcOepE2Coq319vrf2YlALbnfzA80Z6QcMZRkFyX7bTZCyF8hBIgyGkFB3hcKXF2CZUaVuBZUPSAHA9Aq4xYLbFOyAIQhmwwxCqdw2EOnd9p+jx2iCQBoOQB5btNsh6CNWZpR3O/7VYEARhKyR1j98wqBPdUwOIwRB0PJee1yg9I0E8BYIglJOsuqfOGAxCTSIGQ5CwQ67HwG+911VaKkEQahmljHcy02Mg3Qk7DjEYqhHxGgiCUCksC0KNEG4wy1C9dCcIgBgM1YFlQajBfMKNZikNVBCE7cCyjd4JN7oGgngtheyIwVAJ7FCGBS8NVBCEbcIOp15Mwo3StSDkjRgM20GoLt1A2Ep6WkEQhK3g1z/hJglKFLaMPLlKjVK+7gXXQJDYA0EQtgMv/incJPpHKDliMJSKcAPUt0JdqzRQQRDKRzLZmptwzQ6lEq9JrhWhjIjBUAx22BgJ9W3SzVAhlFJZJ1EpBUeOHAHgrrvuKsv5BWEdSrkP/1Q69nXrYhAUjeiNrSGvwoVi2dAYgchV0L4PmjrEWKhyTp48STQarbQYgpCdcCO0Xg7tL4fObqNXdu2FlsuMfmloSwUnirGwbYjeWI8YDPmglPEktO2Bjquh+RIzNjlAPP7443zuc5/jxRdfrLQo287x48eZmZmptBiCkMKyobE9ZRzUt8qLR5UhemM9YjBshFJQ12Qs/44us6xrqrRUW+b73/8+Dz/8MH/+53/Oxz72MR599FHi8XjO49rb2zly5Ai9vb309vYyPT3N4OAg3d3dSdcbGGvcq3P48OG047190WiUkydP0t3dTXd3NydPntz0uoUeNzg4mJRhYmICgOHhYSYmJrjzzjvp7e1Nkzfb+Y4cOUJ7ezv9/f2Mj4/nvD+CUBB1TdC6233x6Kz6IY3f+MY3+KM/+iMeeugh1tbW8j5O9EaNorUO3Kf31TdpffaJ8nxmn9J6aVbrRFxXq3CTjQAAC/BJREFUCsdxSv75zGc+oz/ykY8kP/fcc4++55579Be+8AV9/vz5DWUB9Pj4uNZa64GBAd3T06O11vrChQs6EolorbU+depUslxrre+66y49MjKSPH5gYEBrrfXU1FRavZ6eHn3hwoUNr5vPceYnrJMyZcqmtdZ9fX16amoqub3R+aampnRXV1fyPENDQ8nvEUSASV0F7bVaPmXVG5vqlCe1XjyvdTxW3D+0Anz+85/XH/nIR/RHP/pR/dGPflTfd999m+oLD9Ebtak3xAcGrjehBRp2mdEOFWRubo4//MM/LPt1vLeFb37zm0xPT/PBD35ww7p9fX0A3HrrrXR1dQEQiUSS/Xujo6McOnQoWf/uu++mt7c3GfRz7NgxACYmJohGo/T39yfrzszM0NPTk/W6hR43MzPDvffeSzQaTX4ikci68250vomJCQYGBpLHZDtWEPJCKTO0saEN6pq35ZIvvvgisVispOf02rh33unpaR555BH27t3LG97wBrq7uzc8VvRG7bGzDYZQvdugq2coZFtbGx/4wAdKft777ruPb33rW8ntcNi4Qvfv389rX/vakl5rdnY2bdvfgIaHh/OOHi7kuOnpae68807uv/9+IpEI7e3tm5472/k8d6QgbBk7ZEZNNeza1gRJ8/PzfPzjHy/7dRzHwXEcnnnmGT71qU9t+qJRKKI3qp+dZzAo5TbotqoNXPQe5qXEtu3kudva2rj99tu56aabqKsrPi31oUOHGBwcZGhoiEgkwujoKAMDA+vqDQwM0N/fn6w3PT1NV1dXTos8n+NmZmaSZTMzM2nRzZFIJE0ZbXS+vr4+BgcHGRkZAYwy8b9NCMKGWLYJhq5vrcjlW1tb+fCHP1zy8953331MTU0lt+vq6rAsi9e+9rXceuutRZ1b9Ebw2DkGQ7jBGAr1rTtyaNK1115LLBbjtttuY8+ePSU9d09PDyMjIxw8eBCAAwcOMDo6uq5eV1dXWr2urq6k+3Az8jmur6+Pw4cP093dvc7l2N/fz+DgIAcOHGBsbGzD8/X09DAwMJA8Ry27FoUSoRQ0RMyIhyrxUpaDUChEZ2cnb3jDG7j++uuTLyDFIHojeCgT4xAsDty8X09O/F3uipblGgltMsGTsONQSk1prQ9UWo5qIW+9kS/hRmi+tKZ1y+OPP85jjz3G7bffzu7duystjrANbKY3atPDEG50YxNadqQ3QRCCglJqAIgCEWBGaz2dbx2l1AgwCswCd2itj26L0BXufthOrr32Wq699tpKiyFUCbVjMFgW1O8yhkKVj20WBAGUUl1Av9Z62N0eAwYLqNMDjAPTwJ3bILAJZmzsqOnuB0HYiEAaDFOPPragLn3FDyotR4HsAi5WWIZLgHNlvkY5vmcp5C5Grq1cv9DrleN/s6/E5ys1fRjPQRKlVE+Gl2GzOqNa642z+Ji6Q8CQu7mqLn3FY8WLva2I3tg6oje2xoZ6I5AGA/CDoPXNKqWOaq2HctcsqwyT5b5v5fiepZC7GLm2cv1Cr7cd/5sqJAKc923PAl0Yj0E+dbqUUn1uHbIZD243xVEI5j0WvVHUOUVvlJigGgxB5POVFmCbqNbvud1yVet9qBm01skcw0qpKaXUhNa61mYL2im/o2r9nqI3fIjBsE1orav6h1AqqvV7brdc1XofthO3O2CjMWZHSQUyenQAmbP9ZK3jBkJ2+YyGbN6JwLNTfkfV+j1Fb6QTVINhe6Kha4+g3rdKy70d16/0dyw5uUYtKKVOACO+okiWURJZ6ygz+slvXHRkG2GRQc3d420iqPet0nLXnN4IZB4GQRBqAzcGIYnWesItnwIOaq2jm9Tx0gJ2AdNeuSAI5UEMBkEQBEEQchLULglhC3gJcORNrPIopXqAA363fT5JjARhuxG9UT1UWm8EymDwuSBvBR7KNQZbSKGUigCHgOOVliVfXJnvwPRVd21bJr/U9b0he1FIucJLdN5h4CFfWc4kRsLWEL2xdURvbOn6Nas3ApOuzLWsolrrk1rrw8Ax94ch5McBfD+0gDCmtT7qNrhtnf7N/W11ub+3CUxWwZLgnm88ozhrgqJSXXOnInqjaERvFECt643AGAyYwCb/P98bRiXkQCnVFzR3omtNJxuC1npb37bd8fzDrgUP0FnmS26UoEgoDtEbW0T0RuHUut4ITJeE60Y8CUkrLp9hVDse94ebObY9CHRBmnsvaya/MnMYmHKzqdXmBPc1juiNrSF6oyhqVm8EycPgZwTorbQQAaEH6HH7cW8FbvVZv9WM19gn3AY/XAG5ezCTGnW4w/zKSWaGwmxJjITiEL2RP6I3tk7N6o3AGQzuD3hUay3KNA/cvrSTbuOZwQR9BeHezZD+w49Swv7AXLi/s2n33vUCk5n5AErMCaDbt50tiZGwRURvFIboja1R63ojUAaDe+On3UxvkYBYvFWBGwjTBxwKyH2bIL0vbrvT/mZa6mOUyHJ3f8f9QL+nTNy+zzGlVJ9bNrLZOYT8Eb2xdURvFExN643AJG5yf7hjpFwwXVrr9gqKJJQZ11rvcDdnt7sv0p0LAUwgkYxDDyCiN3YeojfKR2AMBkEQBEEQKkeguiQEQRAEQagMYjAIgiAIgpATMRgEQRAEQciJGAyCIAiCIOREDIZtRuYHEAShUERvCNWAGAzbiDuO+UARx0fc2ciECiP/B2G7EL1ROwT9/yAGQwlxG+boJlUGMJm5tsoYJk95tmvfpZS6aysndeXetvG15bxWMfehQO4NeuMXqgPRG3lfT/RGhRGDoUS4yULuZvOZwjrdzFxbOb83TW/J07O6MnXnrLiNKKUGqnkaYjf9aldAst8JVYrojdIieqO8iMFQItxsYhu+Jbg/kGLmlR8Gjhdx/KZUYZ74Q1T/NMTHMf8XQdgSojdKjuiNMiIGw/Yx4E9RqpQaUUqNK6VOucvxHIFNB8jIie660S4opcYxOcb9+wbcc59y32JQSo156z4Zhtx17SvvU0pN+WSLZDtfJq4sA+6x+R4z5taf8vKju+7ZPuCYf7a3jc6X4z5ccPd714hscn+63O/rffeebGW+00+wjRPbCDsS0RvZjxG9UQm01vIp0Qdj2Y5vsG/Et96HmVUMYCjPc1/wjnG3e4BTvvOMAnf59k356k5hpn3t88uXcbz2fQd/eddG58siowbGNpPBfy133SuLABd85eNAT8b3zfadNrwPPpn63PUxTH/wRuca8d1D73zryjL+36cq/buTT7A/ojdEbwTlE0IoOyrDrajdyUhUYf1YEZ3ej9kHnPSVZe6LuJazR5fWesK1zCOYH+2MXt83OgBMeOVa6xllgoHWnY/ss8DduZkMWY7pUkrdjWl4EaVU5vfc9Dux+X3A/Q7e5C8z7jFdG5xrHDPzWzdGgUxvUOYxS2qSG0EoKaI3Nj1G9EYFEINhexjQWh/JUt4HTOZ5jkKDnkY3uOYJ4A5MsNJGfafZrrXR+dLIaLSbHuO66Y4BB7XWUaXUhRynX3c+tfW55rPKppTqxfxfxpRSw66yXFfmVu/ANH5BKAeiN7IgeqNySAzD9tC5QfmgNlGz+TCb8WYxgbHqPfx9YieBYV+/W49KRQ6PAoNk9I1mHDvgO7YL+NIm59uIzWTwSL6tuNfx74+SboVvdL7N7kNBsimlurTWM1rroxiXYn+2skz587ieIGwF0RuiN6oK8TCUCNdiHQQOKBMQdML9QfdQXJSzhxcoMwNmeI5S6qRS6hTG3ZW00F134GHgfqUU7jF3+o7byC2Y7dhZ93tlPd9GbCZDxnca8X0HP55LbxKjILOeb7P7sAXZelT6GOmDQF+WMo+eLHILQt6I3sh5HtEbVYRygzCEMqGUuisfl1we5+kB7tZaD5ZALKEEKBOJPairb2iZEHBEb9QuQdYb0iVRfjZyKxaE64KMFBjwJJQJVxHPBLHRC4FA9EYNEnS9IQZDGXEb6XjOivkziOkPEyqPvLUJZUH0Rk0TaL0hXRKCIAiCIOREPAyCIAiCIOREDAZBEARBEHIiBoMgCIIgCDkRg0EQBEEQhJyIwSAIgiAIQk7EYBAEQRAEISf/Dx0OevGMm411AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 540x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f = plt.figure(figsize=[7.5, 3])\n",
    "\n",
    "ax = plt.subplot(1, 2, 1)\n",
    "plt.title(\"clean\")\n",
    "print(var_[0])\n",
    "errorfill(1.0 / taus, np.array(mse[0]), \n",
    "          np.array([np.sqrt(a) for a in var_[0]]), label='local', color='#17becf')\n",
    "errorfill(1.0/taus, np.array(mse[1]), \n",
    "          np.array([np.sqrt(a) for a in var_[1]]), label=r\"Ditto ($\\lambda^*$)\", color='#d62728')\n",
    "errorfill(1.0/taus, np.array(mse[2]), \n",
    "          np.array([np.sqrt(a) for a in var_[2]]), label='global', color='#ff7f0e')\n",
    "\n",
    "plt.ylabel('test error')\n",
    "plt.xlabel(r'$1/\\tau$ (device relatedness)')\n",
    "plt.xscale('log')\n",
    "\n",
    "\n",
    "plt.legend(frameon=False, fontsize=12)\n",
    "plt.ylim(-0.08, 0.28)\n",
    "plt.xlim(2, 12)\n",
    "ax.xaxis.set_major_formatter(mticker.ScalarFormatter())\n",
    "ax.xaxis.set_minor_formatter(mticker.NullFormatter())\n",
    "plt.xticks(np.array([2, 4, 6, 8, 10]))\n",
    "\n",
    "\n",
    "ax.arrow(3.3, -0.062, 0.7, 0, width=0.00013, head_width=0.015, head_length=0.25, fc='gray', ec='gray')\n",
    "ax.annotate('more related', (4.5, -0.07), fontsize=13)\n",
    "\n",
    "ax = plt.subplot(1, 2, 2)\n",
    "plt.title(\"adversarial\")\n",
    "errorfill(1.0 / taus, np.array(mse2[0]), \n",
    "          np.array([np.sqrt(a) for a in var_2[0]]), label='local', color='#17becf')\n",
    "errorfill(1.0 / taus, np.array(mse2[1]), \n",
    "          np.array([np.sqrt(a) for a in var_2[1]]), label=r\"Ditto ($\\lambda^*$)\", color='#d62728')\n",
    "errorfill(1.0 / taus, np.array(mse2[2]), \n",
    "          np.array([np.sqrt(a) for a in var_2[2]]), label='global', color='#ff7f0e')\n",
    "plt.ylabel('test error')\n",
    "plt.xlabel(r'$1 / \\tau$ (device relatedness)')\n",
    "plt.xscale('log')\n",
    "\n",
    "plt.legend(frameon=False, fontsize=12)\n",
    "plt.ylim(-0.05, 0.18)\n",
    "plt.xlim(2, 12)\n",
    "ax.xaxis.set_major_formatter(mticker.ScalarFormatter())\n",
    "ax.xaxis.set_minor_formatter(mticker.NullFormatter())\n",
    "plt.xticks(np.array([2, 4, 6, 8, 10]))\n",
    "ax.arrow(3.3, -0.04, 0.7, 0, width=0.0001, head_width=0.01, head_length=0.25, fc='gray', ec='gray')\n",
    "ax.annotate('more related', (4.5, -0.044), fontsize=13)\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "#f.savefig(\"pe_fig2.pdf\")"
   ]
  },
  {
   "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"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "metadata": {
     "collapsed": false
    },
    "source": []
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}