notebooks/partial_ancestral_graphs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Creating and Handling a Partial Ancestral Graph (PAG)\n",
    "\n",
    "First, let's import the PAG class and edge-mark types. In a PAG, three edge-mark types exist:\n",
    "1. Arrow head ...---> X\n",
    "2. Tail ...---- X\n",
    "3. Circle, a variant edge-mark the represents 'head or tail' ...---o X\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('..')\n",
    "\n",
    "from graphical_models import PAG, arrow_head_types as Mark"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "For plotting the graph we need to import:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "from plot_utils import draw_graph"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Creating a PAG\n",
    "For creating a PAG  we need to first define a set of nodes. This set elements can any immutable object that returns a 1-3\n",
    "length discriptive string when str(object) is called. Commonly numeric values or short strings are used.\n",
    "In this example we in use short strings. The graph that we are creating is in Figure 2-(a) by Jaber et al. (2019)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "nodes = {'X1', 'X2', 'Y1', 'Y2', 'Y3', 'Y4', 'Y5'}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Next, we instantiate the PAG and set the graph edges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "pag = PAG(nodes)\n",
    "pag.add_edge('X1', 'X2', Mark.Directed, Mark.Directed)\n",
    "pag.add_edge('X1', 'Y1', Mark.Tail, Mark.Directed)\n",
    "pag.add_edge('X1', 'Y3', Mark.Directed, Mark.Directed)\n",
    "pag.add_edge('X2', 'Y1', Mark.Directed, Mark.Directed)\n",
    "pag.add_edge('X2', 'Y2', Mark.Tail, Mark.Directed)\n",
    "pag.add_edge('Y1', 'Y5', Mark.Directed, Mark.Circle)\n",
    "pag.add_edge('Y2', 'Y3', Mark.Directed, Mark.Directed)\n",
    "pag.add_edge('Y3', 'Y4', Mark.Directed, Mark.Circle)\n",
    "pag.add_edge('Y4', 'Y5', Mark.Circle, Mark.Circle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Calculating PAG Properties\n",
    "\n",
    "In this example we find the visible edges (Richardon and Spirtes, 2002) and dc-components (Jaber et al., 2018)\n",
    "of the PAG."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[{'Y3', 'X1', 'X2', 'Y2', 'Y1'}, {'Y5'}, {'Y4'}]\n"
     ]
    }
   ],
   "source": [
    "dc_components = pag.find_definite_c_components()  # dc-components\n",
    "print(dc_components)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X2 ---> Y2\n",
      "X1 ---> Y1\n"
     ]
    }
   ],
   "source": [
    "pag.find_visible_edges()  # find 'visible' edges\n",
    "for edge in pag.visible_edges:\n",
    "    print(str(edge[0]),'--->', str(edge[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Plotting a PAG\n",
    "\n",
    "Finally, for plotting the PAG, rendering color-coded edges, we use the following method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS4AAAEuCAYAAAAwQP9DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA8uklEQVR4nO2dd5iU1dn/P2dny/RtgIBUpWhQ0B0MoKIGRSTGHpQoplgJvsZXo0QxwfITSyzY8AXEGDQ2YuwGIbYIRBFmQEEpKiwSihTZ3qbcvz+eWdwKW2bn7Mycz3XdF+zT5nvOM893zjnPKUpEMBgMhkQiTbcAg8FgaC3GuAwGQ8JhjMtgMCQcxrgMBkPCYYzLYDAkHMa4DAZDwmGMy2AwJBzGuAwGQ8JhjMtgMCQcxrgMBkPCYYzLYDAkHMa4DAZDwmGMy2AwJBzGuAwGQ8JhjMtgMCQcxrgMBkPCka5bgMFgsFBKKaAXkA2Ega0iUqZXVedEmRlQDQZ9KKXswITc3NzfVlZWDs3KykrLzs4OhkIhtXv3brvD4dgViUTeKSsre0RE1ujW21kwxtVJUUplAfmAAvaJSIVmSYYYopRSNpvtsszMzIeGDx+e9rvf/c594okn0r179/3HBINBvvjiC1599dXQY489FoxEIv7i4uJfishmjdI7Bca4OhFKqaPdbvf/pKenn1pWVtbX5XLVKKUoLy/PdDqdO0VkSUlJyRPAMjE3LmFRSuV6vd5Xe/bsOfzFF190DRs27KDnBINBHnzwwdCdd95ZEwqFptTU1MyPg9TOi4iY0BzAEV6vd3leXl759OnTgx9//LFUVFRILTU1NbJq1Sp54IEHwr169Sr1eDxfAcfr1m2iTfc61+12b5gyZUpVTU2NtJa1a9dKt27dyjMzM6/VnRat+ahbQKpHVlbWdU6ns+LRRx8Nt+SLHIlEZMGCBZKdnV3hdrtnAjbdaWhn7DxooptmZyfQ3qoA0rxe73+uueaa6kgk0sZki2zevFny8/PLgTN0p0lbXuoWkMrhdDrv6d+/f/nXX38trWX37t0yatSoco/H80qCm1d70K29VZGenv7bo48+uiwYDLYz2SLvvvuuOJ3OPUC27nTpCO0CUjXS09Mv79u3b9muXbukrVRWVsqoUaPKXC7XQ7rT045oD7q1tzgAj91uL//iiy9ERCQUCsnAgQNl8uTJjRI1depU6dGjh3z44YcyadIk6d+/v9jtdhkwYIDMnj17/3G/+tWvKh0Ox4O606YlP3ULSMUA+jidzrI1a9ZIe9m1a5dkZ2dXACN1p6uN0R50awfIBQYf7Dil1OTx48eX1RX/9NNPi91ul507f6gtz5o1SzwejwQCAbn11lvlzjvvlKVLl8o333wjd9xxhwCyaNEiERHZuHGjOByOEiBLdz7EPd91C0jFyM7Ofv32228PilhtVkOGDJHrr79e6rJkyRKx2+0ye/ZsWb9+vVx11VUyaNAgcTgc0qdPH7nzzjultp3kxRdfFK/Xu4HoW+IEi/agWzvANYAA/wJ+3Nxxubm5axYuXFhPfDAYlP79+8vNN98sIiKvvfaa2O12Wbx4cbMJttvtMn369P1/jxgxohg4V3c+xD3fdQtItQB6OByOqqKiIqnl+eefF6fTKbt37xYRkfXr10teXp5MmzZNREQeeeQRuemmm+SDDz6Qb775RubMmSNpaWkyZ84cEREJh8PSo0ePsgQtdbUH3doBfgtUAhGgHPgYOKXujwiQkZGRUVNaWtooAXPmzJHs7GxZtGiROJ1OmT9/frOJXb58uQDyj3/8Y/+26dOnh7Oysu7TnQ9xz3fdAlIt0tLSbvzVr35VKXUIh8MycOBAueWWW2Tnzp3Sv39/mTRpkhyIwYMHyy9/+cv9f997771hj8czX3f62hD1mDBhgthsNqlt+7vpppsEEL/f31Q26NYOcGXUsKROlAGfAWdidSA+6tBDDy1pKgHV1dXSu3dvSUtLkxkzZjR1iIhYbxL79+8vZ511ltR9I/n6669Lfn7+Mt35EPd81y0g1SIvL++dv/3tb9KQp59+WrxerxQUFMiYMWOkurq60TG1fPvtt5KZmSkPPvjg/m3Lly+XnJyczbrT14aox2uvvSbA/tJk//795cgjj2wuK9r0mVEzyQAcgCfaTtUN6An0AQ4DBgFDgGGADxgBnBAtTZ0GnAH8DHgUqGpgXLVRBewB7h06dGhRc4n4zW9+I8OHD29ut3zyySfSvXt3mTRpUqPvxX/+8x/Jy8tb1wH3pVOHGWQdZ4LB4DEFBQWNtk+cOJFrr72WkpIS3n//fTIzM5s8//vvv+e8885j6NChTJkyZf/2oUOHUlZW1ksplSUi1S3VEx3Ya8MacJ8R/bduNNzWkmNafK2okexn/Pjx5OXl8fLLL+Pz+di8eTMzZsxoTvt/2qjDBoSaiGAbtnWPXq8pbFjGmBOJRJo5BD7//HNGjhzZ5L5nnnmGyZMnc9ttt/GHP/yh0f5wOAzWgOyUwhhXnKmpqfEecsghjbZfe+21hEIh9u7dS1pa07MNbdy4kbPPPps+ffrw97//Hbvdvn9f9P8K+DRqRi19oG1YX/z2PLxtPS/YMI2ZmZlMmDCBp556irlz56KU4pJLLmkuO29so9awNHTMNqKUmgjMxTKoWiqx8vQe4BEgf+vWrb9s6vxgMMjatWu55ppr6m2PRCLcfPPNzJ49m5deeomzzjqryc/fvHkzwJZ2JyTR0F3kS7XIysoq37Nnj9Tl9ttvl9zcXFm1apX06tVL7r77bmnIokWLJCcnR6699loJhUKN9ouI2O32amAMVvXmR1jVncOwqj89ga5Y1SIPVjUpHf1vIhvx0UcfCSBKKRk9enSTaY2i/X4CPweK+aFqWA7MoE7HUEDZ7fbSbdu2NUpAIBAQQFavXr1/W2lpqfz0pz+Vrl27yrvvvis7duzYH3WHgomIXHPNNdVKqZt050Pc8123gFSL7OzsLYFAQGqZN2+e2O12+eijj0RE5NFHH5WuXbtKeXn5/mNmzpwpDodD5s6dK81RXFwsGRkZNUC67jS2MhoRiUSkb9++9dq6mkG3doBzgRqsUtYjQJemjsvJyVk4e/bsRuN85s2bJ1lZWVK3N/2LL77YVHuZAPLWW2/tPy4cDsuhhx5aCozQnQ9xz3fdAlItcnJyXn7yySdFROTtt9+WzMxMWbBggdRSWVkp3bt3l4ceekiCwaBcdtll4nK55KWXXqr3y9vw1fqHH34oubm5X+hOXxuiPejWTrRE+wDQ8yDHnT5gwIDS9oxRbMiiRYvE6/V+3QlKzfHPd90CUi2AK8aPH1/26aefisvlkoceekga8sADD0jPnj3lk08+afaX9/HHH693zo033ljjdDpn6k5fG6I96Nbemvue5na7Nz333HMxca5gMChDhgwpU0r9SnfadISZjyvOKKU8drv9u6+++srRq1evmFyzurqabt26VZaUlAwTka9ictH40Z4voIqZijiglPqx1+v9cMOGDY66Ewa2hbvuuit0//33rygpKTlBUvAhNotlxBkRKU1PT//rrbfeWhWraz788MMhpdTyBDStlEJEPg2FQg+MGTOmfN++fW2+zosvvij33HNPWUlJyS9S0bQAU1XUEYDX5XLtbjh2rS18+eWX4nQ6y4D+utPVxmgPurW35d4rl8v12OGHH162du3aViU2FArJvffeG3I6nfuAobrTojUfdQtI1QDGeDyeis8++0zayrZt26Rnz57l6enpl+tOTzuiPejW3tZ7r9LT069yOp1lf/rTn4K1Y1SbIxKJyJIlS+SYY44p83q9KxL4RypmYdq4NJKWlnah2+3+60svveQYP358q871+/2cddZZFfv27ZtRWVl5dwdJjAc7gcY9cg/Od1i91hMWpVQfr9d7f01Nzdnjx48Pn3rqqa5jjz2W/Px8QqEQX3/9NStWrAi/8MILlbt37y6prKy8MxQKPSkizXfDTxGMcWlGKXWKy+VacN5553lmzJhh79OnzwGP37NnD/fff3/w8ccfr6qsrLw6Eom8ECephg5CKZUPXOD1ek+y2WzHhcNhr1IqpFT+f8vLt34QDAbfBT4Q87DuxxhXJ0Cpwz2ZmWP/mpY2/4wTTzwxct5557l9Ph99+vRBKcXOnTvx+/288847FW+99VZaZmbm6yUlJdeLyA7d2g0dg1KMA/4J9BFhm249nQ1jXJpRilysSeiOBeUFJmRnZ49VSo2sqanpIiIqMzOzSCm1sqio6H2oOBMcmcA4EUJ61Rs6CqWYB1wOfAiMEWlXt5GkwxiXRpTix8CbWGMIwyJktOCc/wCjgPeBc0QwS7QnGUqRDnyPNaa0Epgiwl+1iupkmH5cGlAKpRS/x/o17YbVkbKlDa619+xEYKVS9Ii9QoNmTq3zfwfwmFIcqktMZ8QYV5xRihxgIXAn1pdy/64WXiIv+m8mcDjwmVIcFTOBhs7AFYC7zt924DmlEmukQEdijCuOKMVwYD3WLJrONl4mu87/04EuwMdKcVr71Bk6A0rh4Icpn2tJB4YDv9EiqhNijCtOKEUe1kIKhwBZ7biUp8HfCuvX+Q2luKId1zV0Dn4GTb50cQGPKkVsBrgmOMa44oQI3wPjgFVYk821mmijrb2Z3TbgcaXqlcgMicdVNP5xqiULa477lMe8VYwz0XaKU4DXsdqpaktfoYO9VVSKLsBW6ptXFVACzASeEmF3rDUb4oNS2IFSoBrrexHEmmL7Q+BL4CvgYxH8ujR2Fsyc83FGBFGKINbqL1cC92JNs9z06hj1qW2Yr8F6C7kZ68tdIJJ6CyYkIdXA6cBeoBCowCqdnylCjUZdnQ5jXHq4DbhbhPeA45TiJOCnLTivds2+vwGPAf8FvsWaW35dB2k1xIloJ9MP6m5Tih1Ab+AbLaI6KaaqGGeU4njgOWBwLH5FleIuwCvC79otztDpUIoPgLuiP3KGKKZxPv7UlrZiVfSfA0xSql6/H0PyUAj006yh02GMK44oxShgMDA/VtcUYSvwb6DZxQcNCU0hxrgaYYwrvtwG3NMBDa1PAFNMz+qkZDPQX7eIzoYxrjihFCOAI4GnO+Dy72ENHzq+A65t0EshpsTVCGNc8aOjSluIEAH+D5gS62sbtFOIMa5GmLeKcSA6fc3LwEARqjvoM3KxqhWDRNjVEZ9hiD9KYcPqz+XtqO9OImJKXPGhtrTVYV88EfYB/wAzXjGZiHYs/i9w4Dm9UwxjXB2MUhwHDAX+EoePmwVMjv5KG5KHQkx1sR7GuDqe24B741HMFyEAbMeaFsWQPGzGGFc9jHF1INH5t44Bnorjxz6BaaRPNgoxXSLqYYyrY5mOVdqqiuNnLgB8SjEgjp9p6FgKMSWuehjj6iCUwgcUAPPi+blRk3wamBzPzzV0KIUY46qH6Q7RQSjF68B7IvGf+E0pDgM+xVqTryLen2+ILdGFMlaKmIVRajElrg5AKY7FmiP8SR2fL8ImYDlwkY7PN8ScHUBudD56A8a4OorpwJ9FqNSowTTSJwnRkRFbMX259mOMK8YoxTHACGCuZinvAF2i/cgMiY/pElEHY1yxpzOUtmp7XM/GlLqShUJMl4j9GOOKIUoxDBiJNblfZ+AvwHlKka9biKHdFGJKXPsxxhVbpgMP6C5t1RJd8edNzEKiyYCpKtbBGFeMUIqhWPNhzdatpQGzgN8qZe51glOIqSrux3yZY8efsEpbna3f1HKgGGvZK0PiUogpce3HdECNAUpxNPAv4HCRtq1S3ZEoxRXA2SKcrVuLoW1Ep+WuAPI74Y9j3DElrtjwJ+DBzmhaUV4ATlDK/GInKtE1F7dgSl2AMa52oxRDgJOwOnx2SqKG+ixwlW4thnZRiDEuwBhXLPgT8FAnLm3V8n/A5UqRpVuIoc2YN4tRjHG1A6X4EfATOnFpqxYRNgCfAz/XrcXQZgoxxgUY42ovtaWtMt1CWogZv5jYFGK6RADGuNqMUhwJjMHqJ5UovAn0jY6nNCQepqoYxRhX2/kT8HAClbYQIYQ1HOm3urUY2kQhxrgA04+rTSjFEcBHWP22SnXraQ1K0QP4EugnQrFuPYaWE+3LVQ50S6QfzI7AlLjaxh+xSlsJZVoAIuwAFgO/1K3F0DqifbkKMaUuY1ytRSkGA+OAx3VraQezgCnRX3BDYmHauTDG1Rb+CDwiQoluIe1gCRDG6sphSCwKMcZljKs1KMUg4AzgMd1a2kO0ymG6RiQmhZguEca4WsmtwKNJ0qj9N+DU6AoyhsTBVBUxxtVilGIg8FOI/3JjHUG0qvsCcKVuLYZWUYgpcZnuEC1FKf4KbBbhDt1aYoVSHIX1hrGvCMHotgFAkQh7tIozNIlSdAE2ipCnW4tOTImrBUQf5jOBR3RriSUirAW+As5XivOVYnn0bzPVc+dlL5ChFNm6hegkXbeABOFWYJYIRbqFxJJoZ9TdWO1dlYAn+m8ytOElJSKIUhQCfbEGzackxrgOglIcDpwFDNCtJVYoxQisIUunAYL1PfBEdwexftUNnZdCrHYuY1yGZplG8pW2ZgCnALYm9gnwfVzVGFpLyr9ZNMZ1AJTiMOAcYKBuLTHmHOBV4ATA2WCfwhhXZ6eQFDcu0zh/YKYB/yfCPt1CYkl0ttbxwDPQaOEFG8a4OjuFpHiXCFPiaobowhLnkXylLQBECGOtt/gVcBfgiO7KxBhXZ6eQFC9xmX5czaAUc4FdIvxRt5aORinOAZ7HqjZGgPTosCBDJ0QpcrH6FObo1qILY1xNEC1t+YFBIqnxhk0pjsOaY8wmQqa1TXUH+mBVH4uAr0QkpE2kAdg/L1cR1pxqSdWM0VKMcTWBUswB9ohwq24t8USp4FGw4yGv9+jycDh8iog4evfuXWWz2di3b59tz549mW63e0NRUdETIvI3EUnpyex0ohSfAb8WYZVuLTowxtUApegLBEih0haAUupor9f7QnZ2dt8bbrjBec4556T169cPpX6YsqukpIRly5bx2GOPlX/wwQcqHA7fEgwGHxeRiEbpKYlSvAbMF+FV3Vp0YIyrAUoxG9gnwi26tcQLh8Pxe5vN9v8efvhh++WXX67qmlVzrF+/nokTJ5Zv3rz5i5KSkp+KSMqYfGdAKR4GtogwU7cWHZjuEHVQij7AhcCDurXEC6fTeVf37t3v/OKLLxxXXHFFi0wL4IgjjsDv97suv/zyY9xu96dKqfwOlmqoTyEp3CXCGFd9bgbmpsrMCGlpaRfl5eVd/8knnzj79u3b6vNtNhsPPvhg5hVXXNHL6/W+pZQy36f4UUgKd4kwVcUoStEbWA0cIcJuzXI6HKXUIU6nc+OHH37oPe6449p1rXA4TEFBQfm6dev+UFNTk0jrTCYsSjEMeFaEobq16MD8Qv7AzcBTqWBaAG63+66rrrrKcdxxxyEiHHXUUdxwww31jlm6dCkOh4M5c+ZQU1PDFVdcwfDhw7Hb7Qwd+sPzYrPZePHFF102m+3PSilXvNOSomwB+qfsgicikvIB0gvke5BuurXEJ7147HZ7xX//+1+p5fnnnxen0ym7d+8WEZH169dLXl6eTJs2TUREtm3bJrfddpu8+uqrMnToUPn1r38tDTn99NNLlVJX6k5fqgTIPpB83Tq0pF23gM4QII+B3K9bR/zSyyWnnnpqidQhHA7LwIED5ZZbbpGdO3dK//79ZdKkSdKQSCQiHo9HHnvssUb73nrrLcnNzV2tO32pEiABEJ9uHToi5ccqRheLuAQ4UreWeOF0Ok8cO3asu+62tLQ0pk2bxnXXXceiRYvo378/Tz31VKNzN27cSGlpKT6fr9G+E044gfLy8iOUUjYRCXdcCgxRCrHeLPo164g7po0L/gD8RYTvdAuJF3a7fWRBQUGjtpGJEycSiUQoKSnhlVdeITMzs9G5fr8fm83GsGHDGu3LyckhNzc3SBJNutjJKSRF3yymdIlLKXoCk0ih0haAiHjz8xt3u7r22msJhULs3buXtLSmf9MCgQBHHnkkTmfDabwscnNzw999911OLPUamqUQGKRbhA5SvcT1B+CvqVTaihIJh+vX5O644w7+8Y9/8PHHH+NyuXj88cebPNHv91NQUNDshUOhEIAZiB0fNpOinVBT1riiC0VcCvxZt5Z4o5Qq/Oabb/b//dRTT3Hvvffy+uuvc8wxxzB16lRmzpxJRUX9OQZFhFWrVjXZvlW7f9u2bQ7g247Ub9hPISlaVUxZ4wKmYg1S3albSLwpKir69/Lly4MA//znP5kyZQrPPPMMo0ePBuDKK6/EZrMxZ84cAL777jtWr17NO++8Q3FxMU6nk9WrV7N3b/3hid988w02m61URHaDNf2KUuQrxXFKMVEpblWKF5TCrxTnxzXRyUkh0C8l+3Lpfq2p6TVy92i/rR66tehJP+Py8/MrP/nkE3G5XPLQQw9JQx544AHp2bOnVFVVydVXXy1Yi2jUi7feeqveOffff3/Y6bz4fZCvQLaCVIJUgxSBlIAEQQSkCuR03fmQDAGyJ1X6H9aNlBzyoxQPYU2Yd51uLfEiOnp6JDAZODsjI8P29ttve8aOHRuT60ciEXr16lW+Y8dZl8Kc+fyw3FlTlAJ5IqYtrL0oxUpgigif6tYST1KuqqgUhwC/Bu7TLCUuKKWylVLXAJ8B84E1wMBQKPQ/119/fXm0Mb3dzJ49O1JeXv41zH0N6y3tdqCpvlwR4GVjWjGjkBRs50o54wJuAv4mwnbdQjoKZXGcUmoe1hf7JOB/gcEi8oCI7BGRZ7/99tvVd999d7sNZNOmTUydOrWqpKTkFyIiImwDRmCtkt3QvBRwvFL8Pjqw3dA+UvLNYkoZV7S0dRlJWtpSSnmUUldh9aR+CfgaOEJELhKR96VOu4CISGlp6S/uu+++kmeffbbN7QU7duzglFNOqQiHwzeLyLofrs9/scxrD1Ypq5Zi4HdYpbLVSrFUKf5HKbq3VUOKU0gKlri0N7LFM0DuB3lMt47Yp4tjgdlYy4q9AowD0lp47hCn07n39ttvD9bU1EhrWLZsmXTv3r3c4XD86QB53g9kF0gYJATyRJ19mSBngjwbHTD8HsiVqTpwuC0Rzb+FunXEPd26BcTxBncD2QvSq4M/a6e0nZ0tTw8urNLjcqx+U38CerYtb+jt9Xo/OvLII8sWLVok4XD4gCK3bNkiU6ZMqXY6nUXABQe/vhwOshukBmRkM8c4QM4HWQBSDPJPkF+CZLclTakSIENA1unWEe9ImbeKSvFnwCXCNR38Ue3N0AP2yVFKHQVcDVwMLAPmAO9IOwc1R9vFJnk8nts8Hs8hF1xwQeaPf/zjzMMPPxybzcbevXsJBAKyePHisuXLl9tsNtsz5eXl0yXaZ+vg12cAcAtwhciB80gp3MBZwETgFOB94EXgLbFW4TZEUQoXVnXcebB8TSZSwriUoiuwARgmwtYO/riYG5dSygFMwDKsfsA8YJ6IxDwt0W4TP05LSxuTk5NzMtBPRNKVUkWVlZXLKisr/wO8JSJxMRClyAHOxTKxkcA7WCb2jghV8dDQ2VGKXcBQSaHO1KliXPcCXhGmxOHjYmZcSqkjsMzqUmAFVunqLUnRRVmjP0DnY5nYMOBNLBN7V4SgTm06UYpPgd+J8IluLfEi6Y1LKbpglbaOFYnLGLp2ZahSyo71cF4NDAb+glW62hwDbUlDdGaPn2OZ2CCslxIvAv8WabL/WNKiFC8Br4nwgm4t8SIVjOseIFeEyXH6yPYa126szqJzgDdEpCYmqpKY6CK+F2KZWE/gZSwT+1iEpF+sVinuA4pEuEe3lniRVMYVrUrsrf2yKkU+sBEoEGFLnGS017gGisjXsRKTaijFQOAiLBPLxurP9iLgT7bGa6WwA32Ba4BjgE+AHwGHA0+LJO/MJ8nWAXUJsCk6E4ENuAFreEm8TKtJLrzwQtLT09m923oBN3XqVJRSBAKBRsca02ofInwlwl0iHAWMByqBF4CvlWKGUhyd6LMpKMVApdgNlAErsbrFHI81KuRMrJ70u/Qp7HiSrcS1FeiFdUOLsX5xh4mwKY4yGmXo66+/zrnnnsucOXO46qqrOOyww7Db7Xz55ZdNnZ/QD1VnJGpUx2KVwi4CyrFKYS+KsFGntragFBnAWmAgTX9fqoFDRCiOq7A4kmwlrlrcwKFAFvChUlwcLYFpYfz48eTl5fHyyy/j9/vZvHkzkyZN0iUn5Yj2WQyIMBWrNHI5kI/13VilFH9QKnGGzUTfoE6AZruDfJjMpgXJZ1wNSzsZQG+shu7NSukZjJqZmcmECRP44IMPmDt3LkopLrnkEh1SUh4RIiJ8LNaURr2B67HMbIVSfKIU/xtd+alTI8LnwP1ARYNdJcDc+CuKL8lWVdwC9GliVwXwOXBGR/4SKaV6isi2pvYtWbKEk046CaUUJ554Ih999FGzl+kofYbmiVa/xmBVJ8/Bmv7nRaw20k65unlU8xqs7iC135sKID/ZO+cmW4mrKSqBmcCJHWFaSqk0pdQ4pdQrWO0OTXLiiSfSt29fRMRUEzshIgRFWCTCb4AewIPAaOArpVisFJcpRa5elfVposoowBvJblqQfCWuQqzXwwBBrEbYCSK8G/vPUocAvwGuAoqwZmd4QURK2nvpdp5viCFK4cR6UzcROA34N1YXizdEKNWprRaluA24FWv6oJ91xPe9s5FsxvUtVrtFJVbp5xwRdsTu+koBP8Ga/ngs8A+s9rOV8kNGdugga4M+lMILnI1lYqOBxVgm9rYIlRp1pQNbsV44OFJh5EDCGJdSKg2rY50P6I71gO8BAsB6EQkrxT6sLhB/Bv4oB5keWCk1HPj8YL3TlVJdsKZ7vgqrWD4H+JuINFX1NMaVAihFHnAelokdB7yFZWKLRaiOrxaVBn+fCN3Oh5PfwZp1dgsQEJGieGqJF53euJRSvR0OxzVKqclOpzPd5/OFBwwYkAWorVu3Vq9cuZK9e/fa0tPTnykv31EAnsdFeLYF1z0beA24SkTmNbFfYf2qXo1VVXgdy7A+lgNnmjGuFCM6s+4FWCY2BOt79RLw/sF+POtc4zhgZWt69yuljnK73dcHg8FfeL3eyLHHHhvu2rWrLRgMyoYNG8Lr1q1zOhyOzcXFxfdjNWMkz5RAuicEay6A9KysrGlOp7NiypQpVV9++aU0x5YtW+Tmm2+u8Xg8FS6X6xHAcZBrj8Bq/xKs0lrdfbnAdcCXwLro//Naob29aM97E20PkF4g14Msj878+gTIySDNzkgL0je6bNuTIOrgn0G2x+N5Nicnp3z69OnBrVu3SlOEQiFZuHChnHrqqaVOp3MvcKbu/IlZPusW0MyNyfV4PP7jjz++bPPmzU3elKbYtWuXnHPOOeVut3sT0LuZaw/AakyvXR+wHKv6OQprFZwi4DmsBSYO+iVqIuIyA6qJzh8gh4HcDLIaZBvITJCRDc0J5A/RtSbLQf7vQOYFHO10OndddtlllcXFxdJS3n//fenevXu52+2eC9h0502781a3gCZujNfj8XwxZcqU6oNNIdwcf/7zn4Mul2sHDaYyBroC27DaAGqNK4w1V/tG4Eagi+48MJF8AXIEyG0g60A2g9wLciyIAtkYLXFJ1Lwebcq8gCEOh6P4+eefj0gbKC4ulpEjR5a53e7n2vij3GlCu4CG4fF4/n7ppZdWRiJtujf7mT59etDr9X5KdNEIwIn1prGmjmnVRjWQrTvtJpI/okY1DORukE3RCNYxrlrzeqiueQFOl8u1bf78+e16MMrKymTIkCFlGRkZ1+rOi3blo24B9cTAWT169CgvKytr7f1oRCgUkmHDhpVlZGRcC6QDHzRjWoI1KPu3utNvIrUiamJzo6sfSRPm9eda83K5XP/385//vEJiwLp168ThcJQDh+nOgzbnnW4B+4WA8nq9X7/99tsiYhnPwIEDZfLkyY0yfurUqdKjRw8pLCyUSy+9VI4++mhJT0+Xs846q95xa9asEYfDUYw1T3lThlU3vk704rOJxAuQLU2YVm1UW9VGejkcjso9e/aISMuejQ8//FBOPvlk6datm2RmZkrv3r3ltttuk9qazPTp00Mej+c53elvc77pFvDDDeT4Qw89tLRuu9bTTz8tdrtddu78ob171qxZ4vF4JBAIyOeffy4zZsyQt99+W7p27Sq33357oxs5cuTIEqwFYG8EpgGPA68C/wG+wuoLVhNt6+qhOx9MpE6AHBUtWTU0rFB0exBEbLYbX7jqqqsqpQ4HezaWLl0qCxYskLVr18rmzZvlscceE0Cee+45ERH57rvvxG63VwK5uvOhTXmnW0Bt2O32h2+//fZ6rfHBYFD69+8vN998s4iIvPbaa2K322Xx4sV1D5Pdu3cLIG+++aY0ZP78+ZKbm3vQBTOBDN15YCK1IvrGsXah3D0gAZDnQW4FuQjkOJD87Ozsb1esWCF1aemzUUs4HBaHwyGzZs3av+2MM84oASbpzoc25Z1uAbWRl5e3atGiRY0yfM6cOZKdnS2LFi0Sp9Mp8+fPb3TMO++8I4Bs37690b7169eLx+PZpTt9Jkw0DJBskEEgWc0fgzsjI6OmqVXGW/JsiIh8//33csMNN8iAAQNk3759+7ffc889EafTOUt3PrQp73QLqA2Hw1HSlPFUV1dL7969JS0tTWbMmNHkjZkxY4b06NGjyX3hcFjS09ODgFN3Gk2YaG0Axw0YMKDJDlsHezYuuugicTqdAsjYsWOlYb+vhQsXSn5+/qe609iW6DTT2oRCoSy3291oe2ZmJqeddhoFBQVMmzatyXMDgQA+n6/JfWlpaWRmZoYBRyz1GgxxwuXxeJocBnSwZ+OBBx5g5cqVvPDCC6xcuZIHH3yw3n6v14uINH7oEoBOY1w2my1UVdX0NEKff/45I0aMaPZcv99PQUFBk/tEhGAwaIP4Dnw1GGJETXV181/dAz0bvXr14sgjj2TixImce+65LFy4sN7+6POWkMvfdRrjcjqdW5paPCIYDLJ27dpmS1Tff/89hYWFze7fsmULmZmZZSJSFlPBBkN8+Grz5s32aLWxHgd7Nuoe9+mnnzJo0KB629evX08oFGp28svOTKcxrurq6qXLly9vdHfWrl1LdXV1oxLV1q1bWb16Na+++ioA6enprF69mtLS+nO7rVixgqysrM86ULrB0GGIyO60tLTyb775ptG+pp6N5cuXc9ddd7FixQoKCwt57733GD9+PDt27OC2226rd/6yZcsqSkpKlnZ4IjoC3Y1sdRohzxg0aFBJw6E+8+bNk6ysLAkGg/W2jxs3rsmOpGvXrq133Pjx40uVUlfrTp8JE20Nj8fzzJ133hmSBjT1bHz00Udy/PHHS3Z2ttjtdunXr59cffXVsmnTpnrnVlRUiNvtriRBe893mvm4lFJpbrd7++LFiw8ZNWpUTK757bffMnjw4MqqqqpuYqqKhgRFKXVsly5dlm7fvt2ZkZERk2v+9a9/5frrr1+6b9++0TG5YJzpNFVFEYlUVlbeevnll5cHg8FYXI8rrriiIi0tbaYxLUMiIyKrampqVt5zzz0tmpTwYOzbt4/f//73lUVFRbfG4npa0F3kqxuA8ng8H02bNq1+vbANzJ07N+LxeL4GMnWny4SJ9gbQ2+FwlPn9fmkPkUhEJkyYUOl2u5/UnaZ25YduAY0EQQ+Xy7XzoYcealSnbykLFiyIOJ3OYuBHutNjwkSsIi0t7YKcnJyKhu24LSUSich1111X7fF4vgDcutPTntAuoElR0M/lcv330ksvrSwqKmrxjamqqpIbb7yx2uFw7AOO0Z0OEyZiHTab7RK3213x7LPPRlozZ93OnTvlzDPPLHe73WuBfN3paG9oF9CsMPC63e75+fn55TNnzgzXHWPVkPLycpk3b5706dOn1Ov1LgQO0a3fhImOCsDndrs3nXLKKWXvvfeeHMjAdu/eLffdd1/Y6/VWuFyuBwG7bv2xiE7zVrE5lFIjs7Ozb66urj5jyJAhVaNHj3b269cvQynF9u3bw0uXLi0LBAL2rKysj4uKiu4FFktnT5TB0E6UUllKqcs8Hs9Up9PZZfTo0WrUqFGuvLw8QqEQGzduDC1ZsqR81apVWZmZmW+UlJTcIyKrdeuOFZ3euGpRSuUBx6WlpQ13Op29gbTKysod4XB45dH/OvqUjLyM84Ah/gK/toU5DYZ4c+zyYwdGqiLPf37K53NcLtePMzMz80UkWFpauiEcDq8A/iMi3+vWGWsSxrgOhC/gewq4DPgYOM1f4K/QLMlg6HB8AV8u1hJ6hwA5/gJ/UwsUJyWdph9XO+kX/dcHvOsL+JwatRgMHY4v4LMD7wJdsFZXP1avoviSLMY1IPpvJtYNfN8X8Lk06jEYOgxfwGfDmn78R4AtGsO1ioozCW9cvoBPAT3qbLIDwzDmZUhCot/3OVgLFtujmzOif6cMCW9cWKbVcCiEHRgKfGjMy5Bk3A78Amud0LoceG6bJCMZjGsATU+GZseqNr4ZXzkGQ8fgC/iuwlqtqqk23K6+gM8TZ0naSBbjSm+wrQzLzD4CZsZdkcEQY3wBXxpwb/TPSBOHVADHxE2QZpLBuAYDLiyzqgaKgSeBbv4C/xh/gd+UuAwJj7/AHwEOBc7DWht0Jz+swg5WKSxlqosNSyqJSAWwDJgHvA5cDIxMpT4thtQg2rl6MbDYF/AVA3nACuB8YAwwUKO8uJIUHVDr4gv4umN1yuvhL/A3vfqGwZDg+AK+T4Bp/gL/+9G/0wDlL/CH9SqLD8lQVayHv8C/E1gNnKFZisHQIfgCvjysPlzLarf5C/yRVDEtSELjivIScKFuEQZDB3EqsMRf4E/ZJfeS1bheAX7qC/jMIrCGZGQcVltXypKUxuUv8O8CVgLjdWsxGGJJtOf86cAi3Vp0kpTGFWUBcJFuEQZDjDki+u8GrSo0k8zG9QpwhhnyY0gyTgcW+Qv8ydUdoJUkrXH5C/x7gE+An+rWYjDEkJRv34IkNq4oCzBvFw1Jgi/gywJGA+/p1qKbZDeuV4HTfQGfW7cQgyEGnAB84S/wJ91UzK0lqY0reoOXAT/TrcVgiAHjSPG3ibUktXFFMdVFQ7JwOqZ9C0gN43odODWV5ioyJB++gO8QrLUVlmuW0ilIeuPyF/j3AUuAs3RrMRjawVjgA3+Bv+FsvylJ0htXFFNdNCQ6pn2rDqliXK8DP/EFfNm6hRgMrSU6Zc1YTPvWflLCuKKTCn4InK1ZisHQFoYCJf4C/2bdQjoLKWFcUUx10ZComLeJDUgl43oTOMkX8OXoFmIwtBIzzKcBKWNc/gJ/CdZQiXM1SzEYWkx0koARwAe6tXQmUsa4opjqoiHROBnw+wv8pbqFdCZSzbjeAk6IztltMCQCKT9pYFOklHH5C/xlwL8w1UVD4mAa5psgpYwrillIw5AQ+AK+3kBXIKBbS2cjFY3rn8BIX8DXRbcQg+EgnA68G13F2lCHlDMuf4G/HKvN4DzdWgyGg2CG+TRDyhlXFPN20dCp8QV8Nqz1E/+lW0tnJFWNayFwnC/g66pbiMHQDMOB7f4C/zbdQjojKWlc/gJ/BVZb1/m6tRgMzWDeJh6AlDSuKGbdRUNnxrRvHYBUNq53gGN9AV933UIMhrpEp18ahjUBpqEJUta4/AX+KuBtTHXR0PkYA3zsL/BX6hbSWUlZ44piOqMaOiNmmM9BSHXjWgwM9QV8PXULMRgAfAGfwkxjc1BS2rj8Bf5qrHm6LtCtxWCIcjhgB9bqFtKZSWnjimI6oxo6E+OAxf4Cv+gW0pkxxmX1TP6RL+A7VLcQgwHTvtUiUt64/AX+GqxVgH6uW4shtfEFfBnAKcC7mqV0elLeuKKYzqiGzsAo4Ct/gX+3biGdHWNcFu8Bg3wBXx/dQgwpjRnm00KMcQH+An8QeBVTXTToxQzzaSHGuH7AvF00aCM6seUg4GPdWhIBY1w/8AFwuC/g66dbiCElOQ34d/RlkeEgGOOK4i/wh4BXgAm6tRhSEtO+1QqMcdXHVBcNcafOMB/TvtVCjHHV599AH1/Ad5huIYaU4kdADfC1biGJgjGuOkSri//AVBcN8cUM82klxrgaYzqjGuKNGebTSoxxNWYJ0NMX8A3ULcSQ/PgCPgdwAvC+bi2JhDGuBvgL/GHgZUx10RAfTgQ+9xf4i3QLSSSMcTWNmRnVEC9MN4g2YIyraZYBXX0B32DdQgxJj+kG0QaMcTWBv8AfwaoumlKXocPwBXw9gF7ASt1aEg1jXM1jOqMaOprTgfei3XAMrcAYV/N8DOT4Ar4f6RZiSFpM+1YbMcbVDNHq4t8xbxcNHYAv4EsDxmKMq00Y4zowC4CLomPJDIZYcgzwvb/Av0W3kETEGNeBWQ64gCG6hRiSDvM2sR0Y4zoA0bFjf8c00htij2nfagfGuA7OS8CFprpoiBW+gM8NDMeajcTQBoxxHZyVQBYwVLcQQ9JwCrDCX+Av0y0kUTHGdRCi1UXTp8sQS8ZhqontwhhXy1iAqS4aYoeZxqadGONqGQFAYb3CNhjaTHQxlhzgM71KEhtjXC3AVBcNMeR04F/RDs6GNmKMq+WYzqiGWGC6QcQAY1wt5zMgBPh0CzEkJr6ALx04FfiXbi2JjjGuFmKqi4YYcBzwrb/Av0O3kETHGFfrMJ1RDe3BDPOJEca4WsdaoBL4sW4hhoTEtG/FCGNcrcBUFw1txRfw5QJHAUt1a0kGjHG1ngXAhOh8SgZDSxkDLPMX+Kt0C0kGzMPXSvwF/i+AUmCEbi2GhMK0b8UQY1xtwyxfZmgx0Zc5pn0rhhjjaht/x1QXDS1nEJAOrNMtJFkwD14b8Bf41wHfA8fr1mJICE4HFkVf7hhigDGutmPeLhpaipnGJsYY42o7C4Cf+wI+m24hhs6LL+DLAk4C3tWtJZkwxtVG/AX+jcB3wIm6tRg6NaOA9f4C/17dQpIJY1ztw1QXDQfDdIPoAIxxtY+/AxeY6qLhAJhuEB2AMa524C/wfw1sw2rDMBjq4Qv4ugGHA5/o1pJsGONqP6YzqqE5TgM+9Bf4g7qFJBvGuNpPbXUxXbcQQ6fDtG91EErE9IlrL76A71Pgb1hVgpVmPvHUJTq852RgNPB74CR/gf9zvaqSD2Nc7cQX8J1kV/aF+Rn5GREiVUWhotLKSOWZ/gL/at3aDPHFF/Clu9Jcr3lsnpPH5Y5zbqjcoFaXrd5XJVUnRkdbGGKEqd60A1/A58lUmW/ff9j9zuOzj0dEMt7Z9477rm/vWuwL+A41bRspx8RDsw495ZnBz7gy0jIAeGnXS7lP7Hjir5jZRGKKaeNqH2cd4z5Gjs+2hiwqpRifN171zOyZhVVdMKQQXpv3F5d0u2S/aQGc3+V8VRmuLPAFfC6N0pIOY1ztw+mxeRrlocfmAXDGX45BJxGJlJSGS+u1vZRHylFKRbBWiDLECGNc7eOdpcVLbdurt+/f8FXlV6yrWJex/tL1ZUqpnyiljldK5WvUaOhAlFIZSqmhSqlTNj2/6T/zts+r2la9DYBgJMjM/86szlJZr/oL/NWapSYVpnG+nYxaNepam7L9eVzOuLTqULXtvb3v2TbdtSmc+1VumdvtprKykk2bNjnS09NLgNfKysoeFpEvdOs2tB2llEcpdUlOTs7/lJWVDTrkkEOqunbtGolEIpSOKHV0+U2XzF4ZvUK71e4waXxcHik/11/gL9atO6kQERPtCCC7i6/Lq72v6B28YPYFofcC70kwGJS6hMNhWbdunfzxj38M5uTkVHi93oXAIbq1m2j1vVbp6elX2+32svHjx5cuXLhQSktLpSGbtm+Su166S4b8ZEip2+3+FjhBt/ZkC+0CEjmAHzudzt2/+c1vKouLi6UlVFVVydSpU6sdDkcJME53Gky0+F5ne73ej4YMGVK2du3a5m5vI1599VXJycmpcLvd9xOt4ZiIwf3QLSBRAxjldDrLXnvtNWkLS5YsEY/HUw78THdaTBz0Xmd7PJ4vJ0+eXNWwNN0Sdu/eLcOGDStzu91PGfOK0T3RLSARAzjE4XAUv/3229Ieli9fLk6nsxwYrDtNJpq918rr9f7ryiuvrIpEIge6nQekuLhYjjzyyLLMzMz/1Z2mZAjtAhItAOXxeBbedNNN1RIDHnnkkbDH4/kMsOlOm4nGoZS69PDDDy+rrm7/7d6wYUPtD9UA3elK9NAuINECGN2zZ8+yqqoqEREJhUIycOBAmTx5sjRk6tSp0qNHDyksLBSXyyVAvbjpppskHA7LsGHDSoGJutNmotG9znQ4HEUrVqxo1b2uZdeuXXLyyScLIOvWrRMRkbvvvjuUnZ39pu60JXpoF5BokZ2d/cbDDz9cr87w9NNPi91ul507d+7fNmvWLPF4PBIIBOSzzz4TQD777DPZsWPH/igvLxcRkZdffllycnJW6U6bifoBXDRixIiS1tzrWoqLi+W6666TKVOmSG5urtRWM4uKisThcFQBPXSnL5FDu4BECsCRkZFRvW/fPqlLMBiU/v37y8033ywiIq+99prY7XZZvHixiIg8+eST0q9fP2mOmpoa8Xg8lXboIzBG3O4nxeV6QmCE7jSnTECepKffJNnZz0la2nUC2bm5uR8899xz9e7Vwe61iEhlZaVceuml8t1338nJJ58s48aNq3eNiy++uAL4nfY0J3BoF5BIAYwcOHBgk/0e5syZI9nZ2bJo0SJxOp0yf/78/fuuvPJK8Xg8csghh0jfvn3liiuukL1799Y7/yc/+UnxCofjTenZs0zuuy8id90Vlvz8crHbb9Sd7qQP6CdO5x6ZOLFCZs8WOeec8ojLtf0wh6O4btWvJfc6GAzKxIkTZdWqVRKJRCQ7O1tuu+22eufPmzdPcnJyXtGe7gQO7QISKYCrL774Yqt+14Dq6mrp3bu3pKWlyYwZM+rtmzt3rixatEjWrFkjzzzzjOTk5Mg555xT75iZkyeHazyeoNQ1tG+/FbHbqwS66k57UofX+3eZPj1U935Errkm+JTTGWrqTWJz9zoSicgvf/lLWbBggYiIfPPNNwLIwoUL653v9/slJyenUHu6EzjMtDatI7tbt24ZTe3IzMzktNNOY82aNUybNq3eviuvvHL//4866ii2bt3KHXfcQSgUIj3dugVH79mTFvnZzyLk5f1wYu/eMHy4jaVL30WpPR2QHgOA03kSv/51vQVP1OWXp5/x3HOilGp0eHP3+vrrr6dPnz5MmDABgFWrVqGUYsSI+jPa5OXlEQqF3B2RlFTBGFfrCAWDwWYHd37++eeMHDnyoBex2+0ApKX9MMa93OFANm6sP3OqCGzdGgSeAsz4xo4iPf1ItmzpQf/+P2zbsoW9SklPaOxcNL7X27dv55FHHsFms3HPPfcAP9RmunbtyiuvvMLZZ58NQCgUqp0xwtBWdBf5EimA80855ZQm27hqamokKytL/vKXvzS1u95xPp9PLrjggnrbr5w0qaLc4SiTJ56ISCgkUl0tcuedIXG7vxZI0532pI60tMlyxBFlUtuetWGDRHr3Lr/YZgtVVFQ0eQ8b3uvy8nJZs2ZNvRg9erSMGzdO1qxZI5WVlfuP/ec//yn5+fkrtac7gUO7gEQKoG92dnZFU+0egUBAAFm9evX+bVVVVXLTTTfJ0qVLZfPmzfLee+/JmDFjpGfPnrJly5Z65w8ePLjoRviFZGd/Li5XlTgcVZKd/bFAb93pTvoAJS7X7WK3l0uXLmXicJSI3X5DTk7ON5988kmL7nVTHHbYYfLHP/6x0fY777wzbLfbH9Ke7gQO7QISKQDldrt3NvVlnjdvnmRlZdWbGWL79u0yZswYyc/Pl6ysLBk0aJBcfvnlsnXr1nrnbt68Wex2ezngEBEEegh0053elAtwCPQXyBIRXC7XozfddFNNS+51QyoqKiQtLU2ef/75RvuOOOKIEmC89vQmcGgXkGiRnp7+h4kTJzb5ZrGtTJ06tcblcs3SnTYT9QMY5PV6K+pW89rLp59+Km63+zvMEK92hZkBtZWEQqF5b7zxRtjv98fkeoWFhcyaNStYXl7+UEwuaIgZIrIRWHr33XfHZNGTSCTCNddcU15dXT1DRMKxuGbKots5EzGUUpMOP/zwsqYabltDKBSSUaNGlWVmZt6iO00mmg7gUKfTWVp3OE9bmTlzZtjr9a7ClLbaf190C0jEwJoh4h9jx44trx1s3VpCoZBceumllR6PZ5n5InfuUEpNyMvLq9iwYcMB7uiBefnll8XpdBZhZoaISWgXkKgBZHg8njeHDx9etmnTJmkNO3fulHHjxlV4PJ7lgFt3WkwcPNLT0y/Pyckpb+0cbMFgUO6+++5g1LSO1Z2OZAntAhI5AFtWVtYtTqez4p577gk3HH/YkNLSUpk1a1bE6/VWOJ3OhwC77jSYaNX9/onL5fruwgsvrDjY9M3hcFgWLVokRx99dJnX6/0Y6KdbfzKFWeUnBiiljvR6vXfV1NT8dOzYsaGTTjrJffTRR1O7ys8XX3zBsmXLyt9++21bRkbGkuLi4ltFZIVu3YbWo5Ry2+32W9LS0q4ZPHiwOuOMM1zDhw+3de3alUgkQmFhIcuXL6954403akpLS78rLS39fyLyjJgHLaYY44ohSqmuwHiXyzXKbrcPFxGnUqq6pqbms9LS0mXAIhHZqlunof0opTKBsRkZGSO9Xu9oEemilIpEIpEtRUVF/xaRj4AVxrA6BmNcBoMh4TD9uAwGQ8JhjMtgMCQcxrgMBkPCYYzLYDAkHMa4DAZDwmGMy2AwJBzGuAwGQ8JhjMtgMCQcxrgMBkPCYYzLYDAkHMa4DAZDwmGMy2AwJBzGuAwGQ8JhjMtgMCQcxrgMBkPCYYzLYDAkHMa4DAZDwvH/AUz/9ohPFg4XAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Figure 1 (b)\n",
    "fig = draw_graph(pag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## References\n",
    "\n",
    "* Jaber, Amin, Jiji Zhang, and Elias Bareinboim. \"Causal identification under Markov equivalence: Completeness results.\"\n",
    "In International Conference on Machine Learning, pp. 2981-2989. PMLR, 2019.\n",
    "* Richardson, Thomas, and Peter Spirtes. \"Ancestral graph Markov models.\"\n",
    "The Annals of Statistics 30, no. 4 (2002): 962-1030.\n",
    "* Jaber, Amin, Jiji Zhang, and Elias Bareinboim.\n",
    "\"A graphical criterion for effect identification in equivalence classes of causal diagrams.\"\n",
    "In Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence. 2018.\n",
    "\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.8.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}