grapl/tutorials/tutorial_05.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GRAPL: Tutorial 05\n",
    "\n",
    "Tutorial illustrating the use of the GRAPL library, for representing, analyzing and processing\n",
    "**acyclic directed mixed graphs** (**ADMG**s) for structural causal modelling.\n",
    "\n",
    "**License**: CC BY-SA 4.0 2021. \n",
    "\n",
    "\n",
    "If you use this code, please cite:\n",
    "\n",
    "Little, M. A., Badawy, R. (2019). [Causal bootstrapping](https://arxiv.org/pdf/1910.09648.pdf), *ArXiV*, arXiv:1910.09648"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import grapl.dsl as dsl\n",
    "import grapl.admg as admg\n",
    "import grapl.algorithms as algs\n",
    "import copy\n",
    "gobj = dsl.GraplDSL()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read an ADMG from a GRAPL text file, render using Graphviz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 7.0.6 (20230106.0513)\n",
       " -->\n",
       "<!-- Title: &quot;Shpitser thesis, Fig 4.1&quot; Pages: 1 -->\n",
       "<svg width=\"283pt\" height=\"296pt\"\n",
       " viewBox=\"0.00 0.00 283.50 296.29\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 292.29)\">\n",
       "<title>&quot;Shpitser thesis, Fig 4.1&quot;</title>\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-292.29 279.5,-292.29 279.5,4 -4,4\"/>\n",
       "<!-- X -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>X</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"100.75\" cy=\"-184.89\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"95.75\" y=\"-182.19\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">X</text>\n",
       "</g>\n",
       "<!-- Y1 -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>Y1</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"22.75\" cy=\"-108.14\" rx=\"23\" ry=\"23\"/>\n",
       "<text text-anchor=\"start\" x=\"13.25\" y=\"-105.44\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">Y1</text>\n",
       "</g>\n",
       "<!-- X&#45;&gt;Y1 -->\n",
       "<g id=\"edge6\" class=\"edge\">\n",
       "<title>X&#45;&gt;Y1</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M88.22,-171.88C77.16,-161.29 60.73,-145.54 47.1,-132.48\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"49.78,-130.2 40.14,-125.81 44.94,-135.25 49.78,-130.2\"/>\n",
       "</g>\n",
       "<!-- W2 -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>W2</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"174.75\" cy=\"-263.59\" rx=\"24.9\" ry=\"24.9\"/>\n",
       "<text text-anchor=\"start\" x=\"163.75\" y=\"-260.89\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">W2</text>\n",
       "</g>\n",
       "<!-- W2&#45;&gt;X -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>W2&#45;&gt;X</title>\n",
       "<path fill=\"none\" stroke=\"black\" stroke-dasharray=\"5,2\" d=\"M150.06,-237C140.41,-226.99 129.49,-215.69 120.39,-206.25\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"147.44,-239.33 156.91,-244.1 152.48,-234.47 147.44,-239.33\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"123.15,-204.07 113.69,-199.3 118.11,-208.93 123.15,-204.07\"/>\n",
       "</g>\n",
       "<!-- W1 -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>W1</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"137.75\" cy=\"-24.7\" rx=\"24.9\" ry=\"24.9\"/>\n",
       "<text text-anchor=\"start\" x=\"126.75\" y=\"-22\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">W1</text>\n",
       "</g>\n",
       "<!-- W2&#45;&gt;W1 -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>W2&#45;&gt;W1</title>\n",
       "<path fill=\"none\" stroke=\"black\" stroke-dasharray=\"5,2\" d=\"M169.17,-227.61C162.2,-182.6 150.34,-105.98 143.35,-60.88\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"165.71,-228.08 170.7,-237.43 172.62,-227.01 165.71,-228.08\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"146.82,-60.43 141.83,-51.08 139.91,-61.5 146.82,-60.43\"/>\n",
       "</g>\n",
       "<!-- Y2 -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>Y2</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"252.75\" cy=\"-108.14\" rx=\"23\" ry=\"23\"/>\n",
       "<text text-anchor=\"start\" x=\"243.25\" y=\"-105.44\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">Y2</text>\n",
       "</g>\n",
       "<!-- W2&#45;&gt;Y2 -->\n",
       "<g id=\"edge7\" class=\"edge\">\n",
       "<title>W2&#45;&gt;Y2</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M186.09,-240.99C199.62,-214.02 222.3,-168.82 237.4,-138.74\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"240.35,-140.65 241.71,-130.14 234.1,-137.51 240.35,-140.65\"/>\n",
       "</g>\n",
       "<!-- W1&#45;&gt;X -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>W1&#45;&gt;X</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M132.06,-49.33C125.42,-78.07 114.46,-125.53 107.42,-155.99\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"104.03,-155.13 105.19,-165.66 110.85,-156.71 104.03,-155.13\"/>\n",
       "</g>\n",
       "<!-- Y2&#45;&gt;W1 -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>Y2&#45;&gt;W1</title>\n",
       "<path fill=\"none\" stroke=\"black\" stroke-dasharray=\"5,2\" d=\"M225.01,-87.5C207.49,-75.09 184.88,-59.08 167.06,-46.46\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"222.85,-90.26 233.04,-93.18 226.9,-84.55 222.85,-90.26\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"169.14,-43.64 158.96,-40.72 165.1,-49.36 169.14,-43.64\"/>\n",
       "</g>\n",
       "<!-- Y1&#45;&gt;W1 -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>Y1&#45;&gt;W1</title>\n",
       "<path fill=\"none\" stroke=\"black\" stroke-dasharray=\"5,2\" d=\"M50.49,-87.5C68.01,-75.09 90.62,-59.08 108.43,-46.46\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"48.6,-84.55 42.46,-93.18 52.64,-90.26 48.6,-84.55\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"110.4,-49.36 116.54,-40.72 106.35,-43.64 110.4,-49.36\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<grapl.admg.ADMG at 0x2402a71c0d0>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "G = gobj.readgrapl(open('grapl/graphs/shpitser_thesis1.grapl').read())\n",
    "G"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute the interventional distribution, $p_X(Y_2,Y_1,W_2)$ using the fixing algorithm, render resulting Latex formula"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$p_{X}(W2,Y2,Y1)=\\sum_{W1'}[p(Y1|X,W1')p(W2,Y2)p(W1')]$$"
      ],
      "text/plain": [
       "<grapl.eqn.Eqn at 0x2402a7d3f90>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "str, eqn, isident = algs.idfixing(G, {'X'}, {'Y1','Y2','W2'})\n",
    "eqn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Intervene on the variable $X$ and draw the mutilated ADMG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 7.0.6 (20230106.0513)\n",
       " -->\n",
       "<!-- Title: &quot;Shpitser thesis, Fig 4.1&quot; Pages: 1 -->\n",
       "<svg width=\"209pt\" height=\"224pt\"\n",
       " viewBox=\"0.00 0.00 209.50 224.29\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 220.29)\">\n",
       "<title>&quot;Shpitser thesis, Fig 4.1&quot;</title>\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-220.29 205.5,-220.29 205.5,4 -4,4\"/>\n",
       "<!-- X -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>X</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"22.75\" cy=\"-191.59\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"17.75\" y=\"-188.89\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">X</text>\n",
       "</g>\n",
       "<!-- Y1 -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>Y1</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"22.75\" cy=\"-108.14\" rx=\"23\" ry=\"23\"/>\n",
       "<text text-anchor=\"start\" x=\"13.25\" y=\"-105.44\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">Y1</text>\n",
       "</g>\n",
       "<!-- X&#45;&gt;Y1 -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>X&#45;&gt;Y1</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M22.75,-173.32C22.75,-164.43 22.75,-153.28 22.75,-142.78\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"26.25,-142.87 22.75,-132.87 19.25,-142.87 26.25,-142.87\"/>\n",
       "</g>\n",
       "<!-- W1 -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>W1</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"100.75\" cy=\"-24.7\" rx=\"24.9\" ry=\"24.9\"/>\n",
       "<text text-anchor=\"start\" x=\"89.75\" y=\"-22\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">W1</text>\n",
       "</g>\n",
       "<!-- W2 -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>W2</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"139.75\" cy=\"-191.59\" rx=\"24.9\" ry=\"24.9\"/>\n",
       "<text text-anchor=\"start\" x=\"128.75\" y=\"-188.89\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">W2</text>\n",
       "</g>\n",
       "<!-- W2&#45;&gt;W1 -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>W2&#45;&gt;W1</title>\n",
       "<path fill=\"none\" stroke=\"black\" stroke-dasharray=\"5,2\" d=\"M131.4,-155.85C124.81,-127.67 115.67,-88.56 109.09,-60.39\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"127.97,-156.55 133.65,-165.49 134.78,-154.96 127.97,-156.55\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"112.55,-59.84 106.87,-50.9 105.74,-61.43 112.55,-59.84\"/>\n",
       "</g>\n",
       "<!-- Y2 -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>Y2</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"178.75\" cy=\"-108.14\" rx=\"23\" ry=\"23\"/>\n",
       "<text text-anchor=\"start\" x=\"169.25\" y=\"-105.44\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">Y2</text>\n",
       "</g>\n",
       "<!-- W2&#45;&gt;Y2 -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>W2&#45;&gt;Y2</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M150.2,-168.76C154.51,-159.77 159.57,-149.19 164.23,-139.46\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"167.26,-141.25 168.42,-130.72 160.94,-138.22 167.26,-141.25\"/>\n",
       "</g>\n",
       "<!-- Y2&#45;&gt;W1 -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>Y2&#45;&gt;W1</title>\n",
       "<path fill=\"none\" stroke=\"black\" stroke-dasharray=\"5,2\" d=\"M155.41,-82.78C145.96,-72.9 135.02,-61.49 125.43,-51.47\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"152.85,-85.16 162.29,-89.96 157.9,-80.32 152.85,-85.16\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"128.24,-49.35 118.8,-44.55 123.19,-54.19 128.24,-49.35\"/>\n",
       "</g>\n",
       "<!-- Y1&#45;&gt;W1 -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>Y1&#45;&gt;W1</title>\n",
       "<path fill=\"none\" stroke=\"black\" stroke-dasharray=\"5,2\" d=\"M46.08,-82.78C55.54,-72.9 66.47,-61.49 76.07,-51.47\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"43.59,-80.32 39.21,-89.96 48.65,-85.16 43.59,-80.32\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"78.31,-54.19 82.7,-44.55 73.25,-49.35 78.31,-54.19\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<grapl.admg.ADMG at 0x2402a712690>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "G_fixX = copy.deepcopy(G)\n",
    "G_fixX.fix('X')\n",
    "G_fixX"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read another ADMG from a file, draw it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 7.0.6 (20230106.0513)\n",
       " -->\n",
       "<!-- Title: &quot;Bareinboim 2020, Fig. 1.8&quot; Pages: 1 -->\n",
       "<svg width=\"264pt\" height=\"332pt\"\n",
       " viewBox=\"0.00 0.00 264.00 332.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 328)\">\n",
       "<title>&quot;Bareinboim 2020, Fig. 1.8&quot;</title>\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-328 260,-328 260,4 -4,4\"/>\n",
       "<!-- A -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>A</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"165\" cy=\"-306\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"160\" y=\"-303.3\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">A</text>\n",
       "</g>\n",
       "<!-- B -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>B</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"91\" cy=\"-234\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"86\" y=\"-231.3\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">B</text>\n",
       "</g>\n",
       "<!-- A&#45;&gt;B -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>A&#45;&gt;B</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M152.11,-292.81C141.07,-282.37 124.97,-267.13 112.06,-254.93\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"114.51,-252.42 104.84,-248.09 109.7,-257.51 114.51,-252.42\"/>\n",
       "</g>\n",
       "<!-- C -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>C</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"91\" cy=\"-18\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"85.5\" y=\"-15.3\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">C</text>\n",
       "</g>\n",
       "<!-- A&#45;&gt;C -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>A&#45;&gt;C</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M160.48,-288.42C148.12,-240.32 113.66,-106.19 98.42,-46.86\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"101.81,-46 95.93,-37.18 95.03,-47.74 101.81,-46\"/>\n",
       "</g>\n",
       "<!-- F -->\n",
       "<g id=\"node6\" class=\"node\">\n",
       "<title>F</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"238\" cy=\"-234\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"233\" y=\"-231.3\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">F</text>\n",
       "</g>\n",
       "<!-- A&#45;&gt;F -->\n",
       "<g id=\"edge8\" class=\"edge\">\n",
       "<title>A&#45;&gt;F</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M177.72,-292.81C188.6,-282.37 204.49,-267.13 217.22,-254.93\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"219.54,-257.55 224.34,-248.1 214.7,-252.5 219.54,-257.55\"/>\n",
       "</g>\n",
       "<!-- D -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>D</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"91\" cy=\"-162\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"85.5\" y=\"-159.3\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">D</text>\n",
       "</g>\n",
       "<!-- B&#45;&gt;D -->\n",
       "<g id=\"edge6\" class=\"edge\">\n",
       "<title>B&#45;&gt;D</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M79.99,-219.17C77.58,-210.25 76.92,-198.4 78.01,-187.88\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"81.41,-188.75 79.73,-178.29 74.52,-187.52 81.41,-188.75\"/>\n",
       "</g>\n",
       "<!-- D&#45;&gt;B -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>D&#45;&gt;B</title>\n",
       "<path fill=\"none\" stroke=\"black\" stroke-dasharray=\"5,2\" d=\"M103.98,-187.84C104.67,-194.41 104.67,-201.51 103.99,-208.08\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"107.48,-187.52 102.27,-178.29 100.59,-188.75 107.48,-187.52\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"100.59,-207.22 102.27,-217.68 107.48,-208.45 100.59,-207.22\"/>\n",
       "</g>\n",
       "<!-- E -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>E</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"18\" cy=\"-90\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"13\" y=\"-87.3\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">E</text>\n",
       "</g>\n",
       "<!-- D&#45;&gt;E -->\n",
       "<g id=\"edge7\" class=\"edge\">\n",
       "<title>D&#45;&gt;E</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M78.28,-148.81C67.4,-138.37 51.51,-123.13 38.78,-110.93\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"41.3,-108.5 31.66,-104.1 36.46,-113.55 41.3,-108.5\"/>\n",
       "</g>\n",
       "<!-- C&#45;&gt;D -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>C&#45;&gt;D</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M91,-36.09C91,-60.25 91,-103.36 91,-132.32\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"87.5,-132.25 91,-142.25 94.5,-132.25 87.5,-132.25\"/>\n",
       "</g>\n",
       "<!-- E&#45;&gt;C -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>E&#45;&gt;C</title>\n",
       "<path fill=\"none\" stroke=\"black\" stroke-dasharray=\"5,2\" d=\"M38.95,-68.91C48.82,-59.45 60.55,-48.19 70.38,-38.77\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"36.6,-66.31 31.81,-75.76 41.45,-71.36 36.6,-66.31\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"72.54,-41.55 77.34,-32.1 67.7,-36.5 72.54,-41.55\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<grapl.admg.ADMG at 0x2402a7340d0>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "G = gobj.readgrapl(open('grapl/graphs/bareinboim_2020.grapl').read())\n",
    "G"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute the Tian factorized joint distribution over variables $F,B,D$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$p(D,F,B)=\\sum_{A,C}[p(D|A,B,C)p(B|A)p(C|A)p(A,F)]$$"
      ],
      "text/plain": [
       "<grapl.eqn.Eqn at 0x2402a7e1fd0>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "str, eqn = algs.admgfactor(G,{'B','F','D'})\n",
    "eqn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Intervene on the variable $C$, render the mutilated graph, and compute the corresponding Tian factorization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 7.0.6 (20230106.0513)\n",
       " -->\n",
       "<!-- Title: &quot;Bareinboim 2020, Fig. 1.8&quot; Pages: 1 -->\n",
       "<svg width=\"188pt\" height=\"260pt\"\n",
       " viewBox=\"0.00 0.00 188.00 260.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 256)\">\n",
       "<title>&quot;Bareinboim 2020, Fig. 1.8&quot;</title>\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-256 184,-256 184,4 -4,4\"/>\n",
       "<!-- A -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>A</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"54\" cy=\"-234\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"49\" y=\"-231.3\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">A</text>\n",
       "</g>\n",
       "<!-- B -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>B</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"18\" cy=\"-162\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"13\" y=\"-159.3\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">B</text>\n",
       "</g>\n",
       "<!-- A&#45;&gt;B -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>A&#45;&gt;B</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M46.01,-217.46C41.63,-208.95 36.09,-198.17 31.09,-188.46\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"34.35,-187.13 26.66,-179.84 28.12,-190.33 34.35,-187.13\"/>\n",
       "</g>\n",
       "<!-- F -->\n",
       "<g id=\"node6\" class=\"node\">\n",
       "<title>F</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"90\" cy=\"-162\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"85\" y=\"-159.3\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">F</text>\n",
       "</g>\n",
       "<!-- A&#45;&gt;F -->\n",
       "<g id=\"edge6\" class=\"edge\">\n",
       "<title>A&#45;&gt;F</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M61.99,-217.46C66.37,-208.95 71.91,-198.17 76.91,-188.46\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"79.88,-190.33 81.34,-179.84 73.65,-187.13 79.88,-190.33\"/>\n",
       "</g>\n",
       "<!-- D -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>D</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"18\" cy=\"-90\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"12.5\" y=\"-87.3\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">D</text>\n",
       "</g>\n",
       "<!-- B&#45;&gt;D -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>B&#45;&gt;D</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M6.99,-147.17C4.58,-138.25 3.92,-126.4 5.01,-115.88\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"8.41,-116.75 6.73,-106.29 1.52,-115.52 8.41,-116.75\"/>\n",
       "</g>\n",
       "<!-- D&#45;&gt;B -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>D&#45;&gt;B</title>\n",
       "<path fill=\"none\" stroke=\"black\" stroke-dasharray=\"5,2\" d=\"M30.98,-115.84C31.67,-122.41 31.67,-129.51 30.99,-136.08\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"34.48,-115.52 29.27,-106.29 27.59,-116.75 34.48,-115.52\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"27.59,-135.22 29.27,-145.68 34.48,-136.45 27.59,-135.22\"/>\n",
       "</g>\n",
       "<!-- E -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>E</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"18\" cy=\"-18\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"13\" y=\"-15.3\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">E</text>\n",
       "</g>\n",
       "<!-- D&#45;&gt;E -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>D&#45;&gt;E</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M18,-71.7C18,-64.41 18,-55.73 18,-47.54\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"21.5,-47.62 18,-37.62 14.5,-47.62 21.5,-47.62\"/>\n",
       "</g>\n",
       "<!-- C -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>C</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"162\" cy=\"-162\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"156.5\" y=\"-159.3\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">C</text>\n",
       "</g>\n",
       "<!-- C&#45;&gt;D -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>C&#45;&gt;D</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M146.04,-153.24C121.39,-141.26 73.8,-118.12 44.09,-103.68\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"45.92,-100.68 35.39,-99.46 42.86,-106.97 45.92,-100.68\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<grapl.admg.ADMG at 0x2402a7e2c90>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "G_fixC = copy.deepcopy(G)\n",
    "G_fixC.fix('C')\n",
    "G_fixC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$p(D,F,B)=\\sum_{A,C}[p(D|A,B,C)p(B|A)p(A,F)p(C)]$$"
      ],
      "text/plain": [
       "<grapl.eqn.Eqn at 0x2402a7f4510>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "str, eqn = algs.admgfactor(G_fixC,{'B','F','D'})\n",
    "eqn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read a DAG and draw it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 7.0.6 (20230106.0513)\n",
       " -->\n",
       "<!-- Title: &quot;Markov independence example&quot; Pages: 1 -->\n",
       "<svg width=\"116pt\" height=\"188pt\"\n",
       " viewBox=\"0.00 0.00 116.20 188.40\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 184.4)\">\n",
       "<title>&quot;Markov independence example&quot;</title>\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-184.4 112.2,-184.4 112.2,4 -4,4\"/>\n",
       "<!-- D -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>D</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"18.2\" cy=\"-162.4\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"12.7\" y=\"-159.7\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">D</text>\n",
       "</g>\n",
       "<!-- G -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>G</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"18.2\" cy=\"-90.2\" rx=\"18.4\" ry=\"18.4\"/>\n",
       "<text text-anchor=\"start\" x=\"12.2\" y=\"-87.5\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">G</text>\n",
       "</g>\n",
       "<!-- D&#45;&gt;G -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>D&#45;&gt;G</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M18.2,-144.04C18.2,-136.84 18.2,-128.29 18.2,-120.18\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"21.7,-120.35 18.2,-110.35 14.7,-120.35 21.7,-120.35\"/>\n",
       "</g>\n",
       "<!-- I -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>I</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"90.2\" cy=\"-162.4\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"88.2\" y=\"-159.7\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">I</text>\n",
       "</g>\n",
       "<!-- I&#45;&gt;G -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>I&#45;&gt;G</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M77.66,-149.17C67.05,-138.83 51.63,-123.8 39.15,-111.63\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"41.77,-109.29 32.16,-104.82 36.88,-114.3 41.77,-109.29\"/>\n",
       "</g>\n",
       "<!-- S -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>S</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"90.2\" cy=\"-90.2\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"85.2\" y=\"-87.5\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">S</text>\n",
       "</g>\n",
       "<!-- I&#45;&gt;S -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>I&#45;&gt;S</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M90.2,-144.04C90.2,-136.73 90.2,-128.03 90.2,-119.81\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"93.7,-119.86 90.2,-109.86 86.7,-119.86 93.7,-119.86\"/>\n",
       "</g>\n",
       "<!-- L -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>L</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"18.2\" cy=\"-18\" rx=\"18\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"13.7\" y=\"-15.3\" font-family=\"Helvetica,sans-Serif\" font-style=\"italic\" font-size=\"14.00\">L</text>\n",
       "</g>\n",
       "<!-- G&#45;&gt;L -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>G&#45;&gt;L</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M18.2,-71.48C18.2,-64.28 18.2,-55.77 18.2,-47.72\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"21.7,-47.98 18.2,-37.98 14.7,-47.98 21.7,-47.98\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<grapl.admg.ADMG at 0x2402a664990>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "G = gobj.readgrapl(open('grapl/graphs/markov1.grapl').read())\n",
    "G"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute all local Markov conditional independences, render in Latex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$(L \\perp\\!\\!\\!\\perp S,D,I \\mid G),(I \\perp\\!\\!\\!\\perp D),(D \\perp\\!\\!\\!\\perp S,I),(G \\perp\\!\\!\\!\\perp S \\mid D,I),(S \\perp\\!\\!\\!\\perp D,L,G \\mid I)$$"
      ],
      "text/plain": [
       "<grapl.condind.CondIndSet at 0x2402a70f3d0>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ciset, cistr, isdag = algs.localmarkov(G)\n",
    "ciset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test whether $D$ is d-separated from $S$ given $G,I$, render in Latex if yes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isdsep, ci, cistr, isdag = algs.dseparate(G,{'D'},{'S'},{'G','I'})\n",
    "isdsep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$(D \\perp\\!\\!\\!\\perp S \\mid I,G)$$"
      ],
      "text/plain": [
       "<grapl.condind.CondInd at 0x2402a681bd0>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ci"
   ]
  }
 ],
 "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.9.7"
  },
  "vscode": {
   "interpreter": {
    "hash": "5238573367df39f7286bb46f9ff5f08f63a01a80960060ce41e3c79b190280fa"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}