notebooks/causal_discovery_with_latent_confounders.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Causal Discovery in the Presence of Latent Confounders\n",
    "\n",
    "We learn a causal graph from observed data, in the possible presence of latent confounders and selection bias.\n",
    "Here we assume that the true underlying causal diagram is a directed acyclic graph (DAG), $\\mathcal{D}(\\mathbf{O}, \\mathbf{L}, \\mathbf{S})$, where $\\mathbf{O}$, $\\mathbf{L}$, and $\\mathbf{S}$ are observed, latent confounders, and selection variables, respectively.\n",
    "For this DAG, we assume causal Markov and faithfulness assumptions.\n",
    "In this setting, a maximal ancestral graph (MAG) defined over $\\mathbf{O}$ is often used to represent the DAG.\n",
    "\n",
    "The ICD algorithm learns a partial ancestral graph (PAG),\n",
    "which is an equivalence class of the true causal graph (an equivalence class of MAGs).\n",
    "\n",
    "For demonstrating the usage of ICD we follow these steps:\n",
    "1. Setup experiment (create training data)\n",
    "2. Define a conditional independence function to be used by the ICD algorithm\n",
    "3. Learn a causal graph using the ICD algorithm\n",
    "4. Calculate the accuracy of the learned graph\n",
    "5. Compare the performance of the ICD algorithm to the FCI algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initially, the required classes and methods are imported."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('..')\n",
    "\n",
    "import random\n",
    "import numpy as np\n",
    "from causal_discovery_algs import LearnStructICD, LearnStructFCI\n",
    "from causal_discovery_utils.cond_indep_tests import CondIndepParCorr  # import a CI test that estimates partial correlation\n",
    "from experiment_utils.synthetic_graphs import create_random_dag_with_latents, sample_data_from_dag\n",
    "from causal_discovery_utils.performance_measures import calc_structural_accuracy_pag, find_true_pag\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Experiment Setup\n",
    "\n",
    "Initially, we will create the data that will be used in this experiment. Here, we will randomly create a small causal DAG\n",
    "and create a probabilistic model by assigning linear relations and normally distributed noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "rand_seed = 42  # arbitrary seed (ASCII code of the Asterisk symbol)\n",
    "random.seed(rand_seed)\n",
    "np.random.seed(rand_seed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Setup simulation parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "num_nodes = 15\n",
    "num_records = 1000\n",
    "connectivity_coeff = 2.0\n",
    "min_lin_coeff = 0.5  # minimal 'strength' of an edge\n",
    "max_lin_coeff = 2.0  # maximal 'strength' of an edge\n",
    "alpha = 0.01"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate a random DAG, $\\mathcal{D}(\\mathbf{O},\\mathbf{L},\\mathbf{S=\\emptyset})$, to serve as the\n",
    "true underlying causal diagram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "dag, observed_set, latents_set = create_random_dag_with_latents(\n",
    "    num_nodes,\n",
    "    connectivity_coeff\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find the PAG (equivalence class) corresponding to the true DAG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "true_pag = find_true_pag(dag, observed_set)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sample a dataset from the DAG by quantifying the edges and adding normally distributed noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "dataset = sample_data_from_dag(\n",
    "    dag,\n",
    "    num_records,\n",
    "    min_edge_weight=min_lin_coeff, max_edge_weight=max_lin_coeff\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define a Conditional Independence Test\n",
    "\n",
    "Constraint-based algorithms, such as ICD and FCI, employ a statistical hypothesis test for estimating\n",
    "if two variables are statistically independent, conditioned on a set of variables.\n",
    "Several such conditional independence (CI) are implemented.\n",
    "Here, we use a CI test based on partial correlation, which is suitable for our linear-Gaussian case.\n",
    "A class instance is created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "par_corr_icd = CondIndepParCorr(\n",
    "    dataset=dataset,\n",
    "    threshold=alpha,\n",
    "    count_tests=True,\n",
    "    use_cache=True\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We enabled *count_test* so the number of unique CI tests is recorded per conditioning set size.\n",
    "It is not necessary to for the normal operation of the CI test and the ICD algorithm\n",
    "Also, We enabled caching to make sure that the statistic for the same CI test is not calculated more than once.\n",
    "Nevertheless, ICD does not perform the same CI test more than once."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learn a PAG using the ICD algorithm\n",
    "\n",
    "First an instance of the ICD class is created. The CI-test to be used during learning is provided here.\n",
    "Then, the method *learn_structure* is called for learning the graph.\n",
    "At this stage, the learned graph is stored within the class instance and\n",
    "the number of CI tests that were performed is stored in the CI test class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "icd = LearnStructICD(observed_set, par_corr_icd)  # instantiate an ICD learner\n",
    "icd.learn_structure()  # learn the PAG\n",
    "learned_pag_icd = icd.graph"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate the Accuracy of the Learned Graph\n",
    "\n",
    "Here, we calculate the precision and recall of the skeleton (ignoring the edge-marks),\n",
    "and the accuracy of the orientation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "icd_result = calc_structural_accuracy_pag(\n",
    "    pag_tested=learned_pag_icd,\n",
    "    pag_correct=true_pag\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare the Performance the ICD Algorithm to the FCI Algorithm\n",
    "\n",
    "First, we learn a PAG using the FCI algorithm. We repeat the same steps used for learning with ICD."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# Learn the PAG\n",
    "par_corr_fci = CondIndepParCorr(dataset=dataset, threshold=alpha, count_tests=True, use_cache=True)  # CI test\n",
    "fci = LearnStructFCI(observed_set, par_corr_fci)  # instantiate an ICD learner\n",
    "fci.learn_structure()  # learn the PAG\n",
    "learned_pag_fci = fci.graph\n",
    "\n",
    "# Calculate structural errors: FCI algorithm\n",
    "fci_result = calc_structural_accuracy_pag(pag_tested=learned_pag_fci, pag_correct=true_pag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For easy plotting, we define the following function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def print_structural_accuracy(structural_accuracy: dict):\n",
    "    print('Edge precision: {:.2f}'.format(structural_accuracy['edge_precision']))\n",
    "    print('Edge recall: {:.2f}'.format(structural_accuracy['edge_recall']))\n",
    "    print('F1 Score: {:.2f}'.format(structural_accuracy['edge_F1']))\n",
    "    print('Orientation accuracy: {:.2f}'.format(structural_accuracy['orientation_correctness']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we print the accuracy of the learned PAG, and plot the number of CI tests per conditioning set size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "ICD performance\n",
      "---------------\n",
      "Edge precision: 1.00\n",
      "Edge recall: 0.73\n",
      "F1 Score: 0.84\n",
      "Orientation accuracy: 0.52\n",
      "Total number of CI tests:  246\n",
      "\n",
      "FCI performance\n",
      "---------------\n",
      "Edge precision: 1.00\n",
      "Edge recall: 0.50\n",
      "F1 Score: 0.67\n",
      "Orientation accuracy: 0.36\n",
      "Total number of CI tests:  662\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA57klEQVR4nO3dd3hUZdr48e9NCCWUAIoBQQgiRXBXIKBiQSLq2lg7qFiw8fr+dNVdXERZdV+Vta7uKq6uBVFBIhZcRLFngBWRIkhX6QSkSkkIJST374/nJBlCEibJzJyZ5P5c17ky58w5z7kJydw5TxVVxRhjjAGo5XcAxhhjYoclBWOMMUUsKRhjjCliScEYY0wRSwrGGGOK1PY7gKo48sgjNTU1tdLX7969mwYNGoQvoCiJ17jBYveLxR59sRz33Llzt6pq89Lei+ukkJqaypw5cyp9fSAQoG/fvuELKEriNW6w2P1isUdfLMctImvKes+qj4wxxhSxpGCMMaZIxJKCiIwWkc0isijo2DsiMt/bVovIfO94qojsCXrvpUjFZYwxpmyRbFMYA4wC3iw8oKoDC1+LyN+BnUHnr1DVbhGMxxhjypWXl0dWVhZ79+6tclnJycksXbo0DFFVXr169WjdujWJiYkhXxOxpKCq00QktbT3RESAAcBZkbq/McZUVFZWFo0aNSI1NRX3MVV52dnZNGrUKEyRVZyqsm3bNrKysmjXrl3I10kkJ8TzksJkVT2hxPE+wDOq2jPovMXAT8Au4C+qOr2MMocAQwBSUlLSMjIyKh1fTk4ODRs2rPT1JR0zfjzZnTuzo3v3omNN5s2j0bJlrLv66rDdJ9xxR5PF7g+LPTTJycm0b9++ygkBID8/n4SEhDBEVXmqyooVK9i5c+dBx9PT0+cWfv6WelGkNiAVWFTK8ReBoUH7dYEjvNdpwDqg8eHKT0tL06rIzMys0vWH+Ppr1SOPdF9L2w+TsMcdRRa7Pyz20CxZsiRsZe3atStsZVVFaf8mYI6W8bka9d5HIlIbuAx4p/CYqu5T1W3e67nACqBjtGOrsvR0mDABBgyABx90XydMcMeNMSYO+NEl9WxgmapmFR4QkeYikuC9PhboAKz0IbaqS0+H//1feOQRaN8eTjvN74iMMXEkISGBbt26FW2rV68GYNasWfTp04dOnTrRvXt3brnlFnJzcxkzZgx33HFH2O4fsYZmERkP9AWOFJEs4CFVfQ24Chhf4vQ+wMMikgcUALep6q+Rii2ipkyBF16A00+Hb7+FzEz43e/8jsoYE25PPgm9eh1cE5CZCbNnw7BhlS62fv36zJ8//6BjmzZt4sorryQjI4PevXsD8N5775GdnV3p+5Qlkr2PSm1ZVdXBpRx7H3g/UrFETWYmXHUV7NrlfmB27IBrr4WxYyEvDy66yO8IjTHh0qvXwVXEmZnF+2H2wgsvcMMNNxQlBIArrrgi7PcBG9EcXrNnw0svwT/+ASeeCOef735ARo2C/v2hRPY3xsS4vn1hzBj3Oi/P7Y8d6/ZPPhmOPhouvti1IV55pdvfvh0A2bbNnf/RR+78jRtDuuWePXuKqo4uvfRSABYtWkRaWlrY/lnliesJ8WJOaY+M6emuKumzz6BbN3ds926I0dkTjTEV0LQpHHeca0P8859h1qwqF1la9VE02ZNCOOXnu0fInJyDjycmFlcd/fQTtGtX/NeDMSZ2BQIweLB7nZjo9q+91u0nJcFDD8G0afDAA/D6627/sssA0COOcOf37+/Ob9Gi0mF07dqVuXPnVvr6irCkEE7LlsFZZ8EHH5R9TnKyOydogJsxJg4FtyE8/HBxd/TMzLDf6o477uCNN97gu+++Kzr2wQcfsGnTprDfy5JCOKWmut5H5fU2SkmBjAxo3RpUYfjwsDxyGmOibPbsg8chFY5Tmj077LdKSUkhIyODe+65h06dOnH88cfz2WefRWQaDWtTCKcGDeC880I/f+tWlyCSkuCkkyIXlzEm/MpqQ6ziYNWcktXPnt69ezN9+qGz/wwePJjBhVVcYWBJIZwyMqBnT9fwFIrmzV2PpMJsv2SJO9a81FXyjDEm4qz6KFyys+Gaa4q7q4WqSRNISICCAhg40DVKRXCSQmOMKY89KYRLw4awfDnUq1e562vVcgllzx4QcYkhPx9q23+RMSZ67EkhXETg2GPd4JXKOvFEOOUU9/rZZ10vpV27whOfMcaEwJJCuIweDZ9+Gr7yWrRwvZl8XKTDGFPzWFIIl0cfhbffDl9511wDb77pnkA2b4Z+/dyo6GCZmW6OJWOMCRNLCuGybBk880xkyp40CaZPh6uvLh4YUzhwplevyNzTGOOL4FXmfvrpJy644AI6dOhAjx49GDBgAJs2bSIQCJCcnEz37t3p1KkTffr0YfLkyWG5v7VihkudOnDkkZEp+5Zb3IC45cthwAA6nHIKzJxpC/gY47Nx42DECFi7Ftq0gZEjYdCg8JS9d+9eLrzwQp555hn6e1NlBAIBtmzZAsAZZ5xRlAjmz5/PJZdcQv369enXr1+V7mtPCuHw2mvw9NORvccxx7gE0L8/rSZPhlNPtYRgjI/GjYMhQ2DNGtdZcM0atz9uXHjKf/vtt+ndu3dRQgDo27cvJ5xwwiHnduvWjQcffJBRo0ZV+b6WFMLh669dFU+kZWbCRx/xa8+e8M03bt/GNBgTMeXNnH3ffZCbe/D5ublw993u9bZtUpmZs4tUdLrsHj16sGzZsordpBSWFMJh3Dj46quil6mpbthBamr4/moInnxrwVNPwbvvuv0ePdzsjMaYqMrKKv34tm3RjaOQhukPRGtTCJfExKLHycK/HgofJyEM9YzBk28FAu7r6NFwzz1u7iRjTNgFAsWvC2fOLtSmjfsdL6lNG/f1iCP0oPMrOnN2165dmTp1asjnz5s3j+OPP75iNymFPSlU1fjxcN11sGcPI0aU/jg5YkQY7jNs2KFtCP37u15PAwe6/Y8/dg3QxpiIGzny0L/HkpLc8XC45pprmDFjBh9//HHRsWnTprFo0aJDzl2wYAGPPPIIt99+e5XvG7GkICKjRWSziCwKOvZXEVkvIvO97YKg9+4TkeUi8qOIxM9K95s2weLFUK8ea9eWfkpZx8NCxH0tKHALfdxzj7UzGBMFgwbByy9D27bu17BtW7cfrt5H9evXZ/LkyTz//PN06NCBLl268K9//Yvm3oSZ06dPL+qSevvtt/Pcc89VuecRRLb6aAwwCnizxPFnVfWgrjoi0gW4CugKHA18KSIdVTU/gvGFx913F7UsHe5xMqJq1XIN3jk57ic0N9e9PuqoKNzcmJpp0KDwJYFCwVNnd+7cmU9LmSkhJSWFnTt3hvfGnog9KajqNODXEE+/GMhQ1X2qugpYDsTdAgORfpw8rCZN3OI94NaL7d4dIvSDY4ypnvxoU7hDRBZ41UtNvWOtgHVB52R5x2Lb55/DySfDihVA8eNkSop7u1Wr8D5OVshtt8G997rlP40xJkQSrm5MpRYukgpMVtUTvP0UYCugwCNAS1W9SURGATNVdax33mvAFFV9r5QyhwBDAFJSUtIyMjIqHV9OTs5BQ8orqumcObQdO5aFjz1Gfv365OYmcPfd3bjxxlX07h3qQ1LFVTTuBqtW0f6FF/jx3nvZ5/MCPlX9nvvJYvdHNGNPTk6mffv2SGFbXRXk5+eTkJAQhqgqT1VZsWLFIVVN6enpc1W1Z5kXRWoDUoFFh3sPuA+4L+i9z4Dehys/LS1NqyIzM7NK15e0YoVqv36qgYDbLygIa/FFKhz3Rx+pduigunFjROKpiHB/z6PJYvdHNGNfuXKlbtmyRQvC8Mu7a9euMERUeQUFBbplyxZduXLlIe8Bc7SMz9WojlMQkZaq+ou3eylQ2DNpEvC2iDyDa2juAMT2avaFT1hBf1Eceyx8+aV7ffbZ0LIlvPWWD7GVdNFFbu3o2rVd3K+/7mZhreyCQMZUU61btyYrK6tofqGq2Lt3L/V8/h2rV68erQvbGUMUsaQgIuOBvsCRIpIFPAT0FZFuuOqj1cD/AKjqYhGZACwBDgC3a6z3PPrxRzf/0Ntvuw9c3DD4xET39llnQdOm5VwfbYUruH37Ldx8s0sON9/sb0zGxJjExETatWsXlrICgQDdu3cPS1nRFLGkoKpXl3L4tXLOHwlEq59O1SUkwBVXQPv2AKxbB126uCUQLr0U7r/f5/jKcuqp8N//Qu/ebn/r1sjN7mqMiTs2ormyOnRwXYs6dADc2LHrr4fgUeb5+W6LOaed5sY17NgBaWlhGnJtjKkOLClU1vbtB+22bQsvvACdO7v96dPdGIUZM3yILVQNG8KNN8Ill7gV3AoX8ClkK7sZU+NYUqiM3Fxo3rzoA7Og4NCRzMcd5wY6x/SA4tq14a9/dau39erl5lK6807X3mAruxlTI9ksqZVx4AA8/ribXB2YNw969oT334fLLnOntGwJTzzhX4gVduaZcOKJrkosORleeslWdjOmBrKkUBmNG7uJ5zytW8M//wlnnHHwaQUFrtq+WbPohlcptWq5Oq+//AUefdRNrmcJwZgax6qPKmPZsoPmyE5JcbUuJQcLDxgAffpEObaqmDoVXnnFTZHxwguHtjEYY6o9SwqVcfbZcMstgMsNn38Oe/YcetoNN8DQoVGOrbIK2xD+/W+39e/v9i0xGFOjWPVRRanCiy/CEUcA7o/rCy6Azz6Dc889+NSg9bZjX/DKbhMmwOmnu6w2e7ZVIxlTg1hSqCiRgz7tzzwTpkw5tD0BXP7YuBHq13ezWse0YcOKX19xhfvaooUlBGNqGKs+qqjZs2Hp0qLdpCQ3y0X9+oeeumEDHH20W7Ez7kyeDB984HcUxpgos6RQUUOHugFfuJU4n3sONm8u/dSjj4ZRo4p6rsaXp5+Gv//d7yiMMVFm1UcV9eqrrp8prj3hrrtcD6PSBqmJQBjW0fbH2LHFqwUZY2oMSwoV1bFj0cvCAb9t25Z9+vbtsHKlm2IorlRwul1jTPVg1UcVMXMmTJx40Cx37dq5cV9lefppt2JnXl4U4gu3995z6y5EcHU+Y0xssaRQEf/+N/zP/0CtWixZ4qqG1q4t/5JBg2DSpOiEF3Zbtrh1I7zqMmNM9WdJoSJefBGmTQMRliyBN94oXrumLF26uHEMhYvvxJXbboO5c2NstSBjTCRZUqiIevWK5sa+4grXXnD00eVfoup6sS5bFoX4wq1wqdH8fKtCMqaGsKQQqvnz3bSnQesohPrX/7nnuq6rcWnqVDfl6+LFfkdijIkCSwqhmjrVrbFZqxaZmW76o9WrD3+ZiGubDh4wHFc6d3YLThtjagRLCqG66y7Ytg2Sk8nJcQ8MoS6g07cvpKZGMrgISkmBjAw44QS/IzHGREHEkoKIjBaRzSKyKOjYUyKyTEQWiMhEEWniHU8VkT0iMt/bXopUXFXiTWDUv79rf01KCu2yDRvcHHN790YutIjbssVtxphqLZJPCmOA80oc+wI4QVV/C/wE3Bf03gpV7eZtt0UwropbvRpuvhmWLatUm+u0aTBwICxfHpHoIm/nTjeY7fnn/Y7EGBNhEUsKqjoN+LXEsc9V9YC3OxOIj2Gzq1bBf/4D+/czfjy0aQPr1oV++bnnwg8/HDQYOr4kJ7tFd666yu9IjDERJhrBroYikgpMVtVDKqRF5CPgHVUd6523GPf0sAv4i6pOL6PMIcAQgJSUlLSMjIxKx5eTk0PDhg1DO9n7Ps3/oSlTprTg3nuXlTuSOZIqFHeMsdj9YbFHXyzHnZ6ePldVe5b6pqqWuwHtgbre677AnUCTw13nnZ8KLCrl+AhgIsVJqS5whPc6DVgHND5c+WlpaVoVmZmZVbq+IqZMUf3kk/CUFc24DzJrlvuHVIFvsYeBxe6PeI09luMG5mgZn6uh/K37PpAvIscBLwPHAG9XPDc5IjIYuAgY5AWHqu5T1W3e67nACiA2Klv27nXdhz75hJ07XfV6ZYwcCY89FtbIom/48DjuW2uMCUUoSaFAXTvApcDzqvpnoGVlbiYi5wHDgN+ram7Q8eYikuC9PhboAKyszD3CbvNm2L8fCgp47TW3CmdZ6yeUZ+xYt25NXCuc5sMYU22FMnV2nohcDdwAFK5DedixvCIyHlfddKSIZAEP4Xob1QW+EDeFwkx1PY36AA+LSB5QANymqr+WWnC0tWkDM2YA0O8HN6g51PEJwcqbXjtuxG1LuTEmVKEkhRuB24CRqrpKRNoBbx3uIlW9upTDr5Vx7vu4aqqYduKJbquMX36Bt9+Gyy+P44FsAF995abU/te/iudGMsZUG6FUH52jqneq6ngAVV0FxPMwrIo5+WR48klWrYJ586CgoHLFbNsG99zjJseLa8uXw0cfVa4OzRgT80JJCjeUcmxwmOOITXl5bu7rFi34979dfsjNPfxlpencGbZuhSuvDG+IUXfjjW4RCVuq05hqqczqI68d4RqgnYgELxPTmBKD0qqtxER4/XUA7t7oOiFVtttx7dqukTru1anjdwTGmAgqr01hBvALcCTw96Dj2cCCSAYVM/LyiubHbtECzis5aUcFffgh/Pwz/PnPVQ/NV3PnwvXXu0aSyjayGGNiUpnVR6q6RlUDwNnAdFWdiksSrYGa0cJ4/vlwxRV8952bKLSq6yx/9plrn417bdq46qM9e/yOxBgTZqH0PpoGnCEiTYHPgdnAQGBQJAOLCf37Q/36jBkD77xT9faAf/yjmtS+NG8OX3/tdxTGmAgIJSmIquaKyM3Av1T1SRGZH+G4YsNddwEwKt/1HEpIqFpxdeuGIaZYsnev27wpxY0x8S+U3kciIr1xTwYfe8eq+PEYB7Ztg337AJcM2revepG7dsHdd7uu/nEvN9ctUP30035HYowJo1CSwt24kcgTVXWxNw1FZkSjigX33gvt2vH+e8r//R8cOHD4Sw6nfn0YM6aaLHeclAQjRlS99d0YE1MOW33kNTBPFZEkb38lbqbU6u2qq+Dkk/lmhjBlCjz0UNWLTEx0y3hWm4HAQ4f6HYExJswO+6QgIr1FZAmwzNs/UUSqQx+a8p19Ntx6K888A/Pnh6/YapMQCq1ZA9NLXfrCGBOHQqk++gfwO6BwausfcBPYVV8bN7oBBd7COuFsIJ40yfViiuDaRtF1yy1uqdJq8w8ypmYLae0wVS25+GR+BGKJHa+/Dh078tTDe7j++vB+3m3eDEuXVn5dhpjz1FPw6afV8BHImJoplC6p60TkVEBFJBG4C1ga2bB8NnAgtGvHnp+S2L07vJ93t9zitmqjWze/IzDGhFEoSeE24J9AK2A9bgDb/4tkUL479lg49lge9DuOePH992467ZEj7YnBmDgXSvVRJ1UdpKopqnqUql4LHB/pwHyzdStMnkzBzuyI3eLaa93o5mpj7lx47jk3e6oxJq6FkhSeD/FY9fDVV9C/P7fduI9+/SJzi507YffuyJTti0GDYMuWarK8nDE1W3lTZ/cGTgWai8ifgt5qTHUe0dy/P0ydSo+FTWkZoXVkPvooMuX6JinJ7wiMMWFSXptCHaChd06joOO7gCsiGZSvkpKgTx9uq96dbsNv2TK4/Xb4+9+t8dmYOFZmUggayTxGVddEMSb/7N8PL7zA5tMupVmPVGqH0gxfCZmZboK9Dz6oRjUuzZvDpk2uGskYE7cO26ZQlYQgIqNFZLOILAo61kxEvhCRn72vTb3jIiLPichyEVkgIj0qe99KW7QI/vQnbr29Dr16Re42DRvCUUdVs+UIjjjCff/OOcfvSIwxVRDS4LUqGAOUnDFtOPCVqnYAvvL2Ac4HOnjbEODFCMd2qB494JdfuGVYM+65J3K36dULpkxx6zZXO6rVLNsZU7NENCmo6jQOXc/5YuAN7/UbwCVBx99UZybQRERaRjK+UrVoQf8r6zGo+i8hFH7790OHDvDww35HYoypJNEy5nAQkeeBMid4UNWQZkoVkVRgsqqe4O3vUNUm3msBtqtqExGZDDyuqv/13vsKuFdV55QobwjuSYKUlJS0jIyMUMIoVU5ODg0bNizab/fqq3x71O9I6H0czZvvq3S5oXj88c7k5QkPPFDxweEl444lqa+/TnbHjmw77bRS34/l2A/HYvdHvMYey3Gnp6fPVdWepb1XXlPqnHLeCwtVVRGp0MxCqvoy8DJAz549tW/fvpW+fyAQoOj6nBy45BKuSn4Cvm3Ft99WutiQfPONW/O5b9+UCl97UNyx5jBxxXTsh2Gx+yNeY4/XuMtLCu8AjVT1oO4kItIcqMpw300i0lJVf/GqhwpHA6wHjgk6r7V3LDoaNoRff+Xl7/PI3h/5240YEfl7+GbXLli50rqmGhOHymtTeA44o5TjpwPPVuGek4AbvNc3AP8JOn691wvpFGCnqv5ShftUXK1a/KZnXU49NXq3rJYzTl99NVx+eTX9xxlTvZWXFNJU9YOSB1V1IiGupyAi44FvgU4ikiUiNwOPA+eIyM/A2d4+wCfASmA58ArRnnTv3nuZOOSTqK0Xs3ix65Y6eXJ07hdVDzwAb73ldxTGmEoor/qovLkLQl2H4eoy3jpkViF1Ld63h1JuRHz/PcPn/InOm+CM0p6PwqxVK7j4YmjRIvL3irpTTvE7AmNMJZWXFDaLyEmqOiv4oIj0AqrfsNUvvmDOLuXX7dG5XZMm8Mor0bmXL1asgIkT3TrONp22MXGjvKTwZ2CCiIwB5nrHegLXA1dFOC5fNGosNGoc3Xvu2xfe5T5jxtSpMGyYm2CwUye/ozHGhKjMaiDvCeEkQIDB3ibAyar6XTSCi5oHHuCJnu8yfnx0b3vnnW6sV7U0YABs2GAJwZg4U+6Ub6q6GXgoSrH4RhUyVp1Er0zXcSZazjkHjj7a3b/a1bA0bOg2Y0xcifTcR3FBHn2E77e25bnnonvf/v1h+PBqmBAKrVnjFuCZN8/vSIwxIbKk4PWlF4F69aJ/+717ITtyK3/6q1Ej+PprN5DNGBMXykwKIvKW9/Wu6IXjg2ee4Yamk/jn01EYxlxCbi40aACjRkX91tHRrBmsX+8Gshlj4kK5g9dE5GjgJhFp6q2DULRFK8BIK2jdhl8btiF7b52o3zspCZ54AtLTo37r6Knl/Yjl5/sbhzEmJOU1NL+EW+/gWFyX1OCab/WOx71aA6/ko4H+3T+S6zbEhPx86NMHTj/dZUBjTEwrr0vqc6p6PDBaVY9V1XZBW7VICJKXh+Yd8DWG/fvhp598DSGyEhLg5JOta6oxcSKU5Tj/V0ROFJE7vO230QgsGo6YMYPudZfw19v9G6D9z3+6z8sdO3wLIfKeeQZuusnvKIwxIThsUhCRO4FxwFHeNk5E/hDpwKIh56hWnHnidjqe1MS3GC68EN58ExITfQshOg4ccFNfGGNiWrmD1zy34EYx7wYQkSdwM58+H8nAomHv8cfxz3ln+hpDly5uq/YGDYJZs2D0aL8jMcaUI5RxCgIEdx3J5+BG5/h04AC5P25HC/yf83/lymr+R/STT0Lv3q4aqaDAHcvMdMeNMTEllKTwOvCdiPxVRP4KzARei2hUUZA/fyHX3nYm91zsfytvv37w4IN+RxFBvXrByJFuatiEBJcQBgxwx40xMeWw1Ueq+oyIBHArrgHcqKpxP29BXvOj+UPfrzj7Zv8HCbz4YjVdV6FQejpMmABXXskJHTvCzz+7/Wo9QMOY+BRKmwKq+j3wfYRjiY4nn2TczgsZMa4ra9dewRurhJGzFzMo+WM31bMPzjvPl9tGV3o69OnDkRMnwo03WkIwJkbVuLmPxu28kCF/a8eaNaAqrFkDQ/7WjnE7L/Qtph074MsvISfHtxAiLzMTpk9n/UUXwUcfuX1jTMypcUlhxLiu5JZYaTSXJEaM6+pTRDBjhptGe/5830KIrMI2hAkT+HnoUFd1NGCAJQZjYlC5SUFEEkSkWv3mrl1bsePR0Lu3m0z0xBP9iyGiZs8uakOQ/fth8mS49VZ33BgTU8pNCqqaDxSISHK4biginURkftC2S0Tu9no3rQ86fkG47hmsTZuKHY+Gpk1dFXujRv7FEFHDhhW1IWhiIkyb5uYq96kNxxhTtlAamnOAhSLyBbC78KCq3lmZG6rqj0A3cE8iwHpgInAj8KyqPl2ZckM1ctBihvyt3UFVSEnkMnLQKsC/KqTZs2HnTjj7bN9CiA4R+PZbqB1SHwdjTJSF8pv5gbdFQj9ghaqukSgtPzYo+WO4v7D3kdKmjTBy0Cp33Mek8PDDrgrrhx98CyF6ChPC3r3+rGxkjCmTqB5+RK+I1AfaeH/lh+/mIqOB71V1lDcwbjCwC5gDDFXV7aVcMwQYApCSkpKWkZFR6fvn5OTQMEbWEV63rj6JiQW0aLHvsOfGUtwVVRh76wkTOOadd5g5fjxaJ/prWVRGdfi+x6N4jT2W405PT5+rqj1LfVNVy92A/sCPwCpvvxsw6XDXhVBuHWArkOLtpwAJuHaOkbgpu8stIy0tTasiMzOzStf7JV7jVg2KPTNTdehQ1e3bfYymYqrF9z0OxWvssRw3MEfL+FwNpUvqX4GTgB1eEplPeBbYOR/3lLDJK3eTquaragHwinfPGmPHDnj9dVi+3O9IoqRvX3j6aTf1hTEmZoSSFPJUdWeJYwVhuPfVwPjCHRFpGfTepcCiMNwjbmRnuyUHvvzS70iibP582LTJ7yiMMZ5QksJiEbkGSBCRDiLyPDCjKjcVkQbAORzcgP2kiCwUkQVAOvDHqtwj3rRq5aYEuvlmvyOJonXroHt3eOUVvyMxxnhC6X30B2AEsA/3l/1nwCNVuam6tRmOKHHsuqqUGe9q1YLjjvM7iig75hh47z2bB8mYGBLKLKm5wAhvcR1V1ezIh1UzTZsG8+bBXXf5HUkUXX653xEYY4KEshxnLxFZCCzADWL7QUTSIh9azTNlCgwfDvn5hz+3WpkyBf7xD7+jMMYQWpvCa8D/U9VUVU0FbsctvGPC7L77XC+khAS/I4mySZNg1Ci3jrMxxlehJIV8VZ1euKOq/wXstzcCGjeGunX9jsIHjz0Gy5bZ1BfGxIAyfwtFpIf3cqqI/BvXyKzAQCAQ+dBqnrw8ePxxOPlkOPdcv6OJosKxCqpubiRjjG/Ke1L4u7edCHQEHsINZDseb0I7E161a7u17adO9TsSH8ybB7/9LSxc6HckxtRoZT4pqKr1E4wyETeOK06mAgqvY45x9We7dvkdiTE12mErcUWkCXA9kBp8vlZy6mxTvhqZEACOPBK++cbvKIyp8UJpaP4ElxAWAnODNhMB//2vm+5i716/I/HJvn2QleV3FMbUWKF096inqn+KeCQGcGsqfPopbNwIqal+R+ODU06Bli3hk0/8jsSYGimUpPCWiNwKTMZNdQGAqv4asahqsKuvhmuu8TsKH91/v82caoyPQkkK+4GncPMfFa7Io4Rn+mxTQo3vkXnllX5HYEyNFkqbwlDgOG9Ecztvs4QQQffe67qm1libNsELL0BBOGZoN8ZURChJYTmQG+lATLGlS2HlSr+j8NEXX8Add8D33/sdiTE1TijVR7uB+SKSycFtCtYlNUImTfI7Ap9dfjn06gWdOvkdiTE1TihJ4UNvMyY66te3hGCMT0JZT+GNaARiis2fD8OGwbPPQteufkfjk7174U9/gt694boavf6SMVEVyojmVRT3Oipijc2Rk5joptDOrsnLGdWtC3PnujELxpioCaX6qGfQ63rAlUCzyIRjwD0dzJrldxQ+E4Fvv3XrlBpjouawv3Gqui1oW6+q/wAurOqNRWS1iCwUkfkiMsc71kxEvhCRn72vTat6HxPHChPC7t3+xmFMDRLKcpw9graeInIboT1hhCJdVbupauHTyHDgK1XtAHzl7ddI990HAwf6HUUMeOwxaN8e9u/3OxJjaoRQPtz/HvT6ALAaGBCRaOBioK/3+g3cYj73RuheMa1xY5vtAYBTT4XcXDdRXo2dQtaY6BHVQ9qQo3Nj14C9HdeI/W9VfVlEdqhqE+99AbYX7gddNwQYApCSkpKWkZFR6RhycnJo2LBhpa/3S7zGDRa7Xyz26IvluNPT0+cG1dAcTFXL3YC6wDXA/cCDhdvhrguh3Fbe16OAH4A+wI4S52wvr4y0tDStiszMzCpd75d4jVu1krEXFKjOmKG6cWPY46mIGvd9jxHxGnssxw3M0TI+V0Pp2vEfXLXOAdzo5sKtSlR1vfd1MzAROAnYJCItAbyvm6t6n3i1fj0cfzxU4UGo+lizxlUjjR7tdyTGVHuhtCm0VtXzwnlTEWkA1FLVbO/1ucDDwCTgBuBx7+t/wnnfeNK8OXTpAs2s869bWGLSJOjb1+9IjKn2QkkKM0TkN6oazhXVU4CJrtmA2sDbqvqpiMwGJojIzcAaItegHfPq1IH33/c7ihjSv7/fERhTI4SSFE4HBnsNw/sAAVRVf1vZm6rqSuDEUo5vA/pVttzqSNXWWCgycSKsWAH33ON3JMZUW6G0KZwPdMBV8fQHLvK+mgj729+gRQuXGAzw2Wfw+uuQn+93JMZUW6FMiLcmGoGYQ514opsLbv9+NxVQjffUU9CggU19YUwEhWtksomACy90m/E0auS+Wp2aMRFjf3LFOFU3mNd4Zs2Cjh1h0SK/IzGmWrKkEMMKClzX1Acf9DuSGNKuHbRpA3v2+B2JMdWSVR/FsFq14I9/hB49/I4khjRvDl995XcUxlRblhRi3IgRfkcQo3bvhq1boW1bvyMxplqx6qMYpwpbtriqJONRhV694Pbb/Y7EmGrHkkKMe/11OOooyMryO5IYIgKPPgrDa+xyG8ZEjFUfxbjTToNnn4WkJL8jiTGXXeZ3BMZUS5YUYlynTm4zpcjKggkTXGu8jVswJiys+igO7NgB69b5HUUM+vprNw/SggV+R2JMtWFJIQ707Qu33eZ3FDHoyith1So3H4gxJiys+igOPPwwxOiqfv6qX9+6pBoTZvakEAd+/3s46yy/o4hRublw443w5pt+R2JMtWBJIQ7s3Qvz5sHOnX5HEoPq14dly2DTJr8jMaZasOqjOPD9965r6uTJNmvqIURgxgzrfWRMmNiTQhz4zW/gvffcIF5TiqeegsxM2LWr+FhmJjz5pH8xGROn7EkhDjRqBJdf7ncUMaxXL7joIrew9aZN8M03MGCAG8NgjKkQe1KIEz//DNOn+x1FjEpPh5Ej3RJ1Dz5YnBDS0/2OzJi4E/WkICLHiEimiCwRkcUicpd3/K8isl5E5nvbBdGOLZY9/DBce63fUcSwu++GoUPhiSfcoA5LCMZUih9PCgeAoaraBTgFuF1EunjvPauq3bztEx9ii1nDh8PEiX5HEcMyM+HFF2HIEHjsMXjjDb8jMiYuRb1NQVV/AX7xXmeLyFKgVbTjiDddu/odQQzLzCyuMmrXDgIBNx9Smzb2xGBMBYmq+ndzkVRgGnAC8CdgMLALmIN7mtheyjVDgCEAKSkpaRkZGZW+f05ODg3jZKjw/v21mD27KW3a5NK06Za4ibukSHzPjxk/nuzOndnRvbs7oEqT+fNptGwZ6666KmzdVePp56Ukiz36Yjnu9PT0uaras9Q3VdWXDWgIzAUu8/ZTgARcldZIYPThykhLS9OqyMzMrNL10bRzpyqoPv54fMVdUlRj/9vfVG+4QbWgICzF2ffdH/EaeyzHDczRMj5XfemSKiKJwPvAOFX9AEBVNwW9/wow2Y/YYlXjxjB7NnTs6AazmRAcOAB5ee5rYqLf0RgTF6KeFEREgNeApar6TNDxluraGwAuBRZFO7ZY17P0hz1Tlr/8xX0VcUt42qhnYw7Lj95HpwHXAWeV6H76pIgsFJEFQDrwRx9ii2kLFsBLL/kdRRwRcdvmzdCvH3z3nd8RGRPz/Oh99F+gtD/ZrAvqYUyZ4rqmTp6c4Hco8aVWLdiyxSbNMyYENqI5jtx6q/ujNykp3+9Q4suRR8L8+W4OcmNMuSwpxJFmzaB5c6sar5QE7+lq0iS3YtuBA/7GY0yMsqQQZ156CWbNaup3GPFr0yZYvfrgGVWNMUUsKcSZxx6Dr79O8TuM+HXrrW79hWbN/I7EmJhkSSGOjBsHBQXw+ecppKa6fVMJiYmwbx/cfDN8Yv0bjAlmSSFOjBvn5nrLygJVYc0at2+JoZIOHHCNz4tsOIwxwWyRnTgxYoRboz5Ybq47PmiQPzHFtQYNXDVS3bp+R2JMTLEnhTixdm3FjpsQFCaExYvdLKt79vgbjzExwJJCnGjTpvTjzZrBc8+5tgZTScuXu6eGNWv8jsQY31lSiBMjR0JS0sHHkpKgQwf48EMbu1AlF18MP/0EnTv7HYkxvrOkECcGDYKXX4a2bUFEadvW7c+YUZwUduxwyWPvXr+jjUNJSW7SvGeegbfe8jsaY3xjSSGODBrkxl19/fVUVq92+yJuWm1wyeGhh2DpUh+DjGf5+fDxx26SKR8XnzLGT9b7qBoZPBhOPx2OO87tf/gh9O0LTZr4F1NcqV0b/vMf99Rg9XGmhrInhWqmMCFs3gxXXeWeHEwFNGzoZlXdts0Nbtt+yIqwxlRrlhSqqaOOgm++gQcfdPtZWbB+vb8xxZUff4R334U5c/yOxJiosqRQjaWlwRFHuNd/+AOcfLKb3cGE4NRTXRfVefMgM/Pg9zIz4ckn/YnLmAizpFBDPP00vPBC8XitDRv8jScuNG0KvXrBpZe6rAouIQwY4I4bUw1ZQ3MN0b692wA+/dR1zf/sM9cQbcqRnu4Wx37xRdpt3+6+aePHu+PGVEP2pFAD9egBd94JvXu7fRvXcBgffQRDh9J23Dj31HDppfDll+69HTvcpHq2aI+pJiwp1EBHHQVPPeWqkvbvh1NOgYcfdjOupqa6zjc2NXeQmTNh9GhWX3cdvPcenHOOG0oO7snhN7+BhQvd/sKF8OqrkJPjX7zGVEHMJQUROU9EfhSR5SIy3O94qrv8fFeFtHOnm4p7zRo3bivcU3MXJpyzzjozIgknYgktM5Nxv88gtdZajh37BqmJ6xn3RXNYudK936cPjB0LXbu6/cmT3UI+hZNRTZjg9vfvd/vBg+KefJJxIxYfHPeIxeFrxA4qv+j7Hq7yLfbDlh32uEuUH5HvO4CqxswGJAArgGOBOsAPQJeyzk9LS9OqyMzMrNL1folE3G3bqrpPrIO3pk1V+/UrPu/dd1Xvvbd4f/p01XfeKd5ftUp12bLi/fx81bFjVZOSDi43KckdD4dIlj924CRNqnvg4LLrHtCxAyeVfkF+vurKlcX7Tzyh2rVr8f5dd6medpor+/5FmsTug8tmt469f1HVA49w+RZ79MsOZ/nAHC3rc7isN/zYgN7AZ0H79wH3lXW+JYXwESk9KYDqRRcVnzdsmGr79sX7N92k2qpV8f5116mmphbvX3WVamJi6eW2bat6663unEKDB7syCg0a5O5RaMAA1dtuK96/9FLVRo3KLv/cc1WHDy8+Pz1d9YEHivdPO031kUeK9086SfXxx4v3y4s9L0+1SxfVF1905+7e7fZfe83t79jh9t980+1v2aLapeWvmtHfZatWrUovu3m9naqq+vPPql0ardFP+49SVdUlS1S7NFyjX13+gqqqzp+v2qXBap0+6F+qqjp7tmqXBqt05o0vqapqSkrp5bdstEtVVb/4QrVLveW6dNhoVVX9+GPVLnV/1hUPjlFV1Q8/OKBd6v6sax99Q1VV3x23T7vU/Vl/eXpsmX9E1Ga/7nptgqqqvvLUdu1S92fdO/4DVVX916PbtEvdnzV/4n9UVfWfD2zR39RdpvrJJ6qq+tSfN2lavYWqX35ZZvl1Za/qt9+qquqIwev0rKRvVefMUVXVYVev1d81mK76ww+qqnr3ZWu0f8OvVZcuVVXV2/uv1ssafa5tj95fatkN6ubpdckfqm7Y4H72TlulNzV5z/3HqeqAk1bpbU3Hq+50/z+Xdl+lf2j2lmpurqqqXnjCKm1UK7vMn5f049boAy3+Xfyz13atPtL6peKfvVZZ+njqv4r2u6Vs0L8fV7zf9chftGntnWWWXxHlJYVY633UClgXtJ8FnBx8gogMAYYApKSkEAgEKn2znJycKl3vl0jEfdRRp7BpU71Djqek7GXo0JkU3u78891WuH/JJbU555wEAgE3AOK00xpywgl1CAR+BaBDh+bk5XUBDp02Yu1apaBgFbVqJRAIrPKOtgWEQGA1ALVqpVJQoAQCblrrhIRjOXDgAIGAW0giMbE92dmtyyy/W7cN7N6dSyDgRu41aNCR7OwcAgHXJ7dx407s3LmLQOAXAJo06cz27TsIBDYCkJd3ZpllT506lebNu7Bp02YCga3s31+L5s07s2HDJgKBbezZk0Dz5p3IytpIIPAr2dm1ad6xI2vPOIFAIMCGDaWXvXVvIwKBAJs31yU1WdhVy/1///JLPVKTlR1kEwgEWLeuPu2S89mW7/ZXr06iXeM8NuftJBAIsGlT6eVvzG5IIBDgp58acWzjHLKyt7AxEGDFisa0b7yLNTs3szYQYNWqRrRvlM3KHZtZEQiwdl1D2jfKYfm2jaxdq6WWfYDaLN6Qxd5AgM3bkmjfKJcFa1azOxBg2/b6tG+0h3krV5IdCLB9Zz3aNc7l++XL2RUIkJ1Tl9TGOcz9aTlr155Vavn7tA6zlyxh99695OXV4ZjkLcxatIPc7GzyC+rQKnkz3y1YzZ5ffwXq0DJ5KzPnbWTvxo3Ukjq0SN7G2vWlf+zt3pdA85QdfDNrFnnJySTVrUPT5J38d+ZMDjRsSIN6dWnWZBfTvvmGgvr1aZRUl6ZNspk6fTpapw5NGtQnu6BtqWWvXaucdOJGGiTsKfrdTWncgPr5uUX7LZMbUjfo/dZNGlEnqfj9Nk2SWby1W5nlBwJTS32vwsrKFn5swBXAq0H71wGjyjrfnhTCJ5JVMGX91VfRv278KD9ey450+RZ79MsOZ/mU86QQaw3N64FjgvZbe8dMhB08NTdFU3OHY6nPstaCGDmy6mVHuvyIlj1oMUkcvMZqErmMHLS46oVHuHyLPfplR6N8iL3eR7OBDiLSTkTqAFcBk3yOqcYonJq7oICiqbnDVW5pa0FEpvzwJrRIxj4o+WNevn/VwXHfv4pByR9XvfBDytewlm+xRz/uQ8sP//cdiK3qI/dUwwXAT7heSCPKO9eqj+KPxe4Piz36Yjlu4qihGVX9BPjE7ziMMaYmirXqI2OMMT6ypGCMMaaIJQVjjDFFLCkYY4wpIq4hOj6JyBZgTRWKOBLYGqZwoile4waL3S8We/TFctxtVbV5aW/EdVKoKhGZo6o9/Y6jouI1brDY/WKxR1+8xm3VR8YYY4pYUjDGGFOkpieFl/0OoJLiNW6w2P1isUdfXMZdo9sUjDHGHKymPykYY4wJYknBGGNMkRqZFETkPBH5UUSWi8hwv+MJlYgcIyKZIrJERBaLyF1+x1RRIpIgIvNEZLLfsVSEiDQRkfdEZJmILBWR3n7HFAoR+aP3s7JIRMaLyKHL68UIERktIptFZFHQsWYi8oWI/Ox9bepnjGUpI/anvJ+XBSIyUUSa+BhiyGpcUhCRBOAF4HygC3C1iHTxN6qQHQCGqmoX4BTg9jiKvdBdwFK/g6iEfwKfqmpn4ETi4N8gIq2AO4GeqnoCkIBboyRWjQHOK3FsOPCVqnYAvvL2Y9EYDo39C+AEVf0tbjmA+6IdVGXUuKQAnAQsV9WVqrofyAAu9jmmkKjqL6r6vfc6G/fB1MrfqEInIq2BC4FX/Y6lIkQkGegDvAagqvtVdYevQYWuNlBfRGoDScAGn+Mpk6pOA34tcfhi4A3v9RvAJdGMKVSlxa6qn6vqAW93Jm4lyZhXE5NCK2Bd0H4WcfTBWkhEUoHuwHc+h1IR/wCGAQU+x1FR7YAtwOte1derItLA76AOR1XXA08Da4FfgJ2q+rm/UVVYiqr+4r3eCKT4GUwV3ARM8TuIUNTEpBD3RKQh8D5wt6ru8jueUIjIRcBmVZ3rdyyVUBvoAbyoqt2B3cRuNUYRr/79YlxSOxpoICLX+htV5XkrhsVdH3oRGYGr+h3ndyyhqIlJYT1wTNB+a+9YXBCRRFxCGKeqH/gdTwWcBvxeRFbjquzOEpGx/oYUsiwgS1ULn8rewyWJWHc2sEpVt6hqHvABcKrPMVXUJhFpCeB93exzPBUiIoOBi4BBGieDwmpiUpgNdBCRdiJSB9fwNsnnmEIiIoKr116qqs/4HU9FqOp9qtpaVVNx3/OvVTUu/mpV1Y3AOhHp5B3qByzxMaRQrQVOEZEk72enH3HQQF7CJOAG7/UNwH98jKVCROQ8XHXp71U11+94QlXjkoLX8HMH8BnuF2SCqi72N6qQnQZch/sre763XeB3UDXEH4BxIrIA6Ab8zd9wDs97snkP+B5YiPt9j9mpF0RkPPAt0ElEskTkZuBx4BwR+Rn35PO4nzGWpYzYRwGNgC+839WXfA0yRDbNhTHGmCI17knBGGNM2SwpGGOMKWJJwRhjTBFLCsYYY4pYUjDGGFPEkoKp9kRkjIhc4b1+tXASQRG5v8R5M6pwj6Jyo0lEBovI0RU4//fxNDOwiT7rkmqqPREZA0xW1fdKHM9R1Yb+RBUeIhIA7lHVOX7HYqoHe1IwMUVErvfmn/9BRN7yjqWKyNfe8a9EpI13fIyIPCciM0RkZdDTgIjIKG/NjC+Bo4LKD4hITxF5HDd76HwRGee9lxN0/VPeGgQLRWSgd7yvd33hugrjvJHCReUWliMiI71/w0wRSfGOt/f2F4rIo4X3K/HvbyAiH3vXLgq6d5qITBWRuSLymYi09P69PXGD6uaLSP0SZd0pbu2NBSKS4R0bLCKjvNfzg7Y9InKmd//RIjJL3OR/cTGDsAkjVbXNtpjYgK64eeeP9PabeV8/Am7wXt8EfOi9HgO8i/vjpgtuSnSAy3Bz2SfgJoLbAVzhvRfArS8AkFPi/jne18uDrk/BTRfREugL7MTNl1ULN4L19FLKVaC/9/pJ4C/e68nA1d7r20reP+jerwTtJwOJwAyguXdsIDC65H1LKWsDUNd73cT7OhgYVeK8/sB07z5/A64tvMb7/2jg98+GbdHb7EnBxJKzgHdVdSuAqhbOT98beNt7/RZwetA1H6pqgaouoXha5T7AeFXNV9UNwNcVjOP0oOs3AVOBXt57s1Q1S1ULgPlAainX78clAIC5Qef0xiUxgv49JS3ETevwhIicoao7gU7ACXjTJQB/IbS5+RfgniKuxc3SeQgR6QA8BQxQN2neucBw7z4BoB7QJoR7mWqitt8BGFNF+4JeS5Tvl0/pv0N5qqqHOadUqvqTiPQALgAeFZGvgInAYlWt6BKgF+ISZH9ghIj8JvhNcVOwTwBu1eI1CwS4XFV/rOC9TDVhTwomlnwNXCkiR4Bbn9c7PoPiZSQH4ao6yjMNGChuPeiWQHoZ5+WJm4q8pOlB1zfHfbDOqsC/oywzcdVDUMaymF5PolxVHYv7C74H8CPQXLx1oUUkUUS6epdk4yZdK1lOLeAYVc0E7sVVQ5VsVB8NvK6qwd/Pz4A/BLWVdK/wv9LENXtSMDFDVReLyEhgqojkA/NwdeB/wK169mfcCmg3HqaoibiqqCW49oBvyzjvZWCBiHyvqoNKXN8b+AHXPjBMVTeKSOfK/cuK3A2MFbfoyqe49omSfgM8JSIFQB7wv6q632tUfk7c0qC1cavYLca1q7wkInuA3qq6xysnwbtXMu6v/+dUdYf3WY+ItAWuADqKyE3eNbcAj3hlL/ASyyrcegCmhrAuqcZEiYgkAXtUVUXkKlyjs/XuMTHFnhSMiZ40YJRXNbMD15PKmJhiTwrGGGOKWEOzMcaYIpYUjDHGFLGkYIwxpoglBWOMMUUsKRhjjCny/wGeGDGS4s7Z9wAAAABJRU5ErkJggg==\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print('\\nICD performance')\n",
    "print('---------------')\n",
    "print_structural_accuracy(icd_result)\n",
    "print('Total number of CI tests: ', sum(par_corr_icd.test_counter))\n",
    "\n",
    "print('\\nFCI performance')\n",
    "print('---------------')\n",
    "print_structural_accuracy(fci_result)\n",
    "print('Total number of CI tests: ', sum(par_corr_fci.test_counter))\n",
    "\n",
    "plt.figure()\n",
    "x_range = range(len(par_corr_icd.test_counter))\n",
    "plt.plot(x_range, par_corr_fci.test_counter, 'xr:',\n",
    "         x_range, par_corr_icd.test_counter, 'ob:')\n",
    "plt.legend(('FCI', 'ICD'))\n",
    "plt.xlabel('conditioning set size')\n",
    "plt.ylabel('number of CI tests')\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}