18 - Encoding expert knowledge.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "70efc90b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import networkx as nx\n",
    "\n",
    "from castle.algorithms import PC\n",
    "from castle.common.priori_knowledge import PrioriKnowledge\n",
    "\n",
    "from castle.datasets import IIDSimulation \n",
    "from castle.common import GraphDAG\n",
    "\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "624933fd",
   "metadata": {},
   "source": [
    "# Encoding expert knowledge in causal discovery\n",
    "\n",
    "A part of [CausalPython](https://causalpython.io) series on causality.\n",
    "\n",
    "<a href=\"https://causalpython.io\"><img src=\"img/CausalPython.io__flat.png\" width=150 align=\"left\"></a>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0ba361e",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2311a11b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a graph\n",
    "dag = np.array([\n",
    "    [0, 1, 1, 0, 0],\n",
    "    [0, 0, 0, 1, 0],\n",
    "    [0, 0, 0, 1, 0],\n",
    "    [0, 0, 0, 0, 1],\n",
    "    [0, 0, 0, 0, 0]\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "69f86e02",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAA46ElEQVR4nO3deVhU9f4H8PfMMMMMIpuoqAgIKBiIIiqiXs201Gvr7aY3csnuTdJKy9KyMnG91xZvVi7YbbNEq9uvRTNNc0tREVwAFUkWERUVBRGZYWbOnN8fpDfCZQZm5szyfj1PjwVz5nzsAd58vue7yERRFEFEROQm5FIXQEREZE8MPiIicisMPiIicisMPiIicisMPiIicisMPiIicisMPiIicisMPiIicisMPiIicisMPiIicisMPiIicisMPiIicisMPiIiciseUhdAROSqtHoB+eXVqKo1QC+YoFLI4d9CheigllArFVKX57YYfEREViKKIvYVX0J6ZimyT1biXLUOGqUCkAGiCMhkAERAaxDQ1keNhFB/JPcJQWKnAMhkMqnLdxsynsdHRNQ8dUYBazJLsXJnEaq0Bmj1Asz5wSoDoFEp4KdRYuLAcDzaJwSeHuwEbY3BR0TUDDllVXg6/QAqavTQGoQmv49GqUCgtwpLk3siLtjPegVSIww+IqImEEwiFm3Mx6q9JdAZTFZ7X7VSjnF9w/DS8Ggo5Bz+tAUGHxGRhfRGEyatzkZG4cVmdXk3o1Eq0D+yFZYlJ0Dlwcn31sb/o0REFhBMIiatzsbuwgqbhB5QP/ll14kKTE7PhmBib2JtDD4iIgss2piPjMKLVh3evBGdwYTdJy7ijU35Nr2PO+JQJxGRmXLKqjBq5R6zQk+lkCNlUDge7N4Bwf4aVOuM2FNUgcWbC1Bysdbse6qVcnw5MYkTXqyIwUdEZIY6o4Ahi3egrFJ729eqFHJ8/vdE9OkUgHPVOuwvuYRgfw16dPRHTZ0Rf1u5B3lnqs2+d0d/DbZMG8SlDlbCoU4iIjOsySzFxRq9Wa+dfGcE+nQKQPbJSgx+ezueWXMQDy7LwLwfjsLb0wNLRsfDkgmbFTV6fLH/VBMrpz9i8BER3YYoili5s8isySwechkm9O8EAHj9uzzU6v93zYe7inHsbDUi2njj7q5tzb6/1iAgbWcROEBnHQw+IqLb2Fd8CVVag1mv7R0WAF+NEiUXr+LI2cbDmRvyzgIAhloQfABQWatHZskli66hG2PwERHdRnpmKbR685YudG3XEgCQd/ryDT+fd7r6t9f5WFSDVi8gfV+pRdfQjTH4iIhuI/tkpVl7bwJAez8NAKC8WnfDz5dfrp8c08FfY1ENIoDs0kqLrqEbY/AREd2CVi/g3E1C7Ea8VB7Xr7uR2t+eE7ZQWX44ztnLOuhstGjenTD4iIhuIb+8uv5oITNdm6x5sw7xf5+3fKKKRqlAfvkVi6+jhhh8RES3UFVr+F9ameGq3ggA8FLdOCyvhWitmc8Mf08mq5/kQs3D4CMiugW9YIIlqwjOVNU/wwvyUd/w80G+9c/2TpuxEP6PRLF+g2xqHgYfEdEtqBRyWHI4+rGz9UORsR18b/j52A4+v73O/J1brpHJwNMarID/B4mIbsHPS3nzB3Y3sL/kEi5rDQhr1QIxN1iy8OfYdgCALfnnLa5FFAF/L5XF11FDDD4ioluIDvKx6Pgho0nEJxnFAIA598c0mBjz9wGd0LWdDwov1GDLsXMW16I1CIgOamnxddSQ5fNpiYjciEalQFsfNU5Xmf9Mbum2QvSLCETvsABsf/FO7C+5hA5+GsSH+ONqnRFT1x5s0jl77XzVUFsww5RujB0fEdFtJIT6WzKxE3rBhMf+sw+LNx9HTZ0Rd3dti44BXlh3+AzufX+XRSczXCMDkBDib/F11BiPJSIiuo29RRfxxKf7m7QEwVq8VAp8/HhvJHZqJVkNroIdHxHRbSR2CoCfRilpDf5eKvQJC5C0BlfB4CMiMkNv78vw9LBkwNN6NEoFUgaGQ2bJugq6KQYfEdEt5Ofn46677sIvn/wLARpp5gMGeqswundHSe7tihh8REQ3oNPpMHv2bAwYMAAPPvggMvdmIG1cH6iV9v2xqVbKseyxBHh6cDantTD4iIj+YOvWrYiLi0Nubi4OHTqEqVOnwsPDA3HBfhjXN8yiTaubQ6NUYHxSGLrdZBcYahrO6iQi+s2FCxfw4osvYtu2bXjvvffwwAMPNHqNYBKR8nkWdp2ogM5gu30z1Uo5BkQGIm1MLyjkfLZnTez4iMjtiaKIjz/+GLGxsWjVqhWOHDlyw9ADAIVchmXJCRgQGWizzk+jVGBAZCCWJScw9GyAHR8RubVjx47hqaeeQm1tLdLS0tCzZ0+zrhNMIt7YlI9P95RYtfNTK+UYnxSGGcOiGXo2wuAjIrek0+mwcOFCLFu2DLNnz8bkyZOhUFjeweWUVeHp9AOoqNFbtKfnH2mUcgR6e2Jpck/EBfs1+X3o9hh8ROR2fv75Z0yaNAndunXDkiVLEBwc3Kz3qzMKWJt5Cit/KUJlrR5avWDWgQ4y1O8Fqq26gL91D0TquHs4e9MOGHxE5DYuXLiAF154ATt27MD777+P++67z6rvL4oiMksuYfW+UhworUT5ZV39c0BZ/ZFCst/+1BkEBPmqkRDij+TEEPy6eyPS0lZg586dXKRuBww+InJ5JpMJH3/8MV555RWMHTsWqamp8Pb2tvl9dQYB+eVXUFmrh95ogspDDn8vFaKDWjY4ZcFoNCI6OhoffvghBg0aZPO63B2Dj4hc2tGjR/HUU09Bp9MhLS0N8fHxUpd0Qx999BHS09OxZcsWqUtxeVzOQEQuSavVYtasWRg0aBBGjRqFPXv2OGzoAcCYMWPw66+/Ys+ePVKX4vIYfETkcrZs2YK4uDjk5+fj0KFDeOaZZ5o0Y9OeVCoVXn75ZSxYsEDqUlwehzqJyGWcP38e06ZNw65du/D+++/j3nvvlboki+h0OkRERGDdunVmrycky7HjIyKnZzKZ8J///AexsbFo164djhw54nShBwBqtRrTp09n12dj7PiIyKkdPXoUKSkp0Ov1SEtLQ48ePaQuqVmuXr2K8PBwbN26FTExMVKX45LY8RGRU9JqtXjttdcwaNAgPProo8jIyHD60AOAFi1a4Pnnn2fXZ0Ps+IjI6WzevBmTJk1Cz5498c4776B9+/ZSl2RV1dXViIiIwO7du9GlSxepy3E5DD4ichrnzp3DtGnTsHv3bixduhQjR46UuiSbmTNnDkpKSvDxxx9LXYrL4VAnETk8k8mEDz74AN26dUOHDh1w5MgRlw49AJgyZQq+//57lJSUSF2Ky2HHR0QO7ciRI0hJSYHRaERaWhq6d+8udUl288orr6CyshLLly+XuhSXwuAjIoek1Woxb948fPDBB5g7dy4mTpzo8IvQre3ChQuIiopCbm4uOnToIHU5LoNDnUTkcDZt2oTY2FgUFhYiJycHkyZNcrvQA4DWrVtjwoQJePPNN6UuxaWw4yMih1FeXo5p06Zhz549WLZsGUaMGCF1SZI7e/YsYmJicOzYMbRt21bqclwCOz4ikpzJZMLKlSsRFxeHkJAQHDlyhKH3m3bt2iE5ORmLFy+WuhSXwY6PiCSVl5eHlJQUmEwmpKWlIS4uTuqSHE5paSni4+NRUFCAVq1aSV2O02PHR0SSqK2txcyZMzF48GCMHTsWu3fvZujdREhICP7yl79gyZIlUpfiEtjxEZHdbdy4EZMnT0ZiYiL+/e9/IygoSOqSHF5hYSESExNRWFgIX19fqctxagw+IrKb8vJyPP/889i3bx+WLVuG4cOHS12SUxk3bhyioqLw6quvSl2KU+NQJxHZnMlkwooVKxAXF4dOnTohLy+PodcEM2fOxJIlS1BTUyN1KU6NHR8R2VRubi5SUlIgk8mQlpaG2NhYqUtyaqNGjUJiYiJeeOEFqUtxWgw+IrKJ2tpazJ07Fx999BHmz5+Pf/zjH5DLOcjUXDk5ORg+fDgKCwuh0WikLscp8auQiKzuxx9/RGxsLEpLS5GTk4OJEycy9KwkLi4OvXv3xocffih1KU6LHR8RWc3Zs2fx3HPPISsrC8uXL8c999wjdUkuaf/+/fjLX/6CEydOwNPTU+pynA5/BSOiZjOZTFi+fDni4uIQGRmJvLw8hp4N9e7dGzExMVi1apXUpTgldnxE1Cw5OTlISUmBQqFAWloaYmJipC7JLezevRtjx45FQUEBPDw8pC7HqbDjI6ImuXr1KmbMmIGhQ4fiiSeewM6dOxl6dtS/f3+EhoYiPT1d6lKcDoOPiCy2YcMGxMbG4syZM8jNzcWTTz7JySsSeO2117Bw4UIIgiB1KU6FX6lEZLYzZ85g1KhRmDJlCtLS0vD555/zqBwJ3XXXXQgICMDXX38tdSlOhcFHRLclCAKWLVuG7t27o0uXLsjNzeXkFQcgk8nw2muvYf78+TCZTFKX4zQYfER0S4cPH0b//v2xZs0abN++HfPnz+fCaQcyYsQIKJVKrFu3TupSnAaDj4hu6OrVq5g+fTruvvtuPPnkk9ixYwcnrziga13fvHnzwEn65mHwEVEjP/zwA2JiYlBeXo68vDz8/e9/5+QVB/bAAw+grq4OmzZtkroUp8B1fER03ZkzZzB16lQcOnQIy5cvx9ChQ6Uuicy0du1avPfee9i1axdkMpnU5Tg0/gpHRBAEAUuXLkX37t0RHR2NnJwchp6TeeSRR1BRUYHt27dLXYrDY8dH5OYOHTqEiRMnQq1WIy0tDV27dpW6JGqiTz/9FKtWrcLPP/8sdSkOjR0fkZuqqanBiy++iGHDhuGpp57C9u3bGXpOLjk5GcXFxcjIyJC6FIfG4CNyQ+vXr0dsbCzOnz+P3NxcPPHEE5y84gKUSiVefvllzJs3T+pSHBqHOoncyOnTpzF16lQcPnwYK1aswJAhQ6Quiaysrq4OkZGR+Pbbb5GQkCB1OQ6Jv+IRuQFBEPDee++hR48eiImJQW5uLkPPRXl6emL69OmYP3++1KU4LHZ8RC7u4MGDmDhxIry8vLBixQo+x3MDWq0W4eHh+Omnn9CtWzepy3E47PiIXFRNTQ1eeOEFDB8+HJMnT+bkFTei0Wgwbdo0LFiwQOpSHBKDj8gFrVu3DjExMaioqEBeXh4mTJjARc1uZtKkSdi6dSuOHz8udSkOh0OdRC6krKwMU6ZMQV5eHlasWIG77rpL6pJIQvPnz8eJEyfwySefSF2KQ2HHR+QCBEHAu+++ix49eiAuLg45OTkMPcIzzzyD9evXo6ioSOpSHAo7PiInd+DAAaSkpKBFixZYsWIFoqOjpS6JHMisWbNw/vx5pKWlSV2Kw2DwETmpmpoavP7661i9ejUWLVqE8ePH8zkeNVJRUYEuXbrg8OHD6Nixo9TlOAQOdRI5oe+++w533HEHKisrceTIETz++OMMPbqhwMBA/P3vf8ebb74pdSkOgx0fkRM5deoUpkyZgqNHj2LFihUYPHiw1CWREygvL8cdd9yBo0ePIigoSOpyJMeOj8gJCIKAJUuWID4+Hj169EBOTg5Dj8wWFBSEMWPG4O2335a6FIfAjo/IwWVnZyMlJQUtW7bEihUrEBUVJXVJ5ITKysoQFxeHgoICBAYGSl2OpNjxETmoK1eu4LnnnsPIkSMxZcoUbN26laFHTRYcHIxHHnkE77zzjtSlSI7BR+SAvv32W8TExKC6uhp5eXkYN24cJ69Qs7300ktYsWIFqqqqpC5FUhzqJHIgp06dwrPPPov8/HysWLECd955p9QlkYt5/PHHERERgVmzZkldimTY8RE5AKPRiHfeeQfx8fHo2bMnDh8+zNAjm5g5cybeffddXLlyRepSJMOOj0hiWVlZSElJgZ+fH5YvX44uXbpIXRK5uL/97W9ISEjA9OnTpS5FEgw+IolUV1dj1qxZ+OKLL/Dmm29izJgxfI5HdpGbm4t77rkHhYWF8PLykrocu+NQJ7k1rV7AwdJKbMs/j01HyrEt/zwOnaqCziDY7J6iKOKbb75BTEwMrl69iiNHjmDs2LEMPbKbbt26oW/fvvjPf/4jdSmSYMdHbkUURewrvoT0zFJkn6zEuWodNEoFIANEEZDJAIiA1iCgrY8aCaH+SO4TgsROARYFU0ZGBjp06IDQ0NAGHy8tLcWzzz6LgoICpKWlYeDAgVb+GxKZJzs7Gw888AAKCwvh6ekpdTl2xeAjt1BnFLAmsxQrdxahSmuAVi/AnC98GQCNSgE/jRITB4bj0T4h8PRQ3PKaiooKhIWFoXPnzsjOzoZcLofRaMS7776LhQsX4rnnnsP06dPd7ocNOZ4///nPeOCBB5CSkiJ1KXbF4COXl1NWhafTD6CiRg9tM4YwNUoFAr1VWJrcE3HBfjd9XXJyMr7++msolUosXrwY8fHxmDhxIlq1aoXly5ejc+fOTa6ByJoyMjLw2GOPoaCgAEqlUupy7IbBRy5LMIlYtDEfq/aWQGcwWe191Uo5xvUNw0vDo6GQNxz+3L17N+6++25otVoAgFKphJ+fH/79738jOTmZz/HI4QwZMgRjx47F448/LnUpdsPgI5ekN5owaXU2MgovNqvLuxmNUoH+ka2wLDkBKo/6OWIGgwFRUVEoLi6+/jq5XI4RI0Zg/fr1Vq+ByBq2bduGlJQUHDt2DArFrYfxXQVndZLLEUwiJq3Oxu7CCpuEHlA/+WXXiQpMTs+GYKr/3XHx4sUoLy9v8DqTyYRt27Zhy5YtNqmDqLnuvPNOtG7dGl9++aXUpdgNOz5yOQs3HMNne0/aLPR+T6NUYFxSKB7tqkHnzp1hNBobfF4ul0OhUCAkJAQnTpyweT1ETbFx40a8+OKLyMnJgVzu+v0Qg49cSk5ZFUat3GPVZ3q3o1bKEfrr/2HzFx8iICAAQUFBaN++PUJDQxEaGop27dqha9eu6Nevn91qIrKEKIro3bs3Xn31VTz00ENSl2NzDD5yGXVGAUMW70BZpfa2r1Ur5RjYuTWGRLdFXLAvgv01UMhlOHmxFj/mncV/dhWjVm9+x9jRX4PNzw+EWunRnL8CkWS+/fZbzJ07F9nZ2S4/Ccv1e1pyG2syS3GxRm/Wax/o3gErx/bC6N4d4SGXYUfBBewvuYRgfw2m3R2F758egFYtVGbfu6JGjy+zyppaOpHk7r//fhiNRvz4449Sl2Jz7PjIJYiiiP6LtuLMZZ1Zr3+4ZwckhPjjw90lKLxQc/3jrVt64uPxvRHbwRffHTqNqV8cMruGDn4a7Jox2OV/WybX9cUXX+Cdd95BRkaGS38ds+Mjl7Cv+BKqtAazX//1gdN45du8BqEHABeu1OH17/MAAMNjgqBUmP/NX1mrR2bJJbNfT+Ro/vrXv6KyshJbt24FAGi1Wnz++ecud4QRg49cQnpmKbQWPJO7laNnqwEAnkoF/LzMH+7U6gWk7yu1Sg1EUlAoFHjllVeQmpqKN954A+3atcP48eNx6NAhqUuzKj6JJ5eQfbLSrL03zRESUH9Mi95owuVa87tIEUB2aaWVqiCyv8uXL6OwsBC7d+9GVlYWdDodfH19XW5hO4OPnJ5WL+BctXnP9swxoV8nAMCOggvQC5Ytizh7WQedQYBa6Vo/KMj1vfXWW0hNTYXJZIIoitDp6r+nRFF0ueDjUCc5vfzy6vqjhazgzqjWGN2rI/RGExZvPm7x9RqlAvnlrvU8hNxDdXU1BEG4vs/s7zH4iBxMVa2h/vygZopo3QLvjOoBuVyGhT8ew7EmBJhMVj/JhcjZzJ07F2vXroW3t3ej3VsYfEQORi+Y0NxFOW19PPHphD7w81JhxY5CfJJR0qT3EcX6Z4NEzuiBBx7A4cOH0blzZ6jVagCAIAguF3x8xkdOT6WQozlLjnw1Snz2RCKC/b2wJrMU/9qY3+T3kslw/bQGImcUHh6OQ4cOYeLEiVi7di20Wu0Ng0+rF5BfXo2qWgP0ggkqhRz+LVSIDmrp8M+4GXzk9Py8lGjqlE4vlQKfTuiNLm1b4ofcs3j129xm1SKKgL8FSyCIHJFarcaqVauQlJSEyZMnX5/wsq/4EtIzS5F9shLnqnX1z9Zl9V/3MhkAsf7kkrY+aiSE+iO5TwgSOwU43GJ47txCTk+rF9BtziYYTZZ9KasUcnz8eG/0jwzEjoLz+MeqLBiE5n07KOQyHEkd5vC/8RKZq+JSFdbnV2LlziJUaQ3Q6gWzfs+UAdCoFPDTKDFxYDge7RMCTw/H+L5gx0dOT6NSoK2PGqerbr859TVyGfDu3+LRPzIQ+4ovIuXz7GaHHgC081Uz9Mhl5JRV4en0A6io0Vt8zJcIoFYvoFYvYNHG4/hwVzGWJvdEXLCfTWq1BIOPXEJCqD/OVGnNHvEcnxSG4bFBAIDKqwbMfyD2hq9bvqMQhReumvWeMgAJIf5mVkDkuASTiEUb87Fqb4lVjvjSGgScqtRi1Mo9GNc3DC8Nj4ZCLt3wJ4OPXEJynxBsOXbO7KOEfDXK6/9+LQBv5L/ZZWYHn0alQHJiiFmvJXJUeqMJk1ZnI6PwotXPtdQZTPhs70kUVdRgWXKCZBPB+IyPXIKlpzPYAk9nIGcnmERM/CwLuwsrbHqYs1opx4DIQKSN6SVJ58d51+QSZDIZJg4Mt9oOLpbSKBVIGRjO0COntmhjvk06vT/SGUzYfeIi3tjU9KVDzcHgI5fxaJ8QBHpLs5Qg0FuF0b07SnJvImvIKavCqr0lFk9iaSqtQcCne0qQU1Zll/v9Hoc6yaXklFVh1Mo9Nv+N9ffUSjm+SumHbh187XZPImuqMwoYsngHyipvPzM6tr0P/tS5NboH+6FHRz8E+dbv8BI284cm3bujvwZbpg2y61IHdnzkUuKC/TCub5jdhjw1SgXGJ4Ux9MiprcksxcUa8/aYnXJXZ7w0PBrDY4Ouh15zVNTo8cX+U81+H0sw+MjlvDQ8Gv0jW0GttO2Xt1opR//IVpgxLNqm9yGyJVEUsXJnkdlDnAdKK7Hk5wI88el+9FqwGUYLj+76I61BQNrOIthz8JHLGcjlKOQyLEtOwOT0bOw+cdEmzyw0SgX6R7bCsuQESdcjETXXvuJLqNKaf+Dyip1FVq+hslaPzJJLSOzUyurvfSPs+MglqTzkSBvTC+OSQq3e+amVcoxLCkXamF7ckJqcXnpmKbRmrn+1Fa1eQPq+Urvdjx0fuSyFXIaZI7piZLd2Td526fc0SgUCvVUOs+0SkTVkn6xs6h7vViMCyC6ttNv9GHzk8uKC/bBl2iCszTyFlb8UobJWb/ZGu4AIL5UH/L1USBkYjtG9OzrMRrtEzaXVCzhXLd2mD7939rIOOoNgl71uGXzkFjw9FBjfLwzjkkKRWXIJq/eV4kBpJcovNz5aRRQBnUFA6xZKlB7Yjv/8ayqSIlpzcTq5nPzyamiUClypM0pdCjRKBfLLr6BHRz+b34vBR25FJpMhsVOr6w/RdQYB+eVXUFmrh95ogspDDn+v/x2m2bPnSxDODocsso3ElRNZX1WtoX53dQcgk9VPcrEHBh+5NbVSccvfMO+9916sW7cOf/rTn+xXFJGd6AUTHGULE1Gs3yDbHjgljegW7rvvPqxfv17qMohsQqWQw1FG8GUy2G2WNIOP6BYSEhJQWVmJwsJCqUshsjo/LyUkn9L5G1EE/L3ss9cug4/oFuRyOUaOHMmuj1xSdJCP3Talvh2tQUB0UEu73IvBR3Qb9957L4OPXJJGpUBbn+bvt2kN7XzVdlnKAHByC9FtDR06FOPGjUN1dTV8fHykLofIqhJC/XGmSmv2iOfgqDaYclfk9f/2UNT3T99M6nf9Y2v3n8IXWeZvPC0DkBDib/brm4sdH9FteHt7o3///vjpp5+kLoXI6pL7hECjMr/TatVChfgQ/+v/XPP7j7Wz8NQGjUqB5MQQi65pDnZ8RGa47777sG7dOvz1r3+VuhQiq1JVnYRaJqDWzNf/90AZ/nugzKo1+Hup0CcswKrveSvs+IjMMHLkSGzYsAGC4BgTAYia6/Dhw3jooYdw3333YWAbvd3OsPwjjVKBlIHhdt0ZicFHZIbQ0FC0b98e+/btk7oUombJycnBww8/jBEjRmDQoEE4ceIEFk28H4He9llK8EeB3iqM7t3Rrvdk8BGZibM7yZnl5eXhkUcewbBhw9C/f3+cOHECzz33HDQaDTw9FFia3NPmhzf/kVopx7LHEuy+8TuDj8hM17YvI3ImR48exejRozF06FAkJibixIkTmDZtGry8vBq8Li7YD+P6htltyFOjVGB8Uhi6dfC1y/1+j8FHZKY+ffrg3LlzKCkpkboUots6duwYHn30UQwePBgJCQkoLCzEiy++iBYtWtz0mpeGR6N/ZCubd35qpRz9I1thxrBom97nZhh8RGZSKBT485//zOFOcmjHjx/HY489hkGDBqF79+4oLCzEjBkzbhl41yjkMixLTsCAyECbdX4apQIDIgOxLDkBCrk0G4Uy+IgswE2ryVEVFBRg7Nix+NOf/oSYmBgUFhbi5Zdfhre3t0Xvo/KQI21ML4xLCrV656dWyjEuKRRpY3rZbUPqG5GJoqMcSkHk+KqrqxEcHIzTp0+jZUv77CtIdCsnTpzAvHnzsGHDBkydOhVTpkyx2g5DOWVVeDr9ACpq9M3a01OjVCDQW4WlyT0RF+xnldqagx0fkQV8fHzQt29fbNmyRepSyM0VFhZiwoQJSEpKQkREBE6cOIHXXnvNqtvqxQX7Ycu0QXh5eDQ6+GngpVKYfW6tDICXSoEOfhrMHBGNLdMGOUToAez4iCz27rvv4tChQ/joo4+kLoXcUHFxMebPn4/vvvsOzzzzDJ577jn4+fnZ/L6iKCKz5BJW7yvFgdJKlF/W1T8HlNUfKST77U+dQUCQrxoJIf5ITgxBn7AAuy5ONweDj8hCRUVFSEpKwtmzZyGXc9CE7KOkpAQLFizAN998g8mTJ+P555+Hv7/9Nnb+I51BQH75FVTW6qE3mqDykMPfS4XooJZ2O2WhqbhXJ5GFwsPDERgYiP379yMxMVHqcsjFlZaWYsGCBfjvf/+LSZMmoaCgAAEB9tvX8mbUSgV6dPSTuowm4a+rRE3A2Z1ka6dOncKkSZMQHx+PVq1aoaCgAPPnz3eI0HN2DD6iJuAuLmQrZWVlePrpp9GjRw/4+vri+PHjWLhwIVq1aiV1aS6DwUfUBElJSSgrK8OpU+Yftkl0K2fOnMGzzz6L7t27o0WLFsjPz8e//vUvBAYGSl2ay2HwETWBQqHAiBEj8MMPP0hdCjm5s2fPYurUqYiNjYWnpyeOHTuGN954A61bt5a6NJfF4CNqomuH0xI1RXl5OZ5//nnExMRAoVDg6NGjeOutt9CmTRupS3N5DD6iJho2bBh++eUXXL16VepSyImcO3cOL7zwAu644w6YTCYcOXIEixcvRlBQkNSluQ0GH1ET+fr6onfv3vj555+lLoWcwPnz5zF9+nR07doVer0eubm5WLJkCdq1ayd1aW6HwUfUDDyclm6noqICL730EqKjo1FbW4ucnBy899576NChg9SluS0GH1EzXFvPZzKZpC6FHMzFixcxc+ZMREVFobq6GocPH8bSpUsRHBwsdWluj8FH1AyRkZHw8fHBwYMHpS6FHMSlS5fw6quvokuXLrh06RIOHDiA5cuXo2PHjlKXRr9h8BE1ExezEwBUVlZi1qxZ6Ny5M86fP4/s7GykpaUhNDRU6tLoDxh8RM3E7cvcW1VVFWbPno3OnTvjzJkzyMrKwgcffICwsDCpS6ObYPARNVO/fv1QVFSEM2fOSF0K2dHly5cxZ84cREZGorS0FPv27cOHH36ITp06SV0a3QaDj6iZlEolhg8fzl1c3ER1dTXmzZuHyMhIFBUVYe/evfj4448REREhdWlkJgYfkRXwOZ/ru3LlChYsWICIiAgUFBRg9+7d+PTTTxEZGSl1aWQhBh+RFQwfPhzbt2+HVquVuhSysitXruCf//wnIiIicPToUezatQufffYZunTpInVp1EQMPiIrCAgIQHx8PLZu3Sp1KWQlNTU1WLRoESIjI5GTk4MdO3Zg9erViIqKkro0aiYGH5GVcHana7h69SrefPNNRERE4MCBA9i6dSvWrFmDrl27Sl0aWQmDj8hKrm1fJoqi1KVQE9TW1uLtt99GREQEMjMz8fPPP+OLL75ATEyM1KWRlTH4iKwkKioKarUahw8flroUsoBWq8W///1vREREICMjAz/99BO++uorxMbGSl0a2QiDj8hKZDIZZ3c6Ea1WiyVLliAiIgI7d+7Ejz/+iK+//hpxcXFSl0Y2xuAjsiI+53N8Op0O7733HiIjI7F161b88MMP+Oabb9CjRw+pSyM7YfARWdGAAQNQUFCAc+fOSV0K/UFdXR2WLl2KyMhI/PTTT/j+++/x3XffIT4+XurSyM4YfERWpFKpcM8993AXFwdSV1eH5cuXo3PnztiwYQO++eYbrFu3DgkJCVKXRhJh8BFZGQ+ndQx6vR5paWno0qULvv/+e3z11Vf44Ycf0Lt3b6lLI4nJRM69JrKqiooKRERE4Ny5c1Cr1VKX43YMBgM++eQTLFiwAFFRUUhNTUVSUpLUZZEDYcdHZGWBgYHo1q0bduzYIXUpbsVgMODDDz9Ely5d8OWXX2L16tXYtGkTQ48aYfAR2cB9993HZQ12YjQa8fHHHyM6Ohrp6elYtWoVNm/ejP79+0tdGjkoDnUS2cCRI0cwcuRIFBcXQyaTSV2OSzIajVi9ejXmzZuHjh07IjU1FYMGDZK6LHICHlIXQOSK7rjjDshkMuTl5aFbt25Sl+NSBEFAeno65s2bh6CgIHzwwQcYPHiw1GWRE2HwEdmATCa7vpidwWcdgiBg7dq1mDt3Llq3bo3ly5fjrrvuYkdNFuMzPiIb4fZl1iEIAtasWYPY2FgsXboU77//Pn755RcMGTKEoUdNwmd8RDZSV1eHNm3a4MSJE2jdurXU5Tgdk8mEr776CnPnzkXLli0xZ84c3HPPPQw7ajZ2fEQ24unpiaFDh2LDhg1Sl+JUrgVeXFwc3n77bbz11lvYs2cPhg0bxtAjq2DwEdkQN602n8lkwtdff40ePXrgjTfewKJFi7Bv3z6MGDGCgUdWxaFOIhs6f/48unTpgvPnz0OlUkldjkMSRRHffvst5syZA4VCgdTUVNx7770MO7IZzuoksqE2bdqga9eu2LlzJ4YOHSp1OQ5FFEV8//33SE1NBQDMmTMH999/PwOPbI7BR2Rj12Z3MvjqiaKI9evXIzU1FUajEampqXjwwQcZeGQ3HOoksrGcnBw8+OCDKCwsdOsf7qIoYsOGDUhNTYVOp0NqaioeeughyOWcakD2xY6PyMa6desGo9GI/Px8dO3aVepy7E4URWzcuBGpqam4evUqZs+ejYcffpiBR5Jh8BHZ2LVdXNatW+dWwSeKIn766Sekpqbi8uXLmD17Nh555BEGHkmOX4FEduBOh9OKoojNmzdjwIABeO655zBlyhTk5uZi9OjRDD1yCHzGR2QHOp0Obdu2RVFREdRqNX7++Wd4enpi2LBhUpdmNaIoYtu2bZg9ezbOnz+P119/HX/729+gUCikLo2oAQ51EtnB+fPnERoain79+uHkyZMwmUwYNmyYywTf9u3bMXv2bJw5cwavv/46Hn30UXh48McLOSZ+ZRLZ0AcffIB//vOfOHv2LEwmE/R6PQDAw8MDiYmJElfXfDt37sTs2bNx6tQpzJo1C4899hgDjxwev0KJbKi2thanT5++HnjXeHl5IS4uTqKqmm/Xrl2YPXs2iouLMWvWLIwdO5aBR06DT5qJbGjq1Kl4/vnn4eXl1eDjRqPRKc/py8jIwN13342xY8ciOTkZx48fx4QJExh65FQ4uYXIxkRRxFNPPYXPP/8ctbW1AOpPbqitrXWaWY579+7F7Nmzcfz4cbz66qsYP3489x4lp+Uc33VETkwmk2H58uW49957oVarAQDh4eFOEXqZmZkYMWIERo0ahb/85S8oKCjAk08+ydAjp+b433lELkAul2P16tUYNGgQACAyMlLiim4tKysLI0eOxMMPP4z7778fv/76K1JSUhh45BI4ME9kJx4eHvjuu+8QFRWFiIgIAIBWLyC/vBpVtQboBRNUCjn8W6gQHdQSaqX9179lZ2cjNTUVBw8exMyZM/F///d/8PT0tHsdRLbEZ3xEdiSKIvYVX0J6ZimyT1biXLUOGqUCkAGiCMhkAERAaxDQ1keNhFB/JPcJQWKnAJtucH3w4EGkpqYiKysLL7/8Mp588snrw7JErobBR2QHdUYBazJLsXJnEaq0Bmj1Asz5xpMB0KgU8NMoMXFgOB7tEwJPD+t1gocPH0Zqair27duHl156CRMnToRGo7Ha+xM5IgYfkY3llFXh6fQDqKjRQ2sQmvw+GqUCgd4qLE3uibhgv9u+PiMjA5s3b8bs2bMbfS43NxepqanIyMjAjBkz8NRTTzHwyG0w+IhsRDCJWLQxH6v2lkBnMFntfdVKOcb1DcNLw6OhkN94+PPChQuIjo5GTU0NcnJyEBUVBQDIy8vDnDlz8Msvv2D69OmYNGlSozWGRK6OszqJbEBvNGHiZ1n4bO9Jq4YeAOgMJny29yRSPs+C3tj4vU0mEx5++GFcuXIFgiDg5ZdfxtGjRzF69GgMGTIEffr0QWFhIV544QWGHrkldnxEViaYREz8LAu7CyusHnq/p1bKMSAyEGljejXo/BYsWICFCxdeXyyvUCjg6+uLGTNm4Omnn4a3t7fNaiJyBgw+IitbuOEYPtt7slnP88ylUSowLikUM0fUH3CbkZGBoUOHQqvVXn+NTCbD8OHDsWHDBpvXQ+QMGHxEVpRTVoVRK/fYtNP7I7VSji8nJqGDRkDnzp1RWVnZ6DUajQb79+9HTEyM3eoiclQMPiIrqTMKGLJ4B8oqtbd/MYDJgyIQH+KHzm1bolULFVQecly4Uoe9RReRtrMIv56vMfveHf01qFg1FflH8qBSqaBSqRqs+9PpdPjrX/+K9PR0i/9eRK6GwUdkJZ9kFGPRxuNmD3EenzscJhHIL6/Gueo6AECXIG+EB3qjzigg5fNsbD9+waz30igVuCvgMuK9r1x/hieTya7/AwBRUVFOfRQSkbUw+IisQBRF9F+0FWcu68y+pneYPw6fugy90HBYdFxSKObeH4vz1TokLdoKwWTet2gHPw12zRhs0x1eiFwBlzMQWcG+4kuo0hosumZ/SWWj0AOAVXtOouTiVbTxUaNzG/NnYFbW6pFZcsmiGojcEYOPyArSM0uh1VtvFqdRqO/ybrRO72a0egHp+0qtVgORq2LwEVlB9slKs/beNMdD8R0QHtgCRRU1KLl41ezrRADZpY1ndBJRQzyWiKiZtHoB56rNf7b3R8/eFYnQAC9oVB6IbO2NqKCWKL+sw5Q1B2Hm473rzl7WQWcQJDnSiMhZMPiImim/vBoapQJX6oxNuv6uqDaID/G//t9llbWY9uVh5J2ptvi9NEoF8suvoEdHvybVQuQOGHxEzVRVa6g/P6iJHlqeAQDwUXugazsfPD+0C75MScKbm45j6fYTFr2XTFY/yYWIbo7P+IiaSS+YYI1FQdU6I/YVX8K4jzJx7Gw1Xri7C+KCfS16D1G0bEIMkTti8BE1k0ohhzWXzukFE9blnIFcLsOQ6LYWXSuTASoPflsT3Qq/Q4iayc9LCatN6fzNpav1w5WtWqgsuk4UAX8vy64hcjd8xkdkobq6OuTk5CA7OxvZ2dnYf+Awau6ZBWv+HpnYqRUA4OSlWouu0xoERAe1tFodRK6IwUd0C3V1dcjNzb0ecllZWcjPz0fnzp2RkJCAXr164R//+Aee31Zj0XZlf+ocCC+lApuPnWuwZEGpkGFs31Dc3709avVGrDt8xqJ62/mquZSB6DYYfES/0ev1jULu2LFjiIyMvB5yTzzxBLp37w6NRtPg2l7FB7Hu8BmzRzw7t/HG6/fG4MKVOuSduYzLWgMCWqgQ3bYl2viooTMImPblYZRbsD5QBiDhd8siiOjGGHzklvR6PY4cOYKsrKzrIXf06FFERERcD7nHH38c3bt3h5eX123fL7lPCLYcO4daM7ct2378Atq0LETvsADEtPOBn5cKesGEssparM89i08ySlBq4TCnRqVAcmKIRdcQuSOezkAuz2AwNAq5I0eOIDw8/HrIJSQkoHv37mjRokWT7tGU0xmsjaczEJmHHR+5FIPBgKNHjzYKubCwMCQkJCAhIQFjxoxBjx49mhxyNyKTyTBxYLhF5/FZk0apQMrAcIYekRnY8ZHTMhqNDUIuOzsbubm5CAkJud7F9erVCz169Lh+OKst1RkFDF28A6fMPIHdmjr6a7Bl2iB4enBiC9HtMPjIKRiNRhw7dux6F5ednY2cnBx07NixUci1bCnddP6csiqMWrkHOoP9dk9RK+X4KqUfunWwbJcXInfFoU5yOEajEfn5+Y1CrkOHDtdDbtSoUYiPj5c05G4kLtgP4/qG4bO9J+0y5KlRKjAuKZShR2QBdnwkKUEQGoXc4cOH0b59+wadXHx8PHx8fKQu1yyCSUTK51nYdaLCpp2fWinHgMhApI3pBYWcz/aIzMXgI7sRBAHHjx9vFHJBQUHXQy4hIQE9e/aEr69zdzB6owmT07Ox+8RFm3R+GqUC/SNbYVlyAvfmJLIQg49sQhAEFBQUNFgMfujQIbRt27ZRyPn5+Uldrk0IJhFvbMrHp3tKrNr5qZVyjE8Kw4xh0ez0iJqAwUfNZjKZbhhyrVu3brBOrmfPnvD3d7+dRXLKqvB0+gFU1Oib1f1plAoEequwNLkn4oL9rFcgkZth8JFFTCYTfv311wYhd/DgQQQGBjYKuYCAAKnLdRh1RgFrM09h5S9FqKzVQ6sXzNreTIb6HVn8vVRIGRiO0b07cskCUTMx+OimTCYTCgsLGywGP3jwIAICAhqFXKtWraQu1ymIoojMkktYva8UB0orUX5ZB41SAcjqjxSS/fanziAgyFeNhBB/JCeGoE9YABenE1kJg48A1P9AvlHI+fr6Nphd2bNnTwQGBkpdrsvQGQTkl19BZa0eeqMJKg85/L1UiA5qyVMWiGzEaYJPqxeQX16NqloD9IIJKoUc/i34A6IpRFFEUVFRgx1PsrOz4evr26CTS0hIYMgRkctx2OATRRH7ii8hPbMU2Scrca668ZAQxPqDN9v6qJEQ6o/kPiFI7MQhod8TRRHFxcWNQq5ly5bXw+1a0LVu3VrqcomIbM7hgq/OKGBNZilW7ixCldZg8SQAP40SEweG49E+IW43CUAURZSUlDRYJ5ednY0WLVo0Crk2bdpIXS4RkSQcKvg47dt8oiji5MmTjUJOrVY3eCaXkJCAtm3bSl0uEZHDcIjgE0wiFm3Mx6q91l/oO65vGF4a7twLfUVRRGlpaaOQU6lUjUIuKChI6nKJiBya5MGnN5owaXU2Mgq5tRNQH3KnTp1qFHIeHh4NJp0kJCSgffv2UpdLROR0JA0+wSRi4mdZ2F3onpv5iqKIsrKyBovBs7OzIZfLG3RxvXr1YsgREVmJpMG3cMMxux/fMnNEV5vf60ZEUcTp06cbhRyAG4YcZ6YSEdmGZMHX1AM7VQo5Nkz5EyLbeMMomBD52o9mX6tWyvHlxCS7THg5c+ZMg8Xg2dnZMJlMjdbJBQcHM+SIiOxIkuCrMwoYsngHyiq1Fl877e4ueObOSMjlMouDDwA6+muwZdqgGy51OHToEB5++GF89tln6Nevn9nvefbs2QbP47KysmA0GhuFXMeOHRlyREQSk+QE9jWZpbhYo7f4ui5tvfHUwAis3V+K5MTQJt27okaPL/afwriksAYf/+ijj/Dss89Cp9Nh27ZtNw2+8vLyRiGn1+uvh9uECRPw/vvvIyQkhCFHROSA7N7xiaKI/ou24sxlnUXXyWTA/z3VD8H+GgxZvAM5s4c1qeMDgA5+GuyaMRgymQx1dXVISUnBV199hdraWgDAPffcg02bNuHcuXMNhiqzsrKg0+kaLQYPDQ1lyBEROQm7d3z7ii+hSmuw+LrHk8IQH+KPZ9ccQLXO2KwaKmv1yCy5hHaKqxgxYgSKi4uh1f5v2HX79u0IDg5GbW3t9ZAbO3Ys3nnnHYSFhTHkiIicmN2DLz2zFFq9ZbM42/uq8cI9UdhRcB7rcs42uwatXsC8Nduwec5jMBgM+GPTazKZ8N1336Fnz54MOSIiF2P3Fd3ZJyvN2nvz9xY81A0KmQyvfZtnlRpEAAWXjJDJZFCpVPD09GzweS8vL1y8eJGhR0TkguwafFq9gHPVlj3be6BHewyOaoN3t/6KU02YBXozgqcPKqtrkJ2djYULF6J///7w9PSEj48PampqsH//fqvdi4iIHIddhzrzy6uhUSpwpc68Z3T+XkrMGnkH8sur8cEvRVatRaNU4Pi5GvSIiUFMTAymTZuGuro6ZGRkYP369ejaVZqF7kREZFt2Db6qWkP9+UFmmnXvHQjwUmHiZ1kwmqw7+VQmq5/k8nuenp4YPHgwBg8ebNV7ERGR47Br8OkFEyxZPDE0ui10RgEzhkXf8PMeCjnWPtkXADB3/VEcPVtt9nuLYv0G2URE5F7sGnwqhRyWzhfxUnmgb3irm37+2ud81Jb9VWQyOMVpDUREZF12DT4/LyUsmdIZN/enm36u5J8jm7yAHajv+Py9VE26loiInJddW57oIB+7nMRgDq1BQHRQS6nLICIiO7Nr8GlUCrT1UdvzljfVzlcNtbLxRtVEROTa7P6QKyHU35KJnTYhA5AQ4i9xFUREJAW7b1mW3CcEW46dQ62F25b9UdjMH5p8rUalQHJiSLPuT0REzsnuHV9ipwD4aZT2vm0D/l4q9AkLkLQGIiKSht2DTyaTYeLAcGgker6mUSqQMjCc+3ASEbkpSRayPdonBIHe0iwlCPRWYXTvjpLcm4iIpCdJ8Hl6KLA0uSfUSvveXq2UY9ljCfD04GxOIiJ3JdnWJXHBfhjXN8xuQ54apQLjk8LQrYOvXe5HRESOSdI9u14aHo3+ka1s3vmplXL0j2x10z0/iYjIfUgafAq5DMuSEzAgMtBmnZ9GqcCAyEAsS06AQs4JLURE7k4mipacl2AbgknEG5vy8emeEugM1jsxQa2UY3xSGGYMi2boERERAAcJvmtyyqrwdPoBVNTom7Wnp0apQKC3CkuTeyIu2M96BRIRkdNzqOADgDqjgLWZp7DylyJU1uqh1QtmHeggQ/2OLP5eKqQMDMfo3h05e5OIiBpxuOC7RhRFZJZcwup9pThQWonyy7r654Cy+iOFZL/9qTMICPJVIyHEH8mJIegTFsDF6UREdFMOG3x/pDMIyC+/gspaPfRGE1Qecvh7qRAd1JKnLBARkdmcJviIiIisQdLlDERERPbG4CMiIrfC4CMiIrfC4CMiIrfC4CMiIrfC4CMiIrfC4CMiIrfC4CMiIrfC4CMiIrfC4CMiIrfC4CMiIrfC4CMiIrfy/3XPV4NWF8CkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the causal graph\n",
    "plt.figure(figsize=(6, 4))\n",
    "nx.draw(\n",
    "    G=nx.DiGraph(dag),\n",
    "    node_size=1200,\n",
    "    arrowsize=17,\n",
    "    with_labels=True,\n",
    "    font_color='white',\n",
    "    font_size=21,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "c4f9a36e",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2023-03-31 10:03:41,720 - C:\\Users\\aleks\\anaconda3\\envs\\causal_book_py38\\lib\\site-packages\\castle\\datasets\\simulator.py[line:270] - INFO: Finished synthetic dataset\n"
     ]
    }
   ],
   "source": [
    "# Generate teh data\n",
    "data = IIDSimulation(\n",
    "    W=dag, \n",
    "    n=1000, \n",
    "    method='linear', \n",
    "    sem_type='gauss'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d78ddb3c",
   "metadata": {},
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "481a6c00",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Train the causal discovery model\n",
    "pc = PC(variant='stable')\n",
    "\n",
    "# Learn\n",
    "pc.learn(data.X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "0d7a3bc2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdoAAADRCAYAAABrVHPYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZpklEQVR4nO3dfbBcdX3H8feHKwgWNNRcLSSEoI0PjFNRroC1KsWnJIqpM04LKFaKjbSi2LFVtD7RllbrQ9ERTa/IoMWRUkANGIvMKFKVh+QKIjFgIwS4gIYIyJOCId/+cc6lm2Xv7t5z9uz53T2f18wOd++e3d93w/3s9zzt+SkiMDMzs2rsUncBZmZmo8yN1szMrEJutGZmZhVyozUzM6uQG62ZmVmF3GjNzMwq5EbbMJLOkvRPdddhZoMj6c2Svld3HdaZG20FJC2VFJIeV3ctZk0naYukl9ddhzWXG+08ooz/n5kNSCorw6nUYdXwh3YfJO0r6XxJd0q6SdI78t8fImmDpHsl/ULSJ/OnXJb/9x5J90t6YZfXHpP0CUnb8tc+sXVrWNKlkk6V9H3gQeBpko6TtEnSfZJulPTWltc7XNK0pPflr7lF0hvaht1b0jfy518p6ekD+8cyS4ik/wCWABfmWXx3nq/jJd0CfHsmM23Pe3QrWNIukk6W9DNJv5R0rqTf7WPsN0m6OX/OB9pe88OSzpN0tqR7gTfnnyeXS7pH0h2SPiNpt5bXC0nvyDO/TdLH2le8JX1c0t35Z8mK8v+CNghutD3kf8gXAj8CFgEvA94p6VXAp4BPRcQTgacD5+ZPe0n+3wURsWdEXN5liL8EVgAHAc8H/qTDMscCq4G9gJuBrcBrgCcCxwH/Jun5Lcv/HrAwr/fPgUlJz2x5/GjgFGBvYDNwatd/BLN5KiKOBW4BjoyIPfn/jL4UeDbwqj5e5h1kuXwpsC9wN3B6tydIOhD4LPAGYB/gSWR5bLUKOA9YAHwZeAT4G7LsvpDss+av257zOmCC7LNiFfAXLY8dCtyQP/9fgS9IUh/vzyrmRtvbC4DxiPiHiHg4Im4EPg8cBfwW+H1JCyPi/oi4osDr/ylZs56OiLuBj3RY5qyI2BgR2yPitxHxjYj4WWS+C3wLeHHbcz4QEQ/lj38jH2fGBRFxVURsJwv4QQXqNpvPPhwRD0TEr/tY9q3A3+cZfQj4MPD6Hrt7Xw9cGBHfi4iHgQ8C7ReWvzwivhYROyLi1xExFRFX5DnfAvw7WXNv9dGIuCsibgFOI1tpnnFzRHw+Ih4BvkjW4J/ax/uzirnR9rY/sG++O+ceSfcA7yP7Az4eeAZwvaT1kl5T4PX3BW5tuX9rh2V2+p2kFZKukHRXXs9KsrXYGXdHxAMt92/Ox5nx85afHwT2LFC32XzWKWez2R/4akv+N5FtfXZrYjvlOiIeBH7ZrQZJz5B0kaSf57uT/5mdc93+nFlznY8HznYS3Gh7uxW4KSIWtNz2ioiVEfG/EXE08BTgo8B5kn6Hx665dnMHsLjl/n4dlnn09SQ9Hjgf+Djw1IhYAKwDWncR7Z3XMWMJcPscajIbJZ3y2Pq7B4AnzNyRNAaMtzx+K7Ci7TNg94i4rcuYO+Va0h7Ak3vU9TngemBZfjjqfeyca9j588G5nifcaHu7CrhX0nsk7ZGfvPQcSS+Q9EZJ4xGxA7gnX/4R4E5gB/C0Pl7/XOAkSYskLQDe02P53YDH52Nsz094eGWH5U6RtJukF5Mdz/2vPmoxG0W/oHsWfwrsLunVknYF3k+WsRlrgFMl7Q8gaVzSqh5jngccKekP8xOaTuGxTbPdXsC9wP2SngX8VYdl/k7S3pL2A04C/rPHa1oC3Gh7yI93HEl2HPMmYBtwBtnJDcuBjZLuJzsx6qiI+E2+2+ZU4Pv57qbDugzxebJjrNcCV5NtnW4na9id6rmP7OSMc8lOyjgGWNu22M/zx24nOwZ7QkRcP7d3bjYy/gV4f77b9/XtD0bEr8hOOjoDuI1sC7f1LORPkWXsW5LuA64gO/FoVhGxEXg7cA7Z1u19ZCcxPtTlaX9Lluf7yD4XOjXRrwNTwDVk5158oVsdlgZ54ve05FuoayJi/4LPPxw4OyIW91jUzIZE0p5ke72WRcRNBV8j8udvHmRtVj1v0dYs3x29UtLjJC0CPgR8te66zKwcSUdKekJ+vsTHgR8DW+qtyurgRjsEktbkX5Zvv60hO25zCtmu3qvJzmj8YJ312txIOlPSVknXzfK4JH1a0mZJ17Z959nmKUlvmCXXG/NFVpEdvrkdWEZ2aMm7EBNXRZ6969isJEkvAe4HvhQRz+nw+Eqy43UryY7tfSoiuh7jM7N6VJFnb9GalRQRlwF3dVlkFVloI7+oyQJJ+wynOjObiyryXMmFrBcuXBhLly6t4qX7MjU1VdvYKTn44IPrLqFWW7ZsYdu2bT0vQZefZNLNRuA3LfcnI2JyDqUsYucLDUznv7tjDq9RC2c5Dc5yf1mGNPNcSaNdunQpGzZsqOKl++LLe2bq/H+QgomJib6X3WWX2Xfu7Nix4zcR0f+LPVanP8h5cczGWU6Dszy3+KWWZ0/NZEblH+jT7HxFn8X4ij5mlUktzz5Ga40nibGxsVlvA7AWeFN+tuJhwK8iIvndxmbzUYp59hatGeXWgCV9BTgcWKhsXtMPAbsCRMQasqt9rSSbkvBBsqkNzawiqeXZjdaM7sd0esknluj2eABvKzyAmc1Janl2o7XGm9nVZGbzX4p5dqM1w2e3mo2S1PLsRmuNJ6nUriYzS0eKeXajNSO9NWAzKy61PLvRmkFyx3TMrLjU8uxGa40nKbk1YDMrJsU8u9GaUe7rAGaWltTy3Fc1kpZLuiGff+/kqosyG7Zddtll1tsocZatCVLLc88tWkljwOnAK8iu8bhe0tqI+EnVxZkNQ4pnKVbBWbYmSDHP/VRzCLA5Im6MiIeBc8jm4zMbGamtAVfEWbZGSC3P/Yw629x7ZiNj5gSKTrcR4ixbI6SW535Ohupr7j1Jq4HVAEuWLClZltnwpLirqSLOso28FPPcTzV9zb0XEZMRMRERE+Pj44Oqz2woUtvVVBFn2RohtTz3s0W7Hlgm6QDgNuAo4JhKqzIbohTXgCviLNvISzHPPRttRGyXdCJwMTAGnBkRGyuvzGyIUgtmFZxla4rU8tzXBSsiYh3ZZLdmIym1YFbFWbYmSC3PvjKUNV6Ku5rMrJgU8+xGa0Z6a8BmVlxqeXajtcaTlNxsH2ZWTIp5dqM1I735K82suNTy7EZrjZfiGrCZFZNint1ozUhvomgzKy61PLvRmpHeriYzKy61PKd1apZZDWZ2Nc126/M1us7zKulJki6U9CNJGyUdN/A3Ymal81xFlt1ozSg320fLPK8rgAOBoyUd2LbY24CfRMRzgcOBT0jabbDvwsygeJ6ryrJ3HVvjDeDkiUfnec1fb2ae19YJ1QPYS1nS9wTuAraXGdTMHqtknivJshutGaW/4N5pntdD25b5DLCWbLacvYA/i4gdZQY1s85K5LmSLFfSaKempmo9GB3xmCk2G6nuEwLmy/+HPi7ZtlDShpb7kxEx2foSHZ7T/uZfBVwDHAE8HbhE0v9ExL0FSm6M+fI3VDVnuX8l81xJlr1Fa0bPrwNsi4iJLo/3M8/rccBHIvvE2izpJuBZwFUFyjWzLkrkuZIs+2Qoa7yZNeASE0U/Os9rflLEUWS7llrdArwsH++pwDOBGwf4NsyM0nmuJMveojWj3DHa2eZ5lXRC/vga4B+BsyT9mGz31HsiYlv5ys2sXdE8V5VlN1ozys/20Wme1zyUMz/fDryy1CBm1peSK84Dz7IbrTVeitdGNbNiUsyzG60Z9Z/VaWaDk1qe3Wit8VJcAzazYlLMsxutGemtAZtZcanl2Y3WGq+PL7ib2TyRYp7daM0of9axmaUjtTy70ZqR3q4mMysutTy70VrjpXjyhJkVk2Kee25fSzpT0lZJ1w2jILM6lJmPdj5xnq0JUstzPzuyzwKWV1yHWa1KXut4PjkL59lGXGp57rnrOCIuk7R0CLWY1WIUt1xn4zzbqEsxzwM7RitpNbB6UK9nNkwjuOVaWGuWlyxZUnM1ZnOXWp4HVk1ETEbERI95O82SlNqupjq1Znl8fLzucszmLLU8+6xja7wUdzWZWTEp5tmN1oz0djWZWXGp5bmfr/d8BbgceKakaUnHV1+W2XCl9nWAqjjP1gSp5bmfs46PHkYhZnVJ8dqoVXGebdSlmGfvOjYjvUu2mVlxqeXZjdaM9I7pmFlxqeXZjdYaL8VdTWZWTIp5dqM1I71dTWZWXGp5dqM1I71dTWZWXGp5dqM1I701YDMrLrU8u9Fa46V4TMfMikkxz2lVY1aTstdGlbRc0g2SNks6eZZlDpd0jaSNkr470DdgZo8qk+cqsuwtWmu8smvAksaA04FXANPAeklrI+InLcssAD4LLI+IWyQ9pVzVZtZJmTxXlWVv0ZpReov2EGBzRNwYEQ8D5wCr2pY5BrggIm4BiIitA30DZvaoEnmuJMtutGb0vDbqQkkbWm7t8y4vAm5tuT+d/67VM4C9JV0qaUrSm6p7N2bNViLPlWTZu46t8frY1bStxzzLnU5xjLb7jwMOBl4G7AFcLumKiPjpnIo1s65K5rmSLLvRmlH66wDTwH4t9xcDt3dYZltEPAA8IOky4LmAG63ZgJXIcyVZ9q5jM0ofo10PLJN0gKTdgKOAtW3LfB14saTHSXoCcCiwaaBvwsyAUnmuJMveorXGK3vWcURsl3QicDEwBpwZERslnZA/viYiNkn6b+BaYAdwRkRcN4DyzaxFmTxXlWU3WjPKX7ItItYB69p+t6bt/seAj5UayMx6KrniPPAsu9Gakd61Uc2suNTy7EZrjddy2r+ZzXMp5tmN1oz01oDNrLjU8uxGa0Z6s32YWXGp5dmN1hpPEmNjY3WXYWYDkGKe3WjNSG8N2MyKSy3PbrRmpHdMx8yKSy3PPauRtJ+k70jalM+9d9IwCjMblpkvuJeZj3Y+cJatCVLMcz9btNuBd0XEDyXtBUxJuqR1fj6z+S61XU0VcZatEVLLc89GGxF3AHfkP98naRPZtEEOp42MUdpynY2zbE2RWp7ndIxW0lLgecCVlVRjVoMUv+BeNWfZRlWKee670UraEzgfeGdE3Nvh8dVA+4TYZvNCal8HqNJcsrxkyZIhV2dWXmp57mv7WtKuZMH8ckRc0GmZiJiMiIkeE2SbJWlmLbjTbZTMNcvj4+PDLdBsAFLLc88tWmWVfQHYFBGfrL4ks+FK8QvuVXCWrQlSzHM/W7QvAo4FjpB0TX5bWXFdZkOV2tcBKuIsWyOklud+zjr+HjBa+8/MWoziLuJOnGVrghTz7CtDmZHeyRNmVlxqeXajtcZLcQ3YzIpJMc9utGak9wV3MysutTy70ZqRXjDNrLjU8uxGa403cxFyM5v/UsyzG60Z6a0Bm1lxqeXZjdYaL8U1YDMrJsU8u9Gakd4asJkVl1qe06rGrCZlr40qabmkGyRtlnRyl+VeIOkRSa8fWPFmtpMyea4iy96itcYre21USWPA6cArgGlgvaS17ROq58t9FLi4RLlm1kWZPFeVZW/RmlF6i/YQYHNE3BgRDwPnAKs6LPd2splztg6ucjNrVyLPlWS5ki3agw8+mA0bNlTx0jYHEVHr+KldnaWbHsd0Fkpq/YOejIjJlvuLgFtb7k8Dh7a+gKRFwOuAI4AXlKvWmsZZnpsSea4ky951bI3Xx66mbT3mWe70KdT+yXga8J6IeGS+fWiZzScl81xJlt1ozSi9xj4N7NdyfzFwe9syE8A5+TgLgZWStkfE18oMbGaPVSLPlWTZjdYabwATRa8Hlkk6ALgNOAo4pnWBiDigZbyzgIvcZM0Gr2SeK8myG60Z5bZoI2K7pBPJzkAcA86MiI2STsgfXzOYKs2sH0XzXFWW3WjNKH+yR0SsA9a1/a5jKCPizaUGM7OuSq44DzzLbrRmzL+zKs1sdqnl2Y3WGi/Fa6OaWTEp5tmN1oz0ro1qZsWllmc3WjPS29VkZsWllmc3Wmu8uUweYGZpSzHPbrRmpLcGbGbFpZZnN1oz0jumY2bFpZbnno1W0u7AZcDj8+XPi4gPVV2Y2TCltgZcBWfZmiK1PPezRfsQcERE3C9pV+B7kr4ZEVdUXJvZUKT4dYCKOMs28lLMc89GG9n8TPfnd3fNb/XO2WQ2YKmtAVfBWbamSC3PfbV9SWOSriGb5PaSiLiywzKrJW2QtOHOO+8ccJlm1So58fu84SxbE6SW574abUQ8EhEHkU0ZdIik53RYZjIiJiJiYnx8fMBlmlVnZlfTbLdR4izbqEsxz3MaNSLuAS4FlldRjFldUlsDrpqzbKMstTz3bLSSxiUtyH/eA3g5cH3FdZkNVWprwFVwlq0pUstzP2cd7wN8UdIYWWM+NyIuqrYss+Ea1S3XNs6yNUJqee7nrONrgecNoRazWozyLuJWzrI1QYp59pWhzEjvSjJmVlxqeXajNSO9XU1mVlxqeXajtcZL8UoyZlZMinlOqxozM7MR4y1aM9I7pmNmxaWW57SqMTMzGzHeojUjvTVgMysutTynVY1ZDbpdrq3fsxclLZd0g6TNkk7u8PgbJF2b334g6bkDfyNmVjrPVWTZW7RmlPs6QH6lpdOBVwDTwHpJayPiJy2L3QS8NCLulrQCmAQOLVGymc2iaJ6ryrIbrRmlv3d3CLA5Im7MX+scYBXwaDgj4gcty19BNnuOmVWgRJ4rybIbrVUmm2e8PhMTE30vW7LRLgJubbk/Tfc13OOBb5YZ0GyY5lOWoVSeK8myG601Xh9fcF8oaUPL/cmImGx9iQ7P6fjJJOmPycL5R3Mu1Mx6KpnnSrLsRmvW27aI6LZKPQ3s13J/MXB7+0KS/gA4A1gREb8cbIlm1qduea4ky260ZpT+OsB6YJmkA4DbgKOAY1oXkLQEuAA4NiJ+WmYwM+uuRJ4rybIbrRnljtFGxHZJJwIXA2PAmRGxUdIJ+eNrgA8CTwY+m4+1vcdWspkVVDTPVWXZjdaM8rN9RMQ6YF3b79a0/PwW4C2lBjGzvpRccR54lt1orfFSnCjazIpJMc9utGakN3+lmRWXWp7daM1I79qoZlZcanlOqxozM7MR4y1aa7w+vuBuZvNEinlOqxozM7MR4y1aM9I7pmNmxaWW576rkTQm6WpJF1VZkFkdys5HO584yzbqUsvzXNr+ScCmqgoxs6Fxls2GqK9GK2kx8GqyiyibjZRua7+jtkXrLNuoSzHP/R6jPQ14N7BXdaWY1WfUGmoXp+Es24hLLc89t2glvQbYGhFTPZZbLWmDpA133nnnwAo0G4bU1oCr4CxbU6SW5352Hb8IeK2kLcA5wBGSzm5fKCImI2IiIibGx8cHXKaZDYCzbFaDno02It4bEYsjYinZ3Hzfjog3Vl6Z2RCltgZcBWfZmiK1PPt7tGakd0zHzIpLLc9zarQRcSlwaSWVmNVk1LZc++Es26hKMc9pXT7DzMxsxHjXsRnp7Woys+JSy7MbrRnpBdPMikstz951bGZmViFv0ZqR3hqwmRWXWp7daM1IL5hmVlxqeXajtcZL8esAZlZMinl2ozUjvTVgMysutTz7ZCizAZC0XNINkjZLOrnD45L06fzxayU9v446zay7KrLsRmtGuWujShoDTgdWAAcCR0s6sG2xFcCy/LYa+Nxg34GZzSia56qy7EZrRumLkB8CbI6IGyPiYbKZcVa1LbMK+FJkrgAWSNpnsO/CzKBUnivJciXHaKemprZJurnESywEtg2qnnk4fgo11D3+IGrYv5+FpqamLpa0sMsiu0va0HJ/MiImW+4vAm5tuT8NHNr2Gp2WWQTc0U+NdRmBLKdQQ93jp1DDULIMpfNcSZYrabQRUWoSS0kbImJiUPXMt/FTqKHu8YdZQ0QsL/kSnVaTo8AyyZnvWU6hhrrHT6GGYY5fMs+VZNm7js3Kmwb2a7m/GLi9wDJmVq9KsuxGa1beemCZpAMk7UY2qfratmXWAm/Kz1g8DPhVRCS929isgSrJcqrfo53svchIjw/111D3+JBGDT1FxHZJJwIXA2PAmRGxUdIJ+eNrgHXASmAz8CBwXF31DlkK/w/rrqHu8aH+Guoevy9VZVkRyR8mMjMzm7e869jMzKxCbrRmZmYVSqrR9rr01RDGP1PSVknXDXvsfPz9JH1H0iZJGyWdVEMNu0u6StKP8hpOGXYNeR1jkq6WdFEd41t5znO9eXaW05FMo+3z0ldVOwso+53KMrYD74qIZwOHAW+r4d/gIeCIiHgucBCwPD+zbthOAjbVMK4NgPMM1J9nZzkRyTRa+rv0VaUi4jLgrmGO2Tb+HRHxw/zn+8j+OBcNuYaIiPvzu7vmt6GeMSdpMfBq4IxhjmsD5TzXnGdnOR0pNdrZLmvVSJKWAs8Drqxh7DFJ1wBbgUsiYtg1nAa8G9gx5HFtcJznFnXl2VlOQ0qNdl5eoq4KkvYEzgfeGRH3Dnv8iHgkIg4iu+LJIZKeM6yxJb0G2BoRU8Ma0yrhPOfqzLOznIaUGq0vUQdI2pUslF+OiAvqrCUi7gEuZbjHuV4EvFbSFrLdjUdIOnuI49tgOM+kk2dnuV4pNdp+Ln010pTN4fQFYFNEfLKmGsYlLch/3gN4OXD9sMaPiPdGxOKIWEr2N/DtiHjjsMa3gXGea86zs5yOZBptRGwHZi59tQk4NyI2DrMGSV8BLgeeKWla0vHDHJ9sDfBYsjW/a/LbyiHXsA/wHUnXkn1YXhIRjT0t34pxnoH68+wsJ8KXYDQzM6tQMlu0ZmZmo8iN1szMrEJutGZmZhVyozUzM6uQG62ZmVmF3GjNzMwq5EZrZmZWof8DsQ/Ub4TLA6EAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x216 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Results\n",
    "GraphDAG(\n",
    "    est_dag=pc.causal_matrix, \n",
    "    true_dag=dag)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a67be9f0",
   "metadata": {},
   "source": [
    "## Add expert knowledge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "eba4493d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instantiate and encode priori knowledge\n",
    "priori_knowledge = PrioriKnowledge(n_nodes=5)\n",
    "\n",
    "# Add required & forbidden edges\n",
    "priori_knowledge.add_forbidden_edges([(2, 0), (1, 0)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "093f4452",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Train the causal discovery model\n",
    "pc = PC(variant='stable', priori_knowledge=priori_knowledge)\n",
    "\n",
    "# Learn\n",
    "pc.learn(data.X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "d6b54209",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdoAAADRCAYAAABrVHPYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZjklEQVR4nO3dfbBcdX3H8feHKwgWNNRcLSaEoI0ow1TUa9BaleJTEsHUGaYFFCrFRlpR7NgqWp9oa6v1oeCIpldg0OJIKaBGjEVmFKkKkkQQjQEbIcAFNERACCgY8u0f51y6Wfbu7j1nz57f3fN5zeyQvXv2/L435LPf87Tnp4jAzMzMqrFb3QWYmZmNMjdaMzOzCrnRmpmZVciN1szMrEJutGZmZhVyozUzM6uQG23DSDpP0j/VXYeZDY6kN0r6Tt11WGdutBWQtFhSSHpc3bWYNZ2kLZJeUXcd1lxutHOIMv5/ZjYgqWwMp1KHVcMf2n2Q9DRJF0u6S9LNkt6W/3yppPWS7pP0C0mfyN9yZf7feyVtl/SiLusek/RxSdvydZ/Sujcs6QpJH5L0XeBB4OmSTpS0SdL9km6S9OaW9R0uaUrSe/J1bpH0+rZh95X0tfz935f0jIH9ZZklRNJ/AIuAr+ZZfGeer5Mk3Qp8czozbe97dC9Y0m6STpP0M0m/lHShpN/tY+wTJN2Sv+d9bev8oKSLJJ0v6T7gjfnnyVWS7pV0p6RPSdqjZX0h6W155rdJ+mj7hrekj0m6J/8sWV7+b9AGwY22h/wf8leBHwILgJcDb5f0auBM4MyIeCLwDODC/G0vzf87LyL2joirugzxl8By4FDgecCfdFjmeGAVsA9wC7AVOBJ4InAi8G+Sntey/O8B8/N6/xyYlHRQy+vHAqcD+wKbgQ91/Uswm6Mi4njgVuCoiNib/8/oy4BnA6/uYzVvI8vly4CnAfcAZ3V7g6SDgU8Drwf2A55ElsdWK4GLgHnAF4BHgL8hy+6LyD5r/rrtPa8DJsg+K1YCf9Hy2mHAjfn7/xU4R5L6+P2sYm60vb0AGI+If4iIhyPiJuCzwDHAb4HflzQ/IrZHxNUF1v+nZM16KiLuAT7cYZnzImJjROyIiN9GxNci4meR+TbwDeAlbe95X0Q8lL/+tXycaZdExDURsYMs4IcWqNtsLvtgRDwQEb/uY9k3A3+fZ/Qh4IPA0T0O9x4NfDUivhMRDwPvB9pvLH9VRHw5InZGxK8jYkNEXJ3nfAvw72TNvdVHIuLuiLgVOINso3naLRHx2Yh4BPgcWYN/ah+/n1XMjba3A4Cn5Ydz7pV0L/Aesn/AJwHPBG6QtE7SkQXW/zTgtpbnt3VYZpefSVou6WpJd+f1rCDbip12T0Q80PL8lnycaT9v+fODwN4F6jabyzrlbCYHAF9qyf8msr3Pbk1sl1xHxIPAL7vVIOmZki6V9PP8cPI/s2uu298zY67z8cDZToIbbW+3ATdHxLyWxz4RsSIi/jcijgWeAnwEuEjS7/DYLddu7gQWtjzfv8Myj65P0uOBi4GPAU+NiHnAWqD1ENG+eR3TFgF3zKIms1HSKY+tP3sAeML0E0ljwHjL67cBy9s+A/aMiNu7jLlLriXtBTy5R12fAW4AluSno97DrrmGXT8fnOs5wo22t2uA+yS9S9Je+cVLh0h6gaQ3SBqPiJ3AvfnyjwB3ATuBp/ex/guBUyUtkDQPeFeP5fcAHp+PsSO/4OFVHZY7XdIekl5Cdj73v/qoxWwU/YLuWfwpsKek10jaHXgvWcamrQY+JOkAAEnjklb2GPMi4ChJf5hf0HQ6j22a7fYB7gO2S3oW8Fcdlvk7SftK2h84FfjPHuu0BLjR9pCf7ziK7DzmzcA24GyyixuWARslbSe7MOqYiPhNftjmQ8B388NNL+wyxGfJzrFeD1xLtne6g6xhd6rnfrKLMy4kuyjjOGBN22I/z1+7g+wc7MkRccPsfnOzkfEvwHvzw75Ht78YEb8iu+jobOB2sj3c1quQzyTL2Dck3Q9cTXbh0YwiYiPwVuACsr3b+8kuYnyoy9v+lizP95N9LnRqol8BNgDXkV17cU63OiwN8sTvacn3UFdHxAEF3384cH5ELOyxqJkNiaS9yY56LYmImwuuI/L3bx5kbVY979HWLD8cvULS4yQtAD4AfKnuusysHElHSXpCfr3Ex4AfAVvqrcrq4EY7BJJW51+Wb3+sJjtvczrZod5rya5ofH+d9drsSDpX0lZJP57hdUn6pKTNkq5v+86zzVGSXj9Drjfmi6wkO31zB7CE7NSSDyEmroo8+9CxWUmSXgpsBz4fEYd0eH0F2fm6FWTn9s6MiK7n+MysHlXk2Xu0ZiVFxJXA3V0WWUkW2shvajJP0n7Dqc7MZqOKPFdyI+v58+fH4sWLq1h1XzZs2FDb2Cl5/vOfX3cJtdqyZQvbtm3reQu6/CKTbjYCv2l5PhkRk7MoZQG73mhgKv/ZnbNYRy2c5TQ4y/1lGdLMcyWNdvHixaxfv76KVffFt/fM1Pn/IAUTExN9L7vbbjMf3Nm5c+dvIqL/lT1Wp3+Qc+KcjbOcBmd5dvFLLc+emsmMyj/Qp9j1jj4L8R19zCqTWp59jtYaTxJjY2MzPgZgDXBCfrXiC4FfRUTyh43N5qIU8+w9WjPKbQFL+iJwODBf2bymHwB2B4iI1WR3+1pBNiXhg2RTG5pZRVLLsxutGd3P6fSSTyzR7fUA3lJ4ADObldTy7EZrjTd9qMnM5r4U8+xGa4avbjUbJanl2Y3WGk9SqUNNZpaOFPPsRmtGelvAZlZcanl2ozWD5M7pmFlxqeXZjdYaT1JyW8BmVkyKeXajNaPc1wHMLC2p5bmvaiQtk3RjPv/eaVUXZTZsu+2224yPUeIsWxOklueee7SSxoCzgFeS3eNxnaQ1EfGTqoszG4YUr1KsgrNsTZBinvupZimwOSJuioiHgQvI5uMzGxmpbQFXxFm2Rkgtz/2MOtPce2YjY/oCik6PEeIsWyOklud+Lobqa+49SauAVQCLFi0qWZbZ8KR4qKkizrKNvBTz3E81fc29FxGTETERERPj4+ODqs9sKFI71FQRZ9kaIbU897NHuw5YIulA4HbgGOC4SqsyG6IUt4Ar4izbyEsxzz0bbUTskHQKcBkwBpwbERsrr8xsiFILZhWcZWuK1PLc1w0rImIt2WS3ZiMptWBWxVm2Jkgtz74zlDVeioeazKyYFPPsRmtGelvAZlZcanl2o7XGk5TcbB9mVkyKeXajNSO9+SvNrLjU8uxGa42X4hawmRWTYp7daM1Ib6JoMysutTy70ZqR3qEmMysutTyndWmWWQ2mDzXN9OhzHV3neZX0JElflfRDSRslnTjwX8TMSue5iiy70ZpRbraPlnlelwMHA8dKOrhtsbcAP4mI5wCHAx+XtMdgfwszg+J5rirLPnRsjTeAiycenec1X9/0PK+tE6oHsI+ypO8N3A3sKDOomT1WyTxXkmU3WjNKf8G90zyvh7Ut8ylgDdlsOfsAfxYRO8sMamadlchzJVkeyUYb8ZgpNhup7gsC5sr/hz5u2TZf0vqW55MRMdm6ig7vaf/lXw1cBxwBPAO4XNL/RMR9BUpujLnyb6hqznL/Sua5kiyPZKM1m60eh5q2RcREl9f7mef1RODDkX1ibZZ0M/As4JoC5ZpZFyXyXEmWfTGUNd70FnCJiaIfnec1vyjiGLJDS61uBV6ej/dU4CDgpgH+GmZG6TxXkmXv0ZpR7hztTPO8Sjo5f3018I/AeZJ+RHZ46l0Rsa185WbWrmieq8qyG60Z5Wf76DTPax7K6T/fAbyq1CBm1peSG84Dz7IbrTVeivdGNbNiUsyzG60Z9V/VaWaDk1qe3Wit8VLcAjazYlLMsxutGeltAZtZcanl2Y3WGq+PL7ib2RyRYp7daM0of9WxmaUjtTy70ZqR3qEmMysutTy70VrjpXjxhJkVk2Kee+5fSzpX0lZJPx5GQWZ1KDMf7VziPFsTpJbnfg5knwcsq7gOs1qVvNfxXHIezrONuNTy3PPQcURcKWnxEGoxq8Uo7rnOxHm2UZdingd2jlbSKmAVwKJFiwa1WrOhGME918KcZZvrUsvzwKqJiMmImIiIifHx8UGt1mwoUjvUVCdn2ea61PLsq46t8VI81GRmxaSYZzdaM9I71GRmxaWW536+3vNF4CrgIElTkk6qviyz4Urt6wBVcZ6tCVLLcz9XHR87jELM6pLivVGr4jzbqEsxzz50bEZ6t2wzs+JSy7MbrRnpndMxs+JSy7MbrTVeioeazKyYFPPsRmtGeoeazKy41PLsRmtGeoeazKy41PLsRmtGelvAZlZcanl2o7XGS/GcjpkVk2Ke06rGrCZl740qaZmkGyVtlnTaDMscLuk6SRslfXugv4CZPapMnqvIsvdorfHKbgFLGgPOAl4JTAHrJK2JiJ+0LDMP+DSwLCJulfSUclWbWSdl8lxVlr1Ha0bpPdqlwOaIuCkiHgYuAFa2LXMccElE3AoQEVsH+guY2aNK5LmSLLvRmtHz3qjzJa1veaxqe/sC4LaW51P5z1o9E9hX0hWSNkg6obrfxqzZSuS5kiz70LE1Xh+HmrZFxES3VXT4WbQ9fxzwfODlwF7AVZKujoifzqpYM+uqZJ4rybIbrRmlvw4wBezf8nwhcEeHZbZFxAPAA5KuBJ4DuNGaDViJPFeSZR86NqP0Odp1wBJJB0raAzgGWNO2zFeAl0h6nKQnAIcBmwb6S5gZUCrPlWTZe7TWeGWvOo6IHZJOAS4DxoBzI2KjpJPz11dHxCZJ/w1cD+wEzo6IHw+gfDNrUSbPVWXZjdaM8rdsi4i1wNq2n61ue/5R4KOlBjKznkpuOA88y260ZqR3b1QzKy61PLvRWuO1XPZvZnNcinl2ozUjvS1gMysutTy70ZqR3mwfZlZcanl2o7XGk8TY2FjdZZjZAKSYZzdaM9LbAjaz4lLLsxutGemd0zGz4lLLc89qJO0v6VuSNuVz7506jMLMhmX6C+5l5qOdC5xla4IU89zPHu0O4B0R8QNJ+wAbJF3eOj+f2VyX2qGmijjL1gip5blno42IO4E78z/fL2kT2bRBDqeNjFHac52Js2xNkVqeZ3WOVtJi4LnA9yupxqwGKX7BvWrOso2qFPPcd6OVtDdwMfD2iLivw+urgFUAixYtGliBZsOQ2tcBquQs26hLLc997V9L2p0smF+IiEs6LRMRkxExERET4+Pjg6zRrHLTW8GdHqPEWbYmSC3PPfdolVV2DrApIj5RfUlmw5XiF9yr4CxbE6SY5372aF8MHA8cIem6/LGi4rrMhiq1rwNUxFm2Rkgtz/1cdfwdYLSOn5m1GMVDxJ04y9YEKebZd4YyI72LJ8ysuNTy7EZrjZfiFrCZFZNint1ozUjvC+5mVlxqeXajNSO9YJpZcanl2Y3WGm/6JuRmNvelmGc3WjPS2wI2s+JSy7MbrTVeilvAZlZMinl2ozUjvS1gMysutTynVY1ZTcreG1XSMkk3Stos6bQuy71A0iOSjh5Y8Wa2izJ5riLL3qO1xit7b1RJY8BZwCuBKWCdpDXtE6rny30EuKxEuWbWRZk8V5Vl79GaUXqPdimwOSJuioiHgQuAlR2WeyvZzDlbB1e5mbUrkedKsuw92hEWEbWOn9rdWbrpcU5nvqT1Lc8nI2Ky5fkC4LaW51PAYa0rkLQAeB1wBPCCctVa0zjLs1Miz5Vk2Y3WGq+PQ03bImKi2yo6/Kz9k/EM4F0R8chc+9Aym0tK5rmSLLvRmlF6i30K2L/l+ULgjrZlJoAL8nHmAysk7YiIL5cZ2Mweq0SeK8myG6013gAmil4HLJF0IHA7cAxwXOsCEXFgy3jnAZe6yZoNXsk8V5JlN1ozyu3RRsQOSaeQXYE4BpwbERslnZy/vnowVZpZP4rmuaosu9GaUf5ij4hYC6xt+1nHUEbEG0sNZmZdldxwHniW3WjNmHtXVZrZzFLLsxutNV6K90Y1s2JSzLMbrRnp3RvVzIpLLc9utGakd6jJzIpLLc9utNZ4s5k8wMzSlmKe3WjNSG8L2MyKSy3PbrRmpHdOx8yKSy3PPRutpD2BK4HH58tfFBEfqLows2FKbQu4Cs6yNUVqee5nj/Yh4IiI2C5pd+A7kr4eEVdXXJvZUKT4dYCKOMs28lLMc89GG9n8TNvzp7vnj3rnbDIbsNS2gKvgLFtTpJbnvtq+pDFJ15FNcnt5RHy/wzKrJK2XtP6uu+4acJlm1So58fuc4SxbE6SW574abUQ8EhGHkk0ZtFTSIR2WmYyIiYiYGB8fH3CZZtWZPtQ002OUOMs26lLM86xGjYh7gSuAZVUUY1aX1LaAq+Ys2yhLLc89G62kcUnz8j/vBbwCuKHiusyGKrUt4Co4y9YUqeW5n6uO9wM+J2mMrDFfGBGXVluW2XCN6p5rG2fZGiG1PPdz1fH1wHOHUItZLUb5EHErZ9maIMU8+85QZqR3JxkzKy61PLvRmpHeoSYzKy61PLvRWuOleCcZMysmxTynVY2ZmdmI8R6tGemd0zGz4lLLc1rVmJmZjRjv0ZqR3hawmRWXWp7TqsasBt1u19bv1YuSlkm6UdJmSad1eP31kq7PH9+T9JyB/yJmVjrPVWTZe7RmlPs6QH6npbOAVwJTwDpJayLiJy2L3Qy8LCLukbQcmAQOK1Gymc2gaJ6ryrIbrRmlv3e3FNgcETfl67oAWAk8Gs6I+F7L8leTzZ5jZhUokedKsuxGa5XJ5hmvz8TERN/Llmy0C4DbWp5P0X0L9yTg62UGNBumuZRlKJXnSrLsRmuN18cX3OdLWt/yfDIiJltX0eE9HT+ZJP0xWTj/aNaFmllPJfNcSZbdaM162xYR3Tapp4D9W54vBO5oX0jSHwBnA8sj4peDLdHM+tQtz5Vk2Y3WjNJfB1gHLJF0IHA7cAxwXOsCkhYBlwDHR8RPywxmZt2VyHMlWXajNaPcOdqI2CHpFOAyYAw4NyI2Sjo5f3018H7gycCn87F29NhLNrOCiua5qiy70ZpRfraPiFgLrG372eqWP78JeFOpQcysLyU3nAeeZTdaa7wUJ4o2s2JSzLMbrRnpzV9pZsWllmc3WjPSuzeqmRWXWp7TqsbMzGzEeI/WGq+PL7ib2RyRYp7TqsbMzGzEeI/WjPTO6ZhZcanlue9qJI1JulbSpVUWZFaHsvPRziXOso261PI8m7Z/KrCpqkLMbGicZbMh6qvRSloIvIbsJspmI6Xb1u+o7dE6yzbqUsxzv+dozwDeCexTXSlm9Rm1htrFGTjLNuJSy3PPPVpJRwJbI2JDj+VWSVovaf1dd901sALNhiG1LeAqOMvWFKnluZ9Dxy8GXitpC3ABcISk89sXiojJiJiIiInx8fEBl2lmA+Asm9WgZ6ONiHdHxMKIWEw2N983I+INlVdmNkSpbQFXwVm2pkgtz/4erRnpndMxs+JSy/OsGm1EXAFcUUklZjUZtT3XfjjLNqpSzHNat88wMzMbMT50bEZ6h5rMrLjU8uxGa0Z6wTSz4lLLsw8dm5mZVch7tGaktwVsZsWllmc3WjPSC6aZFZdant1orfFS/DqAmRWTYp7daM1IbwvYzIpLLc++GMpsACQtk3SjpM2STuvwuiR9Mn/9eknPq6NOM+uuiiy70ZpR7t6oksaAs4DlwMHAsZIObltsObAkf6wCPjPY38DMphXNc1VZdqM1o/RNyJcCmyPipoh4mGxmnJVty6wEPh+Zq4F5kvYb7G9hZlAqz5VkuZJztBs2bNgm6ZYSq5gPbBtUPXNw/BRqqHv8QdRwQD8Lbdiw4TJJ87sssqek9S3PJyNisuX5AuC2ludTwGFt6+i0zALgzn5qrMsIZDmFGuoeP4UahpJlKJ3nSrJcSaONiFKTWEpaHxETg6pnro2fQg11jz/MGiJiWclVdNpMjgLLJGeuZzmFGuoeP4Uahjl+yTxXkmUfOjYrbwrYv+X5QuCOAsuYWb0qybIbrVl564Alkg6UtAfZpOpr2pZZA5yQX7H4QuBXEZH0YWOzBqoky6l+j3ay9yIjPT7UX0Pd40MaNfQUETsknQJcBowB50bERkkn56+vBtYCK4DNwIPAiXXVO2Qp/D+su4a6x4f6a6h7/L5UlWVFJH+ayMzMbM7yoWMzM7MKudGamZlVKKlG2+vWV0MY/1xJWyX9eNhj5+PvL+lbkjZJ2ijp1Bpq2FPSNZJ+mNdw+rBryOsYk3StpEvrGN/Kc57rzbOznI5kGm2ft76q2nlA2e9UlrEDeEdEPBt4IfCWGv4OHgKOiIjnAIcCy/Ir64btVGBTDePaADjPQP15dpYTkUyjpb9bX1UqIq4E7h7mmG3j3xkRP8j/fD/ZP84FQ64hImJ7/nT3/DHUK+YkLQReA5w9zHFtoJznmvPsLKcjpUY7022tGknSYuC5wPdrGHtM0nXAVuDyiBh2DWcA7wR2DnlcGxznuUVdeXaW05BSo52Tt6irgqS9gYuBt0fEfcMePyIeiYhDye54slTSIcMaW9KRwNaI2DCsMa0SznOuzjw7y2lIqdH6FnWApN3JQvmFiLikzloi4l7gCoZ7nuvFwGslbSE73HiEpPOHOL4NhvNMOnl2luuVUqPt59ZXI03ZHE7nAJsi4hM11TAuaV7+572AVwA3DGv8iHh3RCyMiMVk/wa+GRFvGNb4NjDOc815dpbTkUyjjYgdwPStrzYBF0bExmHWIOmLwFXAQZKmJJ00zPHJtgCPJ9vyuy5/rBhyDfsB35J0PdmH5eUR0djL8q0Y5xmoP8/OciJ8C0YzM7MKJbNHa2ZmNorcaM3MzCrkRmtmZlYhN1ozM7MKudGamZlVyI3WzMysQm60ZmZmFfo/NyHQOq+iHAMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x216 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Results\n",
    "GraphDAG(\n",
    "    est_dag=pc.causal_matrix, \n",
    "    true_dag=dag)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5f7c0d4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:causal_book_py38]",
   "language": "python",
   "name": "conda-env-causal_book_py38-py"
  },
  "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}