notebooks/causal_discovery_with_a_perfect_oracle.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Causal Discovery using a Perfect Oracle\n",
    "\n",
    "A perfect oracle is a conditional independence (CI) test that always yields the true answer.\n",
    "For that, the oracle requires access to the true underlying graph from which it can read-off the true conditional independence relation.\n",
    "Although this is impractical, access to the true graph is granted in a simulated environment for evaluating the number of CI tests required by an algorithm in the large sample limit.\n",
    "\n",
    "For demonstrating how to use a perfect oracle we follow these steps:\n",
    "\n",
    "1. Create a true underlying DAG\n",
    "2. Learn a PAG using ICD and a perfect oracle\n",
    "3. Plot the results\n",
    "\n",
    "Let's start by importing the required classes and methods.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('..')\n",
    "\n",
    "from causal_discovery_algs import LearnStructICD  # import ICD algorithm\n",
    "from causal_discovery_utils.cond_indep_tests import DSep  # import the perfect oracle\n",
    "from graphical_models import DAG  # import a DAG class\n",
    "from plot_utils import draw_graph  # for plotting the graphs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Next, let's create the true underlying DAG, a graph with 7 nodes, and define which nodes are observed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS4AAAEuCAYAAAAwQP9DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxZ0lEQVR4nO3deZgU1dn38e893bP29CwosohsbiiIsisuEcQVEYwoYDY1migxaiKuARM1i48aNRoVvXxwCcaoUXhVNA8iKm4oDiBLQBAVFFkHGGbf+n7/qEZGnKVherqquu/PddUF9FRX3z30/ObUqVPniKpijDF+kuZ2AcYYs7csuIwxvmPBZYzxHQsuY4zvWHAZY3zHgssY4zsWXMYY37HgMsb4jgWXMcZ3LLiMMb5jwWWM8R0LLmOM71hwGWN8x4LLGOM7FlzGGN+x4DLG+I4FlzHGd4JuF2CM8RYR2Qh0aKPDb1LVjq09iNjUzcaYhkSkTUNBVaW1x7BTRWOM71hwGWN8x4LLGOM7FlzGGN+x4DLG7LXrr78eVf3OVldXRygUSsjr21VFY8x3xHJV8YUXXqCsrIzrr7/+28cikQhbtmxp8fjxuKpo47iMMXtt8ODB/O53v2PTpk2uvL61uIwx39FSi6tjx45s2LCBr776iuzsbD7//HP+9Kc/8dJLL8V0fBvHZYxJuIKCAi699FLGjBnDyJEjWbduHS+88AL9+/dPWA3W4jLGfMfejpxPT09n27Zt3Hnnndx+++0t7m8tLmOM62pra6mvr6e+vj5hr2nBZYxplTFjxhAOh3n55ZcT9pp2VdEYE7Nx48YRDAYpKipCRDjjjDO49dZb+fOf/8zSpUsTVocFlzEGEckATgPGN7dfhw4dmDhxIl27dqWsrIyVK1cyceJEpk+fnphCo6xz3pgUJSIB4CRgAvBDYCXwDPD3tnzdeHTOW3AZk0JERIBBOGF1AbAJJ6yeVdV10X08Px+XnSoakwJEpA9OWI0H6nDCariqfupqYfvIgsuYJCUiPXGCagJQAPwLOB9YpD4/1bJTRWOSiIh0wjkFnAD0BJ7HaV29r6qRGI/h+VNFCy5jfE5E2gHn4YRVP+D/4YTVG6patw/Hs+AyxsSfiOQC5+CE1UnAbJywelVVq1p5bM8Hl/VxGeMTIpIJnIETVmcC7+GE1Y9UdWccX2oTbbg8WTwOYi0uYzxMRILAMJywGgMswQmrF1R1awLreAJYrqp3Jeo1m2PBZYzHRMdaHYcTVucDX+GE1XOq+rUL9ewHfAYcmsiwbI6dKhrjAdGwOhpn+MJ4oAInrE5U1dVu1gZcDLzkldACCy5jXCUih+K0rCYA2ThhNRpY4oWxViKSBlwBXOh2LQ1ZcBmTYCLSBRiHE1ZdgOeAS4D5XgirPZwG7AA+crmO77DgMiYBRGR/YCxOWB0FzABuBN7al7FWCTQReMhrgWqd88a0ERHJw7kSOAEYCryGcyr4H1WtdrG0mIhId6AIOEhVK1wu5zusxWVMHIlINnAWTlidCswD/gGcr6plbta2D34BPOW10AJrcRnTaiKSDpyCE1bnAAtxWlYvquo2N2vbV9HBruuAk7w4g4S1uIxvRCe+OzkYDJ6Yl5d3EtAx+vimHTt2zKuvr38P5/68Nl+1IXq17XicsBoLfI4TVjeq6oa2fv0EOA/nyqbnQgusxWV8QEQygsHgVVlZWdd17tw5e/To0dmDBg0KHnjggQB88803LFiwoO6VV16pWLt2bU1NTc1fa2tr7421H0lELgAGqur1LewnQH+csBoHbMeZKuZfqvp5a96j14jIu8BfVXWG27U0SlVts82zG3BMOBxePWzYsPKPPvpIW/Lxxx/r6aefXp6bm/sFMCCG45+HM9izCgg3sU8v4FZgFbAG+CPQ2+3vTRt+z4/GGa0fdLuWJmt0uwDbbGtqA07PyckpnzZtWiQSiWisIpGITp8+PRIKhcqBUc0cf1Q0tBQoBX7c4GtdgeuBRcA3wL3AYKJnKcm8AVOBKW7X0dxmp4rGk0RkaCgUen327Nk5Q4cO3adjLFiwgOHDh1eWlZWdqapv73H8M4EXcEar7zIfeBrnVPBw4EWcfqt5moB+My8QkXzgS+BI9XBfnQWX8RwRCYVCoc+efvrpjqNHj27VsV577TXOP//8reXl5QdrdOoXERkBvMR3Q2uXfwNPArNVtaZVL+5DInIlzv2R49yupTm2krXxnFAo9Jezzz47v7WhBXDmmWdy3nnn5YbD4b8CiMgPcGYIbSy0KoA5qvpKioaWEB0p73YtLXL7XNU22xpuQF5WVlbF+vXrtTm/+c1v9PDDD9dQKKThcFiHDRumixcvbnTfTZs2aVZWViXwM6Aep0+rqW2h298DF7/3JwPL8UE/nrW4jKeIyIWnnHJKpHPnzs3u1717d6ZNm8ayZcuYO3cu5eXlnH/++Y3ue8ABBzBy5MgIzoj2FThXBstwlukqBXYCu1pYx4jIgXF6O37jyfsSG2MDUI2nFBYWjh43blyopf2uuuqqb//evXt3xowZw7333tvk/uPGjct588032xcXF/fZ9ZiIhIEeDbYjcFbGyd33d+BPItIZ5xalS92uJRYWXMZTamtr+w8cODDm/SORCPPnz2fq1KnccsstTe43YMAAqqurj2n4mKqW4kyFvGQfy00ml+IMpI3n3PVtxq4qGk9JS0urr6qqSsvIyGh2vyVLljB06FCqq6vJzMzkiSeeYOzYsU3ur6oEAgFV1XRNkaENsYrea/kFcJaq+iLELbiMZ0SvakXq6+tJS2u++7WiooJvvvmGbdu2cd999zF79mxWrFhB+/btm3xOMBiM1NfX56gPppRJJBH5IfAbVT3R7VpiZcFlPCUrK6vi66+/zt5///1jfs66devo1q0bs2bN4qyzzmp0n9LSUtq1a1dXV1eX4YfO50QSkTnA/6rqM27XEiu7qmg8JRQKrVq0aNFePefdd98F4NBDD21yn0WLFhEOh1dbaH2XiPQC+uDcJeAbFlzGUyorK9+eN29ek31Qt9xyCzNmzGD16tUsW7aMe+65h8svv5wrr7yy2eB65513IjU1Ne+2SdH+djlOa8tXp892qmg8RUT67b///u9u2LAhJxj8/kXvyy+/nP/85z9s2LCBcDjMYYcdxuWXX85Pf/rTJo8ZiUQ48MADyzdu3HiKqn7YlvX7iYiEcCYL7K+qa92uZ2/YcAjjKaq6qKCgYM1zzz131IUXfn9FrKlTp+71MWfMmEFFRcV6PLZSjQdMAN71W2iBtbiMB4nIcfn5+W+sXr06u7mrhLHYtm0bhxxySOX27du/N0NEKotewV2IM2Pr/7ldz96yPi7jOar6QV1d3aNjx46tqKnZ93uda2trGT9+fEVtbe1TFlrfMwQIA6+7Xci+sOAynlReXj5p0aJF80aOHFlRVrb3i+OUl5czZsyYyvnz539YVlb26zYo0e8mAg+rasTtQvaFBZfxJFWtKy0tPefDDz984dBDD62YO3duzM995513OOywwyrefffdl0tLS89Q1do2LNV3oovTjgKecLmUfWZ9XMbzRGRkTk7OE0ceeWTmtddeGz711FPZb7/9vrPPtm3beOONN/jrX/9aunTp0tqKiopLVPX/uVSyp4nI9cARqnqx27XsKwsu4wvR++lGFxYW/qa8vLx/fn5+fceOHesANm7cGCwpKQmGQqFF27dvvw+YoSk4EWAsoku8rQbGqeoCt+vZVxZcxleiy9qvBH7M7jm0NgOf+bW/JpFE5CzgVlUd5HYtrWHjuIxvRC/hPwd0Anqq6mMul+RH/piauQXWOW/85BrgB9G/n+diHb4kIj2AY4Fn3a6ltexU0fiCiBwHvMHuRS4qgTxVrXOvKn8RkTuAdFW91u1aWsuCy3he9PL9CqDhXDc7cSa+e8+dqvxFRLJw7ks8XlVXu11Pa9mpovG06FWwGUDeHl/KAc5OfEW+NRZYlAyhBRZcxvv+APQD9pzLOQicm/Bq/CspOuV3sVNF41kiciSwDGe9w8Z+yVYDB6nqloQW5jMi0g9nEdyeydInaC0u42WrgDHAn4G3cFaajgAlQC0gOEtqmeZdATySLKEF1uIyPiIiDwFfAwuA/sCJwGOqOtPNurxMRApwVvA5QlU3ulxO3FhwGd8QkfnA9ao6z+1a/EJErgKOU9UJbtcSTxZcxheiVxd3Ap38smip26J3GqwAfpFsYW99XMYvDgc2WGjtlWE4fYHvuF1IvFlwGb/oB+zdumVmIvBQMi7JZsFl/KIfzhzpJgYiciAwHJjudi1twYLL+EV/rMW1Ny4DnlHVUrcLaQvWOW88L9rJXIxzSX+T2/V4XXTSxS+B01V1mcvltAlrcRk/6AZUWWjFbDTOxIpJGVpgwWX8wTrm905S3ZfYGAsu4wfWMR8jETkCOAJnRo2kZcFl/MA65mN3Bc5tUEm9WIjNOW/8wE4VYyAiuTiLiBztdi1tzVpcxtNE5ACcSQO/dLkUP7gQeFtVv3K7kLZmwWW8rh/OzJ02bqcZ0SEjvyLJO+V3seAyXmenibE5DmchkTfcLiQRLLiM11nHfGwmAg+nyqK4NnLeeJqIrAbGqOpyt2vxqmg/4Kc4UzNvd7ueRLAWl/EsEckDOuP8UJqmXQK8mCqhBTYcwnjb0cDSZJorPd6iEyxeToqt7G0tLuNl1jHfsjOBTapa5HYhiWTBZbzMOuZblvT3JTbGgst4mbW4miEiBwODgOfcriXR7Kqi8SQRyQR2AIWqWuVyOZ4kInfi/Axf53YtiWad88ar+gCrLbQaJyLZwMXAsW7X4gY7VTReZaeJzTsf+FhV17hdiBssuIxXWcd881KyU34XCy7jVdbiaoKIDAA6Aa+6XYtbLLiM50QHVR4FLHa5FK+6ApiqqvVuF+IW65w3XnQYsFFVS9wuxGtEpBBnlPzhbtfiJmtxGS+y08Sm/Qx4VVU3u12Im6zFZbzIOuYbISJpOJ3yF7tdi9usxWW8yFpcjRsOVALvu12I2yy4jKdEpyC25cgaNxF4yKaxtlt+jMeISDfgA1Xt7HYtXiIiXYAlQFdVLXO7HrdZi8t4jZ0mNu4XwNMWWg7rnDdeYx3zexCRDOAy4BS3a/EKa3EZr7EW1/eNAVaq6n/dLsQrLLiM16R0x7yIZIrIahG5T0R6Rh9O6fsSG2Od88YzRKQ9sApol6pXzqK3O9UBtUA9sBzoAXRW1Wo3a/MSa3EZL+kHLE7V0AKI3n9YAaQDWcAAIBf4SkR+F12KLOVZcBkvsY55R+ke/84A2gNTgPUiclTiS/IWCy7jJSndv9XAjiYerwemASnfSW/BZbzErig6tjXyWCVwi6r+MpWns9nFxnEZTxCRMHAgtmo1wJY9/l0B/FhVZ7hRjBdZcBmvOBpYZqtWA7Ap+qcCJcAZqvqhi/V4jgWX8QrrmN9tZ/TPb4CTVPVzN4vxIgsu4xX9gA/cLqKtRWe/GJaenn5KXl7eyXV1dd1UNRAIBEojkciCkpKSeUAmsB7oq6qN9XelPAsu4xX9SOLR4SKSFggELg2Hw7e0b98+f8KECdmDBg0KHHzwwaSnp1NcXNxx8eLFh86dO3f0rFmzAhkZGa/s3LkzTOMd9SnPRs4b10VXrd6OM2I+6RaAFZHueXl5z/Xo0ePIBx54IHTCCSfgNLwaV1xczAMPPFB31113VdfU1FxbV1f3aCoPym2MBZdxnYj0B55U1aQbWCkiR+Xk5MybMmVK+LrrrgsEAoGYn7t8+XLGjBlTvmHDhifKy8t/beG1m43jMl6QlB3zInJQTk7O24899lj+jTfeuFehBdC7d28WLFgQ6tat20U5OTm3t1GZvmTBZbwg6UbMi4jk5eU9c8MNN4QnTJjQ9HlhCwoKCpg7d24oGAz+VkQGx7NGP7PgMl6QjCPmLzzwwAOPufnmm1t9AaxDhw48+uij2bm5uc9GV/pJefZNMK5KxlWrRUTy8/P/cM8994SCwaZza9q0aZx22mm0b98eEWHhwqYbnRdccAGdOnXaDxjRBiX7jgWXcduhwOYkW7V6QHZ2dqfTTjut2Z127tzJ+PHjueSSS8jIyKBPnz5N7isiTJo0KbegoGBSvIv1Iwsu47b+JFn/FnDyueeeG0xLa/7H65prruGSSy6hvLyco446ioyMjGb3HzVqlFRVVQ2V5sZSpAgLLuO2pOvfKiws/MGQIUMyY92/qKiIAQMGtLhfp06dyMrKEqBba+pLBhZcxm1JF1xpaWndevToEdO+9fX1fPLJJzEFF0DXrl1rgC6tKC8p2C0/LYguDZWPc6f+Dpu9IH4arFrt6+ASkXRgf+AAoH0oFCqsq4vtY7JixQoqKytjDq60tDQB9m5AWBKy4NpD9IdpcDgcviIQCJwUDAYPyszMrAOoqakJtGvXbk1NTc3s8vLyh1TV5o5qnYOAGlXd6HYhDUWvdO5HNIiifzb8+55/hoFiYDOwpaamJvOrr76K6bWKiopIT0/nqKNiu2mguLg4Def2qJRmwdWAiAzOy8t7Micn56Crr746e/jw4Wl9+/YlKysrCM5VoMWLF/eaNWvWwY888shlBQUFH5eUlFysqmvcrt2nEtIxHx37VEjjodPYY4U44bAFJ4w2N/j7kkYe266qkV2vl5mZecdnn312HTF0xRQVFdGnT58WO+YBysvL2bBhQxawMqY3nsTsXkWcD3YoFLojEAhc+fDDD2ePHz+elq4IVVdXc//999f/4Q9/qK6trf1tTU3NIwkqN2mIyK1AQFUn7+XzBOf0vbnwafj3/XAWoGgsiBp7bFtrugRE5JyBAwdOX7BgQbixr9fW1rJ8+XIALrvsMjp16sRtt91GKBTi0EMPbfK4s2fPZty4cSu2b99+5L7WlixSPrhEJC0cDj/Tq1evs19++eWcDh067NXzP/30U0aMGFFRXFx8Z0VFxa1tVGZSEpGXgCeBF3GW4IqlNXQATn9SFc2HT8PHtqpqbQLfV2Z2dvaWRYsWhQ8//PDvfX3+/Pkcd9xx33t87NixPP/8800e9+yzzy5/9dVXr4tEIg/HtWAfSvngys3Nvat3795XvPnmm6GcnJx9OsbGjRsZOHBgxcaNG6+oq6t7Ks4l+pKI5NB0n9Cuv4/ACZh2QITYWkNbgC1en/4mOzv7f84///yrnnrqqax4HG/lypX069evoqqqqqOq7rl8WcpJ6eCK9mm9tXr16uwDDmjdOpuLFy9m6NChZZWVlb1UdX0Lr1sIhFT161a9aAJF58xqKYga/hmk+fCpAh7Fud1ni6qWJ/DttDkRyc/JyVkzc+bM/U499dRWHau+vp7+/fuX//e//72htrb2wTiV6G+qmpIbIHl5eSv++c9/aqwef/xxzc7O1nHjxjX69ZtvvrkmLy/vheZeE7gYp7/lNZfffzrQCeiL0/K5ELga+BNOoMwE3gdW4yzYUIMznfBiYDbwNHAvcBNwKXAOcBxwMJBH9JdiM69/GvCmm9+DBHyPTyssLCxfs2ZNo5+XWEQiEb3qqquq8vLyPgDS3H5PXtlS+ari4FAodNC4ceNi2vnpp58mMzOTyspKjj322Eb3mTRpUvq99957logcoKqbG34tuvrwk8BhQAhnKa64aXAJP9Z+ol2X8Bs7LStq5LESjf40xonvx2+1RFVnZ2Zm3nDsscf+z9y5c3OauxexMbW1tfz2t7+tfuKJJ74qKysbqQ2uXKa6lA2ucDg88eqrr85u6eohwMsvv8zGjRu/HSTYVHAVFhZy7rnn6jPPPHMhcB+AiOTitGIuA7JwWl3ghEyTGlzCj/X0rBCnZdTYadnSRh7b5vIPQj9glouvnxDV1dV/z8jIKBkyZMjDkydPzpw0aVIwPT29xectWbKE8ePHl69fv35hWVnZaFVN+bFb3+F2k8+traCg4IuPPvpIW/L222/rz3/+c1VVveeeezQjI0Orqqqa3P/JJ5/UwsLCV3AC6jycVk0lzsj7hlsV8Avgd8DfgGeAOTjjhDYAtThjiT4F3sW58jYVuB34NTAeGA70AToAQbe/p3uzRd9XH7frSOD77Z6fn/9OYWFh+ZQpU2oXLlyoNTU13/nsbNiwQZ9//nk96aSTSrOyskoDgcBltHDKnapbSnbOi0hmMBgsKy0tDWZlNX3RZ9GiRUyaNIlZs2aRlZXFT3/6U1atWsX8+fObfM7y5csZOnTohp07d1YCPdjdwmrMNJzFPxu7krZVVWv2/t15X3TV6o1AvqbYLVQiclQ4HL46LS3ttIqKio7t27evCgaDWlJSkl5VVaXZ2dnLduzY8SDwvKpWul2vV6XqqWJ+VlZW3a4R8Y1ZtWoVF1100behBU6QnXLKKc0euGPHjlRXV+cDTwE/whlzlM737y+rBCar6oZWvA+/StlVq1V1Kc7FDEQk95tvvjkI57NRApRUV1fvACottJqXqrNDtNjMvOuuu1i6dCndu3cnGAwSDAZZtmwZDzzwAAMHDmz6wKqIiKrqjap6EHACMB0nqBpe8q/F6ZtKRUnfMR8LVS1T1RWqukxVvwIG4Hw2/yEiTX/ITMoG147q6upAWVlZkzv87ne/Y8mSJSxevJjFixczffp0AGbMmMGLL77Y5PPWrVtHZmZm8a5/q2qRql6E04l+DU7fTgWQHX0sFVlwNW4QTnBlAq+LSE+X6/GslAwuVa3Nzc39fPHixU3u0717d/r06fPtVltbSyAQ4IwzzqBr165NPq+oqAjgo0Zes0xVH1PVXsCJwCPA2ta+F5+y4GrcMHZ3KeQB80Rkfxfr8ayUDC6A6urq12fNmhVzH8t///tfDj744Bbv4p85c2Z5SUnJ7Ob2UdWFqvprVV0d6+sni+gI/MNxhmiY72o4KVcaTlfC3OjtU6aBlLyqCCAiR+Tn5xdt3rw5O5YpRWKxfv16DjnkkMqqqqoOaveTNSqZV61uDRFpD3wN7PlhrATeAc5S1fqEF+ZRKdviUtUVwMK///3vcfswTJkypSo9Pf1JC61m2Wli4wbgjO3bUzbwA+CBxJbjbSkbXAAlJSUXTZkypXr16tafsb3++us8++yz5aWlpTfGobRkZsHVuCHAnqeEZUA1u8f3maiUDi5V/ay+vv66ESNGVGzatGmfj7N8+XLGjh1bWVFRMV6Ta33AttAfC67GDIv+WQ7siP75F6CnqnZV1d+7VZgXpXRwAVRVVT20devWvw4aNKhiyZIle/38OXPmcPzxx1eUl5dfqqpz2qDEpBG9EbwvSbRqdRy9AFyH0/JqR3Tsn6p+42pVXuX2PUde2QKBwEXZ2dllkydPrt2+fbu2ZP369XrJJZdUhkKhrcAIt+v3wwb0Ata4XYcfNmAs8KrbdXh1S9mrio0RkYNyc3Mfrq+vP+W8887j9NNPzxowYAAdOnRAVVm3bh1FRUXMnDmz/I033kgLBoP/KCsru17t9DAmIjIBOE9Vx7pdi9eJSDvgS2B/TdJ7VlvDgisqugDDn3AmxushIj8sLCw8ta6u7pjq6up8EdHMzMxiEVmwY8eO/wP+pao73a3aX0TkTpx5vf7kdi1+ICILgGtVdZ7btXiNBRffzpn1DHA6zsjl9qq6zd2qko+IzAH+qqqvuV2LH4jIX4BaVb3F7Vq8JuU750XkMJxR3KfizOJQw/cHAZpWSpZVqxNsDs7n0uwhpYNLREbjTFPcDefGVoB6LLjagidXrfa494A+IpLvdiFek5LBJSIBEbkD+CfOen4NJ/tTLLjagrW29pI6S7B9AJzscimek5LBhTPX+TV8f6QyWHC1FQuufTMHZxUm00CqBtdLOHNiNTYhVwQLrrbQH1jodhE+ZP1cjUjJ4FLVh3DWFPwNzoIUe95onfm9J5nWshbXvlkM7C8iB7ldiJekZHABqGo1zmIVG3HGb32Kc39YNtbiiqvoZHhh4Au3a/EbdZaQewNofrGDFJOywRU1Cme+oz8ARwBjgLk4y4OZ+OkHLFYbNLivrJ9rDyk7ADU6rugD4G5V/bfb9SQzEbkB6Kiqv3G7Fj8SkR44n9VOFv6OVG5x/QDnLvwZbheSAvphHfP7TFW/wLmQ1MftWrwilYPrRuBOtelwE8E65lvPThcbSMngEpF+wFHAP9yuJdlF7wPtAqx0uxafs+BqICWDC6e1dU/0yqJpW0cDyzUFV62OszeBE0XErniTgsElIocCw4FH3a4lRdhUzXGgqsXAKuBYt2vxgpQLLmAS8LDaSjyJYh3z8fM6droIpFhwiUhn4HxsqadEso75+LF+rqiUGscVnYEzU1WvdruWVBDtj9kB7KeqlS6X43sikoWzTFmXVJ8uPGVaXCJSCFwK3ON2LSmkN/C5hVZ8RKe5mY9Nc5M6wQVMBF5W1bVuF5JCbEaI+LN+LiDodgGJICI5wFXsXnTTJIb1b8XfHOBpt4twW6q0uC4BPlDV/7pdSIqx4Iq/xUB7EenidiFuSvrgEpF0nCEQd7hdSyqxVavbRnSam7mk+Oli0gcXMA74QlXnu11IijkU2KKqO9wuJAmlfD9XUgeXiKTh3N5jra3Es4GnbWcOMCI6NVNKSurgAkbirJM42+1CUpD1b7WR6DQ35aTwNDdJG1zR30Y3AXfY5GuusOBqWyk9ij5pgws4EWgPvOB2IanGVq1OiJTu50rm4LKJAt1zEFCnqjZ3f9tJ6WlukjK4ROQY4BjgKXcrSVnWMd/GUn2am6QMLuAG4F6bKNA1dpqYGCnbz5V0wSUiB+Os/PuI27WkMAuuxEjZVa6TLriA64CpqrrT7UJSmM16mhjvAn1EJN/tQhItqW6yFpFOwAXA4W7XkuyiVw77AycXFhb+IC0trbuqBqIzy7YHjhKR9dGpWEwbUNUqEdk1zc3/c7mchEqqiQRF5H+AkKpe6XYtySoaWBPy8/N/n5WVdeC5554bHDJkSGbPnj0JBAJs2rSJmTNnRj799NOypUuXCvBoZWXlH+3Wn7YRXWy3i6r+2u1aEkpVk2IDCoBioLvbtSTrBnTJy8t7+/DDDy979dVXtb6+XpuzatUqveiii6pycnKKgbPcrj8ZN2AAsMLtOhK9JU2LS0RuAo5Q1Z+6XUsyEpHeOTk586699trwlClT0tPT02N+7ty5c7ngggsqysrKJldVVd3bhmWmnOj9uJuBY1T1a7frSZSkCC4RyQa+AEao6jK360k2ItI1Jydn8dSpUwt+8pOf7NONvWvXrmXIkCEVxcXFv66trZ0W7xpTmYg8B7yqqk+4XUuiJMtVxYuBjyy04k9EJC8v71833nhjeF9DC6Bbt2689dZbORkZGQ9Eh6yY+Em58Vy+Dy4RCeIMgfiL27Uko0AgcFHXrl373nTTTa2+At2rVy9uvfXWzLy8vOnxqM18K+WmufF9cOFMFLhOVT9wu5BkIyISCoV+/7e//S0UDDaeW6rKH//4R4YOHUq7du0oLCxkzJgxrFu3rtH9r7nmmkBGRkbf6G1ZJg5U9XOcaW56u11Lovg6uKK/YW7EWltt5YSCgoL9hg1reo2RtWvXUlRUxG9/+1vef/99XnvtNVauXMn48eMb3T8YDHLVVVdlhsPha9qo5lSVWqPo3b6s2cpLwWfjzGkubteSjFsgELjt+uuvr9W99Je//EXT0tK0rq6u0a8vXbpUw+HwBrffXzJtOCu0z3K7jkRtvm5xEZ2WWaP/cya+8vPzfzB48OC97tuaN28evXv3JhAINPr1Xr16UV1dvZ+IFLS2RvOtuaTQNDe+DS4ROQHoCPzb7VqSVSQS6XbwwXt3AfC2227jzTff5KGHHmpyn2AwSKdOnSpx5u0ycaC7p7kZ4nYtieDnexVvAu5S1Tq3C0lWqhpoqlN+T7W1tVx55ZW8+OKLzJ49mxNOOKHZ/aOtMT9//rxoVz/XO24X0tZ82eISkb44U6c86XYtfiUimS3tEwgEdm7btq3FY23dupURI0bw/vvv89FHH3HiiSe2+JySkpIAUBJTsSZWKTOey5fBhTNR4H1qMw/sExE5BKgSkWIReU9E7hORn4nIQBEJ7dqvpqZmweLFi5s91tKlSxk0aBD5+fm8//779OjRo8XX37p1K6WlpQHgy9a9E7OHd3Fm5Uj6aW58F1wi0hM4HZjqdi0+9gVQBrQDhgJXAX8H3gB2iMhWEVlVVlaWMWfOnLKmDvLSSy8xdOhQTj75ZKZOnUp5eTkbN25ky5Ytzb74e++9RygUWq7OqsyIY38R6SMip4rIT0XkBhF5IBqyJgbRX+TzgR+4XUubc/uy5j5c9n0I+LPbdfh9A/4FaDNbPfCLrKysyk2bNmljevXq1ehzBw4c2Oj+u5xyyimlOBdVtgE7gDqgCtgZ/XcpUAtEcG4edv375ZcN52zkfrfraPP36XYBe/mf0iH6Ye/gdi1+33DuOChvIrR2AseqKuFw+Ombbrppr8dyNWX58uWalZVVBhwfDavmwnO1298nv22kyDQ3vpodQkT+AuSp6q/crsWPRKQDzqDdc4BhQC7Q8P42xekwP1lVP4k+p0t2dvbK+fPnh/r27duq16+rq6N///7lK1asuK62tvZhETkap18mt5HdK4AbVPXvrXrRFJMq09z4po8r2uH4C+But2vxi2jfUW8RuUlEPgA+BU4DngO644TGLhGc1uyxu0ILQFW/rqmpuWrMmDHlsVxhbIqqct1119V8+eWXS+rq6qZGH/sEGIkTUnvKBjqLSNd9ftEUpE6/4VzgFLdraUu+CS7gCuA1Vf3C7UK8TETSRWSYiNwLfAa8CnQGbgEOUNVxqvq0qm4DnsY5XawDNgEDVfXTPY9ZX1//+ObNmx8bOnRo+YYNe7/Ga319PZMmTap57LHHvi4tLR2lDZr5qjoPuBCobPAUBeYB+cAiEfk/EbkgliEcBkiFYRFun6vGeN6eDWwAjnK7Fi9uONNWjwf+idNq+giYDPSlmfs4gQNxOuG/ADq38BqSk5Pz+3A4XDF9+vRIJBLRWKxcuVL79etXlpeX9xGwXzPHv4zdfW47cSaF3PV/fyHOD+NW4H7gaLe/517egJ7Rn5ekvYfX9QJi/I+4HHjZ7Tq8tAE9cIYxzIn+oL+CcyrdbAA1cpxfAvvvxf4Dw+HwmkMOOWTn1KlTI+vXr9c9lZWV6ezZs/Xss88uy8rKKk9PT78aSIvh2JNxOuy3NLZ/9D3fCqwDioBfAYVu/194cQPWAH3crqOtNs93zkcnClwF/ERV33O7HrdEO10H4nSsn4Nzn+YrwEvA66panuBaTikoKLi2qqrqhMzMzLSuXbvWBAIBiouL0zZs2JAVDoc/KykpeSgSiTylMa5xGZ2m6H+Ar1X1/mb2C+D04VwCnIFzOjwNmKvRsWGpTkQewbm6eJ/btbQFPwTXBGCiqrZ8H0mSic6lfwpOUI0CtuME1UvAh6pa72J5wLdh0xXnhukAzjisFapak6DXbwf8CPg5zinz48ATqro2Ea/vVSJyPvAzVT3b7VragqeDK/pDsRi4SVVfdbmchGgwZGEUMBxYiBNUL6vqajdr8zoR6Y/TCpuAcyo5DZipKXhrmIjsB3wOtE/UL5FE8npwnYUzu+kx6uVCWyEazkey+xTwCOD/cMLqNXWu/pm9ICJZwLk4IdYPeAb4X1Vd7GZdiSYiHwO/UdWkmy3C68E1D3hYVZ9xu5Z4EpF04AR2h1WQ3aeAbyfjb0i3iEh34CKclaCKcVphT6vqdhfLSggRuQOoVtXfu11LvHk2uETkeOAfwGGaBHNuRQfQnoETVGfiXPXZFVZLkrVF6RW7LijgtMLOBF4D/pck7tAXkRHArap6vNu1xJuXg+tlnEUuH3a7ln0V/W0/CieshuBM8PYS8IqqrnextJQW7dC/ECfE2gFPAI8nW4d+9OLOZqCLqibV3GeeDC4ROQqYDfTwU8dqE0MWZrF7yEKTU8QYd4hIP3Z36C/COZWc4afPXXNE5HXgAVV9ye1a4smrwTUdWKaqd7hdS0v8MGTBtCzaoT8GJ8T640z7M01VF7pZV2uJyA3Agap6ldu1xJPngktEegAfAz292ryNDlkYiRNWNmQhyYhIN3Z36G9nd4e+767wisgA4B+qeqTbtcSTF4Pr70Cpqt7kdi272JCF1BQ99R+OM7j1TOA/OCE2xy8d+tG7DDbh3N+ZNP2qngquaEtmJXCEqm50uRYbsmC+Fe3Qn4BzKrk/uzv0v3SxrJiIyHM4i8UmzeIyXguuPwHtVPUKl17fhiyYFonIMTgBdiHOnR3/izNCv7KZp7lGRH4BnKiqP3G7lniJa3CJyEac6ZXbwiZV7Rjvg0aHLOzqWLchCyZm0Q790TghNpBohz6w0Eu/5KILzLyHM3OIZ+pqjXgHV5t+U1RVWt6reS3MsjDHhiyYfRGdqfUinA79EnZ36Be7WdcuIvI5MEpVl7tdSzykRHDZkAWTKNFfjMNwWmEjcS7i7OrQd+1zlmzT3CRtcNmQBeM2ESlkd4f+Aezu0E/49OPJNs1N0gRXgyELu26xORIbsmA8Irqi0a4O/SU4rbAXE9WhLyL741xs2l9VaxPxmm3J18FlQxaM30QX/NjVoT8IeBbnqmSbd+iLSBEwBagB2uPM778JWKSqO9rytePNteDq3bs3t9xyC8OGDSMcDrNmzRrGjBnDZ5991uRzVFVsyIJJFtEO/Z/hdOiXsrtDf2ucX6d9enr6pSJyUyQSye7bt295ly5d0gDWrl0bWbFiRU5WVtY3ZWVl90UikSd8EWKxTk4fy0bzqxJ/u40YMUI3b96s1157rfbp00d79uyp5513nubm5rb03NdpxcIQttnmxQ1nmcBhwHScqa+fw/nlHGjlcQMZGRnXZmdnl//oRz+qmD9/vtbX1+ue6urq9O2339Yf/vCHZdnZ2WWBQOBSPL5CULz/A1oMrXbt2umWLVv0zDPPjCnk9tjGACG3v2m22dZWG868+VcAC4CvgNtx7tvd2+MUhsPhBQMHDixbvXq1xuqTTz7RXr16leXl5c3x8s9awheEveKKK1i3bh3Dhw/nyy+/ZO3atdx///1kZra81qeqztQErmZjTKKp6g5VfVhVB+FcFQ8DH4rIXBH5cXRoT7NEJD83N/f9n/3sZ0d9+OGHoUMOOSTm1+/bty9LliwJnXXWWceHw+G5sbyeK+L826LFVtPChQu1pqZGH3zwQR04cKD++Mc/1urqar3ttttafK7bKW+bbW5sQCZwPs6srduAh3E69r93OgdIOBye9fOf/7wq1kV7G1NXV6ejR4+uyM3Nneb2+29sS2jnfDAYpKKigpkzZ3LBBRd8+/izzz5Lp06dOOmkk5o9vsZh5LwxfiYiB+F06F8ClOF06E/XaIe+iFzQrVu3aStXrgxlZWW16rV27NjBIYccUlFcXHyOqr7R2trjKaGnigUFBaSnp/Pvf//7O49HIhHKy+0M0JiWqOpXqvpH4BCclcwHAJ+JyPMiMjIUCj3w1FNPNRtaqsrdd99Nz549yc3NZfjw4Y1ezS8oKODRRx/NycvLe7Ct3s++SmhwlZSUUFtbizNW1JGXl8fw4cN57bXXElmKMb6mqhFVfUudGR+6AXOAezIyMgpOPLH5tZMnTpzIfffdx4MPPkhRURGBQIALL7yw0X1Hjx5NVlZWFxEZFPc30RpxPhdvsZ9qxowZunDhQj366KN18ODB+tZbb+mSJUs0KyvL+rhss60VW35+/vMPPvigNmfu3Lmalpam8+fP//axN998UwHduHFjo8+55ZZb6rKzs//m9vtruCX8quJll13GypUreeONN5gxYwarVq1i2LBhVFUlxdoExrhpyHHHHdfsDnfffTejRo1iyJAh3z7Wrl07ADZv3tzoc44//vhAdnZ28x3QCRZM9Atu3bq1yWapMWbfiEh6IBDo3Lt37yb3qaysZPbs2Tz++OPfeXxXYBUWFjb6vGOOOYaKiopD41dt6yW8xWWMaRPZwWAwkpGR0eQOy5Yto66ujv79+3/n8Y8//ph27drRpUuXRp+Xl5dHbW1t6y5RxpkFlzHJoaauri5NtekRSdu3bwegffv233n8xRdf5Jxzzmn6wDU1pKWleWrOOgsuY5KAqlZlZmaWfvnll03uc8ABBwC7Awxg1qxZFBUVcc011zT5vJUrVxIKhTw1jbkFlzFJIisr65OPP/64ya/37t2bbt26ceutt7JmzRqef/55JkyYwOTJkzn66KObfF5RURGRSGR+W9S8z+J5iZK9v2l6rza3L8HaZpuXt7S0tCtHjRpVps34+OOPddCgQZqdna1HHHGEPvbYY83trqqqgwcP3gmMc/v9Ndx8PZGgMWY3EcnPysrasGbNmuzOnTvH5ZjLly9n8ODBJRUVFQeohybmjPep4qY4Hy9RxzbG91S1JBAIPPKrX/0qLtNBqyq//OUvyyORyB1eCi2Ic3CpakdVlTba4r6mojHJpry8/OY5c+YUP/vss60+1kMPPRRZsmTJ2qqqqrvjUFpceWola2NM64nIgJycnLdnzpwZOvXUU/fpGM8995xefPHFOysqKgar6qo4l9hqdlXRmCSjqkUVFRVnjB49uuzOO++sr6+PfQhWTU0NkydPrr344otLKioqTvZiaIG1uIxJWiLSMy8v79nu3bsfcfvtt4dGjhxJIBBodN+amhpmzpzJzTffXL5ly5aPdu7c+WNV/SbBJcfMgsuYJBZdWfvC/Pz8m4LBYPfhw4fr8ccfH+rUqROqytdff817771XNnfu3ICIrNixY8efgBnq8WCw4DImRYjIkcBxubm5x2VmZh6oqlpdXb2uvLz8A+B99dEK7xZcxhjfsc55Y4zvWHAZY3zHgssY4zsWXMYY37HgMsb4jgWXMcZ3LLiMMb5jwWWM8R0LLmOM71hwGWN8x4LLGOM7FlzGGN+x4DLG+I4FlzHGdyy4jDG+Y8FljPEdCy5jjO9YcBljfMeCyxjjOxZcxhjfseAyxviOBZcxxncsuIwxvmPBZYzxHQsuY4zvWHAZY3zHgssY4zsWXMYY37HgMsb4jgWXMcZ3LLiMMb5jwWWM8R0LLmOM71hwGWN85/8DYpQ3w63ka30AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "nodes_of_dag = set(range(7))\n",
    "dag = DAG(nodes_of_dag)\n",
    "dag.add_edges(parents_set={5, 0}, target_node=1)\n",
    "dag.add_edges(parents_set={6, 0}, target_node=2)\n",
    "dag.add_edges(parents_set={5, 2}, target_node=3)\n",
    "dag.add_edges(parents_set={6, 1}, target_node=4)\n",
    "observed_set = {0, 1, 2, 3, 4}  # latents set is {5, 6}\n",
    "fig_dag = draw_graph(dag, latent_nodes=nodes_of_dag - observed_set)  # plot the DAG"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Now, instantiate a perfect oracle, set it to count the performed CI tests, and turn on caching."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "ci_test_dsep = DSep(true_dag=dag,  # the DAG from which to read-off the true conditional independence relations\n",
    "                    count_tests=True,  # count CI tests per conditioning set size\n",
    "                    use_cache=True)  # make sure the same CI test is not performed and counted more than once"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Then, learn the casual structure using the perfect oracle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS4AAAEuCAYAAAAwQP9DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtc0lEQVR4nO3deXxU9b3/8dcnk0ySyUJAUFQEsUVFCcockE1FrTvKrVBbaaXWpb1qvT+3lurVW0urVqu11OJSb9tbtagtSou4ldYFNQGFM4AIuMtWpSACWWayTT6/P05AlkzWSU5m5vN8PM4DyZw585kkvP1+v/M936+oKsYYk0qy/C7AGGPay4LLGJNyLLiMMSnHgssYk3IsuIwxKceCyxiTciy4jDEpx4LLGJNyLLiMMSnHgssYk3IsuIwxKceCyxiTciy4jDEpx4LLGJNyLLiMMSnHgssYk3Ky/S7AGJN6RKQXcBAQALYBn2g3rkoqtgKqMaYtRGRkcXHxNcApNTU1ffv161cTCAR027ZtwYaGhvq8vLxl27ZtmwX8TVXru7QWCy5jTEtEZGhxcfHs3NzcI6655prcSZMmBY488kiys70Om6ryySefsHDhQmbOnFm5atWqhmg0egXwl65qhVlwGWMSysvLuzoQCPz8zjvvzL3iiiuyAoFAq89ZtGgRU6dOrd62bdsrFRUV31DV6mTXZcFljGlWQUHBT/v163fdSy+9VHDYYYe167k1NTVceumlNfPnz19TWVl5QrLDy4LLGLOPrKysqQcddNDvli5dGurfv3+HrtHY2Mi0adNqnnnmmRd37NhxTjLrs+AyxuxBRPqHQqH3Fi5cWDRy5MhOXaumpoYjjzyyet26dZeo6l+SVKLN4zLG7KmoqOjOyy+/PK+tofXHP/6RUCjEBRdcsM9jeXl5PP744wWhUOhBEclJVo0WXMaYXUSkV319/dd/+MMftilkZs+eTW5uLrFYjDFjxjR7ztixYznqqKOyga8mq04LLmPM7qZ85StfaWjLuNb8+fPZtGkTBx54IEDC4AK49tpri3r37n1Vsoq04DLG7FJUVDThjDPOKGztvFdffZV58+Zx/fXXs2zZMoLBICNGjEh4/oknnkhNTc0IEZFk1Gm3/BhjdsnOzh7dUgABLFu2jBkzZvDss8/u+vuIESPIzc1N+JyDDz6YbG/G6oHAJ52t01pcxphd4vF4r379+iV8/L333uM73/kODz/8MHl5eYAXXC11EwFEhJKSknqgdzLqtBaXMWYXEYnH4/GEj991112sXLmSQw89dNfX4vE4q1ev5vXXX2fp0qUJn9vY2ChA4ou3p06bx2WM2alv375vPvLII6POPvvsZh9fu3YtVVVVu/7+9ttvM3XqVObNm8exxx7LwIEDm31eQ0MDBQUF9XV1dX1VtaKzdVqLyxizS2Vl5cIlS5Y4Z599drPDSLu3tMDrJgYCAc4880yCwWDC665Zs4a8vLwttbW1nQ4tsDEuY8xu6urqXvrLX/7S5vsKV69ezZe+9KUWQwtg/vz5cVV9sdMFNrGuojFmFxEJFBQUfPrKK6/06+ztPjvF43EOPPDA6JYtW05UVTcZ17QWlzEGEe+DP9DBtbUn/+6aa66pTlaj5oEHHmisra19J1mhBdbiMiZjiJADzAC+DOwP7Ic3PaEXEAIagWzQdUVFvaK/+MUvjrj88ss71bj58MMPGT58eCwajYZV9Z1OvoVdLLiMyRAiZAPr8SaBJlIHHAsi+fn5b8ydO7fwzDPP7NDrbdq0idGjR0c3b958YywWu7dDF0nAuorGZAhVGoDLgGiCU2LAnaqsUdXVsVjsjMmTJ1c98MADje1t4CxduhTHcaJbtmy5K9mhBXjrRdthhx2Zc4D+E7QBVPc6PgYN7nkuQ4uKilaNHz++6s0339TWfPLJJ3rdddfV5efnV4jIBV31HqyraEwGEeEA4H7gPGD3G55jwCmqLN73OZIdDAb/XzAYvGHgwIF5559/fsGoUaOyBg8eTFZWFp999hmRSIQFCxZUvfjii4GcnJw5lZWV01X13132Piy4jEl/TeNb3wduBv4PL7SuxBuUrwEeVuXylq8hJwA3B4PBFUVFRRMaGhoOUtVAIBDYUV9f/0ZVVdWrwFOquqNr340FlzFpT4QJwCxgE/BfqrwjQj6wFu/TxS3Al1SpTHwNOQJYDqCq+V1dc2tscN6YNCXCQSI8BjyKNw3idFXeAVAlBnyv6dSLWgmtLwPlQB4QFBELLmNMcokQFOGHwFvAR8BQVZ5UZe/u1dOAo8rzia8lhwGL+GI5mmpgSBeU3S52k7UxaUSEU4HfAB8DY1V5P9G5TUEWSXwtORRYDPThi4F8BY7AC0XfWHAZkwZEOAS4B3CAa4D5zbSw2nE9GYgXWvuxZ8+sADiy45Umh3UVjUlhIuSK8N/AMmAVcLQqT3cytA7CC62+7JsRAaDltZ27gbW4jElRIpwF3AusBo5T5aMkXboUKMSbYV/AvuF1dJJep8NsOoQxKUaEwcCv8ALkalWeS/5rSDZwAvBN4DtAA5CD1+KqAULqY3hYV9GYFCFCvgi3AEuAN4HSrggtAFVtUNWXgel4La9z8Gbcb8GbFlHcFa/bVtbiMqaHE0GAc4GZgAtcr8r67nltuQj4qqqe1/R3AQ5R1W55/URsjMuYHkyELwO/Br4E/Kcq/+jmEqYAc3b+pal76GtogXUVjemRRAiJcCvep3uvAMO7O7REpAg4CZjfna/bFtbiMqYHaeoWTsabk1UGHKPKv3wq52ygTFW3+/T6CVlwGdNDiHAk3vSGg/DuH3zF34qYAjzlcw3Nsq6iMT4ToUiEO4HXgOeAEX6HVtON1GcA8/ysIxELLmN80rSzzgV4E0j7401vmKlKvc+lAZwORFR1i9+FNMe6isb4QIRheDdDlwAXqFLmb0X76LHdRLAWlzHdSoReItwDvAQ8CYzsaaElIkG8Cad/9buWRCy4jOkGTd3CacAavFnnR6tynypxn0trzinAu6rq16eZrbKuojFdTIRj8ZZOzgW+qsqb/lbUqsn04G4iWIvLmC4jQm8RZgF/Bx4BxvT00BKRAPBVLLiMySwiZIlwCV63MAAcpcpDPbRbuLcTgI2q+rHfhbTEuorGJJEII/G6hQqcrZp4aeQeagow1+8iWmMtLmOSQIT9RPgt3n19DwLjUy20RCQLb6PYHt1NBAsuYzpFhIAIl+NNIq3B21Hnj6o0+lxaRxwHVKjqGr8LaY11FY3pIBHG4HULo8Bpqv7ufJMEPXrS6e4suIxpJxH2B+7Au5dvOvBYZzan6AmaFgic0nT0eNZVNKaNRMgW4SrgbeBzvG7h7FQPrSbHNP253M8i2spaXMa0gQjHA/cBW4GTVFntc0nJNgV4ys8NMNrDgsuYFohwIPALvJVArwfmpEkLa29TgEv8LqKtrKtoTDNEyBHhOryt5jfidQv/ko6hJSJDgV7Qs2f1785aXMbsRYST8T4t3Agcr8q7PpfU1SYDc1U1ZaZwWHAZ00SEAcDdwBjgWuBv6djCasYU4Dq/i2gP6yqajCdCUIQf4X2i9h7evYV/zYTQEpHDgAF4y0anDGtxmYwmwul4K5G+B4xW5UOfS+puk4F5qpoKN4DvYsFlMpIIg/C2ADsWuFqVZ/ytyDeTgZ/6XUR7WVfRpB0Rhotwe4LH8kS4GW8r++V4K5FmZGiJyMHAkXjLSKcUa3GZtCJCDt5a7oeJ8IoqC3Z7bCLedvZv4a31vtafKnuM84BnVLXO70Lay4LL9CgiUozXCggBtcB7qrq1HZf4Ad6GqgHg9yJ8CW/weSZwBPB9Vf6e1KJT1xS8IE85kiIz/E0aE5FBeXl53w8Gg9+MxWL7H3roodHCwkJqamr4+OOP87OzsyuAeVVVVb9S1VWJr8OX8VpT+U1figKL8e7DuwuYqUptV7+fVCAi/YAPgP6qGvO7nvay4DK+EZHiwsLCe1X1GxdffLFcfPHFucOHDyc7+4uOQGNjI++99x6zZ89umDVrVn1jY+OrFRUV31HVTXteC8ELKQevtbVTHG+t96Xd8Z5ShYhcBpyuql/3u5YOUVU77Oj2A3AKCgo2T5s2LbZ9+3Zti5qaGr3hhhtqQ6FQBXDWntfTS0GrQHWvox70ab/fb087gOeBb/hdR0cPa3GZbicio0Kh0EuPPPJI4ZQp7V/+qby8nLPOOitaUVHxTVWdJ0J/4H2gMMFTosA5qrzcibLThoiUAOuBg1W10udyOsSCy3QrEemXn5//3hNPPFEyadKkDl9n6dKlTJgwIRqNRkeBPoLXRdxJ8Qb264FGIA94Q5UJnak9XYjIhcDXVbXjPwCf2aeKplsVFxf/7uKLLw51JrQARo4cyV133ZU3ffrNf6mu1gNAVgEfAp8Cm4DP8NbO2nn06O22ulnKLNGciLW4TLcRkTH777//P9euXVuQn5+/z+Oqym233cZzzz3HO++8g6oyYcIE7r33XgYOHLjP+Y2NjRx33HFVruteoap/6o73kOpEpBD4BDhUVT/3u56OspnzptsUFxdfP3369PzmQgtg3bp1uK7LddddR3l5Oc8//zzvvPMOF1xwQbPnZ2VlcfPNNxeWlJRM78q608xZwOJUDi2wFpfpJiKSm5OTU7Fp06Zgnz592vy8O+64g5tuuom6ujoCgcA+jzc0NNC3b9/Yjh07jtYevvtyTyAijwOvqOpv/a6lM6zFZbpL6YABA2raE1oAr776KkcffXSzoQWQnZ3N6NGj64FRSagxrYlIHl6L628+l9JpFlymuwwfOXJk8+mTwE9/+lNefvll7r///hbPGzduXGEwGHRaPMkAnAasUNV/+11IZ9mniqa79Np///1z2nJifX09V111FXPnzmXBggUcf/zxLZ6/3377ZQWDwf2SUmV6mwLM9buIZLDgMt0lXl9f3+qa5p999hlTpkzh888/580332Tw4MGtXrihoYHGxsb6pFSZpkQkBzgX+B+/a0kG6yqa7rJuzZo1LS6fsnLlSkaNGkWvXr0oLy9vU2gBrFmzpiYajb6flCrT10nAB6q6we9CksGCy3QXd9myZcFEn2I//fTTjBs3jpNOOokHH3yQ6upqNm3axJYtW1q9cHl5eS3YTdStSPlJp7uz6RCmW4iIFBUVbXzhhRcOGjdu3D6PDx06lHfeeWefr48cOZIlS5YkvO7GjRsZMmRIrKampp+qVie16DQhIgHgX8DxqvqB3/Ukg7W4TLdQVY1Go7+85557os09vmbNmmZXAWgptAAefPDBhuzs7NkWWi0aB/w7XUILrMVlupGIjAoEAmXPP/98zmmnndbp661fv56jjjoqVl1dPUJV033T1g4TkZnA56qacptiJGItLtPlRCQkIrcBz8fj8Se+973vVcdinVt0s7GxkWnTplXH4/E7LbQSExHB28knbca3wILLdCHx/AewGjgMGA5ctHXr1mfPOeecaE1NTYeu29jYyHe/+93aZcuWra6pqWl2Nx+zy0gghvczSBsWXKZLNO2Q/AxwB3Cpqk5V1U9UVSsrKy9csmTJi+PHj6/+6KOP2nXdLVu2MGnSpNicOXNWVVZWnqqqNn+rZVOApzTNxoQsuExSiUieiNwCvAm8Chyjqi/ufo6q1ldWVp63atWqn5WWlsZ+/vOfN37+ecuLFVRXV/PQQw8xZMiQ2CuvvPK/lZWVx6tqRde9k9TX1E1Mq2kQO9ngvEkaETkLbzv7FcC1qrq+Dc85slevXrfW1tZOPO200xomTJhQWFpays5dflavXk1ZWVl0/vz5WTk5OeXbt2+/SVUXd/mbSQMiMhx4Ghicbi0uCy7TaSIyEG/fwuHAVar6Qgeu0Q84u6CgYGxeXt5IVQ2JSG1dvO6t/W7YL2/tDWunq+q6JJee1kTkJ0Cxql7ndy3JZsFlOkxEgsB1eJuw3gv8QlU7NuKegBNxJgCvAFe6YfeBZF473YnISuByVS3zu5ZkszEu0yEicgpel/B44DhV/WmyQ6vJWLzNL37lRJyWl4kwu4jI4cB+wCK/a+kKFlymXUTkoKZVNP8A3ACcq6rt+2iwfSYCAuQCzzgRZ1AXvlY6mQL8VVVbXZEjFVlwmTYRkRwRuQ5vi/uPgKNUdV5XDvo6ESfAntuOFQIvOhEn0f6J5guTSZO1t5pjwWVaJSInAC5wJjBOVW9S1WbvOUyy4UDDbn8PAAcDTzoRx353ExCRQcChwEKfS+ky9sM3CYnIASLyMPAY8DPgDFV9rxtLOBHYe9XUPOAE4LZurCPVTAaeVtWGVs9MURZcZh8iEhCR7wNvA5vxuoVzfJgLdA5eUO0tBEx3Ik7K7sTcxdLu3sS92dLNZg8iMhq4H6gCTlLVVX7U4UQcwftEcSfF6zY2AtuBMrxdq81uRKQ/MAx4sbVzU5kFlwFARPbDu69wIjAdmO3zbOvBQAFQD3wALMfrOo50w+4mH+vq6c4DnlPVWr8L6UrWVcxwIpIlIt/FWz0ghtct/FMPuEVkHd7KBr3csHsUMA0owgsyk1ha3pu4N5s5n8FEJIzXLVTgSlVd5nNJLXIizj+AX7th9xm/a+mJmlrNHwEHdtOnvr6xFlcGEpESEZkFPA/8LzC+p4dWkzJgvN9F9GCTgH+me2iBBVdGaVrY79vAGrzxzaNU9fcpNLvagqtlGdFNBOsqZgwRKQXuw5tKcKWqvulzSe3mRJxi4BOgjxt2W9yjMdOISDGwERiQCeuUWYsrzYlIkYj8Eu/j8ceB0akYWgBu2K3A+4RxhN+19EATgVczIbTAgittNXULL8DrFvYBhqnqA6oa97m0zrLuYvOmkMb3Ju7NgisNiciRwD+AG4ELVPViVd3sc1nJYsG1FxEJAacB8/yupbtYcKURESkQkZ8Dr+NtVOGo6us+l5VsZcD4ppn1xnMGsERVt/pdSHex4EoDTd3C8/AmkQ4ChqvqzDS9yXY9EMebWW88GfNp4k4WXClORL4EPAvcDlysqt9U1U98LqvLuGFXse7iLiKSizcw/zefS+lWFlwpSkTymzZDeANv3aVjVPUlf6vqNhZcX/gKsEpVM+qGcwuuFCQiZ+MtOTMMGKGqd6pqJs1rsuD6QtovYdMcm4CaQppWtpyJF1hXqerf/a3IH07EyQE+Bw5xw+52n8vxjYhk4y3tMzLTtm6zFlcPISKHtPBYroj8N97yyRGgNFNDC8ANu/XAUmCM37X47ERgXaaFFlhw9QhNi/eta9oJeu/HTsXboGIsMEpVf9ZF24ClGusuZmg3ESy4fCciOXhrugvwexHJa/r6wSLyZ7zVG36gqueq6sc+ltrTlJPBwSUiWXiLBmZkcNkKqJ0gIkcAY0Oh0Oi8vLwBgMZisXWxWOwNoLyN+w1OBw5o+u9ewP+IyOd4s94fwJvikPbLlHTAImCUE3Gy3bCbjvPVWjMG+LybNy/pMWxwvp1ERICvl5SU/Dcw5KSTToqPHz++8MADDwRgw4YNWlZWVrVw4cKcQCCwYvv27bcD85tbUbRpDtZKIH+3LyvezPfLMvWXsq2ciPM2cJEbdl2/a+luTTfOV6nqLX7X4gdrcbWDiBxaVFT0+CGHHDLs1ltvLTz33HPJzt7nWyhAUW1tLXPnzh190003PbZ169ZyEfm2qm7a7VoCPAIE93q+4s0Mf78r30ua2DnOlVHB1fS7Mxn4D79r8YuNcbWRiIzNz89/68Ybbxy5YsWKwvPOO6+50NolNzeXqVOnsmbNmoIrr7xyQigUWt20JtZOFwLH4G1yurssYBTeapamZZk6zjUCb8ejlX4X4hfrKraBiBwbCoVee/LJJwvPOmufD/7a5PHHH9fLLrtsRzQaHYW3vdZavF1sEtkCDLRPEBNzIs6XgZfdsJtwKkk6EpHbgICq3uB3LX6xFlcrRCSvsLDw6YceeqjDoQUwdepUuf3224uLior+irf8SAHerjoxvL0CY8AmvBulFwLPse8uzmZPHwJBJ+IM9LuQriIivUVkoIj0b9qoV8iwtbeaYy2uVhQWFt59yimnXPH000+HOnutxsZGxo8fX71kyZLZ8Xj8GbxW1xZgq6ratlsd4EScucAcN+w+7nctySAiQWBy7969L6urqws3NDQUhkKh+oaGBqmpqcnKz89fW1FRcSBwZKbdn7gHVbUjwQEU5ebmRjds2KAtufbaa/WII47QgoICLSoq0pNPPlmXL1/e7LkrV67U/Pz87UCO3+8vHY6wG/5B2A3P8ruOzh6ABAKBy/Lz87ePGTOm4tFHH9UPP/xQGxsbdaetW7fq3/72N504cWJNfn5+TXFx8Z+AEr9r9+X75XcBPfkQkSsmTpxYpa349a9/rWVlZfrxxx/rkiVL9LjjjtMhQ4YkPH/UqFEVwPktvTZoAPRU0Fy/vw89+Qi74bFhNxzxu47OHMB+xcXFrw4dOrQqEoloW2zfvl2/+93v1oRCoa3ASX6/h27/nvldQE8++vTp89Kf//xnba/bb79d+/Xrl/Dx+++/X3v16jW7udcEDYFeBboJVEFH+/196MlH2A3nht1wVdgNF/ldS0cOoF9hYeFHV199dW1dXZ2214IFC7SwsLAaONvv99Kt3ze/C+jJR0FBwdYPP/xQ2yoej2tZWZkOHDhQf/Ob3yQ874033tCSkpIPd38t0P1BbwetBK1qCq0doMf4/X3o6UfYDb8WdsOn+l1Hew8gUFRUtHz69Om12gmLFy/WUChUjbchiu/vq1u+d34X0FMPIC8QCMR3H2NIZMWKFVpQUKDZ2dlaUFCgc+bMafH8HTt2aE5OTq33OnoE6MOgsaZDdzsqQQ/3+3vR04+wG74z7IZv8buO9h65ubk3jh07tioej2tnPfTQQ41FRUVrMmXs1KZDJJYTCATi3qfPLfvyl7/M8uXLKSsrY9KkSVx++eVs2bIl4fm5ubnE4/GACMvxtg+bBuQ1HbsLAqeLcJ4IZ4hwogijRDhahMNE6C9CL5GMnzZRBozzuwgAEUIiHNf6ebJ/VlbW/8yePbsgKyvxP8M//OEPnH766fTr1w8RIRKJNHveZZddJqWlpYeIyKUdrz512HSIBEQkJxAI1ESj0axgcO+7chJbv349gwYN4tlnn+Xss89u9pzNmzczcODAaG1tzU+Am4HcpmNvjXhriQfwdqDObzpCzfwJ3lyw6F5/dsXXalXpMb84TsTpizenq48bdn3dN1KE/8D7mb0GXKXKW82dl5ube/PXv/71/3700Ufzm3t8p5kzZ1JcXMy7777LzJkzqaysJNHv40svvcR55523tqKi4jBN83/Ydq9iAqpaX1JS8q9Vq1YdMmJE2zdOfv11bzewIUOGJDxnxYoVFBQUvF9Tw10i3At8H5iB9/PYvdWlwDRVWl0doqnVtXuYNRdwzX2tN3BQO87PB3JEkhaErT1W01pIumH3MyfifIq3MuyK1r5XXawAqAKOBxaL8DTwQ1U27H5SXl7eFf/1X//VYmgBXHPNNQBcddVVlJaWJgwtgJNPPpnCwsK+FRUVDt5Ci2nLgqtlbyxevDhhcP34xz9mxIgRDBs2jNraWhYsWMBPfvITrrrqqhaDa/HixY2xWOw1AFVqgXtE+APwY+ByvJ9LDl5Lq023/KhSD9QDXb4FuwgBOhaSBUC/Npy3+2O5ItTSStAVHXd/dmMsf5asZOnej7X0vL0fU6Wxk9+efLw7UqTpvycDk0R4EPipKttFpE9ubm5fx3HafFHXdWntfBHh1FNPDTzyyCOjseDKXDt27Hh41qxZZ1xxxRVFzT2+efNmrr32Wj799FOKioo4/PDDmTVrFt/+9rcTXrOxsZEHHnggFovFZu/+dVW2A9eJMBO4C+8m6+wk/ENKOlXieK2Kqq5+LRGy8FqhLYZk8OB/VeT0+XxY9crhG3Y7r8/e5zX33N3/W4Q6OtdKHMuet2rlNB1XAt8T4Rcw6M2hQ3vHAoFAm8Yg4vE4K1as4OKLL2713DFjxuTPmzdvPHBfW66dqiy4Wvb8unXrasrKyorGj993EYIHH3yw3Rd85plnqK6u/hRvW7F9qLIe+IYIxwBntvsF0kxTcEebjoScyG0rgec//d1l93T0tUQQvLHG9rYk84H+TX8eQfP/rnaOY86A+S/27z+9zTtxr1mzhlgs1mqLC+CAAw4gOzu7f1uvnaosuFqgqvGsrKyrp02b9r9r1qwpyM1tbvy87SoqKrjsssuiFRUVV7U2eKrKCvwfr0kl7wHFTsQ5yA27HdoQt2ksrYY2ds+bI8JPgJF7fbkBrxtfBvwERvRXPXVUW6/pui45OTmUlpa2em7Tr1WPa6Unm02HaIWqPrF169bXr7766rrOfFDT2NjIJZdcUhOLxeZqBu/Q01XcsNtIz1ifa/dhhZ0hOBsYocppqpRB/JO1a9e2+YKu6zJs2LAWB+Z3WrduHbW1tWm/N4G1uFqhqioiF8yePfvN/Pz8Qb/85S+DLc27aU5DQwOXXHJJzYIFC1ZXVVX9ZxeVar5YEXWOjzUU47WwaoDfADNV2bzXOSs++uijUG1tLYla8fX19axatQqARYsWMWDAAJYvX05BQUGLH/yUlZVVVVVVlSfjjfRofs+ATZUD6FNUVOSOHj26qj23Aa1atUqHDRtWVVxc/DJQ6Pf7SOcj7IZPCLvhN/2sAfQc0O+DFrR0XklJyZrnnntOE1m0aJHiTYfZ4/ja176W8Dm1tbXaq1evKDDE759Fl3+f/S4glQ4gkJubOz0/Pz/67W9/O7p06VJt7pageDyu5eXlev7551fn5+dXZWdn/ydNk33t6Loj7Ibzw264OuyGQ37X0toBXHbqqadWahI98cQTWlJS4mtwd9dhM+c7QET6B4PB7wWDwSsbGxtLSktLawYMGJClqqxbt65x1apV+cFgcHM0Gr23oaHhD6q61e+aM4UTcRYDN7hh9xW/a2mJiBTk5+f/68UXX+w1duzYTl+vpqaGoUOHVq9du/ZbqjovCSX2aBZcnSQi+wNhvImVCvwbiFhY+cOJOL8EPnfD7m1+19IaEfnaIYcc8vA777wTCoU6t8Du9OnT637729++XFFRcZZmwD9q+1Sxk1R1s6q+oKqPquqfVPUfFlq+6jE3XLdGVZ/cvn37s+edd160rq6uw9d5+OGH9b777qusqKi4KBNCCyy4TPopB8Y6ESclfrcrKysvXLRo0eunnnpqdNOmTa0/YTfxeJyf//zn8SuvvHJbNBo9QVX/3UVl9jgp8cM1pq3csLsJ2AYM9buWtlDVusrKyomu6846/PDDY7/73e+0pqbl+a+qyhtvvMHIkSOr77jjjmXRaNRR1TXdVHKPYGNcJu04EecR4HU37D7kdy3tISIjS0pK7o3H4yMuuuiirOOPPz547LHH0qtXL+rq6nj33XdZunSpPvLII1UbN26M1tTUzGhoaPitqqb9TPm9WXCZtONEnP8Exrph9zt+19JeTsQ5PB6NP/b2SW/P69Wr14T6+vrShoaG/KysrIbc3Nz1sVjstVgs9hzw90wMrJ0suEzacSJOKTDXDbuJp5j3QE7E+RIQwZt9n++GXdvFPAEb4zLpaBXQz4k4B/hdSFs5EWco8CbevY61QNruzp0MFlwm7TTdcL2IFJkW4UScY/Dq7Y23AGEtMMjXono4Cy6TrlJiPpcTcY7DW5++F15ogbf4gQVXCyy4TLrqCUvctMiJOCcAL7HnUjjgLUg4uPsrSh0WXCZdvQEc40Scvbd86xGciDMGeBFvHf69CXBU91aUWiy4TFpyw2413p6Ve69G2lPE8MJ150Ygezuse8tJLRZcJp2V00PHudywu8INuyfgfXp4I1CHF2I7b1o82K/aUoEFl0lnO1dE7bHcsLsZeB3YiLcv5J3AJ0ChE3HavKFGprHgMumsDBiXAgFwIfCYG3Y/cMPuj4EBwCA37Nrs8ARs5rxJa07EWQec7obdd/2upTlOxAkAG4BT3LD7jt/1pAprcZl019Pnc50MfGKh1T4WXCbd9fT5XN/C277MtIMFl0l3PXaA3ok4+cBXgSd8LiXlWHCZdLcSONiJOPv5XUgzzgWWumH3U78LSTUWXCatuWG3AW+iZ+e30km+bwF/8ruIVGTBZTJBjxvnciJOH2AC8Fe/a0lFFlwmE/TEca7zgb+7YbfC70JSkQWXyQSLgbATcYJ+F7KbC7FPEzvMgsukvaZWzQfACL9rAXAiziC8XYhe8LuWVGXBZTJFTxrn+iYwxw27Hd8FNsNZcJlM0SPGuZrum7RuYidZcJlMUQaM7wE3XB+Dt3hguc91pDQLLpMp1gFx/F8S+Vt4K0Fk7J6IyWDBZTJC0xIxvo5zNa0EMRWbdNppFlwmk/g9znUisMUNu6t9rCEtWHCZTOJ3cNlKEEliwWUyyXJgsBNxSrr7hZt2G5oMPN7dr52OLLhMxnDDbj2wFBjjw8tPBJa7YfdfPrx22rHgMpnGr+6irQSRRBZcJtN0e3A5Eac38BXgqe583XRmwWUyzSJglBNxsrvxNb8G/MMNuzu68TXTmgWXyShu2N0GrMebwd5d7NPEJLPgMpmo27qLTsQ5BCgFnuuO18sUFlwmE3XnONdU4Ck37NZ20+tlBAsuk4m6c69F6yZ2AQsuk4k+BHKdiDOwK1/EiTilQG/gta58nUxkwWUyTtMN193RXfwW8LitBJF8FlwmU3VpcDkRJwtvpVObdNoFLLhMpurqFtfxwHY37K7swtfIWBZcJlNFgMOdiFPURde35Zm7kAWXyUhN0xOWAaOTfW0n4uQCU7CVILqMBZfJZF3VXTwLWOmG3fVdcG2DBZfJbF0VXDZ3q4tZcJlMtggY3bQWfFI4EacXcDrwZLKuafZlwWUylht2twCbgGFJvOwU4KWmm7lNF7HgMpku2d1F6yZ2Awsuk+mSFlxOxDkYGAE8k4zrmcQsuEymKyd5N1xPBf7qht2aJF3PJGDBZTLdu0CxE3EOSsK1rJvYTSy4TEZrugG60ztcOxHnKGB/YGEy6jIts+AyJjnjXN8CHnPDbjwJ9ZhWWHAZ08lxrqaVIKyb2I0suIyBJcDRTsQJdfD544BqYEXySjItseAyGc8NuzFgJXBcBy/xLWB20wKFphtYcBnj6dA4lxNxgsD5wGNJr8gkZMFljKej41xnAmvcsLs2ueWYllhwGeMpA8Y2DbS3hw3K+8CCyxjADbubgG3A0LY+x4k4xXgtrjldVZdpngWXMV9o7zjXecBCN+xu7aJ6TAIWXMZ8ob3jXLauvE8suIz5QptbXE7EORAYCczv0opMsyy4jPnCKqCfE3EOaMO5FwDz3LAb7eKaTDOy/S7AmB7knIKsguw6rdtw0oqT1lTGK3/ght1/7HywaYnnM4HhwKXA9/0qNNOJqk32NcaJOGcXB4rnzBg0I3Rs4bEsqVzCjHUzotWN1We5YfdVJ+IUhbJCr/cP9h88rnhc6I2KN7I21m38INYYG+OG3c/9rj/TWIvLGKAoUHTrzQNvDp1YciIAX+n9Faobq0P3bLznJ8ApuZJ707jicUfcMfiOXBFBVbl1/a2D/r7t77cDl/tZeyayMS5jgAZtGHx0wdF7fO2o0FEoOgQgJyvngmkHTMsVEQBEhAsPuDAoyJTur9ZYcBkDZEv2ivKK8j2+Vl5R3oiyFECQWHW8eo/Hq+PVCGLLNPvAgssYoDJeOf2XG38ZffTfjzaurF7J7z/9feMD6x5oWPX/VkWDweDd659Zv/zXG35dUxmvBCAWjzFz48xYndY96HPpGckG540BRCRnv0n7/fCQsw/5Ud6AvKLgZ8GGEZtGBAYEBmTV19fz1sq3ou+NfC8neHww59DcQ2s3NG4gS7LmVjdWX+SG3Xq/6880NjhvMp6IlBYVFc057NPDDr6p/02FE8+YSHZ2ds5ep4UAFq1exKy/zspa+fuV9dEN0VcbGhoasP/3dztrcZmMlpWV9Y1QKPSHmTNn5l966aWyc/C9NatXr+aCCy6oXrt27cLKysrJqlrbxaWa3VhwmYwVCASmFBcXP/Lqq6+GSktL2/382tpaJk+eHHvttdderaysnKiqtlFGN7HgMhlJRA7Nz89f9frrr4fC4XCHr1NXV8cJJ5xQvWLFip/V1NTcmcQSTQvsU0WTcUREiouLH7v55puDnQktgGAwyBNPPFEQCARuEZEjklSiaYUFl8lEYwsLC4dPnz69xQ+nVJW7776bww47jMLCQk455RQ++OCDfc4bPHgwP/rRj4JFRUX/02UVmz1YcJmMU1xc/IPrr78+Pzu75Q/Vr7zySmbOnMl9992H67oEAgG++c1vNnvu5ZdfHqivr58iIr27omazF1W1w46MOQDJzc2Nbtq0SVvy0ksvaVZWli5evHjX115++WUFNNFzJ0yYsAOY7Pd7zITDWlwm0xxWWFjYeMABLS+5dffdd3PuuecyevToXV/r06cPAJs3b272OSeeeGJhMBjs6N6Mph0suEymOfLII49saOmEWCzGggUL+NrXvrbH13cGVu/ezfcGhw0bllVUVNS50X7TJhZcJtPkFRYWtjjL9O2336ahoYG9P3FcunQpffr0YcCAAc0+LxQKQdMMe9O1LLhMpolWVFS0OHlx27ZtAPTr12+Pr8+dO5dJkyYlfF5VVRWqWp3wBJM0Flwm06xes2bN3vch7mH//fcHvggwgGeffRbXdbnmmmsSPu+tt96KV1RUvJmcMk1LLLhMpllfU1PTuHHjxoQnHH300QwaNIgZM2bw4YcfMmfOHKZOncrNN9/MMccck/B5CxcurG5oaFjSFUWbPVlwmYyiqhoMBuf98Y9/THhfYU5ODk899RTvv/8+paWl3HLLLfzqV79ixowZCa+7ceNGIpFIDvByV9Rt9mT3KpqMIyLhvn37vrZx48ZQbm5uUq75ox/9qOH+++//Y2Vl5XeTckHTImtxmYyjqpG6urrXb7nllqQsALhq1SpmzZpVW1VVdWsyrmdaZy0uk5FEpH9+fv57L7zwQtGJJ57Y4etUV1czatSo6vfff/8H9fX1toxzN7EWl8lIqropFotNnjhxYrSsrKxD16iqquL000+PbtiwYX5DQ8Nvk1yiaYEFl8lYqvrPqqqqyaeddlr1bbfd1tDQ0OKE+j2Ul5dz1FFHVb/11ltPVlVVXajWdelW1lU0GU9EBhUXFz/et2/f0htvvLFw6tSpFBQU7HOeqlJeXs4999wTfeGFF+qj0eilqvqUDyVnPAsuY/AWFwTOKCkp+WE0Gj1hyJAh0eOOOy7Yt2/fYG1tbeNbb70VW7ZsWU5jY+O2aDR6dzwe/z9V3e533ZnKgsuYvYhIMXAsMBwoAuqBjwAXWG/dQv9ZcBljUo4NzhtjUo4FlzEm5VhwGWNSjgWXMSblWHAZY1KOBZcxJuVYcBljUo4FlzEm5VhwGWNSjgWXMSblWHAZY1KOBZcxJuVYcBljUo4FlzEm5VhwGWNSjgWXMSblWHAZY1KOBZcxJuVYcBljUo4FlzEm5VhwGWNSjgWXMSblWHAZY1KOBZcxJuVYcBljUo4FlzEm5VhwGWNSjgWXMSblWHAZY1KOBZcxJuVYcBljUo4FlzEm5VhwGWNSjgWXMSblWHAZY1KOBZcxJuVYcBljUo4FlzEm5VhwGWNSjgWXMSblWHAZY1KOBZcxJuVYcBljUo4FlzEm5VhwGWNSjgWXMSblWHAZY1LO/wfCCJjaHitYXQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "icd = LearnStructICD(nodes_set=observed_set, ci_test=ci_test_dsep)  # instantiate ICD\n",
    "icd.learn_structure()  # learn the structure\n",
    "fig_pag = draw_graph(icd.graph)  # plot the PAG"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Finally, let's plot the number of CI tests per conditioning set size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of estimated CI tests\n",
      "----------------------------\n",
      "Condition set size: 0    1    2    3    4    5    \n",
      "Number of CI tests: 10   24   3    1    0    0    \n",
      "Total number of CI tests:  38\n"
     ]
    }
   ],
   "source": [
    "num_ci_tests = ci_test_dsep.test_counter\n",
    "str_cs_size = 'Condition set size: '\n",
    "str_num_cit = 'Number of CI tests: '\n",
    "max_str_len = len(str(max(num_ci_tests))) + 3\n",
    "for cs_size, num_cit in enumerate(num_ci_tests):\n",
    "    str_cs_size += str(cs_size).ljust(max_str_len)\n",
    "    str_num_cit += str(num_cit).ljust(max_str_len)\n",
    "print('Number of estimated CI tests')\n",
    "print('----------------------------')\n",
    "print(str_cs_size)\n",
    "print(str_num_cit)\n",
    "print('Total number of CI tests: ', sum(num_ci_tests))\n"
   ]
  }
 ],
 "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}