experiments/test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2975f7f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "\n",
    "from mcpdf.nnpdf import data, defaults, theory\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ea482711",
   "metadata": {},
   "outputs": [],
   "source": [
    "# datasets\n",
    "datasets_list = [\n",
    "    {'dataset': 'HERACOMBNCEM', 'frac': 0.75}, \n",
    "    {'dataset': 'HERACOMBNCEP460', 'frac': 0.75}, \n",
    "    {'dataset': 'HERACOMBNCEP575', 'frac': 0.75}, \n",
    "    {'dataset': 'HERACOMBNCEP820', 'frac': 0.75}, \n",
    "    {'dataset': 'HERACOMBNCEP920', 'frac': 0.75}, \n",
    "    {'dataset': 'HERACOMBCCEM', 'frac': 0.75}, \n",
    "    {'dataset': 'HERACOMBCCEP', 'frac': 0.75},\n",
    "    {'dataset': 'CHORUSNUPb_dw_ite', 'frac': 0.75}, \n",
    "    {'dataset': 'CHORUSNBPb_dw_ite', 'frac': 0.75}, \n",
    "    {'dataset': 'NTVNUDMNFe_dw_ite', 'frac': 0.75}, \n",
    "    {'dataset': 'NTVNBDMNFe_dw_ite', 'frac': 0.75}]\n",
    "\n",
    "ndatasets = len(datasets_list)\n",
    "\n",
    "# load data and covmat for all the dataset\n",
    "y = data.values(fit=defaults.BASELINE_PDF, dataset_inputs=datasets_list)\n",
    "cov = data.covmat(fit=defaults.BASELINE_PDF, dataset_inputs=datasets_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "45435f05",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2053,)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "f21821c7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2053, 2053)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cov.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "69745980",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "-- Reading COMMONDATA for Dataset: CHORUSNUPb_dw_ite\n",
      "nData: 607 nSys: 1014\n",
      "-- COMMONDATA Files for CHORUSNUPb_dw_ite successfully read.\n",
      "\n",
      "\n",
      "-- Reading COMMONDATA for Dataset: CHORUSNBPb_dw_ite\n",
      "nData: 607 nSys: 114\n",
      "-- COMMONDATA Files for CHORUSNBPb_dw_ite successfully read.\n",
      "\n",
      "\n",
      "-- Reading COMMONDATA for Dataset: NTVNUDMNFe_dw_ite\n",
      "nData: 45 nSys: 1003\n",
      "-- COMMONDATA Files for NTVNUDMNFe_dw_ite successfully read.\n",
      "\n",
      "\n",
      "-- Reading COMMONDATA for Dataset: NTVNBDMNFe_dw_ite\n",
      "nData: 45 nSys: 103\n",
      "-- COMMONDATA Files for NTVNBDMNFe_dw_ite successfully read.\n",
      "\n",
      "\n",
      "-- Reading COMMONDATA for Dataset: HERACOMBNCEM\n",
      "nData: 159 nSys: 170\n",
      "-- COMMONDATA Files for HERACOMBNCEM successfully read.\n",
      "\n",
      "\n",
      "-- Reading COMMONDATA for Dataset: HERACOMBNCEP460\n",
      "nData: 209 nSys: 170\n",
      "-- COMMONDATA Files for HERACOMBNCEP460 successfully read.\n",
      "\n",
      "\n",
      "-- Reading COMMONDATA for Dataset: HERACOMBNCEP575\n",
      "nData: 260 nSys: 170\n",
      "-- COMMONDATA Files for HERACOMBNCEP575 successfully read.\n",
      "\n",
      "\n",
      "-- Reading COMMONDATA for Dataset: HERACOMBNCEP820\n",
      "nData: 112 nSys: 170\n",
      "-- COMMONDATA Files for HERACOMBNCEP820 successfully read.\n",
      "\n",
      "\n",
      "-- Reading COMMONDATA for Dataset: HERACOMBNCEP920\n",
      "nData: 485 nSys: 170\n",
      "-- COMMONDATA Files for HERACOMBNCEP920 successfully read.\n",
      "\n",
      "\n",
      "-- Reading COMMONDATA for Dataset: HERACOMBCCEM\n",
      "nData: 42 nSys: 170\n",
      "-- COMMONDATA Files for HERACOMBCCEM successfully read.\n",
      "\n",
      "\n",
      "-- Reading COMMONDATA for Dataset: HERACOMBCCEP\n",
      "nData: 39 nSys: 170\n",
      "-- COMMONDATA Files for HERACOMBCCEP successfully read.\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(2053, 9, 50)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# load FK tables\n",
    "fks = theory.theory(dataset_inputs=[ds['dataset'] for ds in datasets_list])\n",
    "fk_bare = np.concatenate([fk.elements[0].table for fk in fks], axis=0)\n",
    "fk_bare.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ca33cde9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2053, 450)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# reshape FK table as a matrix\n",
    "ndata, nbasis, nx = fk_bare.shape\n",
    "fk = fk_bare.reshape(ndata,nbasis*nx)\n",
    "fk.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "be179b38-52fc-4986-ac0b-f7a6eb10574c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def distinfo(values):\n",
    "    print(f\"   max = {max(values)}\",\n",
    "          f\"   min = {min(values)}\",\n",
    "          f\"absmin = {min(abs(values))}\",\n",
    "          sep = \"\\n\")\n",
    "\n",
    "    sns.boxenplot(x=values, scale=\"exponential\")\n",
    "    plt.xscale(\"log\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "9c05135a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   max = 308.39648487355925\n",
      "   min = 1.8638700388327588e-06\n",
      "absmin = 1.8638700388327588e-06\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAD8CAYAAABJsn7AAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAANvklEQVR4nO3db2hV9x3H8c8vN2r8NzWapo21dcPWbQzHaDYKK6M0GYyxrFjm2OijztE/MHMfFWJpfWBLFYRBdQ9qaVawuI2lGFZZYXAbZGvZAxOY/ROtXKRWjU1jop1OY3KT3x7EpPlzvfeem3PO996b9+tZ7j055+uv3re3J+fmOO+9AADxq7IeAAAWKgIMAEYIMAAYIcAAYIQAA4ARAgwARqqDbLxu3Tq/cePGiEYBgMrU09NzyXtfN/vxQAHeuHGjuru7w5sKABYA59zZbI9zCgIAjBBgADBCgAHACAEGACMEGACMEGAAMEKAAcAIAQYAIwQYAIwQYAAwQoABwAgBBgAjBBgAjBBgADBCgAHACAEGACMEGACMEGAAMBLolkSI39GjR5VKpazHCKS5uVktLS3WYwAljwCXuFQqpf98dFJjy2pD3/fw1cuSpJqVa0LbZ+L6kCQRYKAABLgMjC2r1Y1v/jT0/Y5e6JUk+fXfDm2fS0+9E9q+gErHOWAAMEKAAcAIAQYAIwQYAIwQYAAwQoABwAgBBgAjBBgAjBBgADBCgAHACAEGACMEGACMEGAAMEKAAcAIAQYAIwQYAIwQYAAwQoABwAgBBgAjBBgAjBBgADBCgAHACAEGACMEGACMEGAAMEKAAcAIAQYAIwQYAIwQYAAwQoABwAgBBgAjBBgAjBBgADBCgAHACAEGACMEGACMEGAAMEKAAcAIAQYAIwQYAIwQYAAwUm09QKU6evSoUqnUvPeTTqelqpUhTBSfdDqtZDJZ1Pc2NzerpaUl5ImA0kSAI5JKpXTi45OqWlk3r/2M3xyVloY0VEz+d3NUH352KfD3jV8dkCQCjAWDAEeoamWdln3/sXnt4/rxI1JmLKSJ4lHsn/v68SMRTAOULs4BA4ARAgwARggwABghwABghAADgBECDABGCDAAGCHAAGCEAAOAEQIMAEYIMAAYIcAAYIQAA4ARAgwARggwABghwABghAADgBECDABGCDAAGOGecAEVerfjdDotLV4d/UAVZj53VJ7EnZVRLghwQKlUSh/2ntKS1fU5t7s5kpEWxzRUBbkxktHpvstFf//NK/2SuLMyygMBLsKS1fW6+5HHc25zvuuwbpbZ3YxLQSFrm8v5rsMhTgNEi3PAAGCEAAOAEQIMAEYIMAAYIcAAYIQAA4ARAgwARrgOOItcn3ZLp9PSsrUxT4Qgwvg0XSniE36VhwBnkUql9NHJT7Ri7Z1znhsezfABtxI3PJrRp198aT1GqK4Nfi6JT/hVGgJ8GyvW3qnvPvrbOY+f+NvrGsmMG0yEQq1Ye1fW/3bl7MTfXrceARHgHDAAGCHAAGCEAAOAkQV7DjjflQ7VX6uLeSIgN+urO7gKI3yxBXhgYEB1dXUaGBiQJNXV1c14fPZ2UUulUuo99YlW16+f89xIZmzh/suEkjWSGVPf5Wsmx77Sf0FSeFdhTH+dF9uAuFoR5bFi6cwHH3yg1tZWJZNJvfLKK5Kk/fv3S5JaW1u1f/9+bdmyZWq7ya+jtrp+vR55fMecx7sOH9DIKL/LF6Xldn9f49B1+EBo+5r+OpeKa0CcrYjyWJEHOJPJaM+ePZKkAwe++o/48ssvyzknSdq7d6/eeOONqe327t2rQ4cOqbq6+PHy3ToonU5rWZbrfAFkV+gpkFynKqb3YM+ePfLeSwrWgOn7CKMVuUR9rMgD3NnZqcuXJ24xMz7+1fWz/f39SiQSkqShoSHt3r17aruhoSG99NJLU18XI51Oa3RsXHUNG7I+PzrGtbxAEKNj47p09UbObQb6zimdTs948zM9yNN7MHk6Upp4zb/44oszGtDZ2alt27bNOcb0feTaLgxRHyvSAA8ODqq9vV3Dw8NznvPeK5PJSJKGh4f1/vvvTz03PDysY8eOqWbpMjVsuLeoY4+Ne9U1bNC2p57N+nzHwX0a4ZZBQMFyvZ4mdRzcp8GL53X1xogkqe/c2akgZzIZ9fb2Tr3rnXz9SxOv+ffee2/G1+3t7WpqalJtbe3U47ObcrvtwhDHsfIG2Dn3pKQnJemee+4JtPOurq4Z73qDathwr556dldR33tw327uyQYYmP66Pbhvt86cPqkTJ04E3s/4+LjefffdGe84szUl23ZhiONYeQPsvX9N0muS1NjY6IPsvKmpSe3t7UWONvGv58F9u4v+3rV33V30sQEUZ/rrtu/cWS1fvlybNm2a8w44n6qqKjU1Nc14LFtTsm0XhjiOFekpiNraWm3fvj3raQjnnBKJhDKZjGpqavTAAw+op6dHw8PDqqmp0YMPPjivc8CJKqeBvnPqOLgv6/MDfee06o6GovcPLDS5Xk/Tt1mUqNLKpRO/smrz/ffNOAfc0dEx1YPJH2ZNNqCxsVHd3d1TDdi+ffuc/9Wf3ZTbbReGOI4V+Q/htm7dqs7OTvX19amqqmrqLX19fb2cc7p48aJqa2u1a9cuPfHEE+rr61Ntba2ef/75SK+CuJLgQ4BAEIsSVVq3cmnObdZtvj/nVRDTezB5Xe1kA1544YUZDdi6dWvefeTaLgxRHyvyAFdXV6utrU2tra3asWPH1HXAzz33nKSJawDb2tq0ZMmSqe3a2trmfalHS0tLzovGk8mk2UXtQDnatGnT1Ou3WNN7sHPnTknBGzB9H2G0otB5ozhWLB/E2LJlizo6OlRXV6eHHnpI0lefhJt8fPZ2cbjSfyHrBeZX+i9oWS3XCKO03O7va1zHblizOZR9zX6dF9OAOFsR5bFi+8Tt5PCz/xD5vo5Kc3PzbZ+7PpiIZQYgiMXVCTWsWWFy7IY1m3O+ZoKa/jovtgFxtSLKYy3YX3mQ6xRFMpmsuDsqoPyFcQoApYWfRAGAEQIMAEYW7CmIfK4Nfp71PlzXBi9q8ap6g4lQqGuDFyvuHmrXBj+X7lhlPQZCRoCzyPXDhvR/B8Sv8SltNYuqtbHSYnXHqlB/CIbSQICzyPcDutN9xX9CD9Hjh1UoF5wDBgAjBBgAjBBgADBCgAHACAEGACMEGACMEGAAMEKAAcAIH8Qows0r/TrfdTjvNlqxLqaJKkcha5vv+9WwJsSJgOgQ4IAK/Tho+vqgbkQ8SyVaurham+YT0IY1fGQXZYMAB5TvVkeTksmkPvzsUgwTVRY+RoyFhHPAAGCEAAOAEQIMAEYIMAAYIcAAYIQAA4ARAgwARggwABghwABghAADgBECDABGCDAAGCHAAGCEAAOAEQIMAEYIMAAYIcAAYIQAA4ARAgwARrgnXITGrw7o+vEj896HltaGNFE8iv1zj18dkMSdpLFwEOCIhHVn3nT6ir4MZU/xWb5kkTbdU0xI13FHYywoBDgihd49OZ9kMqmeM/0hTBQf7mwMFIZzwABghAADgBECDABGCDAAGCHAAGCEAAOAEQIMAEYIMAAYIcAAYIQAA4ARAgwARggwABghwABghAADgBECDABGCDAAGCHAAGCEAAOAEQIMAEYIMAAYIcAAYIQAA4ARAgwARggwABghwABghAADgBECDABGCDAAGCHAAGCEAAOAEQIMAEYIMAAYIcAAYIQAA4ARAgwARggwABghwABghAADgBECDABGCDAAGKm2HgD5Ja4Paempd0Lfr7t6WZJUc/XT0PaZuD4kqT60/QGVjACXuObm5gj3HkUo6yOeGagczntf8MaNjY2+u7s7wnEAoPI453q8942zH+ccMAAYIcAAYIQAA4ARAgwARggwABghwABghAADgBECDABGCDAAGCHAAGCEAAOAEQIMAEYIMAAYIcAAYIQAA4ARAgwARggwABghwABghAADgBECDABGAt2U0zk3IOlsdOOEZp2kS9ZDlAHWqTCsU2FYp9u713tfN/vBQAEuF8657mx3IMVMrFNhWKfCsE7BcQoCAIwQYAAwUqkBfs16gDLBOhWGdSoM6xRQRZ4DBoByUKnvgAGg5BFgADBCgAHAyIILsHPuYefcv5xzrzrnHraep1Q55751a43ecs49Yz1PqXLOfcM51+6ce8t6llLD2uRXVgF2zv3ROfeFc+6jWY//xDn3iXMu7Zxry7MbL+mapBpJ56Oa1VIY6+S9P+m9f1rSLyX9MMp5rYS0Tme899ujnbR0BFmzhbY2xSirqyCccz/SRDwPee+/c+uxhKTTkn6siaAel/RrSQlJe2bt4jeSLnnvx51z9ZJ+771/PK754xLGOnnvv3DO/VzSM5Le9N7/Ka754xLWOt36vre897+Ia3YrQdbMe9976/kFsTbFqLYeIAjv/T+dcxtnPfwDSWnv/RlJcs79RdKj3vs9kn6WY3eXJS2JZFBjYa2T9/5tSW875/4uqeICHPLfpwUhyJpJ6o15vLJTVqcgbmO9pHPTvj5/67GsnHOPOecOSnpT0h8inq2UBF2nh51z+2+t1TtRD1dCgq7TWufcq5K+55zbGfVwJSrrmrE2+ZXVO+AweO+PSDpiPUep894fk3TMeIyS570flPS09RyliLXJrxLeAV+QtGHa13ffegwzsU6FYZ2CY82KVAkBPi7pPufc151ziyX9StLbxjOVItapMKxTcKxZkcoqwM65P0v6t6TNzrnzzrnt3vuMpN9J+oekk5L+6r3/2HJOa6xTYVin4FizcJXVZWgAUEnK6h0wAFQSAgwARggwABghwABghAADgBECDABGCDAAGCHAAGCEAAOAkf8DXvDvTOMbY0MAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# diagonalize experimental covmat and check its eigenvalues\n",
    "l, u = np.linalg.eig(cov)\n",
    "distinfo(l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "be8a6443",
   "metadata": {},
   "outputs": [],
   "source": [
    "#inv_cov = u @ np.diag(1./l) @ u.T \n",
    "#inv_cov@cov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "1c5e7247",
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute inverse of the posterior covmat according to Eq.46\n",
    "v = u.T@fk\n",
    "posterior_cov_inverse = v.T @ np.diag(1./l) @ v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "888bb234",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(450, 450)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "posterior_cov_inverse.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "fbfa1e3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# add a gaussian prior proportional to regularize\n",
    "_, N = posterior_cov_inverse.shape\n",
    "posterior_cov_inverse_prior = posterior_cov_inverse + 0.001*np.identity(N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "a488fb68",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   max = 290362387.0521207\n",
      "   min = -5.200987252054414e-09\n",
      "absmin = 6.259910166570457e-14\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAD8CAYAAABJsn7AAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAMHklEQVR4nO3dXYxcdR3G8ee3u/QFLC/b3fJSKAtpKTSRCGlMvFBJWI0hrLqIUeIFImiaCIx3WLhovLFy4QWlxAr1pTFSBWKJJejFoEESuXBr0NSywFhK6SLu4kJZKLT78veiQ9nubsvO9pzzzOx8P8kmnTOdc34Mp9/MnpkzJ1JKAgAUr8U9AAA0KwIMACYEGABMCDAAmBBgADAhwABg0lbLX+7o6EhdXV05jQIA89OuXbveSCl1Tl1eU4C7urrU19eX3VQA0AQi4pWZlnMIAgBMCDAAmBBgADAhwABgQoABwIQAA4AJAQYAEwIMACYEGABMCDAAmBBgADAhwABgQoABwIQAA4AJAQYAEwIMACYEGABMCDAAmNR0SSKgaDt37lS5XHaP0RS6u7vV09PjHqOpEGDUtXK5rOd2P6/x09vdo5yy90felCQtWnKOeZLpWg8NSxIBLhgBRt0bP71d711+nXuMUzY6sEeSlJavMU8y3eL+J90jNCWOAQOACQEGABMCDAAmBBgATAgwAJgQYAAwIcAAYEKAAcCEAAOACQEGABMCDAAmBBgATAgwAJgQYAAwIcAAYEKAAcCEAAOACQEGABMCDAAmBBgATAgwAJgQYAAwIcAAYEKAAcCEAAOACQEGABMCDAAmBBgATAgwAJgQYAAwIcAAYEKAAcCEAAOACQEGABMCDAAmBBgATAgwAJgQYAAwIcAAYEKAAcCEAAOASZt7AORv586dKpfL7jHmpFKpSC1L3GM0hUqlolKplOk6u7u71dPTk+k65xMC3ATK5bL+uadfC85a5h6lZkcOj0mL3VM0h0OHx9Q/MJzZ+o4cHJQkAnwSBLhJLDhrmc777E3uMWr2+tPbNXp41D1GU1hwdrb7yOtPb89sXfMVx4ABwIQAA4AJAQYAEwIMACYEGABMCDAAmBBgADAhwABgQoABwIQAA4AJAQYAEwIMACYEGABMCDAAmBBgADAhwABgQoABwIQAA4AJAQYAE64JZ1TU1YorlYq0uD337QBT5XGl5Q/MhysuE2Cjcrmsfz3/gs7svCDX7RweHVcrVxaGweHRcb36xkjm63176DVJjX/FZQJsdmbnBfrUV9fluo1nH92id7myMAzy2r+ffXRL5ut04BgwAJgQYAAwIcAAYEKAAcCEAAOACQEGABMCDAAmfA54FvI6Y61SqWjh2edmvl6gGcz1LLt6OoOOAM9CuVxW/wsvatnyFZmud2x8QgszXSPQPMbGJzT8zvs1PWZwYL+k+jmDjgDP0rLlK3TTd+/KdJ3bH7hX777PGWrAXMzl3+T2B+7NaZq54RgwAJgQYAAwIcAAYFLIMeCivvc2L5VKRe3nXegeA0AGZvr0hOuTEYUEuFwu66WXKlrRdUkRm8vcREruEQBkZCIlHR4dP3Z7/76XNTo6qlWrVunyyy8/try/v19Lly49druzszPzWQr7FMSKrku0/gc/LGpzmdq44W6+TxeYJy7uuvS4Fm3ccLf27NmtdevWqVQqqbe3Vzt27NB999133OM2bdqkK6+8MtNZOAYMAFX333+/RkZGtHnz5mn3bdy4UWNjY5lujwADQNXExIRuvvlmjY+PT7tvaGhIO3bsyHR7BBgAJhkeHp5x+djYmLZu3XrC++fiIwMcEd+JiL6I6BsaGspswwDQaMbGxvTUU09ltr6PDHBK6cGU0tqU0to83gUEgEbR1tama6+9NrP1cQgCACZpb2+fcXlbW5tuu+22E94/FwQYAKpaWlq0bds2tba2Truvs7NTvb292W4v07UBQAO74447tGTJEt1+++3T7lu/fr3a2rI9daKwEzH273tZGzfcXdTmMvXKvr3qOP8i9xgAMvDKvr3HtWj/vpe1Zs0a3XnnncfOhOvt7dUVV1wxP86E6+7uLmIzuWmJcI8AICMtEVp42oeHGFatWqnu7u7jTkOWNO12HgoJcE9PT918A/1clEqlmr95H0B9Wrly5bTTjF04BgwAJgQYAEy4JtwsDQ7sz/x6UoMD+3XG0vMzXSfQLObyb3JwYL/aV1+W00S1I8CzkNebiMOt/AICzFVba4vaP7aopse0r76srj4UQIBnIa83EUulkl59YyTz9QLNoJ7eTJsrXoIBgAkBBgATAgwAJgQYAEwIMACYEGAAMCHAAGBCgAHAhBMxzN4eek3PProl9220nsn1/FC8vPbvt4dekzpWZ77eohFgo6JOiay81aqxQrYEHG/haa26qGNJ9ivuWF1XpxTPFQE2Kup7kkulkvoHhnPfDjDVfDhdOE8cAwYAEwIMACYEGABMCDAAmBBgADAhwABgQoABwIQAA4AJAQYAEwIMACYEGABMCDAAmBBgADAhwABgQoABwIQAA4AJAQYAEwIMACYEGABMuCZckzhycFCvP73dPUbNjrw1KC0+xz1GUzjyVrb7yJGDg9Ly9szWNx8R4CbQyFePrbw3rIPuIZrE6QvbtDLLYC5vb+h9rwgEuAkUdfXlPJRKJe3a+1/3GE2BKxgXj2PAAGBCgAHAhAADgAkBBgATAgwAJgQYAEwIMACYEGAAMCHAAGBCgAHAhAADgAkBBgATAgwAJgQYAEwIMACYEGAAMCHAAGBCgAHAhAADgAkBBgATAgwAJgQYAEwIMACYEGAAMCHAAGBCgAHAhAADgAkBBgATAgwAJgQYAEwIMACYEGAAMCHAAGBCgAHAhAADgAkBBgATAgwAJgQYAEwIMACYEGAAMGlzDwB8lNZDw1rc/6R7jFMWI29KkhaN7PMOMoPWQ8OSznWP0XQIMOpad3e3e4QM1XPgzp1nz3VjiJTSrP/y2rVrU19fX47jAMD8ExG7Ukprpy7nGDAAmBBgADAhwABgQoABwIQAA4AJAQYAEwIMACYEGABMCDAAmBBgADAhwABgQoABwIQAA4AJAQYAEwIMACYEGABMCDAAmBBgADAhwABgQoABwKSmi3JGxJCkV/IbJzMdkt5wD1EjZi4GMxeDmY93cUqpc+rCmgLcKCKib6YrkNYzZi4GMxeDmWeHQxAAYEKAAcBkvgb4QfcAc8DMxWDmYjDzLMzLY8AA0Ajm6ytgAKh7BBgATAgwAJg0RYAj4tKI+FlEPDZp2RkRsS0iHoqIbzjnO5mI+HREbImIrRHxV/c8s1V9fvsi4nr3LLM1035SjyLiy9X99rcR8Xn3PLMREddExDPVffka9zyzERErIuLxiPh5RHw/j23UfYCr//GDEbF7yvIvRMQLEVH5qCcnpbQ3pXTrlMU3SHospfRtSV/MeOwPZsxi9mdSSuskPSFpWx5zTpntlGeuukvSI/lMOV2O+0khapk/pfR4db9dJ+lrjnmrs9XynCdJ70haJOlA0bNOmq2WmT+uo434lqSrchkopVTXP5I+I+lqSbsnLWuV9G9Jl0paIOkfktZUn7Anpvwsm/S4xyb9eb2kT1T//HADzP6IpCWN8HxL+pykr0v6pqTrG3k/qcf9fNL9P5Z0ddGzzvE5b6nef66kXzfIzEsl/VnSnyTdksc8bapzKaW/RETXlMWflFRJKe2VpIj4jaQvpZQ2Sprtr7wHJF0o6Tnl9JtAVrNHxApJB1NKI3nMOVkWM1d/xTxDR3fi9yLiyZTSRL3P7VTL/BHxvKQfSfpDSunvxU76oRqf8z3V+9+UtLC4KY9Xy8ySRiVtqD7mMUm/yHqeuj8EcQLLJb066faB6rIZRcTSiNgi6aqIWF9d/DtJX4mIn0jamduk09U0e9WtyuF/fg1qmjmldE9K6XuSHpb0UN7xPYks9hOnE81/h6RuSTdGxDrHYCcx48wRcUNE/FTSryRttkx2Yid6nv8o6c7qPrEvjw3X/SvgLKSU/qejx8smL3tX0i2eiWqTUtrgnmEuUkq/dM9Qi5n2k3qUUtokaZN7jlqklH6noy96GkZKabekG/PcRqO+Ah6QdNGk2xdWlzWCRpy9EWeWGnfuDzTi/Mxcg0YN8N8krYqISyJigY6+4fN780yz1YizN+LMUuPO/YFGnJ+Za+F6N7KGdy23S/qPjh4QPyDp1ury6yS9qKPvXt7jnnO+zN6IMzfy3I08PzOf+g9fxgMAJo16CAIAGh4BBgATAgwAJgQYAEwIMACYEGAAMCHAAGBCgAHAhAADgMn/ARp1njgRPPfKAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# check eigenvalues without prior\n",
    "ll, uu = np.linalg.eig(posterior_cov_inverse)\n",
    "distinfo(np.real(ll))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "cdc23081",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   max = 290362387.05312204\n",
      "   min = 0.0009999923879911324\n",
      "absmin = 0.0009999923879911324\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAD8CAYAAABJsn7AAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAALrUlEQVR4nO3dX4yddV7H8c+vnULLLiJDy4rsn0rKYkhsdNMoJqgXjMasO641EkQv3IWqJAq9XLpeIBdavPACkISsRfTGqjTpZmvwZtCwN3thMSshLJrCArsF7MCwgELpTPt40XEsHUrntOec75zT1yuZ5HDO6fN8fz3TN9PnnKdP67ouAAzfmuoBAC5UAgxQRIABiggwQBEBBigiwABFJnp58saNG7vNmzcPaBSA8fTUU0+93nXdptPv7ynAmzdvzsGDB/s3FcAFoLX20ofd7xAEQBEBBigiwABFBBigiAADFBFggCICDFBEgAGKCDBAEQEGKCLAAEUEGKCIAAMUEWCAIgIMUESAAYoIMEARAQYo0tMliVaTAwcOZGZmpnoMLjBTU1OZnp6uHoMxMbIBnpmZybef+U6OXzJZPcrYOvrOm0mS9ZdeXjzJ6rD23bkkEWD6ZmQDnCTHL5nMez/++eoxxtb84WeTJN3V1xdPsjpseO7x6hEYM44BAxQRYIAiAgxQRIABiggwQBEBBigiwABFBBigiAADFBFggCICDFBEgAGKCDBAEQEGKCLAAEUEGKCIAAMUEWCAIgIMUESAAYoIMEARAQYoIsAARQQYoIgAAxQRYIAiAgxQRIABiggwQBEBBigiwABFBBigiAADFBFggCICDFBEgAGKCDBAEQEGKCLAAEUEGKCIAAMUmRjGTg4cOJCZmZm+bvPQoUPJmkv7uk04m0OHDmXnzp3VYywzNTWV6enp6jHo0VACPDMzk6effS4XXXZl37Z57P2FZEPfNgcr8u77C3nu8Fz1GB9w7K0jSSLAI2goAU6Siy67Mj/yC7f2bXuvPbk38+/P9217sBIX/XB/v4/74bUn91aPwDlyDBigiAADFBFggCICDFBEgAGKCDBAEQEGKCLAAEUEGKCIAAMUEWCAIgIMUESAAYoIMEARAQYoIsAARQQYoIgAAxQRYIAiQ7km3BtvvJGkDWNXcEFarVdrTlyx+aMMJcBvvvlmsmFyGLuCC9L788fzvdffqR5jmbdnX0niis1nMrSrIgOD80ObfjQ/e/Md1WMs863HHq4eYVVzDBigiAADFBFggCICDFBEgAGKCDBAEQEGKDKUzwHPz8+nbRjGnoDVZlhn6Y3iGXdDCfDCwkLWDWNHwKqzcPxE5v776ED3ceTwy0lG74w7Z8IBA3Xl1Z/OrX/wlYHuY+9DfzbQ7Q+KY8AARQQYoIgAAxQZyjHgruuGsRvgAraST1ustk9KCDAwFk50Xd6fP37Gx19+8buZn5/PDTfckCTZtGnTsufMzs4uPXb67Q97/vnyKQhgLHxm8zXZde+fnvHx3fd8Nc8++0xuvvnmJMkDDzyQrVu3Lj3+9NNP56677kqS7Ny5M/fff/8Hbp/+/H4Y+DHgV199ddC7AOjZ7t27s7CwkOTkuQq7d+9eeuzBBx9cdvu+++5ben6/DDzAt95666B3AdCz2dnZ7N+/P0myf//+pUMOSXLixIllt+fm5pae3y8DDfBDDz00yM0DnLOFhYXs2bMnzz//fPbs2XPWn26PHj2aRx55JHNzc32b4awBbq39XmvtYGvt4Kn/h1iJffv2nfNgAIO2sLCQRx99dMWHFk6cOJEnnniib/s/a4C7rvta13Xbuq7b1uu7gLfccss5DwYwaBMTE7ntttsyMbGyzyOsWbMmN910U9/2P9BDEHfcsfqu0gqQnIzvjh07cs0112THjh1njfD69etz++23Z3Jysm8zDPxNuL179w56FwA927RpU7Zv354k2b59+wc+57tmzZpltycnJ5ee3y8DD/BVV1016F0A9GzXrl1LP/VOTExk165dS4/deeedy27ffffdKz5UsVJOxADGwksvvpDd93z1jI+//OJ3c/311+fee+9NsvxMuK1bt+axxx5beuzGG2/8wO2RPROutTaM3QAXsDWt5eJ1a8/4+LXXbsnU1NRHhvTUx850u58EGBgLW7ZsWTp9eFT45ygBiggwQBFvwgEDdeTwywO/ZtuRwy9n8rrPDnQfgzCUAPf7oxvA6JhYuyaTH18/0H1MXvfZTE1NDXQfgzCUMq5bty79/UfcgFExim+ODYtjwABFBBigiAADFBFggCICDFBEgAGKCDBAEQEGKOIUNRgDb8++km899nD1GMu8PftKsvG66jFWraEE+PLLL8/s0WHsCS5MF69bm09tvLR6jOU2XjeSpwgPy1ACfMUVV2T28NwwdgUXJKf7jibHgAGKCDBAEQEGKCLAAEUEGKCIAAMUEWCAIgIMUESAAYoIMEARAQYoIsAARQQYoIgAAxQRYIAiAgxQRIABiggwQBEBBigytKsiH3vrSF57cm//tveDI8mGy/u2PViJYz/o7/dxPxx760hy9WT1GJyDoQR4EFdFPfTeXN7q+1bho11y8US2rLbYXT3pysMjaigBnp6ezvT0dF+3uXPnzjz1wn/1dZtwNq4+TD85BgxQRIABiggwQBEBBigiwABFBBigiAADFBFggCICDFBEgAGKCDBAEQEGKCLAAEUEGKCIAAMUEWCAIgIMUESAAYoIMEARAQYoIsAARQQYoIgAAxQRYIAiAgxQRIABiggwQBEBBigiwABFBBigiAADFBFggCICDFBEgAGKCDBAEQEGKCLAAEUEGKCIAAMUEWCAIhPVA5yPte/OZcNzj1ePMbbaO28mSda/82LtIKvE2nfnknyiegzGyMgGeGpqqnqEC4DYfNAnfN/RV63ruhU/edu2bd3BgwcHOA7A+GmtPdV13bbT73cMGKCIAAMUEWCAIgIMUESAAYoIMEARAQYoIsAARQQYoIgAAxQRYIAiAgxQRIABiggwQBEBBigiwABFBBigiAADFBFggCICDFCkp4tyttZmk7w0uHH6YmOS16uHGIBxXJc1jQZrOn+f6bpu0+l39hTgUdBaO/hhVx8ddeO4LmsaDdY0OA5BABQRYIAi4xjgr1UPMCDjuC5rGg3WNCBjdwwYYFSM40/AACNBgAGKCDBAkQsqwK21X2ut/WVr7e9ba79UPc+5aq19rLX2N4tr+e3qefphXF6b0y2+Vgdba1+onqUfWmtrWmt/0lp7sLX2O9Xz9Etr7dOtta+31v6qtXb3sPY7MgFe/I050lp75rT7f7m19h+ttUNn+43ruu7rXdf9bpI7ktwyyHl71eP6fj3JvsW1/OrQh12hXta0ml+bU53D9+FXkvzDcKfsTY9r+mKSTyaZT/L9Yc/aix7X9RM5+WfqtiQ/NbQhu64bia8kP5/kc0meOeW+tUmeT3JNkouS/HuS6xd/M//xtK8rT/l1f57kc9VrOo/17Uryk4vP+dvq2fuxptX82pzH6/SLSX4zyZeSfKF69j6t6e4kv7/4nH3Vs/dxXVck+Zck/5zky8OacWIFjV4Vuq77Zmtt82l3/3SSQ13XvZAkrbW/S/LFrut2J1n2V77WWktyX5J/6rru3wY8ck96WV9O/uTxySTfzir+W0wva2qtfSer9LU5VY+v08eTfCwn/4C/11p7vOu6E8OcdyV6XNP3khxbfM7xoQ15Dnpc13ySexZ/zb4kjw5jxpEJ8BlcnZPfEP/n+0l+5iOef2eSqSSXtda2dF338CCH64Mzre+BJH/RWvuVJAcqBjsPZ1rTqL02p/rQNXVd94dJ0lr7UpLXV2N8P8KZXqf7kzzYWvu5JN+sGOw8nWldDyf549babyV5cVjDjHqAe9J13QM5Ga+R1nXd/yT5cvUc/TQur82H6brur6tn6Jeu695Ncnv1HP3Wdd0zSX5j2PtdtX99XaHDST51yn9/cvG+cTGO67Om0TCOa0pW2bpGPcD/muTa1tqPtdYuysk3PL5RPFM/jeP6rGk0jOOaktW2rup3Knt4R3Nvklfz/x9/uX3x/s8n+c+cfGfzj6rntD5rsqbV8TUK6/KP8QAUGfVDEAAjS4ABiggwQBEBBigiwABFBBigiAADFBFggCICDFDkfwH9TaVAXe4OXAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# check eigenvalues with prior\n",
    "llp, uup = np.linalg.eig(posterior_cov_inverse_prior)\n",
    "llp = np.real(llp)\n",
    "uup = np.real(uup)\n",
    "distinfo(llp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "dff03497",
   "metadata": {},
   "outputs": [],
   "source": [
    "# posterior covmat\n",
    "cov_posterior = uup @ np.diag(1./llp) @ uup.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "3ee0b896",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA/YAAAEeCAYAAAAgkI6SAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAzLUlEQVR4nO3dfdBkV2Hf+e/pft5G86p3yZKMhJFJbK9jbxRUECJwymDZy6I49o6E4gR2oZS4zG5VtpIKLm+ZFNmkWG/VppyCxJ5yFBvHvCh27IgyRsASoiiYjYQNxiIGhJClGYTeZ0bz/nT32T/uOfP0c2f66dfbb/f7mZp65nb38zx37ulz7vmdc+7tEGNEkiRJkiQtpsasd0CSJEmSJI3OYC9JkiRJ0gIz2EuSJEmStMAM9pIkSZIkLTCDvSRJkiRJC8xgL0mSJEnSAlvZ6cnjf/fHIsAdv98BYDU0AQgEACJ+VN605GN+rH0agH3NjVnuTu19Z/MYAK9evwqA451zgHWiCuX2pry9O6wCsEnRTjXT8y+lurLRWJ3ezi64Xse6nY7tahoLbvcoi6P5mIfi1LK7sQbAiVQ/dlkWE9OvXnTS112pfnxXYxcAnznxTQC+Z9dV09vZJdevLJ7bfBmAa1b3AaTadCHPH9XK5fL9zf0APBFPAXC8cxaAVmwDsBZ27Bprgl7R3AvAN1ovAVvH3pwxfZ966pNh1vtQpc3nHx/6zbR6xSsX7pjs2Hr99CeK/0/unH3oL58A4G89cgkA7Sr3TNvkxq3ds0ugaTrQLOpADizt86GyCD6ejCanV2c5O5taotOdTQD2NNbT69jx+3Sh8jHK2/l9vVkK+Hm7kY7x1St7t33/qVQ/VoKLwyatX71YLbVFX28dBWBXc21Ke1gf/cri2rUiSL7cLgLk1Su7AXipfQaAtTRpYp9qOv649QKwFSLPdVrA1iBwHhTL7Zqqc6RzEtgagDwdN7dtWyc0MZ16vJsGGpbcjEXn7a5HihH/3NTd9/riJHXwofXJ75kuan9z16x3QWx1ANqlDoDhsTq9jm0rbl9RdC7Nuqw3Vnb8Pg0ult7neXix3PHNxz5rpkDfrHb3aq3X+zuXUTuVSa4fV6dZY01er7Jox+Lx3CadTIOQuX44XD8duXxyoD+/3djeFTbQT0+uG1kuG+uEJi7W4101ULDvNYuTA/19txWzMgcfdCagap5wpJ2VQ6gmr1eAcRBlfvQqi6arJ6au3Ca1e/SpNB3l4+25YnY8l2hqOgb7gZUDvUFfkiRJkjRr0Rn70Rn0JUmSJEkz54z95ORAb8CXJEmSJE2NM/aT50y+JEmSJGlqvCt+9Qz6kiRJkqTKOGM/fS7ZlyRJkiRpOHMV7DNn8iVJkiRJY/PmefPDoC9JkiRJGpYfdzfHXLIvSZIkSerLGfv550y+JEmSJKknZ+wXj0FfkiRJknSeH3e3+FyyL0mSJEk15oz98nAmX5IkSZJqyGvsl5cz+ZIkSZJUA87YLz9n8iVJkiRpiTljXz8GfUmSJElaHjF687zac8m+JEmSJC0wl+IrcyZfkiRJkhaQS/HVi0FfkiRJkhaAM/YalEv2JUmSJGkOdbzGXkNyJl+SJEmS5ogz9hqXQV+SJEmSZqgm19g3Zr0DdXLwwTUOPrh2PuBLkiRJkjQug/0MdM/Y33fbOYO+NEEh/dHsWAbzxzKRJNVW7Az/dwAhhHtDCM+GEP60x/NvDCEcCyF8Kf39xa7nbg8hfC2E8FgI4T2T+G8OtBQ/dwYicRK/U11cri9pWRgcpf7sU0nSlFW3FP/XgQ8AH9rhNf85xviW7gdCCE3gg8CbgMPAwyGE+2OMXx1nZwYK9p58pqffHfbtEEg7s27Mjsd+flk288OymA/lgUjLRVpiFQX7GOODIYQbR/jW1wCPxRgfBwghfBS4A6gu2Pdq9AyX1SsH+qsalwBwY9gA4L+2X5z6Pi2jfu/lXs+/3DkLwL7GxrbHrRPjy8e8TdEIr6YrhjbTdjNt52PdTsulrm7uBuD5zmkAVkLxulZN7oQ6CeVjvxGKU8Tx9hkArlopjvHhzeMAXLuyF4CX47ltr89lko/8aiqLTiozy2R4verFuVh8hM9aaKbn47btfaE4l7wcNwHYSI8f7xRlZps1vH5lsZKOcT62+X3fSN9nH2q2dodVYOucsml7NDNXp77tidQ+nUxfrSOatBhn+nF3rw0hfBn4NvAPYoyPAtcBT3W95jBw67i/aMdg36tCWdGm79nOqeIrxVeX7E9Gv/dyr+f3NtYH+n4NL5Y6we3SdvmYN1NozIE+dwgMj8MrH/scVDYaRUc4h8F9zWJAK3fCyq8v6/W4BterXuT3f7tUL/Ixfz6e3vb42diqdD/roF9ZlNuohjPDcyW3W5q9Z1Lftsw6ookbYcY+hHAPcE/XQ4dijIeG/DF/BLwixngihPATwO8BNw+9MwPy4+4WVL8l+1Jd2SGQJEnSeSNM9qQQP2yQL/+M413//kQI4V+GEK4AjgA3dL30+vTYWAz2C86b70mSJElSDzP6HPsQwjXAMzHGGEJ4DcUn0r0AHAVuDiHcRBHo7wLuHvf3GeyXjEFfkiRJkpKKLs8MIXwEeCNwRQjhMPBeYBUgxvgrwE8DPxtCaAGngbtijBFohRDeDTwANIF707X3YzHYLzmX7EuSJEmqreruiv+2Ps9/gOLj8C723CeAT0xyfwz2NeFMviRJkqTaqckNlQ32NWXQlyRJkrT0ZnSN/bQZ7AW4ZF+SJEnSEjLYq46cyZckSZK0NFyKLxn0JUmSJC0wZ+ylC7lkX5IkSdLCcMZe6s2ZfEmSJElzzxl7aXAGfUmSJEmaDYO9KuGSfUmSJEkz51J8aXzO5EuSJEmaGZfiS5PnTL4kSZKkqTHYS9VxJl+SJElS5WKc9R5MhcFec8GgL0mSJGninLGXZscl+5IkSZLGZrCXZs+ZfEmSJEkj86740vwx6EuSJEkamDP20vxzyb4kSZKknrx5nuZRIAAQqccbdFDO5EuSJEm6gDP20uIy6EvShRwcliTVjsFe88jO2Ghcsl8fBpfZswzml2Uye9aP2fL4SzXkzfM0T/KJKPOENJqqZvLtKMxeuQwsk9nxmM8P68H8sSwkabpipx7trsF+QdgRqMakgr7lM3vlMrBMJOuBVGadkGrIpfhS/bhkf/E5Qzl7rp6YP5bB7FgfJGnGXIov1Zc335OGV75kSPPDEDk/LAtJmrKKluKHEO4F3gI8G2P8gYs8/7eAfwQE4GXgZ2OMX07PPZEeawOtGOMt4+7PjsH+O5vHAHjtxvUA/LfNFwHY3SjCTduT09S83DkLwN7G+oz3pJ7Kgf53fmoVgP/h350C4KqV3QC80D4NwFooqpYduNHlkNimGGVt0gAunPV6sX0SgH2NDQC+fPwJAG6/tGhfj7RfBmA1NKew18shH9tWbAOwq1G830+cb4c2tm1vpPf7tc2iHjTSz3mucwbYOldsphFz68XoBq0X52KreD4Uzz9++lkA3rjnewB4OtUbz+OjG7QsbmleCsBvHv8KALfufSUAz7eL80euP5bFdPyXF/4MgNde/moATnaKgftLm7sAOJfavYYDlZX7ystPAvDWA98PmDO0kH4d+ADwoR7Pfwt4Q4zxpRDCjwOHgFu7nv+RGOPzk9qZHYP9q9evAuDL554Dtho9K9r05dByKp2Ackdb1eq1ZDIH+t/9oaLz/NYvnQBgf7MoJzsGoysf816d5by9Jw123bSyD4C/dNl/D8Afbj4DwIFUd2y3+isHlbU0GJI7vvlY50Cf26H8fv9WqxgMPtoq6sflq3sAOJC+r5PKoBUti2H1qhe9QmU+Xz/fKgL8X9tThMmXUtlZH0Y3bBv18TPfAuC2vTcDW4ONuf60arJEdF7ceVUxKfZ46yiwVVdy/2o9FOXiAGT1/saBYgD+kXPfAeCKNDhs+6SJq+ga+xjjgyGEG3d4/vNdm18Arq9kR5Idg/2x1AHYnxo9Tz2zk08w6w2vnpimXif2PAv2t79SnIR2N4ra8eG/VsxQumR/dP1uglfezmXxVJqBzPaeD5MaVD62jfMBv9jOK1BygM/bOZCcH8BKZXHt2v5trz+aziWZA17D61UPGj0GH19Kq4fyzNeTrePbXp/rjYY3bBt17eoBAJ7tFANe+di3HeCaiSfTwEpux3I5uNJu+r7ZLtqlS5uXAJ6vVaERgn0I4R7gnq6HDsUYD42xF+8E/qBrOwKfCiFE4FfH/NnAgNfY28jNDzvE8yEv6y7PtHht/vT1CjaZ7df4Bv3EgVwWFwR+TV1uo/LgjJeizE65vvRrs1StcrtkOcwPy0KVGWEgNQXtscM2QAjhRyiC/eu7Hn59jPFICOEq4NMhhD+LMT44zu9xyF6q0MEH1zj44Nr5gC9JkiRpijqd4f9OSAjhB4FfA+6IMb6QH48xHklfnwV+F3jNuL/LYC9NQfeM/X23nTPoS5IkSdPQicP/nYAQwncD/x742zHGr3c9vjuEsDf/G3gz8Kfj/j4v2JamzOX6kiRJ0pRUdJPSEMJHgDcCV4QQDgPvBVYBYoy/AvwicDnwL0MIsPWxdlcDv5seWwE+HGP85Lj7Y7CXZiwHegO+JEmSNGEVfY59jPFtfZ5/F/Cuizz+OPCXJr0/BntpTjiTL0mSJE1WrOjj7uaNwV6aUwZ9SZIkaUwVzdjPG4O9tCBcsi9JkiQNqaJr7OeNwV5aMM7kS5IkSQNyxl7SIjDoS5IkST14jb2kReSSfUmSJClxxl7SInMmX5IkSbXnNfaSlokz+ZIkSaodZ+wlLSNn8iVJklQXfo69pFow6EuSJEmLzWAvaRuX7EuSJGlpuBRfUp05ky9JkqSFZ7CXpC0GfUmSJC0c74ovSb25ZF+SJElzzxl7SerPmXxJkiTNq2iwl6ThGfQlSZI0Nwz2kjQ+l+xLkiRpZvwce0maHGfyJUmSNHXO2EtSdQz6kiRJqpzBXvMkEABopO029XiDzrtcLtHyGNu4S/Yti/ll2cxO+dhbFrPjsZ8PlsP8sCw0LTHW4z22Y7C3ws2PHOibofhXO7ZntzOik+pEM9URTc6kZvJtvyZn2HDY63nLZHqCbdPcsh5Ihdy37VgnVDVn7DVP8i0fgo3eXGjYaZ6aUYO+HYTJKR/LUY+tZSJJkqbOYO/SvXmSj32rJktJ5p11YnZ6Ldm3LKan3/t/UgMBGt2gZaPp8djPB8thfnhZqabFz7EH1kITgOPtM8DWMvCNUHybFXJ6rmjsAuCasA7An7aPznBvlseoS4pb6VKIm1b2AfCt1nEAVlKdseMwul7Lvttp3UozLd7Lgf51jcsBeNVm8bqPr54A4MnWMQD2NTa2/Tz1Vj7Wq+lYn4ktAPY0ivbnWPs0APubRbt0Om4CsC89v56+77lO8bq8wiWfU855KdHQBq0X+fn1dJ7ObVU+9t88+zwA161duu31GtywZXEu1Z8rmpcUr0sD9C/Hcxd9vaq1kvqyrViUVy6ffP52Rd70XN/cA8DhdnHe3kztlTlDE1dRsA8h3Au8BXg2xvgDF3k+AL8M/ARwCnhHjPGP0nNvB/6P9NL/M8b4G+Puz47BPgf63Hk7mjpzoVE0ep3UKNoIjq9fwDwazwJbwV6zdW06GT26+QIAB1J4bOMql3H1OmaNHsf0kXgUgAfCKQBWWkWn7YG3FMH/J3+/6CjkgUn1lo9tPtb5EqDm+Y5wcSxXUwc4b+fXn+gUQeXF1FG+LAWZS9Kp5mQaADiHwX5Y/epF2Woqs6c3iwGuq1b3AnD7rpsAeDQNDlsSO7tYWz5sWeR68EzrJAA3pAHho52ij5XLyrIYz6Dn3Rzoc8A/2i7apUubq8DWAI192/H1K5PHW0cBuLK5G4DvtIqAb86Yntrcj6W6j7H/deADwId6PP/jwM3p763AvwJuDSFcBrwXuAWIwBdDCPfHGF8aZ2d2DPYbjaKRy7MrlzSKjnLLijZxg56I/rhVBMncsdZ4+h33Xs8fab8MbAX6zR4zNRpfv2OZ26c8M587yT/5+0Vn7XfeUIRMP0ZvePnY5/d1HrjKM1vl2ZR8U8l8rjjeKQYkT3Bu2/OeO8bX73KHY2lg/vKVYhDybKeoB1/sFDP2nkMGM0hb3q8sTnaKtmhX6lM92ykGIdeclZyoQc+7eVY4t0eXpoGX8sCmxtevTNZDUSeOpXOFOWP66tJfrWopfozxwRDCjTu85A7gQ7G4Lf8XQggHQgjXAm8EPh1jfBEghPBp4HbgI+Psz0A3z6tLoS8CO2PzYbUUbHrNJmt6tpa9bp+hn9Rd9uts0Gvmcz3o1SmzkzY9+f1/fnAmbXsGmb5yffF8MVvlfpTlMDvej0VTM0KwDyHcA9zT9dChGOOhIX/MdcBTXduH02O9Hh+Ld8WXVCsGfUmSJO0khfhhg/xMGewl1Vqvu+xLkiRpCVR3jX0/R4AburavT48doViO3/3458b9ZQZ7ScKZfEmSpGU0w4+7ux94dwjhoxQ3zzsWY3w6hPAA8M9CCJem170Z+Plxf5nBXpIuwqAvSZK0BCqasQ8hfIRi5v2KEMJhijvdrwLEGH8F+ATFR909RvFxd/9zeu7FEMI/AR5OP+p9+UZ64zDYS9IAXLIvSZK0eCq8K/7b+jwfgZ/r8dy9wL2T3B+DvSQNwZl8SZKkBTK7a+ynymAvSWO4WNA35EuSJM2HaLCXJA3r4INrzuJLkiTNC4O9JGkULteXJEmaD87YS5ImwqAvSZI0IwZ7SVIVvMO+JEnSdDhjL0mqlDP5kiRJ1TLYS5Kmypl8SZKkyTLYS5Jmwpl8SZKkCYlh1nswFQZ7SZpzBn1JkqTROGMvSZpLLtmXJEkaTOw4Yy9JmmPO5EuSJO3MGXtJ0kIx6EuSJNWTwV6SlpRL9iVJUt1Fb54nSVoGzuRLkqS6cil+l0AxyhGJle6MJKl6VQd9zxlSb9YPSZoub56HJ595ksuikbbblslcsI7MD8tidJNesm8ZzI9yvbCezI7Hfj7kcsgsj9mxTmhaYk3eYgPN2HdShWtQj9GOebQaiki/FpoAnOicm+Xu1F47renJ5dHx5DRzedBrM5VBMz1imQzugpn8158F4M6HNoDBw2GvMJlZJtW7cDCYi26rerkPtZrKpHy+MNxMVzMUx7td6unb152+3Ic6F4sWyTqgqjhjz1YHoJUav3xM7DBP30Zq/C4P6wCcwGA/S22KYN8MqwB06nLxzhQNOtOYO2ProWjOznRawNZgmAFmePlYH3xoPW0X7ntdEfTf9vki6A96bPO5xFoyvkHrRT7mzVQPznU2AdjVKAZvztiRHtuwqyFyWeSLPW2bZiMH97zysWMdmJnct22l9siBeVXFYM9Wo5dH1MqPO8o8PYc3jwPw5c2XAbhp48pZ7k7t7W0UweZEpwg6ucPmyWhy+nWW8+N5FuxY+zQA+5pF2ZxOQWYltV+WyeB6Hfsc6D/yujMAHPz8+kVflztpK6XZmBDyLHI9TrBVGLRebKZhlDzQtTsF+pdSPdlorFa/s0tu0LJopq8nSyvt1tJgpG1Ttcrlk8/b+fhvpK854OeZfMulen++eRSAvY3iXJLzhjlDk+ZS/C55lqVcsaxo07OnWTR6Gw0/yGAe5OCSA33Dk09lenWe89c867WW6sZmzKspHGQZV/nY587WnSngk7Y/dltasv9g0U6tlAaDz89UamL61Yvzgyfp2Od6seY5ZOJ6tTHlMnIAeDbKx7ncPp0fePQ8PnV5wHGlFOgzy0KT4ow9cKa0dC8vdT2VRp2dCRvfoEv5VlNHIAf7fCLScAY9cfcrl+uauwH44pmnAdiTZon3pgGYXD7OTI6vX1k83zoBwKvXrgBgd2qnHmsdBeCS1H61vFyip0Gvic+Pl2fkc6C/LM26vHc93YTv+EsArKfZ4T3p+dVSx1rDG7TM8kqWfF+Q/c1d6fntYVODG/Y8kpcbv9Qu6sXpWPStclnksrFeVOOCGft2MRD56rXLga0Z+rzK5VQsVrmcTV81vl515onTzwLwvZd8F2DOmIa6DmD5OfZsjaTlkFLuzNXtTVGFQVdB5EgSPOZjGfQ9269cnusUS5EvX90DwK50rX3usDkjMzn9yuKylWKQ5blOEWAeTZetvGq9uFzleFp26SBLb4POOObt8jkgP3/+2vwU6P/jXynOIT/+xeLnGVwmZ9BzRx5MOVGqB/myCA1v2PNIvtltvvxhN/k+B/l+INaLacorIPN9Jk6m8/bWMnAHgSetV5155a6rga18Yc6oXl2PaV3mdga6xr58XWRd3xSzlO/iuk660YgnnpnKHYI1Tz6VG3R0OV8neb5DUJcLqirQK9BvLclPlzuUBrCeaZ8E4PKVYsDrbV8qymJ3o6gv+S77+aZ8Gl2/epFngfMAfbYrhUtXfU1PvnZ7pXRJi4F+Osp15EAa7Moz8/kylXZdev5zqLxyyD6VJq1T0Yx9COF24JeBJvBrMcb3l57/58CPpM1LgKtijAfSc23gK+m5J2OMbx13fwa62M6ZrtnLy4gN9PNh6yZh5TtRO/g1af2OZT7mOahckUJlnoWx/Rpdr1nhXu/zrWuIt8sB5uBD25+577a0ZL/0MXvqr1+9OH83/FQvVks3MtT0rDoAPFeO9/i4YGvG9HkpkKaliqX4IYQm8EHgTcBh4OEQwv0xxq9u/d7497te/78CP9z1I07HGH9okvvkXXQkSTORA70BX5IkVaWim+e9Bngsxvg4QAjho8AdwFd7vP5twHur2JHMYC9JmqlyoDfoS5KkSRnl6swQwj3APV0PHYoxHuravg54qmv7MHBrj5/1CuAm4LNdD2+EEB4BWsD7Y4y/N/xebmewlyTNFYO+JEmalFFm7FOIP9T3hYO5C/jtGLddD/eKGOOREMIrgc+GEL4SY/zmOL/EYC9Jmmsu2ZckSaOq6OZ5R4AburavT49dzF3Az3U/EGM8kr4+HkL4HMX19wZ7SdLycyZfkiTNiYeBm0MIN1EE+ruAu8svCiH8BeBS4A+7HrsUOBVjPBtCuAL4q8AvjbtDBntJ0kIy6EuSpH6quCt+jLEVQng38ADFBwLdG2N8NITwPuCRGOP96aV3AR+N2z+H+S8CvxpC6AANimvse910b2Bhp896fvMNt/vZLJKkhWLAlyRpcJ966pNL/dnAf3Lj/zh0pv3BJz6+cMfEGXtJ0lJxJl+SJGUVXWM/dwz2kqSlZtCXJKm+qliKP48M9pKkWvEu+5Ik1ccon2O/iAz2kqRaciZfkqTl51J8SZJqxJl8SZKWj0vxJUmqIWfyJUlaHs7YS5Ikg74kSQusJpfYG+wlSRqGS/YlSVoczthLkqSeugO9IV+SpPnkNfaSJGkgLteXJGk+dWa9A1NisJckacJcri9J0nyIOGMvSZLG4Ey+JEmz1anJ3fMM9pIkTYlBX5Kk6eo4Yy9Jkqrkkn1JkqrlUnzNpZDemLE2n8goScvPmXxJkjQOg/2CyIG+kbbbs9sVdXGgZX5YFvPLshnepIJ++dhbFrPjsZ8PlsP8sCw0Ld4Vv4sVb/ba6S3ZikUZNENjp5erYtaJ+WNZVG/U971lM75Rl+yXj71lMX2eL+ZLuRxCaYmu5TQ9uW/bqMkyac2OS/EvwpPT7MQU6DdjMVdvsJ8PnVQXPCnNTm6Xcgehmda12E5Nj+eG6blgJv/1Z4vHH1rf8fvK9cQ2a3rK9cLzxmz1O/62Z9OT+7aE4ph77FUVZ+wvwoo2OznIh2BHYJ7YMZu9WOqk2U5Vp9dy7vLj5derOuVA32tGv1xPND3l+uLg42yV64DlMDvlSSrLQlUx2NM9qllonh/x91q9aVsLTWCrDE7H1ix3p/bOpeN/aXMXsFVHXmifBmAtFFXLujG+ftcI53ZqV1gF4OHj3wTgdftvBuBY+wzgKpdRXBhIiu3NHqsjTnSK2ePdje2h8sXWSQD2p/piWYyvV73Igf4jP1i0RW//yh4Anjr30rbXX7t6YNu2RtevjWqllXZbA/MG/FnI5XJZoxgMO5HO4+dS+VgO05dXoW6E7XHEnKFJcyk+8D3NfQD8p1NPAHDt2gFgqwK2rWhjG/TGRvnxb7eOA3Bp85Jp7WIt9TuZ/PmZ5wH47j2vBODluDmdHauhfnUjd47PD0Cm0HjGwa+R9WqXcgd4V6MYRDndKd73K2ngMXu+dQLYKovceeucDzQaV7968de/mNukItDfd2A/ALc9/S0ArlzZCzjIMgm9yiJvP3b6OwD8zb3fB8DXOsV5PNef1WCNmKY8c3c2nSN29QiVub3T8Abt2/5I8yoA/u2prwFw3fqlgDlDk9epR64nnL++5SLefMPtEbyr7jzIx7yZRvxbsS6LSuaTq1fmz1b4LDprdgwmr9+5oFkaAChfOuQy8OkpXxaRZ43zIMzHXl+sZOl3bb7G531A5kuv/lRuv3LvyvKp3tZg7/Y6YZ9q+j711CeX+gT9H665e+g30x3f+fDCHZOBhupj+tNrW9XLx7wVO4b6OdAmbguM1onZy2WwGpqshuYFZaTx9TsX5GPeDA2aoUGDsO2vpieW/uQyydsHH1rn4EPr3HfbufPX5asa+ZjneuD5YrZ69ady+2X5TE+vOmEZaNLiCH8HEUK4PYTwtRDCYyGE91zk+XeEEJ4LIXwp/X1X13NvDyF8I/19+1j/wcTPsZckqaYuuMv+kB+nJ0nSvKtiSjSE0AQ+CLwJOAw8HEK4P8b41dJLPxZjfHfpey8D3gvcQjGO8MX0vS+Ns08Ge0mSBBj0JUnLp1PNp4q9Bngsxvg4QAjho8AdQDnYX8yPAZ+OMb6YvvfTwO3AR8bZIYO9JEm6qBzoDfiSpEVV0YUd1wFPdW0fBm69yOt+KoRwG/B14O/HGJ/q8b3XjbtDBntJkrQjZ/IlSYtqlKX4IYR7gHu6HjoUYzw05I/5OPCRGOPZEMLfBX4D+Osj7M5ADPaSJGkoBn1J0qIY5ePuUojfKcgfAW7o2r4+Pdb9M17o2vw14Je6vveNpe/93PB7uZ3BXpIkjcUl+5KkedWp5pN5HgZuDiHcRBHU7wLu7n5BCOHaGOPTafOtwH9L/34A+GchhEvT9puBnx93hwz2kiRpIpzJlyTNmyqusY8xtkII76YI6U3g3hjjoyGE9wGPxBjvB/63EMJbgRbwIvCO9L0vhhD+CcXgAMD78o30xhFi7P1fffMNt/shkpIkaSIM+pI0fz711CcrmdKeFx+67meGzrR/58i/Xbhj4oy9JEmaCpfsS5JUDYO9JEmaKpfsS5KmZZS74i8ig70kSZopg74kqSp1ubbcYC9JkuaKS/YlSZMyysfdLSKDvSRJmkvO5EuSxuVSfEmSpDniTL4kaVgGe0mSpDnkTL4kaVDRpfiSJEnzz6AvSerFGXtJkqQF5JJ9SVJmsJckSVpgzuRLkvy4O0mSpCVi0Jek+vHj7iRJkpaYS/Ylafm5FF+SJKkGnMmXpOVlsJckSaohg74kLQ+vsZckSdK2JfuGe0laLF5j3yVQHI1Ym/GO+dNJxz7G4mszNGa5O7VnnVAd+b5X3XWH+kFn8a038y2XT2Y5TU/u2zaoSerSzLgU/yI8Oc1ODvTt9NZsYrCXYKtdKtcN2ylpi+fvyRt0ub7HfL70C5PWlenJfVtCccw99qpKXd5ROwZ7K9j8WAlNANZDUWStWJexp/mU64Qj/bOXj3k50Nt+TV75WA56jC2L6evVNlkW1el3h33PF/OhfK6wHGZnLfVtc4/WslBVOjV5b+0Y7J/dPA7AjWuXAXCWNgDnYvHVmbHpybOR1zT2APB0++Qsd2fplTu/5e2/0izqxMfPfAuAq1b3AdaJSRo0gDy3+TIA+1d2bXv8suYlAJzsbA70c7SlvApiNb2vN0urIsqrJA401gE4E1sAnOgUAWe9UZxqTneKx1dTZ07D61cv8vOtuP18/eSZ5wD4qb3fB8Cj7WM7/hz116sscqDPz//O3UXb9JbfKvpUub7kAXvLYLpeaJ0Atjr6+5tF+exqrAJOnEzTp5/5EwB+8tpbADgZi/O1OUMazY7BPoeVZ1OIvHalCJWnU0d5NeTOnSalV0chN26XYId4HKPOVpVf//VY1InLU53Io86ejCan16BKefvAShHgDzQ2AHi6VQT9jWZRJsc5C3gN307Kx7S8CmKzFPBzcN9IK4hyN/jb6difSZ2zK5q7t/0cA/34yvWg1+VZIS1tbabXff/u6wFY9zKusfVrk8rbd/xWESSvTPXhX9xQBPx3PlnUn7OddN7w3jmVKJdHDvJ5wPGyUJw7jsbiXLGSysGAPzm9+l6vu/IvAPB8+xQAl6eyMWdUp66rtupSmwe6xn5PmoV5qX0GgLXUmWvX7E0xDb0qWn780c0XANhII8sazqANWb/X5ZNQHuE30FenfCzL27kMXo7F7PBlKeg/3S4607m9skx669fuNM6HyGK71zlgT6OYqVxPp5Y8ABBT0M/fp/GVy6bXYHCuHyc6RWj5TOswAJemFS0aXr82qbydZ+ZfbJ8G4J1PFuHlw68ttg9+fr2S/VShXB7lwbDnOkU5NENe7VKXCDA9vc4xe5rFez8fc3NG9eraF6rL/3rHXtYFHYVgaJk1A/18yDOP+WTUq3Ot6pXLIG+vusy1Mr2OabsUNhvO0M9MLqOzcfvlDwb66Tu/AqY0E1wO9IPeZV/jyaEx142sFT1XTFt5EMWcoarUZbjO6RNJkqQZG/Qu+5Kk4fg59pIkSZqJfnfZlyQNpqq74ocQbgd+GWgCvxZjfH/p+f8deBfQAp4D/pcY45+n59rAV9JLn4wxvnXc/THYS5IkzSln8iVpPFXE+hBCE/gg8CbgMPBwCOH+GONXu172x8AtMcZTIYSfBX4JuDM9dzrG+EOT3CeDvSRJ0oIw6EvScCq6xv41wGMxxscBQggfBe4Azgf7GON/7Hr9F4CfqWZXCgZ7SZKkBeWSfUnaWUVL8a8DnuraPgzcusPr3wn8Qdf2RgjhEYpl+u+PMf7euDtksJckSVpwzuRL0sWNEutDCPcA93Q9dCjGeGiU3x9C+BngFuANXQ+/IsZ4JITwSuCzIYSvxBi/OcrPzwz2kiRJS8agL0mFUZbipxC/U5A/AtzQtX19emybEMKPAr8AvCHGeLbr5x9JXx8PIXwO+GHAYC9JkqTeXLIvqa4qWor/MHBzCOEmikB/F3B39wtCCD8M/Cpwe4zx2a7HLwVOxRjPhhCuAP4qxY31xmKwlyRJqgln8iXVTRWxPsbYCiG8G3iA4uPu7o0xPhpCeB/wSIzxfuD/BvYA/y6EAFsfa/cXgV8NIXSABsU19l+96C8agsFekiSppgz6kpZdRXfFJ8b4CeATpcd+sevfP9rj+z4P/HeT3h+DvSRJkgCX7EtaPrGapfhzx2AvSZKkbZzJl7QsqpqxnzcGe0mSJO3ImXxJi6qim+fNHYO9JEmSBuJMviTNJ4O9JEmSRmLQlzTv6jFfb7CXJEnShLhkX9K8cSm+JEmSNAJn8iXNC2+eJ0mSJE2AQV/SrPhxd5IkSVIFXLIvaVqcsZckSZIq5Ey+pKo5Yy9JkiRNkUFf0qQ5Yy9JkiTNkEv2JY2rE52xlyRJkmbOmXxJo6pHrDfYL5xAAOpzrYg0KuuKJC2vUYO+5wapfvwce0laQHbapN6sH1pW55fsv/5ssf3Q+kVfV37vWyek5VeX+j1QsC83ejaCs+Mxl3ZmHalev3NAfj6zTOZHO91CqFEqI02PfapqlQN9OeiXj7d1YvasE6qaN88DfqB5AIA/OPU4ANes7QdgLRTfZsWbnt1hFYCXY7HUzBPQaIY9efR6/r+88GcA3HnVLQB8p3MKgHOxDVg+4+hXRuXtVzX3AfBE5wQAZzstAJqhAcBmKpPV0Kx83xddPra5o7uWjtnpziYAuxvFjNiJTtFRviRt5/d9fl07Ft+/p1l0pK9s7AJgM/3c451zFf4vllOvepDLqklj2/OXNzaACzsz+RzSTjcS8vw9vGHL4sXWSQD+8vo1ADzWOrbt9SupnnneGE+v8/WdD22kfxWPv69dtEe/sla8fi2V14upbpxJ7dnZ2Kp0fwX/uLUbgHe2nwLMGdPQq/1adi7FB56IRVjZu1I0grmiWcEmr194yR3iqxqXAPB85/S0d3EpTGo0+LWXvxqAx1tHga26kTtmnoxGVz5m/bYfax/ftp0b7/UU7OvSmE9Sfh/n8JcHRTZjDvwr257Pr89fX7V2ObDVQd6XBia/bbs1sl5tSaNHW3MihZLnWsWA175mEW7ygP2fdYp606rJnYJHdbG2fNiyON8mpQC53ijqTx4AaKV60ratGkuvcik/ngP9P2yeAeAdp4s6kssh15Xczml0/fq2/2a9qAMHNouAb86oXr8+1bKqy/9zx1brxXYR7K9d2QvA6dRRMLRMXr+KlmfEnukUI/95NlKzkWcyc0egXCesG+MbdHXFmTRLvNbYPrhSDp3qr9cM5GoKIJulmfzcLuWAsj91iHOgz69/LhYz/Lm+aHSD1ovybOOu1GH+VjqHRAP9QHZqywcti2tWi1VFj7ZeAmBXoxjoatdmceh8yYNeP3OqWDmRV0x86geLr3c+av9qUvr1bb+R6sQVzSLYnzVnqCJ1aW13DPZ55Ox0qYNgRZsdA/1kDPoe7vW6Vmnm0joxeYOOKm+kTvJ5wQ7BuPKxOz9IUtrO7//yDGUO+ucM8JUZdrYlzz6Wz+Ma36BlketLPn9bP2brWLqU6IqVPdsev/urxddmGgv24/Sql/tQ5YFIz9+atLoMZpsSJUmSpC4HH1zj4INr5wO+JM07LyCSJEmSLqJ7xt5ZfGkx1eV+SwZ7SZIkqY9yoDfoS4uhLtfYuxRfkiRJGlL3cn2X7EvzK47wZxAhhNtDCF8LITwWQnjPRZ5fDyF8LD3//4UQbux67ufT418LIfzYJP6fzthLkiRJI3ImX5pvVSzFDyE0gQ8CbwIOAw+HEO6PMX6162XvBF6KMb4qhHAX8H8Bd4YQvg+4C/h+4LuAz4QQvjfG8e6uarCXJEmSJsSgL82Xiu6K/xrgsRjj4wAhhI8CdwDdwf4O4B+nf/828IEQQkiPfzTGeBb4VgjhsfTz/nCcHTLYS5IkSRXJgd6AL81GRdfYXwc81bV9GLi112tijK0QwjHg8vT4F0rfe924O2SwlyRJkirmTL40G4NeM98thHAPcE/XQ4dijIcmtlMVMNhLkiRJU2bQl6ZjlGvsU4jfKcgfAW7o2r4+PXax1xwOIawA+4EXBvzeoRnsJUmSpBlzyb5UjYqusX8YuDmEcBNFKL8LuLv0mvuBt1NcO//TwGdjjDGEcD/w4RDC/0Nx87ybgf867g4Z7CVJkqQ54Uy+NFlV3BU/XTP/buABoAncG2N8NITwPuCRGOP9wL8GfjPdHO9FivBPet19FDfaawE/N+4d8QHCTiMYb77h9kqGNyRJkiQNzoCvqnzqqU+GWe9Dld54/Y8OnWk/d/gzC3dMnLGXJEmS5pwz+dJoOtUsxZ87BntJkiRpwRj0pcHUI9Yb7CVJkqSF5833pIur4hr7eWSwlyRJkpaEM/nSdgZ7SZIkSQvNoK+6q+jj7uaOwV6SJEmqCZfsS8vJYC9JkiTVjDP5qguX4kuSJEmqBYO+llU02EuSJEmqI5fsa1l4jb0kSZKkWnMmX4vOpfiSJEmS1MWgr0XjjL0kSZIk7cAl+5p3zthrLgUCUJ+bQEijsq5IF7JeSAXrwuQ5k695VZd6brBfMHV5Y0rjsq5IF7JeSNsZ8Ktj0Ne86LgUX5IkSVo+Bvnpc8m+ZqUu9X3HYH8utgBYC8XLmmlUs50OjqOc07M7rAJwoFE0gkfaJ2e5O7V3U3MvAH907hkALm1eMsvdWWrldqa8fbJTdBB2p7pxurMJwDUrewA41jk7vZ1dMuVjnc8B52IbgJXQ3Pb6jbR9Jj3fSl/3NzcAOJHKSuPrVy/KZfXU2RcBuPuSVwPwh50Xt71eo+tXFvsb6wCcTn2q3GatpvqSr/1spO/TZJXLI9eFA6u7Abi0uQuAJg1gq85YN6qTA/0TZ54H4OGf2A/AXZ8tHjdnaNKcsQfWU5hs0wG2OgqavlOpQ9Dq1OONOe+ebJ8AtgK9J5/qlI9peTsH+tw+rTeKZu1k3JzC3i238rHupK/N0Ljo6/Ixz3efzcE/D7ZocvrVi/Md41DUi+/ZuBKAL8fjU9i7eulXFi91zgCwVhoIy+eN7bVJVfvu9cuBC8spB3pNz40bVwDwP3222M415L7XFwPydz60MYO90jKqS/98x2C/NfJfnHYcQZudfKxPp46zI/uz1e7TkdP0tUvtVe6kWVcmp9eMZJaPfT7k+Xm7y9NXDo25fhxPK1jKgzOqTrlN2kirIO1TTUe/4+q5YnZC6ZjnOpED/cdeXwyKHXxofbo7pqXjjH2XfqPRmh5PPNLOcvtkXalOr3OA54b50assDPTTV26THBierfLx9lwxO/3OJeVA77X5GlVd2llvnidJkiRprnmXfY3KGXtJkiRJmkPeZV+DcsZekiRJkuaYM/nqJ8ZO/xctAYO9JEmSpKVg0FddGewlSZIkLSWX7Ksz5aX4IYTLgI8BNwJPAAdjjC+VXvNDwL8C9lF8gNA/jTF+LD3368AbgGPp5e+IMX6p3+812EuSJElaas7k11ec/s3z3gP8vzHG94cQ3pO2/1HpNaeAvxNj/EYI4buAL4YQHogxHk3P/8MY428P80v93BtJkiRJtXLwwbVtoT4HfS2fDnHov2O6A/iN9O/fAP5G+QUxxq/HGL+R/v1t4FngynF+qcFekiRJUq0dfHCN+247Z8BfQjHGof+O6eoY49Pp398Brt7pxSGE1wBrwDe7Hv6nIYQ/CSH88xDC+iC/1KX4kiRJkmrP5frLaZTPsQ8h3APc0/XQoRjjoa7nPwNcc5Fv/YXujRhjDCH03IEQwrXAbwJvj1u37/95igGBNeAQxTL+9/XbZ4O9JEmSJJUY9JfDKJ9jn0L8oR2e/9Fez4UQngkhXBtjfDoF92d7vG4f8PvAL8QYv9D1s/Ns/9kQwr8B/sEg+2ywlyRJkqQ+vMP+YprBzfPuB94OvD99/Q/lF4QQ1oDfBT5Uvkle16BAoLg+/08H+aUGe0mSJEkakDP5i2XaH3dHEejvCyG8E/hz4CBACOEW4O/FGN+VHrsNuDyE8I70fflj7X4rhHAlEIAvAX9vkF8adhrBePMNt0/9KEiSJEnSolnUgP+ppz4ZZr0PVbpi3/cOnWmfP/71hTsmzthLkiRJ0picyZ9Po9w8bxEZ7CVJkiRpwgz682EG19jPhMFekiRJkirmzfdmYwbX2M+EwV6SJEmSpsSZ/Olyxl6SJEmSVCmDfrW8xl6SJEmSNFUu2Z+s6FJ8SZIkSdIsOJM/Gc7YS5IkSZLmgkF/NF5jL0mSJEmaSy7ZVzeDvSRJkiQtKGfyd+Y19pIkSZKkhWLQ386l+JIkjSAQgPqMkEvDsH7MlsdfdVT3JfsGe80lT0jSYKwrs+Mxl3qzfkialV4z+cuuLq1uqMsIhiRJkiRJy6gx6x2QJEmSJEmjM9hLkiRJkrTADPaSJEmSJC0wg70kSZIkSQvMYC9JkiRJ0gIz2EuSJEmStMD+f7dHvR7Jna6CAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1440x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA9wAAABUCAYAAABwWz81AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAJYklEQVR4nO3dfcid510H8O8vL+02my62EbukXVe1k+mQoXRaYVocQt3W9g/FrfjWwcicDhEUdIjtGBMURGR2WoOrtQy7qpORzjiZf0gdUrqgW92Lk+6lNO3mtsa+pHle8/z84zx5eozNc07sc3LInc8HLjj3fV33dV3JSXLON9d1P3d1dwAAAICttW3eEwAAAIAhErgBAABgBgRuAAAAmAGBGwAAAGZA4AYAAIAZELgBAABgBgRuAAahqu6oqt+e9zw2U1X/VFVvm/c8AICzY8e8JwAA06iqryT59iSrSU4k+VySu5Mc6O617v7FOU4PAOD/sMINwLnkhu7eleTKJL+b5DeSfGC+UwIAeH4CNwDnnO5+qrsPJnlzkl+oqldX1V1V9d4kqapvraqPVtU3quq/119ffvL6qrqqqu6vqmeq6h+r6v1V9cGx+hur6rNV9eT6NvBXjdV9pap+vaoeqqqnqureqnrRNOMCAOcXgRuAc1Z3P5jkSJLXnVK1LcmfZ7QS/vIkC0luH6v/yyQPJrk0ybuT/NzJiqp6ZZJ7kvxqkm9LcijJfVV1wdj1P53k+iRXJfm+JLdMOS4AcB4RuAE41z2e5JLxE939RHd/uLuPd/czSX4nyY8mSVW9PMk1SW7t7uXu/kSSg2OXvznJ33X3x7t7JcnvJ3lxkh8ea/O+7n68u48muS/JayaNCwCcfwRuAM51+5IcHT9RVS+pqj+tqkeq6ukk9yfZXVXbk+xNcrS7j49d8ujY671JHjl50N1r6/X7xtp8bez18SQXTTEuAHCeEbgBOGdV1TUZBeFPnFL1a0m+O8kPdvfFSX7k5CVJvprkkqp6yVj7K8ZeP57RlvCTY9R6/WNTTGmzcQGA84zADcA5p6ourqo3JflQkg9297+f0mRXRvdPP1lVlyS57WRFdz+S5HCSd1fVBVV1bZIbxq79qyRvrKrXV9XOjEL0UpJ/mWJqpx0XADj/CNwAnEvuq6pnMtri/VtJ/iDJW5+n3R9mdN/1N5M8kORjp9T/TJJrkzyR5L1J7s0oVKe7v5DkZ5P80fr1N2T0OLLlKeY3aVwA4DxS3T3vOQDAXFXVvUn+o7utSAMAW8YKNwDnnaq6pqq+s6q2VdX1SW5K8pE5TwsAGJgd854AAMzBZUn+NqPncB9J8o7u/rf5TgkAGBpbygEAAGAGbCkHAACAGdh0S/mj17y+k+TYkxcmSVZXtydJLrhwNUlSZXX8bPmeL46eePPFV78qSbK04G6Aedq9dyFJcuQ/dydJLt61mCRZa4/aPduOHL04SfLHLxq9B7+3aylJcuG3rCRJFp/ZmcR7sxXWTvzv/6PdsePE6Pz67+1fH780SbJvdXT8mh1Pj+rXRscXXbR0VuZJ8uyzo8/tj9euJMllo4/t3Pzp9yRJvvy6X5rLvM5Hl37X6PPigQf3JklenNHfm70vPZYkqW2j71K+U50dHz6+J0myZ/Q25BUro8+K3TtGDyE4+Xnu/Zi9y99+ZZLkax/4cpLk2WOjf7fkjLPvlZ//2KC/JK1880sbf5h27vmOs/5rldoAAAAYpF5emOv4AjcAAADDtPjsXIcXuAEAABikXhK4AQAAYMvNO3D7KeUAAAAM0/LCc2WCqrqzqr5eVZ85Tf11VfVUVX1qvdw6qU8r3AAAAAxSn9k93HcluT3J3Zu0+efuftO0HQrcAAAADNPi8ambdvf9VfWKrRzelnIAAACGaXlxo1TV/qo6PFb2/z96vLaqPl1Vf19V3zupsRVuAAAAhmlshbu7DyQ58AJ6+9ckV3b3sap6Q5KPJLl6swuscAMAADBMiwvPlReou5/u7mPrrw8l2VlVeza7xgo3AAAAw7QFQfukqrosyX91d1fVazNawH5is2sEbgAAAAapl5anbltV9yS5LsmeqjqS5LYkO5Oku+9I8lNJ3lFVq0kWkrylu3uzPgVuAAAAhmlhceqm3X3zhPrbM3ps2NQEbgAAAAapF6cP3LMgcAMAADBMZ7ClfBYEbgAAAAapF5bmOr7ADQAAwCD1cSvcAAAAsOXWFlbmOr7ADQAAwCD14om5ji9wAwAAMEhrx9fmOr7ADQAAwCCdOD7f8bfNd3gAAACYjRNLtVEmqao7q+rrVfWZ09RXVb2vqh6uqoeq6vsn9SlwAwAAMEgrC9s3yhTuSnL9JvU/keTq9bI/yZ9M6lDgBgAAYJBWFrdvlEm6+/4kRzdpclOSu3vkgSS7q+plm/UpcAMAADBIy0s7NkpV7a+qw2Nl/xl2ty/Jo2PHR9bPnZYfmgYAAMAgraw8t7Ld3QeSHDib4wvcAAAADNLSylT3bk/rsSRXjB1fvn7utGwpBwAAYJCWVndslC1wMMnPr/+08h9K8lR3f3WzC6xwAwAAMEhLPf0ac1Xdk+S6JHuq6kiS25LsTJLuviPJoSRvSPJwkuNJ3jqpT4EbAACAQVo8g03d3X3zhPpO8stnMr7ADQAAwCAt1nzvoha4AQAAGKSlbTXX8QVuAAAABmmxBG4AAADYcgtzfi6XwA0AAMAgLcx3gVvgBgAAYJiWBW4AAADYegvVcx1f4AYAAGCQFrM21/HnfAs5AAAAzMZSeqNMUlXXV9UXqurhqvrN56m/paq+UVWfWi9vm9SnFW4AAAAGaTEnpmpXVduTvD/Jjyc5kuSTVXWwuz93StN7u/ud045vhRsAAIBBWszaRpngtUke7u4vdfdykg8luemFji9wAwAAMEiLfWKjVNX+qjo8VvaPNd2X5NGx4yPr5071k1X1UFX9TVVdMWl8W8oBAAAYpOV+bkt5dx9IcuAFdHdfknu6e6mq3p7kL5L82GYXWOEGAABgkBZ7daNM8FiS8RXry9fPbejuJ7p7af3wz5L8wKROBW4AAAAG6QwC9yeTXF1VV1XVBUnekuTgeIOqetnY4Y1JPj+pU1vKAQAAGKTxLeWb6e7Vqnpnkn9Isj3Jnd392ap6T5LD3X0wya9U1Y1JVpMcTXLLpH4FbgAAAAZpaW1l6rbdfSjJoVPO3Tr2+l1J3nUm4wvcAAAADNLiGQTuWRC4AQAAGKQzWeGeBYEbAACAQVoWuAEAAGDrLZ0QuAEAAGDLLa0K3AAAALDlltcmPn97pgRuAAAABml5zivc1d1znQAAAAAM0bZ5TwAAAACGSOAGAACAGRC4AQAAYAYEbgAAAJgBgRsAAABmQOAGAACAGfgfXXe4ch0FHlEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1440x72 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "id_test = posterior_cov_inverse_prior @ cov_posterior\n",
    "plt.figure(figsize=(20,5))\n",
    "sns.heatmap(id_test, xticklabels=False, yticklabels=False)\n",
    "plt.figure(figsize=(20, 1))\n",
    "plt.title(\"Diagonal\")\n",
    "sns.heatmap(np.diag(id_test)[np.newaxis,:], xticklabels=False, yticklabels=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "e3865f77",
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute central value posterior, Eq. 47\n",
    "v1 = u.T@y\n",
    "w = uup.T@fk.T@ u\n",
    "cv_posterior = uup @ np.diag(1./llp) @ w @ np.diag(1./l) @ v1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "7c4c4e42",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(450,)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cv_posterior.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "e8027c1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# draw samples\n",
    "samples=np.random.multivariate_normal(cv_posterior, cov_posterior, size=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "e87d25ff",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000, 450)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "samples.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "ad7773e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "f = np.array([samples[:,i*50:(i+1)*50] for i in range(0,9)])\n",
    "flavour_cv = f.mean(axis=1)\n",
    "flavour_std = f.std(axis=1)\n",
    "    \n",
    "flavour_labels=['Sigma','g','V','V3','V8','V15','T3','T8','T15']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "bbfa1fd3",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([\n",
    "2.9699999989970886e-05,\n",
    "4.1246302562114121e-05,\n",
    "5.7279854290892427e-05,\n",
    "7.9543107354176351e-05,\n",
    "1.1045380151089824e-04,\n",
    "1.5336544212678162e-04,\n",
    "2.1292711419771039e-04,\n",
    "2.9557946342023174e-04,\n",
    "4.1023635832702205e-04,\n",
    "5.6921788919482173e-04,\n",
    "7.8952001737071158e-04,\n",
    "1.0945283774603655e-03,\n",
    "1.5163052361650959e-03,\n",
    "2.0985919006607802e-03,\n",
    "2.9006588501226953e-03,\n",
    "4.0020755222241709e-03,\n",
    "5.5083182403814438e-03,\n",
    "7.5568314649709480e-03,\n",
    "1.0322652614357851e-02,\n",
    "1.4022011128953883e-02,\n",
    "1.8911583033933310e-02,\n",
    "2.5280751342878608e-02,\n",
    "3.3434951223099259e-02,\n",
    "4.3670479005884652e-02,\n",
    "5.6244653828324422e-02,\n",
    "7.1348265292604018e-02,\n",
    "8.9087633304453406e-02,\n",
    "1.0948053859041228e-01,\n",
    "1.3246538983220960e-01,\n",
    "1.5791909766096124e-01,\n",
    "1.8567792651252230e-01,\n",
    "2.1555676356233966e-01,\n",
    "2.4736436150447358e-01,\n",
    "2.8091396862547191e-01,\n",
    "3.1602987267868177e-01,\n",
    "3.5255081809432581e-01,\n",
    "3.9033126797342838e-01,\n",
    "4.2924130793198417e-01,\n",
    "4.6916577293070166e-01,\n",
    "5.1000298771219632e-01,\n",
    "5.5166336656838788e-01,\n",
    "5.9406801715516011e-01,\n",
    "6.3714742681882575e-01,\n",
    "6.8084026849030499e-01,\n",
    "7.2509233853930255e-01,\n",
    "7.6985562506376237e-01,\n",
    "8.1508749786096646e-01,\n",
    "8.6075000810050684e-01,\n",
    "9.0680928480266942e-01,\n",
    "9.5323501557082535e-01,\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "ae6057a7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(50,)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "2f5d865a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fecb1851250>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD8CAYAAACSCdTiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABfbElEQVR4nO29d3wc1bn//z4zs03aVZcs23KRuy13bJopNgaMaQZCTQECgYRLuGkkgSS/lC/JTW5uCmn3cgkQyE0hhJCE3ksoptjGNu69SLItWb1sm5nz+2NW0q606ruSLZ83L+OdM2dmzo53P/vMc57zPEJKiUKhUChGJtpwD0ChUCgU6UOJvEKhUIxglMgrFArFCEaJvEKhUIxglMgrFArFCEaJvEKhUIxgjOEeQDwFBQVy4sSJwz0MhUKhOK5Yu3btUSllYbJ9x5TIT5w4kTVr1gz3MBQKheK4Qgixv7t9yl2jUCgUIxgl8gqFQjGCUSKvUCgUI5hjyievUCgUyYhGo5SXlxMKhYZ7KMOK1+ulpKQEl8vV52OUyCsUimOe8vJyAoEAEydORAgx3MMZFqSU1NTUUF5eTmlpaZ+PU+4ahUJxzBMKhcjPzz9hBR5ACEF+fn6/n2aUyCsUiuOCE1ng2xjIPUiJyAshviSE2CyE2CSE+LMQwiuEKBVCvCeE2CWE+IsQwp2KaykUCsVw8IMf/ICysjLmzp3L/Pnzee+99/jMZz7Dli1b+nyO4ajfMWifvBBiLPDvwCwpZVAI8RhwLXAh8HMp5aNCiPuAm4H/Gez1FAqFYuJdz6T0fPt+dFGP+1evXs3TTz/NunXr8Hg8HD16lEgkwgMPPNCv60gph/yJJFXuGgPwCSEMIAM4BJwDPB7b/whwWYqupVAoFEPKoUOHKCgowOPxAFBQUMCYMWNYunRp+yr9Bx98kGnTpnHyySdzyy238PnPfx6AG2+8kdtuu41TTz2VKVOm8Prrr3PTTTcxc+ZMbrzxxvZr3HbbbSxatIiysjK+853vpGzsgxZ5KWUF8BPgAI64NwBrgXoppRnrVg6MHey1FAqFYjg4//zzOXjwINOmTePf/u3feOONNxL2V1ZWcs899/Duu+/y9ttvs23btoT9dXV1rF69mh//+MdceumlfOlLX2Lz5s189NFHrF+/HnDcQWvWrGHjxo288cYbbNy4MSVjH7TICyFygVVAKTAGyAQu6Mfxtwoh1ggh1lRXVw92OAqFQpFy/H4/a9eu5f7776ewsJBrrrmGhx9+uH3/+++/z9lnn01eXh4ul4urrroq4fhLLrkEIQRlZWWMGjWKOXPmoGkaZWVl7Nu3D4DHHnuMhQsXsmDBAjZv3twvX39PpCJO/lxgr5SyGkAI8QSwBMgRQhgxa74EqEh2sJTyfuB+gEWLFqmq4gqF4phE13WWLl3K0qVLmTNnDo888kifj21z8wgh2l8DaJqGaZrs3buXn/zkJ3zwwQfk5uZy4403pmzhVyp88geAU4UQGcKZUVgObAFeA66M9bkB+GcKrqVQKBRDzvbt29m5c2f79vr165kwYUL79uLFi3njjTeoq6vDNE3+9re/JT1Pd9E1jY2NZGZmkp2dzZEjR3juuedSNvZBW/JSyveEEI8D6wAT+BDHMn8GeFQI8f1Y24ODvZZCoVAMB83Nzdxxxx3U19djGAZTpkzh/vvv58orHTt27NixfOMb3+Dkk08mLy+PGTNmkJ2d3eU8tm0nPf+8efNYsGABM2bMYNy4cSxZsiRlYxfDEbfZHYsWLZIqn7xCoejM1q1bmTlz5nAPo0eam5vx+/2Ypsnll1/OTTfdxOWXX57Qp6mpiUAgMKjrJLsXQoi1UspFyfqrFa8KhUKRAr773e8yf/58Zs+eTWlpKZdddlmXPt1Z8ulEJShTKBSKFPCTn/yk1z7DIfLKklcoFIohYjjc40rkFQqFYoiQUg650CuRVygUiiFCibxCoVCMcIbaL69EXqFQKIaI4bDkVXSNQqE4/vhu14VGgztfQ2rP1wPKXaNQKBTHIPfccw/Tp0/njDPO4LrrrutTyGRnlCWvUCgUxyAffPABf/vb39iwYQPRaJSFCxdy0kknDehcSuQVCoXiGOPtt99m1apVeL1evF4vl1xyyYDOo6JrFAqFYoSjRF6hUCiOMZYsWcJTTz1FKBSiubmZp59+ekDnUT55hUKhOAZZvHgxl156KXPnzm2v7JQslXBfUCKvUCgUvTGEIY9t3HnnnXz3u9+ltbWVs846a0ATr91Z8rZto2npcawokVcoFIo+cOutt7JlyxZCoRA33HADCxcuHNB5kol8Oq37lIi8ECIHeACYDUjgJmA78BdgIrAPuFpKWZeK6ykUCsVQ88c//AExSGtbSpk0rUE6RT5Vzwe/AJ6XUs4A5gFbgbuAV6SUU4FXYtsKhUJxXCJNMzXnGWJLftAiL4TIBs4iVsNVShmRUtYDq4C2cuaPAJcN9loKhUIxXMjoCSryQClQDfxOCPGhEOIBIUQmMEpKeSjW5zAwKgXXUigUiuHBjA7q8DYhPx5F3gAWAv8jpVwAtNDJNSOdd5D0XQghbhVCrBFCrKmurk7BcBQKhSL1DNZdczyLfDlQLqV8L7b9OI7oHxFCjAaI/V2V7GAp5f1SykVSykWFhYUpGI5CoVCklmXLlvHCiy8ltN17773cdttt/T7XcRddI6U8LIQ4KISYLqXcDiwHtsT+3AD8KPb3Pwd7LYVCoQCY88iclJ7voxs+6nH/ddddx2N/f4KLrru2ve3RRx/lxz/+cZ/OHx8ff9yJfIw7gD8KIdzAHuDTOE8Jjwkhbgb2A1en6FoKhUIxpFx55ZV86xvfIBKJ4Ha72bdvH5WVlZx55pl9Ot62Q4AbOE5FXkq5HliUZNfyVJxfoVAohpO8vDwWzZnDc889x6pVq3j00Ue5+uqrEUL06XjLakXTXMDx6ZNXKBSKEc9VK1fy6KOPAo6r5rrrruvzsZYVbH+ddDGUbQ1+gN2gRF6hUCh6Qdo2Fy9dxiuvvMK6detobW3tV+4ay2pNsNY7W+7SDKdsrJ1RIq9QKBTdYLdZ2LaNP8PHsqVLuemmm/plxdt2BCkTY+y7ivzgYvB7Qom8QqFQdINtOiIvYy6Wa66+mg0bNvTPVWOHkdJOaslLKYlG6pF2albTJkNloVQoFMcdvYU8pgrbilnyMVG+7OKL+zVJatsS2wrFtrqKvGU1Y0brsWXfJnAHgrLkFQqFohtsK2Zhxyx5mWTStOfj7Vj4JEjZcWybyEej9dhWEEn6RF5Z8gqF4oRAStnnkMc22iz5dnG3+hkFIzsiazqLvJQWptkI2MlzvqQIZckrFIoTAnsAuWc6u2tkP0XeyR8faTtbQrtlhRxfPYBQcfIKhUIxKMxo/yNYrE7uGvrrrrGjtPnik1nybQhNxckrFArFgJHSHpQlP1B3jZRm3Gsr7rUSeYVCoUgZtmV1xLz38zjnRf8mXjtCJDuuacdZ8rZtI+nYpymRVygUioFjW3a/I2Oc4xzxPffSS3np7bcTLPl7772XT3/60yxcuJD58+dTVlbGfffd136ctO0eLXnittH6P7a+oqJrFArFccfWGTNTer6Z27Z2aZO23T7Res3ll/PX557jggsvat/flmr4lFNOwePx0NzczOzZs7n00kspLMh3zpEg7N375NOJsuQVCoUiCbbdsUr1iosu4oU33yQScsIh41MNezweAMLhcHvyMWnZjjVPvCWvRF6hUCiOGWRM5KWU5GVnc9Ls2Tz/6qsACamGDx48yNy5cxk3bhxf//rXGTNmDLZtOZkle7Tk05fKIB4l8gqFQpGEdqs8JvRXr1zJY089BSSmGh43bhwbN25k165dPPLIIxw5cgRp210mV2WnOHlpx4V0pnE1VMpEXgihCyE+FEI8HdsuFUK8J4TYJYT4S6xqlEKhUBwXtE/USgm2zcXnnMPrq1d3m2p4zJgxzJ49mzfffNOJ5rEsiHPX0NmStyJxRx8fuWu+AMTPXvwn8HMp5RSgDrg5hddSKBSKtNJWyEPGRN6fkcFZJ5+ckGq4vLycYNDx09fV1fHWW28xffp0x9VjW50s+bgEZVY0wbJPJykReSFECXAR8EBsWwDnAI/HujwCXJaKaykUCsVQ0NldA3DVBRckpBreunUrp5xyCvPmzePss8/my1/+d2bPLotZ8jbEiXyCJR8NIfuZR2egpCqE8l7ga0Agtp0P1MuOmYVyYGyKrqVQKE5wkoU89sTRg/sxoxGKJ03t8zEd7hobLOf1pcuXY5smQtcBOO+889i4cWP7MZFIDbYdxbZthJQIV7y1Hp9qeGiseEiBJS+EuBioklKuHeDxtwoh1ggh1lRXVw92OIPHtqHywz513fjKC7Q2NqR5QAqFYrDYlgWSfq16lXGWfIIV3kNqA8sOYdthsJJFz8SJ/NAY8UBq3DVLgEuFEPuAR3HcNL8AcoQQbU8KJUBFsoOllPdLKRdJKRcVFhamYDgDxLZgw6Pw36fCI6vA6jm8affa93np/l/x28/fxCsP3UdD1ZEhGqhCoegPUsqOHDRW3y3odneNLRMTk/Wwcta2QkgZAVt2Kc4tEkJo0plcOJFBi7yU8m4pZYmUciJwLfCqlPITwGvAlbFuNwD/HOy10oIVhbWPwK9Ogr9/Fo5uh3ADHHin20Nsy+Jff3gIADMcZv0LT/PQF2/lmV/+F1X79gzVyBUKRR+w4yxvux+pDeIteWnHT5omt+SdtMIhJ7WwBDpPrIrhseTTmdbg68CjQojvAx8CD6bxWv0nGoJ1v4e3fwGN5V3373gBSs9KeuiGl56ltjLxGNuy2Pb2G2x7+w0mzlvI4ks/xvjZ89IxcoVC0Q/iXTTd5q+JtII7I6EpIbpG9m7JOwW7bUfkbTddRD7eej9eRV5K+Trweuz1HuDkVJ4/JURaYM1D8M6voflw9/12PA8rftClOdTSzOrH/9zjJfZtWMe+DesonjyVxauuZNopSwY7aoVCMUDiLe+kIt9QDq01UDgDDE97sx0/8Wr37pNvK/Nn22GQXUU+wV2TxiIhnTmxVryGm+DeufDit3oWeICaXVCzu0vzu0/8hWBTY58ud3j3Tp762Q9pqa8byGgVCkUKSHTXdBLohgpoqXaEvG5/ewUo6N5d011O+bYyf7YdRRMaogchH6LoSeBEE/ldL0Pr0b733/F8wmb9kcOsf/6pfl+2fOumfh+jUChSg92dJd9YCS1VHdvRFmjuCKDomHi1ueCmTzuphuPOce+993Lbbbfxta99jbKyMubPP4Ovfe1HSGkjNInQEkU+QfTF0IVQnliphrc927/+O56H025v3/zXHx/CGkB1mfKtm5l+2pn9Pk6hUCTnN597dYBHHkzaevu3c50XTYfBkwXujA5L3ra5auVK/vrcc5y3ZEm7Jd+Wavjuu+9m48aNNDVv4/zzPs5bb63h7AUrkpjQHSIvjqfomuMGy4SdL/TvmP2rIeS4Zsq3bWbne91H3PREhbLkFYrjBAnNjnUfL/KXn3++k2o4GkXadnuqYZfLRSgUIhRqJRRsIRo1KSrKB93u4q6J3+7JlZNqThyR3/8WhPq2cMkKxxxmdhR2v4KUkjd+/8CAL3304AFCzc0DPl6hUAwhVjghl7y0bSfV8Jw5vPDmm2BZ7amGTzvtNJYtW8bYseOYPn05y5efzvTpk0Czuqn25JxTiXw66IerpmJ1LsFal7Ox4wW2vfU6h3fvHPClpbSp2L55wMcrFIohxAwn+O7bXl998cU8/vzzSNtuTzW8a9cutm7dyt69m9m69WX+9a/3eeedtaCZ3UTQtIn8cZTW4Lhhe99EPljrouWwl8Nrs5ESottf5s0/PzLoy5dvVSKvUBwXSAtpxud6d4T50vNX8Pp77/Hhhg3tqYafeOIJTj31VDIzM/D7MzjvvDN4//0NSN1MOrnaFlXT2ZJPp11/Yoj8oY3QkHzCpTM1W/wAhGrcNOzJYG25l6aafkTkdEOFEnmF4rjBNsPtr9ss+UBONmctXsxnv/719iyU40vG8cYbb2CaYaLRKG+9tYbp0ychtWgPlrxEaJ3DMNMn8yeGyPfRig836jRVeNu3D27O4f2jJSkZwpG9u4mGQik5l0KhSC8yXuTbYud1g6tWrmTjtm3tIv+xS69g8uTJLFx4FkuWXMWcOdNZuXJpt+cVQiKEPZQLXk+QEMptz/SpW83WQEJSia15hUTl4G6R5faiRSPYlknljm1MmDt/UOdTKBRw+33n9KmfbVnt+aQ0wwZcFE2YBFVbweze6Ip317SlBRaGzqXLlxPcvBnvjBlOmw333Xcfwdo9mK7WPoxIoiWdkE0fI9+Srz8Ihzf22i3aotGw3wfA3gkrOTBqJgfzAr0c1TvhohLMzCwAyrepUErFyKHl/feHewi90r4QSoBmyI4slFa022OkDXaCyMcmSw0jtt/GrHNWsUeaWjiyZxehOMu/J4QArYurJr2MGJGvr1+TfEcfXTU12/1gCypGn8He0ovZUlI86LXHps+PFcjBDOQAyi+vGFnUPvJ7rGM8NLhN5DXdWYHqpCiwQHYvtLYpkHGpxttFPq7KU7SiArO2FmG3zaT2fUxK5AfIkapnaG7e0XVHH1w1ZlijfncGtTlT2DZpKWbwXWzzwKDHFC5y/PmWPxsJHNq1A8vs3oJQKI4X7GCQlrffJrxt23APpUfactUI3W5PM2BHe7a6raiWmOOmzSdvJv6ghQ8dQqAhEF1SGHSHQCqRHyh1h5vYtfNniY3Betj/drfHtJoudjXl8fLOKbw1qZR3J+iEm/+MGRrYytZ4ov4c7AwnUkcaLmxfJmYkzOFdA4+3VyiOFVreeQcZChHammaR72v+d5lcZNss+fhcMtKM9HzJqEiMk2+z5GMJyNowNQ0hNDSh99mSF2LoRX7ETLy2HG2hOfoyhw6uZfS4k5zGnS+BbRK1NZqiHppNNzXhDA4Fs6gMBmiIOj54dCADoP95aZIhgUhRYklb05+DHmyhfOsmxs6YlZLrKBTDRdOrTu6Y0Nb+1VrtN5XroGRR7/2irU6aYC1R0hLcNXqbyPfsj5e2wI4LwGgXeTsCuNrbLV1HQK8ZJ+PRNGtIV7vCCBL5SG0V7hzJ+2/ciV5zGdFQPU171tHUdCohy9X7CVKImZ2P7fEltgVy8FRXULFN+eUVxzfStml+/Q0AQv0sqN1v9rwOo+f33s+KOoWAMvMTmtsteV0itLa27kXethxxT3gwaNsQIDSt3cq3Y8W8+2PJ61pqDMn+kIpC3uOEEK8JIbYIITYLIb4Qa88TQrwkhNgZ+zt38MPtnnC0HgD/mAMc3P1P9qz7gOp6a8gFXgpBuHBsl3bb48N2eajYvrX76jQKxXFAcP0GrJoaACI7dyGjaZxnqlgLofre+0kLgl3rNsgEkY9Z8j3Ub5YxkU+w5IELbrqJF15/p10xpRDc9/DDfP7uL3HVDZ9iwsTTufrqzyec67bbvsWcORdwxhlXccYZV7Fx4zaEPvQinwpL3gS+IqVcJ4QIAGuFEC8BNwKvSCl/JIS4C7gLpyRgWrDifq5GL65m99MT0nWpHonmFiFd7qT7zEAOkdojVO3bw6hJU4Z4ZApFamh+rSPNr4xGCe/e3R43nnIq1iYV759ec3FKL/OVn98DdIi87FSE9aqVK/nLky+wbNEyJGBqgn88/TT/+c0f0hpppkU7xMMPP97lvPfc82Uuu+z8uJahddVAagp5H5JSrou9bgK2AmOBVUBb0pdHgMsGe62esN0dZbsCY1vxj25J5+WSIjWdcP7obveb/hxA5bFRHN80vRITeZfzlBzakiaXTf1Bp4hHsD49509Cm7smfk2qlJIrVpzPM6+8RSTiTNrurTzE4SNVnHHK6Sw742wCgcwhG2N/SWl0jRBiIrAAeA8YJaU8FNt1GBjVzTG3CiHWCCHWVFdXD/jahkh8KCk+efD5ZvpLJL8YjO4fjqyMTKSmK7+84rglsm8fkT3OCtLc664FIdLnl69Y6/ydxJJPF23h8wkueSAvO4uT55fxwpv/AuDx557lsosuRgiBpnXvkL/nnl9x+ukf4+67f0w43HNUT7pImcgLIfzA34AvSikTiqBKZ3o66XOKlPJ+KeUiKeWiwsLCAV+/88n9xa34Zhr9fjjyaBmcWngJEzL7FwFjGy6ieUl/xzoQGqY/m3Il8orjlKZXX2t/nXvddbhLSwmny5IfYpGXdoebRnaeSRVw3WUX8NdnnwPg788+y5WXXu7s6mbR5He+8wXWrHmS1177M3V1Ddx770PpG3wPpETkhRAuHIH/o5TyiVjzESHE6Nj+0UBVd8enizELKmktnYXl69uj1BjfZC4ouYkJ/lmcWnQJy4qvI9tV0KdjrcJxSK3322kGcgg2NlBT3resmArFsUTTq68A4Jk1E09pKb758wlt29aRxCuVVKxz/u7LxGsKaI+sSbZTwKoVS3lt9Xus27KFYDDEgrnznF3dGPLFxYUIIfB43HziE5exdu3wpDVJRXSNAB4Etkop41cjPQncEHt9A/DPwV6rvwQCNeSOPUrrhBkER0/E1pO7UnRhcFL++ZxZfCVeveMHocg3nvPH3siCvOW4RPLJVADpziCcndenMZmZWUgEH76wun9vRqEYZsy6OoIfrgcg+6KLAPDNn4fd3Ey0vDy1F7MtqPzQeT1UlrzVnuy96z4E/swMzj51Mbd9+9tcdvFFiDb57Cbu/fBhx/0speSZZ15l5szhCbZIhSW/BPgUcI4QYn3sz4XAj4DzhBA7gXNj20POhAkbQICZU0DL5NlEcgsTfqlz3aM4f8yNTMlakPR4TehMy17EypJbmOifnbRPRskCZB/jZNENrMwAm15/n+3vHuq9v0JxjND8xhtOEWshyFq5EgDf/PlAGiZfq7dBNBY8McQi350lD3DtJRfw0fbtXHbJxYhY4/LLLuaGG+7kjTfeY+bMc3n5ZWeV/Wc+cxennXYFp512BTU19Xz1q7cOwbvoyqBDKKWUb9H9UoDlgz3/YMn011NQuJ+j1RNBNwgXTyCaU4jn8AHmuGdTlnsGutB7PY/P8HNK4UVMDsxjbc1L1Ecc79P4wlPZ7Ok+F4bb3UI06kXKjmuY/hyMqgpeeWQrZtSm7MyucfUKxbFGcyyqxrdgAa4xYwDwTJmC5vcT2raVrBXn93R4/2jzx0PS6Jqv/OXpxIYjm8GKm9gMFBPSsqg/XImvINRunNumwKx34Q9GEbrEk9URt97urkmiZm1tq84/h8bNm2nxuBCx1VUvP/kPor6uQSNPP/1g7+9zCBgxuWt6YsKEDUDHAiTbm0Fw4gxqRhUQFl0XJ0i6X6xU4C3hvDE3sDD/PCYF5lGb2306Yl2PMHvOK4wenZg4zQxkg2zCMht4/U/b2Xj/7+DpL/X/jSkUQ4QdidDytmOhZl14YXu70DR8c+ekfvI1QeT7YMnbnb7HwXqkZTk5a2ICLSUIzVnoJKUj+FY0LlSyPbKmh8dyIdsFX7T3O7YXN54QIp+R0UjRqL1d2nfrh/mrZzUb9f3Y2ASzd3Go7AF2LbuD5oL13Z5PExpTsxZSXLSYQ3ryD6AQNjNnvUFmZgPjxm1Gi1vOLF0eLI8PaVaAhDfXTeCjfx2G+sFnvlQo0kHr6tXYra2g62RdsCJhX9vka0qJiXxDxIPd2ovIS+mExsRjhrCjofZ8NUiwoi4n3XB8t1bnCTshsiaZxrcJu9bxI6DFLHk5xLlo+ssJIfIA48dvTFoh3TaCHJzwNBtPv5MDp3yfxrFvYbtaqJz3a5oLEouNRDE5IurZopfzlrGVN13dWy9Tp64mN/cwAG5PMIk1n4NtdkxWvdP0SVrf+fNg3qJCkTbaQiczTzkZoyAx4sw3fz7mkSOYtbUdjT2kDuiVaNCp3ATsjk7kcHUvCxvt5Fkd7UhrR3phS8O2Y97pONWzLYEVEe2uGnCyS3amXcYFXS35JLpyLDFiEpT1hs/XzKjiXRw+NA2AQKCa4tE7KSzch653/ZBI3aRi3i8JbbyG6vpR1IpmGkWwT4mIxo/fwKjiPQltJeM2cejQVGzbWSVoBfKwayva95vSx5o3WjjrAgu03ucIjleiR47gGtXLegLFMYWUkubXHJHPikXVxOObN89ZFLV1K/4lS5zGNQ/B4s9AH8KKu1C5vt39ctCeQLj6IDk9DrAbkY+G0VyOPFu2gR0r5dlmzbd9lc2gju5xhDrkMogYSb5/cYE3Mub/EULDljbyGBf5E8aSBxg//iPGjNnGwoVPMX/B8xQX704q8O3oJq65j1JfsIFGrW8CX1S0mwkTu5YbdLvDjBnb8UhrebxYWjPS7rBSNjcsoXHNi/16T8cT9Y8/zu5zz6Pqpz/DDgZ7P0BxTBDatAmzqgrhchE477wu+/XsbGdRVHza4XW/h0MfDuyCMVeNlFB+1GJ/g6/n/t1Z8rbd7q6xbQPbduROaLJdqJ3DBWZII+QyCHcW+PaV9G0q3+HuEQhsaTEc+Wj6wwkl8h5PK5OnfECmv77Px+i6xaxZr5OTU9lr3+ycQ0yd9m63+0tKtqDrsQgAITD9OdhmhzVv4+K9p0emXz60fTuHv/8DZDRKzW9/y+4LL6LxhZH7gzaSaMsdn3nmmehZWUn7+ObP7wijrN0LRz5y0gQPhJjIV4X8hEIRDrVmOCX7uqM7S14KND3ekncEXOhdZTmkJRF4QBDz2ceH0Gu0pzKwsZS7ZuhI36+prlvMKnudLZuXUV+fPAFZRkY9M2f+q8dK7C5XhLFjt3HgwFynIasA+0g5untae5+dVZNYsH0PBdMnpfQ9DCdWcwsVX/giMhRqbzMPHaLiC1+g/owzKP7WN3FPnDh8AzyBqXnwQVzjxhFYtgzhSp6Wuznmj4+PqumMb948ah+J5SPc+pTz9+7X4Myv9H9QMZE/YDrlM20pMCM95H1J8gPQarqI2jpGLLWv448XSDtWqk+IhKTxSV004FjyMkxNXQNnrLoJkFRW1aLrOkX5RQRDrWT4PVi2iWlarFp1Lt/4xu39f89pZASJfHpxhP41Nm9eRkMnoXe5gpTNfhWXq/cERGNLtlBRMQPLchPJyMRt7Sf+qyXRefcv67n42yNH5A9/+9tE9u1Luq/lrbfYc8ml5N10EwWf+yyar5dHc0XKaHzhRar+6ycA6Lm5ZF10EdmXX4avrKy9T6S8gvD27Qifj8A5y7o9l2/+fCLf+x52MIjWJvIH33cmUV39+DdtOQr1+wHYHyqGWPFsK5r43Sq/680+na6p/VV9wnZnZ6FxS3H7a4FAFwaa0NE1D4buJnf0KNa8+CKGqOGrP36IQCCTuz77TZqj9USMo2RkG0SjUVasuIHzzjuDxYvn9fktO9dMHyeUu2aw6LpFWdlrZGcfbm/TtChls1/F6+1bamPDiFJSsgVwHgGjPhNphxL67K/MoXJHbbLDjztq//QnGp99tsc+Mhql5n//l90XXUTji8qFMxREj1Rx+Nvfbt+26uqo+8Mf2PexK9lz6SpqHvod5tGjNMdcNYFlS9EyMro9n2fqFLSMDMIfrobyD2InDcP+ftZLjlnxlhRU1HY8FfdoyaeYDCNAljsfvysHn+7GpXkQCKSMpTMXgCZp885nBpz7Eo2aRKNmtwnLeiZ9Mq9Evp/oukXZ7FfJzj4C2MyY+SaBQP8EeczYbRiGs0rWyB6DbXb197/7lwFOWvWCaZps2bKFPz36MPf99v/R3Nzc+0EDJLhpM1U/+s++j63yEBX//gUOfu427FCo9wMUA0JKyaG778ZqaEi6P7xjB1U//jE7ly6j+je/AZJH1cTTtigq9MYTgIS82JPontd6PK4LMZE/HAxg2h3uFNuWWObQVFVyad6k7RIPljScMoBCOpE1SCwryhlnXMWUKUtZtuw0Fi2aOyTj7CtK5AdAm9DPnPUv8vMrej+gE44176QbDmdmYJldkzsdqhDs3ZianPi2bbNr1y4ee/wRHnzoJjZv+SyFBT9gcun/8YdH76Kpqan3k/QTq6mJii99CTkAC6z59depvOvu9GQ2VFD3+9/T8k4fLGzTxG5oQMvKwn/mmb12982fT2hDzIo/+y7wZvd/8jUm8ruPdk07HqpPfw4bl+ZpX+TUFZ2I9COEU0pQSguJRDcEb731V7ZseYl16zaxZcvOfl9XKEv+2EPXTQoKBp4ueMzY7bhcQaKaxHYlF/N3n9iOtAcmdFJKDhw4wBP/fISHHr6eDRuvJy/3HqZMeZu8vEo0zUbTbKaUPs8fH7sz5UJ/6BvfIHpw4Pen6fnnqf75vakbkAKA0I4dVP3s5wDoBQVM+NMfybvheoS7+yyrgXPP7XF/G76ZkwkdqAHDCzMuhAlL4PAmx8/eV2Lphfc35XfZFW5s7NKWatzdWPFtmLbfyTopJDY2Usr2Fa85OVmceebi9gRlxwpq4nWY0HWTknGb2btnEe7sbGRLFCctfwe1h8O8fcO3KTF3g6YhhMA1dgyBlSvxL1mSNBqirq6O1S8/SFXL++Tk7iM7UEN29+l1EEIyZeLL/Pmv/851V/2SQKCHzn2k5uGHaXrp5QEdq2WPx7vw05gVH1Bz//24J0wg52NXDHpMCif/TOVXv4YMh9Gzsxn/4IN4p08jY+FC8m64gepf/ZqGJ590Mk3G0VNUTTw+7yHCDTqydBnCE4CJZ8L2Zx1rfs6V7f2kaSKSVVCr3QPBWkJBjaN0/WybtoUdDqN5PF2PTQECgVvvWeQlbsfmFhIbi6O11WQaEXJysggGQ7z22mq++MWb0jK+gaJEfhgZPXoH5eWzCGVm4Wk4hO4a36XPdm0Ouev/jtYWC7wGGv75JHpODoEVK8i66EK8Cxfy0ZYP2bDpEQK+tWTnVzGuqyHULULA5In/4i9/u42rr/hvsrqJhe4LwfXrqfrpz3rv2HUUuKetxD3jYoRmoGePRcuZwKHv/whXSQmZp5w84DEdl7S5qgY0iZec6p/9nPD27WiBAONiAt+Ga8wYxvzwP8i/+Saq7r2X5ped4iB6fj6Zp53ap/PrB1/GlWERyT4dD0DpWc6OOJFvXfchR374Q8Y/8Fv07OzEE5Q7rpqd+4qwk6yUtYUgXF2Nr6SkX++7r7RNsPaEJkT7HKnUohyuPsIXbvoqtm1h2zaXX76CCy44Oy3jGyhK5IcRXbcYP24Tu3efjFccBLqKfMhXQMWYMxlX8XpCu1Vfz8FnnqK8aTXmzgg5RRWUFPewYKQPlI5fzeP/uJkrVj1ITnb/hd6qr6f8y1+GaLRfx4nMInwnfRo9b3JCu2vMArRAMRV33cOEB3+JZ1Jpv8d0XLLrZXj5e1A8Fy791cBSA3Si5d13qX3kEURGBuP+93/xzS5L2s8zZQrjfv1rghs2UPXTn+GZMhmhd40hDwaDVFZWMn78eFwuF4SbYc9r+AozCIWLHJEfVQYZ+bDnDcCZp6m8806ilZUcuOVWxj/0ILrf33HSmD9+X30u3eUxiLQ04zVNSn7kzBFEQkFqKzrNaQnwZkfQXDbhsA87mI9u2uhWmLDbwO2vRbYIfI1RGn0dTwW9WfFRKTGAu77+b+huG7vZpGzmdN5667EejxtulMgPM8Wjd3KwvAxfdjbUPo/bU4ZHH4tL09EBXQAzVmH4/WiahhWsZ1/ebhqW2GROqiLbsyul45lQso6/P/0pLl75MIV5uQBEq6rQMzMRGRndhodJKan8+l2Yld0UQnG5KLjlFnI/9Uns5mas+nqs+gaCW1sJ7/WATC5kemA0voV3UPH1nzH+/v+HkZubkvd5TFK+Fl7+DuyLxYAf3ujkSL/8vkHlM7IaGqi8626E2824//4NGQuTF8iJxzdvHhN+/0h7+ona2lp27tnHhm27OFxZgdXagACE1891H7uMaeGNYIbwzZxNaNcBssF5CpmwBLY+CTW7OXTPr4lWOpFkoY0bOXjrZxn/wG+d0EzLhB3P0VLlpsrl73ZcUSEw6+pwxepBt9bXd+nj9kfRXE74pW0bCKnjDR8FKYm4CrCljtAtTL3jMycQuLTu3UASSVSCSxO0z8tq5jG/2hWGQOSFEBcAvwB04AEp5bBUiDpW0TSb8eM+Ylf4VAhkAHuBvQgp0BHo6HhcIQ5Mb8afVUVG7kEyAzX0rdjgwBg/ehPPvXAd55/7R8TvHqL2oVgBYiHQMjLQMjOdP3GvZTjcbcSGd+5cRn//HrzTYu6B3Fys7FHUfrCD8J7eJ5aFy4dr0pVUfuP3jP3559C96fHJdsv+1bDhT7DtGSg920m8NXFJ6s5/dCe88r2OlaIA/lEQaYGPHnOE/mMPQjflK7vFjMDOFzn0y39g1dRQ8ptfk3lq31wvR48eZe2m7WzZvov66kPYMoyW2UpubjOlU+vJ9Vfh1Q9RUVfGo3+tZbo7xMX48J11EY1//6jjRKVnwdYnqXvwFzQ9/1LCNYLr1nHwtn9j3P/exxNvfYs/ZoT4zME8GjKcf1/Lm0EktwhhW+25Zixdw6qpxcjLw7IsQi2JIcCuDBPD2/FEK6MGrmgLuuWELOu2jW1pGIaN7svDLaNE7TAuvWdXjSUtrLa0CLHMlmhWSt1p6SKtIi+E0IHfAOcB5cAHQognpZRb0nnd441Rxbs5eHA24bAfISwCgRoCgaMEso4SCBzt80Kr3ohG3IRaMwk1ewk3eDBbNbx5ETLzmsnIakpI1jZ21E5ef+lKpjzWQCArC7uxEaTEbmnBbkkcj+2VWAHQskBrAhFL9CF8Pgq/8O/kXX89Is7l0Lqxmvp/7MJu7Rr3rOd48E1yYdZHCO2PguV8oYTQENkLOfSdpxjznVVoGcmX4A+EcGsr4dZmAvmFHU8q9Qdgw6Ow4c/OhCACimbCln/C5iegaBYsugnmXQueAU5WN1bC6z+ED//YkX+lcCac/nlaJ15AdM+7ZD17M2LLP8CKwlUPg9F7lAvgRLX8/bPUv7ubpvdzGLs0gn//z+HVf8HYk2DsIvB3hCnW1Tfw1odb2L5zN9GWKvwZRyjMPsLM0ha8s+qwogdoW30aT0nuGkpO20JkTzb3tXycswsXk3PwWWof34HdapI57RREg5sjf3N8/Mbo0fjmzqXphRcAaH3vPZ6/dhk/uKQJn+nCU+UhMieXaN4orIyO+yoNF7bHi2ZHiWSGiTZsxbY9aIaNbTqfLd1j4cqMq9sgNbAMvOGOdSyuaCshSyDcApeRiQvHSu8tXNe0TaxYHpt2Xdft9upQvWHbOpo2OHfqQBHpjEUWQpwGfFdKuSK2fTeAlPKHyfovWrRIrlmzZkDXeuLhS8gef/z+djQ1ObZ5ZmZ9j/lvekNKaG3OpbUhQKjRINLgIlyrEa0RyGDHv7UubKYEjnKgNY+gaYCQeApMvKMtMoosMvNDZARqqK8t4oC5hEBjJYHIUbx6C8JnIQ3AkGBILKljmgYubxiPtxVdNzGjbsLRDCJ2gKiVhTT9CDMLfzCX4uqx5LeOx4+z3N0o9OGblYcvsB33nv91fLNCww5ZBD0raI0sIdxaiqXZ1GTtoCZ7CzU5O4Ewlu3CljoWBlK4kZoHzeXDcPlwu/343AFMyyQUcYTcCjdBNISwIghsBBZoIIRAWhIsEJaFYZm4ohYZ7kyyxy2gaPaVjJtQRrhqN+H3fofc8RQhGaHVnUnDqNk0FEwn4nJhmY0gLMACaWFLG2Qs1E6a6HYEzQrhCTWTVXsYfzBETrCFSNYsmvwrCNZm0lpdC24d4TYI11VD61HcRhRvVgDP1FNw5+TgzgkgAm40t47uceFyuzFcBoIooed/Q3j1y7TWGYQibgLTIrgKbaLSjS11p0KS1Ii6cgj6iokaDbh8R8nIrMbMDBHM0NAsiSdkI4IG2RNOx104D7e7IO5PIS5XDjW1b1K+85c0RvdjRFxsOHQh7p15ZM3MYkdmFcsrT8K9t5mmmq20hA5iXbiCxlAU12svYZgmBhouVzat+Xm0ZHio0i3C7g7XlN7cgLAtll91LVOmFuF2B7vY21JqWKYLw5VYhtM0XWhHvXgiHda+FDqNAR9efwh3SxGGSBLFI6Nd2pvNJkzhJ6BLollOEkFpCyItBp5A73NQ0agHl6v7MqFYLrJyZ/R6HoCtW7cyc+bMhDYhxFop5aJk/dPtrhkLxAdLlwOnxHcQQtwK3AowfnzXiccThf6umo0nFMogXDceb80sio8uoqhJo7xlB82tO2mKtsUodwj8mClTOH9WlPzd72FZO9jdlMem+mL2VecSrnbRAGSKImZkXsvkSXmUuuuQ0WlYIT9m0EckFCIYPUwoepiQ1UTE9mCRhaFHsEQUnzeAkS3xZIUg0IydUYftOwL+bRhjawhNNzlg6QSDWbSG8nFFPIzeX8PoYAPFc68l4+rfY7sz2LT9dfbteYVw8J+43ZVkZDg/gK11Yxi1eyWzwho2JpGITjhqELIMoqbhpJWVBshYukEhsYWJNKLoXhOv1yQzw8Lvs/D7JIbXhYjUY5evx7a8BHMW0uCdQoth0BQ9TMu+KnYd/D+2SmcOVEeiZVyIFrTQG03EUQ3drMMtPVjuLKKWhRUOYUWjTvk5QyI0G832oJtZuEyB5oZw7nQigQiHc+uI6hEi5rtYRRkYU/Lx+PJwudy4rEyor8E+sosGCdHgViLRTKzKEMK20AwDw+VBkxq6KTFCEiMqMKbORhMC4YlQ64+ABh5N4sPG0MJIPYhbr8IjDuKN1OMLm3irbXwHLbxhm4ZmH09WzCRouvj01JfIuGoVTL6yy2dvdPFljH73JRp2b+DgyUuYPeEZ7OwCth2cj7u6lNVVb5PX1Io3uwj32LMI7Tdp8diEJ82k0ZXMsnUEviDqJbeumYb6Q7jHV+Jxr8LjduYHLNNwRN0dRtNshLC7CDyAtIwEgQcQ0qKt4mdTpJYMI4BHz4h9QyStZiMRK0SOpyjBfWPaUdDBpkPQndKCfTOSbXv4akQM+8SrlPJ+4H5wLPlhHs5xgWm6aKwfjVY7hdyj85nWPAMfcY/xHsjzFDM37yyao3VUtO6isnUX9RzlzI/fwLzzL3TcEs3fQT/wLtPKP2Ba+Qc07DvA3uplaK5TyPWNQwhBsGUvGbkfUrH9CLXBsUTtYrK9eUwsGkemaw5Gq3DEtBNRK0qkppVIVYiIHSRih4hYIRrtw9SM30jW1CNk5dXg9zurGEM4sxF75O8IvfpX3K4gum6R6YPMWH6r2upiDm+ZytKMhfje3IS95CJkuBWt6Siuo4fQag6j+T24CrJwjcrBNSYP99h8Igf2E1y7ldadh9F8frwlY3EV5mNkBYiGPIQjOohi9AlnoZdOw2PYFDQeIefIPqIHagnvPoDV2IiRn4+MRkFKXKOLcZdOxr1gEi6jGteB52DfWmzpxza82C434SYfDc1ZNIX9tFiZhA0vpsuD6fIQbWkhfGAPvqq9lIwqpvDkU8g45RQyFi9CDwSwrCC2HcEw/Aihw5rfwdNfxAxpRAMLic6+nUhVDdXvraZi+zZqfS6OBtw0GR48wiLLbeP3+fAWn4SoG0eBsQ+P8SrbPLup1HIx7AImREYxJVxETnghLnEYn/42Xm0d2xoLePHQVMzYv+urh8Zz8RO3QKjemY+Ix7Zg+7Nkh8GT/Ssy393Ke1k/ZtrUt9C9zyPDOkaTjW5bNOQU4G3NZ2pdKXnMI9pQwgHq2adXUyUaQMAYK5eF5iTyvBGa5r1L4+gKIv7Kdj+4GdKJNBkIKdH0LDxeL6arAenqWu9BN92xT1YiLmkhEUghaTEbMWUUr+6nJVqPpUUw/CbRcBS35nynbGljY6MBtjARUo8VCpFoRt8ky7KGT2rTfeUKYFzcdkmsTdGJPdvm0VgRQAiJx3Djd2fhMzKdmXxhOaULhY07msXoutnMlHlofVgK7XflMj17MdOzF4NLklFXRPCjo3in56L5i7BKVtBadxLB8JVEIo0UxIUu2zLC6/ufoXF3M8W5LiYuGMek5VdSPGVGgo/dbo0S2lWBq7gILcNA8xmEwxbvvPgCm19bQ7S2ERmpoN2nWwc1G4owcovInhOmYGIVPp9jcQkBPm/c47WEmqqxHNkwgRkLr+Gq710MQPSTFdQ9/jiu0lG4x5+Ea/wEXGNGJ4yrM9K2CW3eTMvqd2lZ/Q7BZ538QN6yMlxjxhBZs4/I3r3YkQjeGTPwzZ+Pb+FJ5N10E+642Gxp20QrK4ns3k149x6a9+wmsg9k8xR00YTV2IzZ3IyrMAf3uGxGTy/CPXUm7hnzcU+aDO5YdbDqGkKbNxHatIngpk00PPUUVkMD3pkzyTjlZPScHKLlFUQrKohWVhItH4+MRNG9B3G99D1cZadTOGcmE0rL8YXeQ+hOit2K1iwqMhdTHili344PyC46SOYpF1LTdDnRimZOmigpLa6iSXuX1bUf8WroPWY3jOa8+rPJC36J2uhBin3bORTcgyVNtjcWMaupiknPfMUpqn3WVztu6r63IFhLZPSV1D5zmJ/N+Ccrj97D9n/eS2aOiXtCGZl5+8kr2ECmUYedU0trzk5aS18EU8NVU8r02jnMqpuCHvGSM3oPjaOfZm92YmU1AD2cjR7KQOghdGHg0TMIhpsItUqE8IIm0TWDTCPgGDGWgaSJzmnIPbogLAXevCiRoIdQOEI4Uo3msvFmRRAahENh3DHDyZQmtXW1XP2pVQhhUXW0Cl3XKCjIRUq48sqVPProUwgBs2ZN5b//+x68ccEBEjGslny6ffIGsANYjiPuHwAfl1JuTtb/RPXJNzflse3JqRjNicu23ZqX0sBcpgTm43elOHRQFxgFPsyq1m5T8WevnIhnloVpQ8aoiQO+lH1oE7se/y1vb2uhvtmNjDQh7Y48JFIIvFMM8qY3kD+qAsMwkVJQfXg8hz8cTcB/Cp/66ufwZPZx0rGv4wqHCX74IS3vrCZaUYF39mx88+bhLZs14FWV5tGjCLe72+IavRGtqCC4aTOhTZuwW1owiotxjS7GVVyMMXo0rs0PIFY7aQkomgWNFRCKSzSWNRYuvhemnQ/A0QP7ePLnP6KuspyswiLmr7gM3V3GttXVtDZFKJmeR0bAxdHmvex5/w+4m4NM9M9mcmA+mUYWla272du8iVZzKzdOWoNLs+G0z8OKHzjXe+ZOrPefoEp7gOfKNiM8owk+tZ5gw3oAMm2d7Jm3UVvlxqfXMPuiLFwFa6hpfATR14lIKSgI/B8zx8xub/qPB346oPvbHbff/kk8no4ExJFmPwErDyE0gmYrQasZzSjE5arlv37xUzL8Xv79C9dTWXmEFStu4P33/4HP5+WGG+7k/PPP5BOfWNV+Lss2CLZm4ff34JI9Xn3yUkpTCPF54AUcZ9tD3Qn8icze3QvQW7rGl0fsENsb3md7w/sU+yYxNWshxb7SHhIo9QNLYh5p7Xa3qzgT/xklCF0wWGnVRs9m2h2/YFrFOqhcR3j6NezfUceGd9/jyK5dRGurieyu4tDOKBWBufinWQQPexAtU1l1xy1Mmj2u94sMZFweD5mnntrnsMK+0LnIdX9xjR2La+xYslacn7xDyXfBaoX3/xeq4o0aASfdAOfdA96OH5iC8RP55H/8jBf/91dsX/0m//rD/fhz81h0yRUUlZ7G7g/rOLB5I1W7/4qwW4kCOxvXsrNxLYW+GUzJXsiZoz5GyGpha+NqZmY9jmv1rx3XzcW/wN76Mkcj3+HQ1CBvyU0setZFsGE9GgIbSYtmYVU+wtwJE9lReTIfPOkCTsOTM4bRJ/8Ob07PD/a26aFpzc1kn+WhMVqHTFO5vXiBB9C8YVqamsnQA9hWuP2atog6OcLjhmBZFsFgGJfLIBgMUVycmFzNtnTSmzG+Z9LuKJJSPgv0nFD8BKamZizNFT583T5RaWjGeGrENOpDBcyufp3iwjIyM9JYDFtAzhVTEHqKP5hjF8LYhXiAaYsCTFvkTLTbtqSmvJndH+1n+4cf0rCnkllnnMp5152JrqfgB22ksfI/IdzkxO4D5E50Vsa2pRHohNuXwcVf/DolM2fz+u9/S3NdLa///gF8WX9l0oLFVO1+FWl3RHRpukHxtOXo/pNYu7+BLYf/QVlWGSUZ53Ekcj5usZnMD17Gd2gVtbXXE3Zl8iXjx9y48SvUNz4CQnDeGeeyt/oQO7ZtItTUxJotmzmz8BlaSu5k894SwvUT2PPCHRTM/B2Fs7fT2W6JBnXqdpxK/d6rMJqijD/DRFgCUwxNumHdiNLigqhsxW0GwWgLnzSRzjIwAMaMGcUdd9zA7Nnn4/V6Oeec01i+/PSEc3W4agTDUQ922CdeT2RsW7B3z0noTZ1TqAo0owTNPR3dNRWhOTOPczKe5awxvwUgWPgZGlqvxaxJfd71zJOL8YwfeP6a/qJpgsLxAQrHz+bUi2b3fsCJjhCw6tcQbXXcM+d8C9zdF/RoY/6KiyiePJWn7v1PGquPEGxsYPMbiYnkxs+ey/KbbydvzFgAGo8G+dMPm3mn/A/kuPMpKzidEk8ZEbOMun0mYPL90vu4ouFWGvc/iW2GWXDBJZRd/xlmAeZPvs+edR9g2zZvHBnPfP3/uObuB3j/P29i21Evlas16ndPZPzSSry5EcINLqo25tOwJYuztmzBE/06Akn0k78moOUS9GUSbE3NupGeEICuR7EsgWXotIuzkNhStm/W1TXyzDOvsXHjc2RnB7jhhjv5y1+e5pprLm4/V5vISzk8a6eUyA8jhw9PJdiaRWbzXgCEPgbdPR3dPQ2hZSb09XktTvbHLDdPNr5Pfglv5hia366k8dUDyHBqFlpofoPslSdIjpjjGU2Hqx/p92HFU6bxqR/9gud+81P2rPugvd0XyOLsT91M2dnLE/pnFfj41P93KX/4/yqpr3qJtyv/Qd1kPx8zllAYmc8ruf9HJDQK+cE2rPBh8saO48xP3IjQdQRw8Zfu4vF7vkXlDqfI9/pKD+u/ejvQURKwtcrH9r+V4h/TStPBTEBQUtOIN9qC8PnIvfpqarJ9uAv9GJmjiJQPTbF7wwhjWW6iHg8iHELHCUV14mocXn/9XSZMKKGgwFnncskly3nvvfVJRb5nS17lkz/mME0XLS3ZvXfs4fgD++ehhVrQLBPDtxRP1rUY3gVdBB7g1I/NxDsmVm1nxQ8guwRhaATOLqH4q4vIPLk4JZ+TnEumoHnVb/9Ixuv3c9nXvs2ZH78RTTeYvew8Pv3z+7oIfBuZOR6u/4/b8PidBHI5e4J8Jf8x7hr/Pe7L3cWijQuwwh+i6QYXfv4ruNwdk9Yut4fLv/4d8kt6XgMjLY2mg37aPsRj65pwjR/PxEcfZdTddzlJ0mwLTdPILkyjqzIOw4gCEqnpSMMgQ6uLJZMS7UlCx40rZs2ajbS2BpFS8sYb7zF9ekd9Zkls5S04vvxuUDVejyEsS+fgwTI+eP9yNqy/gJaWnAGd5+DBMqJRL0ZTA6Chu7ufWS+aEGDmkrFwwX/C1BWw8FMJ+3W/m9wrplJ0xwLcpQP/4fFMyyVjXteKPIqRhxCCk1ddyWfve4QVn/sCvkDP7jlfwM1193wTzchASIsL18ygQrq4cscnidQ7KQpOu/I6Rk2a0uVYr9/Px775/whkdTVekuGNRBm3cBGlf32sIx2y0JyYfMCdhmLvIokHSAiJrjtzANJwIUUEIUDGifWiRXNZtepczjrrGk477QpsW3LjjR2Lxsyohw4JH57J17SGUPaXYzmE0rY1Dh2aysEDc4hGOz5kbncL8xc8j8fTfaRKZ0KhTNZ8sAopdTL2bMZlFeEOfCx5ZwFXfm0Ro0pjX8JIC7h7/rIENx2l/uk9WPU9LKPujA7FX1mMkddzulXFic22d97hmV/8BwBGxnlYkS1Is4Ix02Zyzfd+hNZDtszanRt49DtfI2j1nHeobMwEVvz0VwlrHrZuWMvMkhzI70hHLU2T0LZtPZ5LeLLQfLnYwVp0oxXpKsHQDtLiKaQllsEyKy8fjlRijumqhdGol3DY+b4JaePLbCAUCuDRQ+ie3r9fLS05yNiiskxfI0JPngJBWG4CudN7PR/0P4RSWfK9YNuCQ4emsOaDy9iz++QEgQeIRDLZ9NE5RKN9DzTct28+UuqIaBg9HER3z+y278zTRncIPPQq8AC+2QWM+vJJBM4uiT1e9k7WeROVwCt6ZcbppzN72QoAzNaXkWYFLq+Plbd/uUeBB8ibOo/L5zTiEj3PHy368te6LmqLs+TbmwwD3dPL5zuWfVK4MxAuDyARhoZmOD80Lo8HX04umuZG17o+Ieh6R41iKTSEkEip0Rfb2DTd7QIfG3HvB6UBJfLdIKWg6kgpa9dcyq6dp7X/miejtTWXrVuWYtu9386mpnyqq5yJTcdVY6C5uz7iAngyDE69bHLSfb2huXWyV5Yy6gsL8Uzq2YVjFHgInDl2QNdRnHic8+lbyB09lrZJxKXX30xO8eg+HTt67ulcO3ED2a5g0v2FEydRMG5C1x1CA5kka6mn52R+0owAEqF7kGQiiCIMN3osbXMgv8Apq1lcjNvTdY2Dkxun7cdFxkReYNP7CtZIpNOPRjc+edN0pzWyUol8Epqbc1m39mK2bz+DUKhvoYQNDaPYvm1Jr7/we/ecRNsvutFcj+aahBDJnwIWX1xKRtbgliK5ijIovHUuuddMR/Mne0yW5F41A6Hi0RV9xOXxcuEdd6LpBpNOOpm5yy/o+8ETz6TI28In51Uzad78LrtnnbE0+XFCdLHkATTDQhjdC67hMxFEAIEkw3lteNENA2+mH7fPCT3Vs7IwjCxEEkk0jEhsCI7ACymxuyly04Zl6bFFUAlvoks/29Yxo6mfY4hnBH2zU/coVFVVSmtrTr+PO3p0Inv2JHWLxfaX0NAQiwywLPSWpm4nXPPHZjJnaepqWWYuKKL4zkVknjY64VZlLC7GM2HoYuIVI4PiyVM559O3suKz/96/AyeeCZoL77UPcNnd93D6VZ9oz8kuhMaMbkW+q7sG20Ig0QNd1wgIXcedreHy2QgRt5ZEREB3o7kM/PmJhZCF0HC5u5bj0WN+dCFspBS4rRC96Y0ZdSHsThZfEkveNF3YltHr+QbDCBL51DGnein59sAKQVRWzKS8vKuP3bYF+/ae1L5ttDQghBvNlTwm/cxrpqFpqf2H17wGuaumUHT7fFwlfjS/i5wLJ/V+oEKRhHnnXUhGdk7/DsrMd+L7x52MEILTrryOy+/6Dl5/gPFz5uHP7abmmdAAmSj0sUIrenbid1Xz+XBPntyeP0jEZaIURB1LXjcwXF2fkt3u/C5tumES76pxyzCil0d2K2wgZO91IUzTjZAawkxtXqZ4VEB0J7Soj7ENsynC5kX3Bg5r9f0+x949J+F2Bykq2tfedvjQNILBDovZaG5wVrOKro+aUxcVMXZa+mqZuksCFP3bfMzaEJpPfQQUQ8yMixI2S+efxCd/eC+N1Ue6P0aIDmu+bYI3Jvia24Pm92M3N2Pk52OMGuVM3NoBaKnuJPKRHqtraZobQ3gwZTjumLZQSgm247bpKY28ZbmQlkSjk8h3cvFIKbBtAx0QZvpKWqpveCd89VMRaLjROD88l3/wOo2+/j7wCHZsPx23O0hOzhFM08X+/XM7dkuJ3tyA5lvW5UjDo3P6x6YO7k30ZYSawFWQXl+gQtFXsotGkV3UyyInoccmX2Mi3WbVCx0jPx/y8hKzf7qdxVWaiMasfh0hTNB7tpoNIwszWp3Q1tBwhCuuuAGBpKrqKJpuUFiYA8C8eTN5/vk3KCzM4913/0406kbYkv/6xc/441/+Sl6+83Ty3W99hRUXLW4/p7MSVjj/9RJWOhiUyHfCV+f4yC1psab6WeyWbbhGTyCa079FQlLqbNm8lHnzXqCqahKm2RGeqAeb0WwvmtE1u+KilRPw5w5xoWqF4nhA09uF/ZVX46LOeg6V70qS/svP2R13GR90CmcvKszk9Zf/jhA2P/vRr/AEcvjinR8H4O2313DLLdfyuc99C9vWME03GhGQcOstt/LZ2z4LQKYrAjS1n7NtJayAHlfDDhYl8p3IqJuOaUd5p+ofHAruQQDeQ/sRpkmkoG9hYm1YlptNm5YTjSaKtt7UgO6e3lE0OkbOqAzmn3vilkBUKHokTuTTehmj6xOupttoukCaTlFALS7mccmSRezZU4Vt67S25sbGqIEQnQr7JH7fZRqFPWHsQ3KVISAV90uYXuzaPF479CcOBRMr03iqK/AcOdjNkd0TiWR0WhDhhE4mi6o54+qp6MaI+SdRKFKL0MFOf6phkaS4N4DutpC2iGWSlAkiHf9a2HYsFFPw4MMPce655/LlL3+Z+rrEokAdlnx6xV4pShx6TQkvlf+B2sjhpPvdtUfwVu6lT8vdukFEQhhmBpqR6H+cOLeACWVdZ/YVCkUMzWiPqEknQghEksAYJ28NOAX9Oot8nJTaNiC48ROf4qM31/Diiy9SVFTEN7/z/YTztR3v/D99q6EGJfJCiP8SQmwTQmwUQvxdCJETt+9uIcQuIcR2IcSKQY90CDi4r5Gw3XMOGldDDb7yXbF/yP5jNNWju7vmqFi0cuKAzqdQnDAMkbsGQOtmsZO0RaxkiCTe/RIv8kLaIDQKC0Zh6AaapvGJT3yCdR+uTzxX+zHHtiX/EjBbSjkXp5br3QBCiFnAtUAZcAHw3yJZrOAxRtOhvkWbGM0N+A7sAKv/j45GcwNaJ1fN6CnZiflpFApFV4bIXQMgupmudEQ+mSUv2l7EDECNI9XV7U/9zz33HDNnJH7vEyZe08igRF5K+aKU7Qkl3gXalmiuAh6VUoallHuBXcDJg7lW7wzuVllRQWt130MKjWAzGQd39s+ysEyMcCaanhgDryZbFYo+MJSWvNZNSGPsAV7QIew33fQ1LrroY+zevZuTFi3iT4/9FdC450c/4OTzz+Tc5efyzjvv8MPvfTfhVEM18ZrK6JqbgL/EXo/FEf02ymNtXRBC3ArcCjB+/MDFbrC3q/VIBtj9O4sebMFXvofguMl0KVKZBKO5AaOTFZ9d5KN07uCKPysUJwRahyW//JzdcHQnGB7I6UU3pIQjmxyfflH3GV8TL+UBu2uSeQF8787PYaLTEBPphx76cXtKYRENI0wThOA3P/01GYZOk24ihcQvNKAj/v6YmXgVQrwshNiU5M+quD7fBEzgj/0dgJTyfinlIinlosLCQRSs6GkJWh9oruy9RmYyjJYGvJX7+jQZazQ3dPHHz18+DpHi9AUKxYhE6F3TGvTFCyyEszDK6Pv6k2RhlG2uGgCBM7navk8KkCCstvE50TVSyriSIclDKNP97e/VkpdSntvTfiHEjcDFwHLZUYGkAohf6VMSa0sbcrAif2hgIg/gaqxF6gbh4h4sCmnjDgUQmR0pi71+FzNO61/svUJxwtI5uiY+xUFveALtueX7dClXJkQS26TdUYhb6+STBw0DzSnyHdsWWEjsmLgnTtQ6Md9DUzFqsNE1FwBfAy6VUsaHpTwJXCuE8AghSoGpwPuDuVY6saOC1qrBLfF311XhPnqo2/16azO6kWjFzz5rLIb7mJ+PViiODTr75G3LEf6+4AmA3ndLXmhdUx9IKdoXQcWvUm1zuxhx7l5nGZRAIjsWRMXH0ku9PWpy2C35Xvg14AFeiq3efFdK+Tkp5WYhxGPAFhw3zu1SDkGA6wBpqfIh++mPT4anugKpG0Rzu7qdjOZGdHdHThqp2SlNJaxQjHQkGkJaHa5R2Q9Lvh+uGojFypsg4xTSWQjV4TFoUww75qrR41MLCw0kMXdNm1tGtEfDC6kls+97ZSDlWgcl8lLK5CWNnH0/AH4wmPP3B20Qv4cD9ccnw3N4P1I3MLMSI2hcoWyEp+OD5pkZGnRBEIXiRMHr9VJT30C+lAjb6vCb9Ccyu79CbzuWeDt2J0GOCa60NXQJMpZaONEJ0+GTT8hCKfV2S9+x8HsXbyklNTU1eL39K9M5YnLX+OXAs7g1H+pbFfm+IABv5R6C+lSsTCf2XQsHcYtp7X0kNtPPHsQks0JxglFSUkJ5eTnVh6ugbosj8o1VUKv1mlVyoISbjiBdHeJrhjUME7yxIiItWiO6O4plGdiROjRLYsmoEwOiN6PZEYSwkIZBFAufdCN9dYCTWjhsurCReKUbTUrcLVW9jsnr9VJS0j8PwIgR+YFa8sIGESwEuoZLDRQhJb7yXbSOn47ty8Robk4oDnIgZyvXz7ghZddTKEY6LpeL0tJS+Oty+PhjjlX+16vh39dDXvLCO4Nl/U9vonZBh/AefHMUk4+0sCBWJ+Kf2afjn7eDgwfLcG2YS/RoNTXhCrJaI0TG3sWoI+/j9pVjTjqJjcZ+LgmfRP3yr4BmkbtvBe/tm0St1syKyDwyW0LM/MV1aXkfJ3zumuzcxXz8P35DUenACmZ3h7BtfAd3IsIh3KEshOj4PS2ftJEst1rhqlD0G282BOsg1OBs+9JXXMcjE6tUtVb58GgdK27dsXDJaMRLwA4TtJw0wh7TBGmjWxEwQ7hiCQo1NDTb8Tjo0QBGrBi4T6Y3tfiIEXk5wCiV3NzTCOQXcO33/pPJM+ekdEyaZZJxYAcuu8NVExYHyZ/cP5+aQqGI4c1xRD5Y70xuerPTdimPUdz+Wm+B4FEfHr1D5D2m47aJRL3k0EzQbAbAbVroVhjNjqBFgrhiDhMNgTQdndIjfoyYjz5DpndubsSIPAMW+VMApwL9Zd/9IWXerEFlmeyMbvvQjY74+e35rzE1N/2VnxSKEYk32xH4UIPzWqQvANHr66h/7D/kaIJb6wgS9FmOyEcjXnJFOTKW88BtWhhmEN2KokWDuGmz5AWm6ZynzZIXUuBDiXyfaM4I9d6pE5rmJitrQULb0q9/i3kHq7s5ov8YvrPai4O4w3W8X7qO6Xlds1AqFIo+4MvpcNd4c9J6KW/WRIhpetZh50W8JZ9hOqul9KigwN1RWcplOSLfbsnLNkteIxqNFR+PBDDQ8OFGINBSaFh2ZsSIfEtG31eztZEVmIfeaYGEd9YsZi45m0ll8wY9Js0oSYiNz2h4HdOQTMud1sNRCoWiW+J98mn0xwPogRz0OshoMfHj5JfyxFnyGTKMbWtkRsIUeytwx7LSdljyEbTWBrRY/HzIbMY2HYPPFcrHkHq7q0YokU8POTFXTWfG/PA/WP5vX8DlGYTvXEryIx2ZHXQzSLXnbXyGjwlZEwZ+XoXiRMabA6F6548vJ62X0rMC6DWCoqNRXOf/fwAJE68+Qpimm+xoK0JAcTjmk7dsDCuIbkdwmRFawjUANEfrHJEPexyRRycjNukqjtWiIcc7uTnJsx8Ll4usgiJO/di1Az73uNomxh3Z3r495tBqdo8JMzl7MlofMlYqFIokDKUln5WFUSsonvBJXBMdg9Cd6W/f7yNIKOQnF6es3xjbia5xLPkQmhXBbVo0hJxKc82RWmxTQ1Q763LiRV65a9KAEG6ys0/qsc9JF11Gfkn/0x8LaTPtcC2FR9ejmyGEtCgpf43dxUL54xWKweDLiZt4zUnrpbSsLHzbPGQu/SFGbLW6p7Tj6d9HiGAwQC5OOOdY0YiQElebu8aO4jYtaoPlADRHa5CWQFY4LhtDah3umnS+jzSe+5gmK2sOut6zO0Y3DM69+d/6fe7q3Ho8poVuRyk8up7C6g8Rdi3lBajIGoViMAypTz6Ab3cWaDoujxchNNzTlnYMhTChOJHP9ETJaQ057hoziGZF8EVMZm5eD0BTuAbb1NAPmIQty7HkSW+MPIwgkff0I8McdO+q6UzJrNnMOnNZn88b9MNrCxqJxCI6i4+8z/iDr7BvFEhNMD1XWfIKxYCJj5NPs09ey85Gi83LuTwe3D4fYnKHFmiAFfS2i7zusShsbMVlWugxn7wuJfnNzmr6cOthbEvgOWDTGgkmuGvS+j7SfoUhYnbB7H71z8k9tc99z/7UzXjjfHE94VleRotPsGGS8wCWV7edrKYD7C52tqflqcgahWLA+HJiE6/pt+Q1txs9EHBe6zoZ2dmQPxmy4wIqgjoBnAlX3WMzpr4ZDXDFomvA8bfrlkWguQU7opFRbhIOt2CgpX0hFIwgke+PV0sIFznZC/vcPyM7hyXXXt9rv/Gz53LKWRcB8O6MxPHsGS0YnTlapTNQKAaDNwesCDQfSbtPHsAo6CjNmZkbS3NQenZ7W06wtV1EDY9NRsSJvtFj7po2XJEIua0hwlVu3GGLSLAhFkKpLPm0kBWYja73L73wvHMvoHhy9/50oWksvf4WFo1ahKEZrJkq2l02ALtGK1eNQjFo2oTdjqbdkgcwCjtE3p+b77yYtLS9baxV0/5a99odx5khdLtD5N2RCDmtYcwKx3K3WqpwYaR9tSucoCLfXXx8TwhN49zP3I7oJvxxzrLzKZxQSoYrgzkFcwh6Olw2rR44lKcmXRWKQROfqybNPnkAPd6Sz4n9qEw6uz2P/RgOd/T1xIt8h7sGwBMOk90aRq+K6UdjBQHpTXsRb0iRyAshviKEkEKIgti2EEL8UgixSwixUQjRd9/IAOktUiaevk66dmbUpCnMO39ll3ZPRiZLrv1U+/Ypo50fkTaXzZ5RAoQKn1QoBo1uOEW5YWgs+YKOug/+NneNvwhOvwOAYjpSoBjxIm85aQ3acEciuGyb3BYn/YreeAC/OTSZ3gct8kKIccD5wIG45pU4dV2nArcC/zPY6/TGlMlfZ3bZL/H7Z/bYTwiD7OxFA77OGddeT0Z2TkLbKVdcQ0ZWh4VxSrEj8m0um92xWt0qnYFCkQLarPmh8MkXdoh8u08eYNk3oHAmLjrSHAgNNJcj9J0teXfEeT2q0Ym08YTrES0drp50kgpL/uc4xbzjl2ytAn4vHd4FcoQQo1NwrW4RQmPUqIs45eSnmTv3frKykueeCQTKMIyBV4LyZGSy9FM3t2/nFI9m4cpLEvrMK5qHz/AR9Ag2lgr2jBYqnYFCkSq8OU5RbnfqynZ2R/zEqz9e5A0PXP4/XQqJ6znOD5A70pyQqsAddkQ+O+j87QnXYzcfSdewExiUyAshVgEVUsoNnXaNBQ7GbZfH2pKd41YhxBohxJrq6tRkfywsWM7iRU8wf/4j5HRyzeTm9N8f35mZZy5jXNlcAM7+5M3oRmLpQZfmYuEox0O1eoZgd7FQ6QwUilThzR4SfzyAUZDf/jrBkgcYswBW/TfkxQoOnXo7xljntSbNhK7eTq533Y4eOyIvhHhZCLEpyZ9VwDeAbw9mAFLK+6WUi6SUiwoLU1v3ND/vDE5a+GcWLnyUvLwzAbqI/kBZfvNtlM4/iSmLk8fbn1rstL8/XVCVq/zxCkXK8OUMiT8eEt01/s4iDzDvGrhjLXzmVbjgP9DzkvQBCleuRHQqwC1bhkbke/X8SynPTdYuhJgDlAIbYvnSS4B1QoiTgQpgXFz3kljbsJCbs5jc+Q/T2LiRjIzUlPnLHzuOS7/yzW73t02+ht2xRVDKH69QpAZv9pD44wH07I65NrevG/eQEFDi5MHS87r++HjLyshbvJiMk06i5e2329vtpsNd+qaDAfsPpJQfSSmLpJQTpZQTcVwyC6WUh4EngetjUTanAg1SykOpGfLAycqaOyh/fGcMd/cxrjPyZpDjyWnfViKvUKQIb86QWfL9xUhi7Rd8/nZ8Ph+ZZ5yR0G41HETaVpf+qSZdTuJngT3ALuC3QP+zfB3nCCFYXLy4fVulM1AoUsQQ+uT7S2d3jWfmTALLljkiv+T0xM5WBLsx/Q6OlIl8zKI/GnstpZS3SyknSynnSCnXpOo6xxOnjnb88iqdgUKRQobQJ99f9NychO28G5x0KD6fD++0aRijRiXst2p3k25UuEcaaRN5lc5AoUgh3pwh88n3FyPOktcLC8i+8ELAEXmAzCVLEvpbtXvSPiYl8mlkfNZ4RmeOVukMFIpU4s0+hi35DpHP+/jHEbF5O3fs7+xLLk7oPxQiPzTrak9gThl9igqfVChSiS8HIs3DPYqkGIUFeKZNwygsJOfaruVDM087Dc+MGYS3bQNAtlZjhxvTOiZlyaeZU0afoiJrFIpUcgxb8q7iYiY9+U/GP/gARm7yMeZd7/jphctZRGnVpNcvr0Q+zZw+5nSVzkChSCXHsE++L2RffBH5t3yGKa+9imfqFKwjH5HOKq9K5NNMnjdPpTNQKFLJMWzJ9wXhdlP0la9gFBRQ9LWvYR7egJR27wcOEOWTVygUxxceP9gFvfc7DvCfeSbh88/Dqi5P2zWUialQKI4/jtHFUANBfOGr/H6WctcoFArFiMRj6GzMTl96AyXyCoVCMYx4DA3oe2W7/qJEXqFQKIYRj0tDk0rkFQqFYkTiMXQ0fGk7vxJ5hUKhGEY8hrLkFQqFYsSiRF6hUChGMIau4RbpK0quRF6hUCiGGZ/rGBZ5IcQdQohtQojNQogfx7XfLYTYJYTYLoRYMdjrKBQKxUjFZ6RP5AeV1kAIsQxYBcyTUoaFEEWx9lnAtUAZMAZ4WQgxTUqZ/oKGCoVCcZzhxMqnh8Ge+TbgR1LKMICUsirWvgp4VEoZllLuxan1evIgr6VQKBQjkmNZ5KcBZwoh3hNCvCGEaKtcPRY4GNevPNbWBSHErUKINUKINdXV1YMcjkKhUBx/eAw9befu1V0jhHgZKE6y65ux4/OAU4HFwGNCiEn9GYCU8n7gfoBFixbJ/hyrUCgUIwGPK32WfK8iL6U8t7t9QojbgCeklBJ4XwhhAwVABTAurmtJrE2hUCgUnTiW3TX/AJYBCCGmAW7gKPAkcK0QwiOEKAWmAu8P8loKhUIxIhlWd00vPAQ8JITYBESAG2JW/WYhxGPAFsAEbleRNQqFQpEcdxot+UGJvJQyAnyym30/AH4wmPMrFArFicCx7K5RKBQKxSBRIq9QKBQjGI8rfT55JfIKhUIxzChLXqFQKEYw6YyuUSKvUCgUw4yy5BUKhWIEk84Vr0rkFQqFYphR7hqFQqEYwSh3jUKhUIxglLtGoVAoRjDKXaNQKBQjGOWuUSgUihGMEnmFQqEYwai0BgqFQjGCUZa8QqFQjGBcuhJ5hUKhUAyAQYm8EGK+EOJdIcR6IcQaIcTJsXYhhPilEGKXEGKjEGJhaoarUCgUiv4wWEv+x8D3pJTzgW/HtgFW4tR1nQrcCvzPIK+jUCgUigEwWJGXQFbsdTZQGXu9Cvi9dHgXyBFCjB7ktRQKhULRTwZbyPuLwAtCiJ/g/GCcHmsfCxyM61ceazvU+QRCiFtxrH3Gjx8/yOEoFAqFIp5eRV4I8TJQnGTXN4HlwJeklH8TQlwNPAic258BSCnvB+4HWLRokezPsQqFQqHomV5FXkrZrWgLIX4PfCG2+VfggdjrCmBcXNeSWJtCoVAohpDB+uQrgbNjr88BdsZePwlcH4uyORVokFJ2cdUoFAqFIr0M1id/C/ALIYQBhIj51oFngQuBXUAr8OlBXkehUCgUA0BIeey4wYUQ1cD+uKZsoKEfp+itf0/7k+3rra3z/vjtAuBoL+PtD6m+Fz316Wt7T9ud9x3r9yPVn43O2+qzcfx+V46Hz8YEKWVh0j1SymP2D3B/Kvv3tD/Zvt7aOu/vtG/NsXwveurT1/Ze3n/nfcf0/Uj1Z6On+3Gs34t0fzaOt/txPH82pJTHfFqDp1Lcv6f9yfb11tZ5f3/H2x9SfS966tPX9p6203kvBnL+of5sdN5Wn43+j2GgqM9GHMeUu2YkIYRYI6VcNNzjOFZQ96MDdS8SUfejg3Tci2Pdkj+euX+4B3CMoe5HB+peJKLuRwcpvxfKklcoFIoRjLLkFQqFYgSjRF6hUChGMErkFQqFYgSjRH6YEEJkxgqtXDzcYxlOhBAzhRD3CSEeF0LcNtzjGW6EEJcJIX4rhPiLEOL84R7PcCKEmCSEeFAI8fhwj2W4iOnEI7HPxCcGcg4l8v1ECPGQEKJKCLGpU/sFQojtsWpYd/XhVF8HHkvPKIeGVNwLKeVWKeXngKuBJekcb7pJ0f34h5TyFuBzwDXpHG86SdG92COlvDm9Ix16+nlvrgAej30mLh3Q9VR0Tf8QQpwFNOMURZkda9OBHcB5OLnzPwCuA3Tgh51OcRMwD8gHvMBRKeXTQzP61JKKeyGlrBJCXArcBvyflPJPQzX+VJOq+xE77qfAH6WU64Zo+CklxfficSnllUM19nTTz3uzCnhOSrleCPEnKeXH+3u9wSYoO+GQUv5LCDGxU/PJwC4p5R4AIcSjwCop5Q+BLu4YIcRSIBOYBQSFEM9KKe10jjsdpOJexM7zJPCkEOIZ4LgV+RR9NgTwI5wv9nEp8JC6z8ZIpD/3BkfwS4D1DNDzokQ+NSSrhHVKd52llN8EEELciGPJH3cC3wP9uhexH7wrAA9O9tKRRr/uB3AHTuGdbCHEFCnlfekc3BDT389GPvADYIEQ4u7Yj8FIpbt780vg10KIixhg+gMl8sOIlPLh4R7DcCOlfB14fZiHccwgpfwlzhf7hEdKWYMzN3HCIqVsYZCp2tXEa2pQlbA6UPciEXU/OlD3onvSdm+UyKeGD4CpQohSIYQbuBanOtaJiLoXiaj70YG6F92TtnujRL6fCCH+DKwGpgshyoUQN0spTeDzwAvAVuAxKeXm4RznUKDuRSLqfnSg7kX3DPW9USGUCoVCMYJRlrxCoVCMYJTIKxQKxQhGibxCoVCMYJTIKxQKxQhGibxCoVCMYJTIKxQKxQhGibxCoVCMYJTIKxQKxQhGibxCoVCMYP5/Hjr8jgu3hmYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i in range(0,9):\n",
    "    plt.fill_between(x, flavour_cv[i]-flavour_std[i], flavour_cv[i]+flavour_std[i], label=flavour_labels[i])\n",
    "plt.xscale('log')\n",
    "plt.legend()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}