MarkovainToSemiMarkovianGraph.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyPdwQEzivS2yW8iv8TvYk0j",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/aliash98/BPM_CausalReasoning/blob/main/MarkovainToSemiMarkovianGraph.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "!pip install networkx==2.8"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "YMEjr3LiVCH1",
        "outputId": "4c7199f7-ec5d-4744-871e-e002f9a05d89"
      },
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/\n",
            "Collecting networkx==2.8\n",
            "  Downloading networkx-2.8-py3-none-any.whl (2.0 MB)\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m2.0/2.0 MB\u001b[0m \u001b[31m54.8 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hInstalling collected packages: networkx\n",
            "  Attempting uninstall: networkx\n",
            "    Found existing installation: networkx 3.0\n",
            "    Uninstalling networkx-3.0:\n",
            "      Successfully uninstalled networkx-3.0\n",
            "Successfully installed networkx-2.8\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import networkx as nx"
      ],
      "metadata": {
        "id": "K1aAw5i8sPli"
      },
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# A Markovian Graph 'G' as an input (with nodes being colored 'blue' if observed and 'grey' if unobserved)\n",
        "\n",
        "# Casual graph: Fig. 1 -> C\n",
        "\n",
        "# Let's assume V4 to be a grey variable such that V1 -> V4 -> D (we no longer have V1 -> D)\n",
        "\n",
        "G = nx.DiGraph()\n",
        "G.add_nodes_from([\n",
        "    ('V3', {\"color\": \"grey\"}),\n",
        "    ('V2', {\"color\": \"blue\"}),\n",
        "    ('V1', {\"color\": \"blue\"}),\n",
        "    ('D', {\"color\": \"blue\"}),\n",
        "    ('V4', {\"color\": \"grey\"}),\n",
        "    # ('V5', {\"color\": \"blue\"}),\n",
        "])\n",
        "G.add_edge('V3', 'V2')\n",
        "G.add_edge('V3', 'V1')\n",
        "G.add_edge('V2', 'V1')\n",
        "G.add_edge('V2', 'D')\n",
        "# G.add_edge('V1', 'D')\n",
        "\n",
        "G.add_edge('V1', 'V4')\n",
        "G.add_edge('V4', 'D')\n",
        "# G.add_edge('V3', 'V5')"
      ],
      "metadata": {
        "id": "xiYU4XN6VUqD"
      },
      "execution_count": 24,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Another random example to test\n",
        "\n",
        "G = nx.DiGraph()\n",
        "G.add_nodes_from([\n",
        "    ('A', {\"color\": \"grey\"}),\n",
        "    ('B', {\"color\": \"grey\"}),\n",
        "    ('C', {\"color\": \"blue\"}),\n",
        "    ('D', {\"color\": \"blue\"}),\n",
        "    ('E', {\"color\": \"blue\"}),\n",
        "    # ('V5', {\"color\": \"blue\"}),\n",
        "])\n",
        "G.add_edge('A', 'B')\n",
        "# G.add_edge('B', 'C')\n",
        "G.add_edge('B', 'D')\n",
        "G.add_edge('B', 'E')\n",
        "G.add_edge('C', 'E')\n",
        "G.add_edge('A', 'C')"
      ],
      "metadata": {
        "id": "u5eclK2BsHQW"
      },
      "execution_count": 32,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "G.nodes()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ZQx61mnqW1Tr",
        "outputId": "9c226330-bfae-4d2c-99b0-b6eecdb987f4"
      },
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "NodeView(('A', 'B', 'C', 'D', 'E'))"
            ]
          },
          "metadata": {},
          "execution_count": 4
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "G.edges()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "XHO0wijvXMbV",
        "outputId": "bb61e3f2-400f-4d5f-c960-a3f22b5539a0"
      },
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "OutEdgeView([('A', 'B'), ('B', 'C'), ('B', 'D'), ('B', 'E'), ('C', 'E')])"
            ]
          },
          "metadata": {},
          "execution_count": 5
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "nx.draw(G, with_labels=True, font_weight='bold', node_size=500)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 319
        },
        "id": "3U1Et98Tsxga",
        "outputId": "2c8151fa-0c1f-4e22-f59d-093e6366df20"
      },
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAAt+0lEQVR4nO3deViU9d4/8Pd9zwwwIAiKijoKQoAgDoJoqJALrqU+mdpjWXhkUWyzOj1Wp/Kxjj6a5nFJzRTx0vR3ojha2uKGoeGWiAiIKyKCWSaBKAzDbL8/TM8xFgEH7mHm/bour6Nxc/P2XNGbz718v4LJZDKBiIjIRohSByAiImpJLD4iIrIpLD4iIrIpLD4iIrIpLD4iIrIpLD4iIrIpLD4iIrIpLD4iIrIpLD4iIrIpLD4iIrIpLD4iIrIpLD4iIrIpLD4iIrIpLD4iIrIpLD4iIrIpLD4iIrIpLD4iIrIpLD4iIrIpLD4iIrIpcqkDNIbeYERxqQZavRH2chEqNyXkMnY3ERE1nMUXX2lFNZIzipCSWYzCkgooZCJEQYDRZILOYIRneydMClVhSr9ucHW0kzouERFZOMFkMpmkDlGbar0RK1LPIzG9AIIAVOmMdR7roBBhMgFxET0wO8oPdnJOgUREVDuLLL6rZRpMTTyKX8uroKmn8P5MqRDRycUBW+PC0dVV2YwJiYiotbK44rtapsG4Vem4WamDoQnRZIKAto4K7HwpguVHREQ1WNQ1wWq9EVMTjza59ADAYDLhZqUOzyUehc7Q8GmRiIhsg0UV34rU8/i1vKrJpXeXwWTCL+VarEi9YKZkRERkLSzmUmdpRTXCF6VCq699SiteEwND+fUa/7zz9JWw6+Rd6+fYy0UcezuKT3sSEdE9FvM6Q3JGEQThwccpffpB7tb53p9FR5c6jxWEO+ed+ZiPOSISEZEVsJjiS8ksrveVhbvaBI+Eo9+ABp2zSmdEyoliFh8REd1jEcWnNxhRWFLRoGNvn9qDqis59/7cbviMeo+/XFIBvcHIFV6IiAiAhRRfcakGCpkIncHwwGM1+cfv+/ODik8hE1FcqoGXuxMA4O4tTaEh11WJiMjqWETxafVGiA0sog5PvdPgS50AIAoCbmu02Lv3MJKTk7F9+3bMmDEDCxcubGpcIiJqxSyi+OzlIozN9HBppUaDviFqCLdvoLq6GjKZDJ06dWqWr0VERJbPIopP5aZs8Mvmf77H10Y9EnYdveo8XpAp4GSqQrleD5PJBIPBgLfeegufffYZAgMDERAQgMDAQAQGBsLb2xtyuUX8X0JERM3EIv4rL5eJ8GzvhIvXbz/w2D/f43Po3rve4vPu0AZfXSlEQkICtm3bBr1ej1OnTuHmzZs4c+YM8vLysGHDBuTl5eHnn3+Gr6/vvSK8W4q+vr6ws+O7gERE1sBiXmBfeyAfy1PPN+iVhoZyUIh4fbgfZvzxOsOXX36J5cuXIz09vdaHWyorK3Hu3Dnk5eXd+3XmzBkUFhbCy8vrXiHeLUV/f38olVwPlIioNbGY4iurrMajC+teuaUpzLVyi1arxYULF+4rxLy8POTn56Nr1673FWJgYCB69uyJNm3amOlvQURE5mQxxQcAS3afRdKhgkZtRVQX0aiD8UwqRnSqQq9eveDj44OwsDB07tz5wZ/cQDqdDpcuXapRiOfOnUOHDh1qFGJAQABcXV3N9vWJiKjxLKr4qvVGjFx+AEW/ax5qoWqZIKCtXIeTCyYCRgMUCgUEQcC4ceOQkpJixsS1MxgMuHz58n2XS+/+r4uLy30P1Nz95e7u3uy5iIjIwooPMN9+fDteHISnnxiOI0eOAADs7Oxw5swZeHvXvqB1SzAajSguLq4xIebl5cHOzq7GU6aBgYHw8PDgy/ZERGZkccUHPMwO7DJ4uNhjyx87sGdnZyM8PBwA0LFjR/j7+2PDhg1QqVTNFb1JTCYTfvnllxplePr0aRgMhhpPmQYGBqJbt24sRCKiJrDI4gPuXPZckXoeiekFEATU+7SnUiHCaALiI70xO8oXiv9YlzM6OhpXrlzB3r178eGHH2LlypVYunQpnnvuuVZRHL/99tt9l0vv/rp169Z9RXj3915eXpDJZFLHJiKyWBZbfHeVVVYjOaMIKSeKcbmkAgqZCFEQYDSZoDMY4dXeCZP7qvB0WLdan940GAwwmUz3XkzPyspCdHQ0fHx88Omnn6Jjx44t/Vcyi9LS0vvK8MyZMzh9+jRu3LgBf3//GvcQfXx8bPblfL3BiOJSDbR6I+zlIlRuSi5aTmTDLL74/pO5/gOm1Wrx/vvvY+PGjVi9ejWeeuqpZkgrjVu3buHs2bM1Lpv+/PPP8PHxqVGIvr6+sLe3lzq22ZVW/PEDU2YxCmv5gcmzvRMmhaowpV/tPzARkfVqVcVnbkeOHMG0adPQv39/fPzxx3Bzc5M6UrPRaDQ1Xs7Py8vD5cuX4eXlVeOhGn9/fzg6Okodu9Eac4ncQSHCZALiInpgdpQf7OScAolsgU0XH3BntZa3334b27Ztw/r16zF69GipI7Wouy/n//ke4sWLF9GlS5caD9UEBATA2dlZ6ti1avpDUSI6uThg6x8PRRGRdbP54rtr//79iImJwahRo/DRRx9Z7H/cW4per0d+fn6NB2vOnj0Ld3f3Wp80lXJiNtdrMDtfimD5EVk5Ft9/KC8vx+uvv479+/dj48aNGDx4sNSRLI7BYEBhYWGNl/Pz8vLQpk2bGivVBAYGokOHDs36BK05Fz7o3k6JPa8Nvu/JYCKyLiy+WnzzzTeYOXMm/vu//xsLFizgQtQNYDKZarycf/dJU5lMVuOhmsDAQHTu3NkshWjOpe6UChliI3rgjZH+D30uIrJMLL46lJSU4KWXXsLJkyexefNm9O/fX+pIrZLJZMKvv/5a46GaM2fOoKqqqtZC7NatG0SxYRNXaUU1whfVv7h5VXEeyo98Ce3VszDqqiB3bg+ldxjcomIhyBQ1jjfX4uZEZJlYfA/wxRdf4JVXXkF8fDzee+897stnRjdu3KjxUE1eXh5u3ryJgICAGk+a9ujRo8bL+Q/azqoi7wBu7FwKmIxQdOwB+86+0N+8jqorOej2ylaIDjV30XBQiHhtuB9m/rGdFRFZFxZfA/zyyy+Ij49HUVERNm/eDLVaLXUkq/afmwT/56/r16/Dz8/v3v3DsLAwLD2jxMXfKmo9j1FXhaurp8NYdQtOvYai/djXIAh3Jkld6TXIXdxrnfgAwLdjG+x9jfd4iawRi6+BTCYTNm3ahDlz5uD111/HG2+8YbMroUjl9u3bOHv27L1SvPbrdRzqMhE6Q+3/CmsKTuJ68nsAgC7xa6Fo3/A1WhUyAWfeH80VXoisEIuvka5cuYLp06ejsrISmzZtgp+fn9SRbNblGxV4/OMfUVltqPXjt0//gJKdSwEA3d/YBkHe8MvUjnYyfPdyJLzcncySlYgsB3+cbaTu3btj7969eO655zBo0CCsXLkSRqP5do2nhtPqjRDreSpU5uh67/f6m9cbdW5REOp9YIaIWi8WXxOIoogXX3wRhw8fxueff47hw4ejsLBQ6lg2x14uwljPBQv7rgH3Hl65eTgZJtO/i0x/8zpMBn2dn2s0mWDPJcyIrBK/sx+Cr68vfvzxR4waNQphYWHYsGEDeOW45ajclNAZ6p7KRDsHuI2YCQgiKk7/gGsbZ6Pk+49x/cv3cXXdDJh0VXV+rs5ghMqN728SWSMW30OSyWR48803sX//fqxatQrjxo3DtWvXpI5lE+QyEZ7t678H16bXUHR65v+g9AmDofw33M5Nha6kCM7BoyAo6t6Vwqu9Ex9sIbJSfCzRTHr37o1jx45h/vz56NOnD1asWIEpU6ZIHcvqTQpV1fseHwA4dA+CQ/egBp/TQSFict+GPwFKRK0Ln+psBsePH0d0dDTUajVWr14Nd3d3qSNZrbLKajy6sP6VWxqLK7cQWTdey2kG/fr1Q2ZmJlQqFdRqNXbu3Cl1JKvl6miHuIgeUCrM86+yUiFDfKQ3S4/IinHia2YHDx7EX/7yFwwePBjLly9H27ZtpY5kdbg7AxE1Br+7m9ljjz2G7OxsODg4QK1WIzU1VepIVsdOLmJrXDjaOioga+JuD3f349sSF87SI7JynPha0O7duxEXF4f/+q//wocffggnJ64KYk5N34FdBg8Xe2zhDuxENoE/2ragUaNGITs7G+Xl5ejTpw8OHz4sdSSr0tVViT2vDkbMoB6wl4tweMB9PzsRsJMJiI3ogT2vDWbpEdkITnwS2bZtG1588UVER0fj/fffh4ODg9SRrEpZZTWSM4qQcqIYl0sqoJCJEAUBRpMJOoMRXu2dUJD6TxQd/AK7d2zDkCFDpI5MRC2ExSeh69evIyEhAefPn8fmzZsRGhoqdSSrpDcYUVyqgVZvhL1chMpNCblMxNChQ5GWlgYHBwfMmzcPc+bMMcuO8ERk2XipU0IdO3bEv/71L7z11lsYPXo0PvjgA+h0OqljWR25TISXuxP8PZzh5f7vFVnu/sxXVVWFuXPn4i9/+YuEKYmopbD4JCYIAp577jlkZmbi8OHDGDhwIPLy8qSOZRO0Wu2937u7u2P8+PESpiGilsLisxAqlQrff/894uPjMXjwYCxduhQGQ+37zJF5ODs7w9/fH0OHDsWMGTMwceJEqSMRUQvgPT4LdOnSpXuX3TZu3AgfHx9pA1kpvV4PmUyGc+fOYfDgwbh06RJfMSGyAZz4LJC3tzfS0tIwYcIEhIeHY+3atdzuqBnI5XIIgoCePXsiIiICiYmJUkciohbAic/CnTlzBtHR0WjXrh02bNgAlYq7BjSHjIwMTJgwAfn5+bCz4zqdRNaME5+FCwgIwJEjRxAZGYnQ0FB89tlnnP6aQVhYGAICArBlyxapoxBRM+PE14qcPHkS0dHR8PX1xdq1a9GxY0epI1mVH374AQkJCcjLy4NMJpM6DhE1E058rUhISAgyMjLg5+eH4OBgbNu2TepIVmXIkCFwc3PD9u3bpY5CRM2IE18rdfjwYUybNg3h4eFYuXIl3NzcpI5kFb7++mt88MEHyMjI4CouRFaKE18rNXDgQGRlZcHV1RVqtRq7du2SOpJVGDduHKqqqrBnzx6poxBRM+HEZwVSU1MRExODMWPGYMmSJXB2dpY6Uqv22WefYcOGDUhLS5M6ChE1A058ViAqKgrZ2dmorq5GcHAwDh48KHWkVu2ZZ55BYWEhjhw5InUUImoGnPiszM6dOzFz5kxMmTIFCxYsgFLJPeaaYs2aNdi1axd27NghdRQiMjNOfFZm3LhxyMnJwdWrVxEaGorjx49LHalVmj59Oo4fP47c3FypoxCRmbH4rFD79u2RnJyMefPmYezYsXjvvfdQXV0tdaxWRalUYvbs2Vi0aJHUUYjIzHip08pdu3YN8fHxuHr1KjZv3ozevXtLHanVuHnzJry9vXH8+HF4e3tLHYeIzIQTn5Xr3Lkzdu7ciZdffhnDhg3DokWLuN1RA7Vt2xYJCQlYsmSJ1FGIyIw48dmQwsJCTJ8+HRqNBps2bYKfn5/UkSze9evX0bNnT+Tl5cHDw0PqOERkBpz4bIinpyf27duHZ599FgMHDsTHH38Mo9EodSyL1rFjR0ydOhXLli2TOgoRmQknPht1/vx5TJs2DY6OjkhKSoKnp6fUkSxWYWEhQkNDkZ+fD1dXV6njENFD4sRno/z8/JCeno6RI0ciLCwMSUlJ3O6oDp6enhg7dixWr14tdRQiMgNOfIScnBxER0dDpVJh3bp16Ny5s9SRLE5eXh6GDh2KgoICODo6Sh2HiB4CJz5C7969cezYMYSEhKBPnz5ITk6WOpLFCQwMxKBBg5CYmCh1FCJ6SJz46D4//fQTpk2bBrVajdWrV8Pd3V3qSBbj+PHjmDhxIi5evAg7OzsAgN5gRHGpBlq9EfZyESo3JeQy/jxJZMlYfFSDRqPBu+++i88//xxr167FuHHjpI5kMYYPH46Jz0ZD8BmElMxiFJZUQCETIQoCjCYTdAYjPNs7YVKoClP6dYOro53UkYnoT1h8VKcDBw5g+vTpGDJkCJYtW4a2bdtKHUlS1Xoj3vn8EL46Uw6ZTESVru5XQRwUIkwmIC6iB2ZH+cFOzimQyFLwu5HqNHjwYJw6dQoKhQJqtRqpqalSR5LM1TINRi4/gG/O34bOiHpLD7jzca3eiKRDBRi5/ACulmlaKCkRPQgnPmqQXbt2IS4uDhMmTMCiRYvg5OQkdaQWc7VMg3Gr0nGzUgdDE75dZIKAto4K7HwpAl1duU0UkdQ48VGDjB49Gjk5OSgtLUVISAgOHz4sdaQWUa03Ymri0SaXHgAYTCbcrNThucSj0Bm4Ug6R1Fh81GBubm7YsmULFi5ciKeeegpvvfUWtFqt1LGa1YrU8/i1vKrJpXeXwWTCL+VarEi9YKZkRNRULD5qtIkTJyI7Oxvnzp1DWFgYTp48KXWkZlFaUY3E9AJo6rifZzKZULwmBoWLxqJw0VjobhTVez6NzoD1P15CWSX3RiSSEouPmqRjx47Ytm0b5syZg1GjRuHvf/879Hq91LHMKjmjCIJQ98e1RadhKL9+78+3T+9/4DkF4c55iUg6LD5qMkEQ8PzzzyMzMxPp6ekYOHAgzpw5I3Uss0nJLK736c2K0z8AAOw6+fzx5wMPXO+0SmdEyoli84UkokZj8dFDU6lU2LVrF2JiYhAZGYl//OMfrX67I73BiMKSijo/btLrUHnuEADAbVgsRIc2MJRfh7Yo94HnvlxSAT0fciGSDIuPzEIQBCQkJODYsWPYvn07hg4dikuXLkkdq8mKSzVQ1LP0mCb/OIxVtyE6usK+exCUPv0A/HsKrI9CJqK4lO/1EUmFxUdm5ePjg7S0NIwfPx79+/fHp59+2iq3O9LqjRDrucF3+4+Cc3ykPwRBhKPfAABAxdlDMOl19Z5bFARo9Zz4iKTC4iOzk8lk+Otf/4qDBw9i/fr1GDNmDK5evSp1rEaxl4sw1lHYxqrb0ORnAABuZ+9B4aKx+G37/wEATNoKVF48Vu+5jSYT7LmEGZFk+N1HzSYwMBBHjhzBwIEDERISgi1btrSa6U/lpqzzZfOKs+mAQQfB3hFK3/B7v+RuXe58PLf+y506gxEqN67gQiQVudQByLopFArMnTsXY8eORXR0NLZv345PPvkEHTt2lDpaveQyEZ7tnXDx+u0aH7tbbM59RsNtaMy9f151JQe//r+3obl0AgZNOWRKl1rP7dXeiVsXEUmIa3VSi9FqtZg7dy42b96MNWvWYMKECVJHqtfaA/lYnnr+gQtSN4aDQsTrw/0w4zEfs52TiBqHxUct7vDhw5g2bRrCw8OxcuVKuLm5SR2pVmWV1Xh0YapZH0Sxl4s49nYU9+kjkhCvt1CLGzhwILKysuDq6gq1Wo3du3dLHalWro52iIvoAaXCPN8mSoUM8ZHeLD0iiXHiI0nt27cPsbGxGDNmDD766CO0adNG6kj3qdYbMXL5ART9rnmohaplgoDu7ZTY89rget8PJKLmx+9AktTw4cORnZ2N6upqBAcH48cff5Q60n3s5CK2xoWjraMCsvoW7qzH3f34tsSFs/SILAAnPrIYO3bsQEJCAp599lnMnz8fDg4OUke652qZBlMTj+LX8qo6d2uojVIhg4eLPbbEhXMTWiILwR8/yWKMHz8e2dnZuHLlCkJDQ5GRkSF1pHu6uiqx59XBiBnUA/ZyEQ4PuO+nVIiwl4uIjeiBPa8NZukRWRBOfGSRkpOT8corryAhIQHvvPMO7Ows54GQsspqJGcUIeVEMS6XVEAhEyEKAowmEyqrqtG5jQzTB/fE02Hd+CALkQVi8ZHFunbtGuLj4/Hzzz9j8+bNCAoKkjpSDXqDEcWlGmj1RtjLRWxc9RE0lRVYvHix1NGIqA681EkWq3Pnzti5cydefPFFDB06FIsXL4bBYJA61n3kMhFe7k7w93CGl7sTooYNxQ8/PHiHBiKSDic+ahUuX76M6dOnQ6vVYtOmTfD19ZU6Uq20Wi3c3d1RVFQEV1dXqeMQUS048VGr4OXlhdTUVEyZMgUDBgzAqlWrLHKzW3t7e4SHh+PgwYNSRyGiOrD4qNUQRRGvvPIKDh06hC1btmDkyJG4cuWK1LFqGDqUlzuJLBmLj1odf39/pKenIyoqCn379sXGjRstarujYcOGYf/+/VLHIKI68B4ftWrZ2dmIjo5G9+7dsW7dOnh4eEgdCTqdDu7u7sjPz4e7u7vUcYjoTzjxUaumVqvx008/Qa1WIzg4GF988YXUkaBQKBAREYG0tDSpoxBRLVh81OrZ2dlh/vz52LFjB+bOnYtnnnkGJSUlkmbifT4iy8XiI6vx6KOP4uTJk/Dw8IBarca3334rWZZhw4ax+IgsFO/xkVVKS0vD9OnTMWzYMCxbtgwuLi4t+vUNBgM6dOiA06dPo3Pnzi36tYmofpz4yCoNGTIE2dnZkMlkUKvVLT59yWQyPPbYY7zPR2SBWHxktZydnbFu3Tp88skneP755zF79mxUVla22Nfnaw1ElonFR1ZvzJgxyM7ORklJCUJCQnD06NEW+bp8wIXIMrH4yCa0a9cOW7ZswYIFC/Dkk0/ib3/7G7RabbN+zV69eqG8vNwiV5chsmUsPrIpkyZNwqlTp5CXl4f+/fvj1KlTzfa1RFHEkCFDOPURWRgWH9mcTp06Yfv27fjrX/+KESNGYMGCBdDr9c3ytfhaA5Hl4esMZNOKiooQGxuLsrIybN68GT179jTr+c+dO4cRI0agsLAQgiCY9dxE1DSc+MimdevWDbt378b06dMRGRmJ5cuXm3W7Iz8/P+j1ely6dMls5ySih8PiI5snCAJmzZqFo0eP4l//+heGDRuGgoICs5377msNJpPJIvcQJLI1LD6iP/j4+CAtLQ3jxo1D//79sX79+ofe7ig/Px+iKGLevHlwdXVFbGysmdISUVPxHh9RLU6fPo1p06ahQ4cOSExMRNeuXRt9jhs3bqBTp06Qy+Worq6GXC7H/Pnz8eabbzZDYiJqKE58RLXo1asXjhw5ggEDBiAkJARbt25t9PTn7u6OJUuWQC6XAwAcHBwwYMCA5ohLRI3AiY/oAU6cOIHo6GgEBATgk08+QYcOHRr1+QkJCUhMTAQA3Lx5E05OTs0Rk4gaiBMf0QP07dsXJ06cgLe3N9RqNb7++utGff7q1asREBAAR0dHlh6RBeDER9QIhw4dwrRp0zBo0CCsWLECrq6uDfq8W7duITMzE4MiIlFcqoFWb4S9XITKTQm5jD9/ErUkFh9RI1VUVGDOnDnYuXMnNmzYgBEjRtR7fGlFNZIzipCSWYzCkgooZCJEQYDRZILOYIRneydMClVhSr9ucHW0a6G/BZHtYvERNdHevXsRGxuLsWPHYvHixWjTps19H6/WG7Ei9TwS0wsgCECVru53+BwUIkwmIC6iB2ZH+cFOzimQqLnwu4uoiUaMGIHs7GxUVlaiT58+SE9Pv/exq2UajFx+AEmHCqDVG+stPeBOKWr1RiQdKsDI5QdwtUzT3PGJbBYnPiIz+PrrrzFr1ixMnToVs954B5PWH8fNSh0MTfj2kgkC2joqsPOlCHR1VTZDWiLbxuIjMpMbN27gpVdexQXfp1Gmkzep9O6SCQK6t1Niz2uDoeDDL0Rmxe8oIjNxd3dH32nvotKkeKjSAwCDyYRfyrVYkXrBTOmI6C5OfERmUlpRjfBFqdDqa7+fV7wmBoby63f+IIgQHV1g36Un3KLioHD1qPVz7OUijr0dxac9icyIEx+RmSRnFKEhW+4pffrBOfQJiPZO0Fw4ipLvVtZ5rCDcOS8RmQ+Lj8hMUjKLH/j0JgC0CR6JdiNmwm1oDABA/3txncdW6YxIOVH3x4mo8eRSByCyBnqDEYUlFQ069vapPagqzIbm8kkAgKPfwHqPv1xSAb3ByBVeiMyExUdkBsWlGihkInQGwwOP1eQf//cfZArYeTxS7/EKmYjiUg283LnOJ5E58EdIIjPQ6o0QG3KDD0CHp95B9zd3wiN6KWAyoeT7ldCV/VLn8aIg1PnADBE1HouPyAzs5SKMjXhAWhAE2Hk8AlFhD5iM0NdTfEaTCfZcwozIbHipk8gMVG5K6AwNm8ru3uOr/jUfRm0FBIU97Dp41Xm8zmCEyo0ruBCZC4uPyAzkMhGe7Z1w8frtBx579x6faO8Ee1Ug2kZMhczJtc7jvdo78cEWIjNi8RGZyaRQFZannq/zlQbVC0mNPqeDQsTkvqqHjUZE/4E/RhKZyZR+3WDudZAMBiOeDutm3pMS2TgWH5GZuDraIS6iB5QK83xb2YkmaE5+g/ffeRMVFQ17R5CIHozFR2RGs6P80MnFAbIGvtpQF5kgoKubEzI2vo/ffvsNarUa+/fvN1NKItvGRaqJzOxqmQbjVqWbdT++b7/9FgkJCRgzZgyWLFmCtm3bmjs2kc3gxEdkZl1dldj5UgS6tVM2+rKnUiFD93bKGpvQPvHEE8jNzYUgCAgKCsI333xj7thENoMTH1EzqdYbsSL1PBLTCyAIqHcBa6VChNEExEd6Y3aUb72bz+7fvx/x8fEIDw/HihUr4O7u3hzxiawWi4+omZVVViM5owgpJ4pxuaQCCpkIURBgNJmgMxjh1d4Jk/uq8HRYtwbvu1dRUYH33nsP//znP7FixQpMnjwZwkPeVySyFSw+ohakNxhRXKqBVm+EvVyEyk35UC+nHzlyBLGxsfD398fq1avRpUsXM6Ylsk68x0fUguQyEV7uTvD3cIaX+8OvyDJgwACcPHkSvXr1QnBwMJKSksCfZYnqx4mPyEpkZWUhJiYG7u7uWLduHby8vKSORGSROPERWYk+ffrg2LFjGDp0KMLCwrBq1SoYjdzOiOjPOPERWaGzZ88iNjYWoigiMTER/v7+Ukcishic+IisUM+ePXHw4EFMnjwZgwYNwocffgi9Xi91LCKLwImPyMoVFBQgPj4eZWVlSEpKglqtljoSkaQ48RFZuR49emDv3r2YNWsWoqKiMHfuXGi1WqljEUmGxUdkAwRBQGxsLLKysnDq1CmEhobi2LFjUscikgQvdRLZGJPJhOTkZLz66quYOnUq/v73v8PR0VHqWEQthhMfkY0RBAFTpkxBTk4Orl27BrVajbS0NKljEbUYTnxENm7Hjh144YUXMHbsWCxevBguLi5SRyJqVpz4iGzc+PHjkZubC4PBgKCgIHz33XdSRyJqVpz4iOieffv2IT4+HpGRkVi2bBnat28vdSQis+PER0T3DB8+HDk5OXBzc0NQUBBSUlKkjkRkdpz4iKhWhw4dQmxsLHr16oXVq1fDw8ND6khEZsGJj4hqNWjQIGRlZcHf3x9qtRqbNm3ilkdkFTjxEdEDZWZmIiYmBh4eHli3bh26d+8udSSiJuPER0QPFBoaiuPHjyMyMhKhoaFYs2YNtzyiVosTHxE1Sl5eHmJjY2FnZ4fExET4+vpKHYmoUTjxEVGjBAYGIj09HRMmTMCAAQOwZMkSbnlErQonPiJqskuXLiE+Ph7l5eVISkpC7969pY5E9ECc+Iioyby9vbFv3z7MmDEDw4YNw7x581BdXS11LKJ6sfiI6KEIgoD4+HicPHkSGRkZ6Nu3L44fPy51LKI6sfiIyCxUKhV27tyJt99+G2PHjsX//M//QKPRSB2LqAYWHxGZjSAIePbZZ5GTk4OioiKo1WocPHhQ6lhE9+HDLUTUbL766iu8+OKLePLJJ7Fo0SI4OztLHYmIEx8RNZ8nn3wSubm50Gg0CAoKwq5du6SORMSJj4haxp49ezBjxgwMGTIE//jHP9CuXTupI5GN4sRHRC1i5MiRyM3NhbOzM4KCgrBt2zapI5GN4sRHRC0uPT0dsbGxUKvVWLVqFTp16iR1JLIhnPiIqMVFREQgKysLPj4+UKvV+Oyzz7jlEbUYTnxEJKkTJ04gJiYGXbt2xaeffopu3bpJHYmsHCc+IpLU3ZVeBgwYgNDQUKxdu5ZbHlGz4sRHRBbj9OnTiImJgVKpRGJiIh555BGpI5EV4sRHRBajV69eOHz4MMaPH4/w8HAsXboUBoNB6lhkZTjxEZFFunjxIuLi4qDRaJCUlIRevXpJHYmsBCc+IrJIjzzyCPbv34+YmBgMGTIEH3zwAbc8IrPgxEdEFq+oqAgJCQkoLi5GUlIS+vbtK3UkasU48RGRxevWrRu++eYbvPHGG3j88cfx1ltvccsjajIWHxG1CoIg4Pnnn0d2djby8/PRp08fpKenSx2LWiFe6iSiVmnbtm146aWXMHHiRCxcuBBt2rSROhK1Epz4iKhVeuqpp5Cbm4tbt24hKCgIe/bskToStRKc+Iio1du1axdmzpyJqKgoLF26FG5ublJHIgvGiY+IWr3Ro0cjNzcXSqUSQUFB+Oqrr6SORBaMEx8RWZWDBw8iLi4OISEh+Pjjj9GxY0epI5GF4cRHRFblsccew6lTp+Dp6YnevXtj69at3PKI7sOJj4is1vHjxxETEwNPT0+sXbsWKpVK6khkATjxEZHV6tevH06cOIF+/fohJCQE69ev5/RHnPiIyDbk5OQgJiYGzs7OWL9+PXx8fKSORBLhxEdENqF37944cuQIHn/8cTz66KNYtmwZtzyyUZz4iMjmXLhwAbGxsdDpdEhKSkJAQIDUkagFceIjIpvj6+uLtLQ0PP/884iMjMSCBQug0+mkjkUthBMfEdm0wsJCzJw5E7/++iuSkpIQEhIidSRqZpz4iMimeXp64vvvv8err76KUaNG4W9/+xuqqqqkjkXNiMVHRDZPEARMmzYN2dnZOHfuHEJCQnD48GGpY1Ez4aVOIqI/SUlJwSuvvILJkydjwYIF3PLIynDiIyL6k0mTJiEnJwelpaVQq9XYt2+f1JHIjDjxERHV47vvvkNCQgJGjhyJjz76CK6urlJHoofEiY+IqB6PP/44cnNzoVAoEBQUhB07dkgdiR4SJz4iogZKS0tDXFwc+vXrh5UrV6JDhw5SR6Im4MRHRNRAQ4YMQXZ2Nrp27YrevXvj888/56LXrRAnPiKiJjh27BhiYmLwyCOP4JNPPkGXLl2kjkQNxImPiKgJHn30UWRmZiI4OBh9+vTBhg0bOP21Epz4iIge0qlTpxAbGws3NzesW7cOPXr0kDoS1YMTHxHRQwoODsbRo0cxfPjwew++GI1GqWNRHTjxERGZ0blz5xAXFwej0YgNGzagZ8+eUkeiP+HER0RkRv7+/jhw4ACeeeYZREREYOHChdzyyMJw4iMiaiaXL1/GjBkzcOPGDSQlJaFPnz5SRyJw4iMiajZeXl7YvXs3Xn75ZYwcORLvvvsutFqt1LFsHouPiKgZCYKA6dOn49SpUzh9+jRCQkJw9OhRqWPZNF7qJCJqISaTCV9++SVmz56NKVOmYP78+XBycpI6ls3hxEdE1EIEQcDTTz+NnJwc/Pbbb1Cr1fjhhx+kjmVzOPEREUnkm2++waxZs/D4449j8eLFaNu2rdSRbAInPiIiiYwdOxa5ubkAgKCgIHz77bcSJ7INnPiIiCzA/v37ER8fjwEDBmD58uVwd3eXOpLV4sRHRGQBhg0bhuzsbHTo0AG9e/fGF198wUWvmwknPiIiC3PkyBHExsbC398fa9asQefOnaWOZFU48RERWZgBAwbg5MmT6NWrF4KDg7Fx40ZOf2bEiY+IyIJlZWUhJiYGHTp0wLp16+Dp6Sl1pFaPEx8RkQXr06cPjh07hiFDhqBv375YvXo1tzx6SJz4iIhaibNnzyI2NhaiKGLDhg3w8/OTOlKrxImPiKiV6NmzJw4ePIjJkydj4MCBWLx4MfR6vdSxWh1OfERErVBBQQHi4+NRVlaGpKQkqNVqqSO1Gpz4iIhaoR49emDv3r2YNWsWhg8fjv/93//llkcNxOIjImqlBEFAbGwsTp48iaysLPTt2xc//fST1LEsHi91EhFZAZPJhOTkZLz66qt47rnn8MEHH8DR0VHqWBaJEx8RkRUQBAFTpkxBTk4Ofv75ZwQHB+PAgQNSx7JInPiIiKzQjh078MILL2DcuHH48MMP4eLiInUki8GJj4jICo0fPx65ubnQ6XQICgrC999/L3Uki8GJj4jIyu3btw/x8fGIjIzE8uXL0a5dO6kjSYoTHxGRlRs+fDhycnLg5uaGoKAgpKSkSB1JUpz4iIhsyKFDhxAbG4tevXph9erV8PDwkDpSi+PER0RkQwYNGoSsrCz4+/sjODgYmzdvvrflkUajwdWrVxt8Lr3BiMs3KnDul1u4fKMCekPrWDybEx8RkY3KzMxETEwMOnfujE8//RRz5szBgQMHUFBQAAcHh1o/p7SiGskZRUjJLEZhSQUUMhGiIMBoMkFnMMKzvRMmhaowpV83uDratfDfqGFYfERENkyn02Hx4sVYvHgxNBoN5HI55syZg3nz5t13XLXeiBWp55GYXgBBAKp0dU93DgoRJhMQF9EDs6P8YCe3rIuLLD4iIht369YteHp6orS0FADg4OCACxcuQKVSAQCulmkwNfEofi2vgqaewvszpUJEJxcHbI0LR1dXZbNkbwrLqmEiImpxW7duxc2bN+Hi4gI7OztUVVXhiSeeAHCn9MatSkfR75pGlR4AaHRGFP1+5/OvlmmaI3qTcOIjIrJxBoMBBQUFuHLlCoqKipCVlYXCwkJ8/kUKRi4/gKLfNTA8RFXIBAHd2ymx57XBUMikn7dYfEREVKslu88i6VBBoye92igVMsRG9MAbI/3NkOzhsPiIiKiG0opqhC9KhVZfs/SK18TAUH4dgABBYQ9R6QL7Ln5w6T8B9l3qLjZ7uYhjb0dJ/rSn9DMnERFZnOSMIghC/ccofcLg2DMCglyByrPp+GXLHFScTa/zeEG4c16psfiIiKiGlMziel9ZAIA2wSPh/sSr6BK3Bo4BjwFGA37ftRpGXVWtx1fpjEg5UdwccRuFxUdERPfRG4woLKlo8PGCKINrxDMAAGPVLWiLz9R57OUS6Vd4YfEREdF9iks1jX76UubS8d7vDZVldR6nkIkoLpX21QYWHxER3UerN0J80A2+P7nzsMsdMkfXOo8TBaHWB2ZaEouPiIjuYy8XYWzEA/8mowFl6f8EAIgOzrBXBdR5rNFkgr3ES5jJJf3qRERkcVRuSugacB/u9qk9qLxwDNqrZ6D//SogytBu9IsQFbUvcA0AOoMRKjdply9j8RER0X3kMhGe7Z1w8frteo/T5GdAUNhBVLaFY89IuPR/st73+ADAq70T5BKv3sLiIyKiGiaFqrA89XytrzSoXkhq0jkdFCIm91U9bLSHxnt8RERUw5R+3WDudb1MJuDpsG7mPWkTsPiIiKgGV0c7xEX0gFJhnppQKmSIj/SWfLkygMVHRER1mB3lh04uDpA18tWGP5MJAjxc7DE7ytdMyR4Oi4+IiGplJxexNS4cbR0VTS4/mSCgraMCW+LCLWJLIoC7MxAR0QM0fQd2GTxc7LHFwnZgZ/EREdEDVeuNWJF6HonpBRAE1LuAtVIhwmgC4iO9MTvK12ImvbtYfERE1GBlldVIzihCyoliXC6pgEImQhQEGE0m6AxGeLV3wuS+Kjwd1s0iHmSpDYuPiIiaRG8worhUA63eCHu5CJWbUvKX0xuCxUdERDbF8quZiIjIjFh8RERkU1h8RERkU1h8RERkU1h8RERkU1h8RERkU1h8RERkU1h8RERkU1h8RERkU1h8RERkU1h8RERkU1h8RERkU1h8RERkU1h8RERkU1h8RERkU1h8RERkU1h8RERkU1h8RERkU/4/cfc7MDpGxG0AAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Rule 1\n",
        "# Adding all observed nodes to the new graph\n",
        "\n",
        "G_sm = nx.DiGraph()\n",
        "for node in G.nodes(data=True):\n",
        "  # print(node)\n",
        "  if node[1]['color'] == 'blue':\n",
        "    G_sm.add_node(node[0])\n",
        "\n",
        "G_sm.nodes()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "o-ESKHosYDCj",
        "outputId": "c1eb1b36-a2ed-446c-88e0-65745e24cc19"
      },
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "NodeView(('C', 'D', 'E'))"
            ]
          },
          "metadata": {},
          "execution_count": 11
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Rule 2\n",
        "# Importing edges between two blue nodes\n",
        "for edge in G.edges(data=True):\n",
        "  from_node = edge[0]\n",
        "  to_node = edge[1]\n",
        "  from_node_color = G.nodes[from_node]['color']\n",
        "  to_node_color = G.nodes[to_node]['color']\n",
        "  if from_node_color =='blue' and to_node_color =='blue':\n",
        "    G_sm.add_edge(edge[0], edge[1])"
      ],
      "metadata": {
        "id": "tbE1i1wqanFg"
      },
      "execution_count": 12,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "G_sm.edges()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "QwJ0zyLwflOw",
        "outputId": "b839f445-b299-4e80-f342-b13db2fe628f"
      },
      "execution_count": 13,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "OutEdgeView([('C', 'E')])"
            ]
          },
          "metadata": {},
          "execution_count": 13
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Rule 3\n",
        "# Adding edges for directed paths through grey nodes between two blue nodes\n",
        "\n",
        "\n",
        "grey_nodes = [n for n,v in G.nodes(data=True) if v['color'] == 'grey']  \n",
        "blue_nodes = [n for n,v in G.nodes(data=True) if v['color'] == 'blue']\n",
        "\n",
        "for source in blue_nodes:\n",
        "  for sink in blue_nodes:\n",
        "    if source == sink or G.has_edge(source, sink):\n",
        "      continue\n",
        "    for path in nx.all_simple_paths(G, source=source, target=sink):\n",
        "      counter = 0\n",
        "      for between_node in path:\n",
        "        if between_node in blue_nodes:\n",
        "          counter += 1\n",
        "      if counter == 2:\n",
        "        G_sm.add_edge(source, sink)\n",
        "        break"
      ],
      "metadata": {
        "id": "Pdgs3SM5f13o"
      },
      "execution_count": 14,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "nx.draw(G_sm, with_labels=True, font_weight='bold', node_color='lightblue', node_size=500)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 319
        },
        "id": "9iVFZfWGQhPP",
        "outputId": "8b308717-d221-4081-ab8f-93936bee7ac0"
      },
      "execution_count": 15,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAARzUlEQVR4nO3de2yVdZ7H8c/vPOcc6Gmx9KKDOtwpA6y4iPGyLsENl+CIqBF2DW7MZkzMGPei4Q8RxZidFhkTZJ3Ixl13V//wj9GNzhpZo5GR9YKa4IUd4wIOraEUdRyopZSe0p4+z2//KCVeaGlPz+Fcvu/XXz306cOv/sHb7zm/53mc994LAAAjYoVeAAAA5xLhAwCYQvgAAKYQPgCAKYQPAGAK4QMAmEL4AACmED4AgCmEDwBgCuEDAJhC+AAApsQLvYBiEXmvdCZU6L0C55RKBIo5V+hlAQByzHT4esNIB4+ldeh4Wt19oZxzck7yfiCEVclAU85LafrElJIBwzEAlANn8ekMkffae7RLzR3dp14PfWzs1NA3q6ZS8+onMAUCQIkzF750JtSutnb19IcKR/GbB06qiAdaNLlOqUSQvwUCAPLKVPjSmVA7W48oE3pl80s7SYkgpiVT64kfAJQoMx9cRd5rV1t71tGTJC8pE0ba1dauyM7/LwBAWTETvr1Hu9TTH2YdvUFeUk9/qH1HT+RiWQCAc8zErs7eMFJzR/eQm1juWnKljnx5+Ad/vuW/Xtf0uZf84M9DLx3oOKGG2kp2ewJAiTERvoPH0iM67vK/WKZJU6adfn1ebd3w5+1Ma3Zt1ViWBgA4x0yE79Dx9LCXLAxaumatrlr20xGdM/JSa2cP4QOAElP24Yu8V3dfOKJj33jh1/q/3e+ffn3HA78Y9vjuvn5F3nNtHwCUkLIPXzozcEcWjWAX5kdv/vY7r88WPuec0plQVcmy/88IAGWj7P/FDr2Xc9JItnPet+0/RvxWpyQ5N3B+AEDpKPstiYFzIxn2suL9wPkBAKWj7Ce+VCLQSG9O8/3P+JauXqupP5k75PHee+7gAgAlpuzDF3NOlclAXSPY4PL9z/j+5Mo/GzZ8lck4G1sAoMSYuFfnZ+0ntK+9a0SXNIxUzEnz6idwOQMAlJiy/4xPkqZPTOXlvNOq83NeAED+mAhfMohpVk2lghy9Kxk4qaGmituVAUAJMvMv97z6CaqIBxpr+5wGnss3t563OAGgFJkJX8w5LZpcp0QQyzp+g8/jWzS5jk0tAFCiTGxu+TaewA4AtpmZ+AalEoHmJjPa/+7/KOYGdmcOJzh1zKyaKi2bfj7RA4ASV/bX8X3fe++9p+XLl+vkyZPq6cvoYGdarZ096u7rl3NOg7f19N6rMhnX1OoKTatOsZEFAMqEmbc6vffasmWLHn74YfX09CiRSKivr+/09yPvlc6ECr1X4JxSiYDP8QCgDJmZ+B577DGtX79+yNuXxZzjKQsAYICZ9+/uuOMOrVu3Ts45JZPJEd+/EwBQXsyMOLW1tbr88st1xRVXaPny5Xr11VcLvSQAQAGY+YwviiJdeuml2rJli6677rpCLwcAUCBm3up86aWXVFFRoRUrVhR6KQCAAjIRPu+9mpqa9NBDD8mxUxMATDMRvldeeUVhGGrVqlWFXgoAoMDKPnzeezU2Nmrjxo1MewCA8g/fjh071NXVpdWrVxd6KQCAIlD24WtqatKDDz6oWKzsf1UAwAiUdQ3eeustffXVV7r11lsLvRQAQJEo6/A1NjZqw4YNisfNXKcPADiLsg3f+++/r+bmZt1+++2FXgoAoIiUbfgaGxt1//33K5FIFHopAIAiUpa3LPvwww918803q6WlRePGjSv0cgAARaQsJ75NmzbpvvvuI3oAgB8ou4nvk08+0YoVK9TS0qJUKlXo5QAAikzZTXybNm3SunXriB4A4IzKauLbt2+frr32Wn3++eeqqqoq9HIAAEWorCa+Rx55RPfccw/RAwAMqWwmvpaWFl111VVqaWlRdXV1oZcDAChSZTPxbd68WXfffTfRAwAMqywmvtbWVi1cuFAHDhxQbW1toZcDAChiZTHxPfroo7rzzjuJHgDgrEp+4vviiy80f/587d+/XxdccEGhlwMAKHIlH757771XsVhMW7duLfRSAAAloKTD9/XXX2vu3Ln69NNPddFFFxV6OQCAElDSn/Ft3bpVa9euJXoAgBEr2Ymvvb1ds2fP1p49ezRlypRCLwcAUCJKduJ7/PHHdcsttxA9AMColOTEd+zYMc2aNUu7d+/WjBkzCr0cAEAJKcmJb9u2bVq5ciXRAwCMWslNfF1dXZo5c6befvttzZkzp9DLAQCUmJKb+J588kktWbKE6AEAslJSE186ndaMGTO0Y8cOzZ8/v9DLAQCUoJKa+J566ildc801RA8AkLWSmfhOnjypmTNnavv27Vq4cGGhlwMAKFElM/E988wzWrBgAdEDAIxJSUx8fX19amho0PPPP6+rr7660MsBAJSwkpj4nn32Wc2ePZvoAQDGrOgnvv7+fs2ZM0dPP/20Fi9eXOjlAABKXNFPfM8995wuvvhiogcAyIminvjCMNQll1yiJ554QsuWLSv0cgAAZaCoJ74XX3xR1dXVWrp0aaGXAgAoE0U78UVRpAULFmjz5s1auXJloZcDACgTRTvxvfzyy0okErr++usLvRQAQBkpyvB579XU1KSNGzfKOVfo5QAAykhRhu+1115Tb2+vbrrppkIvBQBQZooufN57NTY2auPGjYrFim55AIASV3Rl2blzp7755hutWbOm0EsBAJShogtfY2OjHnjgAQVBUOilAADKUFGF75133lFbW5tuu+22Qi8FAFCmiip8TU1N2rBhg+LxeKGXAgAoU0VzAfvu3bu1Zs0aNTc3K5lMFno5AIAyVTQTX2Njo9avX0/0AAB5dU4nvsh7pTOhQu8VOKdUIlDMOe3Zs0c33HCDWlpaNH78+HO1HACAQXn/MK03jHTwWFqHjqfV3RfKOSfnJO8HQliVDPTuux/qvg0PED0AQN7lbeKLvNfeo11q7ug+9XqYg32kWCymWTWVmlc/QTFuUwYAyJO8hC+dCbWrrV09/aHCUZw9cFJFPNCiyXVKJbiODwCQezkPXzoTamfrEWVCr2xO7CQlgpiWTK0nfgCAnMvprs7Ie+1qa886epLkJWXCSLva2hUVx5UWAIAyktPw7T3apZ7+MOvoDfKSevpD7Tt6IhfLAgDgtJzt6uwNIzV3dA+7iWX/x7v1m6e26bM9H6m3J63aH03SwsVL9DfrH1bie9fvhV460HFCDbWVSgZFc7khAKDE5awoB4+lh/3+rlde0kO3r9ZHb/5W9RdeqGtvWq0fTZ6i159/Vn0ne4Y+b+fw5wUAYDRyNvEdOp4ectrr7Unr337xoKIw1OIbV+vvf/mr08/a+8Ohg0qOrzjjz0Veau3s0ezaqlwtEwBgXE7CF3mv7r5wyO/v//gDnejskCStueue7zxgdtKUacOeu7uvX5H3XNsHAMiJnLzVmc4M3JFlKJ3ftJ/++vyLfzyqczvnlM4MHVUAAEYjJ+ELvddwA1l1bd3pr498cXhU53Zu4PwAAORCTsIXOKfh2vSTy65QVfVESdIL//IrRVF0+nt//OKw+jOZIX/W+4HzAwCQCzm5c0vkvV7+/R8UDXPM29t/oyfW/4OiKNK0OfM0a/5l6vjj1/rde2/r6Xd/p8rzqs/4czFJN86exGd8AICcyMnmlphzqkwG6hpmg8viVbeobtJFeunf/1m//9+PdLilWXWTLtSyv7xtyF2dklSZjBM9AEDO5OxenZ+1n9C+9q7hn8IwSjEnzaufwOUMAICcydkF7NMnpnJ1qu+YVp2f8wIAbMpZ+JLBwPP0ghy9Kxk4qaGmituVAQByKqdVmVc/QRXxQGNtn9PAc/nm1vMWJwAgt3IavphzWjS5TokglnX8Bp/Ht2hyHZtaAAA5xxPYAQCm5CV80sC1ffuOdulAR/ep10MfG7iBZ/A11FRpbn0Vkx4AIG/yFr5BfWGkg51ptXb2qLuvX845OTdwRxbvvSqTcU2trtC06hQbWQAAeZf38H1b5L3SmVCh9wqcUyoRMN0BAM6pcxo+AAAKjfcWAQCmED4AgCmEDwBgCuEDAJhC+AAAphA+AIAphA8AYArhAwCYQvgAAKYQPgCAKYQPAGAK4QMAmEL4AACmED4AgCmEDwBgCuEDAJhC+AAAphA+AIAphA8AYArhAwCYQvgAAKYQPgCAKYQPAGAK4QMAmEL4AACmED4AgCmEDwBgCuEDAJhC+AAAphA+AIAphA8AYArhAwCYQvgAAKYQPgCAKYQPAGAK4QMAmEL4AACmED4AgCmEDwBgCuEDAJhC+AAAphA+AIAphA8AYArhAwCYQvgAAKYQPgCAKYQPAGAK4QMAmEL4AACmED4AgCmEDwBgCuEDAJhC+AAAphA+AIAphA8AYArhAwCYQvgAAKYQPgCAKYQPAGAK4QMAmEL4AACmED4AgCmEDwBgCuEDAJhC+AAAphA+AIAphA8AYArhAwCYQvgAAKYQPgCAKYQPAGAK4QMAmEL4AACmED4AgCmEDwBgCuEDAJhC+AAAphA+AIAphA8AYArhAwCYQvgAAKYQPgCAKYQPAGAK4QMAmEL4AACmED4AgCmEDwBgCuEDAJhC+AAAphA+AIAphA8AYArhAwCYQvgAAKYQPgCAKYQPAGAK4QMAmEL4AACmED4AgCmEDwBgCuEDAJhC+AAAphA+AIAphA8AYArhAwCYQvgAAKYQPgCAKYQPAGAK4QMAmEL4AACmED4AgCmEDwBgCuEDAJhC+AAAphA+AIAphA8AYArhAwCYQvgAAKYQPgCAKYQPAGAK4QMAmEL4AACmED4AgCmEDwBgCuEDAJhC+AAAphA+AIAphA8AYArhAwCYQvgAAKYQPgCAKYQPAGAK4QMAmEL4AACmED4AgCmEDwBgCuEDAJhC+AAAphA+AIAphA8AYArhAwCYQvgAAKYQPgCAKYQPAGAK4QMAmEL4AACmED4AgCmEDwBgCuEDAJhC+AAAphA+AIAphA8AYArhAwCYEi/0AgAAtkXeK50JFXqvwDmlEoFizuXt7yN8AIBzrjeMdPBYWoeOp9XdF8o5J+ck7wdCWJUMNOW8lKZPTCkZ5PbNSee99zk9IwAAQ4i8196jXWru6D71euhjY6eGvlk1lZpXPyFnUyDhAwCcE+lMqF1t7erpDxWOojyBkyrigRZNrlMqEYx5HYQPAJB36Uyona1HlAm9somOk5QIYloytX7M8WNXJwAgryLvtautPevoSZKXlAkj7WprVzTGeY3wAQDyau/RLvX0h1lHb5CX1NMfat/RE2M6D7s6AQB50xtGau7oPuMmlruWXKkjXx6Wc07jKio0YWKtZl16mW782c81+08XnvF8oZcOdJxQQ21l1rs9mfgAAHlz8Fj6rMcsvHaprrlulRLJpN5/bbs2/vXNeu+17cOft/Ps5x0K4QMA5M2h4+lhL1mQpKVr1upvH/knPf7fb+rPr79JYX+//vXh+9Xbc+a4RV5q7ezJek2EDwCQF5H36u4LR3x8EI/rr/5unSTpRGeH9n/8wZDHdvf1Z73JhfABAPIinRm4I8tonH/Rj09/3flN+5DHOeeUzow8qt9G+AAAeRF6r9HebOXIl4dPf11dWzfkcc4NnD8bhA8AkBeBcxpNm8L+fv3ntq2SpKrqGs1ZeMWQx3o/cP5scDkDACAvUolAI7k52Bsv/Fof7nxd+z/+QF8e/FxBPK6f/+MvNa4iNeTPeO+zvoML4QMA5EXMOVUmA3WdZYPLx2+9oeT48Tqvpk7X/HSVbvzZXWq49LJhf6YyGc/6ptXcqxMAkDeftZ/Qvvaus17SMBoxJ82rn6DZtVXZ/XzulgIAwHdNnzj025VjMa06+/MSPgBA3iSDmGbVVCrI0QPVAyc11FSN6eG0hA8AkFfz6ieoIh5orO1zGngu39z67N7iHET4AAB5FXNOiybXKRHEso7f4PP4Fk2uG/OT2NncAgA4J3gCOwDAnMh77TvapQMd3adeD31s4AaewddQU6W59VVjnvQGET4AwDnXF0Y62JlWa2ePuvv65ZyTcwN3ZPHeqzIZ19TqCk2rTo1pI8uZED4AQEFF3iudCRV6r8A5pRJBzqa7MyF8AABT2NUJADCF8AEATCF8AABTCB8AwBTCBwAwhfABAEz5f0ZfCwRm53IkAAAAAElFTkSuQmCC\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Rule 4 \n",
        "# Extracting the hidden confounders\n",
        "confounding_list = []\n",
        "\n",
        "for source in grey_nodes:\n",
        "  potential_confounding = []\n",
        "  for sink in blue_nodes:\n",
        "    for path in nx.all_simple_paths(G, source=source, target=sink):\n",
        "      red_flag = 0\n",
        "      for between_node in path:\n",
        "        if between_node in blue_nodes and between_node != sink:\n",
        "          red_flag = 1\n",
        "          break\n",
        "      if red_flag == 0:\n",
        "        potential_confounding.append(sink)\n",
        "        break\n",
        "  if len(potential_confounding) > 1:\n",
        "    confounding_list.append(potential_confounding)\n",
        "\n",
        "print(confounding_list)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "x00JStKSVB3n",
        "outputId": "ebbeddc1-88a2-4f72-c786-d81d8e8e1387"
      },
      "execution_count": 18,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[['C', 'D', 'E'], ['D', 'E']]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Let's clean the output\n",
        "\n",
        "# we construct a graph\n",
        "G_confounded = nx.Graph()\n",
        "for item in confounding_list:\n",
        "  for i in range(len(item)):\n",
        "      for j in range(i + 1 , len(item)):\n",
        "        G_confounded.add_edge(item[i], item[j])\n",
        "\n",
        "nx.draw(G_confounded, with_labels=True, font_weight='bold', node_size=500)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 319
        },
        "id": "m0bWRaD8uJwj",
        "outputId": "e87fa62c-d956-487b-97c5-4d3728f26028"
      },
      "execution_count": 21,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAAm5ElEQVR4nO3deXhU9d028PvMZEgmQBaWsGUnCUkmJSWAAgapBAGr5LLWIBSwAgkKDHJpa/sq9SpPbZ9q7atQAoosQhOqLPpieRGt4vYCArJUH7IQQsguW0gC2Wc57x+YFDjZM5nfmTn35x9AJsPNdSE33zO/7zmSLMsyiIiINEInOgAREZEzsfiIiEhTWHxERKQpLD4iItIUFh8REWkKi4+IiDSFxUdERJrC4iMiIk1h8RERkaaw+IiISFNYfEREpCksPiIi0hQWHxERaQqLj4iINIXFR0REmsLiIyIiTWHxERGRprD4iIhIU1h8RESkKR6iA5DrstrsKK2sR6PVDk8PHQL9jfDQ899SRKRuLD7qksraJuw8UYI9p0pRVFELg14HnSTBLsuw2OwIGdgXjyYEYs74IPh59xEdl4hIQZJlWRYdgtSvyWrH2oN52HzoAiQJaLDY23ytl0EHWQZSE8OwMikKfTw4BRKRerD4qENlVfWYt/koLl1vQH07hXcno0GHIT5e2JE6ASP8jL2YkIio81h81K6yqnrMSj+E6joLbN34o6KXJPh6G7DPnMjyIyJV4DUoalOT1Y55m492u/QAwCbLqK6zYP7mo7DYOj8tEhH1FhYftWntwTxcut7Q7dJrZpNlXLzeiLUHzzkoGRFR9/FSJ7WqsrYJE14+iEZr61Na6YZFsF2/rPjvwxb+DX2GhLf6NZ4eOhx7PomnPYlIKK4zUKt2niiBJHX8OuPI8fDwH9byY523T5uvlaSb7/vkvSMdEZGIqFtYfNSqPadK211ZaNYvfjq8oyZ26j0bLHbsOVnK4iMioVh8pGC12VFUUdup19Z8+y80FP9Py48HTFvS7usLK2phtdl5hxciEobFRwqllfUw6HWw2Gwdvrb+/De3/bij4jPodSitrEfooL49ykhE1F0sPlJotNqh68wHfAAGP7Kq05c6AUAnSW0emCEicgZebyIFTw8d7L102Ncuy/DkLcyISCBOfKQQ6G/s9LL5nZ/x9Rs9HX0CQtt8vcVmR6A/7+BCROKw+EjBQ69DyMC+yL9c0+Fr7/yMzyv4R+0WX+jAvjzYQkRCsfioVY8mBGLNwbw2VxoCl23t8nt6GXRIGRvY02hERD3Cf3pTq+aMD4KjP+aTZWD2uCDHvikRURex+KhVft59kJoYBqPBQX9EbE3wLv4aNdeUtzkjInImFh+1aWVSFIb4eEHfydWGtuglCaEBfnh4pAEJCQl45513HJSQiKjreJNqapejn8d34sQJLFiwAPHx8diwYQMGDBjQC6mJiNrGiY/aNcLPiH3mRAQNMHb5sqfRoEfwAONtD6EdN24cTp06hSFDhmD06NH4+OOPeyM2EVGbOPFRpzRZ7Vh7MA/rD+bCYDCgvftXGw062GUgbXI4ViZFwtDG+sKnn36KRYsWITk5GX/5y1/g7e3dS+mJiP6DEx91Sh8PHSb1q4Bh/2r86v5RiAzoB4NegncfPfp5esC7jx4GvYTIgH54ZloUjj2fhF9PH9Vm6QHAtGnT8O2336KyshJjxozB8ePHnfg7IiKt4sRHnTZ37lxMmDABK1euBHDzKQ6llfVotNrh6aFDoL+x28vpu3btwooVK7B06VKsWrUKBoPBkdGJiFqw+KhTvv/+e8TGxuLChQvw8/PrlV+jvLwcixYtQkVFBTIyMhAdHd0rvw4RaRsvdVKnbNy4EXPmzOm10gOA4cOH48CBA1i4cCESExOxbt062O18kgMRORYnPupQU1MTQkJC8MknnyAuLs4pv2ZeXh4WLFgAHx8fvP322wgM5K3OiMgxOPFRh95//33ExMQ4rfQAICoqCocPH8a9997LpXcicihOfNShxMREPPvss3jkkUeE/PpceiciR+LER+06ffo0ioqKkJycLCwDl96JyJE48VG7UlNTER4ejhdeeEF0FABceieinuPER22qqKjAe++9h9TUVNFRWnDpnYh6isVHbdq6dStmzZqFgIAA0VFu4+/vjx07duCll17CrFmzsHr1algsFtGxiMhF8FIntcpmsyEiIgK7du3C+PHjRcdpE5feiairOPFRqz788EMEBASouvSA/yy9L1q0iEvvRNQpnPioVTNmzMD8+fOxYMEC0VE6jUvvRNQZnPhI4ezZs/j3v/+NlJQU0VG6hEvvRNQZnPhI4emnn0b//v3xpz/9SXSUbuPSOxG1hRMf3ebGjRvIzMzEU089JTpKj3DpnYjawuKj22RkZGDq1KkICgoSHaXHjEYj1q5di23btiEtLQ1msxl1dXWiYxGRYCw+aiHLMtLT02E2m0VHcSguvRPRrVh81OLzzz+HTqfDlClTREdxOC69E1EzFh+1WLduHcxmMyRJEh2l18yePRunT5/G0aNHMWnSJOTm5oqOREROxuIjAEBRURG++uorzJ8/X3SUXseldyJt4zoDAQCef/551NfXY82aNaKjOBWX3om0hxMfoaGhAVu2bMHy5ctFR3E6Lr0TaQ8nPsL27dvx7rvv4sCBA6KjCMWldyJt4MSncbIstxxq0TouvRNpAyc+jTt69CjmzZuHvLw86PV60XFUg096J3JfnPg0Lj09HcuWLWPp3YFL70TuixOfhl26dAnR0dEoKCiAv7+/6DiqtWvXLqxYsQJLly7FqlWrYDAYREcioh7gxKdhmzZtQkpKCkuvA1x6J3IvLD6NslgsePPNNzW5wtAdXHonch8sPo364IMPEB4ejvj4eNFRXIYkSVi6dCmOHDmCzMxMzJgxA6WlpaJjEVEXsfg0iisM3celdyLXxsMtGvTdd9/hgQceQGFhIQ9q9BCX3olcDyc+DVq/fj2eeuoplp4DcOmdyPVw4tOYyspKhIeHIycnB0OHDhUdx61w6Z3INXDi05ht27bhpz/9KUuvF3Dpncg1cOLTELvdjqioKGRkZGDixImi47g1Lr0TqRcnPg356KOP4OvriwkTJoiO4va49E6kXiw+DUlPT4fZbIYkSaKjaAKX3onUiZc6NSI/Px+TJk1CUVERjEaj6Dia0/ykd19fX7z99tsYMWKE6EhEmsWJTyM2bNiARYsWsfQEaV56nzx5MsaMGcOldyKBOPFpQG1tLYKDg3Hy5EmEhoaKjqN5XHonEosTnwZkZmZi8uTJLD2V4NI7kVic+NycLMsYPXo0Xn/9dUybNk10HLoDl96JnI8Tn5v76quvYLVakZSUJDoKtWLatGn47rvvUFVVxaV3Iidh8bk5rjCon5+fHzIzM/HSSy9h1qxZWL16NSwWi+hYRG6LlzrdWGlpKUaPHo2ioiL0799fdBzqhPLycixatAgVFRXIyMhAdHS06EhEbocTnxvbuHEj5s2bx9JzIVx6J+p9nPjcVGNjI0JCQvDFF19wanBR586dw4IFC+Dj48OldyIH4sTnpnbv3o0f/ehHLD0XFhkZiUOHDnHpncjBOPG5qQkTJuCFF15AcnKy6CjkAFx6J3IcTnxu6JtvvsHFixfx4IMPio5CDsKldyLH4cTnhp544gnExsbiN7/5jego1Au49E7UM5z43MyVK1fwwQcfYPHixaKjUC/h0jtRz7D43MzmzZvxs5/9DAMHDhQdhXoRl96Juo+XOt2I1WpFeHg49u7di4SEBNFxyEm49E7UNZz43Mi+ffsQFBTE0tMYLr0TdQ0nPjeSlJSE1NRUzJ07V3QUEoRL70Qd48TnJrKzs5GdnY2f//znoqOQQFx6J+oYJz43sWzZMgwePBj/9V//JToKqQSX3olax4nPDVRXV+Odd97Bk08+KToKqQiX3olax+JzA9u3b8f06dMxfPhw0VFIZYxGI9auXYtt27YhLS0NZrMZdXV1omMRCcXic3F2ux3r16/HihUrREchFePSO9F/sPhc3Keffgqj0Yh77rlHdBRSOS69E93E4nNx6enpMJvNkCRJdBRyEbNnz8bp06dx9OhRTJo0Cbm5uaIjETkVi8+FFRQU4MiRI/jFL34hOgq5GC69k5ZxncGFPffcc5BlGX/9619FRyEXxqV30hpOfC6qrq4Ob7/9NpYtWyY6Crm45qX3e++9l0vvpAmc+FzUli1bsHfvXuzbt090FHIjJ0+exPz587n0Tm6NE58LkmW55VALkSONHTuWS+/k9jjxuaDDhw9j4cKFyM3NhU7Hf7tQ7+CT3sld8W9NF7Ru3TosX76cpUe9ikvv5K448bmY8vJymEwmFBYWwtfXV3Qc0ohdu3ZhxYoVWLp0KVatWgWDwSA6ElG3cWRwMW+99Rbmzp3L0iOn4tI7uRMWnwtpamrCxo0bsXz5ctFRSIO49E7ugsXnQt5//33ExMTAZDKJjkIaJUkSli5diq+//ho7duzAzJkzUVZWJjoWUZew+FwIVxhILbj0Tq6Mh1tcxOnTp5GcnIwLFy7Aw8NDdByiFlx6J1fDic9FpKenY+nSpSw9Uh0uvZOr4cTnAioqKhAREYG8vDwMHjxYdByiNnHpnVwBJz4XsHXrViQnJ7P0SPVuXXpPSEjAN998IzoSkQInPpWz2WyIiIjArl27MH78eNFxiDpt9+7dMJvNXHon1eHEp3L79+9HQEAAS49cTkpKCpfeSZVYfCrHFQZyZVx6JzXipU4Vy83NxZQpU1BcXAxPT0/RcYh6hE96J7XgxKdiGzZsQFpaGkuP3AKX3kktOPGp1I0bNxASEoJvv/0WQUFBouMQORSX3kkkTnwqlZGRgalTp7L0yC01L70PHTqUS+/kdJz4VEiWZZhMJmzYsAE/+clPRMch6lVceidn48SnQp999hl0Oh2mTJkiOgpRr+PSOzkbi0+FmlcYJEkSHYXIKfz8/JCZmYmXXnoJDz30EFavXg2LxSI6FrkpXupUmaKiIiQkJKCoqAj9+vUTHYfI6crLy7Fo0SJUVFQgIyMD0dHRoiORm+HEpzJvvvkmHn/8cZYeaRaX3qm3ceJTkYaGBgQHB+Pw4cOIjIwUHYdIOC69U2/gxKci7777LsaOHcvSI/oBl96pN3DiUwlZljFu3Dj84Q9/wIMPPig6DpHqcOmdHIUTn0ocO3YMVVVVeOCBB0RHIVIlLr2To3DiU4n58+cjISEBzz77rOgoRKrHpXfqCU58KnDp0iXs378fCxcuFB2FyCVw6Z16gsWnAps2bUJKSgr8/f1FRyFyGVx6p+7ipU7BLBYLwsLCsH//fsTHx4uOQ+SSysvLsXjxYly9epVL79QhTnyC7d27F+Hh4Sw9oh4YPnw4PvzwQyxevJhL79QhTnyCTZkyBWazGSkpKaKjELkFLr1TRzjxCfTdd98hPz8fDz/8sOgoRG6DS+/UEU58Aj355JMIDAzEiy++KDoKkVvi0ju1hhOfIJWVldi1axfS0tJERyFyW1x6p9Zw4hPktddew8mTJ7Fjxw7RUYg0gUvv1IwTnwB2ux3r16+H2WwWHYVIM7j0Ts1YfAJ89NFH8Pf3x4QJE0RHIdIULr0TwOITIj09HWazGZIkiY5CpEkpKSk4ffo0jh07hkmTJiE3N1d0JHIiFp+T5efn48SJE3jsscdERyHSNC69axcPtzjZM888A09PT7z88suioxDRD7j0ri2c+JyopqYGf//73/HUU0+JjkJEt+DSu7Zw4nOijRs34sCBA9i7d6/oKETUBi69uz9OfE4iyzLS09OxYsUK0VGIqB1cend/LD4n+eqrr2C1WjF16lTRUYioA0ajEWvWrMH27duRlpYGs9mMuro60bHIQVh8TsIVBiLXk5SUxKV3N8TP+JygpKQE8fHxKCoqQv/+/UXHIaJu2L17N8xmM5YuXYpVq1bBYDCIjkTdxInPCTZu3Ih58+ax9IhcGJfe3QeLr5c1NjZi06ZNWL58uegoRNRDty69T548mUvvLorF18t2796N+Ph4REdHi45CRA4gSRKeeuopHDlyBDt27MDMmTNRVlYmOhZ1AYuvlzUfaiEi98Kld9fFwy296JtvvkFKSgrOnz8PvV4vOg4R9RIuvbsWTny9KD09HcuWLWPpEbk5Lr27Fk58veTKlSuIiopCfn4+Bg4cKDoOETnJwYMHsXDhQj7pXcU48fWSzZs345FHHmHpEWkMl97VjxNfL7BarQgPD8fevXuRkJAgOg4RCcKld3XixNcL9u3bh6CgIJYekcZx6V2dWHy9gCsMRNTszqX39PR0Lr0LxkudDpaVlYVp06ahqKgIffr0ER2HiFSET3pXB058DrZ+/XosWbKEpUdEClx6VwdOfA5UXV2N0NBQZGVlYfjw4aLjEJGKceldHE58DrR9+3bMmDGDpUdEHeLSuzic+BzEbrcjJiYGW7ZsQWJioug4RORCuPTuXJz4HOSTTz6B0WjEPffcIzoKEbmY5qX36upqLr07AYvPQZpXGCRJEh2FiFyQn58fMjIy8NJLL+Ghhx7C6tWrYbFYRMdyS7zU6QAFBQW46667UFxczEsURNRj5eXlWLx4Ma5evYqMjAw+z9PBOPE5wBtvvIGFCxey9IjIIbj03rs48fVQXV0dgoODcfz4cYSHh4uOQ0RuhkvvjseJr4feeecdTJw4kaVHRL3i1qX3hIQELr07ACe+HpBlGWPGjMErr7yCGTNmiI5DRG6OS++OwYmvBw4fPoy6ujrcf//9oqMQkQZw6d0xOPH1wJw5czBx4kSsXLlSdBQi0hguvXcfJ75uKi8vx8cff4wnnnhCdBQi0iAuvXcfi6+b3nrrLcydOxe+vr6ioxCRRnHpvXt4qbMbmpqaEBISgk8//RQmk0l0HCIiLr13ASe+bnjvvfcQExPD0iMi1eDSe+dx4uuGe+65B7/61a/wyCOPiI5CRKTApff2ceLrolOnTqGkpATJycmioxARtYpL7+3jxNdFixcvRkREBJ5//nnRUYiIOsSldyVOfF1QUVGB999/H6mpqaKjEBF1CpfelTjxdcGrr76KM2fOYPv27aKjEBF1GZfeb+LE10k2mw0bNmyA2WwWHYWIqFu49H4Ti6+T9u/fj4CAAIwfP150FCKibnPk0rvVZkfh1VqcvXgDhVdrYbW5xvoEL3V20vTp0/H4449j/vz5oqMQETlEd5beK2ubsPNECfacKkVRRS0Meh10kgS7LMNisyNkYF88mhCIOeOD4Ofdxwm/i65j8XVCbm4upkyZguLiYnh6eoqOQ0TkMLIsY+PGjXjxxRfx+9//HsuWLYNOp7wY2GS1Y+3BPGw+dAGSBDRY2p7uvAw6yDKQmhiGlUlR6OOhrouLLL5OePrpp+Hj44M//vGPoqMQEfWK9pbey6rqMW/zUVy63oD6dgrvTkaDDkN8vLAjdQJG+Bl7I3a3qKuGVejGjRvIzMzEk08+KToKEVGvaWvpvayqHrPSD6HkWn2XSg8A6i12lFy7+fVlVfW9EbtbOPF1YP369fj888+xZ88e0VGIiJzi5MmTWLBgAeJG/xgXxyxGWXUjbD2oCr0kIXiAEf96ZgoMevHzlvgEKibLMtLT07nCQESaMnbsWJw8eRKVgRNRdKW6R6UHADZZxsXrjVh78JyDEvaMh+gAavbZZ5/Bw8MDU6ZMER2FiMipGux6FBsjIFnbvrzZUJqN61/vRmNZLuyWBnj0Hwhj+Dj4Jy2GpDfc9tp6iw2b/l8BUhPDhJ/25MTXjuZpT5Ik0VGIiJxq54kStPdXX232l7i043+h/vw30PsMQr+4++DhNxQ3/n0AsqWx1a+RpJvvKxo/42tDUVEREhISUFRUhH79+omOQ0TkVNNe/xL5l2ta/Tm7pQFl6xfC3nADfU33YeBDz0CSbs5Rlsrv4eEzSDHxNYsM6IdPnhF7FY0TXxveeOMNPP744yw9ItIcq82OooraNn++sTQH9oYbAADfSY+1lB4AGPyHtVl6AFBYIf4OL/yMrxX19fXYsmULjhw5IjoKEZHTlVbWw6DXwWKztfrztrqqlu97+AZ06b0Neh1KK+sROqhvTyL2CCe+VuzcuRPjxo1DZGSk6ChERE7XaLVD184HfHpvv5bvW6svd+m9dZKExnYOzDgDi+8Osixj3bp1WLFihegoREROJcsySktLcfzrw2hsbGrzdZ4jYqDzuvkxUPWRnZDl/xSZtfoyZJu1za+1yzI8Bd/CjJc673Ds2DFUVVVh5syZoqMQEfUKWZZx6dIlZGVlISsrC2fOnGn5vpeXF2JNcbCOXwm0MfTp+njB//4nUfF/X0dt1udoulIIz2FRsNVcQ33haQStyISkb/18hMVmR6C/2NuXsfjukJ6ejuXLl7d6k1YiIldz9erVVgsOAEwmE+Li4vDjH/8Y8+bNg8lkwqBBgwC0f6oTAPqZ7oNH/8G4fmwPGstyUVNRAo/+g9A/fgYkQ9s38w8d2Bcegu/ewnWGW1y8eBExMTEoKCiAv7+/6DhERJ1WVVXVasE1NDTAZDK1lFzz94cMGdLujvKbX57HmoN57T6Foau8DDo8Oy0KS+4d6bD37A5OfLfYtGkTUlJSWHpEpFo3btxAdna2ouSqq6sRGxvbUnAPPvggTCYTRowY0a2bcMwZH4TXP81zaHZZBmaPC3Loe3YHJ74fWCwWhIaG4sCBAxg9erToOESkcXV1dcjNzb1tejtz5gyuXLmC6OhoxRQXHBzs8I9oXv04F1sPX+jyUxlaYzTosTgxDL+ePsoByXqGE98P9u7di4iICJYeETlVY2Mjzp49qyi4srIyREZGthRbWloaTCYTwsLCoNfrnZJtZVIU9v/P9yi5Vt/jpzMM9fHEyiR1rIhx4vvBlClTYDabkZKSIjoKEbkhi8WCc+fOKT6DKywsRFhY2G2fv8XFxSEiIgIeHuJnk+bn8VXXWbpVfnpJgq+3AfvMiap5GC2LD8B3332HBx54AIWFhTAY2r7VDhFRR2w2G86fP68ouPz8fAQFBSkKLioqCn36iH1aQUe6/wR2PYb6eCJTZU9gZ/EBWLJkCYKCgvDiiy+KjkJELsJut6OwsFBxyOTs2bMYOnSo4hRldHQ0jEb1/OXfVU1WO9YezMPmQxcgSWj3tKfRoINdBtImh2NlUqQqHj57K80XX2VlJcLDw5GTk4OhQ4eKjkNEKtN8N5PmYmv+NicnBwMGDFAcMomJiXHrm9tX1TVh54kS7DlZisKKWhj0OugkCXZZhsVmR+jAvkgZG4jZ44KEP3evLZovvtdeew2nTp1CZmam6ChEJJAsy7h48aLikEl2djb69u2rKLjY2Fj4+vqKji2U1WZHaWU9Gq12eHroEOhvFL6c3hmaLj673Y7IyEjs2LEDEyZMEB2HiJzkypUrt01vzd/q9XrFZ3AmkwkDBgwQHZkcSPyRIYE++ugj+Pv74+677xYdhYh6QWVlpeKQyZkzZ2CxWG4ruNmzZ8NkMiEgoGuP2CHXpOniW7duHcxmc7fuakBE6nH9+vWWu5ncOsXV1NQgNja2peSSk5NhMpkwbNgw/n+vYZq91Hnu3DlMmjQJxcXFLn3SikhLamtrkZOTo5jirl69ipiYGMVlyqCgIBYcKWi2+J555hl4enri5ZdfFh2FiO7Q0NCA3NxcxarA999/j6ioKMWqQFhYGJ+oQp2myeKrqalBSEgITp06hZCQENFxiDSrqakJeXl5ioIrLi7GyJEjFQU3cuRIVdzNhFybJv8E7dixA/feey9Lj8hJrFYrzp8/rzhkUlBQgJCQkJZimzNnDkwmEyIjI1V/NxNyXZqb+GRZxujRo7FmzRokJSWJjkPkVux2Oy5cuKAouHPnzmHYsGEt01vzt6NGjYKXl5fo2KQxmiu+L774AkuXLkV2djY/9CbqJlmWUVxcrDhkkpOTg0GDBikOmcTExMDb21t0bCIAGiy+Rx99FPfddx+WL18uOgqR6smyjPLycsVncNnZ2ejfv/9t01tcXBxiY2PRv39/0bGJ2qWp4ispKUF8fDyKior4PyfRLWRZxuXLlxUFl5WVhT59+igOmZhMJvj7+4uOTdQtmiq+3/3ud7h+/Tr+9re/iY5CJExFRUWrBWe32xW36jKZTBg8eLDoyEQOpZnia2xsRHBwML788ktER0eLjkPU66qrq1tK7daSq6ura7Xghg4dys+9SRM0s86we/duxMfHs/TI7dTU1CAnJ0dxw+XKykrExsa2FNvMmTNhMpkQGBjIgiNN08zEd/fdd2PVqlVITk4WHYWoW+rr65Gbm6tYFbh06RKio6MVU1xISAjvZkLUCk0U3/Hjx/HYY48hPz8fer1edByidjU1NeHs2bOKGy6XlpYiIiJCsQsXHh7OP9dEXaCJ4vvlL3+JuLg4PPfcc6KjELWwWCzIz89XHDK5cOECQkNDFbtwERERMBgMomMTuTy3L74rV64gKioK+fn5GDhwoOg4pEE2mw0FBQWKQyb5+fkYMWKE4pDJqFGj4OnpKTo2kdty++L785//jPz8fGzZskV0FHJzdrsdxcXFikMmZ8+eRUBAgKLgoqOjeTcTIgHcuvisVivCw8Oxd+9eJCQkiI5DbkKWZZSVlSkOmeTk5MDPz09xyCQ2Nhb9+vUTHZuIfuDW6wz//Oc/ERQUxNKjbpFlGZcuXVIUXHZ2Nry8vFqKbeLEiUhNTUVsbCz8/PxExyaiDrj1xDd16lSkpaVh7ty5oqOQyl29elVxyOTMmTOQJElxyMRkMvHzYiIX5rbFl5WVhfvvvx+FhYV8rhe1qKqqUhwyOXPmDBobGxU3XDaZTAgICOCyN5GbcdtLnevXr8eSJUtYehp148YNZGdnK6a46upqxMbGthTbgw8+iLi4OAwfPpwFR6QRbjnxVVdXIzQ0FFlZWRg+fLjoONSL6urqkJOToyi4K1eutNzN5NYpLigoiHczIdI4t5z4tm3bhhkzZrD03EhjYyNyc3MVlynLy8sRGRnZUmxLliyByWRCWFgY72ZCRK1yu4nPbrcjOjoaW7duRWJioug41EUWiwXnzp1THDIpKipCeHi44pBJREQEPDzc8t9vRNRL3K74Pv74Y/z2t7/F6dOn+ZmNitlsNpw/f15RcOfPn0dwcLCi4KKiovh5LRE5hNsV36xZs/Dwww9j8eLFoqMQbk7ghYWFilOUeXl5GDp0qOIUZXR0NLy8vETHJiI35lbFV1BQgLvuugvFxcW8FZSTybKMkpISxSGTnJwcDBgwQLELFxMTg759+4qOTUQa5FbF1/z0hVdffVVwEvclyzK+//57xSGT7Oxs9O3bV3E/ytjYWPj6+oqOTUTUwm2Kr66uDsHBwTh+/DjCw8NFx3ELV65cUdxwOSsrCx4eHoqCM5lMGDBggOjIREQdcpvjcP/4xz8wceJEll43XLt2raXUbi05q9V62+XJxx57rOVuJkRErsotJj5ZljFmzBi88sormDFjhug4qnX9+nVkZ2crpriamhrFKUqTyYRhw4bxZCwRuR23mPgOHz6M+vp63H///aKjqEJtbS1ycnIUqwLXrl1DTExMS7FNnz4dJpMJQUFBLDgi0gy3mPjmzJmDSZMm4emnnxYdxakaGhpa7mZy6xR38eJFREVFKVYFQkNDebsuItI8lyo+q82O0sp6NFrt8PTQIdDfiMuXLsJkMqGwsNBtTw82NTUhLy9PccikuLgYI0eOVKwKhIeH824mRERtUH3xVdY2YeeJEuw5VYqiiloY9DroJAl2WYbFZoe3rRZD64uw87+fhp+3a9/Zw2q1Ij8/X3HIpKCgACEhIYobLkdGRsJgMIiOTUTkUlRbfE1WO9YezMPmQxcgSUCDxd7ma/voJUiShNTEMKxMikIfD3VfzrPZbLhw4YKi4M6dO4fhw4crDpmMGjWKdzMhInIQVRZfWVU95m0+ikvXG1DfTuHdyWjQYYiPF3akTsAIP2MvJuwcWZZRXFysOGSSm5uLwYMHK05SxsTE8I4zRES9THXFV1ZVj1nph1BdZ4GtG9H0kgRfbwP2mROdVn6yLKO8vFxRcNnZ2fD19VUUXGxsLPr37++UbEREdDtVFV+T1Y7pa75EybX6bpVeM70kIXiAEf96ZgoMesdd9pRlGZcvX1YcMjlz5gw8PT0VpyhjY2Ph7+/vsF+fiIh6TlXF9+rHudh6+EKXLm+2xWjQY3FiGH49fVS3vr6ioqLVgpNlWXGK0mQyYdCgQT3OTEREvU81xVdZ24QJLx9Eo1VZeqUbFsF2/TIACZLBEzqjDzyHR8Hnrp/Bc3jbxebpocOx55PaPe1ZXV2tOGSSlZWFuro6xSGTuLg4DBkyhMveREQuTDXLXjtPlKCjPjGOHAedty8ay3JQl3sIdXlfY1Dyc+gb3fqT1iXp5vs+ee9I1NTUIDs7W1FwlZWViI2NbSm4Bx54AHFxcRgxYgQLjojIDamm+PacKm13ZQEA+sVPh3fURMh2G67u+9+oy/kK1z5af7MQDcrj/g0WO177P0fw8i+n4fLlyxg1alTL5HbffffBZDIhJCSEdzMhItIQVRSf1WZHUUVtp18v6fTwS5yLupyvYG+4gcbSHBjDxrT+3kZ/HPjoY0RGjIRer3dUZCIiclGqGHVKK+u7fPpS7/OfR+PY6qrafJ2nwQNeA0ew9IiICIBKiq/Raoeui5+n3TzscpPe26/N1+kkqdUDM0REpE2qKD5PDx3sXThcKtttqDr0DgBA59UfnoExbb7WLsvwVPktzIiIyHlU8RlfoL8RFlvHU1nNt/9C3bljaCzLgfVaGaDTY8DM5a0ebGlmsdkR6C/+9mVERKQOqig+D70OIQP7Iv9yTbuvqz9/ApKhD3RGX3hHT4bPXQ+3u8cHAKED+8LDgXdvISIi16aK4gOARxMCseZgXqsrDYHLtnbrPb0MOqSMDexpNCIiciOqGYXmjA+Co+8hI8vA7HFBjn1TIiJyaaopPj/vPkhNDIPR4JhIRoMeaZPDXf7htERE5FiqKT4AWJkUhSE+XtD38FZheknCUB9PrEyKdFAyIiJyF6oqvj4eOuxInQBfb0O3y6/5eXyZqRMc+kgiIiJyD6p5OsOtuv8Edj2G+ngiUyVPYCciIvVRZfEBNx9Ku/ZgHjYfugBJQrs3sDYadLDLQNrkcKxMiuSkR0REbVJt8TWrqmvCzhMl2HOyFIUVtTDoddBJEuyyDIvNjtCBfZEyNhCzxwXxIAsREXVI9cV3K6vNjtLKejRa7fD00CHQ38jldCIi6hKXKj4iIqKe4rhERESawuIjIiJNYfEREZGmsPiIiEhTWHxERKQpLD4iItIUFh8REWkKi4+IiDSFxUdERJrC4iMiIk1h8RERkaaw+IiISFNYfEREpCksPiIi0hQWHxERaQqLj4iINIXFR0REmsLiIyIiTfn/q0B/SRJjktAAAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "list(G_confounded.edges())"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "pjNb52LexK-8",
        "outputId": "eade303a-a05e-4297-cfc8-b0bac36acda5"
      },
      "execution_count": 23,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[('C', 'D'), ('C', 'E'), ('D', 'E')]"
            ]
          },
          "metadata": {},
          "execution_count": 23
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Let's set up a function\n",
        "\n",
        "# INPUT: A Markovian graph (with nodes being colored 'blue' if observed and 'grey' if unobserved\n",
        "# OUTPUT: A Graph containing only solid arcs and a set of tuples(V_i, V_j) that are connected with dashed arcs\n",
        "\n",
        "def markov_to_semimarkov(G):\n",
        "\n",
        "  G_sm = nx.DiGraph()\n",
        "\n",
        "  # Rule 1 ------- Adding all observed nodes to the new graph\n",
        "  for node in G.nodes(data=True):\n",
        "    if node[1]['color'] == 'blue':\n",
        "      G_sm.add_node(node[0])\n",
        "  \n",
        "  # Rule 2 ------- Importing edges between two blue nodes\n",
        "  for edge in G.edges(data=True):\n",
        "    from_node = edge[0]\n",
        "    to_node = edge[1]\n",
        "    from_node_color = G.nodes[from_node]['color']\n",
        "    to_node_color = G.nodes[to_node]['color']\n",
        "    if from_node_color =='blue' and to_node_color =='blue':\n",
        "      G_sm.add_edge(edge[0], edge[1])\n",
        "\n",
        "  # Rule 3 ------- Adding edges for directed paths through grey nodes between two blue nodes\n",
        "  grey_nodes = [n for n,v in G.nodes(data=True) if v['color'] == 'grey']  \n",
        "  blue_nodes = [n for n,v in G.nodes(data=True) if v['color'] == 'blue']\n",
        "\n",
        "  for source in blue_nodes:\n",
        "    for sink in blue_nodes:\n",
        "      if source == sink or G.has_edge(source, sink):\n",
        "        continue\n",
        "      for path in nx.all_simple_paths(G, source=source, target=sink):\n",
        "        counter = 0\n",
        "        for between_node in path:\n",
        "          if between_node in blue_nodes:\n",
        "            counter += 1\n",
        "        if counter == 2:\n",
        "          G_sm.add_edge(source, sink)\n",
        "          break\n",
        "    \n",
        "  # Rule 4 ------- Extracting the hidden confounders\n",
        "  confounding_list = []\n",
        "\n",
        "  for source in grey_nodes:\n",
        "    potential_confounding = []\n",
        "    for sink in blue_nodes:\n",
        "      for path in nx.all_simple_paths(G, source=source, target=sink):\n",
        "        red_flag = 0\n",
        "        for between_node in path:\n",
        "          if between_node in blue_nodes and between_node != sink:\n",
        "            red_flag = 1\n",
        "            break\n",
        "        if red_flag == 0:\n",
        "          potential_confounding.append(sink)\n",
        "          break\n",
        "    if len(potential_confounding) > 1:\n",
        "      confounding_list.append(potential_confounding)\n",
        "    \n",
        "  # we construct a graph for confounding arcs\n",
        "  G_confounded = nx.Graph()\n",
        "  for item in confounding_list:\n",
        "    for i in range(len(item)):\n",
        "        for j in range(i + 1 , len(item)):\n",
        "          G_confounded.add_edge(item[i], item[j])\n",
        "  # nx.draw(G_confounded, with_labels=True, font_weight='bold', node_size=500)\n",
        "\n",
        "  return G_sm, list(G_confounded.edges())\n"
      ],
      "metadata": {
        "id": "i0Xa9WyixUQJ"
      },
      "execution_count": 30,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Let's test the function\n",
        "\n",
        "final_graph, confounded_vars = markov_to_semimarkov(G)\n",
        "nx.draw(final_graph, with_labels=True, font_weight='bold', node_size=500)\n",
        "print(confounded_vars)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 337
        },
        "id": "R0YrkeEazvb6",
        "outputId": "97d63b5f-8efd-4656-c510-b9951c7e8bbb"
      },
      "execution_count": 33,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[('C', 'D'), ('C', 'E'), ('D', 'E')]\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAAPSUlEQVR4nO3dW4xcd33A8d+Zy87O+rb2Qi5kTWIjGyJokZKYupUDAreBiiLRNgmpgisVnFairgxSWxX1IapUqbyg4IiiPlh5SlRF9VtARVGtqpJpkzRJi1BJMKBcvBbEwvEFr/c2lz5Yaxl7r7PrzFn/Pp+X2d2cGR9npfn6P+ec3ym63W43ACCJSr93AADeScIHQCrCB0AqwgdAKsIHQCrCB0AqwgdAKrV+7wCLa7U7MXZmIqZanWjUKjG6uRm1qn+zAPRC+ErqzPh0PP3iiTjy8li8cXo86tVKVIoiOt1uzLQ7cfvIurj/rtF4aNfWGB4a6PfuAqwZhckt5TLd6sSho8fj8LHXoigiJmc68247WK9Etxuxf8+2OLh3ZwzUrAIBFiN8JXLy7EQ8fPi5eOv8ZEwsELyrNeuVuHnjYDy1f3fcNty8jnsIsPYJX0mcPDsRn/nmsTh3cSbaPfxKqkURm4bq8cyBPeIHsACfjZXAdKsTDx9+rufoRUS0u904d3EmPn/4uZhpL321CJCN8JXAoaPH463zkz1Hb1a7242fn5+KQ0d/vEp7BnDj8VFnn50Zn47dXzsaU625V2lj3/pCtM+fiogiinojKs2N0XjPztj4kd+PxnveP+dzGrVKPP/Vvc72BJiDFV+fPf3iiSiKxbdrvu+eGPrAnihq9bj46rH4+ZN/HeOvHptz26K49LoAXEv4+uzIy2MLXrIwa/2H74t3ffrL8Z7934qhOz8a0WnH29/9x+jMTF6z7eRMJ468NHY9dhdgzRO+Pmq1O/HG6fFlPaeoVGN4zx9FRERn8pcxNfbKnNu9fno8Wk5yAbiG8PXR2JmJqPcweqy68abLX7cvnp1zm3q1EmNnJnrdNYAblvD10VSrE5WlHOC7yqWTXS6pDg3PuU2lKOY9YQYgM+Hro0atEp1lnlTb7bTj7LF/joiIyuCGaIzeOed2nW43GkaYAVzDkOo+Gt3cXPLF5he+/2xc/PHzMXXylWi9fTKiUo0tn/rzqNQH59x+pt2J0c0muABcTfj6qFatxO0j6+Inpy4suu3ET1+Moj4QleamGPrAvbHxI5+d9zq+iIg7Rta5dRHAHISvz+6/azS+cfT4vJc0jH7piWW/5mC9Eg/cPbrSXQO4IVkS9NlDu7bGas/O6XYjHrxn6+q+KMANQvj6bHhoIPbv2RbN+ur8Kpr1ajxy73bjygDmIXwlcHDvzrh542BUe7i04UrVoohbNjbi4N4dq7RnADce4SuBgVolntq/OzYN1XuO3+z9+J7cv7uni+IBsnB3hhLp/Q7s1bhlYyOedAd2gEUJX8lMtzpx6OjxOHzstSiKWHCAdbNeiU434pF7t8fBvTus9ACWQPhK6uzF6Xj6xRNx5KWxeP30eNSrlagURXS63Zhpd+KOkXXxwN2j8eA9W53IArAMwrcGtNqdGDszEVOtTjRqlRjd3HRxOkCPhA+AVCwbAEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASEX4AEhF+ABIRfgASKXW7x0AYG1qtTsxdmYiplqdaNQqMbq5GbVq+ddTwgfAkp0Zn46nXzwRR14eizdOj0e9WolKUUSn242ZdiduH1kX9981Gg/t2hrDQwP93t05Fd1ut9vvnQCg3KZbnTh09HgcPvZaFEXE5Exn3m0H65XodiP279kWB/fujIFauVaBwgfAgk6enYiHDz8Xb52fjIkFgne1Zr0SN28cjKf2747bhpvXcQ+XR/gAmNfJsxPxmW8ei3MXZ6LdQy6qRRGbhurxzIE9pYlfudafAJTGdKsTDx9+rufoRUS0u904d3EmPn/4uZhpL321eD0JHwBzOnT0eLx1frLn6M1qd7vx8/NTcejoj1dpz1bGR50AXOPM+HTs/trRmGrNvUob+9YXon3+1DU/v/VPHo+Bm7fP+ZxGrRLPf3Vv38/2dDkDANd4+sUTURSLb9d8366obb718veVoY3zblsUl173zz76vtXYxZ4JH0BS3W43Xnjhhdi1a1dUKr965OvIy2MLXrIwa/2H74uhnb+5pD9vcqYTR14aE77lWKtTAgDK6M0334zdu3fH6OhoPProo7Fv375oNBrRanfijdPjS3qNC99/Nibf/MHl77f89p8uuP3rp8ej1e709b279OG7EaYEALyTut1uTE1NxeTk5IKPJ0+ejIGBgRgbG4sDBw7EwYMHY9++ffE3f//1qFcrMdNuL/pnTfz0v3/l+8XCV69WYuzMRNzxrnUr+juuRGnDN9+UgKt/ET85dSG+cfR4PPZvx0s7JQDIodPpxNTU1JKis5zH5T5neno6Go1GNBqNGBwcnPexKIpotVoRETE9PR1FUcSPfvSjmGp1orKUA3wR8e4/+Nslf9QZEVEpinlPmHmnlDJ8V04JWMr/oNkoPvG91+I7P/hZ6aYEANfXbHBWMza9RGdmZmbR2CzlcXh4uOfnDw4OxsDAQBRLCNf4+Hhs2LAhms1mfOITn4jHHnssduzYEa//Yjw61+mE/063G40+L05KF76VTAmYmOnEibcvPb9MUwLgRtVut6/LCme58ZkrOL1EY/PmzSuK1lKDUxbNZjMOHDgQ+/bti127dl3++ejm5pIvNr/6GN/6X78vBm66Y97tZ9qdGN3c3/fmUoVvtacEPPuVj0XdyS/cgGaDcz1js5T4tFqtRVcfS4nGli1bVrRCWmvBKYtKpRKPP/74NT+vVStx+8i6+MmpC4u+xtXH+Abf+2sLhu+OkXV9PymxVOG7HlMC/vK+96/S3sGl4LyTsZkvPu12e1kffc3330ZGRla0wqnX64Jzg7r/rtH4xtHj817SMPqlJ5b9moP1Sjxw9+hKd23FSjO5ZbEpARERk2M/jPP/9S8xdfLV6MxMRm3DSDS33xOb934ximr9mu3LMiWAlWu1Wn2LzZWP7Xa75+Muq3HsZ/ZRcLjezl6cjt/4h4Xfk5erLO/JpVnxLTYlYPyH/xG/eObrEd1O1G/aFkO37ojWuVPxy//91xj+6OfnDF9ZpgSsZXMFpx/R6XQ6PQfkyq83bNiwoujUajXBIYXhoYHYv2dbPPG915Z1K6L5NOvV+OKebX2PXkSJwrfQlIDOzGS8/ew/RXQ7se6DH4+R3/tKFMWlz4hnzvwsinpjzueVZUpAL1qt1qqe3tzrY7fbXZVVyqZNm1b0GoID77yDe3fGd37wszjx9sSKDkFViyJu2diIg3t3rOLe9a4U4VtsSsDU2CvRmfxlRERs+q3PXY5eRET9ihlxc1nOlIBut7usj9SuZ3QiYlU+IrsyOL2ucICcBmqVeGr/7lW5H9+T+3eX5mTDUryrjZ2ZWHBKQPvi2ctf1zbdtKzXbs9Mx6c/98fRPX9qSdEpimLFsWk2mzE8PLyiYz2CA5TBbcPNeObAnh7vwF6NWzY24smSXVtdinfXxaYEVIeGL3/dOncq6iNLPyuoXqvFZ//wgdi+ZfH4CA7AtW4bbsazX/7YnNO05tKsV6LTjfjinm1xcO+O0qz0ZpXiXb5Rqyw4JaBx251RGVwfnckLce4/n/6VY3ytc6eiun5LFNW5/yrVWi1+93c+3te5cABr3UCtEn/1yQ/EI/duvzQ/+aWxeH2O+cl3jKyLB+4ejQfvKe/85FJcztBqd+LOR78bM+35d+XC//17nP72Y5fP6mzcujPaF96Oidf/J7b+xZNRGVw/5/Pq1SJe+btP9f2CSYAbzVq9Y04pVnxLmRKw/oMfj9qGd8f554/E1MlX48LpE1Hb8K7Y8OFPzntWZ0Q5pgQA3Ihq1cqa/DStFOGLWHxKQETE4Hs/FIPv/dCSX7MsUwIAKI/SLIUe2rU1VvtD12434sF7tq7uiwKwppUmfLNTApr11dmlZr0aj9y7vbQHVwHoj9KEL+LSlICbNw5GdYUTOso2JQCA8ihV+GanBGwaqvccvzJOCQCgPEpxOcPVrrwD+40wJQCA8ihl+CIu3ZR2uVMCHrl3eymnBABQHqUN36yzF6fX/JQAAMqj9OG70lqdEgBAeayp8AHASlkuAZCK8AGQivABkIrwAZCK8AGQivABkIrwAZDK/wOSxRITBn40ygAAAABJRU5ErkJggg==\n"
          },
          "metadata": {}
        }
      ]
    }
  ]
}