Results_on_synthetic_linear_graphs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "bdKrXxMHM34Z"
   },
   "source": [
    "Code for generating the results for synthetic linear Gaussian causal graphs (Fig. 10)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "8Xqc_6q-WNW1"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import networkx as nx\n",
    "import numpy\n",
    "import pandas\n",
    "import networkx\n",
    "from itertools import combinations, permutations\n",
    "import logging\n",
    "from causallearn.utils.cit import CIT\n",
    "from sklearn import linear_model\n",
    "import collections\n",
    "import matplotlib.pyplot as plt\n",
    "import pickle as pkl\n",
    "from datetime import datetime\n",
    "import copy\n",
    "\n",
    "from causal_discovery.utils import *\n",
    "from causal_discovery.pc_alg import PCAlgorithm\n",
    "from causal_discovery.mb_by_mb import MBbyMBAlgorithm\n",
    "from causal_discovery.sd_alg import SequentialDiscoveryAlgorithm\n",
    "from causal_discovery.ldecc import LDECCAlgorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "pnj3tN45cbZq"
   },
   "outputs": [],
   "source": [
    "def get_params(num_conf, num_med, p_link_conf, p_link_conf_med, \n",
    "               p_link_med, p_link_x, p_link_y):\n",
    "\n",
    "  def generate_links(max_links, prob_link, atleast_one=True):\n",
    "    while True:\n",
    "      links = np.random.binomial(n=1, p=prob_link, size=(max_links,))\n",
    "      if not atleast_one or np.sum(links) > 0:\n",
    "        return links\n",
    "  \n",
    "  # covariates + treatment + outcome\n",
    "  A = np.zeros(shape=(num_conf + num_med + 2, num_conf + num_med + 2))\n",
    "\n",
    "  # generate edges between the covariates\n",
    "  for i in range(1, num_conf):\n",
    "    A[i, :i] = generate_links(i, p_link_conf)\n",
    "  \n",
    "  # connect some covariates to the treatment.\n",
    "  A[num_conf, :num_conf] = generate_links(num_conf, p_link_x)\n",
    "\n",
    "  # connect the covariates to some mediators.\n",
    "  for i in range(num_conf + 1, num_conf + 1 + num_med):\n",
    "    A[i, :num_conf] = generate_links(num_conf, p_link_conf_med, atleast_one=False)\n",
    "\n",
    "  # connect the treatment to some mediators\n",
    "  A[num_conf+1, num_conf] = 1  # connect treatment to first mediator.\n",
    "  if num_med > 1:\n",
    "    A[num_conf+2:num_conf+1+num_med, num_conf] = generate_links(num_med-1, p_link_conf_med)\n",
    "  \n",
    "  # generate edges between the mediators\n",
    "  for i in range(num_conf + 2, num_conf + 1 + num_med):\n",
    "    A[i, num_conf + 1:i] = generate_links(i - (num_conf + 1), p_link_med)\n",
    "\n",
    "  # add edges between covariates and outcome.\n",
    "  A[num_conf+num_med+1, :num_conf] = generate_links(num_conf, p_link_y)\n",
    "\n",
    "  # add edges between mediators and outcome.\n",
    "  A[num_conf+num_med+1, num_conf+num_med] = 1 \n",
    "  if num_med > 1: \n",
    "    A[num_conf+num_med+1, num_conf+1:num_conf+num_med] = generate_links(num_med - 1, p_link_y)\n",
    "\n",
    "  # disconnect treatment from outcome.\n",
    "  # A[-1, num_conf] = 0\n",
    "\n",
    "  # A *= np.random.uniform(low=.1, high=1, size=A.shape)\n",
    "  A *= np.random.uniform(low=.25, high=1, size=A.shape) * (2 * np.random.binomial(1, 0.5) - 1)\n",
    "\n",
    "  labels = [\"W%d\" % d for d in range(1, num_conf + 1)]\n",
    "  labels.append(\"X\")\n",
    "  labels += [\"M%d\" % d for d in range(1, num_med + 1)]\n",
    "  labels.append(\"Y\")\n",
    "\n",
    "  # noise_vars = np.random.uniform(low=0.3, high=1, size=(num_conf+num_med+2))\n",
    "  noise_vars = np.random.uniform(low=0.1, high=0.2, size=(num_conf+num_med+2))\n",
    "\n",
    "  B = np.linalg.inv(np.eye(A.shape[0]) - A)\n",
    "  noise_cov = np.diag(noise_vars)\n",
    "  cov_matrix = B @ noise_cov @ B.T\n",
    "\n",
    "  return A, noise_vars, labels, cov_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "hWBTdTWgNYAf"
   },
   "source": [
    "Below we demonstrate the various (local) causal discovery algorithms on a randomly generated graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "0qSr_-irdHwp"
   },
   "outputs": [],
   "source": [
    "NUM_COVARIATES = 10\n",
    "NUM_MEDIATORS = 3\n",
    "A, noise_vars, labels, cov_matrix = get_params(num_conf=NUM_COVARIATES,\n",
    "                                                       num_med=NUM_MEDIATORS,\n",
    "                                                       p_link_conf=0.1,\n",
    "                                                       p_link_conf_med=0.05,\n",
    "                                                       p_link_med=0.1,\n",
    "                                                       p_link_x=0.1,\n",
    "                                                       p_link_y=0.1)\n",
    "graph_true = adj_matrix_to_nx_graph(A, labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 516
    },
    "id": "NhhHfF9PdQ1w",
    "outputId": "f1a5fae8-fd66-4507-e9d9-9f76b7d1d8f2"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApQAAAHzCAYAAACe1o1DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd1iTZxcH4F8GCg6U4d6jKiJKcSKoqKgo7j3qwi1Wba2rVuuqrbOfrVpwb+vWulBEtAIuEMWtgAucCMhGkjzfH08BkRXIeJNw7uvKJWS87yFgcvKMc0SMMQZCCCGEEEIKSSx0AIQQQgghRL9RQkkIIYQQQlRCCSUhhBBCCFEJJZSEEEIIIUQllFASQgghhBCVUEJJCCGEEEJUQgklIYQQQghRCSWUhBBCCCFEJZRQEkIIIYQQlVBCSQghhBBCVEIJJSGEEEIIUQkllIQQQgghRCWUUBJCCCGEEJVQQkkIIYQQQlRCCSUhhBBCCFEJJZSEEEIIIUQllFASQgghhBCVUEJJCCGEEEJUQgklIYQQQghRCSWUhBBCCCFEJZRQEkIIIYQQlVBCSQghhBBCVEIJJSGEEEIIUQkllIQQQgghRCWUUBJCCCGEEJVQQkkIIYQQQlRCCSUhhBCiIsaEjoAQYUmFDoAQQgjRBzIZcOcOEBgIBAUB164BoaFASgq/zcgIMDYGrKyAli2Bpk2BZs3492IaviEGTsQYfa4ihBBCchMRAWzcCPz1FxAVxa+TSnkSmRsjIyAtjX9dvTowZQowejRgaan5eAkRAiWUhBBCSA4CA4GlS4ETJwCRCJDLC38ssRiQSIBBg4CffgLq11dfnIToAkooCSGEkM+kpAALFwIrVvAkMK+RyIKSSnlyumwZ8N13/PiEGAJKKAkhhJD/XL8OfPMNEBYGKBSaPVfz5sDOnUCDBpo9DyHaQMuECSGEEABbtwKtWgHh4ZpPJgEgOBiwtQVOndL8uQjRNEooCSGEFHlr1wJjxvDyP6qslSwImQz49Ano1Qs4eFA75yREU2jKmxBCSJG2aRMwfrxw5xeJ+Kad48cBV1fh4iBEFZRQEkIIKbK8vYEuXYQvTC4S8VJDgYGAjY2wsRBSGJRQEkIIKZI+fuRFx9++1c6ayfxIpYC1NXDjBk8uCdEntIaSEEJIkfTDD8C7d7qRTAJ8TWVICLB8udCREFJwNEJJCCGkyPH2Bjp3FjqKnEmlwM2bNPVN9AuNUBJCCClSGANmztTd/tqM8W46hOgTGqEkhBBSpFy/DrRsKXQUeROJgBcvgKpVhY6EEOXo6OczQgghRDPWrePTyrpMLAY8PYWOghDl0QglIYSQIiMqCqhcGUhLEzqS/FlYAK9eAcWKCR0JIfmjEUpCCCFFhpdXwZLJkSP5mkbGgHbtcr7Pkyf8dl9f/n3FisCSJUBAAPD+PS9PFBgIjBtXsHWbHz4AV64of39ChEQJJSGEkCIjMLBwNR7j4nhrxi+1awfUrctvT9e0KTBiBODjw//t1w+4dAn46y/elUdZYjEQFFTwWAkRgo6vIiGEEELU5/r1wk13798PDBsGuLsD8fGZ148Zw0ciTU0zr/P3B+rU4XUl050/z6eup0wBfv4ZiIjI/5wiESWURH/QCCUhhJAiQS4HgoML99h9+/i/Q4ZkXmdqykcft27Net/Y2KzJZLrr1/m/yu7clsuBq1cLHishQqCEkhBCSJHw4gWQklK4x8bFAYcOAW5umdcNGcK77Ozfr9wxOnTgo6OPHyt/3qdP9WMDESGUUBJCCCkSEhJUe/zWrbx+ZcOG/Hs3N+DgQeWO26kTMHw48OefQHS08udkDEhKKly8hGgTJZSEEEKKhMKOTqa7dAkIDeWJZKNGQIsW2ae7c/L118CBA3z6eu7cgp9X1bgJ0QbalEMIIaRIKMzu7i9t2wZMnQoYGwOPHgF+fnnf39aW9w1/8gTo1g349Kng56Q6lEQf0AglIYSQIsHERPVjbN8OWFoCEyfy5DIvtrZ8d/fz50DnzllLCxWEsXHhHkeINtEIJSGEkCKhalUGkQhgTFToY7x6BaxcCTRoAOzYkfv9mjThyWREBF8/GRtbuPNZWFBCSfQDJZSEEEIMTlxcHO7cuYOQkJCMy507d8DYNQBWKh07v3WQ9erxZBIA5s0Dvvoq6+1hYbwFZH5EIr5OU1T4/JcQraGEkhBCiN6SyWQIDQ3NkjSGhITg2bNnAACpVIoGDRqgcePG6NmzJ86fL42LFxlkMs1lafb2fFocAE6ezH77qFF5j26mk0iA5s3VGhohGiNijDGhgyCEEELy8/79+2yjjvfu3UPKf9ugK1WqhMaNG2e5NGjQAMU+29Wybh3fVKMv73zHjwM9ewodBSH5o4SSEKL7UlOBO3d4I+Zbt3ghv5QUPhdobAxUrAjY2QHNmvHFbRKJ0BETFaSmpuLhw4dZEseQkBC8efMGAGBsbIxGjRplSRxtbGxgmT4smIf79wFra03/BOohlQJv3vB1lIToOkooCSG6KTYW2LmTb6W9e5f3shOJ+LusTJY5xCQW8wQyvZ2IsTGfJxw/HhgwACheXLifgeSJMYbIyMhsieOjR48g+693Ya1atTISxvTksW7dupCo8KHB0RG4coV3udFVUikwaBCwe7fQkRCiHEooCSG65fZtYP16YNcuPjIJFHx+Uizm2ULZsry+y8SJQI0a6o+VKC0xMRF3797NttYxJiYGAGBqapolaWzcuDEaNWoEU1NTtceyfz8weLDaD6t2AQF8PSYh+oASSkKIboiPB2bOBDw9M0ch1UEi4SObixfz40tpL6ImKRQKhIeHZ0kaQ0JCEBYWBsYYxGIx6tWrl22tY/Xq1SHS0nbmT5+AypWBDx+0croCE4t5e8eQENrhTfQHJZSEEOH5+AAjRwKvX2tuHlIk4pWmd+3Sn0V0Oi4mJibbJpm7d+8iMTERAGBpaYkmTZpkSRytrKxgoo4K4ypauRKYPVt3N+f8/Tef8iZEX1BCSQgRjkIB/PAD8PvvmdPUmpQ+Ovnnn3wanCglLS0Njx8/zrbWMSIiAgBQrFgxNGzYMMsGmcaNG6NChQpaG3UsKJmM13gMCQHkcqGjySSVAl278t3dOvrUEZIjSigJIcJISwNGjwb27BHm/EuX8qrTJANjDG/fvs1WDPz+/fv49F8T6mrVqmVb61ivXj0YqaNRtpbdu8cHrdW1ukJVIhFgago8fMgLFxCiTyihJIRon0IBjBgB7N0r7Jzjr78Cc+YId34BJScn4/79+9nWOr5//x4AUKJEiWyJo42NDczMzASOXL1+/ZV/rtCVd8I9e4ChQ4WOgpCCo4SSEKJ9330H/O9/QkfBbdkCuLkJHYXGMMbw/PnzbGsdHz9+DIVCAZFIhDp16mTbJFOrVi2IxWKhw9c4uZyvVTx6VPgyQtOnA2vW0FQ30U+UUBJCtOvcOaBLF6GjyGRszOtc1qkjdCQqi4uLy1KaJ330MS4uDgBgZmaWbZ2jtbU1SpUqJXDkwkpNBXr1Ary9hUsqR48GNm/mS4kJ0UeUUBJCtCcujneyeftW+OGgdFIp0LIl8O+/BX43j4iIgFgsRuXKlTUUXM7kcnmW/tXply/7V385ZV2lShWd3SQjtNRUYMgQ4Ngx7U9/T5kCrF1LySTRb5RQEkK0Z9w43vlGl7bVpvvjD+Dbb5W++86dOzFhwgS0bt0aPj4+GgsrKioq24jj3bt3M/pXV6xYMcf+1cWpQ1CByeXAkiV8v1b695oilfLLqlXA5Mk0zU30HyWUhBDtCAzkLRF1lbExEBkJmJvnebfExES4u7tjx44dAIAyZcogJiZG5ZG/9P7VX651fP369X/hGcPa2jrbJply5cqpdF6SXXAw8M03wIMHmhutdHAAduwwiJUWhACghJIQoi2jR/PGxLpSo+VLIhEfLvr++1zvcu/ePfTt2xehoaFQfDZlHxkZqfS0N2MMr169yjZd/fDhw4z+1TVr1sy21rFu3bqQUpcfrXn6NBJz5sTi1ClrJCXxPw9VVmlIJHzE09wcWLgQcHenKW5iWCihJIRo3ocPvNfdf7UMdVaNGkB4eLZ3esYYtm3bhsmTJ0Mmk0H+xVyol5cXuuSw0SgxMRH37t3Lljym968uXbp0lqQxvX91mTJlNPczknydPXsWvXv3xqdPnxAXJ8fevXyN4717fJpaoVAuuRSLeSIqlwOtWjHY2fljxYpWKFmSPhgQw0MJJSFE81avBmbN0p2NOHnx8sqyCz0hIQHjx4/Hvn37cry7RCLBsmXL0K9fv2xrHUNDQzP6V3/11VfZ1jrWqFGDNsnokOTkZMyaNQvr1q0DAJiYmCApKQkAn/q+epUXKQgMBK5dA/4r2ZmjypWBVq2AZs0AV1fgzZtz6NKlC1q2bIkLFy6gRIkS2viRCNEaSigJIZpnbQ3cv6/6cU6cABwdARsb4L+2fxnMzPgQUng40KZN4Ra/SaXAgAG84DqAkJAQ9OnTB8+fP882Kvk5iUSScbuFhUW2/tUNGzbUif7VJHfBwcEYNGgQwsLCMpYzlC9fHm/fvs31Ma9fA0+eAMnJfPC9eHGgRAleyMDSMut9T58+DVdXVwCAnZ0dTp8+jQoVKmjs5yFE22jcnRCiWQkJvJecOowdy2tGbt4MuLhkvW3dOqB0aWDkyMLvpJDJAH9/MMbg4eGBqVOnQqFQZFkvmRMLCwvs3LkTjRs3RsWKFWnUUY/I5XKsXLkS8+fPB2Msy+/a2Ng4z8dWqsQvykhf5gAAt2/fRtOmTXHu3Dk0bNiwUHETomtoSTAhRLNu3VLfVPfbt7zGSpcuwPjxmdf37s371c2cCYSFqXaOFy8wtn//jPWS+SWTABAdHY0OHTqgUqVKlEzqkWfPnqFt27aYO3dujmtj80soCyI6Ojqj85BcLsebN28ypr8JMQSUUBJCNCsoSL3bWQ8eBPbt4zuya9Tg22Y9PPjiNg8PtZxizNdfw87OLmOaOr/d1TKZDI8ePVLLuYnmMcawa9cuWFtb49q1a7neT53LFKKjoyGRSDK+l8vlSEpKQufOnTNKUBGizyihJIRoVnCw+uujuLsD8fHA1q3Ahg1AsWLq68ctFqO1iQmCgoIQFxeHS5cu4fvvv8dXX30FABCJRDn2uA4JCVHP+YlGxcTEYODAgRgxYgSSkpLyXBurzoTy8ynvdAqFAnK5HKNGjcqYcidEX1FCSQjRrOho9deejIkBxowBOnQABg3iHW4iI9VzbImExww+Mtm2bVssX74cjx8/Rnh4OP744w907NgRRkZGGfcBKKHUF5s2bcKhQ4eUuq86d2JHR0fnuXxi6dKlGDZsGFJTU9V2TkK0iRJKQohmJSZq5rheXsCVK8Djx8CePeo9dnJyjlfXqlULU6ZMwblz5xATE4Pjx49j1KhRqFSpEpWB0RNTpkzB0qVLUaJEiSxT0DlR5wjlhw8f8hwNBYD9+/ejY8eOiP7vAw0h+oQSSkKIZmlyk0pqqmaKpSsRc8mSJdGzZ09s2rQJr169woIFC9QfB1G7EiVKYN68eXj27Bnc3d0hlUpzXMIgEonUmlBGRUXlex/GGPz9/bF9+3a1nZcQbaGEkhCiWSVLCh1BwVHNSINXrlw5rF27FleuXMlIKD8fsRSLxWrd5Z1XQpm+bKJRo0bYvHkzJk+erLbzEqItlFASQjSrXDleMFxfyOU8ZlIk7NmzByYmJvDx8YGTkxMAnlgyxtSaUOa0KQfgieugQYMQEBCA27dvY8yYMWo9LyHaQgklIUSzmjblSZq+UCh4zMTgPX36FOvXr8fs2bPRoUMHnD9/Ht7e3rC2toZCoUBJNY2uM8YQHx8PIHMUtHLlyrCzs0PlypWxY8cO2NvbUw1TotcooSSEaFbTpoXvXCMEkQiwtRU6CqIF8+fPh4WFBaZPn55xnbOzM4KDg3Hq1Cl8//33ajlPamoqZP9VOnBycsLx48fx4sULeHp6IiIiAqdPn1bLeQgREvXyJoRoVkoKUKqU/oxS1q3LGzQTgxYcHAw7Ozt4eHhgwoQJGj/fyZMnUa9ePdSrVy/L9S1atIC5uTm8vLw0HgMhmkQJJSFE85o25QXOdf3lRioFhg/nBdOJQXNxccHTp09x7969fDshadKOHTswatQoPH78OKN4PiH6iKa8CSGa93nfbV0mk/GC6cSg+fj44OzZs/j1118FTSYBYNCgQbCwsICHmtqGEiIUGqEkhGheQgJQoQKQlCR0JLkTiQArK+DuXc3WziSCUigUaNGiBYyMjBAQEKATG2Fmz56NjRs3IjIykgrkE71FI5SEEM0rVQoYPVr3ywdNnUrJpIE7ePAggoKCsHz5cp1IJgFg4sSJ+PjxI/bt2yd0KIQUGo1QEkK048EDwNpaN9dRikSAqSkQEcGTX2KQPn36hIYNG8LKygonTpwQOpwsunfvjlevXiEoKEhnEl1CCoJGKAkh2mFlBcyerZsjgIwBGzZQMmngNm7ciPDwcPz6669Ch5KNu7s7goODce3aNaFDIaRQaISSEKI9KSlAkyZAWJjulBGSSIBu3YDjx3Uz2SVqER8fjzp16sDV1RXbtm0TOpxsFAoFvvrqK7Ru3Rq7du0SOhxCCoxGKAkh2mNsDOzaxbvR6AKRiI9KbtxIyaSBW716NeLi4rB48WKhQ8mRWCzGpEmTcODAAbx7907ocAgpMEooCSHa1aIF5L/8InQUmXbsACpWFDoKokFv3rzBqlWrMHXqVFSrVk3ocHI1evRoiMVibKU6qEQPUUJJCNG4jx8/4tixY5g8eTLKli0Lox9/hFxNbe1UsnUr0KuX0FEQDVuyZAmMjIwwZ84coUPJk4WFBQYPHgwPDw/IdWVJCCFKojWUhBC1k8lkuHHjBry9vXH69GncuHEDCoUCIpEIjDFUrlwZkRERwJw5wIoV2g0ufWp761Zg1Cjtnpto3ZMnT9CwYUMsW7YMM2fOFDqcfAUGBqJ58+Y4fvw4evbsKXQ4hCiNEkpCiFqEh4fD29sbZ8+ehbe3NxISEiCRSHIcadm4cSPGjRvHv/n9d2DWLP61TKbZICUSoEQJYOdOoHdvzZ6L6IRBgwYhICAAjx8/homJidDhKKVly5YoW7Yszp49K3QohCiNEkpCiEoiIyPh5OSE0NBQiEQiiMXiPKfrjI2N8e7dO5QuXTrzyjt3eA/tkBDN1KkUiQDGkNapE2abm2PW//6HirRuUn98+gRERfEqAXI539xVpgyvHZqHGzduoEWLFti6dStGjx6tpWBVt3PnTowcOZL6exO9QgklIUQlsbGxsLa2xuvXr5Hfy4lUKsWoUaOwadOm7DempfHp74UL+S5wdewET5/eLlUKWLcOJ83M0KNnTxgZGeHnn3/G9OnTUbJkSdXPQ9Tr3j3A3x8IDASuXQPu38959LpaNaBVK6BpU6B5c6Bt24xuTIwxdOjQAe/fv8ft27chkUi0/EMUXkpKCqpWrYoRI0ZgzZo1QodDiFIooSSEqOzBgwdo1aoVEhISoMgnEbx+/TqaN2+e+x1evQI2bwbWrwfevePT1AXdoCCV8gSkXj3eTnH4cMDUFO/fv0f58uUBACKRCBYWFli6dCnc3NxgZGRUsHMQ9UpJAQ4dAv74A7hxg38YkEjyXwaRnijK5bxf/OTJwLhx8Lp9G127dsWJEyfQvXt3zcevZtTfm+gbSigJIWoREBCA9u3bIy0tLceRSpFIBGtra4SEhCjXWi4tDThxgpf1uXIFeP+eXy+VZq8ZqVBkJp3Vq/ORqjFjgHbtstyXMQZTU1MkJCRkeXjt2rWxYsUK9O3bl9reaVtKCvDbb8DatUBsLCAWqzY6LRaDiUTwLlUKW776Cn9fv66Xv9OnT5+iTp062LRpE8aMGSN0OITkixJKQojaHD9+HH369Mk1oVy3bh0mT55cuIO/fg0EBQG3bgHR0UByMk8+TEyA8uX5tKedHWBmludhmjZtips3b2a5TiwWQ6FQoGnTpli9ejXatWtXuBhJwVy9ykePw8PVXuxeJhJBXKwYxKtXA5Mm8b8VPUP9vYk+oYSSEKIWjDGsWLEi11p/xYsXx9u3b1GmTBktR5bViBEjsHfv3hw3DqXvSu/atSuWL18OGxsbASIsAlJSgAULgFWreKKn6ZqLbdoA27cDtWtr9jxqdubMGXTr1g0BAQGwt7cXOhxC8qR/H9kIITpHLpdj2rRpmDNnDubPn4+ffvopy+1SqRSDBw8WPJkEgHr16uU62pOeZHp7e6NJkyaYMmVKvhuNSAF9/Ah07AisXs139GujgPeVK3z0+vp1zZ9Ljbp06YI6depg/fr1QodCSL4ooSSEqCQlJQWDBw/G+vXr4eHhgcWLF2Px4sVwc3PLSNxkMhkmTJggcKRcvXr1IMtno4dMJgNjDNevX6eEUp1iYvj61mvXtNvPXSYDEhIAJyfAz09751VRen/vgwcPUn9vovMooSSEFFpsbCxcXFxw8uRJHD58OCNpFIlE8PT0hIuLCwCgfv36aNWqlZChZqhfv36et4vFYohEIsycOROXL1+GWA/X3umkxETAxYWXBBKiraBcDqSm8hi+WEOry9L7e2/ZskXoUAjJE62hJIQUSmRkJFxcXBAZGYkTJ07AwcEh230SExMxYsQIDBs2DH379hUgyuwSExNRqlSpHG8Ti8WoUKEC9u7dCycnJ+0GZuiGDAEOHhQmmfycRAKYmwMPH/J/9YCbmxt8fHwQHh6uV/U0SdFCCSUhpMDu378PFxcXiEQieHl5wcrKSuiQCqRixYp4+/ZtxvfpPcbNzMzw6NEjlCtXTsDoDNCRI0C/fkJHkUkiAQYPBnbvFjoSpQQFBaFZs2bU35voNJrLIYQUiJ+fHxwcHFC2bFkEBAToXTIJIEvMUqkUJUuWxNKlSxEfH4+1a9cKGJkBiooCxo3LXjtUSHI5sGcP8M8/QkeilKZNm6JFixa0OYfoNEooCSFKO3LkCJydnWFra4vLly+jSpUqQodUKFZWVhkbhuzt7XHv3j3MmzcPCxcuxK+//go/Pdq4ofO+/Zbv7Na1yTCxmBe/j4kROhKluLu749y5c3j8+LHQoRCSI0ooCSFK2bBhA/r3749evXrBy8tLJ0oAFVazZs0gkUiwcuVKXLx4EdWrVwcAzJkzB/b29hg+fDji4uIEjtIAPH4M/P238Osmc6JQ8AL5OfWV10EDBw6EhYUF/vrrL6FDISRHtIaSEJInxhh++uknLFu2DNOnT8fq1av1fuezQqFAQkICTE1Ns9329OlTNGnSBP369cO2bdsEiM6AfPcdsG5d/v24hVS1KvDsWWZPcB02e/ZseHp6IjIyEiVLlhQ6HEKy0O93BUKIRqWlpWHMmDFYtmwZVq5ciTVr1uh9Mgnw3dw5JZMAUKtWLfzxxx/Yvn07Dh06pOXIDEhiIrB5s24nkwAQEQF4eQkdhVImTpyIuLg47Nu3T+hQCMmGRigJITlKTEzEgAED4O3tje3bt2PYsGFCh6Q1jDEMGDAAvr6+CAkJ0du1ooLasgUYO1boKPInkQCdOgFnzggdiVJ69OiBiIgI3Lx5k/p7E52i/0MNhBC1e//+Pdq3b4/Lly/j9OnTRSqZBDILsxcvXhyjR4+GQptdXQzF6dN844uqtmzhvb8bNcp+2+zZfLNP9+6FP75cDpw/D3z6VPhjaJG7uztu3bqFK1euCB0KIVnQCCUhJIuwsDC4uLggPj4ep0+fhp2dndAhCebcuXPo0qUL1q5di6lTpwodjn6pWhWIjFT9OKVLA3fuAB8+AC1bZk6hN2oEBAYCe/cCbm6qn+fmTeDrr1U/joYpFArUq1cPrVq1wm49qaNJigYaoSSEZAgKCkLr1q0hEokQEBBQpJNJAOjcuTOmTp2KWbNm4d69e0KHoz8+fFBPMgkA8fG8vI+tLfDTT/w6qRTYtQt4+xaYPl095wkMVM9xNIz6exNdRQklIQQAcPbsWbRr1w61atVCQEAAateuLXRIOuG3335DnTp1MGzYMKSmpgodjn4IClLv8Xx8AA8P4McfATs7YOFCnmCOGQOoo7yTVKr+mDUovb/35s2bhQ6FkAyUUBJCsGvXLnTv3h1OTk7w8fGBpaWl0CHpDBMTE+zZswf379/H/PnzhQ5HPzx6pP7OODNnAi9eAIcO8bWTf/3F1z6qg0wG3L+vnmNpgbm5OYYMGQIPDw/IdbHGJymSKKEkpAhjjGH58uUYMWIERowYgWPHjlF9uxzY2tpi6dKlWLVqFS5evCh0OLovOVn9dR2TkviUd61awPv3PMFUp8RE9R5Pw9zd3fHy5UucPHlS6FAIAUAJJSFFllwux7Rp0zBnzhzMnz8fmzdvhlQqFTosnTVjxgy0bdsWI0aMQGxsrNDh6LZPn9Q/QikS8TaOcjlQvjzQpIl6j69nyxmaNm2Kli1bUn9vojMooSSkCEpJScHgwYOxfv16eHh4YPHixVTTLh8SiQQ7d+5EXFwc3N3dhQ5HtxUrpv7e3T/8ALRuDQwdCjx5AmzdChgbq+/4xYur71ha4u7uDm9vb+rvTXQCJZSEFDGxsbHo0qULTp48icOHD2PChAlCh6Q3qlevjg0bNmDv3r3UrSQvJibq7d9tZQUsXgzs2AEcOACMGgXUrQv88ov6zqGHSz0GDBgAS0tL6u9NdAIllIQUIREREWjTpg3u3LmD8+fPo3fv3kKHpHeGDh2KwYMHY9KkSXjx4oXQ4eim+vXVN0IpkfBEMioKmDaNX3ftGrBmDf++dWvVzyGVAg0bqn4cLTM2NsaYMWOwbds2JOrZGlBieCihJKSIuHfvHlq3bo24uDj4+/vDwcFB6JD01oYNG1C6dGmMHDmSuujkpGlT9R1r7lygeXPexvHjx8zr589X39S3TKbemLUovb/33r17hQ6FFHGUUBJSBPj5+cHR0RFly5bFlStXYGVlJXRIes3MzAw7d+7EpUuXsGbNGqHD0T0WFoA6+p83bswTx40bgbNns96Wmqreqe9mzVQ/hgBq1qwJV1dXbNiwAdT4jgiJWi8SYuCOHDmCoUOHwt7eHseOHUOZMmWEDslgzJw5E2vXrsWNGzfQRN27jvVdv37AsWOAPjVjsCcAACAASURBVIzgSqW8bFCxYkJHUiheXl7o2rUr/P390VodSwAIKQQaoSTEgG3YsAH9+/dHr1694OXlRcmkmi1duhRWVlYYNmwYUlJShA5Ht3Trph/JpEQCODvrbTIJ8BahdevWpRJCRFCUUBJigBhjmDdvHtzd3TFt2jTs27cPxfWwLIquK168OHbv3o3Q0FDMnTtX6HB0y5AhQKlSQkeRP7kcmDJF6ChUQv29iS6ghJIQA5OWlgY3NzcsW7YMq1atwu+//w6xmP6ra4qNjQ1+++03/O9//4O3t7fQ4eiOEiX4RhpdL5ZftSrg4iJ0FCobNWoUJBIJ9fcmgqE1lIQYkISEBAwcOBDnz5/Htm3bMGzYMKFDKhIUCgW6dOmC+/fvIyQkBBYWFkKHpBsePQIaNBA6ityJxcCyZbw3uAEYM2YMvL29ER4eTl2viNbRsAUhBuLdu3do3749Ll++jNOnT1MyqUVisRjbt29HcnIyJk6cSLtt09WvDwwerP6+3mrARCIwMzNg/HihQ1Gb9P7ep06dEjoUUgRRQkmIAQgLC4ODgwNevnyJf//9F87OzkKHVORUqVIFnp6eOHToEHbt2iV0OLrjzz+BMmXU39tbRSLGMEEqhc/Nm0KHojZ2dnbU35sIhhJKQvRcYGAg7O3tIRaLceXKFXz99ddCh1RkDRgwACNGjMCUKVPw9OlTocPRDZaWwKZN6u/trQqJBHE9euBxgwZwdnaGm5sbYmJihI5KLai/NxEKJZSE6LGzZ8/CyckJtWvXhr+/P2rVqiV0SEXen3/+CQsLCwwfPhxydfaz1md9++rO1LdEApibw3T7dly4cAEbN27E4cOHYWVlhYMHD+r9coX0/t4bNmwQOhRSxFBCSYie2rlzJ7p374727dvDx8cHlpaWQodEAJiammLnzp0ICAjA8uXLhQ5Hd2zezNsbCplUisW8TaOXF2BuDrFYjHHjxuHBgwewt7fHwIED0adPH0RGRgoXo4qMjY0xduxYbN++nfp7E62ihJIQPcMYw2+//YaRI0di5MiROHr0KEqWLCl0WOQzbdq0wZw5c/Dzzz8jKChI6HAElZKSgpMnT6J99+5oER0NWFsLk1RKJEDx4jyZtLPLclPlypVx9OhRHDp0CNeuXUPDhg3h6empt33aJ0yYQP29idZR2SBC9IhcLsf06dOxbt06LFiwAAsXLoRIxzY7EO7Tp0+wt7dHYmIibt68iRIlSggdktZ8/PgRp0+fxpEjR3Dq1CkkJycDACwtLfE+NJR30bl6VXuddKRSoGRJwNsbaN48z7vGxMRg5syZ2LJlC9q2bYuNGzeifv362olTjXr27IkXL14gODiYXiOIdjBCiF5ITk5m/fr1Y2KxmHl6egodDlHCgwcPmImJCZs8ebLQoWjcmzdvmKenJ+vUqROTSqUMAJNIJAxAxmX16tX8zsnJjM2axZhIxJhEwhjfsqO5S5s2jIWHF+jn8fHxYXXq1GHFixdnv/zyC/v06ZMGnjXN8fLyYgCYv7+/0KGQIoJGKAnRAzExMejduzeuX7+O/fv3o2fPnkKHRJS0fv16TJkyBadOnUK3bt2EDketwsPDcfToURw8eBDXr18HYwxisTjHqWKRSISIiAhUrlw588qrV4Hhw4HwcPWPVkql/LJqFTBpEl8/WUBJSUlYtGgRVq9ejUaNGmHz5s1o1qyZeuPUEIVCgfr166NFixbYs2eP0OGQIoASSkJ0XEREBFxcXPD69WucOHECrVu3FjokUgCMMXTr1g3BwcG4c+cOypUrJ3RIavHs2TPUrVs3I3nM661ELBbD0dERly5dyn5jcjLw22/AH38AsbE88VMluRSLec3Lvn35cWvXLvyx/nPz5k2MGTMGISEh+O6777Bo0SK9WLe8Zs0azJkzBy9fvkSFChWEDocYONqUQ4gOu3fvHuzt7REfHw9/f39KJvWQSCTCtm3bIJfLMW7cOL0vS5OuWrVq6N27N4C8k8n024cPH57zjSYmwKJFwOvXwK5dfCc4wJNCZdoHSiSZm3wqVAAWLgRevgQOHFBLMgnwguHXr1/HsmXLsH79etjY2OD8+fNqObYmjR49GlKplPp7E62gEUpCdNTly5fRs2dPVK9eHWfOnMk6VUj0zrFjx9CnTx9s2rQJY8eOFToctUhLS0O/fv1w6tSpPHdES6VSvHv3DmZmZsod+N49wM8PCAoCrl0D7t8HZLLs96tWDWjZEmjWjF/atVMuCVXBkydPMH78eFy8eBGjR4/GqlWrYG5urtFzqoL6exNtoYSSEB10+PBhDBs2DK1bt8bRo0dRpkwZoUMiajB27Fj8/fffuHXrFurWrSt0OGqRmpqKHj16wMfHJ8ekUiKRoGvXrjhx4oQqJwGio/n0uFzORzVNTflFAIwxbNmyBT/88AOMjY3x559/on///jq5m/rmzZto2rQpjh49mjGiTIgmUEJJiI5Zt24dpk6dikGDBmH79u0oXry40CERNUlISICtrS0sLS3h5+dnMCNGiYmJqFevHl69epXj7X///TcGDRqk5ag079WrV5gyZQqOHj2Knj17YsOGDahSpYrQYWVjb2+PUqVKwdvbW+hQiAGjNZSE6AjGGH788Ud8++23+O6777Bnzx5KJg1MqVKlsHv3bgQGBuKXX34ROhy1SK+N+urVK9SuXRuSL4qWGxsbo0ePHgJFp1mVK1fGkSNHcPjwYVy/fh0NGzaEh4eHzhVEd3d3x/nz5/Ho0SOhQyEGjBJKQnRAWloaRo8ejV9//RWrVq3C6tWrIS5EmROi+1q1aoWffvoJS5YswdWrV4UORyUymQyjR4/G1q1bsXPnTgQHB+Prr7/OSCqlUin69u1r8EXd+/btiwcPHmDgwIGYNGkSnJycdCp569+/PywtLfHXX38JHQoxYDTlTUgOGOObThMT+fItIyO+bKtyZfWv+U9ISED//v1x4cIFbN++HUOHDlXvCYjOSUtLg6OjIz58+IBbt26hVKlSQodUYGlpafjmm29w+PBh7N27FwMHDgQAxMbGol27drh79y4UCoVB1t/Mi6+vL8aPH4+XL19iwYIFmDlzJoyMjIQOC3PnzsVff/2FyMhIvSh5RPQPJZSEgK/39/bmm0qvX+f/JiRkv1+xYoCNTebG0g4dgBo1Cn/et2/fwtXVFY8fP8bRo0fRsWPHwh+M6JUnT57A1tYWw4YNw8aNG4UOp0BSU1MxaNAgnD59GgcOHMi22ePDhw9o06YN3r17h9evX+tEQqVNycnJWLRoEVatWgVra2ts2bJF8ILoz58/R61ateDh4YHx48cLGgsxTJRQkiItMBBYvx7Yuxf49ImPRKal5f+49PuJREDXrsCUKUCXLgVrxhEaGgoXFxckJibizJkzsLW1LfwPQvTSpk2bMH78eBw7dgy9evUSOhylJCcno1+/frhw4QKOHDmS6+hjQkICoqOjUb16dS1HqDtu3ryJsWPH4vbt25g+fToWL14s6Ohgz5498fz5c9y6dUsnd6QT/UYJJSmSzp0D5swBgoP5FHZOJe6UJZHwSibVqgHz5wNjx/JEMy+BgYHo1q0bzM3N4eXlhZo1axY+AKK3GGPo3bs3AgICcOfOHVSsWFHokPKUmJiI3r17w9/fH//88w+cnZ2FDknnyWQyrFmzBj///DMqVaqEjRs3Cva8nT17Fi4uLvDz84ODg4MgMRDDRav+SZESF8cTvi5dgNu3+XWqJJMATyYB3pxj/Hg+Df78ee739/LygpOTE+rUqQM/Pz9KJoswkUiETZs2QSwWw83NTae76MTHx6Nbt264evUqvLy8KJlUklQqxaxZs3Dnzh3UrFkTnTp1wujRoxEdHa31WDp16oS6deti/fr1Wj83MXyUUJIiw9sbaNAA2L6df6+pyh5+foCVFeDpyTf3fG7Hjh3o0aMH2rdvDx8fH1haWmomCKI3ypcvj61bt+LMmTPw8PAQOpwcxcbGonPnzrh16xbOnTuHtm3bCh2S3qlbty58fHywefNmHD16FFZWVjhw4IBWP0SIxWJMnjwZhw4dwtu3b7V2XlI0UEJJigQPDz4q+fZt5oiipshkvKHHxInApEn8fIwxLFu2DKNGjcKoUaNw9OhRgy+lQpTn6uqKSZMmYcaMGXj48KHQ4WQRHR0NZ2dnPHr0CD4+PrC3txc6JL0lEokwZswYPHjwAG3atMGgQYPQq1cvREREaC2GUaNGUX9vohG0hpIYvNWrgR9+EObcIhEwaJACZctOhYfHeixcuBALFiygBfEkm6SkJNjZ2aFkyZK4cuUKihUrJnRIeP/+PZydnfHq1SucP38eTZo0ETokg3L06FG4u7sjISEBy5cvx4QJE7RSf3bs2LE4e/Ysnj59ajDdmojwKKEkBs3Tk48UCosB2ARPTxHGjx8ndDBEhwUGBsLe3h6zZs0SvJPO69ev4ezsjOjoaPj4+KBhw4aCxmOoYmNjMWvWLGzatAmOjo7YtGkTGjRokO/j5HJ5tq5Eykrv733kyBH06dOnUMcg5Es05U0M1uXLfMpZeCIA45GaSskkyVuzZs2wcOFC/Pbbb/Dz8xMsjoiICLRr1w4fP37EpUuXKJnUoLJly2Ljxo3w9fXFmzdv0KRJEyxduhSfPn3K9TFhYWGoWLEidu3aVahz2tnZoVWrVtiwYUNhwyYkGxqhJAYpKQmwtuY7rzW9ZlJZxsbA3btAnTpCR0J0mVwuR7t27RAZGYnbt2/D1NRUq+d/9uwZOnToAIVCgQsXLqB27dpaPX9R9nlB9IYNG2LLli1o3rx5lvswxtChQwdcvHgRZcuWRVhYGMzNzQt8rt27d2P48OF4+PAh6tevD4C/Vj56xBs7hIQAHz8CKSm8NJqJCVC1KtC0Kb/QfkLyJUooiUH67jvgjz80t5O7MKRS3mHn338LVgCdFD1Pnz5FkyZN0LdvX2xPL0ugBaGhoejQoQOKFSuGCxcuFOmi5EIKDg7G2LFjcevWLUybNg1LlizJKIi+bds2uLm5AQAkEgkmTJhQqDJAqampqFatGvr2HQUbmxXYs4fX5U1J4bcbGfEqFYzxteAiEX89Tf+AXqkSL5E2YQLg6Jh/7V1i+CihJAYnIIC/wOnqX/YffwDffit0FETXbd++HaNHj8bBgwfRv39/jZ/v4cOH6NChA8qUKQMfHx9UrlxZ4+ckuZPJZPj999+xYMECVKxYEZ6enmjcuDHq16+PuLi4jPuJRCLcvHmzwJ227twBRo68huDgxhCJTAAU/DUzvSmElRUwdSowbBhQunTBjkEMByWUxOC4uADnz+vOVPeXypUDIiP5CAAhuWGMYcCAAfD19UVISAiqVKmisXPduXMHzs7OKF++PM6fP48KFSpo7FykYEJDQzF+/Hj4+vqiRo0aiIiIgPyzFzeJRILmzZsjICBAqeoRHz/yGZxt2wCplEEmU31oMf20ZmbApk1A374qH5LoIZp4IwYlLIy3VdTVZBIA3r8Hjh8XOgqi60QiETw9PVG8eHGMHj0aCg2t37h58yacnJxQpUoV+Pr6UjKpY9ILok+ZMgXPnz/PkkwCfM3t1atXsWfPnnyP5eXFmzvs3Mm/V0cyCWROjcfEAP36AYMHA1FRajk00SOUUBKD4ump++sTJRLgzz+FjoLoAwsLC2zfvh3e3t74UwN/NNeuXUPHjh0zkhbq3KSb4uLicODAgVxHIEUiEb777rssU+Gfk8l4+bSuXYF37zT3gTt9vvPQIaB+fcDXVzPnIbpJx996CVFeSgqwcaNuj04CPL5//wUePBA6EqIPOnfujKlTp2L27Nm4d++e2o7r5+eHTp06wdraGt7e3jAzM1PbsYl6zZ49Gx8+fMi1TSNjDNHR0Vi0aFG221JS+Kjhxo38e21sVJTLgdhY3p3sn380fz6iGyihJAYjMJCvD8pPv378k/TAgdlvu3WL39a5c/bbQkN5OY3SpYEff+Sfvl+/BuLjeYmNWbOA4sWVi1UsBs6eVe6+hPz222+oU6cOhg0bhtTUVJWPd+HCBXTp0gXNmjWDl5eX1ksTEeVdvnwZnp6e2aa6v6RQKLB27Vrcv38/47q0NKB/f+DkSe1vUlQo+Mho377AqVPaPTcRBiWUxGAEBio33X3xIn+xa98+6/VmZoCNDZCQkP22KlV4/UhfX6B6dWD6dODmTWD8eKBnTz7Fs3Ahf+FWhkjE4yVEGSYmJtizZw/u37+P+fPnq3Sss2fPwtXVFW3atMGpU6dQqlQpNUVJNGHcOOUbIjDGMHnyZDDGwBjg5gacPi1c+TTG+Ln79gX8/YWJgWgPNfEkBiMoSLlaaB8+8ALjTk5Zr2/Xjn+i3rIle0KZ/r2vL/D0KVCzJi+ens7XF0hMBFatAhwc8n/xlMuBq1fzj5WQdLa2tvjll18we/ZsdOvWDU5f/gEr4cSJE+jfvz+6dOmCAwcOwNjYWP2BErXq1q0bzMzM8PLlS7x9+xYymSzjNpFIBKlUCsYYZDIZFAoFLl26hDVr1qBixRnYvVvAwP/DGH9dHTyYL/Ohzy+Gi8oGEYPx1Vd8WloZ//sfMG0aL8775k3mdc2aAYsX8ykaMzM+WgkAmzcDI0cCFhZALuve0aYNXxs5ZAjw99/5xyAS8Sl6qttGlCWXy9GxY0eEh4cjJCQEZcuWVfqxhw8fxuDBg9GrVy/s3bsXxYoV02CkRBMUCgWioqIQERGByMhIREREZHz9/PlzhIaG4tWrV2jVqjfu3j2M+HjdqccrFvONQYWowU70BCWUxCAwxovsKju106sXcOxY1uTv9m0+Zb1sGS9/0asXcOYMvy0sjJfBaNky92P+/DOf9m7cmBcNVkZwMFDAesSkiHvx4gUaN24MV1dXpUrFMAb8/vs/+OGHfbC2HoCOHXsjLU0MIyPeDrR2bd5Kz8YGoBxT/zHGX7tOn9bNDYoXLmSfASKGgaa8iUFISyvYOqFLl/iLrZMTTyjNzYFGjYCZM/nU9c2b/EXvzBnev7Z2beDgwdyPZ2PDN+UcOaJ8MgkAycnK35cQAKhevTo2bNiAYcOGwdXVFUOHDs12n48fea3BQ4eAa9fSkJraE0BPPHzI8OhR1nUhaWn8X6kUaNiQ78ydMIF6zuurU6eAEyeEjiJnYjEwZgyfSdL18m6k4OhXSgzCp08Fu39sLB+RTF+G1q4dTzDT1z5eupT5Kfrz9ZM5qVGDj2y+fAmMHVuwONSwYZcUQUOHDsWQIUMwefJkvHjxIuP627f5tGLFinxJx7//MqSmZrZkkslESEtDlkvmbbxawZo1QN26PLE8cUI3R7lI7n7/nde61UUKBV+DThUuDBMllMQgKFuu53O+vrz4bqVKPGkMCuKjkwBPKL/+GjA15belpQF+ftmPUb06P45MBnTsyKfKNR03IQCwfv16lC5dGiNHjsS7dwoMHsyXT2zZwmsP8sVMBe+Ekp5A+vjwCgY2NnzEnui+x4/5lLIufwiQSoF164SOgmgCJZTEIEilBf9Unj7i6OTEL5cuZd6Wnjy2bctvu3EjM9lMV706L0EkEvGkMzKy4HGXKFHwxxACAGZmZti5cycuXjRHrVqpOHSIX//ZJmCVpCcljx8DzZsD8+cXfCaAaJeHh+6OTqaTyfhSomfPhI6EqBsllMQgiER8mq4g/v2Xv7j17w9YW/PkMF1cHC9yPnIkUKtW9unuatX4/SUSoEMH4LNZR6WJxQWPmZB0KSnApk3tARxGcrKxxkal5HI+VfnLL3wEVNlKCkS7FApg61bdHp1MJxZn9hMnhoMSSmIwbGwSIJEoX7QgPp5P5fXuzV+Mv6wdeekSL8gLZE0oy5Xj31eqBMydC5Qvz3d/p1+qVFHu/HXrAiVLKh0uIRni43k3p/37+feMFXxqu6AY46OVrVrxtZZEtzx5olynsHQjR/LfKWN8DXlux2Qs6+vf8OHAvn3Aw4c8eX36tOCxKhRAQEDBH0d0GyWURC/J5XIEBwdj3bp1GDx4MKpVq4ZDh+ZCLi9YFSxfX/5pOTiYv0l/7tIlfltqatYXv4YN+Q5YY2Ngzx5eoPzzizIbc6RS/sZMSEElJfFkMiBA+x1Q0ns0t20LqLGtOFGDoKDCPS4uju+8/lK7dvxD75d1d4cP5zM616/zcmqFwRhfRkRFCw0LlQ0ieiEhIQFXr16Fv78//P39ceXKFSQkJMDIyAjNmjXDkCFDYGHRHXPmFOwz0pw5/JKTf/7JufPOpUvKdeTJi0LBa/8RUhByOdCnD38zF6qdnlye2Z40MJCvJSbCCwoCjIyy7txXxv79wLBhgLt71g/VY8bwDy1ftnnv0iUzETxxgpdbK4zoaL7uvGrVwj2e6B5KKIlOioyMhL+/P/z8/ODv74/bt29DLpfD3NwcrVu3xrx58+Do6IhmzZpltI/79AlYsYK/UOk6hQLo2lXoKIi+WbMGOHdO6Ch4UhkTA4wYwXcVU01B4d24UfBkEuDT18OG8SYPGzfy60xNgX79gKlTgenTs95fnaOKQUGUUBoSSiiJ4ORyOe7du5cx+ujn54fnz58DAOrUqQNHR0dMnDgRDg4OaNCgAcS5vHsVK8Zr8C1frtsL0yUSPrrz1VdCR0L0yYMHwLx5QkeRSSbjo/WensCkSUJHQ96+Ldzj4uJ4AXw3t8yEcsgQ/qF3//7sCaU6RUVp7thE+4pcQpmczIv/BgUBjx7xUjCpqTwZKVGCv8k3a8Z3M9KGCc1ISkrC9evXM0Yfr1y5go8fP0IqlcLOzg59+/aFo6MjWrdujYoVKxbo2OPHA7/+qqHA1UQuB6ZMEToKok/kcj4aKNQ0d16+/56PttesKXQkRZsqXbe2buVVKxo2BO7f58nlwYN8aYOmiES8UgExHEUioXz+nH+KPnaMJ5EKBf9jlkozd7mJRPwik/HvxWK+8aJHDz7qRaNJhffmzZsso4/BwcGQyWQoU6YMWrdujZkzZ8LBwQEtWrRACRULM9aoAbi68jpnujlKyWBpmQZXV2qaTJS3cSNfr6iLZDL+AenkSaEjIYV16RIvB+XmBmzfDrRoAcyYofnz0qYcw2KwCaVCwdcarVsHnD7NE8TPEwzG8l5volDwkglr1/J1S87O/EXT1ZUnoiRnCoUCDx8+zBh99Pf3R9h/WwFr1qwJBwcHuLm5wcHBAdbW1rlOX6vi5595P1vdJEJUlDtGjUrGsmXLUJ12NJB8MAasXs0/8OriG7BMxv+/hYVR/291S0pKwvv37xEVFYWoqKg8v3716gSAwv8Ctm3jayaNjfnAS06dwdSJMcDERLPnINplkKnR48d8eujaNb5ejbHCj1alP87XFzh/nu9o270baNJEffHqs5SUFNy4cSNj9DEgIAAxMTGQSCSwtbWFq6srHBwc4ODggCrKFmhUUbNmwOzZfC2lLr0BS6VAmzYKDBrUHAsWzMfhw4fx/fffY86cOShdurTQ4REd5etb+PIs2iKR8FmgFSuUf8ypU6ewePFirFy5Em3bttVccDpCJpMhOjpa6QQxKioKyTnMY5coUQLlypWDpaUlLC0tUb16ddjZ2SEpyQjPnzMUpt0mwEcmFy/mM3LaWqtrYaGd8xDtMKiEUi7nI4pz52auNVLXtGf6cR484OVeFizg5zEyUs/x9cX79+8REBCQMQIZGBiItLQ0lC5dGq1atcL06dPh4OCAli1bolSpUoLFuXAhcOQIfyPWlanvYsWA7dvFqF59PIYOHYLly5dj9erV2Lx5MxYvXowxY8ZASsPf5Avr1wMSCYNcrvni5YUll/Np+UWL8h91evbsGaZOnYoTJ04AAPz9/fUuoWSMIT4+vkDJYUxMTLbjSCQSWFhYwNLSMiNJrFWrVpaEMf3rcuXKwcLCIs9lQevXF26nNwC8egWsXAk0aADs2FG4YxSUnZ12zkO0w2Devd694x1PrlzR7HnSk5OFC/nOuJMnDbcOG2MMjx8/zrL+8fHjxwCAqlWrwtHREUOHDoWjoyNsbGwg0aEmssWL85FkXSoevnZt5t9K6dKlsXTpUkyYMAHz5s3DxIkT8eeff2LVqlVwcXERNlCiMz58AI4dY1AodDeZTPfxI6/dOmhQzrenpqZi1apVWLJkCWT/NRw3MjLC+/fvtRhlzlJTU/Hhw4cCJYhpOWRupqamWRLBr776Cvb29jkmiJaWlihbtqzalv00bVr4ZDLd3Ln538fKim/eAYCKFflm1n79+Pf37/NBF2WULctb2BLDYRAJ5cuXgJMT33yjLYzx/zitWvEpqfr1tXduTUlNTcXNmzezrH+MioqCSCRC48aN0alTJyxcuBAODg56sfaveXP+SXv4cKEj4Tthc+pGUa1aNezcuRPTpk3DjBkz0LVrV3Tu3BmrVq2CjY2N9gMlOsXXNx4KhfLLIUaO5FOXAH9NvHQp+32ePOEdUC5e5OWrAGDTJv5aVrUq/zD28iVw/DhfNvLhg3LnNjLiH+hzSijPnTuHiRMn4tmzZ2CfrUNRKBSIUnPtGIVCgdjYWKWSwvSv479skwWe7H6ZCFpZWWVLCtO/trCwQPHixdX6sxSEtholDBzIB1Q+d+gQ/3fhQj5KnR+RiL8+q9oggugWEWO6tMqs4F6/5i+Er17xxeHaJpUCZmb8hVTfFqRHR0cjICAgY/Txxo0bSE1NRYkSJdCyZUs4OjrC0dERrVq1gumX7RL0iIeHsHXyxo3j68vye/FkjOGff/7BzJkzERYWBjc3NyxZsqTApZOIYXj//j1sbP7G27eToOxn//SEMi6OJ4QjRmS9vV07nkjGxfE+9ukJ5d69/DUsNJSXcmnWjK+ji4gAvv5a+ZGvVq2yzhK9fPkS06ZNw9GjRyEWi6HIoe5Rly5d4OXllesxk5KSlB41jIqKwocPHyDPYZ2LmZlZriOFnyeF5cqVQ7ly5VC6dGmI9Cjjkcv5e1EOubHOkUj435cyHAwlOQAAIABJREFUySfRH3qdUMbFAS1b8hdBIZLJdBIJULkyL+tRvrxwceSFMYbw8PAso4/3798HAFSqVAmOjo4Zm2eaNGkCIwNbHLprFzBqFE/qtLGmMn1H7owZfF1SQd6XPn36BA8PDyxatAipqamYPXs2ZsyYoXJJJaI/Xr9+DWdnZ4SFrURaWlelp7zTE8pNm3j3k4oVsyYYO3fyD76mpryodHpCmZOJE4G//gI6dOCzMMowNua1C+XyT/j999+xcOFCpKWl5ZjgpatatSomTpyYa4KYlJSU7TEmJia5Joc5JYvm5uZFYn3y1Kn8dybk+6EyRCI+Uq5vgzAkb3qdUI4bx0sd6MKmC4mE16w8cqRgyUN8fDzGjx+PevXqYZEaP66lpaUhODg4S/vCt2/fQiQSwdraGg4ODhlJZM2aNfXqk3hhBQQA33zDl0ZoskC0RAKUKsUXyA8dWvhpnZiYGPzyyy/4448/UL58eSxbtgzffPONRkotEd0RERGBDh06ICkpCWlpz/DunfKJUHpC2aEDX9/93XdZ2+m9fp3ZTi+/hLJfPz6V2bYtcPmy8vGvX++Ln37ql+MmlNykjwrmNXr4+df04SpnDx/yNY66TCLhZfjyGJQm+orpqbNn00uS69bl77+V/xnCw8NZgwYNGABWvnx5plAoCv18xMTEsNOnT7N58+YxJycnZmJiwgAwY2Nj1q5dO/bjjz+yU6dOsejo6EKfwxAkJjL2/feMiUSMSSTq/d2LRPzfHj0Ye/1afTGHhYWxAQMGMADMzs6O+fr6qu/gRKc8ffqU1apVi9WoUYOFhYWx4sUL9jc4ciQ/TtOmjO3YwdjVq5m3TZjAWHw8Y6VKMXbnDmO+vtkfL5EwVqIEY61bM3b/PmP//suYWFywGLp3X8kAKH0xNjYW8ik3OO3aqf+1Td2XEyeEfpaIJuhlQhkby1ilSgV/odP0RSRirGxZxt68yf9nuHTpEitbtiyTSCQZL6zPnj1T6udXKBTs6dOnbPfu3WzixInMxsaGiUQilp6Y9unTh61evZpdvXqVpaamqvhsGyZ/f8acnPjvTSpV7fee/ngrK8b27WNMhc8FefLz82MtW7ZkAFjPnj3Zw4cPNXMiIognT56watWqsTp16rBnz54xhaLgf4ufJ5Tt2vGvGzbkt127xtjWrfzrnBLKli2zxnPyJE8+CxrDmTP8w/LcuXOZubk5A5DldS6nS1JSkhafacN29Kjw74W5XcRixqpXZ0wmE/pZIpqglwnlDz/o7icwiYSxUaPyjn/z5s1MKpUysVic5UV1z549Od4/LS2NBQYGsrVr17KBAweyypUrZzzGysqKjR07lm3bto09efJEpVHOoujBA8amTct845RI8v+gIhIxZmSUmUwOGcLYsmX/surVa7DQ0FCNxqtQKNi+fftYjRo1mFQqZVOmTGHv37/X6DmJ5j148IBVqlSJ1a9fn0VERDDGmMoJJcDYkyeMrVrFWKNG/HpHx9wTyhIl+OPatGHs228Zi4xk7MoVxkxMChbD6dOZP1dqaio7cOAAa9euXZ6J5YsXL7T6fBsyhYIxFxfVPyhr6nLunNDPENEUvUsoExMZK11a+P8UeV2MjBiLisoee1paGps6dWqOL6hSqZRNmjSJMcZYXFwcO3fuHFuwYAHr2LEjK1myJAPAihUrxhwcHNjs2bPZP//8w6JyOgkplIQE/ka4ZAljvXoxVqFCzr9bU1PGOnZkbM4cxg4dYuztW/74H3/8kQFgJUuWZJcuXdJ4vMnJyWz58uXM1NSUlSlThq1cuZKlpKRo/LxE/UJCQlj58uVZo0aN2JsvpjeKFVMtofzxRz5jsm4dYw8fZt4vtynvzy8tWvBjTZ9esBhyW5Hx6NEjNmPGDFamTBkGIMsH6ps3b2rq6S2SIiL4h+T0ZTi6cJFIGBs7VuhnhmiS3iWUW7YI/x8jv4tYzNjKlVnjjomJYR07dsyYms7pYmZmxmxtbTNeaC0sLFjPnj3Z8uXLmb+/PyUMWhYXx9dCPnvG2KtXjEVH5z6dvWLFiozfo0QiYRs2bNDKaPG7d++Yu7s7k0gkrFatWmz//v00Sq2jrl69ygICAlhycnLGdUFBQczc3JzZ2trmONJcrpxqCWXlynx6USZjbPbsgiWUYjF/3IYNBYvh1q28n4fk5GS2e/duZm9vn/F/5vz58+p7ogljjLFt24R/L/z8b6lyZcY+fhT6WSGapFcJpULBWOPGurd2MqdL9eqMyeU87kePHrE6derku44IABs2bBjbvHkze/DgASUGemTRokXZfr/jxo3T2hrWBw8esB49ejAAzN7engUEBGjlvER5lSpVYumzEba2tqx///6sRIkSzMbGhn348CHHx3TrVrBRpi8TSoCxX3/l6+oqVixYQtm+PT/W998rf34jI8YK8id/9+5dtmLFChYfH6/Sc0uyUygYGzhQ+PdLkYhPv1+8KPQzQjRNrxLKkBDhE8WCXHx9Gfvnn3+YiYlJtvWSuV3OnDkj9NNMCmH27NnMyMgoy+9SLBazVq1aZZvG1CQfHx9ma2vLALCBAwey8PBwrZ2b5K179+65zlCULl2aOTs7swULFrBTp05ljFYuWFCwtXA5JZQ5XT5PKF1dGTt2jDE3N76cw8WFsZ9+4st2Hj/myzyUPf/XXwv29JIcpKYy5uwsXFIpEvHLkSNCPxNEG/SqoN3Vq8rdr18//uc8cGD2227d4rd17pz9ttBQICgo+/XGxsCjR/xxM2YoF4NEAhw69BI9e/ZEcnJyjh0iviSVSuHv76/cCYhOyan4skKhQGBgIGxtbREYGKiVODp06IDAwEBs27YNfn5+aNCgAWbNmoXY2FitnJ/kztnZOdfb4uPj4ePjg2XLlsHV1RXlypVDyZIl8fz5EY0XqQ4NBT59AubPB06c4LV0hw4FtmzhjSPi4pQ7jpERvz/RHcWK8Y5JnTppv82hWMzfBw8cAPr00e65iTD0KqEMCuKtDvNz8f/snXdYFFcXxt8tdEGaimLHEmMhotiDNYq9JvYWTeyJCtZgiS3YTWLBYDdWYhKxRE1EAUFARAFBFNSoGFFEUToLe74/7rcLSNtF2Fng/p7nPrpz78y8O7PMnr33lKsscfX7SXvNzICWLVklh/f7rK1Z1v6CKkKsXg0YGamv9/nzWpgxYwZ69eoFc3Nz5fbCKjZkZWXBx8dH/RNxBCc1NRVElG97VlYW4uPj0blzZxw5ckQjWiQSCSZNmoT79+/ju+++w44dO9CoUSNs374dMlVr6HFKnR49ehT4GVFARMjKZT2mpaWhRYs0tc5x8CAzHAr6YZybli1znoH37rEf3w0aAIaGrH38MbBoEaBGbnJkZQHt2qkll6MBDA0BT0+W9B5ghl5ZI5GwRPpnzwIjRpT9+TjaQbkyKAMCVCsplZAA3LkDdOuWd3vXrmz/vXvzG5SK1+8blPb2wJw5wLffqqc1OxsICpJg586d+Pvvv/Hq1StERUXB3d0do0ePhrW1tXJsbgMzMDCQf+mXQ1JTUwudhc7OzkZmZibGjRsHZ2fnIsvQlSZGRkZYvnw5YmJiMGTIEHzzzTdo2bIlPD09izRsOGXDw4cPVRonkUhgZmaG69evw9l5LHr2ZF/Q2o6eHjBsmNAqOAWho0Po1Mkd48Z5wNy87D5PilnQAQNY1Z4+fcrmPBwtRdgVd9WRydTzJdq2je2X2xF92zaia9eIevdmx8udtHfPHrYtt7+Qjg7zNdq8mahePXY8Jyf1fEiKKkzz5MkTOnLkCE2bNo0aN26s9Kfifm/lj8GDB6vkIwuA1qxZI4jG27dvU69evQgAde/enadq0RBxcXE0d+5cEolExQbmSSQSsrGxoQcPHij31+ZE1YomlbJKPBztw8/Pj9q1a0cAyNzcnBISiMaOzblvpekraWpatsUdONpNuTEoExLU+4APHsz2GzUqZ1toKNHatURGRkSZmUR9++b0PXjAKknkPsbq1UQPH7KEvyU1KKOjVX+P8fHxFBQUxKO7yyEKQ60oQ0EkEtGgQYMoWp0PRSkjl8vp3Llz1KxZMxKJRDRx4kR6+vSpYHoqIq9fv6Y//viD5syZQ82bN1d+BkxMTJQJvgtqYrGYOnXqlC/iWyZjlcGENhqLa8WlC+JolvDwcBowYAABUAaDDRs2TNl/4wbR5Mks16nCICzJDwmAqGFDNmHz5o2Ab5gjOOXGoHz2TL0Puqkpy6Hm5sZem5uzND69e7PXAQFEGzaw/9euzc7h6pqzv60ti5BTjC+pQRkeLtgl42iQDh06FGgkiEQi0tXVpenTp9P9+/eFlqlEJpPRrl27qFq1amRgYEDLli3jqVtKSFJSEp0/f56cnZ3Jzs5O+eXdsGFDcnBwIIBF3GdlZdH9+/cLNShHjRpVaK7ZdeuET/9SWJNIiDp21PBF5xTKo0ePaPz48SQSiUgqlVLuHyxr167NNz4hgWjLFhahr6gABrD/KyqHSSTMeMw9o2lhQTRiBNE///AZSQ6j3BiUsbHqP+hu3sypDjF0KJuVNDJir9evZ7/QAKLx49k5+vTJeUDevEl06FDOsUpqUIaFCXbJOBqkZcuWeYxIgFU2MjMzo+fPnwstr1Devn1LS5YsIT09PbKysiJ3d3fK4oV2iyQtLY28vLzIxcWFOnXqpPzSrlWrFo0bN4727dtHjx49Ijc3NwJAX3/9NWX/PymtXC6nGjVq5DMmlyxZohxTEKmpRDY22llyViIhunVLU1efUxhxcXE0e/ZskkqleQzJ3O3UqVNFHiMjgygkhMjdnWj2bKIJE4g+/5yVl/3yS5bG6vRp9n3MjUjO+5Qbg1LdJW+A1bAlYstFP/3E6tIq+vr1YzOYJiZE+/blNTadnNjUvY0NUdWqrLVsyY713XfstaqzBQKubnI0SMOGDZUP7YYNG9Lu3bvJ19eXANBvv/0mtLxi+ffff2nMmDEEgFq2bEmXeMFdJZmZmeTv709r1qyhHj16kJ6eHgGsktWIESNo586dFBUVlcdVZdu2bQSAvvnmm3wuLBMmTFC6QIjFYnJ3d1dJh7+/dpXSA5ielStL9XJz1CQxMZGWLVtG+vr6xfroRkRECC2XU4EpNwZlZqb6DsT9+7N9R49mM4W5l7RNTJhBOWAA85P088vp27+/eD22tqppKCooh1NxGDBgAHXs2JFOnz6dZ6apW7du1L59+3LjFxsYGEhdunQhANS3b1+6c+eO0JI0TnZ2NoWEhNCmTZuoX79+VKVKFVL4QA4cOJC2bt1KoaGhhc4orl+/ngDQggULCrzvBw4cIABkYGBAFy5cUEubk5P2LH1LJETNm6tXGYdTuly7do2qVq2qUuEMsVjMy/dyypRyY1ASsbKL6jzwjI2ZQ/upU8x/0tExb39wMJGHBzv2mjU525s2JeraNW8bOZKN27mTvVbMZhbVrK0FuUwcLeLs2bMEgK5duya0FJWRy+V06tQpsrGxIbFYTNOmTdNotR9NI5fLKTIykrZv307Dhg0jMzMzpcH32Wef0Q8//ECBgYEkk8mKPc73339PAGjZsmWF/ohISEigUaNG0e0SRLGkprLnYGlF55a0icUsWJG79AiLn58fGRoaqlTWt379+kLL5VRwypVB+dVX6j9IAwOZMSmTMQMzd9/mzTn1tnv2LPo46vpQisXZZGsbTWFhYUX6RnEqNtnZ2dSsWTMaMmSI0FLUJiMjg7Zs2UJmZmZkbGxMa9eupdTUVKFlfTByuZwePHhA7u7uNHr0aLKysiIApKOjQ59++imtWLGCvL291ZrNkcvltHTpUtJEWqiXL5k7jlBGpVjMIoN5bWbt4MGDB2Rvb19oWU+FX/eAAQOElsqp4JQrg3L3bvV9iFxd2b5BQfn7Bg1ifenpRAYGpWtQikTZJBazLxgzMzMaOHAgbdy4kQICAigzM1OQ68cRhj179pBIJKJ79+4JLaVEJCQk0Ny5c0lHR4fq1KlDhw8fLnc/kmJjY+nw4cM0efJkqlevHimWANu1a0eLFi2iixcvUnJycomOLZfLaf78+QSANikct8uY2FiiJk00H6QjkRDp6xNxF1vtQiaT0fLlyws1KnV0dGjBggVCy+RUcMqVQXn7tjC/yEvazp9PoytXrtD3339PvXr1IkNDQwJAhoaG1LNnT/r+++/pypUrFWLWh1M4aWlpVKNGDZo+fbrQUj6I6OhoGjZsGAGgtm3bko+Pj9CSCiU+Pp48PDxoxowZ1LRpU+UXa6tWrWju3Lnk6elJiYmJH3ye7OxsmjlzJgGg7du3l4Jy1YmPZ+l6NGlMWlqy4CCOdvLdd98pfyy9b1Tu3btXaHmcCk65MijlcqIWLbQv0rGgVrs2C/rJTWZmJgUEBNDGjRtp4MCBZGpqqvz12KlTJ1q0aBGdO3euVL7oONrF6tWrSV9fn16+fCm0lA/Gx8eH7O3tCQANHTpUK/JrJiYmkqenJ82bN49sbW2VX6JNmjSh6dOn08mTJ0v92mdlZdHUqVNJJBLRL7/8UqrHVl0Dy2ahq1t2S+CKIKCxY1m2DY52EhMTQ1WqVKFRo0YpMzbkbv78lwCnjBEREaEc4e4OfP210CqKRiwG1q0DFi0qepxcLkdERAR8fHzg6+sLHx8fPH/+HCKRCLa2tnBwcMCnn36KTz/9FDVq1NCMeE6ZkJCQgDp16mDRokVYsWKF0HI+GLlcjmPHjmHJkiV4/vw5Zs2aheXLl8Pc3LzYfY8ePYrOnTujXr16JT5/amoq/Pz84OXlBS8vLwQHB0MuF8HKqidsbL6AgUFnAPVAZACZDDA0BKpWBVq2BNq0Ya1atRKfHllZWfjyyy9x5MgR7N+/HxMmTCj5wUqBe/eACROAoCBWp7k0ysUr6jJbWBAWLYqGs3OTDz8op0yQyWRwcHDAy5cvcevWLZiYmODo0aP4+uuvkZqaCiJCQkKCSn+fHE6JEdigVZvkZNUirIVsUilznFcXuVxO0dHRtG/fPpo0aVKe3IZNmzalqVOn0qFDh+jff/8t/QvLKXNmzZpFlpaWFcrFITU1ldatW0fGxsZkZmZGW7ZsoYwi8shcvXqVALb0rI4vcUZGBvn4+NDKlSvJwcGBdHR0CABVr25NnTtvo6ZNn5GOjlz5N5i74kfumbbcPoc1axI5OxPFxKj3njMzM+mLL74giURCx48fV2/nMiQri9X97tUrZ4m6pM8vgKhRI6IffyRav54laO/Xrx+95nnQtJLly5eTRCLJNwv56NEj6tChA9nY2AikjFOZKHcGJRHRvHnaWTFC8RBXVN4pDWJjY+nYsWM0c+ZMatGihdLArFOnDo0dO5bc3NwoMjKy3OQ5rMzExMSQSCQiNzc3oaWUOnFxcTR9+nQSi8VkY2NDv/32W77PZHZ2NrVq1YrEYjGJRCJaWURGbJlMRkFBQeTq6kq9e/dW+h+bmprS0KFDac2afTRjxkuqVk3+QcaTRMJcaPr0ITpzJifrQ2Gkp6fT4MGDSUdHh37//ffSuHRlwv377DlZu3be9/r+dRKL8y6Vm5iwcnpeXjmVUM6dO6d87lSrVo1Onz4t7Jvj5OHatWskFovp+++/L7BfLpfz/JMcjVDulrwBIDER+OgjID4ekMuFVpODSAQYGwNRUUDNmmVzjoSEBPj5+cHHxwc+Pj4ICQlBdnY2LC0t8emnnyqXyW1tbSGVSstGRCUlMREICWHLi2lpQGYmoKvLllObNgXs7NiyalGMGDECYWFhiIqKglgs1oxwDRIREYEFCxbgr7/+QpcuXbB582a0a9cOAHDo0CFMnDhROVYikSAwMBBt2rRRun8olrC9vb3x9u1bGBkZwcHBAT169ECPHj3QqpUt9uyRYP58dv1LY2mXaWHH6toVOHAAqF8//5i0tDQMHz4cXl5eOHXqFPr37186Jy9jXr9mn9ubN4E7d4CUFPb51dUF9PUBG5scN4B69XKWuhUEBASgY8eOAACRSAQiwujRo/Hzzz/DwsJCgHfEUfD27VvY2tqidu3auHr1Kn/mc4RFYIO2xJw/L/xsZEHt1181ex2SkpLo0qVLtGzZMuratauyLJyxsTE5OjrS2rVrydfXl/9CLQGpqaxq0hdfENWtm3OPRSI20yOV5sxwKfrq1WNJ8A8dIkpLy3/M69evEwD6888/Nf12NMqlS5eU9c1Hjx5Nd+/epRo1auRJayIWi6lGjRo0fPhwqlatGgEgPT096tGjB61evZr8/PzyLIs/esSKCpTl369UytLi7NyZd7YyOTmZevXqRQYGBpWuLGVkZKTynimaRCIhCwsL+uOPP4SWV6kZO3YsmZiY0KNHj4SWwuGUzyVvBZMmac/St0RCNHBgzjKRUKSnp9O1a9do3bp15OjoSMbGxsovagcHB3JxcaGLFy/Su3fvhBWqxURHs3yjJiY591bdzwLAar4vXMhKe+amc+fO1KVLF2HenAbJysqiPXv2kJWVVZGVPGrVqkUuLi7k5eVFaQVZ4cR8Aw0MNJvMu1cvonfviN69e0cODg5kZGREVythNu9nz54VeN8UPw6++OILio+PF1pmpePXX38lAPSrpmcxOJxCKNcG5Zs3zHFc6DJkEglRrVpEz58LfUXyk5WVRTdv3qRt27bRsGHDyNLSUjnD0LZtW5o/fz79+eef9OrVK6GlCk5CAtG4cSUzIov6bIhERBMnss8rEdEff/xBAOj69euCvt+yJi4ujo4dO0YjR44s1JhUGCZFlaY8cIBdQ02nC5NIiFq1klGbNr3JxMSE/Pz8NHj1tIekpKQi759EIiEzMzPyUNSx5ZQ5Dx8+JBMTExo7dqzQUjgcJeXShzI3sbFAhw7AixdAVpbmzy+RML+5gACgcWPNn19diAj37t1T+mD6+Pjg6dOnAIDmzZsrfTAdHBxgbW0tsFrNcfo0MGUK85MsLb+83EgkgIUFsH8/0KdPNpo1awZbW1t4eHiU/skE4s2bN/D29lb6QUZERAAATE1N8fbtWxT2qBGLxahTpw4iIiJgZGSUp+/oUWDcOGbiCYFIlA2xOBL//CNDt252wogQGCKCVCqFvAiHdYVv5fDhw3HgwAFUqVJFgworF1lZWejWrRuePXuG27dvo2pxjtscjoYo9wYlAPz7L9CtGzMuy8IYKAypFDA3B65cAT7+WHPnLW0eP36cJxfmvXv3AAANGjSAg4OD0shs1KgRRO977JdzUlKAadOAI0dY/tCyDPJSHH/iRMDOzh3z5k3H/fv3YWNjU3YnLUOSk5Ph6+urNCBv3boFIoKNjY0yiKZGjRro2bNnocakArFYjOnTp2PHjh3KbX//DTg6Ch94J5EQOncW4fJl9jdfGTExMUFSUlKRYxRG5dWrV9G1a1cNKat8rFq1Ct9//z18fHzQuXNnoeVwODkINDNa6vz3H1HbtppcFpNTs2byfP5xFYG4uDj67bff6JtvvqHWrVsrfaWsrKzo888/p59//plCQ0PLXT3n93nzhqh9e8374YrFRJ06ZZG5eX2aPXu20JdBZdLS0sjLy4tcXFyoU6dOJJVKCWA+kOPHj6f9+/fny5Has2fPIv0n329///03ERG9fk1UvXpOlRZtaBs2CHHVtYNatWoVe+9atmxJ3t7eQkut0Pj7+5NEIqHly5cLLYXDyUeFmKFUkJUFbN4MuLjkvC5tJBJALicQrcGOHTUxc+bU0j+JlpGYmAh/f3/lDOaNGzcgk8lgamqKLl26KGcw27RpAx0dHaHlqsTbtyxFzJ07mp3VViCRANWrP0NiYms8fXpXK9OvyGQyBAcHK2cg/fz8kJGRAUtLS3Tv3l05C9m4ceMCZ64DAwPRoUMHZXqkopZMFVSpUgUPHz6Es3M1HDkizL0pDB0dIDQUaNZMaCWap3nz5oiMjMy3XSQSwczMDBs2bMCkSZMgkUgEUFc5ePfuHT755BPUqFEDvr6+PEUQR+uoUAalgshIVobs5s3SK0OmOE6zZsCvvwK//DIdhw8fxq1bt9CkSeUqSZaWlobAwEDlMrm/vz9SU1NhaGiIDh06KJfJ27dvD0NDQ6Hl5iM9HejVi/m9CmmwiMUEIh+sWOGPFSuWCCfk/8jlcoSGhioNSB8fHyQnJ8PExATdunVTGpEtWrRQKYdmfHw8Nm7ciDdv3iA5ORlJSUlISkrCu3fvkJSUhJSUFKSkpCA1NRXZuW7E+PHHcfjwyLJ8qyVCKgU++YR9biqb3dS5c2f4+/srX0ulUmRnZ8PQ0BD3799HrVq1BFRXOZg4cSJ+//13hIaGomHDhkLL4XDyI+wEadmRlUXk6ckqYCjyBpZkmUsRQe7gQOThQaRIi5ecnEyNGzcme3t7tUrIVUQyMzMpICCANm7cSAMHDiQzMzMCQDo6OtSxY0dauHAhnT17lt4owpxLgQ0bNlDz5s3p3r17au+7YIE2LaXKychoc6HpcsoSuVxOkZGRtH37dho2bBiZm5sTADIwMKDevXuTq6srBQUFkUwmK3MtGRkZ9Pr1a3r06AnVri3XovuTv+3ZU+aXQ+vo378/ASx3KAAaNGgQXb16lfT09MjFxUVoeRWeY8eOEQA6ePCg0FI4nEKpsAZlbmJimBFhY5PjYykWs3q/EklOjV+pNK+hUb8+0Zw5RJGRBR83KCiIJBIJf6C+R3Z2NoWFhdH27dtp5MiRVLNmTQJYehhbW1uaM2cOeXh4UFxcXInP8cknnxAAMjExocuXL6u83/Xrmk8/U3zLomXLyj5BtFwupwcPHpC7uzuNHj2arKyslIa/g4MDrVy5knx8fARNgu/pKfS9KLqJxUQtWwqfb1bTjBs3jgDQRx99RP/8849y+9KlS0lfXz+f7yyn9Hj8+DFVrVqVRo4cyUvscrSaCrnkXRRJScDt22w5PCoKSE0FMjJYGTIDA6BJE1aCTJUyegCwZs0arFixgkfcFQER4eHDh8prwEQQAAAgAElEQVQ0Rb6+vnjw4AEAoEmTJso0RQ4ODqhXr16xkeRJSUkwNTWFXC6HWCyGSCTCrl278NVXXxW5X1oa0LIlywqgTb55IlE2dHQe4e3bhtDXL91yjM+ePcOVK1eUy9iPHz+GWCxG27ZtlT6QnTt31hrXhN69AS8v7bo/BREQALRvL7QKzREaGorw8HCMGjUqj+9ecnKy8m/4xIkTAiqsmGRnZ6N79+54/PgxQkNDYWpqKrQkDqdwBDZoyz0ymYw6depEDRo0oLdv3wotp9wQGxtLx44do5kzZ1KLFi0IYJGitWvXpjFjxtCuXbsoIiKiwF/kFy9eVI7P3ebNm0dZWVmFnnPZMm1a6n6/ZdP48VEffF3j4+PJw8ODZsyYQU2bNlVeG1tbW5o7dy6dPn2aEhMTP/g8ZUFMjND3QLUmlRKNHy/01dIeDhw4QADIx8dHaCkVjrVr15JYLObR85xyATcoS4EHDx6QsbExTZo0SWgp5ZZXr17R6dOnycnJidq1a6dMNWNpaUlDhgyhLVu2UHBwMMlkMnJxcVGmrMndRCIR9evXr8CykmlprBSi0MZIUU1H5w1lZKh33RITE8nT05PmzZtHtra2ymvRtGlTmjFjBnl4eNDLly9L6S6VLd9/X7yv8/DhbOwXX+Tvu32b9fXunb8vJobo5k1Wa70o/vpL1XtFJKBngFaRnZ1N9vb21Lp16yJ/0HHUIzAwkKRSKX333XdCS+FwVIIblKWE4lc6Lz9WOiQlJdGlS5do2bJl1LVrV9LX1ycAVKVKFTI1Nc1nTCqaRCKhZs2a5fPpOnhQeINRlfbDDw+KvC4pKSl06dIlWrx4MbVr104ZJFG3bl2aPHkyHT58mGJjY8vy1pQZ/foV799qYUGUnU20a1fe7WZmbHtSEtEPP+Tts7Zmx9+0iUhXl+Uefb/98AMb8/XXqt+rGzeEu1bahr+/PwGgPZUxYqkMSEpKIhsbGx70ySlXcIOylJDL5TRixAgyMzMrt1/o2kx6ejpdu3aNVq1apUy0XpRRaW5unqdWdtu22rzcrWgysrS8m+d9Z2RkkI+PD61cuZIcHBxIR0eHAFCNGjVo9OjR5O7uTpGRkRXiS8fSUrXrFBpKdPdu3m1DhhBlZBBt20YUEJC3b9w4dvz+/Qs/ppcXUXIykbGxahpEIiI3N2Gvl7YxZswYql69Onf9KQW+/PJLMjIyoujoaKGlcDgqU7oRAJUYkUgENzc3GBgYYNKkSSolceaojp6eHjp37ozu3buDiIocm52djTdv3qBLly5wd3dHeDgQHCx8Cb/ikeLVq4+wa9dlrF+/Hn369IGZmRkcHBzw448/wsLCAlu2bEFERASeP3+OPXv24NWrV2jXrh1cFNn8yyn//Qe8eqXa2CtXgI8+AqyscrZ16wbcuAGcP8+C6nKXku7WjRU58PUt+HgNG7Ik9ydPsqA9VZBIWGAfJwdXV1ckJSVh7dq1Qksp1/z222/Yt28ffv75ZzRq1EhoORyO6ght0VY0Ll26RABo69atQkupkKxbt06tUn4A6Ouvbxe7lPqhvnkAW051diYKD2ezXXFxROfPE3XsqM4spZyAyWRkZER9+/aljRs30s2bN/P4pslkMnJ3d6fq1asrZ2vHjBkjwN0oPc6dU/0aDR7M9hk1Ku+s5dq1REZGLFds3745fQ8eEAUGFn68NWvY8Tp1Um9G+ZNPhLte2srKlStJV1eXYmJihJZSLnny5AmZmprSiBEjeIogTrmDz1CWMp999hnmzp2LxYsXIzw8XGg5FY6rV68WOvubO52JlZUVBgwYgFmzZiE9vXmxlU2uXmUzmN27591uZsZSDSUn5++ztgZsbNiMGQC4uwOursCffwIDBwKzZgHVqgHe3oC9vWrvTyTKhkTSHo8ePcL58+fh7OwMOzs7SCQSEBH+/PNPNGvWDF999RVevnwJIoJIJEJ6erpqJ9BSXr9Wfay3N0sr1K0be21uDrRowbanpAAhITn3qnZtNgOpuEfvIxYDEycCd+8CuQrBqERCgnrjKwMLFixAjRo14OzsLLSUckd2djYmTJiAKlWqYPfu3cWmT+NwtA6hLdqKSFpaGrVo0YJatmwpSAWUikpWVhYZGhoSwCK6FTOVYrGYmjdvTrNmzaKjR4/S48eP8+z38cdl75unq0skkxEdOpR3jJUVG7Ntm+ozXyJRMK1fvz7Pe/D19SV7e3vl+8V7s7B9+/bVxC0oVeRyOaWnp9Pbt29p8+a3as0O3rxJFBXF/j90KJuVNDJir9evZwEzQE56nz59Cj5Ov36s38lJvdlJgPl8cvKjqOqiTsEBDpGrqyuJRCLy8vISWgqHUyJ4dfkyQF9fH0eOHIG9vT1cXFywadMmoSVVCN68eYPU1FTo6+sra4Z37twZHTp0gImJSYH7pKcD9+6pdvwrV4Bvv2W+eXFxbFtu37xZs5hvXnJyTp/CN08uZ+3t27zHfPeOzaapM4EoEtli27ZhmDt3LqKjo7Fo0SKcO3cOkv9PsxY0Q5uamqr6Cf6PXC5HRkYG0tPTC2xF9ZXWPjl8CWCvytqvXAGcnICaNdls5M2bbHYSYDOVTk6AiQnrk8mAa9cKPs6UKUBmJnDokNqXDyqUM6+UjBw5Etu3b8fcuXMREhKSZ+WAUzA3b96Ei4sLFi5ciO7vL4VwOOWESlcpR5Ns2bIFTk5O+Oeff9CzZ0+h5VQI4uPjYWZmpvKX1NOnQN26qh178GC2XD16NHD8ONsWGgqcPQusWwe8ecPG/PUX63vwgAWSKCqmbN3KDJTx41m1F3Nztp+jI9CuHRuvOtXRokUNREREQCQSFRvkVa1aNTg6Oqpl5GVmZqojCABzK9DX14e+vj709PSU/1elFTU+MNAGmzfbqayjf392X8aMAZYsYQb/4sWsz8SELaEPGQL89BPw/DlQUBGratWA2FjA0xP4/HO1LwXq1gUeP1Z/v8pAcHAw7O3tsXPnTsyYMUNoOVpNSkoK7OzsYGxsDH9/f+jq6goticMpEfynYxkyd+5cnDt3DhMnTkRYWBjMzc2FllTuqVatmlrj09JUH5vbN+/48RzfvAUL8vrm/fVXjm+eh0fO/vPmsRnKU6eg9Nl8/Bjo0UNdYxIwMqqGO3fuAECxUe0AK4H3+PFjpdFmampaIkNPT08PBgYGhY4rq9kmCwtg82bVx/v4sNnhESOA5s2BhQtz+t69Y+VVJ04EGjQAjh4t+BgTJrCSq3tVnxjNQ61aJduvMtC2bVtMmjQJy5Ytw6hRo2BmZia0JK1l3rx5iI2NRUhICDcmOeUbgZfcKzxPnz4lU1NT+vzzz3nUngBERqrnF/chvnnffceiu11ciLp2JRo4kOjiRaKXL1lEsDo63NxYtgATExOlzyhQeCR7o0aNBL3OH8rr1+r7MAYGsmTmMln+/JGbN7M+IqKePQvePzKS6PHj4pOpF9R0dIi+/VbYa6bt/Pfff1SlShWaO3eu0FK0lt9//50AkLu7u9BSOJwPhnsBlTG1a9fG7t274eHhgV9//VVoOZUOAwP1xl+5AjRtWrhvXuvWBfvmffQRsGoVsGIFsGYNG3vmDFuaTUwEtmxRT8fAgb1gZ2cHe3t77Nu3T5mPTlJIuHpGRoZ6J9AyzMyAOnXU2+fKFebHeOtW/vyR3t6sLyOj4Ojtjh2BZs2A/fuZiaguMhnLd8kpnJo1a2Lp0qXYvn07oqKihJajdTx79gxTp07F0KFDMWXKFKHlcDgfDDcoNcAXX3yBCRMmYNasWfj333+FllOpqFpVvfGK9DLdurHm7Z3TpzAeHRxygnUUxqatLTNgbtzIe7ysLOaH2aKFejpMTERwdnbG5cuXYWdnh6ioKHh6eqJdu3YA8huW5T1tEAB06IBi0zvlZvFiQCRi/qnv4+nJ+vT1C3Z7uH6d9a9cWWK53KBUgXnz5qFOnTpwcnISWopWIZfLMXHiROjr68Pd3Z2nCOJUCLhBqSF+/vlnWFhYYPz48cjOzhZaTqXBzCxvRZXieN837+rVnL73ffNy5zb87z/2b4cOeY+nqwvY2bHgD1WpX59Fk48YMQJ169bF5s2bIRaLMXDgQPj7+8PPzw/9+vUDkJN7syQBNtpGp04lmy0UAhMTNpPNKRp9fX1s2rQJ58+fx4ULF4SWozVs2bIFly9fxqFDh2BhYSG0HA6nVOBR3hrE19cX3bp1w+rVq7F06VKh5VQahgxhy8+qll4MDATatmXjzc3zLqdu3gzMnctmI3v1Ai5fZttFIiAggCVB37CBGaZVqwJz5rDl8XHjgCNHij+3RMIijo8dY6+3bt2KhQsX4tGjR6hdu3aesVFRUdi4cSMOHTqEKlWq4M2bN6q9QQ0QHh4Od3d3mJiYwNzcHObm5jAzM8v3f319feU+8fEs0CUrS0DhKiCRsPRS6gQRVWaICD169EBcXBzCwsKgo6MjtCRBuXXrFtq3b49vv/0WGzduFFoOh1N6COzDWelYunQpSaVSunHjhtBSKg2rVxNJpaoHXLi6sv2CgvL3DRrE+tLTiQwM8vaZmLBzRUTklF708iJydFT93GIx0caNOdrfvXtHVatWpQULFhT6/p49e0a3bt0qwyuoPrt37yYAJJFIiiyVqaenR5aWlmRlZUU1a9akIUOS1LpXQrX794u/BllZWZSamkpv3ryhFy9eVOqgvNu3b5NYLKYff/xRaCmCkpKSQh999BF98sknlJ6eLrQcDqdU4TOUGiYzMxOdOnVCUlISQkJCYGRkJLSkCs/FiywXZHnhypWcsoIAsGjRIri5ueHp06eFJnDXNhITE1GzZk21fDv19PRw6tQzDBigvUuAEgnQpYsMUVEsgsjQ0BAymQyZmZl5/s3KysqXO3TZsmVYtWqVELK1gmnTpuHkyZOIjo6GpaWl0HIEYcaMGTh48CBu3ryJZs2aCS2HwylVuA+lhtHV1cWvv/6Kp0+fYsGCBULLqRR0786WrssDVlZAly55t33zzTdITU3Fnj17hBFVAkxNTTF+/HiV81YaGBggKCgI/fpZwM5OveAcTZKdzfKNpqam4sWLF3j06BFiY2Px8uVLvHnzBikpKcjMzCwwEX27gqKHKhGrV6+GXC7HihUrhJYiCJ6ennBzc8OWLVu4McmpmAg9RVpZ2bVrFwGgs2fPCi2lUvDdd0QSifBLpcUtd69eXbD+CRMmUJ06dSgzM1OzF+4DuH37dqFL3bmbWCymCxcuKPe7dUs775VEQjR8ONMYFRVFBgYGKr+/jz76iLIViTErMZs2bSKxWEzh4eFCS9Eo//33H1laWtKgQYMqtesDp2LDl7wFgogwcOBA3LhxA+Hh4ahevbrQkio0T56w6Glt/rRLJCwavKCo9LCwMNja2uLIkSMYM2aM5sWVkJYtWyIiIgJFPWZ2796Nr7/+Os+2779nTVvul0gEmJoCUVGA4k/1zz//xNChQ1Xa/9dff8XYsWPLUGH5IDMzEy1atEC9evVw6dKlSpEuRy6Xo2/fvggLC0NYWJja1b44nHKDsPZs5SYuLo6qVatGAwYM4L9aNcCgQeoF52iySaVEX3xRtP7evXtT69attf6zkpqaSvv376e2bdsWO3u3cOHCAo+RmUnUooV2zVSePJlf59KlS4utYlS7dm2SyWRlfNXLD2fOnCEAdPr0aaGlaIStW7cSALp48aLQUjicMoUblALj6elJAMjNzU1oKRWee/eIdHWFN0wKavr6RA8fFq3/0iVWjvGff/7RzAVTkwcPHpCzszOZm5sTAOrTpw+dOnVK+Tp3E4vFNGLEiCKXgSMiWOS80EalSET01VcFa8zKyqJevXoVGcnerFkzunr1qtb/ENAUcrmcevfuTY0aNarwkc63b98mXV1dmjdvntBSOJwyhxuUWsC0adPIwMCAoqKihJZS4dm8WXjjsaC2fXvx2uVyObVq1Yr69u1b9hdKRbKysujs2bPUr18/EolEZGpqSvPnz6f7ufLqfPfdd3kMLolEQu3ataPU1NRijx8QwNIzCWVUikREn39OlJVVuMaEhASqU6dOgUalkZERNW/enABQx44d6fTp09yXkogiIiJIIpHQxtw5sioYqamp9PHHH1OrVq0oLS1NaDkcTpnDDUotIDk5mZo0aUJt27YtV0EX5ZGsLKL27bVn6VsqJfr0UyJVbYxDhw4RAMGDGl69ekUbNmygBg0aEABq3bo17dmzh1JSUvKNffLkiXJZWCKRUN26denly5cqn8vfn8jYWBijcswYtvxeHLdu3SJdXV16fxZ27dq1JJfL6ezZs9SlSxcCQB9//DEdPHiw0v+tz549m4yNjSkuLk5oKWXC7NmzSV9fn+7cuSO0FA5HI3CDUku4ceMGSaVS+u6774SWUuG5d4+oShUWVS2kMSmRsCXdBw9U156RkUHW1tY0efLksrtARRAUFEQTJ04kPT090tXVpXHjxtH169eLXc4dNGgQASATE5MSzcSHhhI1bKiZeyaRsJnJRYtUN/SJiA4fPpzHoDQ0NKQ3b97kGePr60v9+/cnAFS3bl366aefCjTCKwOvXr0iMzMzmjp1qtBSSp2zZ88SANquytIDh1NB4AalFrFmzRoSi8Xk6+srtJQKz7VrRHp6whmVYjFbyg0MVF/7hg0bSEdHh/7777/SvzAF8H6QTb169eiHH35Qa5bx2rVrVK1aNfL29v4AHUROTszYK6vZSrGYGa5+fiXTOHv2bBKJRCQWi4usbhQaGkpjx44liURClpaWtGrVKnr9+nUJr0z55aeffiKRSEQhISFCSyk1FMGW/fr1436znEoFNyi1iKysLOrcuTPVr1+f3r59K7ScCs/lyywYRtNLqRIJkaEhUUltq8TERDI2NqYlS5aU7gV5j4KCbDw9PSmrKIfCIiitL1c/PyIbm5xrWZqzkvPnM8O1pGRkZFDHjh1VNvgfPnxIM2fOJH19fapSpQo5OTlRbGxsyQWUMzIzM6lZs2bk4OBQIYwvuVxOffv2perVq9OLFy+ElsPhaBRuUGoZDx8+JGNjY5owYYLQUioF168TmZtrzqiUSIgsLIiCgz9M9/z588nU1JSSkpJK50L8H0WQTd++ffME2URHR5fqeT4UmYzozz+JevYsuWGpmJ02N2eJ7x8/Lh1t7969U9vHNS4ujpYuXUpVq1YlXV1dmjp1Kt27d690BGk5Fy5cIADk4eEhtJQP5qeffiIAdP78eaGlcDgahxuUWsiBAwcIAJ0sKPEdp9R5+ZJoxAhmXIhEZWNIKo47cGASTZmykN69e/dBmh8/fkwSiYR+/PHHUrkG8fHxtH79emWQjZ2dHe3du7dc+Pfdv0+0YAGRnR2Rjk7ONdfRydtyG501ahANHkx09CiRNmWuSUxMpPXr15OVlRWJRCIaMWIEBX/or49ywIABA6h+/frlOho6PDyc9PT0aM6cOUJL4XAEgVfK0UKICCNHjsQ///yD8PBwWFtbCy2pUuDhAUybBrx9CxRQirnEiMWAmRmwZw8QEbEWLi4uMDIywpo1azB9+nTo6+uX6Lhjx46Fv78/oqOjVa6Z/T5BQUHYsWMHTpw4ofzczZw5E+3bty+XVUxkMiAyErh5E3j8GMjIYNsMDABjY6BVK6BNm5xqN9pKeno6Dh06hA0bNuDBgwf47LPPsHjxYnTv3r1c3pfiuH//Plq0aIGVK1di6dKlQstRm/T0dNjb2wMAbty4UeK/aQ6nXCOsPcspjISEBKpVqxb17NmT563TIPHxRKtWsRmsD/HRU+xXsybR2rVECQns+Ldu3SKARQGLRCKysrKi3bt3lyiFzM2bNwkAnThxQq39UlNTad++fR8UZMPRDDKZjI4fP06ffPIJAaB27drR77//XiGfCfPnzycjIyN69uyZ0FLU5ttvvyU9PT0KCwsTWgqHIxjcoNRi/v77bwJAW7ZsEVpKpUMmI/r9d6Lu3XN87USiwvNXSqU5y9oSCVGvXkSnT+dPiJ2amlpgqb569erR4cOH1Q546dGjB9nb26sU0BATE5MnyMbR0fGDgmw4mkMul9OFCxeoa9euBICaNm1K+/bto4yMDKGllRpv3rwhS0vLcuc//tdffxEA2rZtm9BSOBxB4QalljNv3jzS1dXlv3wFJDWVVWzZsYPoyy9ZYvSWLYmaNGH/duxINGUK0c6dLA1QcW5g9erVy2dQKozMJk2a0KlTp1SOeD1//jwBKDQdT3kJsuGojr+/Pw0ePJgAVid869atpR6cJRRubm4EgAJLkk9LAF68eEE1atQgR0fHChGlzuF8CNyHUstJT09Hu3btADB/N+6bU/4ZPHgwzpw5g4L+9MRiMeRyOVq1agVXV1c4OjoW6TNHRGjRogVsbGzg6emp3P7q1Svs27cPbm5uePToEezs7DBr1iyMGjUKhoaGZfK+OJolMjIS69evx9GjR2FiYoI5c+Zgzpw5sLCwEFpaicnOzoadnR2MjIzg5+en1f6iRIRBgwYhICAA4eHhsLKyEloShyMsgpqzHJUICwsjXV1dmj9/vtBSOKXAkiVLSEdHJ98sZe6mqAutyszHvn37CABFRUVRYGAgTZgwQVnJZvz48RQQEMBnTyow//77L82ZM4cMDAzI0NCQ5s6dS0+ePBFaVonx8vIiAHTkyBGhpRTJzp07CQCdOXNGaCkcjlbADcpywpYtWwgA/f3330JL4Xwgv/76a5HGZO6mij/ZmzdvyMTEhKpVq6b0x3R1deVBNpWMly9f0rJly8jU1JSkUilNmjSJIiMjhZZVIoYNG0a1a9em5ORkoaUUSEREBOnr69PMmTOFlsLhaA18ybucIJfL0bt3b9y9exfh4eEwNzcXWhKnhNy+fRutW7cutF8ikUAqlWL79u2YMmVKoct+Dx48gJubG/bt24fXr19DJBLh4MGDGDNmDCQSSVnJ52g5SUlJ+OWXX7BlyxY8f/4cQ4YMweLFi5WuM+WBhw8folmzZliyZAlWrlwptJw8ZGRkoH379sjMzERwcDB3IeFw/o9YaAEc1RCLxThw4ADS0tIwbdq0Av3vOOWDpk2bQiwu+E9PJBKhadOmuHXrFqZOnZrPmMzOzsa5c+fQr18/NG7cGHv37sXkyZNx48YNGBoaIiYmhhuTlRxjY2M4OTnh4cOHcHd3x507d9C+fXv06NEDf//9d7l4djRs2BDz58/Hhg0b8OTJE6Hl5GHp0qW4e/cujh49yo1JDic3wk6QctTl5MmTBIAOHjwotBTOB1C3bt0Co7xNTEwKrAGsSiWbb775hiwsLMpFdRuO5sjKyiIPDw9q06aN8rNz8uRJrU8X9e7dO7KysqJRo0YJLUXJpUuXCABt3rxZaCkcjtbBDcpyyIQJE8jY2JgePnwotBROCRk4cKDSiJRIJFS1alXavXs3GRsb07hx45Tj1AmyefjwIYnFYtq5c6cm3wqnnCCXy+nvv/+mnj17EgBq3Lgxubu7U7o21Z58D0XA2bVr14SWQvHx8VSzZk367LPPKmRieQ7nQ+E+lOWQd+/ewdbWFtbW1rh69WqJy+5xhGPp0qVwdXUFEcHBwQFHjx6FtbU1jh49irFjx2Lq1Km4ffs2goODUa9ePcyYMQNffvklqlWrVuRxR44ciZCQEERFRfGlb06hBAUFYf369fjjjz9gZWWF+fPnY9q0aTA2NhZaWh7kcnmetGmFuYqUNUSEoUOH4tq1awgLC0OtWrUE0cHhaDXC2rOckuLr60tisZjWrFkjtBROCfD09CSJRELr1q1TLj3GxMSQk5MT6erqEgD69NNP6cyZM2otTQYFBREA+v3338tKOqcCcffuXZo8eTJJpVIyNTUlFxcXrcsO4OvrSwBo//79gmnYvXs3AaA//vhDMA0cjrbDZyjLMS4uLli/fj38/f1hb28vtByOGhARkpOTYWhoiL/++gs7d+7EhQsXYGpqinHjxsHT0xOWlpbw9/eHrq6uWsfu2rUrsrKy4OfnV0bqORWNp0+fYuvWrdi9ezeICFOmTIGzszPq1asntDQAwOjRo3H16lXcv39f47OoUVFRsLOzw/jx47F7926NnpvDKVcIa89yPoTMzExq27YtNWnSRGvztXEKRhFkU79+fQJAbdq0oX379ikDaoKDg0lHR4ecnZ3VPranpycBID8/v9KWzangvHr1ilauXEnm5uYkkUho/PjxdOfOHaFl0ePHj8nAwICWLFmi0fNmZGSQnZ0dNW3alD9jOZxi4DOU5Zx79+4pfz27ubkJLYdTBESEGzduYMeOHThx4gQA5vM4c+ZMtGvXLl+KoC1btsDJyQkXLlxAnz59VD6PXC7Hxx9/jObNm+PUqVOl+h44lYOUlBTs2bMHmzZtQmxsLAYNGoTFixejY8eOgmlasWIF1q9fj8jISDRs2LDAMcnJwK1bwM2bwJMnQHo6kJUF6OsDJiZAy5ZAmzaAjQ2gSlXHRYsWYevWrQgICICdnV0pvyMOp4IhsEHLKQXc3NwIAHl6egothVMAqamptG/fPmXaFkUlm/j4+CL3y87OJkdHR6pevTrFxcWpdc5ffvmFRCIRRUdHf4h0TiUnIyOD9u/fTx999BEBIAcHB/rrr78EKeWZnJxM1tbWNGzYsDzbQ0KIpk0jsrEhEomIACKxmEhHhzWpNOf/AGtVqhB17Urk5kaUlFTw+S5fvkwikYjWr19f9m+Ow6kAcIOyAiCXy2nAgAFUrVo1tQ0PTtmhCLIxMzNT1uVWN8gmLi6OatSoQX369FErVUlaWhpVr16dl4bjlArZ2dn0+++/k729PQEgW1tbOnbsGMlkMo3qUJQtvXjRmw4fJrK3ZwaiVJpjLKraRCLWDA2JZs8miojIOc+rV6/I2tqaevTowVMEcTgqwg3KCsKLFy+oevXq1L9//wJnD65fv06XL18WQFnlIisri86cOUOOjo4EgMzMzMjJyYliYjSXLh8AACAASURBVGJKfMyLFy8SANq0aZNa+61atYoMDAyKnQnlcFRFLpeTl5cXffbZZwSAGjZsSG5ubpSWlqax8zdrNod0dF4oZyLVNSQLagqDdOpUosREOQ0fPpzMzMzo6dOnGnlfHE5FgBuUFYgzZ84QANq1a5dyW0ZGBi1evJhEIhHVqlVLQHUVm/j4eHJ1dS00yOZDcXZ2Jh0dHQoODlZLk4GBAa1atapUNHA4uQkODqbPP/+cRCIR1ahRg1xdXent27dldr5374i+/loxu5hdKobk+00iITI1TSHgM/rtt9/K7L1wOBURHpRTwZgxYwYOHjyIkJAQACzoIzw8XFm/NyUlhdefLSWICEFBQdi5c2eeIJtZs2bB3t4+X5DNh5CZmYlOnTrh7du3CAkJUTl1yqxZs+Dh4YEnT55AX1+/1PRwOAru37+PjRs34uDBgzA0NMTMmTPx7bffokaNGqV2jhs3gCFDgLg4QC4vtcMWiEgkB5EY06cDP/8M8LoRHI5qcIOygpGSkoLWrVsjPT0dL1++RFZWFrKzs5X9oaGhaNWqlYAKyz9paWk4duwYdu7ciZs3b6J+/frKSjaWlpZldt7o6GjY2dlh+PDhOHDggEr7xMTEoEmTJti9eze++uqrMtPG4Tx79gzbtm2Dm5sbZDIZvvzySyxYsAANGjT4oONeuQL07w9kZgK5HmVljkgEDBgAnDzJosQ5HE7RCFPHilNmJCcnw8LCAk+fPkVGRkYeYxJgaYY4JSMmJgbOzs6wtrbGlClTUL16dZw5cwYxMTFYuHBhmRqTANC4cWPs3LkTBw8exJEjR1Tap1GjRhg6dCg2b94MeVlP7XAqNdbW1ti4cSOePHkCFxcXeHh4oHHjxhg7dizCwsJKdExvb8DREcjI0KwxCbBF8HPngGHDAJlMs+fmcMolQq63c0oXT09PMjMzI4lEQgDyNalUyks1qklZBNl8KOPGjSNjY2OVNfj7+/O0UhyNk5KSQj/99BPVrVuXAFC/fv3I19dX5f1v3WIR2KUVeFPSJhIRjR5NJECmJA6nXMGXvCsAKSkpmD9/Pn755ReIRCIUdkslEgnGjBmDQ4cOaVhh+ePVq1fYu3cv3Nzc8O+//6JNmzaYNWsWRo0aBQMDA0G1vXv3DnZ2djA3N8e1a9dUKs3YuXNnSKVSeHt7a0Ahh5ODTCbD8ePH4erqisjISHTu3BmLFy9G//79C/UzzsgAWrUCHjzQ/MxkYRw4AEycKLQKDkd74UveFQBXV1f88ssvAFCoMQkA2dnZiIiI0JSscgcRITAwEBMmTEDt2rWxYsUKODg4IDAwEMHBwZg8ebLgxiQAmJiY4NixY7h16xaWL1+u0j7Ozs7w8fFBUFBQGavjcPKio6OD8ePHIzw8HKdPn0Z2djYGDhwIW1tbHDlyBFlZWfn2+f57IDpae4xJkQiYPRt49kxoJRyOFiPo/CinVIiLi6PJkyeTSCQiqVRa4HK3olWpUkWQKhfaTEpKCu3du1dZyaZ+/fq0fv16rc/fuGHDBgJAly5dKnZsVlYWNWrUiL744gsNKONwCkcul5O3tzf17dtX+fe2Y8cOSk1NJSKioKCcijfa1CQSor59+dI3h1MYfMm7AhEREYHFixfj7NmzkEgk+QJyFLx48QLVq1fXsDrtIyYmBrt27cL+/fuRmJgIR0dHzJo1C46OjpBIJELLKxa5XA5HR0eEh4cjNDS02Hu6a9cuzJ49GzExMR8cecvhlAa3b9/G+vXrcfLkSVhYWODbb+fi6NFFuHdPojWzk+9z6hQL1OFwOHnhBmUFxM/PD05OTggMDIRYLM4X3evr64suXbrk248IePQIuHkTCA5mud8iI4G0NJayQ0cHMDAAPvoIaN8eaNOGNRsbtiRUHsjOzsb58+exc+dOXLhwAebm5vjyyy8xffp02NjYCC1PbeLi4tCqVSu0bdsWZ8+ehVhcuBdLamoq6tati7Fjx+LHH3/UoEoOp2gePHiAjRs3Yu/ee8jKuiK0nEKRSICOHQFfX6GVcDhaiLATpJyyQi6X0+nTp6lJkyb5lr337NmTZ+ybN0Q//kjUqFHO8o6OTtHLP7n7GzQg2rKF6PVrgd6sCrx8+bLASjaKZbbyzPnz5wkAbd26tdixy5cvJyMjI0pISNCAMg5HPQYPTiWJpGyq4JRmu3NH6CvF4Wgf3KCs4MhkMtq7dy9ZWVkpDcpJkyYREVF0NNFXXxHp6zOfpZL6LSn21dMjmjyZKCpK4Df9f+RyOQUEBND48eNJT0+P9PT0aOLEiRQYGCi0tFJn3rx5pKOjQyEhIUWOe/HiBenp6dG6des0pIzDUY24uJya2trcpFKiWbOEvlocjvbBl7wrCWlpadi6dSuWLVuGVq0+wcSJN7F4MYuiLCDIssRIpWz5e80awMmJLRFpmtTUVBw/fhw7duxASEiIxirZCElGRgY6duyIlJQU3Lx5E1WqVCl07LRp0+Dp6Yl///0Xenp6GlTJ4RTO1q2As3PZl1YsDQwNgcRE5gbE4XAYPG1QJcHAwABLly5FcPBb6OoGYd48luutNI1JgB1PJgMWLQI6dADu3i3d4xdFTEwMnJycULt2bUydOhU1atTA2bNnNVbJRkj09PRw/PhxPHv2DN9++22RY+fPn4+4uDgcPXpUQ+o4nOLx9y/eF3v4cDZP+MUX+ftu32Z9vXvn74uJYb7hCgwNWWqie/eA9HTg1SvAywto1Eg1rampmn22cTjlAW5QViLOnwc6daqCkBDNTBvevg3Y2gJ//ll258jOzsaZM2fQt29fNG7cGAcOHMCUKVMQHR2N8+fPo3///uUiYrs0aNKkCbZv3459+/bh+PHjhY5r2rQpBg0ahE2bNoEvUHC0hYCA4vNOXr3KZjC7d8+73cwMaNkSSE7O32dtzQIHr/w/1sfIiB1nyhTg55+ZATp5MhAYyAxNVQkOVn0sh1MpEHjJnaMhTp5kedQ0nd9NJGKl0w4fLt33836QTdu2bWn//v0VIsjmQ5DL5TR69GgyMTGhhw8fFjrOx8eHAND58+c1qI7DKZiEBNWfKaGhRHfv5t02ZAhRRgbRtm1EAQF5+8aNY+fo35+93rqVKCmJBROW9Lmmo0M0c6aw14zD0Tb4DGUl4Nw5YPRo9ste0xNSROy8Eyaw/G0fdqz8lWy6du2KoKAg3LhxA5MmTdKKSjZCIhKJsGvXLlhYWGDMmDGQyWQFjuvSpQvatWuHTZs2aVghh5OfW7dUH3vlCktdZmWVs61bN5bm7Px5lsostwtxt27MFcfXl6U9mzoV8PBgKdJKikwG8KJTHE5euEFZwQkNZUl4hTAm32fUqJI9hFNTU7Fv3z60bdsWHTp0wLVr17BmzRrExsbiwIEDsLe3L32x5ZiqVavi2LFjCA4OxsqVKwscIxKJ4OzsDC8vL4SEhGhWIIfzHi9fqj5WsXTdrVvOtu7dAW9vwM+PPec+/TRvX0gI8O5djrEZHQ3s3Am8fs18yW/cAPr1U0/zixfqjedwKjrcoKzAyGTAuHHML0loY1KxWDRuHHOCV4X3g2ysrKxw9uxZREdHY8GCBRU6yOZDad++PVavXo0ffvgBXl5eBY4ZOnQoGjRogM2bN2tYHYeTl7Q01cd6e7NnmsKgNDcHWrRg21NSmPGo8KOsXRto2DDHCLW2Zv8uWsR8LidMAIYOZcbmmTMFB/QUhqrPMQ6nssANygrMunVARETxju6aIjsbePCARVcWPoYF2Tg6OuYLsjl37lylCrL5UBYuXIju3btj3LhxePXqVb5+qVSKefPm4cSJE3jy5IkACjkVkRMnTqB69eoYPXo0jhw5goSEhGL3UecHb2IiW3lRGJRdu7Jni58fe+3tnWNQKv5VGJSKQlKZmUDfvsDZs2yZfMAA4PlzYNky1XUI/SOdw9E2uEFZQQkNBVav1r6HnlwOrF+ff+k7Pj4erq6usLGxwaBBg5CQkID9+/cjNjYWGzduLJdlEYVGLBbj8OHDyMzMxOTJkwuM6J48eTJMTEx4KUZOqZGRkYH4+Hh4eHhg3LhxqFatGjp06IAffvgBYWFhBX4O1XV9vnIFaNoUqFmTGY03b7LZSYAZlK1bAyYmrE8mA65dY30K29bfn0WEK0hLY/vZ2amuoZK7a3M4+eAGZQVl8WKhFRSOWAwsXMiCbAICAjB+/HjUrl0bK1euRLdu3XiQTSlSq1YtHDhwAGfPnsX27dvz9VepUgUzZszAL7/8gsTERAEUcioarVq1AsBWG4CcYDoXFxfY2tqiVq1amD59Os6ePYvU1FQAgIWFeufI7UfZrRszBhUojEcHh5xgHYWxGRZW+DFFIvWSqqurmcOp8AgaY84pEx480Hx6oJK0Zs2GEQBq0KABbdiwgeLj44W+dBWWb775hnR1den27dv5+v777z/S1dWlDRs2CKCMU9FIT08nsVhMAAptUqmUAJCOjg59/PHHtGqVm1rPDmNjIpmM6NQpouxsIkfHvP3BwUQeHkzPmjV5+/z8iOLj2TEU2wwMiJ49I/r7b9XTBn31lbDXmcPRNnjpxQrIokXA5s3a4ztZECJRNurWPYsdO6RwdHTkfpFlTHp6Ojp06ICMjAwEBwfDyMgoT/+UKVNw8eJFPHz4ELq6ugKp5FQUmjZtivv376s8vmXLloiPD0NcnOrnCAwE2rZls4rm5kBSUk7f5s3A3LlsNaRXL+Dy5Zy+jh3ZDGdICHO/IWJlYjt0YP6YAQHFn1skIri5ifD116rr5XAqOnzJu4KRng7s3q3dxiQAEEnw6tVgdOvGg2w0gb6+Po4fP44nT55g3rx5+frnz5+PZ8+e4cSJEwKo41QUXr9+DW9vb1SpUgWiYuooKvoXLFiAkJAQtG+fEzSjCleusPG3buU1JgG2BC4Ws5RA/v55+65fB3r2ZH1HjgBHjzI/y27dVDMmAYBIBFfXz/HNN9/g9OnTePv2rerCOZwKCp+hrGBcugT06VP8uOHDgd9+A0aOBE6ezNunKJnYpw87Xm5iYoC3b1k+tytX8uaCU3DhAougVIU//wQGD1ZtLOfD2bt3L6ZOnYqTJ0/i888/z9PXv39/xMbG4vbt28UaA5zKjUwmw/379xEWFpanxcbGAgAkEonSh7IgJBIJTExMcPz4cfT+f64eV1fAxUX7fwwDgFSajfHjZ+Pq1Yt49OgRxGIx7O3t0bNnT/Ts2ROdOnWCvr6+0DI5HI3CDcoKxrp1wPLlxT+ULSxYMuFffgFmzMjZbmYGvHoFpKYC27cDS5bk9FlbA7GxbDnJ2ZkZlHXqAGPH5j12YiJw717xWqVSdvxVq1R/f5wPg4gwatQoXLx4EaGhoahXr56y78qVK+jRowcuXbqEzz77rNBjZGezqFiRCNDXB/gEc8WFiPDixYt8huPdu3eRmZkJAKhTpw5atWqVp0VHR2PIkCEFHlMkEqFLly44fvw4atWqpdz+5AlQv772ZaZ4H6mU5dPdv5+9fvjwIS5fvozLly/Dy8sL8fHx0NfXR5cuXZQGpp2dHV+J4VR4uEFZwRg2DDh9WrVoxdBQQFcXaNYsZ9uQIcCJE8CuXcynqEOHnL5x44DDh1nOtnPnmEFpackSBJcEkYglEr5woWT7c0pGYmIiPvnkE1hbW8Pb2xtSqRQAMx7atm0LS0tLXLx4EQArTxcYyNKyBAay5cXc6VYAlp6lTRvA3p75tLVvD9Stq+l3xflQ0tPTERkZmc94jI+PBwAYGhqiZcuWeQzHli1bwszMLN+xnj59irrvfQjEYjGICC4uLlixYkWBBtagQcBff7FSidrMjRvss/4+crkc4eHhSgPT29sbKSkpMDU1Rffu3ZUGZtOmTfkqAKfCwQ3KCkatWixBryps2wZ8+y3L5aZwht+2jT0oV61iRqOZWY4BsWcPMHEim9189+7DDUqAOdO/esWMS47m8Pf3h4ODA5YuXYpVuaaIjx07hjFjJsLVNRqenvWU/mc6OszPrChyj+nWDZgzhxkI/7dXOVoCEeHJkydKgzE8PBxhYWG4d+8e5HI5RCIRbGxs8hiNrVq1QsOGDSH+X3t3Hhd1nf8B/PWdGQ6VVLzPhMEbwQM1zAs0TXMzTTu27G5zLbvstN2fdtq9abm2W7napqXmWqt5bKWIpJkHpqLmAaICad6gcQzM9/fH2xkYhmOGOb5fhtfz8eABzPmhRnjN53i/XdzkqKoqGjVqhLwrmxuNRiMaN26MJUuW4Lrrrqv0fq5u2dGKwSDbgVztVlpUVIRt27bZA+bWrVthsVjQtm1be7gcPnw42tpa+BDVYgyUASQvT2aLXHXTTbKH8Y9/BJYskct275buEbNmAefPy23WrpXr0tMl/F1zjXyflCRf5+fL8x47Jo/z6qvutSU7dQpo0cL125N3vPrqq5gxYwY2bNiAhIQEXL4MvP56CV5//SKs1iYwGmu+n81235YtgcceA6ZNk+Vx8q+8vDykpaU5zDju3bvXfoikcePGTsvV0dHRCAsL8/i5Bw0ahM1X2tcMHToUS5YsQatWraq8j9UKdO4sM+Pu1IT0p08/lZaNNXHp0iWkpKTYA+bPP/8MAOjatas9XCYkJFQ460ukdwyUAeTXX2WG0lWNG0tA/OQT4M9/ltnC06flQM2338qJx02bpAh5u3bAiRNSZsNWNP2VV4DsbOCXX6RrxOjR8jhbtkiHCldfWenp0m+3KqqqIjU1FUuXLkVKSgqWLVuG9u3bu/7DkpOSkhIMHz4cR44cwYcf7sPUqY2QleX9P+SKAnTsCCxaBPTv793HJlFSUoKMjAyn5eqMjAwAMkPYtWtXpyXrdu3a+Wzpddq0aZgzZw5mzpyJv/zlLy7vIdywQU5h643JJK/flBT3TqNX5fTp00hKSrIHzPT0dBgMBsTFxdkD5sCBA9nggWoFBsoAcvRo9cGsvJ07gQYNgK5dgfHjZf9keLh0lnjzTWDYMNkbd9ddwL//DYwaBVzZXlehadPk0M748TL76Yp9+4Du3Su+Li0tDUuWLMHixYuRmZkJRVGgqip27dqFXr16uffDkpNDh7IRG/sNCgsnezQjWR2jUYLqc88BM2dyttITZ8+etS9T2z7S0tKQn58PAGjZsqXTrGPXrl39fur4woULOHXqFLp06eL2fadMkQODepqlDAkB0tLkzZGvZGZm2sPl+vXr8dtvvyEkJAQDBw60B8y4uDj7vmciPWGgDCAnTrh/GOKdd6Sob5s2cuK6Xz8p/AsAN9wArFwpM5ezZ8uhHFvYrEyLFrKEXXYmszqHDgGdOpV+f/DgQSxduhSLFy/GoUOHKixBcvHiRTR0Z32fnJw6JTNBBw6osFr9s4lVUeQ1tm6dvJaochaLBQcPHnSadczOzgYABAcHIzo62umQTMuWLTUeuecuXZLDgjk5+gmVtj3n/qKqKtLS0hwO+OTl5aFhw4YOB3y6devGAz6kCwyUAeTMGaB5c/fuM2aM7Jm84w4JlGvWlAbBhg2Bc+fk5Pf778uS+sCBVT+eLVC+/jrwwguujeH4caC4+CiWLVuGRYsWIS0trco6do0bN8b58+fd+CmpvOxsYPBgeRPi7xO1RiPQpQuwcaP7r9dApKoqTp486bTPcf/+/bBcOeV09dVXOwXHzp07B/RMVUqKrJCUlGhbSshgkDde69Z5b6m7JiwWC3bs2IHvv/8e69evx5YtW2CxWNC6dWuHAz7cCkRaYaAMIFYrEBYmh2RcddVVEhpXrpTgOGaMYxmfHTtkKX3iROC116TwcFWeeQZ46y05zLNyZfXPHxQEdOnSB2lpu1wec1xcHHbs2OHy7cnRmTPAtdfK/1etyrOYTEB0tHQ0adRImzFoIT8/v8LSPGfOnAEANGjQoMLSPI0bN9Z45Nr4+mtpwmDrou1vRqPMqH//vWwN0pPLly/jhx9+sM9g7tq1C6qqolOnTrjuuuswfPhwJCYmokmTJloPleoIBsoAM3Cgc6ux6tSkJ+6gQcBf/gJ89RWQkSF74kaPBh56SELCiBGu/QHo3Ru45pop+Oc//wlXXoqKomDo0KH46KOP0KFDB/addlNJicxMbtumfUcSo7F05ifQVuxUVcWxY8ec9joeOnTIXpqnY8eODmV5YmNjERkZ6XJpnrriP/8Bbr9dfj/5c/nbYJA6vGvW1I43PWfPnsXGjRvtM5iHDx+Goijo3bu3PWAOGjQI9evX13qoFKAYKAPMtGnS4aa6moFlvfGGHJbYvt35FO7YsVIovbBQ9rzZZj+jooA5c6QmW7NmEh4PH5ayQe++C1xpolGloCDgwQeBefOAPXv24I477sD+/ftdCpaAFEpu3749oqKiYDabERUV5fB1XZ3Vqcp778lrRE/mzwfuv1/rUdRcbm5uhaV5cnNzAQDh4eEVluZpoLcpLx379lvg1ltl/7avZ9UVRX6fjRsnvb5ra/46fvy4wwGfkydPIjg4GAMGDLAHzH79+gX0tglfsVrl793hw9JVrqhIDm01aCDbeczmwHuT7AoGygCzeLEcnqktPvkEeOAB+bqoqAgvv/wyXn/9dSiKUukeyrfffhu9e/dGRkYG0tPT7Z/T09Pt9fUA+UNePmTaPrdt27bOtUI7dEiK0LsS9v2pQQPgwAFp46lnJSUlSE9Pd1quPnr0KIDS0jzlw2Pbtm15aMILTp6UFZBVq0pDn7eZTBIg//53aSkbKP/bVFXFgQMH7LOXGzduRG5uLq666ioMHTrUHjCjo6P5Wq2A1Sqrc+vWyepOaqoEycpcdZV0D+vfXzrLDRoUOK+lqjBQBpgTJ4AOHfTfD9fm8GHnMhzbtm3DnXfeiYyMDFgrWONKTk7GkCFDKny8c+fOOYTMsmEzKyvLPvsZHByMiIiICmc3zWZzwC0LWa2yb3LnTv21tTOZ5PCFnpa+z5496xQc9+3bZy/N06pVqwpL84SEhGg88sCmqsAXXwAPPyxbc7y1BG4yyb+LMWOAjz+W7mGBrLi4GDt37rQHzM2bN6OoqAgtW7bEsGHD7AGzQ4cOWg9VU+fOAQsXAh98AGRmlr5OXGW7fdeu0uBh0iQJm4GKgTIAjRkjS0R6Cw5lGY3AkCFSxLgi+fn5eOGFFzBnzhwoiuIQLLOysmrUqqywsBCZmZkOIbNs8Mwvc5qpVatWTrObtq9btGhR697Fr14t75T1LCVF3sn7U1FRUYWleXJycgAAISEhFZbmacHWTpo6dw5YsED+0B875v4fehuDQd7EjB8PPPIIMHSoft7U+FN+fj42b95sD5g7d+6EqqqIiorC8OHDcd111yExMRHNmjXTeqh+kZ8v9XLnzCndPuZJUrK9pkJDpYrK9Omy5SvQMFAGoLVrpYak3q1YIb/Iq5KcnIxJkyYhJycHVqsVwcHByM/P9/rBBVvplopmNzMyMnDq1Cn7bRs0aOC0hG77Wq8HhUaPBr77TvuDOJUxmaSSwBdf+ObxVVXFr7/+6tS/+sCBA/bSPB06dHDqX92pUyfuMdMxq1Ve1x9+KK1gr2xbhckkAcBqVaGqEhwNBsVegshgkJWRSZNkH3egz0i66/z58w4HfA4ePAgA6NWrlz1gDh48OCD3Af/4o7wuMjN9cwhMUaTCxaJFcgYhkDBQBiCrFYiMlPqOetWqlSzPu/K3Oi8vD0899RQ+/vhjdOzYEYcPH/b9AMu5dOkSMjIyKlxKz8zMRPGV6ZGyB4Uq2rupxUGhjAz546n3f+kmE5CVJf2/PfH7779XWJrn7NmzAICwsDCn0jw9evTgIa5aTlUlBOzcKR/HjwOnT+fhu++SMGzYAERGNkdsrOxt69VLf2WA9CwrKwsbNmywB8ycnBwEBQUhPj7eHjD79++PoFo87WaxSO3kd9+VNxy+fPNt274/c6ZUSwmUwg4MlAHqH/+Q9mV69be/AU8+6d59NmzYAEVRkJiY6JtB1VBJSQlOnDhR4VJ6+YNCTZo0qXR201cHhZ59Vv5763V20sZgkP7wrhbEt5XmKR8cDx8+7FSap+xHREQES/PUEYcPH0bnzp2xceNGDB06VOvhBARVVXHw4EGsX78e33//PZKSknDx4kWEhYVhyJAh9oDZo0ePWvPvLD9f6p2uW+ffN96KIk1FFiwIjCVwBsoAZbXKHsWfftLXXkqTSWYHtm4tfZcW6GpyUKj87KYnB4WiomSWsioTJgDLlwO33QYsW+Z43c8/y9LM9dfL3tyyjhwBLl6UWZ8xY6S0S+/esgk9KMj9/Wh9+sjsUnm5ubkONR1tX+ddKZoaHh6Onj17OgTH7t27B+SSHLkuIyMDUVFRWL9+PYYNG6b1cAJSSUkJUlNT7QFz8+bNKCgoQPPmzTFs2DB7wIyMjNR6qBUqKABGjZI93Fq0+VQUaQTy5ZeurdjpGQNlAEtPB3r0kH8wehEUBOzZI4GDnA8KlZ/lLHtQqHXr1pXOblZ2UCg3F2jcuPp33U2bAr/9Bnz0kePMdni4dNb5/Xepbzp9eul1bdvKEvW77wJPPy0loAYPBnbtkhDbt6/7gTIoSEVq6iEcOOA465iZmQkAMJlMFZbmadOmTa07KEW+d+zYMURERODbb7/FiBEjtB5OnVBQUIAtW7bYA+aOHTtgtVoRGRlpbw85bNgwrx1su3DhAubPn48pU6a4/aa7uFj28a9Zo23PeEUB7rpLZipryaRuhRgoA9wHH0i5Ar146y1pz0jVsx0UquxUemUHhcoGzjNnojFpUjuXnm/3biA4GOjWrfSyceOApUvl0EN8vHzYTJoEfPaZnB5fvdqxNuAHHwBTp9b0xGwvALvRunVrp04yLM1D7sjOzka7du2wdu1ajBo1Suvh1EkXLlxAcnKyPWAeOHAAABAbG2sPmEOGDMFVNaynM3v2bDz55JPo06cP1q5d61ZQfeUV2ceolxQ0b56+t6pVh4EywFmtwN13A59/ru0/GoMBuPlmqQd0ngAAIABJREFU6aRTV5a6fc21g0JPAngHQPVve2fPBh5/XE68njxZelnfvsDLL0toDA8HLl2S6z75BLjnHpndtJ2utal5oFTx1FO/4LnnmqF58+bu3pnIwcmTJ9G6dWusWrUKf9B73aw6IicnBxs2bLAHzKysLJhMJlxzzTX2gBkfH+9ytYwxY8Zg7dq1MBgMaNu2Lb777jt07ty52vvt2SNbdfS0JSw0FNi/Xw7V1koqBTyLRVXHjlVVRVFViZX+/TAYVHXkSFUtLNT6v0TdUVxcrB49elS94YZs1WQqcen/0003yX1vv730st27VfW111S1QQNVLSpS1dGjS69LT1fVn36q+LE++EAey93XSlCQqk6bpt1/Nwosv/32mwpA/frrr7UeClXAarWqBw8eVOfNm6dOmDBBDQ8PVwGo9evXV0eNGqW+/fbbampqqlpSUlLh/YuKitR69eqpAFQAqtFoVBs1aqRu3ry5yuctKlLVmBhVNZm0+ZtY2YfJpKpDhqhqJT+u7tXi1XpylclUeuDC39vMbBuOV66U5VTyD6PRiIiICDRq1AZWq2v/zJOT5SR4QoJ836SJ7MFNTpYeyqmpgO2Afbt20q82Kcn7Y6+qpRmRO2xVEypr40raUhQFnTt3xpQpU7B8+XKcPn0aO3bswMyZM2G1WjFjxgz06dMHLVq0wK233op//vOfOHLkiP0g4/bt2x32mZeUlCAvLw+JiYlYsWJFpc/77rtAWpq+ZicBGc+mTcD8+VqPpGYYKOuIoCDp8/3ii7Lk7OtlZ6NRlrlfeEFODXPbmzYKC13f6nDhguyjtAXKoUMlYG7eLN8nJ5cGSttnbwdKVdXXITKq3Rgoaxej0Yi4uDg8++yz+N///ofz588jKSkJU6ZMQVZWFh555BF06tQJEREReOCBB/Dee+85lVqzWq2wWCyYOHEi5syZ4/QcRUXA22/rZ99keYoCvP66toeEaoqBsg4xGIAZM6QsS7duvputVBQppP3TT8Crr9b+Ugi1mbule5KSgC5dZB9lYqK8Vi5fluuSk6UkUMOGcp3FAvzwg3fHqyicySbvsXU5YqCsnUJCQpCQkIBXXnkFW7Zswblz57By5UqMHz8e27Ztw/Llyyv8f6uqKlRVxRNPPIEnnnjCoXXvihXSulOvVBU4erTytsR6xkBZB/XsKcuXL74I1K8vf8Q9DZe2xwgNlVnJ3bvlMAdpq14998pQ2GYcExLkIzm59DpbeBwyRK7bvr00bHpTvXref0yqm2yzV8V6W9ukGmnYsCFuvPFGzJ49G1u3bnWpEcT777+PiRMn2pfGP/hA/wdDTSYp01bbMFDWUUFBMlt58qSUKrDVhTSZ3AuXttnHTp3kH+qpUzIrySVufYiKcm9pZ9Mm2cczcaL0m924sfS63Fwpcn7PPXIK0Rf7J4uLZcxE3sAl78D1ww8/uPT/VVVVfPXVV+jSpQs2bTqDLVv03zWsuFjOHeTkaD0S9zBQ1nFXXQX8+c/Avn3SKeCZZ4Bhw2RZs7r7JSYCTz0lweKXX4BHHqn+fuRfffu698szL09mr8eNkz08tv2TNsnJUv4JcA6UV18tHXcmTCgNhbbv4+Jce35Vdf22RNVhoAxc69evt29pqIyiKPb2jydOnMDcuWnVTphMmCC/h2691fm6n3+W60aOdL7uyJHSLl9BQcBLL0mHssJC6TE/a5as4LlKVb2/pcjXuLuNAMis5KBB8gHIi/nYMeDgQTl1W1gos4716gGdO8sMFRuT6F9NwllSEtC/v/xyvNLZ0C45GZg2TV4PW7Y4XpeYCCxc6HjZ8uXyeeFC4L77qn9ug0FacxJ5gy1MMFAGnnXr1jltZQgKCoLFYgEANGnSBPHx8YiPj0e/fv3Qr18/PPtsUxiNVZ/u3rhR3kwnJjq2oQ0PB2JipA5vYqJjG9q2beVN9LvvyvdffAHccIPU792+HRgwAPjrX2XV56abXPv5goLkd3BFwVavGCipQooCRETIB9VezZsDrVqVFip3xfPPy0dFVq6s/I3Ep5/Khyc6dZJ9vUTeYjKZGCgDTF5eHtLS0uzf16tXD3FxcRgwYAD69++P/v37o3379k7tWLdurb5U0NmzUlLIVu3CZuhQue/8+aVVLmzKVr245hqZ5Zw2DXjvPbl8/Xq57+uvA9ddB3z/ffU/o8UCbNtW/e30hIGSKMCNHi0tEvV+LsFolHf1RN5kNBp5KCfAhIWF4c0330SzZs3Qv39/dO3atdoDOgUFsuLmiqQk6RpW9s247SDimjWyvSssrLRrWEKC/H5NSQEefFAuW7PG8TG/+UYC5YQJrgVKQGYoVbX2rAZyDyVRgHv4Yf2HSUD2ev75z1qPggKN0WjkDGWAURQFzzzzDO677z5ER0e7dNo7O9v1/eRlq13YJCbKlp/NmyXkDR7seF1qqhxctJU9Kyx0fEzb97Gxro0BkC1H5dva6hkDJVGA69tX6ke6Uz7I34xGOQzmQgteIrcwUBLgXsMET7qG7d8vnwcOdHxM2/mEpk19N26t6fhPDBF5y2OP6bvzQkkJMHWq1qOgQMRASYB7qzSedA1buxY4fBh4803ZL9moEXD99XLKu7jY/d/DV84Y1QoMlER1wG23SVkfPRb0NRqlDuqNN2o9EgpEJpOJeyjJ7drINe0aZrHIvvXjx4HvvpNwuny5BMrz52Xp3R21qdEDAyVRHVCvnpzA1uNEjaoCixaxRSf5BmcoCZDaye7wpGtYejpw7bVSTigmBmjRAvjyS6m6sWmTe+OoTVUvGCiJ6oiEBFlW1tNeSkWRVp0sZk6+wkBJANCmjXuh0htdw3JypARRfr40Dbl0ScoOuSoionbNUHJOgKgOeeMNYNUqICtL+9lKo1EO4fzf/2k7DgpsDJQEyJvXfv2ADRtcu70rXcOeeEK+Lh8on3lGyg0dPw60bCnFyceNA+66y/V2ikYjEB/v2m31QkdzFUTkaw0aAKtXSw01LfdTGo3SeWLVqtIyG0S+wMLmZNOvn3SgcVVSkqzo7NpVcdcwg6HirmGhocCMGcC6dcA//iHL1gkJ0kHHVbWxDa2iqqqq9SCIyL+2bZMyPQUF/p+pNBpl6SklRUpxEPlSp06dcPPNN+PNN9/UeiikseXLgVtu0XoUrtuwwbkrj55xhpKoDurfX/YINWzo38MwJhPQrJm8o2eYJH/gkjfZjBxZe/YkNmvmXMtS7xgoieqoPn2AHTuk96y/JCTIc3br5r/npLqNgZJsGjYE7r1X/xUljEZgypTatx2IgZKoDjObZaby/fdl348vftGaTLJ385NPgG+/lc4SRP7CQEllTZmi/1a0qgo89JDWo3AfAyVRHWcwAI8+KuUtrr9eTkN648CO0SiPfdNNwIEDwAMPyGMT+RMP5VBZMTFSI1KPTR4AeQN+44218403AyURAQCiooBvvpGivE8/DTRuLJe7OmupKKW3bdZM6kseOyYb4du3982YiapjNBrZKYccvPOOvlvRzpql9QhqRuc7CYjI3yIjpV7liy8Ca9bIifCffpLWY+VLZ9iEhORj8OB66N9faqddf33t2/9DgYlL3lTegAHAtGnAe+/pK1gqCvDKK0D37lqPpGZYNoiIXKKqQGYmcPasdH4wGGTf5bPP3of69c9i1aqVWg+RyMnAgQPRuXNnLFiwQOuhkI7k58vyd2am9k0eAFnd6dkT2LpV/4eGKlNLh01E/qYoMnsZGel4effuYUhK2q7NoIiqwRlKqki9esCiRcDgwfK7TcupNYNBQuS//117wyTAPZRE5CGz2YyMjAxwsYP0yGQycQ8lVSg+Hvj8c23HoCgSKFetqr1L3TYMlETkEbPZjPz8fJw6dUrroRA54QwlVeWWW4CFCyXY+bsKhcEgp83/8x/guuv8+9y+wEBJRB4xm80AgKNHj2o8EiJnDJRUnbvvllBnMvlvydlkAkJCpN/32LH+eU5fY6AkIo9EXtlUmZGRofFIiJwxUJIrxo8Hdu+WgzH+mKkcOBDYvx8YPtz3z+UvDJRE5JGwsDC0aNGCgZJ0iYXNyVXduskp6zfeAIKCvD9baTJJZYx584ANG4CICO8+vtYYKInIY5GRkQyUpEssbE7uMJmAZ58F9u6VpfCQEM/2V9ruW78+MHmydA2bMkX2TwaaAPyRiMjfbCe9ifSGS95UE126APPnA7/+Cvztb47l0oKCKr9f2Y5hANC1q8xInjwJzJ0beLOSZdXiikdEpBdmsxkpKSlaD4PICQMleSI8HHjiCeDxx6UI+s6dwI4dwPbtMttYUABYLBIy69UDevQA+vUD4uLko317/58e1woDJRF5zGw2Izs7G4WFhQgJCdF6OER2DJTkDWUbO0ycqPVo9IlL3kTkMbPZDFVVcezYMa2HQuSAh3KI/IOBkog8xtJBpFc8lEPkHwyUROSxdu3awWQyMVCS7nDJm8g/GCiJyGNGoxEREREMlKQ7DJRE/sFASURewdJBpEcMlET+wUBJRF5hNpvZz5t0x2QycQ8lkR8wUBKRV9i65aiqqvVQiOw4Q0nkHwyUROQVZrMZubm5OHfunNZDIbJjoCTyDwZKIvIKs9kMgKWDSF8YKIn8g4GSiLyCgZL0iIGSyD8YKInIKxo3bozw8HAGStIVHsoh8g8GSiLymsjISJ70Jl3hDCWRfzBQEpHXsBYl6Q0DJZF/MFASkdcwUJLeMFAS+QcDJRF5jdlsxvHjx2GxWLQeChEACZTcQ0nkewyUROQ1ZrMZJSUlOHHihNZDIQIgh3I4Q0nkewyUROQ1LB1EesMlbyL/YKAkIq+5+uqrYTAYeNKbdIOBksg/GCiJyGuCgoLQvn17zlCSbjBQEvkHAyUReRVPepOe2AKlqqpaD4UooDFQEpFXMVCSnphMJgCA1WrVeCREgY2Bkoi8ioGS9MRoNAIAl72JfIyBkoi8KjIyEufOncOFCxe0HgoRAyWRnzBQEpFX2UoH8aQ36YEtULK4OZFvMVASkVcxUJKecIaSyD8YKInIq5o1a4awsDDuoyRdsB3KYaAk8i2T1gMgosCiKAoP5pAmli9fjo0bN6KkpATFxcUoKSlBZmYmAODBBx+0t2G0Xf/www9jzJgx2g6aKEAwUBKR10VGRiI9PV3rYVAds3DhQqxevRpBQUH2upO2ckErV66EqqoO9SgTEhIYKIm8hEveROR1ZrOZeyjJ755//nkAgMViQXFxMYqLi+2B0mq1OoTJsLAwPPTQQ5qMkygQMVASkdeZzWYcO3aM+9bIrwYNGoQRI0bYD+JUxmAw4Mknn0SjRo38NDKiwMdASUReZzabUVRUhJycHK2HQnXMrFmzqn0jExwcjMcff9xPIyKqGxgoicjrbKWDeDCH/K1v374YO3as/XR3eUajEVOnTkXTpk39PDKiwMZASUReFxERAYCBkrTx6quvVjpLaTQa8dRTT/l5RESBj4GSiLwuNDQUbdq0YaAkTcTExOCWW25xmqU0Go146KGH0KpVK41GRhS4GCiJyCdYi5K09PLLLzvNUiqKgmeffVajEREFNgZKIvIJlg4iLXXp0gV33323/cS30WjEvffei/bt22s8MqLAxEBJRD7BGUrS2syZM+21J1VVxfTp0zUeEVHgYqAkIp8wm804deoULl++rPVQqI6KjIzEhAkTAAAjRoywVx8gIu9joCQin4iMjAQALnuTpl555RU0bdoUL730ktZDIQpo7OVNRD5RthZljx49NB4N1QW//Qbs3Ckf27fL59zcLigqOoOEBCA0FGjdGrjmGqBvXyAuDujZE6hXT+uRE9V+DJRE5BOtWrVCaGgoZyjJp4qLgW++Ad5/H0hKksuMRkBVgSttvO0KCoALF4DDh4FPP5XbhIYC994LPPwwEBPj9+ETBQwueRORTxgMBkRGRvJgDvnEpUvAa68B7dsD48cDmzaVXldS4hwmyyouljAJSMj85BMgNha49lpgxYrS64jIdQyUROQzPOlNvpCcDHTvDsyYAZw8KZdV0767SsXF8nnbNmDCBOCmm0ofl4hcw0BJRD7DGUrypsuXgalTgYQEIDu76lnImrCF0rVrga5dgc8/52wlkasYKInIZ2zFzVX+VSYPZWUBvXsDH34o33s7TJZVXAzk5gJ33glMnuzZ7CdRXcFASUQ+YzabkZ+fj1OnTmk9FKrFMjKA+Hjg6FHfBsmybO+BPvkE+OMfS5fFiahiDJRE5DNlSwcR1UR2tixxnzypTahTVWD5cuC++/wXZolqIwZKIvIZW3FzBkqqicJCYORI4NdftV12VlVg0SJg5kztxkCkdwyUROQzYWFhaN68OQMl1chLLwEHDuhnufm116RgOhE5Y6AkIp9i6SCqie3bgTfe0Ncpa4MBuOsumTklIkcMlETkU7aT3kSuKiwEJk2SAKcnJSXAoUMyc0pEjnT2z5WIAg1nKMld//63BDc9lutRVeDttwEWLiByxEBJRD5lNpuRnZ2NgoICrYdCtYCqSl9uRdF6JJWzWoF//UvrURDpCwMlEfmU2WyGqqo4duwYAODixYvYtWsXVqxYgR07dmg8OtKbrVuBtDR97Z0sz2oF5s7V5wwqkVZMWg+AiAKLxWLB8ePHcfToUWRkZGDXrl0AgNGjR+Ps2bPIzc2133bQoEFISUnRaqikQ3//O2Ay6edkd2VycoA1a4Abb9R6JET6wBlKIvKaJUuWIDQ0FB07dsSIESMwefJkzJ8/HwBw9OhRhzCpKApGjx6t1VBJh0pKgBUrHMPkhAkyW3nrrc63//lnuW7kSOfrjhwBdu6Ur8eMAT79FNizBygqqnr202QCZsyQrjwFBVK2aOrUim+3fLl7Px9RIGOgJCKv6dSpk9NlFoulwtuqqorbbrvN10OiWuTgQSA/3/GyjRtliTkx0fHy8HAgJga4dMn5urZtgagoIClJvh8/Xlo37t8P7N5d9RjmzQOmT5eZ0uuvB776CpgzRy4rq7hYlueJSDBQEpHXxMXFYfr06VCqOVGhKAr69OmDqKgoP42MagPbjGJZZ8/KnsqEBMfLhw6VUDd/vnOgtH1vC5R/+hPQpQtw++1Vh8Du3YEHHpCyQO+8AyQnAy+8IP28//pXCbFlHTkCXL7s1o9IFLAYKInIq2bMmIHu3bvDaDRWehtFUTBp0iQ/jopqg507gaAg58uTkoCuXYFWrUovS0iQ4udr1gBxcUBYmON1xcWAbXuuqwd8xo2T2pcLFjhevmABUL8+MGqU4+VWa/UznkR1BQMlEXlVcHAwFi9eXOUspaqquLWiTXFUp+3eDVS0Q8I201h2ljIxUWYQN2+WwDh4sON1qalAmS27LunRA/jtN+cak3v2lF5flqIAe/e69xxEgYqBkoi8rmfPnnjxxRcrDJUGgwEDBw5E27ZtNRgZ6dnFixVfnpwsB3ZsgbJJEwl3ycmy5JyaWrrM3a4dYDaXhlB3NG0KnDvnfPnvv0v3nqZNHS83GGQPJxExUBKRjzz33HPo3bu309K3qqq48847NRoV+UtKSgpSUlJQ7Eb9n8p6ZF+4ILOXtkA5dKgEzM2b5fvk5NJAWX7/pLuqWh4vf52iyElwImKgJCIfMZlMWLRoEQzlGjIrioIJEyZoNCryl7vuugtDhgxB06ZNcccdd+CLL77AuYqm/8owVVEZOSlJDta0bi2hcefO0gMxyclA795Aw4ZyncUC/PCD+2M+e9Z5FhKQ/ZMhIRXPXla055OoLmKgJCKf6datG2bNmmX/XlEUDBs2DM2bN9dwVOQPthJSubm5+PLLL3HHHXegWbNmGDhwIN566y3s378farkpv/r1K3+8svsoExIkRNrYwuOQIaWHdWpy+nrvXqBFC6BlS8fLY2Lkc1qa4+VWKxAa6v7zEAUiBkoi8qknn3wS8fHxAGS5m6e764ZevXoh6Mr0nW3ZW1VV/Pjjj5g+fTqio6Nx9dVX49FHH8X//vc/FBQUoGNHoLLiAJs2ycntiROB6GipT2mTmytFzu+5B4iMrPly93//KyHxnnscL7/3XtlHuW6d4+VWqzwfEbH1IhH5mNFoxGeffYYuXbpAVVWMGzdO6yGRH8TExFRY1F5VVfvMZFZWFubNm4e5c+dCURTcc88uqGrPCh8vL08O34wbJ0HOtn/SJjkZeOIJ+bp8oLz6aqBfP/naVvrUtusiM7O0/uX+/VLX8qWXZI/m9u3Sheehh6QO5fnzzuOKi6vuvwRR3cBASUQ+17FjR/zlL/+HHTsuIDW1EfLz5Q92vXpAo0Yy41TVcifVPrGxsS7dzmq1AgBatWqFG25oiYULK79tUhLQv78EwLw8x+uSk4Fp0+Rgz5YtjtclJsLpcW1tExcuBO67r/Tyhx8GsrOBRx+VupeZmcDjjwNz5zqPp1kzoE0bF35IojpAUctvYiEi8gKLBVi9WopLb9sms0u//17xbQ0GoGNHYMAAoG9fmT1q3dq/4yXPXL58GWlpadi7dy/27NmD3bt3Y9OmTVXeR1EUqKqKp59+GrNmzUJhYRAaNnS9ELmWFAUYPVpe40TEQElEXpaVBXz8MfDhh8Dp03IKtpJ23k5MJlnOVBQJlY88IgWrq+nkSH5UUlKC9PR0e3C0fc7IyICqqjAYDOjUqRNiY2Oxfv36Sk92m0wmhIWFYdGiRRgzZoz98r59gV275HWgZwYD8NZbwFNPaT0SIn1goCQirzh1SvawLVsmAbCkxLPHM5nkEEbXrrLcOHy4d8ZJrjt9+rRTcNy3bx/y8/MBAC1atEBsbCxiYmLsn7t374569eoBAG6//XYsX74cJeVeDIqiID4+HsuWLUO7du0criu/BK1XQUHAr79WXGaIqC5ioCQij6gqsHQpMGWK7GvzNEiWZzDIbNXkycDbbwNXXeXdxyegoKAABw4ccAiOe/fuxcmTJwEAoaGhiI6OdgiOMTExaFm+vk45s2bNwowZM+yB0mAwQFVVTJ8+HS+99BJMFRSezM+XvYvutk30J5MJuOMO4NNPtR4JkX4wUBJRjZ05A/zpT8DXX8uspC9/mxgMsq9y0SLHns7kOlVVcezYMafgeOjQIXvoM5vNDsExNjYWHTt2dOp45IpVq1Zh7NixAOS0f6NGjbBkyRKMGDGiyvs9/TQwe7b335x4008/yQEhIhIMlERUIydOyOnZzEz//eE3GCS4fvYZ8Mc/+uc5a6sLFy5g7969DsFx7969yLtyPDo8PNwpOEZHR+MqL04BZ2ZmIvJKocahQ4fiiy++QGsXTltlZspWh8paMWrJZALi46UuJvf2EpVioCQitx0/Dlx7reybdKNVs1fY/ogvWOBcgLouslgsOHjwoNNexxMnTgAAgoKC0LVrV6e9jm3btoXi40SkqipGjRqFQYMG4YUXXnBrlnPOnNK6knoSEiIdczp21HokRPrCQElEbvntNynvc/y4/8NkWYoitQRvvlm7MfiTqqrIyclxCo4HDhywFxBv166dU3Ds0qULgoODNR69+6xWOeG/bZu2r7PyZs+WupRE5IiBkohcpqrSOWTjRu3/yCsKEBwcmLNFly5dwr59+5z2OtpK8ISFhdkPxpQ9JBMeHq7xyL3ryBGgRw+gqEj72pQmk+yZTEmRrRdE5IiBkohc9vHH0oZOL2r7H3lbTcfywTE9PR2AnIru3Lmz017HDh06wFAbf+Aa+M9/gFtu0TZQmkxAy5bA1q1AuSpHRHQFAyURueT4caBbt8q73Whpzhzgsce0HkXVTp8+7RQcy9Z0bNmypVNw7Natm72mY122YAHwwAPahEqjEWjSRNo5BtpMOJE3MVASkUtGjpReylovdVckJAT45RcgIkLrkUhNx/379zvtdTx16hSA0pqO5fc6tmjRQuOR69vixcDdd8vX/uqiYzJJqaqkJCAqyj/PSVRbMVASUbV+/hno3VvrUVTOaJSDEu++67/nVFUVmZmZTsHx8OHDDjUdywfHmtZ0JGD9euCuu+RgmD9KVY0cKZ172FeeqHoMlERUrcmTgX/9S5+zkzYNG0orvPr1vf/YtpqOZYNjWlqaQ03H8sHR2zUdSVy8KP2z588v7aLkTSaTzHh/8AFw772sNUnkKgZKIqrSxYvSCq+gQOuRVO9f//KsD7StpmP5vY5V1XSMjY1FmzZtfF7TkRx9+y3wyCNyEtzW990TRqOE07FjpXc8D98QuYeBkoiq9Pe/A48+qn3ZluoYDEDPnkBqavW3tdV0LB8cy9Z0bN++vVNZntpa0zFQqarsb5w7F/jvf+UyRXF9OTwoCLBYgEaNZBZ+8mTAbPbdeIkCGQMlEVVp/Hhg5cqqlxYnTJAi47fdBixb5njdzz9L0Lv+eplVKuvIEZkBjYsDXn0VuOEGoEMHWbbOyQG+/x547TU5Ye6qCxckINiUrelYNkCeP38eQGlNx7Kzjj169Ai4mo6BLjtbXns7dgA//ggcPVr5bYODgZgYaaE4cKC8xkND/TdWokDEQElEVWrTRvYmVqVpUzko8dFHwJQppZeHhwNnzkipoblzgenTS69r2xbIypKDNE8/LdcfOwYcOADk5QHduwN//avMPEZHA1dqelfrxReTYbVuYE3HOi43F9i7V95gFBTIknZoqLyeu3eXZXIi8h4GSiKq1OnTgKvVbHbvlpmfbt1KLxs3Dli6FPjwQ5kNio8vvW7SJOCzz4A//AFYvbrixxw1Cli7Frj/fqlFWD0rgGfQsuVi1nQkIvIjvkcjokrt3On6bZOSpHRPq1bAyZNyWUICsH07sGaNHKAICwMuXSq9rrhYutxU5vRp+ezqgQujUcGNN76Gr77yY/0gIiIC13mIqFJHj7peNiUpST4nJJRelpgIJCcDmzfLAYrBgx2vS02VpcmybEuTvXoBs2cDBw8CK1a4NoaSEgUnTnAzHBGRvzFQElGlCgpc75GdnCyna22BskkToEcPufzyZQmPiYlyXbsSKaeOAAAEj0lEQVR2cprWFkJtWraU2cj8fGDXLtnnlpgo93fVlU6GRETkRwyURFQpi8X1GcoLF2QfpS1QDh0qAXPzZvk+Obk0UNo+lw+UZ84AffvKydsHH5RQmpQky+iuKix0/bZEROQdDJREVKmQEPc6kSQlAV26SKu6xETZg2mbXUxOlvaNDRvKdRYL8MMPjvcvKZH7bNkinVCGDZOZzOefd30MPHdDROR/DJREVKmwMPcDJSCzlAkJEiJtbOFxyJDSwzrVLWVnZ0s9ys6dXR9Dw4au35aIiLyDgZKIKhUd7d7tN22SPZATJ8p9N24svS43V4qc33MPEBnpvNxdkago2W955Ihrzx8UBMTGujdmIiLyHMsGEVGlevaUQzmuzlLm5cnhm3Hj5D62/ZM2ycnAE0/I12UDZUwM8N570m0nI0PuGxMDPPkkcPYs8M47rj2/xSJdd4iIyL84Q0lElapXz73lZkCCosEgp7Tz8hyvS06W6woLZZ+kzalTsrT91FNSIuibb6Ru5TffSEB0p/UiAyURkf+xUw4RVen++6WjjavFxbUUFCSF04ODtR4JEVHdwhlKIqrS+PG1I0yaTMCNNzJMEhFpgTOURFSlkhKgQwc5ca13GzaU1rgkIiL/4QwlEVXJaASmTnW9Y44WFAXo2NGx7SMREfmPjv9EEJFe3H+/vgMlADz6qOtdfYiIyLt0/ieCiPSgRQvg6af1GdiMRlmSf/BBrUdCRFR3cQ8lEbmksFCKhqeny75KvVAUKUc0eLDWIyEiqrs4Q0lELgkJARYtAvT0FtRgAB57jGGSiEhrDJRE5LJ+/YAXXtB6FMJkAiIigFmztB4JERFxyZuI3GK1AnffDXz+uXazlUYj0KwZsHWrhEoiItIWZyiJyC0GA7BggRQ81+KQjskENGkiLR4ZJomI9IGBkojcFhQELF0K3HuvfO+vYGk0Am3ayMxkt27+eU4iIqoeAyUR1YjJBMyfD3z6KRAWJt/7iq0G5t13A7t3A2az756LiIjcxz2UROSxnBypA7l2rYQ/q9V7j23bL7lgATB6tPcel4iIvIczlETksTZtgNWrZRm8Tx+5zJMZS6NRPoeHA889B/zyC8MkEZGecYaSiLwuNRWYN0/qVhYWyqylolReEF1RJEQWF8v3AwZIK8Wbb5b6l0REpG8MlETkM5cuAbt2ATt2ADt3Aj/+CJw8KSHTapWw2KABEBcnNS7j4oC+fYH27bUeORERuYOBkog0oar67A1ORETu4x5KItIEwyQRUeBgoCQiIiIijzBQEhEREZFHGCiJiIiIyCMMlERERETkEQZKIiIiIvIIAyUREREReYSBkoiIiIg8wkBJRERERB5hoCQiIiIijzBQEhEREZFHGCiJiIiIyCMMlERERETkEQZKIiIiIvIIAyUREREReYSBkoiIiIg8wkBJRERERB5hoCQiIiIijzBQEhEREZFHGCiJiIiIyCMMlERERETkEQZKIiIiIvIIAyUREREReYSBkoiIiIg8wkBJRERERB5hoCQiIiIijzBQEhEREZFHGCiJiIiIyCP/DwzVSGlHyr6lAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "draw_graph(graph_true, highlight_nodes={\"X\", \"Y\"})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 226
    },
    "id": "xoFVii0MdmFC",
    "outputId": "dc55c177-95b0-4a53-9312-3e2af0e5b44c"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>W1</th>\n",
       "      <th>W2</th>\n",
       "      <th>W3</th>\n",
       "      <th>W4</th>\n",
       "      <th>W5</th>\n",
       "      <th>W6</th>\n",
       "      <th>W7</th>\n",
       "      <th>W8</th>\n",
       "      <th>W9</th>\n",
       "      <th>W10</th>\n",
       "      <th>X</th>\n",
       "      <th>M1</th>\n",
       "      <th>M2</th>\n",
       "      <th>M3</th>\n",
       "      <th>Y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-0.524526</td>\n",
       "      <td>0.082740</td>\n",
       "      <td>0.511285</td>\n",
       "      <td>-0.428790</td>\n",
       "      <td>0.206569</td>\n",
       "      <td>-0.124775</td>\n",
       "      <td>-0.274969</td>\n",
       "      <td>0.975348</td>\n",
       "      <td>0.223864</td>\n",
       "      <td>-0.261680</td>\n",
       "      <td>-0.238098</td>\n",
       "      <td>0.389241</td>\n",
       "      <td>-0.567587</td>\n",
       "      <td>-0.837017</td>\n",
       "      <td>0.502616</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.292501</td>\n",
       "      <td>0.016717</td>\n",
       "      <td>0.320840</td>\n",
       "      <td>-0.879261</td>\n",
       "      <td>-0.427722</td>\n",
       "      <td>0.035111</td>\n",
       "      <td>-0.687620</td>\n",
       "      <td>0.801554</td>\n",
       "      <td>-0.022219</td>\n",
       "      <td>0.166778</td>\n",
       "      <td>0.362314</td>\n",
       "      <td>-0.447789</td>\n",
       "      <td>0.124390</td>\n",
       "      <td>-0.199582</td>\n",
       "      <td>-0.220629</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-0.648133</td>\n",
       "      <td>0.562219</td>\n",
       "      <td>0.385289</td>\n",
       "      <td>-0.137838</td>\n",
       "      <td>0.348370</td>\n",
       "      <td>0.444976</td>\n",
       "      <td>-0.771743</td>\n",
       "      <td>0.694962</td>\n",
       "      <td>-0.790796</td>\n",
       "      <td>0.027081</td>\n",
       "      <td>0.493028</td>\n",
       "      <td>-0.227073</td>\n",
       "      <td>-0.071652</td>\n",
       "      <td>-0.614490</td>\n",
       "      <td>0.402219</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.344170</td>\n",
       "      <td>-0.317331</td>\n",
       "      <td>-0.397835</td>\n",
       "      <td>-0.059340</td>\n",
       "      <td>0.439778</td>\n",
       "      <td>-0.009194</td>\n",
       "      <td>-0.915783</td>\n",
       "      <td>0.124095</td>\n",
       "      <td>-0.216393</td>\n",
       "      <td>0.328786</td>\n",
       "      <td>0.710924</td>\n",
       "      <td>-0.776587</td>\n",
       "      <td>0.095804</td>\n",
       "      <td>0.129171</td>\n",
       "      <td>0.715350</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.534461</td>\n",
       "      <td>0.242387</td>\n",
       "      <td>-0.178900</td>\n",
       "      <td>-0.957031</td>\n",
       "      <td>-0.057134</td>\n",
       "      <td>-1.591982</td>\n",
       "      <td>2.103615</td>\n",
       "      <td>-2.068360</td>\n",
       "      <td>0.725326</td>\n",
       "      <td>0.228350</td>\n",
       "      <td>-1.610583</td>\n",
       "      <td>2.057864</td>\n",
       "      <td>0.274302</td>\n",
       "      <td>1.125149</td>\n",
       "      <td>-0.491953</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         W1        W2        W3        W4        W5        W6        W7  \\\n",
       "0 -0.524526  0.082740  0.511285 -0.428790  0.206569 -0.124775 -0.274969   \n",
       "1  0.292501  0.016717  0.320840 -0.879261 -0.427722  0.035111 -0.687620   \n",
       "2 -0.648133  0.562219  0.385289 -0.137838  0.348370  0.444976 -0.771743   \n",
       "3  0.344170 -0.317331 -0.397835 -0.059340  0.439778 -0.009194 -0.915783   \n",
       "4  0.534461  0.242387 -0.178900 -0.957031 -0.057134 -1.591982  2.103615   \n",
       "\n",
       "         W8        W9       W10         X        M1        M2        M3  \\\n",
       "0  0.975348  0.223864 -0.261680 -0.238098  0.389241 -0.567587 -0.837017   \n",
       "1  0.801554 -0.022219  0.166778  0.362314 -0.447789  0.124390 -0.199582   \n",
       "2  0.694962 -0.790796  0.027081  0.493028 -0.227073 -0.071652 -0.614490   \n",
       "3  0.124095 -0.216393  0.328786  0.710924 -0.776587  0.095804  0.129171   \n",
       "4 -2.068360  0.725326  0.228350 -1.610583  2.057864  0.274302  1.125149   \n",
       "\n",
       "          Y  \n",
       "0  0.502616  \n",
       "1 -0.220629  \n",
       "2  0.402219  \n",
       "3  0.715350  \n",
       "4 -0.491953  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = get_dataset(mixing_matrix=A, noise_vars=noise_vars,\n",
    "                   num_samples=1000, labels=labels)\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "_REfaO6ID9N5"
   },
   "outputs": [],
   "source": [
    "sample_cov = data.cov().values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "-Q_8Qc_jdsm6",
    "outputId": "8378d8dd-33c5-4dcb-b304-c454971a1d3d"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total CI tests done: 4341\n",
      "Estimated ATE set: {0.053902203097293586}\n"
     ]
    }
   ],
   "source": [
    "pc_alg = PCAlgorithm(use_ci_oracle=True, graph_true=graph_true)\n",
    "cpdag_pc, _, _ = pc_alg.run(data)\n",
    "print(\"Total CI tests done: %d\" % pc_alg.ci_test_calls[\"total\"])\n",
    "print(\"Estimated ATE set: %s\" % str(\n",
    "    set(get_ATE_using_cov(g, sample_cov) for g in generate_local_dags_from(cpdag_pc))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "ocEqP4jDejM1",
    "outputId": "c9c1dc15-d5e3-4dd2-fa1b-2e62b944e6bd"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total CI tests done: 1466\n",
      "Estimated ATE set: {0.053902203097293586}\n"
     ]
    }
   ],
   "source": [
    "mb_by_mb_alg = MBbyMBAlgorithm(use_ci_oracle=True, graph_true=graph_true)\n",
    "result_mb_by_mb = mb_by_mb_alg.run(data)\n",
    "print(\"Total CI tests done: %d\" % mb_by_mb_alg.ci_test_calls[\"total\"])\n",
    "print(\"Estimated ATE set: %s\" % str(\n",
    "    set(get_ATE_using_nodes_and_cov(list(data.columns),\n",
    "                                    list(result_mb_by_mb[\"tmt_parents\"]) + par,\n",
    "                                    sample_cov)\n",
    "    for par in get_all_combinations(result_mb_by_mb[\"unoriented\"], \n",
    "                                    result_mb_by_mb[\"non_colliders\"]))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "AhlAd4GUjgA5",
    "outputId": "c9973d1c-b023-4d57-93e4-9dbf2c8e6d5f"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total CI tests done: 2711\n",
      "Estimated ATE set: {0.053902203097293586}\n"
     ]
    }
   ],
   "source": [
    "sd_alg = SequentialDiscoveryAlgorithm(use_ci_oracle=True, graph_true=graph_true)\n",
    "result_sd = sd_alg.run(data)\n",
    "print(\"Total CI tests done: %d\" % sd_alg.ci_test_calls[\"total\"])\n",
    "print(\"Estimated ATE set: %s\" % str(\n",
    "    set(get_ATE_using_nodes_and_cov(list(data.columns),\n",
    "                                    list(result_sd[\"tmt_parents\"]) + par,\n",
    "                                    sample_cov)\n",
    "    for par in get_all_combinations(result_sd[\"unoriented\"], \n",
    "                                    result_sd[\"non_colliders\"]))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "TV5TXUenuTbq",
    "outputId": "aff14a65-8d6d-4ebc-e57f-25dec4874599"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total CI tests done: 584\n",
      "Estimated ATE set: {0.053902203097293586}\n"
     ]
    }
   ],
   "source": [
    "ldecc_alg = LDECCAlgorithm(use_ci_oracle=True, graph_true=graph_true, ldecc_do_checks=True)\n",
    "result_ldecc = ldecc_alg.run(data)\n",
    "print(\"Total CI tests done: %d\" % ldecc_alg.ci_test_calls[\"total\"])\n",
    "print(\"Estimated ATE set: %s\" % str(\n",
    "    set(get_ATE_using_nodes_and_cov(list(data.columns),\n",
    "                                    list(result_ldecc[\"tmt_parents\"]) + par,\n",
    "                                    sample_cov)\n",
    "    for par in get_all_combinations(result_ldecc[\"unoriented\"], \n",
    "                                    result_ldecc[\"non_colliders\"]))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "3gjXpu5JNihj"
   },
   "source": [
    "The code below evaluates the various algorithms on randomly generated linear graphs.\n",
    "We use `ipyparallel` to execute the runs in parallel. The following command starts the `ipyparallel` engines:\n",
    "\n",
    "```\n",
    "ipcluster start -n <num_engines>\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "R74t7mjIfRDj"
   },
   "outputs": [],
   "source": [
    "import ipyparallel as ipp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "dsMrRbQdhDnY"
   },
   "outputs": [],
   "source": [
    "# Verify that ipcluster is running and import the necessary Python packages.\n",
    "\n",
    "parallel_client = ipp.Client(debug=False)\n",
    "dview = parallel_client[:]\n",
    "# Execute an identity map in parallel.\n",
    "ar = dview.map(lambda x: x, (i for i in range(0, 20000, 2)))\n",
    "assert ar.get()[0] == 0\n",
    "\n",
    "# Import the required Python packages.\n",
    "with dview.sync_imports():\n",
    "  import numpy\n",
    "  import os\n",
    "  import pandas\n",
    "  import networkx\n",
    "  from itertools import combinations, permutations\n",
    "  import logging\n",
    "  from causallearn.utils.cit import CIT\n",
    "  import collections\n",
    "  import copy\n",
    "  import ipyparallel\n",
    "\n",
    "  from causal_discovery.utils import orient_colliders, get_dataset, get_connected_component_with_node, remove_directed_only_edges, apply_meek_rules, adj_matrix_to_nx_graph, dag_to_cpdag, get_neighbor_types, get_all_combinations, get_ATE_using_nodes_and_cov, get_ATE_using_cov, generate_local_dags_from, to_dict    \n",
    "  from causal_discovery.pc_alg import PCAlgorithm\n",
    "  from causal_discovery.mb_by_mb import MBbyMBAlgorithm\n",
    "  from causal_discovery.sd_alg import SequentialDiscoveryAlgorithm\n",
    "  from causal_discovery.ldecc import LDECCAlgorithm\n",
    "  \n",
    "  try:\n",
    "    from cPickle import dumps, loads, HIGHEST_PROTOCOL as PICKLE_PROTOCOL\n",
    "  except ImportError:\n",
    "    from pickle import dumps, loads, HIGHEST_PROTOCOL as PICKLE_PROTOCOL\n",
    "\n",
    "# Make sure ipyparallel is still able to execute functions.\n",
    "dview = parallel_client[:]\n",
    "ar = dview.map(lambda x: x, (i for i in range(0, 20000, 2)))\n",
    "assert ar.get()[0] == 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "mMjn9PkChGOp"
   },
   "outputs": [],
   "source": [
    "def execute_iteration(params, num_samples, seed, run_oracle=False, \n",
    "                      progress_log_filepath=None):\n",
    "  np = numpy\n",
    "  nx = networkx\n",
    "  \n",
    "  np.random.seed(seed)\n",
    "\n",
    "  # linearly scale down alpha.\n",
    "  alpha = (32000 / num_samples) * 0.01\n",
    "\n",
    "  result = {\n",
    "      \"graph_true_sample_cov_ate\": None,\n",
    "      \"graph_true\": None,\n",
    "      \"oracle\": {\n",
    "          \"ate\": None,\n",
    "          \"neighbors\": None,\n",
    "          \"parent_sets\": None,\n",
    "      },\n",
    "      \"oracle_sample_cov\": {\n",
    "          \"ate\": None,\n",
    "      },\n",
    "      \"pc-alg\": {\n",
    "          \"ci_tests\": None,\n",
    "          \"ate\": None,\n",
    "          \"neighbors\": None,\n",
    "          \"parent_sets\": None,\n",
    "      },\n",
    "      \"ldecc\": {\n",
    "          \"ci_tests\": None,\n",
    "          \"ate\": None,\n",
    "          \"neighbors\": None,\n",
    "          \"parent_sets\": None,\n",
    "      },\n",
    "      \"ldecc-checks\": {\n",
    "          \"ci_tests\": None,\n",
    "          \"ate\": None,\n",
    "          \"neighbors\": None,\n",
    "          \"parent_sets\": None,\n",
    "      },\n",
    "      \"mb-by-mb\": {\n",
    "          \"ci_tests\": None,\n",
    "          \"ate\": None,\n",
    "          \"neighbors\": None,\n",
    "          \"parent_sets\": None,\n",
    "      },\n",
    "      \"sd-alg\": {\n",
    "          \"ci_tests\": None,\n",
    "          \"ate\": None,\n",
    "          \"neighbors\": None,\n",
    "          \"parent_sets\": None,\n",
    "      },\n",
    "      # oracle CI results.\n",
    "      \"pc-alg-oracle\": {\n",
    "          \"ci_tests\": None,\n",
    "      },\n",
    "      \"ldecc-oracle\": {\n",
    "          \"ci_tests\": None,\n",
    "      },\n",
    "      \"sd-alg-oracle\": {\n",
    "          \"ci_tests\": None,\n",
    "      },\n",
    "      \"mb-by-mb-oracle\": {\n",
    "          \"ci_tests\": None,\n",
    "      },\n",
    "  }\n",
    "\n",
    "  A, noise_vars, labels, cov_matrix = params\n",
    "  graph_true = adj_matrix_to_nx_graph(A, labels)\n",
    "  data = get_dataset(mixing_matrix=A, noise_vars=noise_vars,\n",
    "                     num_samples=num_samples, labels=labels)\n",
    "  sample_cov = data.cov().values\n",
    "\n",
    "  result[\"graph_true\"] = graph_true\n",
    "  result[\"graph_true_sample_cov_ate\"] = get_ATE_using_cov(graph_true, sample_cov)\n",
    "\n",
    "  cpdag_oracle = dag_to_cpdag(graph_true)\n",
    "  result[\"oracle\"][\"ate\"] = [get_ATE_using_cov(g, cov_matrix) \n",
    "                             for g in generate_local_dags_from(cpdag_oracle)]\n",
    "  result[\"oracle\"][\"neighbors\"] = get_neighbor_types(cpdag_oracle, \"X\")\n",
    "  result[\"oracle\"][\"parent_sets\"] = [list(g.predecessors(\"X\")) for g in generate_local_dags_from(cpdag_oracle)]                     \n",
    "  \n",
    "  result[\"oracle_sample_cov\"][\"ate\"] = [get_ATE_using_cov(g, sample_cov)\n",
    "                                        for g in generate_local_dags_from(cpdag_oracle)]\n",
    "\n",
    "  pc_alg = PCAlgorithm(use_ci_oracle=False, alpha=alpha)\n",
    "  cpdag_pc, _, _ = pc_alg.run(data)\n",
    "  result[\"pc-alg\"][\"ci_tests\"] = to_dict(pc_alg.ci_test_calls)\n",
    "  result[\"pc-alg\"][\"neighbors\"] = get_neighbor_types(cpdag_pc, \"X\")\n",
    "  result[\"pc-alg\"][\"parent_sets\"] = [list(g.predecessors(\"X\")) for g in generate_local_dags_from(cpdag_pc)]\n",
    "  result[\"pc-alg\"][\"ate\"] = [get_ATE_using_cov(g, sample_cov)\n",
    "                             for g in generate_local_dags_from(cpdag_pc)]\n",
    "  \n",
    "  ldecc_alg = LDECCAlgorithm(use_ci_oracle=False, alpha=alpha, ldecc_do_checks=False)\n",
    "  result_ldecc = ldecc_alg.run(data)\n",
    "  result[\"ldecc\"][\"ci_tests\"] = to_dict(ldecc_alg.ci_test_calls)\n",
    "  result[\"ldecc\"][\"neighbors\"] = {\n",
    "      \"parents\": result_ldecc[\"tmt_parents\"],\n",
    "      \"children\": result_ldecc[\"tmt_children\"],\n",
    "      \"unoriented\": result_ldecc[\"unoriented\"],\n",
    "      \"non_colliders\": result_ldecc[\"non_colliders\"],\n",
    "  }\n",
    "  tmt_par = list(result_ldecc[\"tmt_parents\"])\n",
    "  result[\"ldecc\"][\"parent_sets\"] = list(list(p) + tmt_par for p in get_all_combinations(result_ldecc[\"unoriented\"], result_ldecc[\"non_colliders\"]))\n",
    "  ate_ldecc = [get_ATE_using_nodes_and_cov(list(data.columns), tmt_par + par, sample_cov)\n",
    "               for par in get_all_combinations(result_ldecc[\"unoriented\"], result_ldecc[\"non_colliders\"])]\n",
    "  result[\"ldecc\"][\"ate\"] = ate_ldecc\n",
    "\n",
    "  ldecc_alg = LDECCAlgorithm(use_ci_oracle=False, alpha=alpha, ldecc_do_checks=True)\n",
    "  result_ldecc = ldecc_alg.run(data)\n",
    "  result[\"ldecc-checks\"][\"ci_tests\"] = to_dict(ldecc_alg.ci_test_calls)\n",
    "  result[\"ldecc-checks\"][\"neighbors\"] = {\n",
    "      \"parents\": result_ldecc[\"tmt_parents\"],\n",
    "      \"children\": result_ldecc[\"tmt_children\"],\n",
    "      \"unoriented\": result_ldecc[\"unoriented\"], \n",
    "      \"non_colliders\": result_ldecc[\"non_colliders\"],\n",
    "  }\n",
    "  tmt_par = list(result_ldecc[\"tmt_parents\"])\n",
    "  result[\"ldecc-checks\"][\"parent_sets\"] = list(list(p) + tmt_par for p in get_all_combinations(result_ldecc[\"unoriented\"], result_ldecc[\"non_colliders\"]))\n",
    "  ate_ldecc = [get_ATE_using_nodes_and_cov(list(data.columns), tmt_par + par, sample_cov)\n",
    "               for par in get_all_combinations(result_ldecc[\"unoriented\"], result_ldecc[\"non_colliders\"])]\n",
    "  result[\"ldecc-checks\"][\"ate\"] = ate_ldecc\n",
    "\n",
    "  mb_by_mb_alg = MBbyMBAlgorithm(use_ci_oracle=False, alpha=alpha)\n",
    "  result_mb = mb_by_mb_alg.run(data)\n",
    "  result[\"mb-by-mb\"][\"ci_tests\"] = to_dict(mb_by_mb_alg.ci_test_calls)\n",
    "  result[\"mb-by-mb\"][\"neighbors\"] = {\n",
    "      \"parents\": result_mb[\"tmt_parents\"],\n",
    "      \"children\": result_mb[\"tmt_children\"],\n",
    "      \"unoriented\": result_mb[\"unoriented\"],\n",
    "      \"non_colliders\": result_mb[\"non_colliders\"],\n",
    "  }\n",
    "  tmt_par = list(result_mb[\"tmt_parents\"])\n",
    "  result[\"mb-by-mb\"][\"parent_sets\"] = list(list(p) + tmt_par for p in get_all_combinations(result_mb[\"unoriented\"], result_mb[\"non_colliders\"]))\n",
    "  ate_mb = [get_ATE_using_nodes_and_cov(list(data.columns), tmt_par + par, sample_cov)\n",
    "            for par in get_all_combinations(result_mb[\"unoriented\"], result_mb[\"non_colliders\"])]\n",
    "  result[\"mb-by-mb\"][\"ate\"] = ate_mb\n",
    "\n",
    "  sd_alg = SequentialDiscoveryAlgorithm(use_ci_oracle=False, alpha=alpha)\n",
    "  result_sd = sd_alg.run(data)\n",
    "  result[\"sd-alg\"][\"ci_tests\"] = to_dict(sd_alg.ci_test_calls)\n",
    "  result[\"sd-alg\"][\"neighbors\"] = {\n",
    "      \"parents\": result_sd[\"tmt_parents\"],\n",
    "      \"children\": result_sd[\"tmt_children\"],\n",
    "      \"unoriented\": result_sd[\"unoriented\"],\n",
    "      \"non_colliders\": result_sd[\"non_colliders\"],\n",
    "  }\n",
    "  tmt_par = list(result_sd[\"tmt_parents\"])\n",
    "  result[\"sd-alg\"][\"parent_sets\"] = list(list(p) + tmt_par for p in get_all_combinations(result_sd[\"unoriented\"], result_sd[\"non_colliders\"]))\n",
    "  ate_sd = [get_ATE_using_nodes_and_cov(list(data.columns), tmt_par + par, sample_cov)\n",
    "            for par in get_all_combinations(result_sd[\"unoriented\"], result_sd[\"non_colliders\"])]\n",
    "  result[\"sd-alg\"][\"ate\"] = ate_sd\n",
    "\n",
    "  # Results when a CI oracle is used.\n",
    "  if run_oracle:\n",
    "    pc_alg = PCAlgorithm(use_ci_oracle=True, graph_true=graph_true)\n",
    "    cpdag_pc, _, _ = pc_alg.run(data)\n",
    "    result[\"pc-alg-oracle\"][\"ci_tests\"] = to_dict(pc_alg.ci_test_calls)\n",
    "    \n",
    "    ldecc_alg = LDECCAlgorithm(use_ci_oracle=True, graph_true=graph_true)\n",
    "    result_ldecc = ldecc_alg.run(data)\n",
    "    result[\"ldecc-oracle\"][\"ci_tests\"] = to_dict(ldecc_alg.ci_test_calls)\n",
    "    \n",
    "    sd_alg = SequentialDiscoveryAlgorithm(use_ci_oracle=True, graph_true=graph_true)\n",
    "    result_sd = sd_alg.run(data)\n",
    "    result[\"sd-alg-oracle\"][\"ci_tests\"] = to_dict(sd_alg.ci_test_calls)\n",
    "\n",
    "    mb_by_mb_alg = MBbyMBAlgorithm(use_ci_oracle=True, graph_true=graph_true)\n",
    "    result_mb = mb_by_mb_alg.run(data)\n",
    "    result[\"mb-by-mb-oracle\"][\"ci_tests\"] = to_dict(mb_by_mb_alg.ci_test_calls)\n",
    "\n",
    "  if progress_log_filepath is not None:\n",
    "    os.system(\"echo '%d' >> %s\" % (seed, progress_log_filepath))\n",
    "\n",
    "  np.random.seed(None)\n",
    "  return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "7ii8I2ZPlWYh"
   },
   "outputs": [],
   "source": [
    "def async_results_to_results(async_results, results):\n",
    "  for N, ar in zip(results[\"sample_sizes\"], async_results):\n",
    "    if N not in results:\n",
    "      results[N] = []\n",
    "    for r in ar.get():\n",
    "      results[N].append(r)\n",
    "  return results\n",
    "\n",
    "def execute_estimation_in_parallel(results, start_graph, end_graph, iters_per_graph=1):\n",
    "\n",
    "  num_threads = len(parallel_client.ids)\n",
    "  dview = parallel_client[:]\n",
    "\n",
    "  num_graphs = end_graph - start_graph\n",
    "  print(\"Executing %d iterations across %d engines\" % (\n",
    "      num_graphs*iters_per_graph, num_threads))\n",
    "\n",
    "  dview[\"get_dataset\"] = get_dataset\n",
    "  dview[\"get_noise_matrix\"] = get_noise_matrix\n",
    "  dview[\"execute_iteration\"] = execute_iteration\n",
    "\n",
    "  # The random seed below replicates the results in the paper.\n",
    "  RANDOM_SEED = 716697\n",
    "  np.random.seed(RANDOM_SEED)\n",
    "  new_param_list = []\n",
    "  for _ in range(start_graph, end_graph):\n",
    "    NUM_COVARIATES = 20\n",
    "    NUM_MEDIATORS = 3\n",
    "    A, noise_vars, labels, cov_matrix = get_params(num_conf=NUM_COVARIATES,\n",
    "                                                    num_med=NUM_MEDIATORS,\n",
    "                                                    p_link_conf=0.10,\n",
    "                                                    p_link_conf_med=0.05,\n",
    "                                                    p_link_med=0.1,\n",
    "                                                    p_link_x=0.1,\n",
    "                                                    p_link_y=0.1)\n",
    "    params = (A, noise_vars, labels, cov_matrix)\n",
    "    graph_true = adj_matrix_to_nx_graph(A, labels)\n",
    "\n",
    "    for graph_iter in range(iters_per_graph):\n",
    "      new_param_list.append(params)\n",
    "      results[\"params\"].append(params)\n",
    "      results[\"ground_truth_ate\"].append(get_ATE_using_cov(graph_true, cov_matrix))\n",
    "      \n",
    "  sample_sizes = [8000, 16000, 24000, 32000]\n",
    "  results[\"sample_sizes\"] = sample_sizes\n",
    "  async_results = []\n",
    "  for N in sample_sizes:\n",
    "    print(\"%s: Starting iterations %d for N=%d\" % (\n",
    "        str(datetime.now()), len(new_param_list), N))\n",
    "    \n",
    "    def execute_estimation_iteration(i):\n",
    "      return execute_iteration(new_param_list[i], N, (RANDOM_SEED + N + i),\n",
    "                               run_oracle=(N==sample_sizes[0]))\n",
    "    \n",
    "    async_result = dview.map(execute_estimation_iteration, range(len(new_param_list)))\n",
    "    async_results.append(async_result)\n",
    "\n",
    "  np.random.seed(None)\n",
    "  return results, async_results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "rvMPJwZfodBA"
   },
   "outputs": [],
   "source": [
    "results = {}\n",
    "results.update({\n",
    "    \"params\": [],\n",
    "    \"ground_truth_ate\": [],\n",
    "})\n",
    "\n",
    "results, async_results = execute_estimation_in_parallel(results, \n",
    "                                                        start_graph=0,\n",
    "                                                        end_graph=250,\n",
    "                                                        iters_per_graph=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Vfbe49lusOGd"
   },
   "outputs": [],
   "source": [
    "# Monitor progress of the tasks.\n",
    "for ar in async_results:\n",
    "  ar.wait_interactive(timeout=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Fg8iX-g7sYBv"
   },
   "outputs": [],
   "source": [
    "results = async_results_to_results(async_results, results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "CiIL04UYOh4o"
   },
   "source": [
    "The code below plots the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "lDvbGe-Y1cOP"
   },
   "outputs": [],
   "source": [
    "# For the color map:\n",
    "# https://gist.github.com/AndiH/c957b4d769e628f506bd\n",
    "\n",
    "# Tableau 20 Colors\n",
    "tableau20 = [(31, 119, 180), (174, 199, 232), (255, 127, 14), (255, 187, 120),  \n",
    "             (44, 160, 44), (152, 223, 138), (214, 39, 40), (255, 152, 150),  \n",
    "             (148, 103, 189), (197, 176, 213), (140, 86, 75), (196, 156, 148),  \n",
    "             (227, 119, 194), (247, 182, 210), (127, 127, 127), (199, 199, 199),  \n",
    "             (188, 189, 34), (219, 219, 141), (23, 190, 207), (158, 218, 229)]\n",
    "             \n",
    "# Tableau Color Blind 10\n",
    "tableau20blind = [(0, 107, 164), (255, 128, 14), (171, 171, 171), (89, 89, 89),\n",
    "             (95, 158, 209), (200, 82, 0), (137, 137, 137), (163, 200, 236),\n",
    "             (255, 188, 121), (207, 207, 207)]\n",
    "  \n",
    "# Rescale to values between 0 and 1 \n",
    "for i in range(len(tableau20)):  \n",
    "    r, g, b = tableau20[i]  \n",
    "    tableau20[i] = (r / 255., g / 255., b / 255.)\n",
    "for i in range(len(tableau20blind)):  \n",
    "    r, g, b = tableau20blind[i]  \n",
    "    tableau20blind[i] = (r / 255., g / 255., b / 255.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "U9mLTMRc1kkK"
   },
   "outputs": [],
   "source": [
    "SMALL_SIZE = 8\n",
    "MEDIUM_SIZE = 8\n",
    "BIGGER_SIZE = 8\n",
    "\n",
    "plt.rc('font', size=SMALL_SIZE)          # controls default text sizes\n",
    "plt.rc('axes', titlesize=SMALL_SIZE+6)     # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=MEDIUM_SIZE+6)    # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=SMALL_SIZE + 4)    # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=SMALL_SIZE+4)    # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=SMALL_SIZE+2)    # legend fontsize\n",
    "plt.rc('figure', titlesize=BIGGER_SIZE+20)  # fontsize of the figure title"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "2MAKzoHd1mJj"
   },
   "outputs": [],
   "source": [
    "plot_prop = {\n",
    "    \"pc-alg\": [\"dashed\", tableau20blind[0], \"o\"],\n",
    "    \"ldecc\": [\"dashdot\", tableau20blind[1], \"^\"],\n",
    "    \"ldecc-checks\": [\"solid\", tableau20blind[5], \"s\"],\n",
    "    \"mb-by-mb\": [\"dashdot\", tableau20blind[7], \"D\"],\n",
    "    \"sd-alg\": [\"dotted\", tableau20blind[3], \"v\"],\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "NmlANMMf2HlC"
   },
   "outputs": [],
   "source": [
    "def get_hausdorff_distance(arr_1, arr_2):\n",
    "\n",
    "  arr_1 = np.array(arr_1)\n",
    "  arr_2 = np.array(arr_2)\n",
    "\n",
    "  def get_distance_point_from_array(point, arr):\n",
    "    return np.min(np.abs(arr - point))\n",
    "\n",
    "  dist_1_to_2 = np.max([get_distance_point_from_array(p, arr_2) for p in arr_1])\n",
    "  dist_2_to_1 = np.max([get_distance_point_from_array(p, arr_1) for p in arr_2])\n",
    "\n",
    "  return max(dist_1_to_2, dist_2_to_1)\n",
    "\n",
    "def plot_mse(results):\n",
    "  plt.figure(figsize=(6, 4))\n",
    "  samples = results[\"sample_sizes\"]\n",
    "\n",
    "  dict_to_name = {\n",
    "      'pc-alg': \"PC\", \n",
    "      'ldecc': \"LDECC\", \n",
    "      'ldecc-checks': \"LDECC-checks\",\n",
    "      'sd-alg': \"SD\",\n",
    "      'mb-by-mb': \"MB-by-MB\",\n",
    "  }\n",
    "\n",
    "  errors = collections.defaultdict(list)\n",
    "\n",
    "  for N in samples:\n",
    "    for k in dict_to_name.keys():\n",
    "      errors[k].append([])\n",
    "\n",
    "    for i, res in enumerate(results[N]):\n",
    "      for k in dict_to_name.keys():\n",
    "        errors[k][-1].append(get_hausdorff_distance(\n",
    "            res[\"oracle\"][\"ate\"], res[k][\"ate\"]\n",
    "        )**2)\n",
    "\n",
    "      \n",
    "  samples = np.array(samples)\n",
    "  mses = {}\n",
    "  stds = {}\n",
    "  for k in dict_to_name.keys():\n",
    "    errors[k] = np.array(errors[k])\n",
    "    mses[k] = np.mean(errors[k], axis=-1)\n",
    "    stds[k] = np.std(errors[k], axis=-1)\n",
    "\n",
    "    plt.plot(samples, mses[k], label=dict_to_name[k],\n",
    "             color=plot_prop[k][1], linestyle=plot_prop[k][0],\n",
    "             marker=plot_prop[k][2])\n",
    "    ci = 1.96 * stds[k] / np.sqrt(errors[k].shape[-1])\n",
    "    plt.errorbar(samples, mses[k],\n",
    "                 yerr=ci, ls=\"none\", color=plot_prop[k][1])\n",
    "  \n",
    "  plt.title(\"MSE vs Sample size\")\n",
    "  plt.xlabel(\"Sample size\")\n",
    "  plt.ylabel(\"MSE (Hausdorff distance)\")\n",
    "  plt.legend()\n",
    "  plt.show()\n",
    "  \n",
    "plot_mse(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "7q0qKVkeAZd8"
   },
   "outputs": [],
   "source": [
    "def plot_accuracy(results):\n",
    "\n",
    "  dict_to_name = {\n",
    "      'pc-alg': \"PC\", \n",
    "      'ldecc': \"LDECC\", \n",
    "      'ldecc-checks': \"LDECC-checks\",\n",
    "      'sd-alg': \"SD\",\n",
    "      'mb-by-mb': \"MB-by-MB\",\n",
    "  }\n",
    "  \n",
    "  def is_correct_adj_set(r, key):\n",
    "    return set(frozenset(p) for p in r[\"oracle\"][\"parent_sets\"]) == set(frozenset(p) for p in r[key][\"parent_sets\"])\n",
    "\n",
    "  N = len(results[results[\"sample_sizes\"][0]])\n",
    "  print(\"Total = %d\" % N)\n",
    "\n",
    "  same_count = collections.defaultdict(lambda: 0)\n",
    "  correct_count = collections.defaultdict(lambda: collections.defaultdict(lambda: 0))\n",
    "  \n",
    "  for s in results[\"sample_sizes\"]:\n",
    "    # for i in tqdm(range(N)):\n",
    "    for i in range(N):\n",
    "      for k in dict_to_name.keys():\n",
    "        if is_correct_adj_set(r=results[s][i], key=k):\n",
    "          correct_count[k][s] += 1\n",
    "        \n",
    "  plt.figure(figsize=(6, 4))\n",
    "\n",
    "  accuracy = {}\n",
    "  samples = results[\"sample_sizes\"]\n",
    "  for k in dict_to_name.keys():\n",
    "    accuracy[k] = np.array([correct_count[k][s] * 100 / N for s in samples])\n",
    "\n",
    "    plt.plot(samples, accuracy[k], label=dict_to_name[k],\n",
    "             color=plot_prop[k][1], linestyle=plot_prop[k][0],\n",
    "             marker=plot_prop[k][2])\n",
    "  \n",
    "  print(accuracy)\n",
    "  plt.title(\"Accuracy vs Sample size\")\n",
    "  plt.xlabel(\"Sample size\")\n",
    "  plt.ylabel(\"Accuracy (%)\")\n",
    "  plt.legend()\n",
    "  plt.show()\n",
    "\n",
    "plot_accuracy(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "FOmTtLHaBBX2"
   },
   "outputs": [],
   "source": [
    "def plot_recall(results):\n",
    "\n",
    "  dict_to_name = {\n",
    "      'pc-alg': \"PC\", \n",
    "      'ldecc': \"LDECC\", \n",
    "      'ldecc-checks': \"LDECC-checks\",\n",
    "      'sd-alg': \"SD\",\n",
    "      'mb-by-mb': \"MB-by-MB\",\n",
    "  }\n",
    "\n",
    "  N = len(results[results[\"sample_sizes\"][0]])\n",
    "  \n",
    "  correct_count = collections.defaultdict(lambda: collections.defaultdict(lambda: 0))\n",
    "  \n",
    "  for s in results[\"sample_sizes\"]:\n",
    "    for i in range(N):\n",
    "      for k in dict_to_name.keys():\n",
    "\n",
    "        if frozenset(results[s][i][\"graph_true\"].predecessors(\"X\")) in set(frozenset(p) for p in results[s][i][k][\"parent_sets\"]):\n",
    "          correct_count[k][s] += 1\n",
    "  \n",
    "  plt.figure(figsize=(6, 4))\n",
    "  \n",
    "  coverage = {}\n",
    "  samples = results[\"sample_sizes\"]\n",
    "  for k in dict_to_name.keys():\n",
    "    coverage[k] = np.array([correct_count[k][s] * 100 / N for s in samples])\n",
    "\n",
    "    plt.plot(samples, coverage[k], label=dict_to_name[k],\n",
    "             color=plot_prop[k][1], linestyle=plot_prop[k][0],\n",
    "             marker=plot_prop[k][2])\n",
    "  \n",
    "  plt.title(\"Recall vs Sample size\")\n",
    "  plt.xlabel(\"Sample size\")\n",
    "  plt.ylabel(\"Recall (%)\")\n",
    "  plt.legend()\n",
    "  plt.show()\n",
    "\n",
    "plot_recall(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "aDhNjpg8g_sG"
   },
   "outputs": [],
   "source": [
    "def plot_number_of_tests(results):\n",
    "\n",
    "  plt.figure(figsize=(6, 4))\n",
    "\n",
    "  dict_to_name = {\n",
    "      'pc-alg': \"PC\", \n",
    "      'ldecc': \"LDECC\", \n",
    "      'ldecc-checks': \"LDECC-checks\",\n",
    "      'sd-alg': \"SD\",\n",
    "      'mb-by-mb': \"MB-by-MB\",\n",
    "  }\n",
    "\n",
    "  samples = results[\"sample_sizes\"]\n",
    "  samples = np.array(samples)\n",
    "  \n",
    "  tests_count = collections.defaultdict(lambda: collections.defaultdict(lambda: 0))\n",
    "  for k in dict_to_name.keys():\n",
    "    for s in samples:\n",
    "      tests_count[k][s] = np.mean([t[k][\"ci_tests\"][\"total\"] for t in results[s]])\n",
    "\n",
    "  for k in dict_to_name.keys():\n",
    "    plt.plot(samples, [tests_count[k][s] for s in samples], label=dict_to_name[k],\n",
    "             color=plot_prop[k][1], linestyle=plot_prop[k][0],\n",
    "             marker=plot_prop[k][2])\n",
    "\n",
    "  \n",
    "  plt.title(\"Number of CI tests\")\n",
    "  plt.xlabel(\"Sample size\")\n",
    "  plt.ylabel(\"Number of CI tests\")\n",
    "  plt.legend(loc=(0.67, 0.555))\n",
    "  plt.show()\n",
    "\n",
    "plot_number_of_tests(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "-EQQ352JhMt9"
   },
   "outputs": [],
   "source": [
    "def plot_test_statistics(results, range_start=0, range_end=30000):\n",
    "  plt.figure(figsize=(6, 4))\n",
    "\n",
    "  res = [r for i, r in enumerate(results[8000]) if i%5==0]\n",
    "  total_tests_pc = [r[\"pc-alg-oracle\"][\"ci_tests\"][\"total\"]\n",
    "                        for r in res]\n",
    "  total_tests_ldecc = [r[\"ldecc-oracle\"][\"ci_tests\"][\"total\"]\n",
    "                        for r in res]\n",
    "  total_tests_sd = [r[\"sd-alg-oracle\"][\"ci_tests\"][\"total\"]\n",
    "                        for r in res]\n",
    "  total_tests_mb = [r[\"mb-by-mb-oracle\"][\"ci_tests\"][\"total\"]\n",
    "                        for r in res]\n",
    "  \n",
    "  plt.hist([total_tests_pc, total_tests_ldecc, total_tests_sd, total_tests_mb],\n",
    "           label=[\"PC\", \"LDECC\", \"SD\", \"MB-by-MB\"],\n",
    "           range=[0, 20000],\n",
    "           color=[plot_prop[\"pc-alg\"][1], plot_prop[\"ldecc\"][1], \n",
    "                  plot_prop[\"sd-alg\"][1], plot_prop[\"mb-by-mb\"][1]])\n",
    "  plt.legend()\n",
    "  plt.title(\"Distribution of CI tests\")\n",
    "  plt.xlabel(\"Number of tests\")\n",
    "  plt.ylabel(\"Number of graphs\")\n",
    "  plt.show()\n",
    "\n",
    "plot_test_statistics(results)"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}