notebooks/cgm-examples.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# An Introduction to `CausalGraphicalModels`\n",
    "\n",
    "`CausalGraphicalModel` is a python module for describing and manipulating [Causal Graphical Models](https://en.wikipedia.org/wiki/Causal_graph) and [Structural Causal Models](https://en.wikipedia.org/wiki/Structural_equation_modeling). Behind the curtain, it is a light wrapper around the python graph library [networkx](https://networkx.github.io/).\n",
    "\n",
    "This notebook is designed to give a quick overview of the functionality of this package."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CausalGraphicalModels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from causalgraphicalmodels import CausalGraphicalModel"
   ]
  },
  {
   "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 2.38.0 (20140413.2041)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"164pt\" height=\"260pt\"\n",
       " viewBox=\"0.00 0.00 164.25 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>%3</title>\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-256 160.246,-256 160.246,4 -4,4\"/>\n",
       "<!-- sprinkler -->\n",
       "<g id=\"node1\" class=\"node\"><title>sprinkler</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"42.2463\" cy=\"-162\" rx=\"42.4939\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"42.2463\" y=\"-158.3\" font-family=\"Times,serif\" font-size=\"14.00\">sprinkler</text>\n",
       "</g>\n",
       "<!-- wet -->\n",
       "<g id=\"node3\" class=\"node\"><title>wet</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"85.2463\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"85.2463\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">wet</text>\n",
       "</g>\n",
       "<!-- sprinkler&#45;&gt;wet -->\n",
       "<g id=\"edge4\" class=\"edge\"><title>sprinkler&#45;&gt;wet</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M52.4367,-144.411C57.7064,-135.832 64.2434,-125.191 70.0921,-115.67\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"73.173,-117.341 75.425,-106.988 67.2085,-113.677 73.173,-117.341\"/>\n",
       "</g>\n",
       "<!-- rain -->\n",
       "<g id=\"node2\" class=\"node\"><title>rain</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"129.246\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"129.246\" y=\"-158.3\" font-family=\"Times,serif\" font-size=\"14.00\">rain</text>\n",
       "</g>\n",
       "<!-- rain&#45;&gt;wet -->\n",
       "<g id=\"edge3\" class=\"edge\"><title>rain&#45;&gt;wet</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M119.262,-145.116C113.749,-136.345 106.783,-125.264 100.594,-115.416\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"103.476,-113.425 95.191,-106.821 97.5495,-117.15 103.476,-113.425\"/>\n",
       "</g>\n",
       "<!-- slippery -->\n",
       "<g id=\"node5\" class=\"node\"><title>slippery</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"85.2463\" cy=\"-18\" rx=\"38.9931\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"85.2463\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">slippery</text>\n",
       "</g>\n",
       "<!-- wet&#45;&gt;slippery -->\n",
       "<g id=\"edge5\" class=\"edge\"><title>wet&#45;&gt;slippery</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M85.2463,-71.6966C85.2463,-63.9827 85.2463,-54.7125 85.2463,-46.1124\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"88.7464,-46.1043 85.2463,-36.1043 81.7464,-46.1044 88.7464,-46.1043\"/>\n",
       "</g>\n",
       "<!-- season -->\n",
       "<g id=\"node4\" class=\"node\"><title>season</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"85.2463\" cy=\"-234\" rx=\"34.394\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"85.2463\" y=\"-230.3\" font-family=\"Times,serif\" font-size=\"14.00\">season</text>\n",
       "</g>\n",
       "<!-- season&#45;&gt;sprinkler -->\n",
       "<g id=\"edge2\" class=\"edge\"><title>season&#45;&gt;sprinkler</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M75.2731,-216.765C70.084,-208.317 63.623,-197.799 57.7952,-188.312\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"60.6804,-186.322 52.4639,-179.633 54.7159,-189.986 60.6804,-186.322\"/>\n",
       "</g>\n",
       "<!-- season&#45;&gt;rain -->\n",
       "<g id=\"edge1\" class=\"edge\"><title>season&#45;&gt;rain</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M95.4513,-216.765C100.838,-208.195 107.564,-197.494 113.595,-187.9\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"116.738,-189.476 119.097,-179.147 110.812,-185.751 116.738,-189.476\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x7f10a219fb38>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sprinkler = CausalGraphicalModel(\n",
    "    nodes=[\"season\", \"rain\", \"sprinkler\", \"wet\", \"slippery\"],\n",
    "    edges=[\n",
    "        (\"season\", \"rain\"), \n",
    "        (\"season\", \"sprinkler\"), \n",
    "        (\"rain\", \"wet\"),\n",
    "        (\"sprinkler\", \"wet\"), \n",
    "        (\"wet\", \"slippery\")\n",
    "    ]\n",
    ")\n",
    "\n",
    "# draw return a graphviz `dot` object, which jupyter can render\n",
    "sprinkler.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P(season)P(sprinkler|season)P(rain|season)P(wet|rain,sprinkler)P(slippery|wet)\n"
     ]
    }
   ],
   "source": [
    "# get the distribution implied by the graph\n",
    "print(sprinkler.get_distribution())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# check for d-seperation of two nodes\n",
    "sprinkler.is_d_separated(\"slippery\", \"season\", {\"wet\"})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('sprinkler', 'rain', {'season'}),\n",
       " ('sprinkler', 'slippery', {'wet'}),\n",
       " ('sprinkler', 'slippery', {'rain', 'wet'}),\n",
       " ('sprinkler', 'slippery', {'season', 'wet'}),\n",
       " ('sprinkler', 'slippery', {'rain', 'season', 'wet'}),\n",
       " ('rain', 'slippery', {'wet'}),\n",
       " ('rain', 'slippery', {'sprinkler', 'wet'}),\n",
       " ('rain', 'slippery', {'season', 'wet'}),\n",
       " ('rain', 'slippery', {'season', 'sprinkler', 'wet'}),\n",
       " ('wet', 'season', {'rain', 'sprinkler'}),\n",
       " ('wet', 'season', {'rain', 'slippery', 'sprinkler'}),\n",
       " ('season', 'slippery', {'wet'}),\n",
       " ('season', 'slippery', {'rain', 'wet'}),\n",
       " ('season', 'slippery', {'sprinkler', 'wet'}),\n",
       " ('season', 'slippery', {'rain', 'sprinkler'}),\n",
       " ('season', 'slippery', {'rain', 'sprinkler', 'wet'})]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# get all the conditional independence relationships implied by a CGM\n",
    "sprinkler.get_all_independence_relationships()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# check backdoor adjustment set\n",
    "sprinkler.is_valid_backdoor_adjustment_set(\"rain\", \"slippery\", {\"wet\"})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "frozenset({frozenset({'sprinkler'}),\n",
       "           frozenset({'season'}),\n",
       "           frozenset({'season', 'sprinkler'})})"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# get all backdoor adjustment sets\n",
    "sprinkler.get_all_backdoor_adjustment_sets(\"rain\", \"slippery\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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 2.38.0 (20140413.2041)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"172pt\" height=\"268pt\"\n",
       " viewBox=\"0.00 0.00 172.25 268.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 264)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-264 168.246,-264 168.246,4 -4,4\"/>\n",
       "<!-- sprinkler -->\n",
       "<g id=\"node1\" class=\"node\"><title>sprinkler</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"122\" cy=\"-166\" rx=\"42.4939\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"122\" y=\"-162.3\" font-family=\"Times,serif\" font-size=\"14.00\">sprinkler</text>\n",
       "</g>\n",
       "<!-- wet -->\n",
       "<g id=\"node3\" class=\"node\"><title>wet</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"76\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"76\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">wet</text>\n",
       "</g>\n",
       "<!-- sprinkler&#45;&gt;wet -->\n",
       "<g id=\"edge1\" class=\"edge\"><title>sprinkler&#45;&gt;wet</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M111.562,-148.208C105.598,-138.614 98.0078,-126.404 91.3728,-115.73\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"94.3326,-113.862 86.0807,-107.217 88.3876,-117.557 94.3326,-113.862\"/>\n",
       "</g>\n",
       "<!-- rain -->\n",
       "<g id=\"node2\" class=\"node\"><title>rain</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"31\" cy=\"-166\" rx=\"27\" ry=\"18\"/>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"31\" cy=\"-166\" rx=\"31\" ry=\"22\"/>\n",
       "<text text-anchor=\"middle\" x=\"31\" y=\"-162.3\" font-family=\"Times,serif\" font-size=\"14.00\">rain</text>\n",
       "</g>\n",
       "<!-- rain&#45;&gt;wet -->\n",
       "<g id=\"edge2\" class=\"edge\"><title>rain&#45;&gt;wet</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M42.8232,-145.557C48.3487,-136.471 54.9925,-125.546 60.8787,-115.866\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"64.0234,-117.431 66.2287,-107.068 58.0424,-113.794 64.0234,-117.431\"/>\n",
       "</g>\n",
       "<!-- slippery -->\n",
       "<g id=\"node5\" class=\"node\"><title>slippery</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"76\" cy=\"-18\" rx=\"38.9931\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"76\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">slippery</text>\n",
       "</g>\n",
       "<!-- wet&#45;&gt;slippery -->\n",
       "<g id=\"edge3\" class=\"edge\"><title>wet&#45;&gt;slippery</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M76,-71.6966C76,-63.9827 76,-54.7125 76,-46.1124\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"79.5001,-46.1043 76,-36.1043 72.5001,-46.1044 79.5001,-46.1043\"/>\n",
       "</g>\n",
       "<!-- season -->\n",
       "<g id=\"node4\" class=\"node\"><title>season</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"122\" cy=\"-242\" rx=\"34.394\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"122\" y=\"-238.3\" font-family=\"Times,serif\" font-size=\"14.00\">season</text>\n",
       "</g>\n",
       "<!-- season&#45;&gt;sprinkler -->\n",
       "<g id=\"edge4\" class=\"edge\"><title>season&#45;&gt;sprinkler</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M122,-223.837C122,-215.058 122,-204.163 122,-194.267\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"125.5,-194.07 122,-184.07 118.5,-194.07 125.5,-194.07\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x7f10a219f550>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# get the graph created by intervening on node \"rain\"\n",
    "do_sprinkler = sprinkler.do(\"rain\")\n",
    "\n",
    "do_sprinkler.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Latent Variables "
   ]
  },
  {
   "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 2.38.0 (20140413.2041)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"90pt\" height=\"228pt\"\n",
       " viewBox=\"0.00 0.00 90.00 227.60\" 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 223.6)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-223.6 86,-223.6 86,4 -4,4\"/>\n",
       "<!-- y -->\n",
       "<g id=\"node1\" class=\"node\"><title>y</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">y</text>\n",
       "</g>\n",
       "<!-- z -->\n",
       "<g id=\"node2\" class=\"node\"><title>z</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"55\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"55\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">z</text>\n",
       "</g>\n",
       "<!-- z&#45;&gt;y -->\n",
       "<g id=\"edge2\" class=\"edge\"><title>z&#45;&gt;y</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M48.3644,-72.411C45.0866,-64.2164 41.0558,-54.1395 37.3819,-44.9548\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"40.5516,-43.4548 33.588,-35.4699 34.0523,-46.0546 40.5516,-43.4548\"/>\n",
       "</g>\n",
       "<!-- x -->\n",
       "<g id=\"node3\" class=\"node\"><title>x</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"55\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"55\" y=\"-158.3\" font-family=\"Times,serif\" font-size=\"14.00\">x</text>\n",
       "</g>\n",
       "<!-- x&#45;&gt;z -->\n",
       "<g id=\"edge1\" class=\"edge\"><title>x&#45;&gt;z</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M55,-143.697C55,-135.983 55,-126.712 55,-118.112\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"58.5001,-118.104 55,-108.104 51.5001,-118.104 58.5001,-118.104\"/>\n",
       "</g>\n",
       "<!-- Unobserved_0 -->\n",
       "<g id=\"node4\" class=\"node\"><title>Unobserved_0</title>\n",
       "<ellipse fill=\"black\" stroke=\"black\" cx=\"27\" cy=\"-217.8\" rx=\"1.8\" ry=\"1.8\"/>\n",
       "</g>\n",
       "<!-- Unobserved_0&#45;&gt;y -->\n",
       "<g id=\"edge4\" class=\"edge\"><title>Unobserved_0&#45;&gt;y</title>\n",
       "<path fill=\"none\" stroke=\"black\" stroke-dasharray=\"5,2\" d=\"M26.7296,-215.901C25.5161,-211.815 20.5645,-194.562 19,-180 13.8723,-132.275 15.4079,-119.865 19,-72 19.6366,-63.5179 20.8384,-54.3361 22.1208,-46.0356\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"25.5836,-46.5472 23.7517,-36.1119 18.6763,-45.412 25.5836,-46.5472\"/>\n",
       "</g>\n",
       "<!-- Unobserved_0&#45;&gt;x -->\n",
       "<g id=\"edge3\" class=\"edge\"><title>Unobserved_0&#45;&gt;x</title>\n",
       "<path fill=\"none\" stroke=\"black\" stroke-dasharray=\"5,2\" d=\"M27.486,-215.866C29.2799,-212.419 35.7477,-199.992 41.9617,-188.052\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"45.0705,-189.66 46.5825,-179.174 38.8611,-186.428 45.0705,-189.66\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x7f10a219fa58>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dag_with_latent_variables = CausalGraphicalModel(\n",
    "    nodes=[\"x\", \"y\", \"z\"],\n",
    "    edges=[\n",
    "        (\"x\", \"z\"),\n",
    "        (\"z\", \"y\"), \n",
    "    ],\n",
    "    latent_edges=[\n",
    "        (\"x\", \"y\")\n",
    "    ]\n",
    ")\n",
    "\n",
    "dag_with_latent_variables.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "frozenset()"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# here there are no observed backdoor adjustment sets\n",
    "dag_with_latent_variables.get_all_backdoor_adjustment_sets(\"x\", \"y\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "frozenset({frozenset({'z'})})"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# but there is a frontdoor adjustment set\n",
    "dag_with_latent_variables.get_all_frontdoor_adjustment_sets(\"x\", \"y\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# StructuralCausalModels\n",
    "\n",
    "For Structural Causal Models (SCM) we need to specify the functional form of each node:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from causalgraphicalmodels import StructuralCausalModel\n",
    "import numpy as np\n",
    "\n",
    "scm = StructuralCausalModel({\n",
    "    \"x1\": lambda     n_samples: np.random.binomial(n=1,p=0.7,size=n_samples),\n",
    "    \"x2\": lambda x1, n_samples: np.random.normal(loc=x1, scale=0.1),\n",
    "    \"x3\": lambda x2, n_samples: x2 ** 2,\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The only requirement on the functions are:\n",
    " - that variable names are consistent \n",
    " - each function accepts keyword variables in the form of `numpy` arrays and output numpy arrays of shape [n_samples] \n",
    " - that in addition to it's parents, each function takes a `n_samples` variables indicating how many samples to generate \n",
    " - that any function acts on each row independently. This ensure that the output samples are independent\n",
    " \n",
    "Wrapping these functions in the `StructuralCausalModel` object allows us to easily generate samples: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>x1</th>\n",
       "      <th>x2</th>\n",
       "      <th>x3</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>1.004020</td>\n",
       "      <td>1.008056</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>0.902045</td>\n",
       "      <td>0.813686</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>1.150323</td>\n",
       "      <td>1.323242</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>0.984634</td>\n",
       "      <td>0.969504</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>1.023118</td>\n",
       "      <td>1.046770</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   x1        x2        x3\n",
       "0   1  1.004020  1.008056\n",
       "1   1  0.902045  0.813686\n",
       "2   1  1.150323  1.323242\n",
       "3   1  0.984634  0.969504\n",
       "4   1  1.023118  1.046770"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds = scm.sample(n_samples=100)\n",
    "\n",
    "ds.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7f109cf51518>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3Xd0VNXawOHfSTLJpPdeSEgBQodE\negcBASliwYKAXgH1+qko2Ati9yKWa+EqqCjiFZEqSO8lJJQQAiGF9F4myWQykynn+yNxwphCLJdQ\n9rMWCzhnz8w+yVrnnbP3ft8tybKMIAiCILSFVXt3QBAEQbh2iKAhCIIgtJkIGoIgCEKbiaAhCIIg\ntJkIGoIgCEKbiaAhCIIgtJkIGoIgCEKbiaAhCIIgtJkIGoIgCEKb2bR3B/5uXl5ecmhoaHt3QxAE\n4ZqSkJBQKsuy9+XaXXdBIzQ0lPj4+PbuhiAIwjVFkqSstrQTw1OCIAhCm4mgIQiCILSZCBqCIAhC\nm4mgIQiCILSZCBqCIAhCm4mgIQiCILSZCBqCIAhCm4mgIQiCILSZCBqCIAhCm4mgIQiCILSZCBqC\nIAhCm4mgIQiCILSZCBqCIAhCm4mgIQiCILSZCBqCIAhCm4mgIQiCILSZCBqCIAhCm4mgIQiCILRZ\nuwUNSZKCJUnaI0lSsiRJZyVJ+r9m2kiSJH0oSVKaJEmJkiT1aY++CoIgCPXac49wA7BAluUTkiQ5\nAwmSJO2QZTn5kjbjgciGP/2ATxv+FgRBENpBuz1pyLJcIMvyiYZ/VwPngMDfNZsMfCPXOwq4SZLk\nf4W7KgiCIDS4KuY0JEkKBXoDx353KhDIueT/uTQNLEiS9JAkSfGSJMWXlJT8r7opCIJww2v3oCFJ\nkhPwE/C4LMtVf+Y9ZFleLstyjCzLMd7e3n9vBwVBEASzdg0akiQpqA8Y38myvK6ZJnlA8CX/D2o4\nJgiCILSD9lw9JQFfAudkWV7aQrONwMyGVVT9gUpZlguuWCcFQRAEC+25emoQcB9wRpKkUw3HngNC\nAGRZ/gz4BbgFSAM0wOx26KcgCILQoN2ChizLBwHpMm1k4JEr0yNBEAThctp9IlwQBEG4doigIQiC\nILSZCBqCIAhCm4mgIQiCILSZCBqCIAhCm4mgIQiCILSZCBqCIAhCm4mgIQiCILSZCBqCIAhCm4mg\nIQiCILSZCBqCIAhCm4mgIQiCILSZCBqCIAhCm4mgIQiCILSZCBqCIAhCm4mgIQiCILSZCBqCIAhC\nm7Xndq+CIAhtotfrqampwWQyYTKZkGUZe3t7HBwcsLIS332vJBE0BEFod7W1taSlpZGamkpubi75\n+fnk5+dTVlZGVVUVtbW1zb5OkiScnJxwd3cnODiYoKAggoOD6datG506dUKhUFzhK7n+iaAhCMIV\npdVqOX/+PGfOnCExMZELFy6Qn5+PLMsA2NjYEBAQQEBAAB07dsTV1RVnZ2ecnJywtrZGkiQkSaK2\ntha1Wk11dTUlJSXk5uZy8uRJNBoNALa2tnTu3Jl+/foxcuRIIiIikCSpPS/9uiD99ou6XsTExMjx\n8fHt3Q1BEBrIskxGRgZHjhzh8OHDnDx5Er1eD0BgYCDR0dGEh4cTERFBZGQkfn5+WFtb/+nPKi4u\n5syZMyQlJZGYmEhSUhImk4ng4GBGjx7Nbbfdhp+f3995idcFSZISZFmOuWw7ETQEQfi7GQwGTp06\nxYEDB9i3bx+5ubkAdOzYkQEDBtC3b1+6deuGh4fH/7wv5eXl7Nu3j927dxMXF4csy4waNYqZM2fS\npUuX//nnXytE0BAE4YqSZZmUlBS2bNnCr7/+Snl5OQqFgtjYWIYPH87AgQPb/Rt+YWEhP/zwA+vW\nraOmpoZRo0bx+OOP4+/v3679uhqIoCEIwhWhUqnYtm0b69evJy0tDYVCwZAhQxg7diz9+/fH0dGx\nvbvYhFqt5vvvv+err74CYPbs2cyaNQsbmxt3mlcEDUEQ/mdkWeb06dP88MMP7N27F71eT5cuXZg8\neTJjxozB1dW1vbvYJoWFhSxdupTdu3cTHR3NW2+9RUBAQHt3q12IoCEIwt/OYDCwZ88evvvuO5KS\nknBxcWH8+PFMnjyZqKio9u7en7Zr1y5ee+01rKyseP311xkwYEB7d+mKE0FDEIS/jU6nY+PGjXz9\n9dcUFhYSFBTEPffcw6RJk1Aqle3dvb9FTk4OCxcuJCMjg+eff55bb721vbt0RbU1aNy4A3iCIFyW\nVqtl3bp1rFq1ipKSEnr06MHTTz/N4MGD//Sy2N8zGI2UqcopLiuhpKwUtaaGmloNtdpaTLKMrUKB\nwkaBvVKJt4cXvp7e+Hh64+Lk/Ld8/m+Cg4P54osvePrpp1m8eDE6nY7bb7/9b/2M60G7Bg1JklYA\nE4FiWZa7NXN+OLABuNhwaJ0sy4uvXA8F4cZUV1fHhg0bWLFiBSUlJfTt25fFixcTExPzlxLkispK\nOH0+idTMdDJyMsnIySS/uBCTyfSH38vbw5Mu4Z3oEt6J3l2607NLd2z+YiBzdHTkgw8+YNGiRbz9\n9tu4ubkxZsyYv/Se15t2HZ6SJGkooAa+aSVoPCXL8sS2vqcYnhKEP89kMvHLL7/w+eefU1BQQK9e\nvZg/fz59+/b9U+9XUani8Mk4jpw6zulzZygsLQbA2tqaEP8gwoNDCQkMxtfTBx9PL7w9vHBxdMLB\n3gEHpT2SlRV6fR11ej01tRpKykspLiulsLSI1Mx0zqVfICs/B1mWcXNxZfhNgxk9cBix3fv8pZpU\nOp2O+fPnk5KSwueff063bk1uT9eda2ZOQ5KkUGCzCBqC0L7i4+NZtmwZ58+fJzo6mocffph+/fr9\noScLWZZJz8lk56G9HDpxlJSLaciyjIebO3269KBH52707NyNyA4d/7a6UGpNDXGJCew+sp8D8UfQ\naGsJDQxh5pS7GD909J9eRltRUcGsWbPQarWsWbMGd3f3v6W/V6vrKWj8BOQC+dQHkLPNtHsIeAgg\nJCSkb1ZW1v+wx4JwfcnLy2Pp0qXs27cPX19fHn30UcaOHfuHvqnnFxeyec82dh7ex8XcLKysrOjR\nqSv9esYwqE8/OoVFXJFqtLq6OvYeO8DX69eQmpmOr5cP/7h9JpNGjvtTn5+amsrMmTMZNGgQ7777\n7nVdu+p6CRougEmWZbUkSbcAH8iyHNna+4knDUFom7q6OlatWsWKFSuwsrJizpw5zJgxo82roXR1\ndew5doDNe7YRl3gCgD7RPRk9cBgj+g/B0+1/XyKkJbIsc/hkHF/+uIozF5LpHd2DV//5DP4+fzwj\n/dtvv2XZsmUsXryYW2655X/Q26vDdRE0mmmbCcTIslzaUhsRNATh8uLj43njjTfIzs5m5MiRLFiw\nAF9f3za9tkxVzk+/bmTtto1UVKnw9/Zl4oix3DpyPH7ebXuP5hiMBjQ6LdZW1tja2GBjbfOXv9mb\nTCY27d7G+199ApLEsw89ztgho/7we8yZM4fCwkLWrVuHg4PDX+rT1eq6WHIrSZIfUCTLsixJ0k3U\n7zRY1s7dEoRrllqt5oMPPuDnn38mKCiIjz76qM2JbNn5uXyzYQ1b9+2gTq9ncN/+3DXhNmK7927T\n0I9GV0taQRYX8i6SW1ZIYUUJeeVFlFSVU6PVoNPXWbSXJAk3B2f8PXzwc/cmxCuA7qGd6BHaGVeH\nti23tbKyYvLoW4jp3puXP3yTF5a9TkZOJvNmzGlzQLKysuLJJ59kzpw5fPPNN8ybN69Nr7tetffq\nqe+B4YAXUAS8DCgAZFn+TJKkR4H5gAGoBZ6UZflwa+8pnjQEoXlHjx7ltddeo6SkhLvvvpt58+a1\naSjqQmY6X61bza4j+1DY2DBx+FjumngboYEhrb6uvFpFQnoSx1MTOXUxmaySxj0z7BS2+Lt7E+jh\nh4+bJ45KBxyVDjjYKjGajNQZ9Oj0dahqqsgvL6agopi8siKMJiMAHX2DGdG9PxNjRxLo2bYhJ4PR\nyNvLl7F+5xbum3wn/7zvoT/0JLNw4ULi4uLYsmXLVVlP66+6Zoan/m4iaAiCJa1WywcffMCPP/5I\naGgoL7/8Mt27d7/s67Lzc/n0+xXsPLwXR3sHbht7K3dPmt7iXIUsy2QW57IvKY7dZ45wPjcdAEel\nA73DoukSHEGnwDA6BXbEx9XzDw89afU6krPTOJ15jvi0M8SnnUGWZfqGd2P6wPGM7DHgsu8pyzLv\nfPEha7dt4L7Jd/LYzLlt/vykpCRmzZrFggULmDFjxh/q+7VABA1BEEhNTeX5558nIyODGTNm8Oij\nj2JnZ9fqa0oryvniv1+zftcvKGwU3D1pOvdMur3FDOwiVSm/JOzll/i9ZJXkAdAtJIohXWO5KbIH\nnQLD/3LSXUufuyV+D5uP7ya3rJCeoZ15YvIDRAdHtPq6SwPHs3OfZNrNbV7Rz+zZs6murubHH3+8\n7lZSiaAhCDcwWZb56aefWLp0Kc7Ozrz66qv079+/1dfo9XrW/LKOL39chbZOx7QxE5kz/T683Js+\nWRhNRg6dS+Cnw9s4duE0JtlEn45dGd1zEEO6xuLr5vW/urRm+7Lp+G4+2/odFTVVTIwZwf9Nmo2L\ng1PLrzEaefyN50hIOsWXb3xEl/C2FVtcv349S5YsYeXKlW16WruWiKAhCDcorVbL66+/ztatWxk4\ncCCvvPLKZXfIi0s8wbtffEhmXjaD+vTjydmPEBIQ1KSduraGjcd38eOhX8grK8LH1ZOJsSOZGDOC\nIK/mNzKSZZmC6lIulueTWZ5PlqqA0hoVqlo1ldpqNHod1pKEJFlhY2WNl6MrPk6e+Dl70tEjkN6B\nnfB39rrsN3u1VsOKnT/y/f5NBHr6suyBF1rsE4CqqpJ7nnoIB6U9q//1nzYlG9bU1DB27FgmTZrE\nokWLLtv+WiKChiDcgHJzc1m4cCGpqanMnTuXOXPmtLqyqbK6imVff8rmPb8S6OvPk7MfYWjswCbt\nqjRq1hzYxPcHNlOj1dAzrAt3Dp7A8G79mww9ybJMelkux3OTOZV/gdP5F6iorTKfd7d3xsfJAzd7\nZ1yVTtgr7EAGo2zCYDJQUqOiqLqcYnU5dcb6vcR9nTyICYpmfOeB9A3qgpXU8jWdzEhm4VdvAfDe\n7GfpGdbylq4H4o/w5JvP8+i9/+D+qW2bp1iwYAHnz59n8+bN19UQlQgagnCDSUhIYOHChZhMJpYs\nWcKgQYNabb/v+CHe/Ox9VFUqZk6ZwQO334edra1FG7VWw3d717Pm4BZqtBpGdO/PrFHT6RIUbtHO\naDJxIu8c+zNOciDzJAVV9alU/i5e9PKPood/JOGeQYS6B+Bq3/Kw0aVMsomMsjxO5qdwKi+FYzlJ\nVOs0BLv5ckePMdwaPRSlovn5mZzSAh7/4jWKVKW8/8ALxEb2aPFznnrrRY4lJrDx09W4u7pdtl+/\nDVGtWbOGiIjW50+uJSJoCMINZMOGDbz55psEBQWxbNkygoKaDi39pqZWw3tffsTmPb8SGRrOy48u\nolOY5c3PYDSyMW4nn29bTUVNFSO69+eBMXcQFRBm0S5bVciWcwf55fwhitXl2FkriA3uyuCwXgzs\n0ANfZ8+/7Rq1hjr2psfzY+JOkgrT8Xf24smh9zC0Y59m26tqqpj3yQsUV5ax8rF36OAT2Gy7zNxs\n7nh8NrOn3c38ux+4bD9yc3OZMmUKzzzzDNOnT/9L13Q1EUFDEG4AsizzySefsHLlSvr168dbb72F\ns3PLiW/n0i/w3NLF5BcXMmvqDB68fWaTsfyEtCTeW/8f0guz6RnWhSdunWOxIskkmziadYbvT/1K\nXM5ZrCSJ/iE9mNhlMINCe7b47b9Gr6VYU0GJRoXGoMNgMmIwGbCWrHG1c8Rd6YyH0hkPpctlh31O\n5J3n3b3fkFGex8iIWJ4dMRsXZdPcifzyYmZ/8DTODk58+8RSlLbN923hOy8Tn3SSLct/wF5p3+pn\ny7LM2LFjGTBgAK+++mqrba8l10VGuCAILTMYDCxZsoTNmzczdepUFi1a1GJFV1mW+XHbBpZ99Snu\nrm58vvh9enWxXP1Tqanmw81fsyluF/7uPrx539MWuQ8Go4Et5w+y+uQ2MisK8HZ0Z37/6UzoMhhv\nJ8sKsCqdmrOlmVwozyGlPId0VR5qfW2brstd6Uy0ZyhdPTsQ49eZYBefJm36BHZm1V2L+fbkVpYf\n+5kLJVl8OPlpAl0t2wZ4+LDk3gU88vnLLN++hscm3t/sZ86YeBt7jh1gz9ED3DL85lb7J0kSUVFR\nZGRktOl6rjciaAjCNUir1bJw4UIOHz7M3LlzefDBB1v8dq7V6Xjjs3+xdf9OBvftz8v/XISbs6tF\nm12nD/Puz8up1FQzc8RUHhxzp/lbudFkYkfqUf5z7GdyK4vp5N2BxTfPY1RELDbW9bcQWZa5WFnA\nsYJzxBWcI6U8BxkZG8maMDd/hgT1wN/JEx8HN3wc3HFUKLGxskFhZY3BZESlU1OhraaktpKU8myS\nSzM5lHeG5YmbifYMZXLEIAYFdsPaqnHS3cbahlkxk+gT2JkFm9/nH2uX8Nm05whxt8wQj43swZR+\nY1i9byM39xpM59/NxwD06tKdQF9/Nu/dftmgARAaGsqGDRuQZfm6mgxvCxE0BOEao9Vqeeqppzh2\n7BjPPfcc06ZNa7FtaUU5T7/9Ikmp55h712zm3HaPxWqqGm0t7/68nF8S9tI5KJwPH3rZYt7ieM5Z\nlh1YTVpZLlFeIfxr4hMMCu1pvlGqtGr2ZJ9ge1Y8mZWFAES5B3Nv9Bh6+0YS7haArfXll7L6O106\n91E/gV+iUbEv5xS/ZBzjzWPf4efowaO9p9HXzzKnood/JJ/f9jwPr3uTJzcv5cvbX8JVaTnZ/s+J\n97PnzFE+3/Y97z/4QpPPlySJMYNGsGr9GtQ1apwcW5+sDwwMpLa2lsrKStzcLj95fj0RQUMQriFa\nrZYnn3yS48eP8+KLL3Lrrbe22DY9+yKPv/4squoq3n76FUb2H2pxPjk7lRe+W0p+eTEPjrmDOaPv\nMC+fLVFXsOzg9+xMPUagizdLxj7MqMhY81LXDFU+P6bs5UBuIkbZRCf3YB7pPZVBgd1wVzY/p2KS\nZTQGHRq9FkmSUFrborRWYGNl3ey3dW8HN6Z3Gs7UqKEcy09mZdJWXjj4BTeHxjKv163Y2zTOT3T0\nCOTtWx7jkZ/f5tmtH/PhrU+Zn4IAnO0duWfYZD7Z+i3nc9Obfdro3yuWr9atJuHsaYbd1PrKMy+v\n+uTFsrIyETQEQbg6GQwGnnnmGY4fP87LL7/MxIktl784fT6JJ954DjtbO/6zZBmdO1p+O1935Ffe\nW/8FXs7ufPbwa/QKiwbqh5l+TtrDR4d+wGgy8o9+U7mvzy3Y2dQvxT1bmsma87uIL0zB3saWSeED\nGRd2Ex1cG4eEZFmmsKaCi1WF5FSXkqcuo6Cmghq9FpmmC2+U1raEuHgT6uJDBxdfunqGYHfJ04m1\nZMXAwG7E+HXiu+SdrE3ZS2pFLq8MmoWPQ+NcSs+AKJ4bOZtXd/6HH8/sYkavsRafM33gOFbuWsvP\nR7fz7PT5TfrRIyoaW4WCU+fOXDZo/BYoVCpVq+2uRyJoCMI1wGQy8eqrr3Lw4EGeeeaZVgPGoRPH\nWPTuK/h6evPRS+8QcMnGQ3qDnvfWf8HPR7czsHMfXr37cXOZ8dIaFUt2fcmRrERuCu7KohH3E+Ra\nvz9GbnUJK878wpH8s7jaOXJ/13FMCO+Ps2393hJGk5Hk8hzOlGZytiyLcq0aADtrGwKdvOjlHYaz\nrQNOCiWOCiUmZHSGOrRGPeXaarKrStiVfRqjbEJprSDWL4qhgd0Icm4sR2JrrWB29/H08O7IG0e/\n45l9y3l3+Dw87RvnZ8Z3HsT2C/XzLzdH9cfTofGck70jw7v1Y+fpQyyY8iC2NpbDZgqFgo7BoaRm\npV/292FvX7/CSqvVXrbt9UYEDUG4Bnz66ads3bqV+fPnt5obcCD+CAvffZnw4FA+fPFtPFwbv4mr\na2tY9M07HE9N5P6R05g37m7zxPKRrERe3bEcjV7HU8PuY3r3UUiShEav5bvkHWxIO4SttYL7u45j\nSuRglA1PHsUaFYfzz3Gk4DxVdRrsrBV0dg9iXGgMUe6BeNu7YtXGiWK9ycjFykIO55/jaMF5DuSd\nJdY3krs6DcPhkmW8ff068fqQB3nuwHKWHFnFu8PnY9NwHZIk8eTQe7lr9XOsStjC40PutviMsb2H\nsPXEPhLSzjCgc9P8jvCQjsQlXn7Jvm1DEqROp2vTtV1PRNAQhKvc2rVrWblyJVOnTmXOnDkttjuY\ncJRF775CZIdw/v3yuzhfMplbWlXO/33xGhmFObx812NMiBkB1OdcrDi+kf8c+5kIzyA+HfcIYR4B\nAJwousAHCT9RolExNiyWmV3HmucrsqqK2XLxOGdKM5GQ6ObVgUEB0XT1DDHfwP8ohZU1Ue6BRLkH\ncnvUYHZmn2J71gnSKwuZ03UM4W6NdaQ6e4bwf32n89ax1XyXvIP7u40znwtx92NURCybkg/wUL9p\nONg27hnSJ7x+BVZCelKzQcPf25fSinIMBkOLy5eh/skPuCL7nl9tRNAQhKtYXFwc77zzDoMHD2bR\nokUtLu88de4Mi959mYgOYXz80jsWAaO4soz5n75IaVUF7z/wAv079QJAZ6hj8Y7/sDMtjvGdBvLM\nyNkobWzRGfV8mbiFTemHCXL25r3h84n2CgUgt7qUjRnHOFOaiYONHRPCYhkcEI3bJauVZFlGpdNQ\noFFRrlNTrq2hsk6DSa6f0ZBlGWeFEi97Z7yUzgQ4uuHxu9VOjgolk8P7090rlBVJO1h64mdmdx1D\njG+kuc2w4F6cKErlh/N7GBTYnQj3xozv23uMZvuFo+xKi2NSdOMCAKWtHdHBEZy+eK7Zn6O3hxey\nLFNWWYGvp3eLvxeDwQCA9f+g5PvVTgQNQbhKlZaW8sILL9ChQwfeeOONFr/5Zufn8tTbL+Ln7cuH\nL75tse9FaVU5D3/2EuXVKj6e+wrdO3QCoFKr5unNH3C64AL/HHQn9/QejyRJ5KvLWHLkGy5WFjA1\ncgj3dxuHnbUCraGOTRlx7MlJRGmjYFLHfowI7oF9wzCVSTaRWVVKelUxF6tKqKzTACABrnYOuNk6\nYC1ZmYNeZV0t2SVlGOX6b+xBTh708upAlKsf1pd8e+/o6sfz/e7gk9NbWJG0AyeFks4ewebzD/Wc\nyOG8JFYlb+fVQbPNx7v7ReDj5MGhzNMWQQMgwr8DuxOPNPuzdGyYq6itbT0RUaOpv77rdb/w1oig\nIQhXIaPRyAsvvIBGo+Gzzz5r8eZUpa7miTefA2DZc29aJO1Vaqp55PNXKKks58N/vGQOGBW1VTz6\n89tkVRTy+riHGR3ZD4BTxWm8fmQVAIsHzSHWvzMASaVZrD6/F5VOzeDArkwO74+jon7Ip9ZQx+nS\nbE6VZlGt16KwsibEyZNYnzCCnDxxt3NocbjKJJuo0NWQXlnMqdJsNmeexNHGjls69CTUpfFbvr2N\nHY/0nMg78Wv5Imk7r/S/Gyfb+pu7o8KeqZFDWJW8nZyqYnP2uCRJDAjpzs60OIwmk0UgCvbyp1JT\nTZVG3WTPDaVd/XVp61qf4Far6yf6nZzaVnzxenLjDcgJwjVg5cqVxMfH8/TTT9OxY8dm2xiNRp5/\nfwn5xYW8u3Axwf6NwzNavY6nVrxBbmkB/5rznLk8eGWtmkd/foecymLev/VJc8DYlZXAiwe+xNPe\nhY9G/x+x/p0xmoysvXCQf5/ejNJGwVMxt3F35+E4KpRoDXr25CXzWdIuDhSk4G7nyOSwvjzafQzT\nwmPp7R2Kt71zq/MbVpIVnkpnbvIN5x/Rw5keHou9jS1r049zujTboq3SxpYHuo1Fo9ex5eJxi3Pj\nwm4C4GDeGYvj3fwjqKmrpaC6xOK4p3PDctmaKn7P1PDkI7VSeh3q8zMAPD3/voKM1wrxpCEIV5mT\nJ0+yfPlyxo0b12ry3ierv+ToqeM8P38BvaMbS3/Lsswr339AYlYKS+5ZQExEfY2pmrpa/rnhHbJV\nhbw38XFig7sCsCH1IJ+d3khP7wheGHAfTrb2VNVpWJ64jfTKAoYHdWda5CAUVtaYZJnTpVkcLLiA\n1qinq0cQsT4d8bZvuUhiW0iSRJiLDwGO7my6eJLtOWdQ67UM8m/MLwl08mRwYDT7884yIrgHPg71\nN38Pexc6uQdzND+ZGV1GmduHuddP6GeU5ZmXDgM4NRQ2VNfWNOmHXl+/f4eilUlwgOLiYqytrXF3\nd2+13fVIBA1BuIrodDpee+01/P39efbZZ1uc+E44e5pVG35g6piJTBk9weLc9wc2sTvxCP+cMJMx\nveqT1AwmI89v+4S00hzenfg4/UK6AbDt4jE+O72RgQHdeKb/3SisbCisqeDjU5upqtMwp+vNxPrV\nTz5X6jRszjpJfo2KECdPRgRF42PvYvHZJlmmXKehoLaayjodWqOeWqMBvcmIg40tTg1//Byc8FU6\nNbk+O2sF08Jj2JZ9hsOFqfg6uBJxyQ3/lrBYDuYlc6wghUnh/czHe/lE8OOFfehNBhRW9bc1v4ay\n7KWaSovP+G1uSG80NPm5qhvmKhwdmlbMvVROTg5+fn5iIlwQhPb11VdfkZ2dzccff4yjY/M3Lk1t\nLYs/fptAX3+emGWZ2ZyUdYGPNn/DsG79uHf4FPPxDw+u4UhWIotGzGJQaE8A9uec5sOEdcT4dmJR\nv/qAkVtdygcnNwAST/SZTFhDpndaZRG/ZJ1ClmFCh150cQ+wuOGrdLWkVJVSUFuN3mRCAlwUSuxt\nbHC1VaKwskZj0KPW11GsVZNaXYaLwo6ubj4EOVoWT7SSrLg5uBsltVXsyDlDiJMntg0lQdzsHIl0\nD+BEcToTO95k7kOIqy8m2UR+dak5O92xYd6jps5yUtvQECwUNk1rYpVXVgDg4dJ6aZDMzEzCwsJa\nbXO9EkFDEK4ShYWFfP3114wbN47+/fu32O6T1V9QUFLE8teWWez9UKvT8vL3y/Bx9eDFOx4131D3\npMXzw+nt3NnzZqZ1q8/PSCksEiVPAAAgAElEQVTP4b3jPxDt2YHnB9yHrbUNpbVVfHRqEzZW1jze\nZwq+Dm7Iskx8yUX25p3D196VW8P64GbXOClfrddxVlVMTk0lCisrAh1c8Ld3xkfphG0L38KNsomc\nmkpSKks5UpJDB001sV6BFkHIxsqakYHRrEk7ygVVId08GzeV6u4Vyk+ph6iuq8WloS/+jvVPFcUa\nlTlo/BZo9A1bxl76cwKa3fejuKwENxfXVvcLr6urIzMzk4EDm26LeyMQQUMQrhLLli1DkiQeffTR\nFtucPp/Ef7eu547xU5rsh/HRlq/JLSvk03mLzauCiqrLeGP3Crr4hPHPQXcC9XtdvH5kFe5KZ14c\neD9KG1uq62r56ORGDCYjT/WdZg4YBwtSOFqUTpSbH7d06IWiYWLbJMucVRWRUlmKlSTR2dWbTi5e\nLQaKS1lLVoQ6uRPi6EayqphzlSXYWdvQ08OypHmQkwfOCiUXVAUWQcOvod5Uca3KHDR+q1VVZ2oM\nEFpDHQBKG8vgUFpd/zTh7erRpG9Z+TmE+Le86yHAhQsXMBgMdO3a9bLXej0SQUMQrgIJCQns3LmT\nefPm4efn12wbk8nE2//5AF9Pbx6++0GLc2eyUlh7eBszhk6iT3j9fIUsy7y992v0JgOvjZ2PwtoG\nWZb5MOEnVDo1S0c8jKudIybZxJdJ26nQqfm/3pPxd6q/mZ4oyeRoUTo9PUMYE9zN/CRQZzRypCSb\nYm0NoU5udHfzRdnMUM/lWEkS3dx9qTMZuVBVSgcnV9xsG5+cJEkiwtWXM2U5FvtWeDZMulc01LeC\n+kAE9XM3v1E3DEs5KBozwgEKK0pRKuxwUlouY5ZlmYycTIbGtP4EkZiYCHDDBg2x5FYQ2pksy3z0\n0Ud4e3tz7733tthu24FdpGam8+i9/8DB3t7i9e9vXIGnsztzx84wH99/8SSHMk/zj5umEuxWP5l8\nIDeRI/lnua/rzUS413+j3pl9ipSKXO7sNNRcqiOjspg9eclEuPpaBAy1vo7dhRmUaDXEegYS6xX0\npwLGpbq5+WAlSaRXlTc5527niEE2UXvJEJOEZP7Xb6oakgldbBvngfKr6pfa+rs0Fj0ESC/MIsw3\nuMkkfF5RAaqqSqIjOrXa32PHjhESEtJicL/eiaAhCO3s0KFDJCUl8dBDD6FUKpttYzQaWf7DV3QK\ni2DMoBEW5/YmHSUp6wLzx9+Ng119MNEbDby//zsiPIO4s+cYAGoNOj47tYFI90CmRQ4BoEijYmP6\nMXp7hzPQvz6Xo7pOy+ask3jbuzChQ69LnjAM7C3MQGs0MMwvlFDnv2e5qa21DcEOrmTXVCLLlqXT\nfytUqNE3FgbUmxpKeFySsFeurc+5cLNrTLbLUdVvChXo2pgoKMsyqfmZRPiHNOlHYspZALpFdmmx\nr1qtlhMnTtCvX78W21zv2jVoSJK0QpKkYkmSklo4L0mS9KEkSWmSJCVKktS0wpggXONWr16Nj48P\nkyZNarHNvuOHyCsqYM70ey2K5MmyzFe71hHk6WcuQgiwKXk/BdWl/HPQXebNiNanHqRCp2Z+rynm\n6rY/px7GxsqauzoNNQeHXblJGE0mbg3rY55MlmWZ+LJ8tEYjQ31D8Va2viT1j3JR2GGQTZh+FzTq\nGlY62V2yoVJJbX2A8FI2Lve9qCrAWrIi8JJS6meLMnC2cyDgkuzyjMJsVDVV5mTHSx09dRw3F1ci\nQ5tu0GRuc/QotbW1DBs27A9e4fWjvZ80vgLGtXJ+PBDZ8Och4NMr0CdBuGLS0tKIi4vjjjvuaLWq\n6rcb/0ugrz/DYi03Bzp64RTnctOYOXKaORAYTEa+it9Ed78Icz5GjV7LTxf20d8/mi6eHQBIVxVw\nuvQiY0P7mCeUs6vLSK0sYoBfJO52jYEhp6aSPE0V3d198LCz5+/221ODviEj+zdqff1KJ4dLJrML\n1PXDWL6OjctiL1Tk0sHF12Jr2cSCVLr5hpt3GwQ4nlY/H9E33HIRgcFo5OipePr1jGm1cu3u3btx\ncXEhJibmD13f9aRdg4Ysy/uBpgOZjSYD38j1jgJukiT5t9JeEK4pmzZtwsbGhilTprTYJjkthTMp\nycyYeFuTZLIfD/6Ch7MbE/oONx87nHmaInU59/W5pXHZbfZJavRa7uoy0txuT04i9jZ2jAruaT6W\nUHIRBxtbYnwacxBkWSa5sgRXhZKo380P/F2q9DpsJCvsfld2JFddjpfS2WIo6lx5DsHO3uYVUxq9\nljOlGfTwbnxCyKss5mJ5vjlo/mZ34lHCfIMJ9PS1OH488QQVVSpG9hvSYh/1ej0HDhxgyJAhrQb4\n6117P2lcTiCQc8n/cxuOCcI1T6fTsWXLFoYMGdLqPtNrf92AvVLJhGE3WxwvUpVy+PwJJt802iJR\nbWPyPjwdXM1JfABbM44S7hZAlHt9hdhKXQ0nSzIYGNDZ/O28qq6W9MoiengGW9SMKqitplqvo7Or\nV4sZ6n+FLMuU1KrxUjpYvH+d0UBeTQWhlww5VdfVklFZSI+GUu0AxwvPYzAZGRjYGCD2pNdvpDQs\nvK/5WJGqlFMXk7m51+Amfdi6fwfOjk4M6ttyfkxcXBzV1dWMGjWqxTY3gqs9aLSJJEkPSZIUL0lS\nfElJyeVfIAhXgQMHDqBSqbjttttabKOrq2P3kf2MGTgCJ0fLiqq7Eg9jkk1MjG18etDUaTmalcTN\nUf3NcxkF6jIyKgsY1aGv+aZ8tiwbk2yiv19n82szq0qQgWgPy+9lRbVqrCWrJpnbf5dibQ3VhjoC\nHCxLkpwtz8Mom4h0a1yldCg/GRmZPj6NTxXbLsbhZe9q3vNDlmU2nt1PV99wi/mMDcd2IEkSY3tb\nPk2oa9TsPnqA0QOHY9ewI19zNm/ejIuLyw09CQ5Xf9DIA4Iv+X9QwzELsiwvl2U5RpblGG/vljdO\nEYSryZYtW/Dx8SE2NrbFNsdOx1NTq2H0wOFNzu1PiiPcL4Rgr8YR2/jcZPQmA4NDe5mPJRRdACD2\nkgCRUpGLk8KeAKfGKq15NRXY29jiYWcZnFR1WtxslW3etvWPkBuSBO2tbQh1anzaMsom4osz8HNw\nJdCxfpWWwWRkb+4ZungEm/udVVXEqeI0JnQcYM7VOJ6TTJaqgOk9Gp8IDEYjG+N20T+qF0FeliPc\n2w/tRVenY/Ko8S32s6Kigj179jBhwgTs7Jpmkt9IrvagsRGY2bCKqj9QKctyQXt3ShD+KrVazZEj\nRxg7dmyrRe8OnTiGo70Dsd17WxzX1uk4nXmeQV36Whw/mZeCnbWCngGN1WHPlWXhqXQh0KlxmCe7\nqoSOrr4WgaBcp8Zb6dxkCKrOZGwy1/B3SakspUxXS7Sbj/mmDxBXlI6qTsMAv0hzf/bmnqFSV8Po\nkMaAuDp5B3bWCsZ1rC+PLssyK45vwMPBlVERjcF45+lDFFeWMW3AWIvPN5lMrNnyE5Gh4URHdKYl\nP/74IwaDgalTp/4t130ta+8lt98DR4BOkiTlSpL0gCRJ8yRJmtfQ5BcgA0gD/gM83E5dFYS/1bFj\nxzAYDAwdOrTVdglJp+gd3aPJxGtyThpGk5GeoZZLR88VXyTSOwTFJUtUc6qLCXHx/V2BQTUeSsty\n5nqTybzE9lJOCluqG0py/J0y1RWcURUR7OhKmFNjzkexporDhal0dgswV7it0mnYknGcbp4diPas\nz7FIKc9mf24i06KGmvMzjmUncTI/hTmxt2L3266CJhMrd62lo28wQ6Itn+oOxB/hYm4WM6fc1eJ8\njUajYc2aNQwdOrTFvU1uJO26BECW5RmXOS8Dj1yh7gjCFXP06FGcnJzo3r17i21UVZVk5ec0O2xy\nPi8dgK4hkRbH08tyGR15k8WxoppyooIaJ8UNJiNaox5n26ZLZ02/W/IK4KZQkq+pprJOi6tt88mH\nf4Qsy1yoKiWxoggfpaNFsUKNoY4NFxNQWtsyumG/D5Mss+rcbgwmA9OjBpuv4eMTP+Nm58T0qPqc\nCYPRwIeH1uDv4sWUrsPNn7clfg8Xi3JYcs+TFstpjUYjn//wFYG+/s0O//3m+++/p6qqitmzZ7fY\n5kZytQ9PCcJ16dSpU/Tq1avVpZupWfWBISosssm53NICnJQOuDs1Tk7rDHVU6WrwcbIsxKcz6s17\neUN9nSYrSUJ/SZ0mAC+lEyW11U0+K8LFE4WVNafKC5pkbP9RepOR+LI8EiuKCHJwYbBPB/OwlN5k\nZH1GPGq9lqkd+5r7vD3rBEllWUyPGoxvw8ZLa87tIk2Vx6N9ppprS60+9SvpZbk8OeQe85NWdW0N\nH2/5hu4dOjG6p2WOy8bdW81lWWxaGCJUqVR8/fXXDB8+vNUAfyMRQUMQrrDq6mouXrxIjx49Wm2X\nnp0JQESHpvs25JcXE+jpZznk1HDDd//dxkh6k9FivkCSJOytbanRW+6D7evgSrVeS4XOckc7O2sb\nurn5UKyt4UxF0Z8KHEbZRGpVGVtzL5CpVtHF1Zv+3sHm/AuNXsePacfIq6nglg69CGiY/D5ZnM7G\n9GPE+EYytGFJ7aniNL4/v5uRIX0YFFh/I08vy+WLuPUMDevD0I6NhSM+2vw1Kk01T0/9h8VThqq6\nkk9Xr6BXl+6MGtBydve3336LRqNh3rx5Lba50dy4GSqC0E7S0+ufIDp1ar0wXmFpEXa2dni4Nq3x\nVF1bg6tD81us/n6Vk6utA5V1loEg0NmLrKpii2Nd3APYn3+ehOKLjA62TIoLd/agSq8jpaqUirpa\nurr54NWGUiJqvY6cmkoyqivQGPV4Kx3p4e6LxyV7cpTWVrMu4zg1eh2TQnvTyb1+ddO58hxWJG0n\nzNWXe7uMQJIk8tVlvHH0W4KcvXm4d31CpKZOy3NbP8bR1p5nRtxvft8jKSdZf2wH9w6fQucgy9Ig\nn3z3JVXqKhY++FiLcxkFBQWsXr2a8ePHExERcdlrvVG0GjQkSfIDkGW5UJIkb2AIkCLL8tkr0TlB\nuB6lpaUBEB7eco0jgOKyUnw8mk+o0+hq8XC2TAj8rYzIbwX9fuPl4EqxRmVxrKOrH9uzTlCj1+LY\nMLzjpFDS1T2QM2U59PUJsygjIkkSvT38cVHYkawqZk/hRdxslXjZOeBiq0RpZYOttTU6o4Eagx6N\noY5SnQZVXf3TjJedA329Aiy2eJVlmbPleezKPVtf/ypyAP4NpUFOlWSwImkHvo7uPNxzAnbWClRa\nNa8cWgnAywNn4ahQYpJNvL7rS7JVhXw0eSGeDa+vUFfy+n//TZhvsEXlX4BT586wfucWZkyc3mqd\nqd/2N3nkETGteqkWg4YkSXOBZ+r/Kb0NzAKSgDclSXpHluUvr0wXBeH6kp+fj0KhwMfHp9V2NZoa\nnJ2cmj1nZWWF/LtJazelE9aSFSXqCovj4W6BHMhNxCibzMNUfX0i2JaZwP7cJMaHNdZRGugfRWpl\nEeszErgzsj8Ol8yFSJJEhIsnoU7uXFRXkK1WkaVWNakXBWAjWeFmq6SHuy/Bjq4W7wNQUlvNztwk\nctXlBDi6MSm0Dy629siyzI7sk6xPO0IHFx/m95yAo0KJSqfm2f3LKdZUsGTwg+Y8jU+PrGVnWhyP\nDryDmOBooD4n47lV71FZU817s5/FTtH42VXqal5c9joBPn48dOf9tGT79u3s2rWLhx9++IYtgd6S\n1p40HgW6AvZAFhDR8MThDuwBRNAQhD+hqKgIHx+fVgvjQf22pHa2zSeSKRV2aOssl8HaWNvg5+xJ\nTqXlsFNPnwi2XYwjtTyXzg3LVYOcvejqGcLunNMMC+puLkHuYmvPpNA+/JxxnB9SjzAptA9e9pbD\nYDZWVkS6eBLp4olJNqEzGtEZDega8jkcbGxRWFk1+4RUoashofgip0qzsbO2YWxwd7p71u9toTXU\nsSZlP8cKU+jrE8HM6FHYWttQVFPBS4dWUKgu45VBs+nmXT/H82PiTr5J2MLUbiO4t88tQP3Tyweb\nVpKQnsQrM/7PYlhKlmVe//RflFSU8eXrH+Jo79CkfwClpaW8/fbbdO3alZkzZ7b4+7lRtRY09LIs\nawCNJEnpsiwXAsiyXCFJ0l9bQiEIN7CKigo8PT0v37AV7k6u5JY2zXON9ArhbFG6xU53fX2jsLWy\n4dfMOHPQAJjUsR/vxP/EqnO7eaj7OHP7UBcvbguPZWPmSb5JOchAv0j6+oSZt3q9lJVkhb2NFfat\nbMRkkk1kVpVysjSTjKoSrJDo7hnMkIBO5ieQ5LJsvju/hwqtmglhsdwSFouVJHG2NJPXjnyNwWRk\n8eA59PSpn1v47+kd/Gv/twwJ681Tw+4z9/27fRv44eAWZgydxC2XFHEEWPnTd+w+up/HZs6lawt7\nZsiyzOLFi9Fqtbz66qs3dGHClrT2E5ElSVLIsqwHJvx2UJIkJWLVlSD8aWq1utUChb+xs7Wjuqbp\nEliAIE8/jl04ZREcAGKCo9mbkUBuZbF5tz5nWwdGdejLzqwE7u86Djdl/ZBXBxcfpoT3Z13aYXbn\nnGbUJZnWIc5ezO48lJ25SRwoSCGuOJ1o90A6ufvj7+BmUdDw90yyTFVdLXk15VysKuFiVQlaox5H\nGzsG+kXS0ysEp4Z5lCqdhp/Tj3C04Dx+Du48FXMbHV39kGWZXy8e5+OT6/C2d+PVQbMJdvFBlmW+\nTtjMp0fWMjSsD2+Mf8Tcly3xe/hw89eM7jmI/5s4y6JPhxKO8tmalYwdMop7b72jxb6vW7eOw4cP\n8/TTTxMaGtryL+cG1lrQmAogSVK0LMvJlxz3AJ76n/ZKEK5jOp2uTfWLnB2dyCvMb/ZcmG8wOn0d\nWSV5hPoEmY8P6tCT91jFrymHebBfY8mLaVFD+DXzOF+e2cKC2DvNx0eH9CK9soC1qYfQGQ2MC+1j\n3n/CUWHH5LC+5KjLOF2aTWJZDidLs4D6FVludg44KZTIsoxJNmGQTah0Gip0NRgb5jkcbGyJcPUl\n3NWHcBdf8xLbSl0NO7JOsj/vLEbZxLjQvtwSFovCyhqVVs2HJ37iSP5ZenqH81z/e3Gxc6TOqOet\n3V+x5fxBxkYN4KXRDzZuMHVsB2+u/ZTYyB68MuP/LIb+ElPO8sy/FhMZGs4L8xe0uFrqzJkzvPfe\ne/Tr14/bb7/9sr+fG1WLQUOW5WwASZL+K0nSKuAdQAk8C8QAO65IDwXhBuXn5cO+uIOYTKYm8x99\nwuuzpRPSkiyCRoCrNwNCuvPz2b3MiplkvqkGOftwR6fhrDm/m5v8uzAkqD5HRJIkHug2llXJu9iU\ncYwLFbnMjB5lUWIk2MmTYCdP6owGctRlFNSoqNBpUNXVUK5VYyVJWElWWEtWuNjaE+rihaedEz4O\nrvjau1jcpPPVZRzIO8uh/GQMJhM3+UVxS1gMPg5uyLLMwdwzfHRiHRqDlge6T2Bq1BCsJSvKalQ8\ns/VjEgtS+cdNU3jgpinm911zYBNLN6xgQKfevD1rEbaXDJWlZV/kiTeew8vdgw+efwulXfMZ7cXF\nxTz11FP4+Pjw+uuvX3a+6UbWlgG7fsDbwGHAGfgOGNTqKwRBaJG1tTUGg+Gy7QJ8/KjT6ykuK8HP\n23LToCBPP/zcvDmQfJzbBlpufnlHzzE8sWkpPyXt5s6ejXtw3BM9hpPFabwXtwYnhT29feszzRVW\n1szuOobOHkGsSTnAK0dWMyK4B4MCuuDj0DiMZmttQ7irL+Guln25nHJtNYklFzlakEJWdTHWkhWx\nflGMD+1rfv+0ily+SPyF0yVpRLgF8lTsnXRwrV+1tC/jBG/sWkGtQccb4x5hVEOZFIPRyAebVvLD\nwS0M79aPJfcuaBIwHn55AXa2tnz80rt4uXs07Rz1taUef/xxamtr+fe//92mocMbWVuChh6opX4V\nlRK4KP9+rZ8gCG3m5ORETU3NZdt1Ca+vVJucltIkaEiSxPi+w/h69zqKK8vwcW2cWB/QoQf9Qrrx\n2ZGfGBkei3dDMUAbK2teGTSLZ/cv55VDK3l50Cz6+EaZ329gQDSd3IP4Oe0IO7JOsj3rBB1d/ejn\n14kIN398HdwtdtBrjkk2UVpbRZ66jFRVPufKcijU1C8BDnLy4vbIwcT6RZnrXuWry/gueQe7s0/g\nYuvAvJ63MiF8ADZW1tTqdbx/4Ds2nN1HlFcIr46dR8eGvT7UWg0vfPsvDp8/wYyhk3hs4v3mPBWA\n1Mx0Hn71KRQ2Nnz66lICfZvf8NNoNPL888+TlpbG+++/L5L42qAtQeM4sAGIBbyAzyRJuk2WZTHo\nJwh/gqurK1lZWZdtFxkajq1CwanzZxg5oGk13EmxI1m5ay1rD23l4VvuNR+XJImnh83k7tXPs2TX\nlyyd9KT5Zu9m58RbQ+fy7P7lvHxwJbd1GsZdnUeibFjF5GnvwoPdx6LSqjlWeIGjBef5PmUfUJ97\n4WHvgo+9K3bWCuysFVhJEhqDjhq9jhp9LUWaSnNyocLKmgi3AAYFRhN9yR4YsiyTVHqR9akHOJx3\nFoWVNXd2GsHtnYfjqKjP1diddpwPD66hsLqMmX0m8FD/aeZ6Uil5GTz/7b/IKyvi2enzmdrfckfD\nE2dP8/Q7L6G0tePTV5cSEhBEc0wmE2+88QYHDhxg4cKFDBokBlDaoi1B4wFZluMb/l0ATJYk6b7/\nYZ8E4brm7+/P4cOHm6x8+j1bhS039ejLvrhDPDHr4SZtg7z8ubn3EL4/sInpg8ZbPG0Eu/ny5NB7\neGvPV7y37xsWDr/f/HpXO0feHjaXz09v4ofzu9mXc4r5vSZzk3/jMlQ3pRNjQ/twc4feFNZUkF1d\nQp66jJLaSsq11eiMeuqMBoyyEQcbJQ4KO9yVznRyDyLAyZNAJ0/8HT0sSq2Xa6vZl3OK3VknSFPl\n4aSw545Ow5kUMRBP+/rCixlluSzd/x3Hc5OJ8Azi02nP0juwvtyKLMusObCZj7d8g6ujM/+e+wp9\nwi3Lnfx6YBevfvwOAT5+fPDCWy0+YZhMJt588002bNjAgw8+yB13tLyiSrB02aBxScC49Niq/013\nBOH6FxwcjE6no7i4GF/f1ucHRvQbwsGEo5xNPU+3qKa5BfPH3cPuxCN8tPkbXrvnCYtzU7uNIK+y\nmFUnfkGSrHhy6D3m5anOtg48FXsnYzr05d8n1/PyoZV08ejA6NC+DA3qiVPD8JEkSfg7eeDv1Px8\nQGtkWSarqoiTRanEFZzjdHEaJmTC3QJ4pPcURneIMT/hZKsK+er4JralHMbR1p6nht3H1G4jzP0t\nrCjh9R8/4diFUwyJjuXFOx/FzbGxMKPBaOTzNSv5at1qekf34N2Fi3F1dmm2X0ajkTfffJP169fz\nwAMPMHfu3D98bTcykbkiCFdYVFTDXEVy8mWDxsgBQ1n61Sd8v3ktrz/5YpPzgZ6+zBk1neXb19C/\nUy8mxIywOP/IwDuQgW9P/EKOqpDXxz2CyyWFBnv6RPDvMY+zOf0I2y7G8dGJdXx2aiN9fCPp4tmB\nTh4hRLoHmetTtcQom6jQVpNVWUi6Kp+MynySSi5Spq2q76eTF3d2HsnwkF6EuDRec3pZLt8kbGH7\nhSMorBXc3nM0s2Nuxa0hC91gNLB6/ya+2P4DAAunPsRtA8dZPHWVV1bw/PtLiD9zkimjJ7DwwcdQ\nKJpPNtTr9bz00kvs2LGDOXPmMG/evFaf9oSmRNAQhCusc+fO2NjYkJiYyIgRI1pt6+TgyNTRE/h+\n81rmF84hyC+wSZvZo6dzPC2Rd9YtJyogjMiAUPM5SZL456A76eDuz9t7vmLmmpf45+C7GBkeY75Z\nKqxsmBo5hCkRg0mtyGVnVgInii5wrOCc+X0cFUrc7JxwUtjj0FAo0Cib0Bn1lGurqdBWW2zg5OPg\nTrRXKL19IunjG4mvY+OTit5oYE96PD+d2cWp/AsobWy5q9dY7u1zC54OjfuDxKUmsnT9F2QU5TC0\n600smPwA/h6W9bqOnjrO4n+/U19T6pGnuXVky/t819TU8Oyzz3L48GEee+wxUSLkT5L+6qYqV5uY\nmBg5Pr7JiJogXFXmzp1LRUUF//3vfy/btqS8lGmPzqR/z768u+i1ZtsUqUp54KNn0Bv0fPbwEsJ8\ng5u0OVOYxpu7V5Jelkt3vwgeG3wXPfybbvD0m+o6DRfKc0hT5VFWW4VKp0ZdV4vWUGfOzbC1tsHD\n3gVPpQue9q6EuPgQ5uqPs61lXSeDycjJvPPsTI1jT3o8lVo1gS7eTO0+klu7DMXVvrEw49nsC3y6\ndTVxqafxd/dhwZQHGNrVcjfCWm0tH65aztptGwgL6sBrjz9Pp7CWVz4VFRXxxBNPkJ6ezqJFi5g2\nbVqLbW9UkiQlyLIcc9l2ImgIwpW3Zs0a3nvvPdauXdumchUr163mk+++YNlzbzCob/9m22SV5DH3\nkxcA+NfsZ+kaEtWkjdFkYsu5A3x+bB2lNSo6+4QyscsQRkfe1GTzpr+qtEZFXHYSx3LOciz7DBW1\n1dgr7Bgc2otbOg+if4fu5uxzgOTsVFbsWsv+s3G4Obowa9Rt3DZgnEWVWqjfN/2Nz5aSXZDL3ROn\nM//uB1C2kmGfnJzMggUL0Gg0vPnmmwwcOPBvvc7rhQgagnAVKykpYeLEicyYMYPHH3/8su3r9HXc\n+9RcKqurWPXuZ/h4ejfb7mJRDo9/8RrFleU8PP4e7hk2udns5lq9jk3J+9lwdi9pZbkAdPeLoId/\nJD39I+nsE4qX4+XzMqD+KaJYXU5eZTEpJVmcL87kfEkmOaoioH4nwdjgaEaExzCwQw+UisYbvCzL\nHD5/gm/3richPQknpQP3DJvMXUMm4ai03MO8TFXOB19/xtb9O/n/9u47PKoq4eP49yST3nslBAIJ\nTQgtAaRIB0UQRRdFVhfQFdauu1bW14Kia2FdyyplUVlExELvTTqEBAgpJCSEdNL7ZJKZOe8fCSwC\nmUwoCQnn8zx5mGTOvaxp5nMAACAASURBVHNugPnNuaf5e/sxd86L9Lutd4P1klKyevVqPvjgA9zd\n3VmwYIGah2GCCg1Fucm9+uqr7Nu3jw0bNuDg0PgueGcyz/LIS7Pp1L4jX735SYOdvWVVFbz74xfs\niD1AZGg4r06ZfVlfwHlSSk4XZrA7NZpD6bEknEu7MM/C0sISb0c3vB3csbOywcrSEitLK/QGA5W1\nWqpqqimpLievvOjCWlMAPo7udPEO5jbfTkQG9aCTZ7vftSgACsqKWHdkJ2uPbCejIAdvFw8eGjaR\nSRGjLwuLmtoaVm1ew8IfvkFXU8Mf75nKo/c+ZLJ1odVqef/991m3bh2RkZHMmzdPzfRuhAoNRbnJ\nxcXF8cgjj/DEE08wa9Yss47Zsm8nr338NiMHDuWd5+aisbzyarNSSn45uIVPVi/BII1MjBjJn0ZO\nwcfV0+T5dfoaEvPSSCnKIresgNzyQvIri9Hpa6g16Kkx6LGy1GBvZYuDtR3Otvb4OXni7+yFn7Mn\nnTzbNXibq0JbyZ74I2w7vp/9iUcxGI307tiNeyJHMzp88IV1ss4zGo1s3beTL5YvITsvhwG9+vHC\nzCcJDgi64vnPi4uLY+7cuWRkZDBr1ixmzZqFZQO/J+V/VGgoSivwwgsvEBUVxerVq83+JLx87Y98\nsvRLhvW/nXnPz8XG2rrBsudKCli6/SdWH96GAO7sdwcT+o2gZ3CXGz7UVEpJWl4mh5NPcDAxhkPJ\nx9Eb9Hi7eDC291AmRY4iyMv/suOMRiO7Du1l8U/LSDpzms7BITw9/XEGhPc3+Xp6vZ6lS5eycOFC\nPD09efPNN+nXr9H3QKWeCg1FaQVSU1OZOnUqEydO5PXXXzf7uJUbfuEfi/9F985dmP/CG5etTXWp\n7KI8lm5fxeaYPWhrqvF19WJIt34MCOvNbcFhv5sod7Vq9bUkZ6cRl5HMyfRkjp6OJa+0EKhbYPGO\nHpGM6DmIbu06XbGfpba2lq37d7H05+WcyTxLkF8gM6c8zLihoxpddTYxMZG33nqLpKQkxo4dy0sv\nvYSz8/Xt2G/rVGgoSivx6aef8u2337JgwQIGDx5s9nE7Dv7GW599gEaj4f+eepnBDYyquliVTsuu\nk4fYemwvMalxVOmqAejgHUhH3yACPHwI8PDFy9kdexs7HGztsLW2wWg0UmvQozcYqNBWUlRRSnFF\nKfmlhaQXZHM2L5vMwlwMRgMA7k6u9AruwoCw3kR07kWAR8Ohdq4wn1+2rOPXbespLCkiJKgDM+6b\nxsiBwxq9rVRVVcXChQtZvnw5rq6uvPTSS4wYMcLs36HyPyo0FKWVqKmpYfr06ZSUlPDtt982Okv8\nYmezM3j5H29yOj2VUQOH8dQf/4y/t69Zx1bX6kjISOH4mQROpCWSXpBNdlEeekPjy7afZ62xop2n\nH0Fe/gR5+RMW0JEeQaH4uHqavP2l1+vZF3OIdTs2sSfqAEYpub1PJFPGTWJgeP9GWxZSSrZu3cqC\nBQvIy8tj0qRJPPPMM6p1cQ1UaChKK5KSksKMGTMICAhg4cKFZo2mOq+mtoZvf/2BpT8vR0ojU++6\nj4funoKHa9PXizIYDeSXFlFQVkxVTTVVOi3VNdVYCAusNBo0FhocbO1xd3LBzdEFZztHs/tGDAYD\nsUnx7Dy4h017t1NUUoy7qxt3DRvDvWPuJtD38v6NK4mLi+Of//wn0dHRhIWF8be//Y1evXo1+VqV\n31OhoSitzIEDB3j22Wfp27cvH330EXZ2do0fdJFzhfl8tmwhm/dsx0qjYczgEdw/7h66hoS22PpK\nldoqok7GsO/oIXYf3kdRaTEajYbBfSK5e8R4BvWOQKMxbzWjjIwMvvjiC7Zu3YqrqytPPPEEkydP\nViOjrpNWERpCiHHAPwFLYJGUcv4lzz8K/APIqv/RZ1LKRabOqUJDac3Wr1/Pm2++SY8ePfjkk09w\ncXFp/KBLpGdnsnzdKjbs3oK2upp2vgEMi7idYRG3c1totxv6JltWUU5cciInkxOIOhnDiVNx6PV6\n7GxtGdx3IHdEDmZQ7wgc7c1vSWVmZrJ06VLWrl2LlZUV06ZNY/r06Tg6OjZ+sGK2mz40hBCWQBIw\nGsikbrOnB6WU8ReVeRToJ6V80tzzqtBQWrudO3fy6quv0q5dOxYsWIC/v3m3bS5VUVnBln072Xlo\nL1EnY9Dr9TjY2dO1Uxg9OnclNLgTHQKDCPDxw87W/FaNlJLyygpy8nLJyM3i9NkzpGScITU9jfSc\nutnlQgg6B4cwoFc/Bob3p2eX7lhbNTw0+EpSU1P5z3/+w+bNm9FoNEyaNImZM2fi6Wl6rolydVpD\naAwE/k9KObb++1cApJTvXVTmUVRoKLegqKgoXnjhBaBu5vjYsWOv6XwVVZXsjzlMTPwJ4pITSEpL\nwWAwXHjew9UddxdXnBydcLJ3xNbGhvPvDUajkQptFRWVFZRXVlBQXEilturCsRYWFrTzDaBDu/Z0\nCwmje+eudOsU1qTWxHlSSo4cOcKyZcvYv38/tra23H///UybNk2FxQ3WGkJjCjBOSjmr/vvpQOTF\nAVEfGu8B+dS1Sp6TUmaYOq8KDaWtyMrKYu7cuZw4cYIJEybw4osvXrdbMrqaGs5kppGenUlmbjZZ\n53IoKS+lvLKCisoKqnW6un4QAQKBo4MDTvaOONo74OHmgZ+XD35ePvj7+BIc0N7kkh7mqKioYOPG\njfz8888kJyfj7u7OAw88wJQpU9TyH82krYSGB1AhpdQJIf4M/EFKedkgbCHE48DjAEFBQX3N2X9Z\nUVoDvV7PokWLWLJkCa6ursyaNYvJkyc3uO5UayKlJC4ujjVr1rBx40a0Wi2hoaFMnTqVsWPHYnON\nQaQ0TWsIjUZvT11S3hIoklKa7BlULQ2lLYqPj2fBggVER0cTEBDA7NmzGTVqlNkjj24WUkpSUlLY\nvHkzW7ZsISsrCxsbG8aMGcN9991Hjx49Gj+JckO0htDQUHfLaSR1o6OOAA9JKeMuKuMnpcypfzwZ\neElKaXLaqwoNpa2SUrJ//34+//xzkpKS8PX1ZcqUKUyaNAk3N7eWrl6DpJQkJSWxfft2tm/fztmz\nZ7GwsKB///6MGTOG4cOHq0l5N4GbPjQAhBB3AguoG3K7REo5TwjxFhAlpVwjhHgPmAjogSJgtpQy\n0dQ5VWgobZ3RaGT37t38+OOPHD58GCsrKwYPHszw4cMZMmQITk5OLVo/KSXZ2dlER0dz+PBhDh8+\nTGFhIZaWlvTt25cRI0YwfPhwPDw8WrSeyu+1itC4EVRoKLeS1NRUfvrpJ7Zv305BQQGWlpb07NmT\nyMhI+vTpQ/fu3W9434BWqyU5OZnExERiYmKIiYmhoKAAAHd3d/r3709ERARDhw69qVtEtzoVGopy\nCzEajcTHx7Njxw6OHDlCYmIiUkqsrKzo0KEDISEhdOzYkXbt2uHl5YW3tzceHh5Ym1hW/Ty9Xk9R\nURGFhYUUFhaSk5NDZmYmmZmZpKWlkZ6efmF4rre3N3369CE8PJzw8HBCQkJabDa60jQqNBTlFlZW\nVkZMTAzHjh0jJSWFlJQUzp07d1k5Kysr7O3tcXBw+N1Mcb1eT1VVFVqtlpqamsuOs7GxISAggKCg\nIEJDQwkLCyMsLAwfHx8VEq2UuaHRuoZeKK1CRUUFeXl5lJSUUFZWRklJCVqtFp1Oh06no7a29kJZ\nIQTW1tbY2tpia2uLvb09rq6uuLq64ubmhpeXV5sYXtrcnJ2dGTZsGMOGDbvws4qKCnJycsjLyyM/\nP5/CwkKqqqqoqqqisrISo/F/W7ZaWlpib2+PnZ0d9vb2uLm54eHhgbu7O76+vnh6eja6Eq3SNqnQ\nUK5KaWkpycnJpKenk5GRQXp6OpmZmeTm5lJZWWnyWI1Gc+HTqNFo/N3M5EsJIfD09MTPz+/CJ9t2\n7doRHBxMhw4d1Fj+JnB0dKRz58507ty5pauitGIqNJRG6XQ64uPjiYmJITY2lqSkpN/d6rCysiIw\nMJB27drRt29ffH198fb2xs3NDRcXF1xcXHBwcMDGxgZra+vLbl8YDAZ0Oh3V1dVUVFRQWlpKSUkJ\nhYWFnDt3jtzcXHJzc4mOjmbjxo0XjrO0tCQ4OJjQ0FB69epFeHg4HTt2VJ+AFeUGUqGhXEZKSWpq\nKvv372ffvn2cOHHiwn3tDh060KdPHzp37kxoaCjBwcF4eXld08qp52+F2Nvb4+5ueg+I6upqsrKy\nOHPmDElJSSQlJXHkyJELYeLi4kLfvn0ZNGgQgwYNwtvb+6rrpSjK5VRHuALUBUVCQgJbtmxh27Zt\n5ObmAhASEsKAAQPo3bs34eHhN+U6QFJKsrKyiImJITo6moMHD5Kfnw9A586dGTVqFGPGjKFdu3Yt\nXFNFuXmp0VOKWfLy8li9ejXr168nMzMTjUbDwIEDGTZsGAMGDMDX17ytQ28mUkpOnz7NgQMH2L17\nN8ePHwegW7du3HXXXdx5550tPgFOUW42KjSUBkkpOXr0KCtXrmT37t0YDAYiIiIYO3Zsm1zSITc3\nl23btrFx40ZOnTqFra0t48eP54EHHlCdwopST4WGchkpJQcOHGDRokWcOHECFxcXJk6cyH333Udg\nYOB1fZ3CkiLyiwopKK77qqisQKurRltdTU1tTX1nuMDCQmBna4ejvQMOdg64Ojvj4+mNj4cX7i5u\n171TOyEhgVWrVrFp0yZ0Oh1Dhgxh1qxZdO/e/bq+jqK0Nio0lN+JjY3l448/JjY2Fh8fHx599FEm\nTpx4zUNW84sKSEhJIiEliZT0M2TkZJGRm4WuRnfF8jbWNthYWyOlREqJ0WhEq6vmSv8ONRoN7XwD\naB/Qjvb+7Qjr0JnbQrvi4+l9zRPIysrKWLlyJcuXL6esrIyhQ4fy7LPPEhQUdE3nVZTWSoWGAkBh\nYSGfffYZa9euxcvLi8cee4y77777qifMFZUWc/hENAePRXEkNpq8wroO5/O7t7XzCyDIL5BAvwC8\n3T3xdHPHw9UdJ0cnbK1trjjK6nxwVFRVUlJawrnCfM4V5pObf470nEzOZmWQnpN5YT6Hh6s7vbr0\nYFDvCAb27o+3h9dV/34qKyv54Ycf+Oabb6ipqWH69OnMmDEDW1vbqz6norRGKjRucVJK1qxZwyef\nfEJ1dTXTpk1j5syZ2NvbN/lcBcWFbNu/m637dnLiVN3K9S6OzvS/rTe9ut5G15BQQoNDmrTPdFPV\n1taSfDaVk8kJnEyK5+jJY+QV1S2KF9qhE6MH3cG4ISPx9fK5qvMXFBTw6aefsmHDBvz8/HjjjTfo\n16/R/z+K0mao0LiFlZeXM2/ePLZt20bfvn155ZVXCA4ObtI5jEYj+6IPsXLDLxw6cRQpJZ3bd2Tk\noGEMDI8grEOna5qbca2klKSkn2Ff9CF+O7KfE6fiEELQp1svJo4cz+hBd1xVayo6Opp33nmHzMxM\nHnvsMWbMmNGi16kozUWFxi0qPj6el156iby8PGbPns306dOb9KZXU1vDr1vX8/36n8jMzcbb3ZOJ\nI+9kzODhdAhsf9X1klKirammoroKXW0NFkJgISywsLDA2c4RW2uba+qnyMzNZtOe7WzYtYWM3Cw8\n3Tx4aMJ93Dd2EvZ2TWsBVVVVMX/+fDZs2ED//v2ZP38+Li4mN4xUlFZPhcYtaO/evbz88su4ubnx\n3nvvNWnrTCkl2/bv5rNlC8nOy6FnWHem3nUfwyMHm72lqJSSzMJcEjJOk5R9hqzCc+QU55FTlEdp\nVQVGaWzwWBuNNS4OTni7ehDk6U977wCCvQPoERSGl4vpWeIXMxqNHDwexbLVKzkSG42bsytPPDiD\nSSPHNyk8pZSsXbuW9957j3bt2vH555/j5XX1fSeKcrNToXGLWbNmDfPmzSM0NJQFCxY0aVe0hJQk\nPlj4T04mJ9C5fUee/uOfGRDe36xjs4vy2BsfxYFT0ZxIS6RcW7dYocZSg5+bF/7u3vi5eePh5IqD\nrT0OtvbYWFljNBoxSiMGo5FybSUllWWUVJSRW5JPen42eaWFF17Dz82LnsFdiQztxbAekTjZOZhV\ntxOn4vjsu6+JSYgltEMn/jbraXp1adoe1FFRUTz//PO4ubnx5Zdf4u/v36TjFaW1UKFxC9m2bRuv\nvPIKkZGRfPDBB2Z3dhuNRv679kc+/+8i3Jxdmf3QDO4aNqbRT+QV2ko2Rf/G6sPbOJWVCkA7Tz/6\ndbqNru060TUwhI4+7bDSXP2S5lU6LWfOZXAi7RTHzyRwLC2BovISNJYaBoSGM7bPUEbcNqDR1zjf\ngvrnt/8mrzCfmVOmM+v+pt2yO3nyJE8//TSenp4sWbIER0fHq74uRblZqdC4RRw7dow5c+bQpUsX\nvvjiC7OHihaVFvN//3qfAzGHGR45hNfnvIizo+mlNVJy01n+2xq2xuylulZHl4COjO09lMHd+9He\nK+B6XE6DpJTEZySz7fh+th3fx7mSAjyc3Jhy+3juGzgWVwfTs9grtVX8Y9G/WL9rMxE9+/L2s6/i\n7mL+1qNHjhzhySefZNCgQXz44Yeqc1xpc1Ro3AIKCgqYOnUqzs7OLFmyxOzFBLPzcpn9xgsUFBfw\n/J/+wr1j7jbZCV1cUcpXm7/n14NbsbGyZmzvIUweOJaugSHX61KaxGg0cij5OCt+W8uBUzE42Ngx\na8wf+MPgCWgaeTNfs30jHyz6J97uXnzx5kf4epq/Cu7KlSv54IMPeOaZZ5g+ffq1Xoai3FRUaLRx\nRqORp556imPHjrFs2TI6dOhg1nHp2ZnMefNFtNVaPn19Pt07dzX5Gj8d2MS/Ny2nSqdlyqDxzBz9\nQKOf6s8zGI2kl+SSXJDOuYoiSrTlFGvLqKqpRmNhicbCEmtLKzwcXPFz9sTfyZNgd388HcxfSTcl\nN53P1n3LvsSjhAV05JUps+nWrpPJY06ciuPpd17GxdGZr99ZgI+ZkwOllLz44oscOHCA5cuXN3kY\ns6LczFRotHE///wz7777Lq+++ir33nuvWccUlRYz/a9PUFNbw2d//wdhHRp+c9Xqqpm7/GN+iztC\nZGgvnps4k46+jS8tXlhVyo7kw+xKjSbuXAra2v8tJ2JtaYWrnRMO1rYYjEZqjXpq9LUUa8swXvTv\nMNjNn8igHgwI6kFEUA80FqZbD1JKdpw4wEerF1FUXsozdz/Kg0PvNnlM3OlE/vLmX/Hx8OKb97/A\n1sa823oFBQX84Q9/ICQkhK+//tqsYxSlNVCh0YZVV1czefJk/P39WbRokVnzGwwGA0+/8zLHEk6w\n+N1/0aVjaINlC8qKeWHJu5zKSuW5iX/igcF3mXwNg9HIrtSj/HpyJ1GZ8RilpKN7AP0Cu9HVO5hO\nnkEEuHhhb2V7xfPoDXryKorJLsvnVP5ZDmfEEZOViM5Qi7ejO5N73MGkbsPwaKQFUqGt5J2Vn7Mj\n9gCPjLiXOeMfNlnvQ8ejeOrtl7hn1F28+sTzJs99seXLl/Pxxx+zaNEiwsPDzT5OUW5mKjTasBUr\nVvDhhx/y1Vdf0bdvX7OO+eaX7/ls2UJem/0C94y6q8FyheUlzPzXyxSVlzDv4RcY0t300NuYrETe\n3/UNZ4qy8Xf2YmzoAMaEDqCjx7WtmqvT13DwbCw/ndzBofSTWFpY8kDPUcwZdD/Wlg2PmDIYDfzj\nl4X8fGAzDw69m+cmzjD5Ov/67mu+/XUFH78yjyH9BppVN61Wy4QJEwgPD+ejjz5q0nUpys3K3NBQ\n2722MlJKfvnlF3r06GF2YBSVFrNk1TKG9b+dSSPvbLCcwWjg9WUfUVhWzJez36ZH+4ZbIzp9Df/c\n+z0/xe7A39mLeePmMDykP5bXaSlzG401w0L6MiykL+kluSyL3sD3xzYTnZXIO+PmEOR65c2hLC0s\neeneP2NpYcn3v62lZ/sujOw1qMHXmf3gDHYd3suX3y9hcN8BZrXa7OzsmDBhAitXrqSyshIHB/Pm\njShKW3B9NytQbriTJ0+SkpLCpEmTzD5myapl6Gp0PDn9MZNviou3/sjRlJP87d7HTQZGha6K2T/P\n56fYHTwYPpblD81jVOfI6xYYlwpy9eXVETP44K5nyCkr4JEVb7DnTEyD5YUQPDfxT/QICuWdHz/n\nXElBg2U1Gg0z7nuY5LQUfovab3adhg4dSm1tLQcPHmzStShKa6dCo5XZu3cvFhYWjB492qzyNbU1\nrNu1hbFDRhIc0PBeEcUVpXy36xdGhw/m7oiRDZYzSiNvbPmKxPw05o9/imeHPISd1bXtyWGuYR37\nsOzBt2nv5sfrm77gbHFOg2U1lhrenvY81TU6/rt7jcnzjh0yEg9Xdzbt3mZ2XXr16oW1tTWxsbFm\nH6MobYEKjVYmJiaGsLAws2clHz4RTWVVJWMGjzBZ7oe966nR1zJr9AMmy/3nyFr2ph3jucEPMrxT\n8y8d7uPkwQd3PY21xoo3tnyF3qBvsGyAhw8jew1i7ZHtVFZrGyynsbRkcN8BHDgehV7f8Pl+d4xG\nQ/v27Tlz5kyTr0FRWrMWDQ0hxDghxCkhxGkhxMtXeN5GCPFD/fOHhBDBzV/Lm0tycnKTtiaNOhmD\njbU1Ebf1MVluS8weBob1poNPw8Nqy3WVLI1ay6hOEUzpOcrsOlxv3o7uvDD0YRLyznAw/aTJspMj\nx1BZXcXRFNMtgv639aayqpLUzLNm1yMwMJDs7GyzyytKW9BioSGEsAQ+B8YD3YAHhRDdLik2EyiW\nUnYCPgHeb95a3lxqamooLy/H09PT7GMyc7IJ8PE3ubdEaVU5mYW59O5oOox2p0ZTY6jlod7jrnm7\n1Ws1olN/7Kxs2Jd2zGS5sMCOAJw5l2Gy3PnNmwqKC02Wu5idnR063ZW3tVWUtqolWxoRwGkpZaqU\nsgZYAVzauzsJ+Kb+8SpgpGjpd6sWVF5eDoCzs3kzsgHyivLxaWSpjLRzmQB08jO9X0Z0ViJuds50\n8+lo9uvfKNaWVvT060xsborJco629rg6OJNdlGey3Pl1t8oqys2ug4WFxYUtaBXlVtGSoREAXPzx\nL7P+Z1csI6XUA6WA+Wt+tzHnFyNsyqdbK41Vo/fpz6/XZGq/CwCNhSUWQrR4K+O8Gn0tjtamN1iS\nUlJZXYWjrelhsefDwsXR/EAuKCho0hL0itIWtImOcCHE40KIKCFEVH5+fktX54axs7PDwsKCsrIy\ns4+xt7OnoqrCdBmbujfeskbKOVjbUVZdiU5fY/br3ygGo5HM0jw87E3vqFdUUUqtQY+7k+ly5wrq\nWiIebuZv+JSbm4u3t/kLHipKW9CSoZEFXNzrGlj/syuWEUJoABfgspvOUsqvpZT9pJT92vLuahYW\nFgQFBZGSYvqWzMVCgoJJST9DTW3Db/TtPP1xtnPk6GnTncq3B4dTa9SzPnGf2a9/o+xMOUJ+ZTEj\nOpmesb4n/ggA/TrdZrJc1MljONjZ09HMLW3Ly8s5e/YsYWFh5lVYUdqIlgyNI0BnIUQHIYQ1MBW4\ndED9GuCR+sdTgB2yra170kRdu3YlPj4ec38NPcO6U1NbS/zpUw2W0VhaMiAsnH0JUehMhEvfgC50\n9+nIt1HrKKuubHLdrxdtrY5Fh3+lvasfd4Q0POxXSsnGo7to5+lHqH/DqwAbDAb2RR+if88+Zm9t\ne/z4caSU9OzZs8n1V5TWrMVCo76P4klgM5AArJRSxgkh3hJCTKwvthjwEEKcBp4HLhuWe6uJiIig\noKCAuLg488rf1gc7W1vW7thkstzEyNEUV5axcu/6BssIIXh68IMUVJbwwrpPqG6B21TVtTr+tv6f\nnC3O4dkhD5qchb477jAxqfHcf/t4k/0wOw/tIa8wn3FDGp7UeKmtW7fi4OBAr169mlR/RWntWrRP\nQ0q5QUoZKqUMkVLOq//Z36WUa+ofV0sp75dSdpJSRkgpU1uyvjeD4cOHY21tzcaNG80q7+jgyJ1D\nR7N573aKS0saLBfRuSe3d+nLf7avoqCsqMFy4f6hvDnmz8TmnOaVDZ9RUdPwpLnrraCyhOfWfsyR\njHheGzmTQcENv2GXVpXz4S8L6eTXnimDxjdYzmg08u2vKwjyC+SOiMFm1aOyspIdO3YwZswYs3dK\nVJS2ok10hN9KHB0dGTZsGOvWraOiwnTH9XlT77oPvcHAp999ZbLcMxMfRW8w8Nf/zKe6tuERWiM7\nR/DS8Ec4mH6CPyx7mU2n9pt9u+xqGKWR9Ql7efC/rxKXm8KbY/7MhK5DGiyv1VXz/OJ5FFeU8tr9\nc9BYNnzLadWm1SSkJDFjysNmb+H63//+F61Wy+TJk5t8LYrS2qnQaIUefvhhKisr+fXXX80qHxwY\nxMMTH2Ddzs0ciY1uuJx3IG9Ne474zNO8vuxjavS1DZad3GM4i+//O94Obryx5Sue+PldorMSGx22\n2xS1Bj1bkg7y8PdzeWvbQoLd/Vj24DuMDWt4CXNdbQ2vfPcP4tKTeefh5+ke1PDCixk5Wfxr2UIG\n9o7gzmHmreVVVFTEsmXLGDFiBN26XToXVVHaPrWfRis1e/ZskpKSWLVqFW5ubo2Wr9bpmPbiY1RU\nVrJ0/uf4eV95aXGAlXs38OGvC+kZ3IX5f/wrns4ND0M1SiPr4vfw+f4fKakux8vBleGd+jMipD/d\nfTua3PviSsqqKzmSEceB9Fj2nomhWFtOkKsvj0VOZlTnCCxEw59zsovyeH3ZR5xMT+LVKbO5Z8CY\nBsuWlJUy67WnKS4rYfnHi8za8lVKyd/+9jf27NnDihUr1HavSpuiNmFq41JSUpg2bRpjxozhrbfe\nMuuYtMx0/vTKX/Dx9GbRvE9xtG94wtu24/t4a8W/cLJz4L0//pWewV1Mnltbq+O31Gi2nz7MgbOx\n1Bhq0VhYEuIRSBfvYNq5+OBq54SrnRNWlhoMRiM6fQ3aWh2ZpXmcLc7mbHEOqUVZGKXEycaeyHY9\nmNBtCJFBPUyG+fHyuQAAEKNJREFUBcDO2IO8s/IzjFLy9z88yfDbGm6NVGm1/OXNF0lKO83nb3xI\neFfTw3HP27hxI3PnzuWpp57ikUceafwARWlFVGjcAr788ksWL17M22+/zfjxDXf2XuzQ8aM8885L\ndAkJZcFr7+FqYtJbcnYaf106n5ziPO7qO5zZ46fh5dL45LfKGi2H0+OIz0slMS+NxLw0ynQND9G1\nEAJ/Z286uPvR2bM9A9vfRjefjo3uDQ6QlpfJFxuWsevkIboGduLd6S8Q4NFwKyq/qIDn33uNpLQU\n3n/x/7gj0rzO79OnTzNz5kxCQkJYuHCh2f0fitJaqNC4Bej1eubMmUNcXByLFy+mSxfTrYHzdh/e\nx6sfv4Wfly8LXnuXQN9LV2/5nwptJUu2r+KHPeuwtLBk+vDJTB0yASc783erk1KirdVRUl1Oibac\nGoMeSwsLbDXW2FnZ4O3o3uTbWOn52SzeupLNMXuwsbJmxqgpPDR0Ilaahs9zKjWZ5997jfKqCuY9\nN9fs7V0LCwt59NFHqa2tZenSpfj6NhxKitJaqdC4RRQVFTF9+nQAvvrqKwIDzdub+1hCLC/On4tE\n8trsFxgxYKjJ8lmFufxr/bfsOHEAWysbRoffzj0DxtAjKLTZ1qLS6qrZHXeYTdG7OXjqGNYaK+6/\n/U4evmMSbo4Nt5gMBgMrN/7K5/9dhIuTM5+8+i6hwSFmvWZJSQl/+ctfSEtLY+HCharzW2mzVGjc\nQpKSknjiiSewt7fn66+/xt/f36zjMnKyeO2Td0hIOcXEEeN59tHZODmY3twpKesMq/ZvZHPMHrQ1\n1fi5eTO4Wz8Gd+tH35AeWJv4pH81zpUUcCT5BAeTjrE3/ghVump8Xb0Y33cYDwy+Cw8nV5PHp2Wl\n8/bn/+DEqTgG9x3Aa7NfwNPNvEUGCwsLmTNnDhkZGXz44YcMGtTwXuOK0tqp0LjFJCYmMnv2bBwd\nHfn000/p0KHhZTMuVltby9crv+HbX1fg5uzK7AdnMGH42Ebv2VdWa9l2fB974g9zKOk4utoarCw1\nhPp3ICywI538ggn2DsDPzQtvFw+Tt42klJRpKygoKyazIIeU3HRSctNJzEwho6BuS1d3RxeGdO/P\n+D53EN6hKxaN7EdeWl7Gd6t/YMX6n7CxtuGFGU8yfugos1tFZ8+e5bnnniMvL4+PP/6YiIgIs45T\nlNZKhcYtKDExkaeffhqdTsf8+fMZONC8e/YAialJfLDoU2JPxRMS1IHHHvgjd0QMNqvDt7pWR1Ry\nLNEpJ0nITCExK5UqnfbChD8hBE62DtjZ2GJvbYuVxopagx69QY+utoaSyrLL5oT4u/vQyS+IPiE9\niOjcixDfILPe8CsqK/h+/U/8d+0qqrRVjBk8gmcfecLs1gXU7cP+2muvYWVlxYcffkh4eLjZxypK\na6VC4xaVk5PD888/T0pKCk899RTTpk1r9FP5eVJKdh7cw+fLF5GenUmQXyBTJ9zHXcPGYG9net+K\nixmNRnJLCsgqzCW3OJ+c4nzKtBVodVqqdNXU6GuxstRgpbHCylKDq6Mzns5ueDm74+vmRUefIBxs\nzX89gNSMNH7eso51OzdRqa1ieOQQHp/6KJ2CzGtxQd3OiF9//TVLly4lLCyMDz/8ED8/vybVQ1Fa\nKxUat7CqqireeOMNdu7cSWRkJHPnzm3SiB+DwcCuw3v59tcVxJ8+ha2NLcMibmf8kJFE9upn9kqw\nN9q5wnx2HdrD1n27OJ54EiuNFSMGDuXhiffTpWPDM8GvJDExkTfeeIOUlBTuueceXnzxRbWulHJL\nUaFxi5NS8vPPP7NgwQIsLCx4+umnmTx5stmtjvPnOJmcwLqdm9m2fxdlFeW4ODkzoFc/BvaOILJX\nPzybsGnRtaqtrSXudCJHYmPYH32Ik8kJAHQIbM/dw8cyYfg43FxMd4xfqqKigsWLF7N8+XLc3Nx4\n/fXXGTzYvLkbitKWqNBQAMjKyuLtt98mKiqKbt268eyzz9KnT58mn6e2tpb9xw6z48BvHDwWRVFp\nMQBBfoF07RRG905dCOvQiSD/QDxc3a95GK6upobM3CySz6ZyKjWZU2dOE5sUT7WuGiEEXTqGMjxy\nMMMjhxAcGNTk8xsMBlavXs2///1vioqKmDhxIs8++2yT9l9XlLZEhYZygZSSTZs28emnn5Kfn8+g\nQYOYM2eO2ZMBL2U0GklOS+HIyRiOJ8QSf/oUeUUFF563t7UjwNcfTzcP3F1ccXN2xdHBAWsra2ys\nbdBYatDXd4Tr9XrKKysoqyintLyU/KJCsvNyKSj+3waN1lZWhAR1oGdYd/r2CKdPt164OF3dm7te\nr2fr1q0sXryYtLQ0wsPDee655+jevftVnU9R2goVGsplqqurWblyJUuXLqWsrIzbb7+dhx56iIiI\niGtuGeQXFXD67BkycjJJz8ki61w2hSXFFJcWU1RaTE1twyvmWlpY4OzojIuTM+6ubvh7++Lv7UeA\njx+hwSEEBwRdcz9KZWUlGzdu5LvvviMrK4tOnTrx+OOPM3z48GabnKgoNzMVGkqDysvLWbFiBT/+\n+CNFRUWEhIQwdepURo8ejaOj6cl9V0NKicFgQFdbQ01NDbX6WjQaDRrLui87W9sb8sYtpSQhIYFf\nfvmFTZs2odVq6d69O3/6058YOnRok/p3FKWtU6GhNEqn07FlyxaWL19OcnIyNjY2DB8+nAkTJtCv\n380zSqqpMjMz2bx5M5s3byY1NRUbGxvGjBnDvffeS48ePVTLQlGuQIWGYjYpJSdPnmTdunVs2bKF\n8vJynJycGDRoEEOGDGHAgAG4ujZtVFJzqq2tJTY2lgMHDnDw4EESEupGVYWHhzNu3DjGjRt3Q1pQ\nitKWqNBQropOp2P//v389ttv7N27l+LiulFSnTt3pl+/fvTr14/u3bvj6enZYnWsrq4mMTGRo0eP\nEh0dzYkTJ9BqtVhaWtKjRw+GDBnC2LFj1cQ8RWkCFRrKNTMYDMTHx3P48GGOHj3K8ePH0enq9g73\n9vamS5cuhISE0L59e4KCgggICMDNze269BVIKamoqCArK4vMzEwyMzM5ffo0SUlJpKWlYTTWbSsb\nEhJCnz59iIiIoH///qpFoShXSYWGct3pdDoSEhKIj48nISGBhIQEMjIyMBgMF8pYWlri6emJh4cH\nLi4uODs74+joiKOjI9bW1mg0GiwtLTEYDOj1dUNuq6qqKC8vp7y8nNLSUgoLCyksLKS6uvp3r+/j\n40NoaCihoaF07dqV8PDwm/q2maK0Jio0lGah1+vJzs4mPT2drKwsCgoKKCgoID8/n7KyMsrLy6mo\nqKCiooLaKwy7tbCwwMHBAScnJxwdHXFxccHDwwMPDw88PT3x9/cnMDCQgIAA1YpQlBvI3NBoncNj\nlJuGRqMhKCiIoKDGZ2VLKdHr9RgMBjQaDRYWFmrYq6K0Mio0lGYjhMDKygorq+u7UZOiKM1HfcxT\nFEVRzKZCQ1EURTGbCg1FURTFbC0SGkIIdyHEViFEcv2fbg2UMwghjtV/rWnueiqKoii/11ItjZeB\n7VLKzsD2+u+vRCulDK//mth81VMURVGupKVCYxLwTf3jb4B7WqgeiqIoShO0VGj4SClz6h/nAj4N\nlLMVQkQJIQ4KIVSwKIqitLAbNk9DCLEN8L3CU69d/I2UUgohGpqW3l5KmSWE6AjsEELESilTrvBa\njwOPA2ZNMlMURVGuzg0LDSnlqIaeE0KcE0L4SSlzhBB+QF4D58iq/zNVCLEL6A1cFhpSyq+Br6Fu\nGZHrUH1FURTlClrq9tQa4JH6x48Aqy8tIIRwE0LY1D/2BG4H4puthoqiKMplWio05gOjhRDJwKj6\n7xFC9BNCLKov0xWIEkIcB3YC86WUKjQURVFaUIusPSWlLARGXuHnUcCs+sf7gduauWqKoiiKCWpG\nuKIoimI2FRqKoiiK2VRoKIqiKGZToaEoiqKYTYWGoiiKYrY2t0e4ECIfONvS9TCTJ1DQ0pVoZuqa\nbx234nW35mtuL6X0aqxQmwuN1kQIEWXORu5tibrmW8eteN23wjWr21OKoiiK2VRoKIqiKGZTodGy\nvm7pCrQAdc23jlvxutv8Nas+DUVRFMVsqqWhKIqimE2FRjMSQrgLIbYKIZLr/3RroJxBCHGs/mtN\nc9fzehBCjBNCnBJCnBZCXLYHvBDCRgjxQ/3zh4QQwc1fy+vLjGt+VAiRf9Hf7ayWqOf1JIRYIoTI\nE0KcbOB5IYT4tP53ckII0ae563i9mXHNdwghSi/6e/57c9fxRlKh0bxeBrZLKTsD2+u/vxKtlDK8\n/mti81Xv+hBCWAKfA+OBbsCDQohulxSbCRRLKTsBnwDvN28try8zrxngh4v+bhdd4fnWZikwzsTz\n44HO9V+PA182Q51utKWYvmaAPRf9Pb/VDHVqNio0mtck4Jv6x98AbXXf8wjgtJQyVUpZA6yg7tov\ndvHvYhUwUgghmrGO15s519zmSCl/A4pMFJkEfCvrHARc63frbLXMuOY2TYVG8/KRUubUP84FfBoo\nZyuEiBJCHBRCtMZgCQAyLvo+s/5nVywjpdQDpYBHs9TuxjDnmgHuq79Ns0oI0a55qtaizP29tDUD\nhRDHhRAbhRDdW7oy11OLbMLUlgkhtgG+V3jqtYu/kVJKIURDQ9faSymzhBAdgR1CiFgp5WV7oyut\nzlrgeymlTgjxZ+paWiNauE7K9RdN3f/hCiHEncCv1N2eaxNUaFxnUspRDT0nhDgnhPCTUubUN9Hz\nGjhHVv2fqUKIXUBvoDWFRhZw8afowPqfXalMphBCA7gAhc1TvRui0Wuu37HyvEXAB81Qr5Zmzr+F\nNkVKWXbR4w1CiC+EEJ5Syta6JtXvqNtTzWsN8Ej940eA1ZcWEEK4CSFs6h97ArcDrW1v9CNAZyFE\nByGENTCVumu/2MW/iynADtm6Jw01es2X3MufCCQ0Y/1ayhrgj/WjqAYApRfdom2ThBC+5/vnhBAR\n1L3PtuYPRL+jWhrNaz6wUggxk7qVeB8AEEL0A56QUs4CugJfCSGM1P1jmy+lbFWhIaXUCyGeBDYD\nlsASKWWcEOItIEpKuQZYDHwnhDhNXafi1Jar8bUz85qfFkJMBPTUXfOjLVbh60QI8T1wB+AphMgE\n3gCsAKSU/wY2AHcCp4Eq4E8tU9Prx4xrngLMFkLoAS0wtZV/IPodNSNcURRFMZu6PaUoiqKYTYWG\noiiKYjYVGoqiKIrZVGgoiqIoZlOhoSiKophNhYaiNBMhRLgQ4oAQIq5+KZE/tHSdFKWp1JBbRWkm\nQohQ6laQSRZC+ANHga5SypIWrpqimE21NBTlBhBC9K9vTdgKIRyEEHGAtZQyGUBKmU3dMjJeLVpR\nRWki1dJQlBtECPEOYAvYAZlSyvcuei6CugULu0spjS1URUVpMhUainKD1K9BdQSoBgZJKQ31P/cD\ndgGP1O8xoSithro9pSg3jgfgCDhR1+JACOEMrAdeU4GhtEaqpaEoN0j9/u4rgA6AH/A8sBFYK6Vc\n0JJ1U5SrpVa5VZQbQAjxR6BWSrm8fv/w/dSt5DsU8BBCPFpf9FEp5bEWqqaiNJlqaSiKoihmU30a\niqIoitlUaCiKoihmU6GhKIqimE2FhqIoimI2FRqKoiiK2VRoKIqiKGZToaEoiqKYTYWGoiiKYrb/\nB4Zpnck22NlhAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f10a2183d30>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# and visualise the samples\n",
    "import seaborn as sns\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "sns.kdeplot(\n",
    "    data=ds.x2,\n",
    "    data2=ds.x3,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And to access the implied CGM\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "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 2.38.0 (20140413.2041)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"62pt\" height=\"188pt\"\n",
       " viewBox=\"0.00 0.00 62.00 188.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 184)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-184 58,-184 58,4 -4,4\"/>\n",
       "<!-- x1 -->\n",
       "<g id=\"node1\" class=\"node\"><title>x1</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-158.3\" font-family=\"Times,serif\" font-size=\"14.00\">x1</text>\n",
       "</g>\n",
       "<!-- x2 -->\n",
       "<g id=\"node3\" class=\"node\"><title>x2</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">x2</text>\n",
       "</g>\n",
       "<!-- x1&#45;&gt;x2 -->\n",
       "<g id=\"edge1\" class=\"edge\"><title>x1&#45;&gt;x2</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M27,-143.697C27,-135.983 27,-126.712 27,-118.112\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"30.5001,-118.104 27,-108.104 23.5001,-118.104 30.5001,-118.104\"/>\n",
       "</g>\n",
       "<!-- x3 -->\n",
       "<g id=\"node2\" class=\"node\"><title>x3</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x3</text>\n",
       "</g>\n",
       "<!-- x2&#45;&gt;x3 -->\n",
       "<g id=\"edge2\" class=\"edge\"><title>x2&#45;&gt;x3</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M27,-71.6966C27,-63.9827 27,-54.7125 27,-46.1124\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"30.5001,-46.1043 27,-36.1043 23.5001,-46.1044 30.5001,-46.1043\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x7f10a219f320>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scm.cgm.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And to apply an intervention:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "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 2.38.0 (20140413.2041)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"70pt\" height=\"196pt\"\n",
       " viewBox=\"0.00 0.00 70.00 196.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 192)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-192 66,-192 66,4 -4,4\"/>\n",
       "<!-- x1 -->\n",
       "<g id=\"node1\" class=\"node\"><title>x1</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"31\" cy=\"-166\" rx=\"27\" ry=\"18\"/>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"31\" cy=\"-166\" rx=\"31\" ry=\"22\"/>\n",
       "<text text-anchor=\"middle\" x=\"31\" y=\"-162.3\" font-family=\"Times,serif\" font-size=\"14.00\">x1</text>\n",
       "</g>\n",
       "<!-- x2 -->\n",
       "<g id=\"node3\" class=\"node\"><title>x2</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"31\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"31\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">x2</text>\n",
       "</g>\n",
       "<!-- x1&#45;&gt;x2 -->\n",
       "<g id=\"edge1\" class=\"edge\"><title>x1&#45;&gt;x2</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M31,-143.993C31,-136.061 31,-126.915 31,-118.478\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"34.5001,-118.315 31,-108.315 27.5001,-118.315 34.5001,-118.315\"/>\n",
       "</g>\n",
       "<!-- x3 -->\n",
       "<g id=\"node2\" class=\"node\"><title>x3</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"31\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"31\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x3</text>\n",
       "</g>\n",
       "<!-- x2&#45;&gt;x3 -->\n",
       "<g id=\"edge2\" class=\"edge\"><title>x2&#45;&gt;x3</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M31,-71.6966C31,-63.9827 31,-54.7125 31,-46.1124\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"34.5001,-46.1043 31,-36.1043 27.5001,-46.1044 34.5001,-46.1043\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x7f109cf82518>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scm_do = scm.do(\"x1\")\n",
    "\n",
    "scm_do.cgm.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And sample from the distribution implied by this intervention:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>x1</th>\n",
       "      <th>x2</th>\n",
       "      <th>x3</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>0.258736</td>\n",
       "      <td>0.066944</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>1.142529</td>\n",
       "      <td>1.305371</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>1.927665</td>\n",
       "      <td>3.715894</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>2.878977</td>\n",
       "      <td>8.288507</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "      <td>4.098250</td>\n",
       "      <td>16.795653</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   x1        x2         x3\n",
       "0   0  0.258736   0.066944\n",
       "1   1  1.142529   1.305371\n",
       "2   2  1.927665   3.715894\n",
       "3   3  2.878977   8.288507\n",
       "4   4  4.098250  16.795653"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scm_do.sample(n_samples=5, set_values={\"x1\": np.arange(5)})"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}