notebooks/linear_income.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sanity check: income dataset with linear models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# limit the number of threads used by xgboost models\n",
    "import os\n",
    "os.environ[\"OMP_NUM_THREADS\"] = \"8\"\n",
    "os.environ[\"OPENBLAS_NUM_THREADS\"] = \"8\"\n",
    "os.environ[\"MKL_NUM_THREADS\"] = \"8\"\n",
    "os.environ[\"VECLIB_MAXIMUM_THREADS\"] = \"8\"\n",
    "os.environ[\"NUMEXPR_NUM_THREADS\"] = \"8\"\n",
    "\n",
    "# standard libaries\n",
    "import sys\n",
    "import copy\n",
    "import dill\n",
    "from collections import defaultdict\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xgboost\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import shap\n",
    "\n",
    "# flow specific packages\n",
    "if '../' not in sys.path:\n",
    "    sys.path = ['../'] + sys.path\n",
    "from flow import Node, CreditFlow, Graph, get_source_nodes, topo_sort, flatten_graph, eval_graph, boundary_graph, single_source_graph, viz_graph, save_graph, hcluster_graph\n",
    "from flow import ParallelCreditFlow, GraphExplainer, translator\n",
    "from flow import group_nodes, build_feature_graph\n",
    "from flow import CausalLinks, create_xgboost_f, create_linear_f\n",
    "from flow import edge_credits2edge_credit\n",
    "from on_manifold import OnManifoldExplainer, IndExplainer, FeatureAttribution\n",
    "from linear_evaluation import get_error, get_indirect_effect_flow, get_effect_asv, get_effect_ind, get_effect_manifold, get_direct_effect_flow\n",
    "from linear_evaluation import get_direct_effect_ground_truth, get_indirect_effect_ground_truth\n",
    "\n",
    "np.random.seed(42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dataset and Training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "X, y = shap.datasets.adult()\n",
    "X_display,y_display = shap.datasets.adult(display=True) # human readable feature values\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)\n",
    "model = LogisticRegression(max_iter=10000).fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Explanation setting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "nsamples = 100 # paper: 1000 samples\n",
    "bg = X.fillna(X.mean()).sample(1) # background sample, no need for multiple sample because of the linearity axiom\n",
    "fg = X[:nsamples] # foreground samples to explain"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Build a causal graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "learning dependency for Marital Status: 100%|██████████| 8/8 [00:00<00:00, 85.40it/s]\n"
     ]
    },
    {
     "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=\"1121pt\" height=\"188pt\"\n",
       " viewBox=\"0.00 0.00 1121.02 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 1117.02,-184 1117.02,4 -4,4\"/>\n",
       "<!-- Marital Status -->\n",
       "<g id=\"node1\" class=\"node\"><title>Marital Status</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"864.999\" cy=\"-90\" rx=\"61.1893\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"864.999\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">Marital Status</text>\n",
       "</g>\n",
       "<!-- predicted income &gt;50k -->\n",
       "<g id=\"node2\" class=\"node\"><title>predicted income &gt;50k</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"539.999\" cy=\"-18\" rx=\"92.8835\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"539.999\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">predicted income &gt;50k</text>\n",
       "</g>\n",
       "<!-- Marital Status&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge1\" class=\"edge\"><title>Marital Status&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M819.864,-77.8049C811.604,-75.8226 803.051,-73.814 794.999,-72 733.587,-58.165 663.555,-43.7032 612.44,-33.3879\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"612.972,-29.9248 602.478,-31.3815 611.59,-36.787 612.972,-29.9248\"/>\n",
       "</g>\n",
       "<!-- Age -->\n",
       "<g id=\"node3\" class=\"node\"><title>Age</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"416.999\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"416.999\" y=\"-158.3\" font-family=\"Times,serif\" font-size=\"14.00\">Age</text>\n",
       "</g>\n",
       "<!-- Age&#45;&gt;Marital Status -->\n",
       "<g id=\"edge2\" class=\"edge\"><title>Age&#45;&gt;Marital Status</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M437.479,-149.881C442.437,-147.6 447.807,-145.454 452.999,-144 600.177,-102.787 644.403,-134.093 794.999,-108 800.313,-107.079 805.823,-105.986 811.312,-104.805\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"812.133,-108.208 821.124,-102.603 810.6,-101.378 812.133,-108.208\"/>\n",
       "</g>\n",
       "<!-- Age&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge3\" class=\"edge\"><title>Age&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M390.237,-159.547C303.359,-154.56 33.344,-136.842 7.99855,-108 -2.56321,-95.9813 -2.73974,-83.8612 7.99855,-72 36.2675,-40.7748 291.917,-27.1841 438.283,-21.919\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"438.758,-25.4046 448.628,-21.5542 438.511,-18.4089 438.758,-25.4046\"/>\n",
       "</g>\n",
       "<!-- Hours per week -->\n",
       "<g id=\"node4\" class=\"node\"><title>Hours per week</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"1011\" cy=\"-90\" rx=\"66.8882\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"1011\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">Hours per week</text>\n",
       "</g>\n",
       "<!-- Age&#45;&gt;Hours per week -->\n",
       "<g id=\"edge4\" class=\"edge\"><title>Age&#45;&gt;Hours per week</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M437.068,-149.942C442.131,-147.599 447.652,-145.41 452.999,-144 660.715,-89.2176 722.506,-139.53 934.999,-108 940.91,-107.123 947.046,-106.042 953.153,-104.853\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"954.261,-108.199 963.358,-102.769 952.86,-101.341 954.261,-108.199\"/>\n",
       "</g>\n",
       "<!-- Capital Loss -->\n",
       "<g id=\"node5\" class=\"node\"><title>Capital Loss</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"83.9986\" cy=\"-90\" rx=\"55.4913\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"83.9986\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">Capital Loss</text>\n",
       "</g>\n",
       "<!-- Age&#45;&gt;Capital Loss -->\n",
       "<g id=\"edge5\" class=\"edge\"><title>Age&#45;&gt;Capital Loss</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M391.016,-156.403C343.02,-147.79 237.294,-128.257 148.999,-108 144.408,-106.947 139.649,-105.803 134.889,-104.623\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"135.721,-101.224 125.169,-102.169 134.007,-108.011 135.721,-101.224\"/>\n",
       "</g>\n",
       "<!-- Occupation -->\n",
       "<g id=\"node6\" class=\"node\"><title>Occupation</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"208.999\" cy=\"-90\" rx=\"51.1914\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"208.999\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">Occupation</text>\n",
       "</g>\n",
       "<!-- Age&#45;&gt;Occupation -->\n",
       "<g id=\"edge6\" class=\"edge\"><title>Age&#45;&gt;Occupation</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M393.138,-152.97C359.38,-141.609 297.038,-120.629 254.323,-106.254\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"255.316,-102.895 244.722,-103.022 253.083,-109.529 255.316,-102.895\"/>\n",
       "</g>\n",
       "<!-- Education&#45;Num -->\n",
       "<g id=\"node7\" class=\"node\"><title>Education&#45;Num</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"345.999\" cy=\"-90\" rx=\"67.6881\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"345.999\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">Education&#45;Num</text>\n",
       "</g>\n",
       "<!-- Age&#45;&gt;Education&#45;Num -->\n",
       "<g id=\"edge7\" class=\"edge\"><title>Age&#45;&gt;Education&#45;Num</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M402.287,-146.496C392.915,-137.255 380.557,-125.072 369.809,-114.475\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"372.26,-111.976 362.681,-107.448 367.345,-116.961 372.26,-111.976\"/>\n",
       "</g>\n",
       "<!-- Workclass -->\n",
       "<g id=\"node8\" class=\"node\"><title>Workclass</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"479.999\" cy=\"-90\" rx=\"48.1917\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"479.999\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">Workclass</text>\n",
       "</g>\n",
       "<!-- Age&#45;&gt;Workclass -->\n",
       "<g id=\"edge8\" class=\"edge\"><title>Age&#45;&gt;Workclass</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M430.359,-146.155C438.533,-137.072 449.18,-125.243 458.516,-114.87\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"461.249,-117.065 465.337,-107.29 456.046,-112.382 461.249,-117.065\"/>\n",
       "</g>\n",
       "<!-- Relationship -->\n",
       "<g id=\"node9\" class=\"node\"><title>Relationship</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"600.999\" cy=\"-90\" rx=\"55.4913\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"600.999\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">Relationship</text>\n",
       "</g>\n",
       "<!-- Age&#45;&gt;Relationship -->\n",
       "<g id=\"edge9\" class=\"edge\"><title>Age&#45;&gt;Relationship</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M438.165,-150.586C443.001,-148.324 448.145,-146.009 452.999,-144 486.302,-130.212 524.634,-116.544 553.942,-106.542\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"555.168,-109.822 563.514,-103.296 552.92,-103.193 555.168,-109.822\"/>\n",
       "</g>\n",
       "<!-- Capital Gain -->\n",
       "<g id=\"node10\" class=\"node\"><title>Capital Gain</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"729.999\" cy=\"-90\" rx=\"55.7903\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"729.999\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">Capital Gain</text>\n",
       "</g>\n",
       "<!-- Age&#45;&gt;Capital Gain -->\n",
       "<g id=\"edge10\" class=\"edge\"><title>Age&#45;&gt;Capital Gain</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M437.55,-150.119C442.505,-147.827 447.856,-145.617 452.999,-144 544.168,-115.329 571.567,-128.11 664.999,-108 669.603,-107.009 674.371,-105.908 679.138,-104.756\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"680.005,-108.147 688.866,-102.338 678.317,-101.353 680.005,-108.147\"/>\n",
       "</g>\n",
       "<!-- Hours per week&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge11\" class=\"edge\"><title>Hours per week&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M962.473,-77.4879C953.372,-75.5159 943.916,-73.5926 934.999,-72 830.487,-53.3359 709.357,-38.0054 629.075,-28.7172\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"629.337,-25.2244 619.003,-27.5589 628.537,-32.1786 629.337,-25.2244\"/>\n",
       "</g>\n",
       "<!-- Capital Loss&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge12\" class=\"edge\"><title>Capital Loss&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M124.392,-77.627C132.499,-75.5732 140.989,-73.5843 148.999,-72 251.143,-51.7969 369.833,-36.9005 449.411,-28.1314\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"449.842,-31.6053 459.403,-27.04 449.082,-24.6467 449.842,-31.6053\"/>\n",
       "</g>\n",
       "<!-- Occupation&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge13\" class=\"edge\"><title>Occupation&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M246.962,-77.7302C254.248,-75.7132 261.838,-73.7115 268.999,-72 333.962,-56.4725 408.486,-42.1476 462.92,-32.3157\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"463.62,-35.7461 472.843,-30.5329 462.382,-28.8564 463.62,-35.7461\"/>\n",
       "</g>\n",
       "<!-- Education&#45;Num&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge14\" class=\"edge\"><title>Education&#45;Num&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M384.327,-75.1703C414.176,-64.3998 455.9,-49.345 488.638,-37.532\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"489.832,-40.8221 498.051,-34.1357 487.457,-34.2376 489.832,-40.8221\"/>\n",
       "</g>\n",
       "<!-- Workclass&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge15\" class=\"edge\"><title>Workclass&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M493.915,-72.7646C501.271,-64.1822 510.459,-53.4625 518.692,-43.858\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"521.47,-45.9943 525.321,-36.1239 516.156,-41.4387 521.47,-45.9943\"/>\n",
       "</g>\n",
       "<!-- Relationship&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge16\" class=\"edge\"><title>Relationship&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M586.542,-72.411C579.056,-63.8197 569.766,-53.1594 561.46,-43.6275\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"563.983,-41.1959 554.775,-35.9562 558.706,-45.7948 563.983,-41.1959\"/>\n",
       "</g>\n",
       "<!-- Capital Gain&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge17\" class=\"edge\"><title>Capital Gain&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M694.706,-75.9976C665.501,-65.2378 623.636,-49.8137 590.819,-37.7232\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"591.979,-34.4209 581.386,-34.248 589.559,-40.9893 591.979,-34.4209\"/>\n",
       "</g>\n",
       "<!-- Sex -->\n",
       "<g id=\"node11\" class=\"node\"><title>Sex</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"488.999\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"488.999\" y=\"-158.3\" font-family=\"Times,serif\" font-size=\"14.00\">Sex</text>\n",
       "</g>\n",
       "<!-- Sex&#45;&gt;Marital Status -->\n",
       "<g id=\"edge18\" class=\"edge\"><title>Sex&#45;&gt;Marital Status</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M512.427,-152.708C521.835,-149.62 532.839,-146.309 542.999,-144 653.323,-118.93 683.899,-129.375 794.999,-108 800.087,-107.021 805.363,-105.913 810.631,-104.745\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"811.736,-108.082 820.703,-102.439 810.174,-101.259 811.736,-108.082\"/>\n",
       "</g>\n",
       "<!-- Sex&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge19\" class=\"edge\"><title>Sex&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M468.54,-149.799C463.582,-147.521 458.206,-145.397 452.999,-144 406.37,-131.493 52.6908,-143.523 19.9986,-108 9.16368,-96.2269 9.27376,-83.8735 19.9986,-72 47.4901,-41.564 295.483,-27.6838 438.833,-22.1549\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"439.109,-25.6471 448.969,-21.7712 438.844,-18.6521 439.109,-25.6471\"/>\n",
       "</g>\n",
       "<!-- Sex&#45;&gt;Hours per week -->\n",
       "<g id=\"edge20\" class=\"edge\"><title>Sex&#45;&gt;Hours per week</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M512.105,-152.532C521.563,-149.381 532.688,-146.072 542.999,-144 714.526,-109.536 762.163,-135.151 934.999,-108 940.903,-107.073 947.034,-105.958 953.137,-104.75\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"954.253,-108.093 963.339,-102.643 952.837,-101.238 954.253,-108.093\"/>\n",
       "</g>\n",
       "<!-- Sex&#45;&gt;Capital Loss -->\n",
       "<g id=\"edge21\" class=\"edge\"><title>Sex&#45;&gt;Capital Loss</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M468.504,-149.93C463.547,-147.647 458.181,-145.487 452.999,-144 322.224,-106.462 282.686,-133.273 148.999,-108 144.235,-107.099 139.305,-106.043 134.388,-104.906\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"134.914,-101.432 124.371,-102.48 133.266,-108.235 134.914,-101.432\"/>\n",
       "</g>\n",
       "<!-- Sex&#45;&gt;Occupation -->\n",
       "<g id=\"edge22\" class=\"edge\"><title>Sex&#45;&gt;Occupation</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M468.417,-150.212C463.464,-147.917 458.121,-145.682 452.999,-144 373.827,-118.01 350.043,-127.374 268.999,-108 264.971,-107.037 260.807,-105.982 256.638,-104.885\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"257.529,-101.5 246.962,-102.269 255.703,-108.258 257.529,-101.5\"/>\n",
       "</g>\n",
       "<!-- Sex&#45;&gt;Education&#45;Num -->\n",
       "<g id=\"edge23\" class=\"edge\"><title>Sex&#45;&gt;Education&#45;Num</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M467.898,-150.671C446.424,-140.159 412.608,-123.606 386.035,-110.598\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"387.49,-107.414 376.97,-106.161 384.412,-113.701 387.49,-107.414\"/>\n",
       "</g>\n",
       "<!-- Sex&#45;&gt;Workclass -->\n",
       "<g id=\"edge24\" class=\"edge\"><title>Sex&#45;&gt;Workclass</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M486.82,-144.055C485.829,-136.346 484.631,-127.027 483.517,-118.364\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"486.966,-117.747 482.22,-108.275 480.024,-118.64 486.966,-117.747\"/>\n",
       "</g>\n",
       "<!-- Sex&#45;&gt;Relationship -->\n",
       "<g id=\"edge25\" class=\"edge\"><title>Sex&#45;&gt;Relationship</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M507.998,-149.125C524.343,-138.91 548.299,-123.937 567.804,-111.747\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"569.784,-114.637 576.409,-106.369 566.074,-108.701 569.784,-114.637\"/>\n",
       "</g>\n",
       "<!-- Sex&#45;&gt;Capital Gain -->\n",
       "<g id=\"edge26\" class=\"edge\"><title>Sex&#45;&gt;Capital Gain</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M512.878,-153.32C522.25,-150.348 533.111,-146.957 542.999,-144 588.85,-130.288 641.292,-115.502 678.95,-105.041\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"679.975,-108.389 688.676,-102.344 678.104,-101.643 679.975,-108.389\"/>\n",
       "</g>\n",
       "<!-- Country -->\n",
       "<g id=\"node12\" class=\"node\"><title>Country</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"591.999\" cy=\"-162\" rx=\"39.7935\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"591.999\" y=\"-158.3\" font-family=\"Times,serif\" font-size=\"14.00\">Country</text>\n",
       "</g>\n",
       "<!-- Country&#45;&gt;Marital Status -->\n",
       "<g id=\"edge27\" class=\"edge\"><title>Country&#45;&gt;Marital Status</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M621.531,-149.956C627.599,-147.852 633.974,-145.761 639.999,-144 707.882,-124.164 726.191,-124.345 794.999,-108 799.694,-106.884 804.565,-105.719 809.452,-104.544\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"810.552,-107.879 819.453,-102.131 808.911,-101.074 810.552,-107.879\"/>\n",
       "</g>\n",
       "<!-- Country&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge28\" class=\"edge\"><title>Country&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M620.584,-149.512C626.891,-147.363 633.601,-145.366 639.999,-144 688.727,-133.595 1053.21,-144.62 1087,-108 1097.85,-96.2409 1097.75,-83.8467 1087,-72 1057.71,-39.7376 791.688,-26.5313 641.704,-21.613\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"641.661,-18.1099 631.554,-21.287 641.436,-25.1063 641.661,-18.1099\"/>\n",
       "</g>\n",
       "<!-- Country&#45;&gt;Hours per week -->\n",
       "<g id=\"edge29\" class=\"edge\"><title>Country&#45;&gt;Hours per week</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M620.884,-149.623C627.114,-147.494 633.716,-145.477 639.999,-144 768.577,-113.771 804.836,-130.444 934.999,-108 940.717,-107.014 946.655,-105.881 952.577,-104.677\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"953.419,-108.077 962.488,-102.599 951.983,-101.226 953.419,-108.077\"/>\n",
       "</g>\n",
       "<!-- Country&#45;&gt;Capital Loss -->\n",
       "<g id=\"edge30\" class=\"edge\"><title>Country&#45;&gt;Capital Loss</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M559.499,-151.689C548.609,-148.855 536.349,-145.971 524.999,-144 359.599,-115.272 314.291,-137.333 148.999,-108 144.03,-107.118 138.884,-106.049 133.762,-104.881\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"134.538,-101.468 123.997,-102.535 132.903,-108.274 134.538,-101.468\"/>\n",
       "</g>\n",
       "<!-- Country&#45;&gt;Occupation -->\n",
       "<g id=\"edge31\" class=\"edge\"><title>Country&#45;&gt;Occupation</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M559.104,-151.895C548.296,-149.107 536.185,-146.199 524.999,-144 412.259,-121.838 381.422,-131.716 268.999,-108 264.827,-107.12 260.52,-106.109 256.217,-105.027\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"256.811,-101.564 246.25,-102.405 255.03,-108.334 256.811,-101.564\"/>\n",
       "</g>\n",
       "<!-- Country&#45;&gt;Education&#45;Num -->\n",
       "<g id=\"edge32\" class=\"edge\"><title>Country&#45;&gt;Education&#45;Num</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M559.332,-151.705C518.737,-140.153 448.885,-120.277 400.213,-106.427\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"401.048,-103.026 390.472,-103.655 399.132,-109.758 401.048,-103.026\"/>\n",
       "</g>\n",
       "<!-- Country&#45;&gt;Workclass -->\n",
       "<g id=\"edge33\" class=\"edge\"><title>Country&#45;&gt;Workclass</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M569.603,-147.003C553.255,-136.785 530.843,-122.778 512.524,-111.328\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"514.239,-108.273 503.904,-105.941 510.529,-114.209 514.239,-108.273\"/>\n",
       "</g>\n",
       "<!-- Country&#45;&gt;Relationship -->\n",
       "<g id=\"edge34\" class=\"edge\"><title>Country&#45;&gt;Relationship</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M594.223,-143.697C595.215,-135.983 596.407,-126.712 597.513,-118.112\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"600.996,-118.469 598.799,-108.104 594.053,-117.576 600.996,-118.469\"/>\n",
       "</g>\n",
       "<!-- Country&#45;&gt;Capital Gain -->\n",
       "<g id=\"edge35\" class=\"edge\"><title>Country&#45;&gt;Capital Gain</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M617.632,-147.998C638.551,-137.387 668.413,-122.239 692.095,-110.227\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"693.829,-113.272 701.164,-105.626 690.662,-107.029 693.829,-113.272\"/>\n",
       "</g>\n",
       "<!-- Race -->\n",
       "<g id=\"node13\" class=\"node\"><title>Race</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"677.999\" cy=\"-162\" rx=\"28.6953\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"677.999\" y=\"-158.3\" font-family=\"Times,serif\" font-size=\"14.00\">Race</text>\n",
       "</g>\n",
       "<!-- Race&#45;&gt;Marital Status -->\n",
       "<g id=\"edge36\" class=\"edge\"><title>Race&#45;&gt;Marital Status</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M702.063,-151.992C731.583,-140.942 782.258,-121.973 819.273,-108.116\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"820.706,-111.317 828.844,-104.534 818.252,-104.762 820.706,-111.317\"/>\n",
       "</g>\n",
       "<!-- Race&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge37\" class=\"edge\"><title>Race&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M706.571,-159.607C797.932,-154.863 1078.54,-137.949 1105,-108 1115.59,-96.0097 1115.77,-83.8302 1105,-72 1074.52,-38.5317 796.629,-25.8022 642.256,-21.2835\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"641.911,-17.7723 631.815,-20.9849 641.711,-24.7694 641.911,-17.7723\"/>\n",
       "</g>\n",
       "<!-- Race&#45;&gt;Hours per week -->\n",
       "<g id=\"edge38\" class=\"edge\"><title>Race&#45;&gt;Hours per week</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M705.04,-155.578C751.913,-146.145 851.264,-125.998 934.999,-108 940.369,-106.846 945.948,-105.63 951.536,-104.401\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"952.549,-107.761 961.558,-102.185 951.038,-100.926 952.549,-107.761\"/>\n",
       "</g>\n",
       "<!-- Race&#45;&gt;Capital Loss -->\n",
       "<g id=\"edge39\" class=\"edge\"><title>Race&#45;&gt;Capital Loss</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M656.955,-149.766C651.859,-147.489 646.339,-145.374 640.999,-144 428.665,-89.3599 365.293,-143.888 148.999,-108 143.95,-107.162 138.722,-106.112 133.526,-104.948\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"134.164,-101.502 123.626,-102.592 132.544,-108.312 134.164,-101.502\"/>\n",
       "</g>\n",
       "<!-- Race&#45;&gt;Occupation -->\n",
       "<g id=\"edge40\" class=\"edge\"><title>Race&#45;&gt;Occupation</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M656.538,-149.681C651.554,-147.467 646.186,-145.4 640.999,-144 480.63,-100.723 432.23,-138.764 268.999,-108 264.631,-107.177 260.122,-106.181 255.629,-105.088\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"256.436,-101.682 245.878,-102.572 254.686,-108.46 256.436,-101.682\"/>\n",
       "</g>\n",
       "<!-- Race&#45;&gt;Education&#45;Num -->\n",
       "<g id=\"edge41\" class=\"edge\"><title>Race&#45;&gt;Education&#45;Num</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M656.464,-149.939C651.484,-147.711 646.136,-145.573 640.999,-144 547.102,-115.245 519.424,-126.592 422.999,-108 417.315,-106.904 411.408,-105.707 405.507,-104.471\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"406.132,-101.025 395.622,-102.366 404.674,-107.872 406.132,-101.025\"/>\n",
       "</g>\n",
       "<!-- Race&#45;&gt;Workclass -->\n",
       "<g id=\"edge42\" class=\"edge\"><title>Race&#45;&gt;Workclass</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M655.869,-150.341C651.013,-148.14 645.867,-145.913 640.999,-144 602.898,-129.031 558.665,-114.68 526.214,-104.696\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"527.128,-101.316 516.542,-101.742 525.083,-108.01 527.128,-101.316\"/>\n",
       "</g>\n",
       "<!-- Race&#45;&gt;Relationship -->\n",
       "<g id=\"edge43\" class=\"edge\"><title>Race&#45;&gt;Relationship</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M662.416,-146.834C652.035,-137.396 638.157,-124.78 626.22,-113.929\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"628.466,-111.24 618.712,-107.103 623.757,-116.42 628.466,-111.24\"/>\n",
       "</g>\n",
       "<!-- Race&#45;&gt;Capital Gain -->\n",
       "<g id=\"edge44\" class=\"edge\"><title>Race&#45;&gt;Capital Gain</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M689.539,-145.465C695.926,-136.867 704.021,-125.97 711.288,-116.187\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"714.245,-118.076 717.398,-107.962 708.625,-113.902 714.245,-118.076\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.files.Source at 0x7f72b0533310>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "causal_links = CausalLinks()\n",
    "categorical_feature_names = [] # a linear model doesn't handle categorical features\n",
    "display_translator = translator(X.columns, X, X_display)\n",
    "target_name = 'predicted income >50k' # 'predicted hazard'\n",
    "feature_names = list(X.columns)\n",
    "f = lambda x: x.dot(model.coef_.reshape(-1, 1)).ravel() + model.intercept_ # logit is linear log(p / 1-p)\n",
    "\n",
    "# build causal links\n",
    "causal_links.add_causes_effects(feature_names, target_name, create_linear_f(feature_names, f))\n",
    "A = ['Age', 'Sex', 'Country', 'Race']\n",
    "D = list(set(feature_names) - set(A))\n",
    "causal_links.add_causes_effects(A, D)\n",
    "\n",
    "causal_graph = build_feature_graph(X.fillna(X.mean()), causal_links, categorical_feature_names, display_translator, target_name, method='linear')\n",
    "\n",
    "causal_graph.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Explain\n",
    "\n",
    "Flow, ASV, and SHAP converge in 1 run b/c the model is linear."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "bruteforce sampling: 100%|██████████| 1/1 [00:00<00:00, 58.57it/s]\n"
     ]
    }
   ],
   "source": [
    "cf_flow_adult = GraphExplainer(causal_graph, bg, nruns=1).shap_values(fg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "manifold bg samples: 100%|██████████| 100/100 [01:20<00:00,  1.24it/s]\n"
     ]
    }
   ],
   "source": [
    "bg_on_manifold = X.fillna(X.mean()).sample(100) # used to learn conditional dependence\n",
    "explainer = OnManifoldExplainer(f, pd.concat([bg, bg_on_manifold]), nruns=100, sigma_sq=0.1)  # need > 1 run b/c conditional expectation is non-linear\n",
    "cf_on_manifold_adult = explainer.shap_values(fg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1/1 [00:00<00:00, 1088.58it/s]\n"
     ]
    }
   ],
   "source": [
    "explainer = IndExplainer(f, bg, nruns=1)\n",
    "cf_ind_adult = explainer.shap_values(fg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compare baselines"
   ]
  },
  {
   "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=\"1047pt\" height=\"239pt\"\n",
       " viewBox=\"0.00 0.00 1046.53 239.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 235)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-235 1042.53,-235 1042.53,4 -4,4\"/>\n",
       "<!-- Education&#45;Num -->\n",
       "<g id=\"node1\" class=\"node\"><title>Education&#45;Num</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"652.742\" cy=\"-105\" rx=\"92.0841\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"652.742\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">Education&#45;Num: 16.00</text>\n",
       "</g>\n",
       "<!-- predicted income &gt;50k -->\n",
       "<g id=\"node2\" class=\"node\"><title>predicted income &gt;50k</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"451.742\" cy=\"-18\" rx=\"113.98\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"451.742\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">predicted income &gt;50k: 1.91</text>\n",
       "</g>\n",
       "<!-- Education&#45;Num&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge1\" class=\"edge\"><title>Education&#45;Num&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"#ff0051\" stroke-width=\"4.95694\" d=\"M615.876,-88.4097C583,-74.507 534.691,-54.078 498.99,-38.9803\"/>\n",
       "<polygon fill=\"#ff0051\" stroke=\"#ff0051\" stroke-width=\"4.95694\" points=\"500.64,-34.9692 489.741,-35.069 497.262,-42.9588 500.64,-34.9692\"/>\n",
       "<text text-anchor=\"middle\" x=\"577.242\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#ff0051\">2.13</text>\n",
       "</g>\n",
       "<!-- Education&#45;Num noise -->\n",
       "<g id=\"node3\" class=\"node\"><title>Education&#45;Num noise</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"652.742\" cy=\"-213\" rx=\"88.2844\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"652.742\" y=\"-209.3\" font-family=\"Times,serif\" font-size=\"14.00\">Education&#45;Num noise</text>\n",
       "</g>\n",
       "<!-- Education&#45;Num noise&#45;&gt;Education&#45;Num -->\n",
       "<g id=\"edge2\" class=\"edge\"><title>Education&#45;Num noise&#45;&gt;Education&#45;Num</title>\n",
       "<path fill=\"none\" stroke=\"#ff0051\" stroke-width=\"4.83127\" d=\"M652.742,-194.969C652.742,-178.378 652.742,-152.883 652.742,-133.431\"/>\n",
       "<polygon fill=\"#ff0051\" stroke=\"#ff0051\" stroke-width=\"4.83127\" points=\"656.97,-133.341 652.742,-123.341 648.515,-133.341 656.97,-133.341\"/>\n",
       "<text text-anchor=\"middle\" x=\"665.242\" y=\"-155.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#ff0051\">2.08</text>\n",
       "</g>\n",
       "<!-- Hours per week -->\n",
       "<g id=\"node4\" class=\"node\"><title>Hours per week</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"96.7422\" cy=\"-105\" rx=\"91.784\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"96.7422\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">Hours per week: 60.00</text>\n",
       "</g>\n",
       "<!-- Hours per week&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge3\" class=\"edge\"><title>Hours per week&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"#ff0051\" stroke-width=\"1.71873\" d=\"M146.553,-89.8693C184.118,-79.4108 236.988,-65.1116 283.742,-54 312.782,-47.0984 344.823,-40.2381 373.068,-34.4439\"/>\n",
       "<polygon fill=\"#ff0051\" stroke=\"#ff0051\" stroke-width=\"1.71873\" points=\"373.791,-37.8687 382.89,-32.4406 372.392,-31.0099 373.791,-37.8687\"/>\n",
       "<text text-anchor=\"middle\" x=\"296.242\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#ff0051\">0.72</text>\n",
       "</g>\n",
       "<!-- Hours per week noise -->\n",
       "<g id=\"node5\" class=\"node\"><title>Hours per week noise</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"87.7422\" cy=\"-213\" rx=\"87.9851\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"87.7422\" y=\"-209.3\" font-family=\"Times,serif\" font-size=\"14.00\">Hours per week noise</text>\n",
       "</g>\n",
       "<!-- Hours per week noise&#45;&gt;Hours per week -->\n",
       "<g id=\"edge4\" class=\"edge\"><title>Hours per week noise&#45;&gt;Hours per week</title>\n",
       "<path fill=\"none\" stroke=\"#ff0051\" stroke-width=\"1.22583\" d=\"M89.1882,-194.969C90.5969,-178.378 92.7616,-152.883 94.4132,-133.431\"/>\n",
       "<polygon fill=\"#ff0051\" stroke=\"#ff0051\" stroke-width=\"1.22583\" points=\"97.9112,-133.601 95.2699,-123.341 90.9363,-133.009 97.9112,-133.601\"/>\n",
       "<text text-anchor=\"middle\" x=\"106.242\" y=\"-155.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#ff0051\">0.51</text>\n",
       "</g>\n",
       "<!-- Relationship -->\n",
       "<g id=\"node6\" class=\"node\"><title>Relationship</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"387.742\" cy=\"-105\" rx=\"93.6835\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"387.742\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">Relationship: &#160;Husband</text>\n",
       "</g>\n",
       "<!-- Relationship&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge5\" class=\"edge\"><title>Relationship&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"#ff0051\" stroke-width=\"5\" d=\"M400.694,-86.799C409.889,-74.5865 422.383,-57.9921 432.71,-44.2776\"/>\n",
       "<polygon fill=\"#ff0051\" stroke=\"#ff0051\" stroke-width=\"5\" points=\"436.29,-46.7958 438.81,-36.1754 429.3,-41.5326 436.29,-46.7958\"/>\n",
       "<text text-anchor=\"middle\" x=\"436.242\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#ff0051\">2.15</text>\n",
       "</g>\n",
       "<!-- Sex -->\n",
       "<g id=\"node7\" class=\"node\"><title>Sex</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"243.742\" cy=\"-213\" rx=\"50.0912\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"243.742\" y=\"-209.3\" font-family=\"Times,serif\" font-size=\"14.00\">Sex: &#160;Male</text>\n",
       "</g>\n",
       "<!-- Sex&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge6\" class=\"edge\"><title>Sex&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"#ff0051\" stroke-width=\"1.27165\" d=\"M240.771,-194.905C237.229,-168.647 234.602,-118.161 259.742,-87 260.318,-86.2866 340.783,-57.9394 397.35,-38.0746\"/>\n",
       "<polygon fill=\"#ff0051\" stroke=\"#ff0051\" stroke-width=\"1.27165\" points=\"398.631,-41.3342 406.907,-34.7191 396.313,-34.7294 398.631,-41.3342\"/>\n",
       "<text text-anchor=\"middle\" x=\"272.242\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#ff0051\">0.53</text>\n",
       "</g>\n",
       "<!-- Sex&#45;&gt;Relationship -->\n",
       "<g id=\"edge7\" class=\"edge\"><title>Sex&#45;&gt;Relationship</title>\n",
       "<path fill=\"none\" stroke=\"#ff0051\" stroke-width=\"1.54108\" d=\"M264.681,-196.587C288.999,-178.686 329.327,-149 357.153,-128.517\"/>\n",
       "<polygon fill=\"#ff0051\" stroke=\"#ff0051\" stroke-width=\"1.54108\" points=\"359.377,-131.226 365.356,-122.479 355.228,-125.588 359.377,-131.226\"/>\n",
       "<text text-anchor=\"middle\" x=\"350.242\" y=\"-155.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#ff0051\">0.65</text>\n",
       "</g>\n",
       "<!-- Relationship noise -->\n",
       "<g id=\"node8\" class=\"node\"><title>Relationship noise</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"387.742\" cy=\"-213\" rx=\"76.0865\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"387.742\" y=\"-209.3\" font-family=\"Times,serif\" font-size=\"14.00\">Relationship noise</text>\n",
       "</g>\n",
       "<!-- Relationship noise&#45;&gt;Relationship -->\n",
       "<g id=\"edge8\" class=\"edge\"><title>Relationship noise&#45;&gt;Relationship</title>\n",
       "<path fill=\"none\" stroke=\"#ff0051\" stroke-width=\"3.37734\" d=\"M387.742,-194.969C387.742,-178.378 387.742,-152.883 387.742,-133.431\"/>\n",
       "<polygon fill=\"#ff0051\" stroke=\"#ff0051\" stroke-width=\"3.37734\" points=\"391.242,-133.341 387.742,-123.341 384.242,-133.341 391.242,-133.341\"/>\n",
       "<text text-anchor=\"middle\" x=\"400.242\" y=\"-155.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#ff0051\">1.45</text>\n",
       "</g>\n",
       "<!-- Age -->\n",
       "<g id=\"node9\" class=\"node\"><title>Age</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"492.742\" cy=\"-159\" rx=\"50.0912\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"492.742\" y=\"-155.3\" font-family=\"Times,serif\" font-size=\"14.00\">Age: 40.00</text>\n",
       "</g>\n",
       "<!-- Age&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge9\" class=\"edge\"><title>Age&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"#ff0051\" stroke-width=\"1.02694\" d=\"M494.556,-141.011C495.578,-126.405 495.789,-104.923 490.742,-87 486.498,-71.9257 478.19,-56.664 470.328,-44.429\"/>\n",
       "<polygon fill=\"#ff0051\" stroke=\"#ff0051\" stroke-width=\"1.02694\" points=\"473.222,-42.4598 464.747,-36.101 467.407,-46.3566 473.222,-42.4598\"/>\n",
       "<text text-anchor=\"middle\" x=\"508.242\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#ff0051\">0.42</text>\n",
       "</g>\n",
       "<!-- Marital Status noise -->\n",
       "<g id=\"node10\" class=\"node\"><title>Marital Status noise</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"840.742\" cy=\"-213\" rx=\"81.7856\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"840.742\" y=\"-209.3\" font-family=\"Times,serif\" font-size=\"14.00\">Marital Status noise</text>\n",
       "</g>\n",
       "<!-- Marital Status -->\n",
       "<g id=\"node11\" class=\"node\"><title>Marital Status</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"900.742\" cy=\"-105\" rx=\"137.576\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"900.742\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">Marital Status: &#160;Married&#45;civ&#45;spouse</text>\n",
       "</g>\n",
       "<!-- Marital Status noise&#45;&gt;Marital Status -->\n",
       "<g id=\"edge10\" class=\"edge\"><title>Marital Status noise&#45;&gt;Marital Status</title>\n",
       "<path fill=\"none\" stroke=\"#008bfb\" stroke-width=\"1.34228\" d=\"M850.382,-194.969C860.048,-177.893 875.053,-151.384 886.172,-131.74\"/>\n",
       "<polygon fill=\"#008bfb\" stroke=\"#008bfb\" stroke-width=\"1.34228\" points=\"889.246,-133.415 891.126,-122.988 883.154,-129.967 889.246,-133.415\"/>\n",
       "<text text-anchor=\"middle\" x=\"893.242\" y=\"-155.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#008bfb\">&#45;0.56</text>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.files.Source at 0x7f7366d68d50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<AGraph <Swig Object of type 'Agraph_t *' at 0x7f7360436480>>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sample_ind = 20\n",
    "cf_flow_adult.fold_noise=False\n",
    "cf_flow_adult.draw(sample_ind, max_display=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdoAAAD9CAYAAAAbOraVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deZhcRdn+8e8NCASQRUgAgRD5+SIo7sMSQEVRQNlkF3BBkbyyKRBFEZGwiOwgii8EBVzYI0siEBIDYQkEmLAqIGvAsIY1ZE/M8/ujqsnJSc/WMz3Tk7k/19VXp+vUOadO2/JM1alTjyICMzMzq4+leroBZmZmSzIHWjMzszpyoDUzM6sjB1ozM7M6cqA1MzOro2V6ugHWWNZYY40YNGhQTzfDzKxXmTRp0usR0b/aNgdaW8SgQYNobm7u6WaYmfUqkp5vaZuHjs3MzOrIgdbMzKyOOhRoJQ2TFC28vtnKfiMkje90a2sgaTtJR1Qpv1RSt46RSpqcv6t9S+Ur5fIDurM9ZmZWf7Xco30H2KFK+dOdbEu9bAfsCZxbKj8J6Nf9zQHg55KujL6y/qV27+kWmJm1La6ty2FrCbTzI2Jil7ekm0XEMz106vHANsCuwPU91AYzM+smXX6PVtJ6km6SNCsPlX6/Sp3Fhm0lDcrDpzsVypaWdIykJyXNkTRF0qWF7TtKGivpNUnTJE2UtF1h+zBgKLB+YYj70lba8ClJ4yTNlPSWpMskrVmljXtLulDSO7lNJ0hq73c5CRgNHNtapfzdnVkqOyCff6X8eZv8eVtJN0iaIempPFy+tKQzJL0u6UVJR7WzfWZm1oVqCrSSlim/crmAG4BNgAOBo4AfAYNrbN+FwAnA1cBOpKC5QmH7h4BRwLeAPYC7gZslbZW3/wG4HHglt2Ewaci42jX1J/U2VwD2Aw4HvgCMlbRsqfrpwHTSkPRfgV/mf7fXyUCTpGpD8LW4ELgL2A14HhgB/A54P+laRgBnSdq8i85nZmbtVMvQ8erAvHKhpA8BHwU+DWwREffm8knAM8BTHTmJpI1IwfpHEXFeYdNVlX9ExO8K9ZcCbgM+lvebEBFTJL0MzGnHcPfQ/L59REzLx3wKmEgK4lcU6t4REZX6Y3PA3J30B0GbImKCpNtJvdrR7dmnDX+JiDNym6cA/wI+EhFfymX/APbJbby3vLOkIcAQgIEDB3ZBc8zMrKKWHu07wKZVXi8BmwGvVoIsQEQ8Txou7agv5vdLW6ogaV1Jf5L0IjCf9AfAdsCGNZxvM2BMJcgC5OuYDGxdqjum9PkxYN1Cu4q9/aVbON/JwNaSvlBDW8vGFf5dmZR2a6UgIhYAzwLrVNs5IoZHRFNENPXvX3VhEzMzq1Gtk6GqPhYjaS3gtSqbXiMNY3bE6sCMYuArnWspYGQ+7i9JAWYGcCIwoIPnAlib1BMsexX4QKns7dLnucDyuV2DgOcK254HBpUPGhH/kHQvqVfb2Wm577UnIuamEfyW22hmZt2nq5dgfIXqQW4AMKvweTZQvu+5WunzG8CKklZuIdh+mDRM/dWIeG/4VVKtj+y8TPW2r0nHeuQvkXr4FXNaqfsr0h8Lm1XZ1p7vyMzMGlxXzzq+H1izOOlG0kDgM6V6U4BBkoo9rO1KdSpDn99u4VyVgPpeIJO0PrBVqV57e3L3AttLeq/nLWlTUm/0rnbsD6QeZUQ0F16PtlJ3FPAw8Isqm6cAG5fKyt+RmZk1uFp6tMtI2qJK+X+Am0iB4xpJPyUFwRNYfDj5etIQ7x/y4zafBr5XrBAR/5Y0nDRbdgBwB7AqsGdEfAN4ghSMzpJ0HGkI+QTgxdK5niAF/wOAfwKvR8TkKu0/GzgYuEXSacBKwKnAo8DfWvtCOukUChO8Cq4Dfivp56Q/YPYgTfTqfer0ELiZWW9QS492FeCeKq/v5pWOdiFNDroYOIf0mMk9xQNExD9JgXUwaej0C8B3q5zrEFLw/CYpiJ8LzMzHmEO6tzmf9PjKScCvgdtLx7iaNKHqdFLAGlbtoiJiKmkC1mzSDOPzgTuBr0TE3Na/kk4ZQfpjoGw46Xp/SLqGOaQJVGZm1ouor6wCaO3T1NQUTpNnZtYxkiZFRFO1bc7eY2ZmVkcOtGZmZnXkQGtmZlZHDrRmZmZ15EBrZmZWR20GWknDJL3ewrbFUs1Zz6uWctDMzHpGVy/BaLY4dXYpZ7MaebEUawC9Zug4JzIvr/3b4zqxtrKZmfUBXR5oJX1K0jhJMyW9JekySWsWtm+ThzU3Ke03XtKIwudLJTVL+rqkf5FWbNpc0qqS/iDpJUmzJb0g6aI22jRe0ghJQyRNljRL0o2S1inVW17S6ZL+I2mOpIclfa1UZ7KksyQdl3O/tpRd6ARJTxY+ryhpnqQHCmVrSFog6SuFss9Juj1/f29Iuqi4/nKuM1DSlZLezPVukfSRNr6DL0p6V9IprdUzM7Ou1e6hY0nV6qpUpz8wHngc2I+F6wWPldRUw1KGg0hLJ55Iygz0HGlN4i2BI3PZesDn23GswcBHgKNISQZOI625XMy0M4KUSed4UrL6vYGRue0PFertR0qpdwgtf4d3Ar+UtGZEvJrbPB/4ZCEj0eeABeQlKiVtBfwjt2tPUqrAU0lZe/bMdT5ASnLwBvAD0pKUPwP+IWnDiChmSSLvsz1p7eRfR8RJ7fiuzMysi7Q30K5OSqpeTTGF3ND8vn0ltZ2kp4CJpEXxr+hg+1YHvlwMcpI2A86PiOJC/H9tx7EGAIMj4oV8nOeBuyTtEBGjJW0L7AhsExGV9ZLHSNqQlDN2r9LxdoqI2a2c7x5SYP0cKYB/jrRe82BS0B2dyx6MiOl5n1OBuyNin8L1vgiMk7RJXiP6SGBF4FMR8WauM4GUoP57pDWaKey/C2mt5F9ExJnt+J7MzKwLtXfo+B1Sz6/8+nup3mbAmGL+2Ii4lxQEtq6hfS+WepIADwE/kXRIDoLt9UAlyOZ2TSBlFarkgv0yqYc8QdIylRcwDiivXzmujSBLRMwAHiAFU0i97jtIPd1i2Z0AklYgBeGrS+e/i/RHzmcL7RwLTCvUeZf0B0+5nXsA1wBDWwuyeUi9WVLz1KlTW7ssMzProPYG2vmlHKvNEdFMGr4sWht4tcr+rwIfqKF91Y51GGlo9ZfAvyU9Jekb7ThWOVVfpWzt/O81gLVIQa34GkYanm6rXdXcCXwuT+LaPH+ulL0f+FT+DGl4eGng96XzzwHeV2jDGsA+Vdr5xSrt3AV4kzRs3KKIGB4RTRHR1L9//3ZempmZtUdXP97zMmmItmxNFg4xV3qC5RnEqwHl53UXSy0UEW+TUsf9UNIngKOByyQ9EhGPtdK2au0akNsMKSC9CHy9lWO02K4W3Eka6t2WlID+IeC/wJmkwLg0C5PKv52PO4w0xFz2UqGdI0lpAcveLX0+nHRPeoykL0RE+Q8jMzOrs64OtPcCB0t6f0S8CyBpU9KkpkpAmZLfNyYNrSJpPWAj4KmOnCwiHpH0E2D/vH9rgfYzkgYW7tFuRQq09+Xt40j3mKdHRLX8sLW4kzRh7GfAhIhYIOlRYFY+1xM5Dy4RMUPSROAjEXFiK8ccR5qk9a9qE59KpgHbk3L03iLpS8VhfTMzq7+uDrRnAweT/qN+GgtnHT8K/A0gIqbk1aROkjSTNHz9c1JPrU2S7iINhf6T1AM8CJjBwoDZkqnAjZKOZ+Gs4wciYnTePha4hTRD+jTSrOKVScO7y0fEMe1pX1FEvCnpMdK92GNy2YI8eWlHoPxY0tGkiU8LSBOo3gUG5rrHRsSTpO/4m8Ctkn5L6oWvCXwBuCsiFplwFhFv5MeH7gT+nid/zezotXSKFw0wsz6sS5+jzb2zL5KGh68gzYC9E/hK6dGefYEXSLOFTyE9vvPvdp7mHuAAUiC6mnTP8qsRMaW1nYC7c3vOBf5ICtTvDRNHRAC7AxcDR5CC7oWkCUp3lQ/WAZV7sHdUKVvkuBFxFyko9wf+AowiBd//kO8LR8TrwBbAE8A5wBjSI1CrAI9Ua0BEvEwavh4EXNuIC3+YmS2plOLLkk3SeOD1iNizp9vS6JqamqK52ctXm5l1hKRJEVF+8gPoRUswmpmZ9UYOtGZmZnXUJ7L3RMQ2Pd0GMzPrm9yjNTMzqyMHWjMzszpyoDUzM6ujhrlHK2kYKT1dxatAM/DziKj6fGgLx7kU2KSladYt7DOAlPLu0oiYXCjfBrgN+HjOnNNt2nvuWq6322n3nm6BNQIvXGJ9VMME2uwdYIf870GkhSzGStq4khKuTgaQgvx4UqahigdIC1Y8U8dzd9ZJQL+eboSZmVXXaIF2fkRMzP+eKGkyaSWoHYDLu7sxeV3giW1W7EER0ch/BJiZ9XmNfo/24fz+Xvo3Sd+X9C9JcyQ9L+no1g4gaW1JF0t6VtIsSU9KOrmyDKGkQaS1mAFukxSSIm/bJn/epHC8FSSdJ+kVSbMl3S9pu9I5x0saIWk/SU9LmibpZknrluodk7fPlvSqpNGS1ipdwhqSrpE0PV/DIaVjXJrXjq58PiC3eVNJdxauebfWviczM6uPRg+0A/P7cwA5U8//kfLR7pT/fZKkw1o5xhqkhAVHkXrGZwDfBX6bt79Myv4DcChpqHhwK8e7KO//K2A30jrEN0oqJ7bfnJQ7dygwBPgMMLyyUdK3SckUziZl2DkYeBpYscr5Hs7nGg+cL2kz2nYVcANp/eZHgWskfbId+5mZWRdqtKFjJFXatD7wO1IO1xskrUy6j3pyRJyQ64yVtALwC0n/FxH/LR8vIh4Fflw4/gRStp+LJR0eEXMkVSZbPVYYuq7Wto1JCRG+GxF/ymW3kBbzP44UMCtWBnaMiLdyvbWAcyT1y+ntNgPGRMTvC/tUmy1yRUScnI8xHtiZFDzbylb0h4g4s9DGx0gZhL5R5bqGkP4YYODAgeXNZmbWCY3Wo10dmJdfTwOfBnaPiDmkXuaKpJ7ZMpUXcCspTdy61Q6o5AhJj0malY99GbAcC3vM7bUpKb/sNZWCiFiQP5d7tPdXgmxWyZW7Tn5/CPiapBMkbSZp6RbOOaZwrnmknL1Vr7XkulIbbyAF98VExPCIaIqIpv79+7fj0GZm1l6NFmjfIQWzLYD/BZYFLpe0FGkIGFKe2HmF1225fD2qOwI4kxR4diUFm0PztuU72L61SYnhy/lcXwVWkLRcoeztUp1KmsDKOS8mDR3vDdwLvJrvHZcDbrXjtKfdr1X5vHY79jMzsy7UaEPH8yOiMrHn3twD/TOwFwsTw+9Ezs1a0lI+272AERFxbKVA0kdrbN/LwEqSVigF2zWBmbnn3S65l3kOaTh5PdJ94l8BU4ALamxf0QDgjdLnl7vguGZm1gGNFmjL/gr8NL++BMwCPhgRN3bgGP2AcgDcv/S53Ntsyf1AAHuS/gBAkvLnmpPDR8R/gFMlfReo9Y+Ast2AxwHyiMCutH1ftz68UIGZ9WENHWgjIiSdQrqn+llgGPAbSesDd5CGvjcEvhgRLT2+Mhb4oaR7SQtP7A98uFTnBVIQ/46kd4B5hZ51sT2PS7oC+J2k9+fjHQRsRJo13G6SLiT10ieShsy/CPwP6Y+KrvB9SXOBfwLfJ13zvl10bDMza6eGDrTZVaQAe3REbC/pJeBI0mMzs4Enc52WnAj0B07On68FfgiMqlSIiNmSDiLNar4deB9p0lM1BwGnAb8EViU9OrNTRHS0R3tPPtb/knrSTwMHRcT1HTxOS75BGpo+mfQI0j4R8WAXHdvMzNpJEdHTbbAuJOkA4BLg/RExvaP7NzU1RXPzYp15MzNrhaRJLa0532izjs3MzJYoDrRmZmZ15EC7hImISyNCtQwbm5lZ13OgNTMzqyMHWjMzszqq2+M9koaRHpd5OiL+p8r2p0jPdp4QEcO64HyDSFl+do6Iv+eyo4H7ImJ8DccbD7weEXu2Uud9wOHA94ANgJmkZ2uvi4hTc50Ngf2AcyOivJxiW20YABwCXBoRkzt6DQ1Du/d0C5ZcXgzErOHVu0c7G/iQpEWmPEvaFBiUt3eVl0mJB4rPsx4NbNOF5yj7Hek53ctIS0MOIT2Hu3OhzoakPzhWreH4A/K+gzrVSjMz6zH1XrBiBvAAafGE4sOZ3yBl3flsV5xE0vIRMZu0ylK3yOn5vgscGxFnFDZdm5dlNDMz65Z7tFcCe1eCT37fO5cvQtJgSSMlvSxphqSHJO1fqnOApMip5cbnxAM/kTQol++U600mpd07PpeHpG3ytqGS7pf0jqRXJY2SVF6WsS0rklaQeqW8IfIqIPl8lRWonsttmJy3rS3pYknPSpol6cmcvWfZvH0QadUpgNsq11D6DlYqfTeTJZ1Z+Ly1pDslTcuvhyTt1cHrNDOzTuiOQHstKbtNJV/r50hLIla7ubQ+MAE4kDT8+jfgEknV1ui9ghTEvgb8vcr23UhrCP+RNKQ8mNS7hpTP9XekhfYPApYG7pa0SnsvKiKmkpY2HCZp97z2cdkDLEw6v3tuQ2VN5jVIax0fBewAnEHqIf82b3+ZhckPDi1cQ7tIWpn0vTwL7EFKfPAXahvCNjOzGtV9reOIeFvSaNJw8Z35fXREvFMeYY2I93q5ued7BykoHkQKrEXnRcRvCvUHlY71oKT5wJSImFjadmRhv6VJiQdeIwXeP3fg8g4g9cz/BiyQ9GD+fF5EzI2IaZIq6fseLE5oiohHWRiEkTSBNNR+saTDI2KOpEfy5sfK19AOGwKrAIdFxLu5bEy1ipKGkO4vM3DgwA6exszMWtNdj/dcCeyZE6PvSZVhYwBJq0k6T9LzLEzsPoQUNMo6kiqvfJ4tJI2V9AYwnzRbeKUWztOiiLgV+H+kPx4uJg1VnwHcmlPTtdYGSTpC0mN5+HseaVLVckBXRLtngOnA5ZJ2ldRiTzYihkdEU0Q09e/fvwtObWZmFd0VaEeSAtmvSPc2R7VQ71JgH1Kw2g7YlBTAquWJrZb8vU2SBpJ6diJlztkqn+e1Fs7Tqoh4NyKuioiDSI/4nJSPuXPre3IEcCZwHaknvRlpiJha2lGlXW8BXyHdR74amCrpRkkbdPbYZmbWft2SJi8iZkj6Oym93TURMaNcR9LypEdkDo2ICwrlLf0xUGvaoR2AFYBdK+2QtAzwgRqPt7BBKX/uGcBxpBy1N7RSfS9gREQcWymQ1N6k75XHopYtla9Was9EYAdJ/YAvA2cDlwNbtPM8ZmbWSd2Zj/b/SMOiF7SwfTlSD3tOpSBPMNqF2oPqXBbvHfYDFpCGjCv2poPfRV6sYsUqi1BUFueo9Ljn5vdq7ZhTKtu/9Lmlfafk941Jk8eQtDmwcrW2RsQsYJSkTYBjqtWpKy+qYGZ9WLcF2rw60/hWtr8j6X7gl5KmkYLhz0gzh6sGkHZ4AtgxT8aaDvyb9Pzu0qTZzH8EPkaalNShVZtIE42elPQn4Lbczo+QAtmLpCFh8jkB/lfSlcDMPBFqLPBDSfeS7qfuT1opq+gFYBbwHUnvAPMiohm4L5/jPEnHkXrjRwPTKjtK2pG0YtX1+TjrkIbKb+3gdZqZWSc02lrH+5EeR/kz8BvSbN6OzAIu+wlpJu+NwP3AZ3OQOwDYnPT4y36kYdx3OnjsacDppEU3LgJGk4LdLcAWEfEOQEQ8Twrku5N6n5X70yeSZlKfnN/nAj8sniAvwnFQPsft+RqIiLmkx4QWACOAocDBwFuF3Z8mjQScQronfXpu4/c6eJ1mZtYJymsrmAHQ1NQUzc3NbVc0M7P3SJoUEU3VtjVaj9bMzGyJ4kBrZmZWRw60ZmZmdeRAa2ZmVkcOtGZmZnXUnQtWWDtIOgA4nLTu8nxgMnBbRBzVg83qHO3e0y1ofF7Uw2yJ5R5tA5F0DPAH0rO4uwPfJi3juEtPtsvMzGrnHm1jOQy4MCJ+XigbJemEnmqQmZl1jnu0jWVV4JVyYZRWFZG0vKTTJf1H0hxJD0v6WmH7TyTNLiYpkLSZpPmSDqrrFZiZ2SIcaBvLA8Dhkr4jafVW6o0gLSN5Cikd3/3ASEmfytvPApqBP0laJmdG+hMwJiIuqlvrzcxsMR46biyHkpIAXAqEpMdJ6z2fGRHTACRtC+wIbBMRt+f9xkjaEDgW2CsiFuRJVQ+RkhysBqwJbFvtpJKGAEMABg7sipzzZmZW4R5tA4mIR0ip73YBfk9KTn8c0CxppVzty6Th5Qm5t7pMzqc7DmgqHOtp4Kd5/x8Bh0fESy2cd3hENEVEU//+/et0dWZmfZMDbYOJiDkRMSoiDouIjwLfJ+W4PTBXWQNYC5hXeg0D1isd7m/5/U3gmjo33czMqvDQcYOLiD9KOh3YKBe9ScpF+/V27H4BKRftAFIg/nmrtc3MrMs50DYQSQMi4rVSWX9SkvlXc9E4Uv7Z6RHxRCvH+jawE7AV8Ang/yRdHxH31aXxrfFiDGbWhznQNpZHJd1AStT+GrA+KWn8TNKsYYCxpAUtxko6DfgXsDLwKWD5iDhG0rrAb4AzIuJe4F5Je5BmIX86J5Q3M7Nu4Hu0jeVEYBBwHinYnkQKpJtFxHPw3jO1uwMXA0eQgu6FwGDgrnycPwBTgOMLxz6QdG/3V/W+CDMzW0iltRCsj2tqaorm5uaeboaZWa8iaVJENFXb5h6tmZlZHTnQmpmZ1ZEDrZmZWR050JqZmdWRA62ZmVkd+TnaKvIzp4cCnwH6Ac8Dfyct7l91veA6t2cI8FpEXN/d5+4S2r2nW9AyL6ZhZnXmHm2JpLOAq4FngW8B2wHnkDLfnN9DzRpC+5ZcNDOzBuMebYGknYGjgAMj4uLCptslDScF3YYlqV9EzOrpdpiZ2ULu0S7qSOCBUpAFICL+GxE3A0haQ9KfJL0haaak8ZIWeVBZUkg6rFQ2TNLrhc8H5HoflzRW0gxJT0gLx1oljQc+C3wn142caxZJkyWdJek4SVOAaZK+JmmBpA+Vzv2hXL5rZ78kMzNrPwfaTNL7gC2B0e2ofj2wPWkd4n1I3+Ntkj5c4+kvB0YCuwFPAVfm9YoBDgGeAG4iLbM4GLixsO9+wBdyvX1ISzK+BHyndI4DSOsn34iZmXUbDx0vtDqwHCmtXIsk7UDKiLNNRNyey24FJgM/Af63hnOfU+lFS5pEytSzE3BBRDwmaQYwNSImtrD/TsVEAZIuJfWAT4iIkCRS4P1rRMyvck1DSPeBGThwYA3NNzOzlrhHu7i2Fn/ejDQD+Pb3doiYQZqVvHWN5xxTONYbpJ7nui1XX8S4Ktl4LiZl/tkmf/5i/nxJtQNExPCIaIqIpv79+3ek3WZm1gYH2oXeAOYAbXXp1iYFwrJXgQ/UeO63S5/nAsu3c99XywUR8SwwHvhuLvoucF9E/KvG9pmZWY0caLOImAdMIN17bc3LwIAq5WsCbxY+zwGWLdVZreYGtqylHvgfgD0krUNKq1e1N2tmZvXlQLuoc4EmSeWJREhaKt+fvRcYIOnzhW0rADuyMB8spHywGxf3Jz2LW4uO9HArrs37XUn63/nKGs9tZmad4MlQBRExStLZwB8lbQXcAEwHNgJ+AEyOiN0k3Q1cJelnpCHnH5NWkDqjcLjrgEMlPUha/OL7wMo1Nu0JYHtJ2+fzPZfv5bZ2LbMlXUZa4eqKiCgPT3cfr75kZn2Ye7QlETGU9JjM/5AeuxkLDAXGAQfnal/P5ecC1wACvhQRTxcOdULedjJwKfAQtQ/fngw8Tlqx6n5g53buV1mycbHngs3MrHsooq1JttZbSTod2BvYICIWtGefpqamaG5urm/DzMyWMJImRURTtW0eOl4CSfoI8FFSD/yE9gZZMzPreg60S6YLgc1Jq02d18NtMTPr0xxol0ARsU1Pt8HMzBJPhjIzM6sjB1ozM7M6cqBtMJKey6nwas0EZGZmDcT3aBuIpMHAoPxxX+CknmtNF1qYXrd7eaEMM2sA7tE2ln2BGaRlHvft4baYmVkXcKBtEJKWJi0uMZK0ktPGkj5ZqrONpEckzZZ0v6TNJL0uaVip3q6SmnO9VySdnhPbm5lZN3OgbRxfJGUAuhIYAcyj0KvNWXhuIqXo25P0rOxlpDWWKdTbm5RQ4D5gF9JSkEOAX9f9CszMbDG+R9s49iXlpR0dEXMljQG+IemYSOtkHgHMBHaOiFkAkqYBV1UOIEmkxAZ/johDCuVzgPMl/bpaMgJJQ0jBmIED20rHa2ZmHeEebQOQtCwpZ+x1ETE3F18JrA8Mzp83BcZWgmw2snSoDUmJ66+WtEzlBdxKSrO3SbXzR8TwiGiKiKb+/ft3zUWZmRngQNsovgqsCtwkaVVJqwLjScnjK8PHawFTiztFxGxSGr+KNfL7TaSh58rruVy+Xj0ab2ZmLfPQcWOoBNNrqmzbS9IRwCvAIt1NScsDKxWK3szvQ4AHqxzruSplZmZWRw60PUzSiqT8slcAw0ubPw2cDXyJlIf2u5L6FYaPdynV/zfwIjAoIi6qX6vNzKy9HGh73q7ACsBvIuLe4gZJE4BjST3e44BDgVGSziENJf+MNEFqAUBELJA0FPiLpJWBm4G5wAakZPV7RsTMbrmqIi8cYWZ9mO/R9rx9gafKQRYgIuYBV5MmSr0O7AgMID2+czjwPWBpYFphn6tIwftTpKHoa4FDgAdIQdfMzLqRe7Q9LCJ2bmP7IaRACXAb8InKNklbA8sBD5f2uZnUmzUzsx7mQNuLSDqNNMnpFeAjpOHkR4Dbe7JdZmbWMgfa3mU50oIUawLvAmOAoyJiQY+2yszMWuRA24tExBGkFaLMzKyX8GQoMzOzOnKgNTMzqyMH2m4iaZikKLxekfR3SZ9oe28zM+utfI+2e70D7JD/PQg4ERgraeOIeLPFvXo77d595/LiGGbWYBxou9f8iJiY/z1R0mTgHlLwvbzHWmVmZnXjoeOeVVloYj0ASYMljZT0sqQZkh6StH95J0nrS7pC0uuSZkp6RNJ+he3LSzpd0n8kzZH0sKSvddM1mZlZgXu0PauSZb2SVWd9YAJwATAb2Aq4RNKCiLgCQNIAUi94JvBj4D+kPLPFFHgjgM2A44FngL2BkZKaIuKhul6RmZktwoG2mwweBawAABmcSURBVOVE7JCC6u+Ah4AbACLiykI9AXcA6wIHkbL7ABwJrAJ8NiJezmXjCvttS1oTeZuIqKwYNUbShqQEBXtVadMQUmo9Bg4cWN5sZmad4EDbvVYnJWKveAPYNCLmAEhaDTiBlBRgHVLCAEip7yq+BIwuBNmyL5OWaJxQCOqQgvEB1XaIiOHkFH1NTU3RgesxM7M2ONB2r3dIgXBp4JPAmcDlkrbKyyheCmwBnAQ8RsrKczAp8FasTspN25I1SCn05lXZ9t9Ott/MzDrIgbZ7zY+I5vzveyXNAv4M7CXpBmAn4NCIuKCyg6TyhLU3gLVbOcebpB7w17uu2WZmVivPOu5ZfwX+BfyUlDBgKWBOZaOk9wO7lPYZB2wvac0WjjmO1KOdHhHN5VeXX4GZmbXKPdoeFBEh6RTgMqCJNCT8S0nTgAXAz0jDzSsXdjsH+DZwp6RfkWYdbwysGBGnA2OBW0gLYZxGCuQrkxLBLx8Rx3TLxRV5EQkz68Pco+15VwFPAUcD+wHPkoaTfwP8Lf/7PRExlfTYz4PAucDfSTOGX8jbA9gduJiU6ecW4EJgMHBX3a/GzMwWofTfZbOkqakpmps9wmxm1hGSJkVEU7Vt7tGamZnVkQOtmZlZHTnQmpmZ1ZEDrZmZWR050JqZmdVRrwq0Sg6QdK+k6ZKmSbpdUnlRh15H0tGStqlSHpIO64EmmZlZF+htC1b8npTJ5vfAL0jt/wZwg6SfRcRpPdm4TjqalM1nfKl8MAvT6DU8nTl/sbL4cW/7mZmZdZ1e819ASV8HfgAcXFwLGLhZ0ivAKZLGRsQDPdPC+oiIiT3dBjMzq11vGjr+EfA0cFGVbacA7wLvDbFK2k3SfZJmSXpD0k2S1i9s/4SkUZLezsPQ90n6St52QB6yXal4EkmTJZ1Z+Dxe0ghJQ/K2WZJulLROab9TJT2azzNF0mWS1ioel5SV5/h83qgMI1cbOpZ0mKSnJM2R9LSkI0vbh0l6XdKnJU2UNFPSg5I+165v2szMukyvCLQ5r+pgYFRELJbqLSLeAW4DPp/rfwu4FngG2Bv4LvAk0D9v3wiYQMqC8wNgN+A6YL0amjcYOBw4CjgQ+ARwfanOANIfAzuSlkXcALi1kJlnN9Kaxn/MxxsMVO2ZSzoI+C0wEtgZuAY4S9LPSlVXAP5EWn5xD1KygmslrVDDNZqZWY16y9DxGqTsNs+3Uud5YIccvE4FrouIfQvbRxb+fTwpsH0uImblsrE1tm0AMDgiXgCQ9Dxwl6QdImI0QER8r1JZ0tLAPcAUYGvgjoh4UNJ8YEprQ8X52oYBl0bE0Fw8RtIqwDGSzo2I2bm8H3BERNya932ZtD7y54HRpeMOIa2XzMCBA2v8GszMrJpe0aPtoI8AHwQuaaXOl4CrCkG2Mx6oBFmAiJgAvAZsVimT9FVJd0t6B5hPCrIAG3bwXOuSru2aUvlVpAw9Hy+UzWXRiVWPFY6xiIgYHhFNEdHUv3//DjbJzMxa01sC7eukoc/1W6mzPinh+er588ut1F29je0d8VoLZWsDSNqU1JueAnyLNCy8Ra63fAfPVUn4/mqpvPL5A4WydyNiQeVDRMyt8ZxmZtYJvSLQRsR80nDrjoX7mu+RtDKwDXAH8EYuXrtcr+CNNrZXhl+XLZWvVqXugBbKKoF8N2AqsE9EjMxDw6+0cu7WVI5ZPmclCfybNR7XzMzqpFcE2uw3pKHW71fZ9jPS0OnvgH+TerbfaeVY44C9JbXUu6sM7W5cKZC0OYsmYK/4jKSBhXpbkQLhfbmoHzAvFs1HuH+V48yl7d7mFOAlYK9S+d7ANODRNvY3M7Nu1lsmQxER10u6ADhf0kdJCc+XAfYBDgCOqTxDK+lo4DJJlwFXAEG6L3tFRDQDJwD3A3dIOovUw/008EZEXEwKki8C50k6jjQkezQpmJVNBW6UdDwpUJ5Gum9bmXA0FjhC0rnAKGBL4JtVjvMEqcc+GpgO/Dsi3i19BwskDQMulPRGPvYXgIOBnxcmQvUYL05hZrao3tSjBTiEtDLUYOAG0qSgDYBdI+LUSqWIuJz0SMtGwAjgz/nfU/P2f5Nm/L4O/IH0aM+e5FnN+X7mbsCCvP9QUjB7q0qb7gbOB84lPZ7zT+DrhbbcBPw0t2ckKTDuVOU4PwFmADeS/gj4bLUvICIuIj1TvBvpj419gaHF6zczs8ahRUc0rSMkjQdej4g9e7otXaWpqSmam5t7uhlmZr2KpEkR0VRtW2/r0ZqZmfUqDrRmZmZ15JkrnRAR2/R0G8zMrLG5R2tmZlZHDrRmZmZ1tMQEWkkH5pRy65bKT8vl3yyVfyWXb9mJcw7Kx6j2uE63HcPMzBrXknSP9u78viVwdaF8S2Bmfv9rqXwOMKlbWtdH6Mz5i5V5EQsz68uWmB4taWWlN0kBFABJ7wOaSAtWlHuuWwKTImJOLSdrZflGMzOz9ywxgTavJXwPiwbUT+f33wObSHo/vJfXdXNS8nck7S3pUUlzJP1H0q9ysnny9gPy8O5mksZLmkVayWkxkr4o6V1JpxTK1pd0haTXJc2U9Iik/Vq6FknflnSXpDclvSXpNklNpTofkzQ615kh6XFJhxa2by3pTknT8ushSeU1ks3MrM6WtDG9u4FhkvrlXLODSUPD/yQlet8c+AfwMWAVYIKk7Uj5XP9MCp6fAE4ipdL7Qen4V5CC9gnA2+WTS9qetJzjryPipFw2gPQHwEzgx8B/gE2A9Vq5jkG5Pc+QMgjtC9wp6WMR8WyuMwp4nLRu8hxSHt6V8zlXJi3PeANwIiBSrtpVWzmnmZnVwZIWaCcA7wM2JaXM2xK4JyJC0sT8+R8s7PXeTQpY4yOiku1ntCSAX0s6OSKmFI5/XkT8pvJB0qDCv3ch3Rv+RUScWdjnSFJQ/2xEVNLcjWvtIiLixMJxlyIlD9iMFFRPlLQG8CHSGs+VjD3FY26Yz3lYITHBmNbOaWZm9bHEDB1n9wPzWRhItyT1JgEmlsqfIt3T/QwpOUHRVaTvZnCp/MYWzrtHPsbQUpCFlDVodCHItknSxpKuk/Qq8F9gHqnHumGu8iapZ3yBpH1yr7noGVIGoMsl7Sqp1Z6spCGSmiU1T506tb3NNDOzdliiAm1EzAQeArbMj/msy8LZyPcAWyh1V7ck9X7XIPWAXy0dqvL5Ay2Ul+1CCn7XVdm2OgsTtrcp30ceQxpaPgr4HKmH/jA5X21ELAC2IyWQvxh4Jd+P/XTe/hbwlXxtVwNTJd0oaYNq54yI4RHRFBFN/fv3b29TzcysHZaoQJtNIPVEtwQmR8Qrufw+4P3ANsCHSQH4dVJvsdwjXDO/v1kqbynV0eGk/LVjJK1e2vYGsHYH2j+Y9AfCNyPisoi4K+fQXWWRhkQ8ERF7kO67fpkUhG/MQ81ExMSI2CFv353UG768A+0wM7MusCQG2rtJPdXvsHDYmIiYBvyLNCEJYEJE/Jc0Wao8G3dvUi7ae2ifacD2+d+35MlIFeOA7SWtufhuVfXL7+89dpQX1RhUrXJEzIuIW4GzSQF91dL2WRExitTz/Wg722BmZl1kSZsMBQuHir9KSpBedA8pcfxbpBm7AMeTguMlwJWk2bknAReVJkK1KiLekPQV4E7g75J2yEPZ5wDfJs0a/hXp3urGwIoRcXqVQ00k3V+9SNLppN7tMFKPGQBJnwDOJN1LfhZYjZRc/uGIeFPSjsD3gOuBF4B1gP8Fbm3v9dTKi1OYmS1qievR5uD4AumRlnKP9J5KeX7ulogYA3yDtLDFKOAI4CzgsBrO/TKwLan3ea2kZSNiKrAV8CBwLumxmyG5jdWO8Sqph70W6fGcI0iPGT1dqPYK6X7xscDNpEeOHifdKybXDeAU0v3e04HRpOBrZmbdSDnemAHQ1NQUzc3NPd0MM7NeRdKkiGiqtm2J69GamZk1EgdaMzOzOnKgNTMzqyMHWjMzszpyoDUzM6ujhgu0kvaQdKukt3PauiclnS3pg3U41zY5/d0m+fOykoZJ+lSNx5ssqbzWcbV6y0k6UtL9OaXeHElPSbpA0sdqOG9I6vDjSGZmVn8NtbqApLNIz41eQlroYRppNaMfkLLV7NbFp3yAtOThM/nzsqQFLCaT1kzucpJWID3b+nHgt8AvgLmk1HkHATvSegq9agYDz3VhM2umM+cvVuZFLMysL2uY/wJK2pm0iP6BEXFxYdPtkoaTFtHvUnlZxoldfdw2/Ar4FLB5RPyrUH6bpPOBAzt6wIjo7mswM7N2aqSh4yOBB0pBFoCI+G9E3Fz5LOlUSY9Kmi5piqTLJK1V3KcyjCvpOEmv5LqXSVqlUGeRoWOgkrv1klwelZyz7TlnW3Jvdgjw+1KQrVzngoi4qFB/RUm/k/RvSTMlPSfp/NJayosNHUsaL2mEpP0kPS1pmqSbc0YjMzPrRg3Ro5X0PlK2nbPaucsA0vKCLwH9gaHArZI2ySnkKvYlLUd4EGnB/dOBP7B4EoGKL5HWAz6ZhblnKynu2nvO1nwWqAwdt8cKwNKkpRankoaUjyXlvt2+lf0ANgc+mNvZD/gNMBz4WjvPbWZmXaAhAi0pZ+tytLD+b1lEvLdmr6SlSWsYTwG2Bu4oVO0H7BgR03PdGcBfJG0cEY+zuPvz+zPl4dgOnLM1lQld/ykW5tR2740uRMT8/D4VOLhQbxnSvdi7JA2MiNa+r5VJ1/5W3nct4BxJ/SJiVjvba2ZmndRIQ8fQcr7XRUj6qqS7Jb0DzCcFPEg5V4vGVoJsdh0pqcCmHW1YB87ZHuXrHEnKizsPmFcYykbStyQ9KGl63n5XO897fyXIZo/l93XKFSUNkdQsqXnq1KkduQ4zM2tDowTaN0j5Vwe2VVHSpqTANAX4FmnG7RZ58/Kl6q8VP+S0ddPpWCL2jp6zNS/l9/K90iNIwf8HpfPuBvyZ1HveK5+zMvO6rfO+Xfo8t6X9ImJ4RDRFRFP//v3bOKyZmXVEQwwdR8Q8SRNI9x1/0Ub13Uj3K/eppLqTtH4LdQcUP+TJSCux8L5re3XknK2ZBMwkzaB+LzdsRDydj7lSqf5ewL0RcUilQNIXajivmZn1kEbp0ULK1dok6TvlDZKWkrRD/tgPmFcJeNn+LRzzK6XgtRtp2LalPHAt9fo6cs4W5R71cOBQSRu3Y5d+pJ5+UYfPa2ZmPacherQAETFK0tnAHyVtRUp6Ph3YiDSkOpmUvHwscISkc0mJ2rcEvtnCYWcBN0o6gzRcfAZwXUQ8Vq1yRMyV9Bywt6R/ArOBRzp4zrYcC2wG3CPpd8Cd+TzrAN8B/pvbTT7v+ZKOBe4lzRjetsbzdgsvTmFmtqiG+q9iRAyVdDdwGHA5qUc3mXR/9Mxc5yZJPwUOJz22cw+wE/BklUNeSXo29o+kIeORFGbxtuAH+Vz/IM2E/lAHz9nWNc6U9CXgUGA/0v3ZZUgzkccBn4yIykpVFwIbAD8i9bLH5n28QIWZWS+hRUdDlxySJgMjIuLHPd2W3qSpqSmam1saWTczs2okTYqIpmrbGukerZmZ2RLHgdbMzKyOGuoebVeKiEE93QYzMzP3aM3MzOrIgdbMzKyOHGjNzMzqqGEDraQ9JN0q6W1JcyQ9KelsSR9se+8On2uRvLSSlpU0TNKnajzeZElntlFnvKQRtRy/19Hui7/MzPqIhgy0ks4CrgaeJS3ivx1wDmlVpPPrcMoHSIkCKgtFLAscD9QUaM3MzCoabtaxpJ2Bo4ADI+LiwqbbJQ0nBd0uFRHT8GpLZmZWB43Yoz0SeKAUZAGIiP9GxM2Vz5JOlfSopOmSpki6LCc4p1BnsqQzJR0n6ZVc9zJJqxTqLDJ0TFq2EeCSXB6SBrX3nF1F0mGSnspD509LOrK0fV1JV0t6TdIsSc9IOqmw/WOSRkt6U9IMSY9LOrQebTUzs+oaqkcr6X2kBfvPaucuA4BTSHle+wNDgVslbRIRCwr19gWeJq1TvDZwOvAHUhq6ar5ESmN3MnBjLquk1mvvOTtF0kHAb4GzgVuALwJnSVouIk7N1f5MWg96CCn/7AakJAwVo4DHSQkQ5gAfAVbuqjaamVnbGirQAquTFvJ/oT2VI+J7lX9LWpq02P8UYGvgjkLVfsCOETE9150B/EXSxhHxeJVD35/fn4mIRYaUO3DOmklaChgGXBoRQ3PxmNwLP0bSuRExm5QFaN+IGJXrjC8cYw3gQ8CuEfFoLh7XwvmGkII1AwcO7IpLMDOzrBGHjiHljG2TpK9KulvSO8B8UsAD2LBUdWwlyGbXAQI27WjDOnDOzlgX+CBwTan8KlKP9OP580PAryUdIKkcId8kZQS6QNI+kga0dLKIGB4RTRHR1L9//665AjMzAxov0L5BGuJss1slaVNS2rsppJnJg4Et8uZy4vbXih9yAvbppGHkduvgOTuj0q5XS+WVzx/I7/uQktifAzwv6SFJ2wLkYeztgFeAi4FXJN0p6dNd2E4zM2tDQwXaiJgHTAC2b0f13YCpwD4RMTIP8b7SQt1FenOSViDlp325evUuOWdnFO8HF62Z398EiIgXI+IA0pD74NyWkZJWz9ufiIg9gFWBL5P+GLgxD02bmVk3aMT/4J4LNEn6TnmDpKUk7ZA/9gPmxaIJdfdv4ZhfkbRS4fNupOHplhKvzs3v5V5qR87ZGVNIk63Kk7X2BqYBjxYLI2JBDvonACsA65e2z4uIW0kTq9YmBd7uE9cu/jIz6yMabTIUETFK0tnAHyVtBdxAGubdCPgBMBkYDYwFjpB0Lml27Zak2bXVzCL15M4gBZozgOsi4rEW2jBX0nPA3pL+CcwGHungOdtjHUl7Vjn/CEnDgAslvZHP+wXgYODnETE7T4y6hTTz+EnSJLKhpF7t45I+AZxJuq/7LLAa8FPg4Yh4sxNtNjOzDmi4QAsQEUMl3Q0cBlxO6klOJt0fPTPXuUnST4HDSY/t3APsRAo6ZVeSno39I2nIeCQpaLXmB/lc/yAFsQ918JztsQWLT3gCUERcJGl54Ef5NQUYGhHn5DqzST3bHwHrATNJi25sFxGzJL1Cuqd7LGli1dvAbaRga2Zm3USLjoIueSRNBkZExI97ui29QVNTUzQ3tzSibmZm1UiaFBFN1bY14j1aMzOzJUZDDh1bz5k0adLrkp7v4sOuAbzexcdc0vg7apu/o7b5O2pdPb+f9VvasMQPHVvPk9Tc0pCKJf6O2ubvqG3+jlrXU9+Ph47NzMzqyIHWzMysjhxorTsM7+kG9AL+jtrm76ht/o5a1yPfj+/RmpmZ1ZF7tGZmZnXkQGtmZlZHDrTWKZI+KmmcpJmSXpJ0oqSl27HfKpIukfSWpHckXVbJOrQkqeX7kTRIUlR5Xdld7e5Okj4s6UJJj0j6r6Tx7dyvT/yGoLbvqC/9jiTtJWmkpBclTZc0SdK+7dhvOUlnSXpN0gxJN0oa1NXt84IVVjNJq5HWgn4M2BX4f8BZpD/gftHG7lcDGwLfBxYApwHXA5+rV3u7Wye/H4Afk9JGViypCxF8DPgaaa3u93VgvyX+N1RQ63cEfeN3dBTwHHAk6fq+BlwuaY2I+G0r+50H7Jn3mwoMA8ZK+nhEzO6y1kWEX37V9AKOAd4CVi6UHU1KcLByK/sNJqUp/HyhbLNc9uWevq4G+H4G5e9ip56+hm76npYq/HsEML4d+/SJ31Anv6M+8zsC1qhSdjnwXCv7rAvMB75dKFuHlCb1+13ZPg8dW2d8FbglIqYVyq4kZVv6Qhv7vRoRd1QKIuI+0l+kX61HQ3tIrd9PnxIRC2rYra/8hoCav6M+IyKq9dIfJGUua8l2+f29BNkR8SJwF138G3Kgtc7YCHiiWBARL5B6bBt1ZL/s8Tb2621q/X4qLsn3416WdLakfvVoZC/VV35DXaGv/o4G03oK042AKRExvVTe5b8h36O1zliNlOe27K28rZb9NuiCdjWKWr+fOcD5wBhgGrANKY/w/yPd67W+8xvqjD77O5K0LfB14HutVKv1/58d5kBr1mAi4mXgsELReEmvAr+X9MmIeLiHmma9SF/9HeVZw5cDN0TEpT3amMxDx9YZbwGrVClfLW/r6v16m668zhH5/bOdatGSo6/8hrraEv07kvQB4GbgeWD/Nqp322/IgdY64wlK9zIkrQesQPX7Zy3ul7V03623qvX7qSZK731dX/kNdbUl9nckaQXg78CypJnWM9vY5QlgPUkrlsq7/DfkQGudcTOwvaT3F8r2AWYBt7ex31qStq4USGoi3Vu7uR4N7SG1fj/V7JnfJ3VFw5YAfeU31NWWyN+RpGWAa4D/AXaIiNfasduY/L5b4TgfJD2H3aW/IScVsJrlBRkeA/5JWixgA+Bs4NyI+EWh3tPA7RFxYKHsFtL/KX7MwsUGXouIJWaxgVq/H0nDgPeTFhmYBnwe+AlwU0Ts0Z3X0B1yT+Rr+eNQYGXg+Pz5poiY2Vd/QxW1fEd96XckaThwEPAj4L7S5gcjYo6kcQARsW1hvwuBPVh0wYrVAS9Y4VfjvICPAreSemkvAycBS5fqTAYuLZWtClxCmvU3jTR5YbGHznv7q5bvB/gG0Ay8Q3p4/mngRGC5nr6eOn1Hg0hDmdVeg/r6b6jW76gv/Y7ytbf1/YyntNAHsBzpj9+pwAzgJuBDXd0+92jNzMzqyPdozczM6siB1szMrI4caM3MzOrIgdbMzKyOHGjNzMzqyIHWzMysjhxozczM6siB1szMrI7+P29XOkKKqtUoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cf_on_manifold_adult.draw(sample_ind)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdoAAAD9CAYAAAAbOraVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd7xcVbn/8c+X3oskAQRC4KcIyuXq9VAieEFRQAGRLmABlagUpSiKqIQi0kEULwQpFnoEIZSQGAy9nVCVXgIEKSGUkEIKeX5/rDVkZzNzyuTMmTnJ9/16zWsya6+999pnlGfW2muvRxGBmZmZNcYizW6AmZnZgsyB1szMrIEcaM3MzBrIgdbMzKyBHGjNzMwaaLFmN8BaS79+/WLQoEHNboaZWZ8ybty41yOif7VtDrQ2j0GDBtHe3t7sZpiZ9SmSnq+1zUPHZmZmDeRAa2Zm1kAtE2glDZUUhdcrkq6TtFE3j3ORpG6NfUoakM8/qFS+VW7Lht05Xk/o6rnruV4zM+s9rXaP9m1gu/zvQcCxwGhJG0TEGw087wDgaGAsML5Qfj8wGHimgeeeX8cBSze7EQsU7dLsFphZM8RVDTlsqwXa2RFxd/733ZLGA3eRgu8lvd2YiJgM3N1pxSaKiFb+EWBmttBrmaHjGh7K72tVCiR9V9K/Jc2Q9LykIzo6gKTVJV0g6VlJ0yU9Kel4SUvk7YOAR3L1f1aGrvO2DwzfSlpG0ll5aPtdSfdJ2qZ0zrGShkvaW9LTkiZLulHSmqV6R+bt70p6VdJISauVLqGfpCslTcnXcEDpGPMMHUvaN7d5Y0m3Fa55547+TmZm1hitHmgH5vfnACT9BPg/4O/ADvnfx0k6qINj9APeAA4j9YxPAfYDfpe3vwzsk/99IGmoeHAHxzsv7/9rYGfgReB6SVuU6m0KHAQcDgwB/gcYVtko6ZvAz4HTgW2BHwBPA8tWOd9D+VxjgbMlbdJB+youB64BdiH9kLhS0n93YT8zM+tBrTZ0jKRKm9YGfg88CFwjaQXSfdTjI+KYXGe0pGWAX0j6v4h4r3y8iHgE+HHh+HcAU4ELJB0cETMkPZw3P1oYuq7Wtg2AvYD9IuJPuewm4GHgl6SAWbECsH1EvJnrrQacIWnpiJgObAKMiog/FPapdoPg0og4Ph9jLLAjKXjeW6ud2R8j4tRCGx8FjgS+VuW6hpB+DDBw4MDyZjMzmw+t1qNdBZiVX08DnwJ2iYgZpF7msqSe2WKVF3AzsCqwZrUDKjlE0qOSpudjXwwsydwec1dtDAi4slIQEXPy53KP9r5KkM0eze9r5PcHgS9LOkbSJpIWrXHOUYVzzQKeosa1llxdauM1pOD+ARExLCLaIqKtf/+qC5uYmVmdWi3Qvk0KZpsB3wOWAC6RtAhpCBjg38wNxrOAf+bytajuEOBUUuDZiRRsDszblupm+1YHpkTEtFL5q8AykpYslL1VqjOzdM4LSEPHewD3AK/me8flgFvtOF1p92tVPq/ehf3MzKwHtdrQ8eyIqEzsuSf3QP8M7E66zwrp3uyrVfZ9osYxdweGR8RRlQJJH6+zfS8Dy0laphRsVwWm5Z53l+Re5hmk4eS1SPeJfw1MAM6ps31FA4BJpc8v98BxzcysG1qtR1v2V1IP9qekx3ymAx+OiPYqr3dqHGNpoBwA9yl9Lvc2a7kPCGC3SoEk5c+3d3o1NUTEixFxImm4vN4fAWXvzzLOIwI70fl9XTMz62Gt1qOdR0SEpBNI91Q/DQwFfitpbeBW0g+F9YDPRUStx1dGAz+UdA9p4Yl9gI+U6rxACuLfkvQ2MKvQsy625zFJlwK/l7R8Pt7+wPqkWcNdJulcUi/9btKQ+eeAj5J+VPSE70qaCfwL+C7pmvfqoWMv2Br00LqZLZxaOtBml5MC7BERsa2k/wCHkh6beRd4Mtep5VigP3B8/nwV8ENgRKVCRLwraX/SrOZbgMVJk56q2R84CfgVsBLp0ZkdIqK7Pdq78rG+R+pJPw3sHxF/7+ZxavkaaWj6eNIjSHtGxAM9dGwzM+siRUSz22A9SNK+wIXA8hExpbv7t7W1hdPkmZl1j6RxEdFWbVur36M1MzPr0xxozczMGsiBdgETERdFhOoZNjYzs57nQGtmZtZADrRmZmYN1K1AK2loJY1cldfXO9hveF4Qv9dJ2kbSIVXK50kv10ttGZ//VnuVypfL5fv2ZnvMzKzx6nmO9m1Surmyp+ezLY2yDWnlpjNL5ceRVo1qhp9Luiz8bFVr0i7NboGZ9aYGL1JTT6Cd3VEqub4iIp5p0qnHAluRlkTsqcUpzMysRfX4PVpJa0m6QdL0PFT63Sp1PjBsK2lQHj7doVC2qKQjJT0paYakCZIuKmzfXtJoSa9JmizpbknbFLYPJa0gtXZhiPuiDtrwSUljJE2T9KakiyWtWqWNe0g6V9LbuU3H5PWEu2IcMBI4qqNK+W93aqls33z+5fLnrfLnrSVdI2mqpKfycPmikk6R9LqklyQd1sX2mZlZD6or0BbzwRbywlYW2L8G2BD4DnAY8CNSLtl6nAscA1xBytpzOLBMYfs6pKUUvwHsCtwJ3Chp87z9j8AlwCu5DYNJQ8bVrqk/qbe5DLA3cDCwJSm5/BKl6icDU0hD0n8lLce4G113PNAmqdoQfD3OJSU12Bl4HhgO/B5YnnQtw4HTJG3aQ+czM7MuqmfouJKcfR6S1iFlnvkUsFlE3JPLx5EW33+qOyeRtD4pWP8oIs4qbHp/XeOI+H2h/iKk3LSfyPvdERETJL0MzOjCcPfh+X3biJicj/kUadH/XYFLC3VvjYhK/dE5YO5C+kHQqYi4Q9ItpF7tyK7s04m/RMQpuc0TSBmPPhYRn89l/wD2zG28pwfOZ2ZmXVRPj7aSnL38+g8pqfqrlSALEBHPk4ZLu+tz+f2iWhUkrSnpT5JeAmaTfgBsQ8ro012bAKMqQRYgX8d4YItS3VGlz48CaxbaVeztlxO5VxwPbCFpyzraWjam8O/KpLSbKwU59+2zwBrVdpY0RFK7pPaJEyf2QHPMzKyi3slQVR+LkbQa8FqVTa+RhjG7YxVgajHwlc61CHBtPu6vSAFmKilbz4BungtgdVJPsOxV4EOlsrdKn2eSc9lKGgQ8V9j2PDCofNCI+EdO3XcUqac5P95vT0TMTCP4tdtYpS3DgGGQkgrMZ1vMzKygp9PkvUL1IDeAlO+14l2gfN9z5dLnScCyklaoEWw/Qhqm/lJEvD/8KqneR3ZepnrbV6V7PfL/kHr4FeWk80W/Jv1Y2KTKtq78jczMrMX19Kzj+4BVi5NuJA0E/qdUbwIwSFKxh7VNqU5l6PObNc5VCajvBzKlhPCbl+rV7MmV3ANsq5TQvXK8jUm90S7nmo2ImRHRXng90kHdEcBDwC+qbJ4AbFAqK/+NzMysxdXTo11M0mZVyl8EbiAFjisl/ZQUBI/hg8PJfycN8f4xP27zKeDbxQoR8YSkYaTZsgOAW0mJ1neLiK8Bj5OC0WmSfkkaQj4GeKl0rsdJwX9f4F/A6xExvkr7Twd+ANwk6SRgOeBEUmL3v3X0B5lPJ1A9cf3VwO8k/Zz0A2ZX0kQva7QGP7xuZguXenq0KwJ3VXntl1c6+gppctAFwBmkx0zuKh4gIv5FCqyDSUOnWwL7VTnXAaTg+XVSED8TmJaPMYN0b3M26fGV44DfALeUjnEFaULVyaSANbTaRUXERNIErHdJM4zPBm4DvhgRMzv+k8yX4aQfA2XDSNf7Q9I1zCBNoDIzsz5EXgXQitra2qK9vVeXgDYz6/MkjYuItmrbnL3HzMysgRxozczMGsiB1szMrIEcaM3MzBrIgdbMzKyBHGjNzMwaqNMFK3JO14Miol+VbRcBG9aa0mzNUVhveceIuK65remDNL9LT1tVXgjEFlLu0ZqZmTVQnwm0khatkoC96eYjiYGZmS0EejzQSvqkpDGSpkl6U9LFklYtbN9KUkjasLTfWEnDC58vyjlSvyrp36SlETeVtJKkP0r6j6R3Jb0g6bxO2jRW0vCcd3W8pOmSrpe0RqneUpJOlvSipBmSHpL05VKd8ZJOk/TLnGS9Vhq/YyQ9Wfi8rKRZku4vlPWTNEfSFwtln5V0S/77TZJ0XjHRQa4zUNJlkt7I9W6S9LFO/gafk/SOpBM6qmdmZj2ry0kFJFWrq1Kd/sBY4DFgb+YuzD9aUlsdawYPIq1RfCwpBd9zpMX/PwMcmsvWAv63C8caDHwMOIyUzeckUnKDYkq74aSUdUcDzwB7ANfmtj9YqLc3KXftAdT+G94G/ErSqhHxam7zbOC/C6n/PgvMIa8FLWlz4B+5XbuRcvKeSEqPt1uu8yFSNqFJwPdJaz//DPiHpPUiopiOkLzPtqQkBb+JiOO68LcyM7Me0tVAuwowq8a2Yq7Ww/P7tpUcspKeAu4mZZ+5tJvtWwX4QjHISdoEODsiihlv/tqFYw0ABkfEC/k4zwO3S9ouIkZK2hrYHtgqIiqJCUZJWo+UnH330vF2iIh3OzjfXaTA+llSAP8sKTHCYFLQHZnLHoiIKXmfE4E7I2LPwvW+BIyRtGFOxnAosCzwyYh4I9e5AxhPStRwdrERkr5CSkrwi4g4tVpDJQ0BhgAMHDiwg0syM7Pu6urQ8duknl/5VZ7RugkwqpioPSLuIQWBLepo30ulniTAg8BPJB2Qg2BX3V8Jsrldd5DS91WSrn+B1EO+Q9JilRcwBijPqh7TSZAlIqYC95OCKaRe962knm6x7DYAScuQgvAVpfPfTvqR8+lCO0cDkwt13iH94Cm3c1fgSuDwWkE2t3VYRLRFRFv//v07uiwzM+umrgba2aVk5u0R0U4avixaHXi1yv6vAh+qo33VjnUQaWj1V8ATkp6S9LUuHKucE7dStnr+dz9gNVJQK76GkoanO2tXNbcBn82TuDbNnytlywOfzJ8hDQ8vCvyhdP4ZwOKFNvQD9qzSzs9VaedXgDdIw8ZmZtYE9SR+78jLpCHaslWZO8Rc6QmWZxCvDLxeKvtADr+IeIuUo/WHkjYCjgAulvRwRDzaQduqtWtAbjOkgPQS8NUOjlGzXTXcRhrq3RqYSeqNvwecSgqMi5J6rABv5eMOJQ0xl/2n0M5rSfl3y94pfT6YdE96lKQtI6L8w8jMzBqspwPtPcAPJC0fEe8ASNqYNKmpElAm5PcNSEOrSFoLWB94qjsni4iHJf0E2Cfv31Gg/R9JAwv3aDcnBdp78/YxpHvMUyKiWiL2etxGmjD2M+COiJgj6RFgej7X4znhPBExVdLdwMci4tgOjjmGNEnr39UmPpVMBrYFbgFukvT54rC+1eCFFcysB/V0oD0d+AHpP+onMXfW8SPA3wAiYoKkduA4SdNIw9c/J/XUOiXpdtJQ6L9IPcD9ganMDZi1TASul3Q0c2cd3x8RI/P20cBNpBnSJ5FmFa9AGt5dKiKO7Er7iiLiDUmPku7FHpnL5uTJS9sD5ceSjiBNfJpDmkD1DjAw1z0qIp4k/Y2/Dtws6XekXviqwJbA7RExz4SziJiUHx+6DbguT/6a1t1rMTOz+vToc7S5d/Y50vDwpaQZsLcBXyw92rMX8AJptvAJpMd3nujiae4C9iUFoitI9yy/FBETOtoJuDO350zgfFKgfn+YOCIC2AW4ADiEFHTPJU1Qur18sG6o3IO9tUrZPMeNiNtJQbk/8BdgBCn4vki+LxwRrwObAY8DZwCjSI9ArQg8XK0BEfEyafh6EHBVKy78YWa2oFKKLws2SWOB1yNit2a3pdW1tbVFe3t7s5thZtanSBpXa93/PrMEo5mZWV/kQGtmZtZAPT0ZqiVFxFbNboOZmS2c3KM1MzNrIAdaMzOzBlooho77Ekn7klZ0Wo+UlGA88M+IOKyJzVq4aJdmt6B5vFiHWY9zj7aFSDoS+CPpGd5dgG8C15DWLDYzsz7IPdrWchBwbkT8vFA2QtIxzWqQmZnNH/doW8tKpFR984jSqiKSlpJ0sqQXJc2Q9JCkLxe2/0TSu5I+XijbRNJsSfs39ArMzGweDrSt5X7gYEnfkrRKB/WGk5ahPAHYEbgPuFbSJ/P204B24E85Z+1SwJ9IuYLL6yubmVkDeei4tRxIyrV7ERCSHiMlYzi1knVH0takJANbRcQteb9RktYDjgJ2z4kL9iWl5TuSlIJwVdJ6xx8gaQgwBGDgwIGNuTIzs4WUe7QtJCIeJqUP/AopAbyAXwLtkpbL1b5AGl6+I/dWF5O0GCl9XlvhWE8DP837/wg4OCL+QxURMSwi2iKirX///g26OjOzhZMDbYuJiBkRMSIiDoqIjwPfBT4KfCdX6QesBswqvYYCa5UO97f8/gZwZYObbmZmVXjouMVFxPmSTiYltocUNF+ikOKvA+eQ0hEOIAXin3dY28zMepwDbQuRNCAiXiuV9Sflmn01F40BDgemRMTjHRzrm8AOwObARsD/Sfp7RNzbkMYvSLxog5n1IAfa1vKIpGtIydxfA9YGfgxMI80aBhhNWtBitKSTgH8DKwCfBJaKiCMlrQn8FjglIu4B7pG0K2kW8qci4t1evSozs4WY79G2lmOBQcBZpGB7HCmQbhIRz8H7z9TuAlwAHEIKuucCg4Hb83H+CEwAji4c+zuke7u/bvRFmJnZXCqthWALuba2tmhvb292M8zM+hRJ4yKirdo292jNzMwayIHWzMysgRxozczMGsiB1szMrIEcaM3MzBrIz9G2GEnPkR7x+Wher9h6m3Zpdgt6hhfeMGsJ7tG2EEmDSUEWYK8mNsXMzHqIA21r2QuYCtyDA62Z2QLBgbZFSFoU2AO4lrTq0waS/rtUZytJD0t6V9J9kjaR9LqkoaV6O0lqz/VekXSypMV77WLMzOx9DrSt43Ok5OyXAcNJqe/e79VKWgO4gbQG8m6kZRcvBpYuHkTSHsBVwL2kvLbHkJK6/6bhV2BmZh/gyVCtYy/gLWBkRMyUNAr4mqQj8/rGh5CSC+wYEdMBJE0GLq8cQJKAU4A/R8QBhfIZwNmSfhMRk8onljSEFIwZOHBgwy7QzGxh5B5tC5C0BClRwNURMTMXX0bK3jM4f94YGF0Jstm1pUOtBwwErpC0WOUF3AwsBWxY7fwRMSwi2iKirX///j1zUWZmBjjQtoovASsBN0haSdJKwFhgBnOHj1cDJhZ3yunuphSK+uX3G0hDz5XXc7l8rUY03szMavPQcWuoBNMrq2zbXdIhwCvAPN1NSUsByxWK3sjvQ4AHqhzruSplZmbWQA60TSZpWWBH4FJgWGnzp4DTgc8D9wH7SVq6MHz8lVL9J4CXgEERcV7jWr2A80IPZtaDHGibbydgGeC3EXFPcYOkO4CjSD3eXwIHAiMknUEaSv4ZaYLUHICImCPpcOAvklYAbgRmAusCXwV2i4hpvXJVZmYG+B5tK9gLeKocZAEiYhZwBWmi1OvA9sAA0uM7BwPfBhYFJhf2uZwUvD9JGoq+CjgAuJ8UdM3MrBe5R9tkEbFjJ9sPIAVKgH8CG1W2SdoCWBJ4qLTPjaTerJmZNZkDbR8i6STSJKdXgI+RhpMfBm5pZrvMzKw2B9q+ZUnSghSrAu8Ao4DDImJOU1tlZmY1OdD2IRFxCGmFKDMz6yM8GcrMzKyBHGjNzMwaqGFDxzl129HA0xHx0SrbnwI+AhwTEUN74HyDSCsf7RgR1+WyI4B7I2JsHccbC7weEbt1UGdx5j5msy7pmdZnSGsWn5jrrAfsDZwZEW91sw0DSDOOL4qI8d29hmbQqbOb3YT5Fj/2HRUz6zmN7tG+C6wjqa1YKGljYFDe3lNeJi3Af3uh7Ahgqx48R9nvgWNJ6ep2IC19eAtppaeK9Ug/OFaq4/gD8r6D5quVZmbWNI3+6T6VtFDC14D2QvnXSBllPt0TJ5G0VF5g/+6eOF4Xz7kMsB9wVEScUth0VU5XZ2Zm1iv3aC8D9qgEn/y+Ry6fh6TBkq6V9LKkqZIelLRPqc6+kkLSJpLGSpoO/ETSoFy+Q643HlgFODqXh6St8rbDJd0n6W1Jr0oaIekj3byuZYHFSc+0ziPnjyWfb0Qufi63YXzetrqkCyQ9K2m6pCclHZ9T5lWGwh/J+/6zcg2lv0ExoQCSxks6tfB5C0m3SZqcXw9K2r2b12lmZvOhNwLtVaTnPrfInz9LykJTbeX2tYE7gO+Qhl//Blwoaa8qdS8lBbEvA9dV2b4z8DZwPmlIeTCpdw2wJmnYdydgf9IyhndKWrGrFxURE4EXgaGSdpG0fJVq9wM/zv/eJbdh5/y5HynbzmHAdqTnY/cDfpe3vwxUfmQcWLiGLslrHV8HPAvsCuwG/IX6hrDNzKxODZ/1ERFvSRpJGi6+Lb+PjIi3yyOsEfF+Lzf3fG8lBcX9SYG16KyI+G2h/qDSsR6QNBuYEBF3l7YdWthvUWA08Bop8P65G5e3L6ln/jdgjqQH8uezImJmREyW9ESu+0BxQlNEPMLcIFxJIDAVuEDSwRExQ9LDefOj5WvogvWAFYGDIuKdXDaqWkVJQ0j3lxk4cGA3T2NmZh3prcd7LgN2k7QkqWf1gWFjAEkrSzpL0vPMTVo+hBQ0yq6vtzGSNpM0WtIkYDZptvByNc5TU0TcDPw/0o+HC0hD1acAN0vq8G+r5BBJj+bh71mkSVVLAj0R7Z4hJYW/RNJOOZl8resYFhFtEdHWv3//WtXMzKwOvRVoryUFsl+T7m2OqFHvImBPUrDaBtiYFMCWqlL31XoaImkgqWcn4HvA5vk8r9U4T4ci4p2IuDwi9ic94nNcPmaHyQJIKzydClxN6klvQhoipp52VGnXm8AXSfeRrwAmSrpe0rrze2wzM+u6XnlgMCKmSroOOBS4MiKmlutIWor0iMyBEXFOobzWj4GosznbkfK/7lRph6TFgA/Veby5DYoISaeQFvtfH7img+q7A8Mj4qhKgaSPd/FUlceiliiVr1xqz93AdpKWBr5ASiJ/CbBZF89jZmbzqTdXhvo/Uk/2nBrblyS1Z0alIE8w+sp8nHMmH+wdLk1KlF5cWWEPuvmjQ9LiNYZjK4tzVHrclRyw1doxo1S2T+lzrX0n5PcNCu3ZFFihWlsjYnpEjCCNDnQ1mJuZWQ/otSVw8upMYzvY/rak+4BfSZpMCoY/I80crhpAuuBxYPs8GWsK8ATp+d1FSbOZzwc+QZqU1K1Vm0gTjZ6U9CdSnti3SanrjgReIg0Jk88J8D1JlwHT8kSo0cAPJd1Dup+6D2mlrKIXgOnAtyS9DcyKiHbg3nyOsyT9ktQbP4JCAnhJ25NWrPp7Ps4apKHym7t5nd3iVZXMzObVamsd7016HOXPwG9Js3m7Mwu47CekmbzXA/cBn85Bbl9gU9LjL3uThnHf7uaxJwMnkxbdOA8YSQp2NwGbRcTbABHxPCmQ70J6dKlyf/pY0kzq4/P7TOCHxRPkRTj2z+e4JV8DETGT9JjQHGA4cDjwA+DNwu5Pk4bXTyDdkz45t/Hb3bxOMzObD8prK5gB0NbWFu3t7Z1XNDOz90kaFxFt1ba1Wo/WzMxsgeJAa2Zm1kAOtGZmZg3kQGtmZtZADrRmZmYN1KcCbV4feF9J90iaklO/3SJpfha1aAmSjqik8SuVh6SDmtAkMzPrAX1tdYE/kJ4r/QPwC1L7vwZcI+lnEXFSMxs3n44gpe4bWyofDDzX661ZmGmXnjlOVMsEaWYLmz4TaCV9Ffg+8IPiWsjAjZJeAU6QNDoi7q9+hL6pjvR4ZmbWQvrS0PGPSKsdnVdl2wnAO8D7Q6ySdpZ0r6TpkiZJukHS2oXtG0kaIemtPAx9r6Qv5m375iHb5YonkTRe0qmFz2MlDZc0JG+bnjPkrFHa70RJj+TzTJB0saTVisclpdg7Op83KsPI1YaOJR0k6SlJMyQ9LenQ0vahkl6X9ClJd0uaJukBSZ/t0l/azMx6TJ8ItDm7zmBgRES8V96elzv8J/C/uf43gKtIawjvAewHPAn0z9vXJy2HuDqpl7wzaW3itepo3mDgYOAw4DvARqT1hYsGkH4MbE9Kj7cu8+as3Zm0BOT5+XiDgao9c0n7A78jpR7cEbgSOE3Sz0pVlwH+BJwL7EpKYHCVpGXquEYzM6tTXxk67kfK7vN8B3WeJ6WEWwQ4Ebg6IvYqbL+28O+jSYHtsxExPZeNrrNtA4DBEfECQE5af7uk7SJiJEBEvL++sKRFgbtIGXi2AG6NiAckzQYmdDRUnK9tKHBRRByei0dJWhE4UtKZeX1kSNmBDsnJ6ZH0MvAA6cfIyNJxhwBDAAYO7Imc82ZmVtEnerTd9DHgw8CFHdT5PHB5IcjOj/srQRYgIu4gJZHfpFIm6UuS7swZeGYzN83det0815qka7uyVH45KcPRfxXKZjLvxKpHC8eYR0QMi4i2iGjr379/N5tkZmYd6SuB9nXS0OfaHdRZm5Q6bpX8+eUO6q7SyfbueK1G2eoAkjYm9aYnAN8gDQtXEq+X88x2ZvX8/mqpvPK5mLz+nYiYU/mQM/7Uc04zM5sPfSLQRsRs0nDr9oX7mu+TtAKwFXArMCkXr16uVzCpk+2V4dclSuUrV6k7oEZZJZDvDEwE9oyIa/PQ8CsdnLsjlWOWz7lqfn+jzuOamVmD9IlAm/2WNNT63SrbfkYaOv09KdH6S8C3OjjWGGAPSbV6d5Wh3Q0qBZI2pXoC+v+RNLBQb3NSILw3Fy1NSthezEe4T5XjzKTz3uYE4D+k/LlFe5Dy4z7Syf5mZtbL+spkKCLi75LOAc6W9HFS0vbFgD1JidyPrDxDK+kI4GJJF5OSqgfpvuylEdEOHENKon6rpNNIPdxPAZMi4gJSkHwJOEvSL0lDskeQglnZROB6SUeTAuVJpPu2lQlHo4FDJJ1JSvr+GeDrVY7zOKnHPhKYAjwREe+U/gZzJA0FzpU0KR97S1LS958XJkLZ/PBCE2bWg/pSjxbgANLKUIOBa0iTgtYFdoqIEyuVIuIS0iMt6wPDgT/nf0/M258gzfh9Hfgj6dGe3cizmvP9zJ2BOXn/w0nB7M0qbboTOBs4k/R4zr+ArxbacgPw09yea2bYxfEAABpuSURBVEmBcYcqx/kJMBW4nvQj4NPV/gARcR7pmeKdST829gIOL16/mZm1Ds07omndIWks8HpE7NbstvSUtra2aG9vb3YzzMz6FEnjIqKt2ra+1qM1MzPrUxxozczMGqjPTIZqRRGxVbPbYGZmrc09WjMzswZyoDUzM2sgB9oqJO0q6eacQm+GpCclnS7pw01qz5Ccj9fMzPoYB9qSvIDFFcCzpLWJtwHOALYmPS/bDEMoPJtrZmZ9hydDFUjakZxXNq8QVXGLpGGkoNuyJC3dQxmJzMysh7hHO69DScsnXlDeEBHvRcSNAJL6SfqTpEmSpkkaK2meB5UlhaSDSmVDJb1e+LxvrvdfkkZLmirpcUm7FOqMJa0S9a1cNyTtm7eNl3SapF9KmgBMlvRlSXMkrVM69zq5fKf5/SOZmVnXOdBmkhYnrUM8srO6wN+BbYEfk9ZaXgT4p6SP1Hn6S0jLM+4MPAVcJqmSN/YA0jrIN5CWnhxMWqaxYm/Sso4H5LbcREo8UE6qsC8pfd/1mJlZr/HQ8VyrAEsCL3RUSdJ2wObAVhFxSy67GRhPWq/4e3Wc+4xKL1rSOFJ+2R2AcyLiUUlTgYk5xV41OxQTCki6iNQDPiYiQpJIgfevOeVg+ZqGkO4DM3DgwPJmMzObD+7RflBniz9vArxWCbIAETGVtMD/FnWec1ThWJNIPc81a1efx5gqWXsuANYm5egF+Fz+fGG1A0TEsIhoi4i2/v37d6fdZmbWCQfauSYBM4DOunSrkwJh2aukdHr1eKv0uSu5aYvnnUdEPAuMBfbLRfsB90bEv+tsn5mZ1cmBNouIWcAdpHuvHXmZlNi9bFXgjcLnGcASpTor193A2mr1wP8I7CppDWAXavRmzcyssRxo53Um0CapPJEISYvk+7P3AAMk/W9h2zLA9sDthV0mABsU9yc9i1uP7vRwK67K+11G+p4vq/PcZmY2HzwZqiAiRkg6HThf0uak5PJTSEnjvw+Mj4idJd0JXC7pZ6Qh5x8DSwOnFA53NXCgpAdIi198F1ihzqY9Dmwradt8vufyvdyOruVdSRcDBwKXRkR5eNrMzHqBe7QlEXE46TGZj5IeuxkNHA6MAX6Qq301l58JXAkI+HxEPF041DF52/HARcCD1D98ezzwGGnFqvuAHbu439/z+weeCzYzs96hiM4m2VpfJelkYA9g3YiY05V92traor29vbENMzNbwEgaFxFt1bZ56HgBJOljwMdJPfBjuhpkzcys5znQLpjOBTYlrTZ1VpPbYma2UHOgXQBFxFbNboOZmSWeDGVmZtZADrRmZmYN1LKBVtKukm6W9JakGZKelHS6pA834Fxb5fRzG+bPS+SUdp+s83jjJZ3aSZ2xkobXc3wzM+s7WjLQSjqN9Mzos8A3SAnXzyCtrHR2A055Pyn93DP58xLA0UBdgdbMzKyi5SZDSdoROAz4TikB+y2ShpGCbo+KiMlArRR0ZmZmdWvFHu2hwP2lIAtARLwXETdWPks6UdIjkqZImiDpYkmrFfepDONK+qWkV3LdiyWtWKgzz9Ax8E5+vzCXh6RBXT1nT5F0kKSn8tD505IOLW1fU9IVkl6TNF3SM5KOK2z/hKSRkt6QNFXSY5IObERbzcysupbq0UpaHPgMcFoXdxkAnAD8B+hPWirxZkkblhZp2At4GtiflObuZFJ2m91rHPfzwM2kpQ+vz2Uvd/Oc80XS/sDvgNOBm0g5ZU+TtGREnJir/Zm0xvIQUqq9dUnrMleMIC3d+HVSNqGPUf96y2ZmVoeWCrTAKsCSwAtdqRwR3678W9KiwF2krDlbALcWqi4NbB8RU3LdqcBfJG0QEY9VOfR9+f2ZiJhnSLkb56xbzvQzFLgor70MMCr3wo+UdGZO9r4JsFdEjMh1xhaO0Q9YB9gpIh7JxWNqnG8IKVgzcGBn6XjNzKw7WnHoGGrnWJ2HpC9JulPS28BsUsADWK9UdXQlyGZXkxIBbNzdhnXjnPNjTeDDpKQERZeTeqT/lT8/CPxG0r6SyhHyDeBF4BxJe0qqlkMXgIgYFhFtEdHWv3//nrkCMzMDWi/QTiINcXbarZK0MWmJwQmkmcmDgc3y5nLu1teKHyJiGin93erdaVw3zzk/Ku16tVRe+fyh/L4n0E6akf28pAclbQ2Qh7G3AV4hZe95RdJtkj7Vg+00M7NOtFSgjYhZwB3Atl2ovjMwEdgzIq7NQ7yv1Kg7T28uJ2pfjrn3XbuqO+ecH8X7wUWr5vc3ACLipYjYlzTkPji35VpJq+Ttj0fErsBKwBdIPwauz0PTZmbWC1rxP7hnAm2SvlXeIGkRSdvlj0sDs2LePH/71DjmFyUtV/i8M2l4ulY+uJn5vdxL7c4558cE0mSr8mStPYDJwCPFwoiYk4P+McAywNql7bMi4mbSxKrVSYHXzMx6QatNhiIiRkg6HThf0ubANaRh3vWB7wPjgZGkxOuHSDqTNLv2M6TZtdVMJ/XkTiEFmlOAqyPi0RptmCnpOWAPSf8C3gUe7uY5u2INSbtVOf9wSUOBcyVNyufdkpT27ucR8W6eGHUTaebxk6RJZIeTerWPSdoIOJV0X/dZYGXgp8BDEfHGfLTZzMy6oeUCLUBEHC7pTuAg4BJST3I86f7oqbnODZJ+ChxMemznLmAHUtApu4z0bOz5pCHja0lBqyPfz+f6BymIrdPNc3bFZnxwwhOAIuI8SUsBP8qvCcDhEXFGrvMuqWf7I2AtYBpp0Y1tImK6pFdI93SPIk2segv4JynYmplZL9G8o6ALHknjgeER8eNmt6UvaGtri/b2WiPqZmZWjaRxEdFWbVsr3qM1MzNbYDjQmpmZNVBL3qPtSRExqNltMDOzhZd7tGZmZg3kQGtmZtZALRdoJe0q6WZJb+X0cE9KOl3ShxtwrnnS40laQtJQSXUlfK+k5OtCvSUlHSrpPknv5Ot8StI5kj5Rx3lD0kH1tNnMzBqrpe7RSjoNOAS4kLR+72Tg46RnWtchrejUk+4nLV34TP68BHA06ZndB3v4XMD7yz+OIiUG+B3wC9JKVBuSns3dnvRcbHcMBp7rwWaamVkPaZlAK2lH4DDgO6Wk77dIGkZaIL9HRcRk0iIPvenXwCeBTSPi34Xyf0o6G/hOdw9YTuVnZmato5WGjg8F7i8FWQAi4r2IuLHyWdKJkh6RNEXSBEkXS1qtuE9lGFfSLyW9kutenJcurNSZZ+iYtHoUwIW5PCQN6uo5O5N7s0OAP5SCbOU650TEeYX6y0r6vaQnJE2T9JyksyWtUDruPEPHksZKGi5pb0lPS5os6UZJa3anvWZmNv9aokcraXHSusGndXGXAcAJpIX3+5PW+L1Z0oY5PVzFXsDTpCHZ1YGTgT/ywcX6Kz4P3AwcD1yfy4qZdLpyzo58mrTo/6gu1l8GWJS0jOJE0pDyUaRlGzvLcLQpaenFw0lLWP4WGAZ8uYvnNjOzHtASgZaU5m1J4IWuVI6Ib1f+LWlR0prDE4AtgFsLVZcGtq8kfZc0FfiLpA0i4rEqh74vvz9THo7txjk7UpnQ9WKxMKete390ISJm5/eJFNZklrQY6V7s7ZIGRkRHf68VSNf+Zt53NeAMSUtHxPQuttfMzOZTKw0dQ0pd1ylJX5J0p6S3gdmkgAewXqnq6EqQza4GBGzc3YZ145xdUb7Oa4FZlVdhKBtJ35D0gKQpefvtXTzvfZUgm1UyFa1RrihpiKR2Se0TJ07sznWYmVknWiXQTgJmAAM7qyhpY1JgmgB8gzTjdrO8uZw/9rXih4iYRkq5t3p3GtfNc3bkP/m9fK/0EFLw/37pvDuT0uDdRRru3oy5M687O+9bpc+1cuwSEcMioi0i2vr379/JYc3MrDtaYug4ImZJuoN03/EXnVTfmXS/cs9KAnZJa9eoO6D4IU9GWo659127qjvn7Mg4Ujq7bUj3ggGIiKfzMZcr1d8duCciDqgUSNqyjvOamVmTtEqPFuBMoE3St8obJC0iabv8cWlgViXgZfvUOOYXS8FrZ9Kwba08cLV6fd05Z025Rz0MOFDSBl3YZWlST7+o2+c1M7PmaYkeLUBEjJB0OnC+pM2Ba0jDvOuThlTHAyOB0cAhks4ERpBmK3+9xmGnA9dLOoU0XHwKcHVEPFqtckTMlPQcsIekf5GSqz/czXN25ihgE+AuSb8HbsvnWQP4FvBebjf5vGdLOgq4hzRjeOs6z2tmZk3QMoEWICIOl3QncBBwCalHN550f/TUXOcGST8FDiY9tnMXsAPwZJVDXkZ6NvZ80pDxtRRm8dbw/Xyuf5BmQq/TzXN2do3TJH0eOBDYm3R/djHSTOQxwH9HRGWlqnOBdYEfkXrZo/M+XqDCzKyP0LyjoQsOSeOB4RHx42a3pS9pa2uL9vZaI+tmZlaNpHER0VZtWyvdozUzM1vgONCamZk1UEvdo+1JETGo2W0wMzNzj9bMzKyBHGjNzMwayIG2l0gaWki9Fzl133WSNmp228zMrHEcaHvX26R1kgeTnp9dDxgt6UNNbZWZmTXMAjsZqkXNLqTfuzs/63sXsB1pgQ4zM1vAuEfbXA/l97UAJA2WdK2klyVNlfSgpA+sbSxpbUmXSnpd0jRJD0vau7B9KUknS3pR0gxJD0lywnczsyZwj7a5KmkBn8vvawN3AOeQ1j/eHLhQ0pyIuBRA0gBSL3ga8GPS0o0bkoN1Npy0nvLRwDPAHsC1ktoi4sGGXpGZmc3DgbaXSar8zdcGfg88SEqgQERcVqgn4FZS7tr9gUvzpkOBFYFPR0Ql3d+Ywn5bA9sDW0XELbl4lKT1SAkNdm/AZZmZWQ0OtL1rFWBW4fMkYOOImAEgaWXgGGAnUjafRXO9lwr7fB4YWQiyZV8AXgHuKAR1SMF432o7SBoCDAEYOHBgtSpmZlYn36PtXW8DGwObAd8DlgAukVT5Hi4C9iSl89sm172AefPjrkLHiev7AauRAnrxNZR5h5ffFxHDIqItItr69+9fz3WZmVkN7tH2rtkRUUmNc4+k6cCfgd0lXUNKvXdgRJxT2aEQhCsmkXLr1vIGqQf81Z5rtpmZ1cs92ub6K/Bv4Kek3LeLADMqGyUtD3yltM8YYFtJq9Y45hhSj3ZKRLSXXz1+BWZm1iH3aJsoIkLSCcDFQBtwH/ArSZOBOcDPSMPNKxR2OwP4JnCbpF+TZh1vACwbESeTksPfRFoI4yRSIF8B+CSwVEQc2SsXZ2ZmgHu0reBy4CngCGBv4FnScPJvgb/lf78vIiaSHvt5ADgTuI40kemFvD2AXUj3dg8hBd1zSatR3d7wqzEzs3ko/XfZLGlra4v2do8wm5l1h6RxEdFWbZt7tGZmZg3kQGtmZtZADrRmZmYN5EBrZmbWQA60ZmZmDeRAa2Zm1kALTKCV9B1JIWnNUvlJufzrpfIv5vLPzMc5B+Vj7NDMY5iZWetaYAItcGd+LwfOz5Byt1YrnwGMa3C7zMxsIbYgBdrHSQvqvx9QJS1OWtrwz1QPtOMqKeq6S9JSndcyM7OF3QITaPPSg3cxb0D9VH7/A7BhXqS/khFnU+CO/HkPSY9ImiHpRUm/LuZylbRvHt7dRNLYnHXnJ9XaIelzkt7JaxhXytaWdKmk1yVNk/SwpL1rXYukb0q6XdIbkt6U9E9JbaU6n5A0MteZKukxSQcWtm8h6TZJk/PrQUlO+m5m1ssWtKQCdwJDJS0dEdNJ6/uOA/5FWpx/U+AfwCeAFUnJ0bchrTf8Z1Lw3Ag4jpT39ful419KCtrHAG+VTy5pW+Bq4DcRcVwuG0D6ATAN+DEpCcCG1MgNmw3K7XmGlLN2L1ISgU9ExLO5zgjgMeDrpCHwj5GTD0hagbQG8jXAsYCA/wJW6uCcZmbWAAtaoL0DWJyUMP1WUu/2rpwl5+78+R/M7fXeSQpYYyPiW7lspCSA30g6PiImFI5/VkT8tvJB0qDCv78CXAH8IiJOLexzKCmofzoiKgnbx3R0ERFxbOG4i5Ay8mxCCqrHSuoHrAPsFBGPVDnmevmcB0XEO7lsVK3zSRpCSkzAwIEDO2qamZl10wIzdJzdB8xmbiD9DKk3CXB3qfwp0j3d/wGuLB3nctLfZnCp/Poa5901H+PwUpAF+DwwshBkOyVpA0lXS3oVeA+YReqxrpervEHqGZ8jac/cay56BpgCXCJpJ0kd9mQjYlhEtEVEW//+/bvaTDMz64IFKtBGxDTgQeAz+TGfNZk7G/kuYDOl7upnSL3ffqQe8KulQ1U+f6hGedlXSMHv6irbVgG6E2SXJ/U+1wIOAz5L6qE/BCwFEBFzgG2AV0jp8F7J92M/lbe/CXwxX9sVwERJ10tat6vtMDOznrFABdrsDlJP9DPA+Ih4JZffCywPbAV8hBSAXyf1Fss9wlXz+xul8lo5BQ8GXgJGSVqltG0SsHo32j+Y9APh6xFxcUTcHhHtpKHguQ2JeDwidiXdd/0CKQhfn4eaiYi7I2K7vH0XUm/4km60w8zMesCCGGjvJPVUv8XcYWMiYjLwb9KEJIA7IuI90mSp8mzcPYA5xf07MRnYNv/7pjwZqWIMsK2kVT+4W1VL5/f3HzvKi2oMqlY5ImZFxM3A6aSAvlJp+/SIGEHq+X68i20wM7MesqBNhoK5Q8VfAn5U2nYXsD/wJmnGLsDRpOB4IXAZaXbuccB5pYlQHYqISZK+CNwGXCdpuzyUfQbwTdKs4V+T7q1uACwbESdXOdTdpPur50k6mdS7HUrqMQMgaSPgVNK95GeBlYGfAg9FxBuStge+DfwdeAFYA/gecHNXr8fMzHrGAtejzcHxBdIjLeUe6V2V8vzcLRExCvgaaWGLEcAhwGnAQXWc+2Vga1Lv8ypJS0TERGBz4AHgTNJjN0NyG6sd41VSD3s10uM5h5AeM3q6UO0V0v3io4AbSY8cPUa6V0yuG8AJpPu9JwMjScHXzMx6kXK8MQOgra0t2tvbm90MM7M+RdK4iGirtm2B69GamZm1EvdobR6SJgLPV9nUjzRLe2Hga11wLUzX62vtXWtHRNWFCBxorUsktdcaFlnQ+FoXXAvT9fpaW4eHjs3MzBrIgdbMzKyBHGitq4Y1uwG9yNe64FqYrtfX2iJ8j9bMzKyB3KM1MzNrIAdaMzOzBnKgtaok7S/pKUnvShonaesu7DNUUlR5bdcbbe6MpI9LGiNpmqT/SDpW0qJd2G9FSRdKelPS25IurpKlqeXUc72SBtX4Di/rrXbXQ9JHJJ0r6WFJ70ka28X9+tx3W8+19sXvVdLukq6V9JKkKfm/Q3t1Yb8lJZ0m6TVJU3OK0EGNb3FtC2JSAZtP+X/M55CSGdwO7EdKlLBxRPyrk93fBsqB9bFqFXuTpJWBfwCPAjsB/4+0pvUiwC862f0KUprB75KyOp1EStjw2Ua1d37N5/VCynJ1R+FzsxcD6MwngC+TknIs3o39+tx3S/3XCn3rez0MeA44lNTOLwOXSOoXEb/rYL+zgN3yfhNJ/x0bLem/IuLdxja5hojwy695XsATwAWFz4sAjwB/7WS/ocDrzW5/jbYdScratEKh7AhgWrGsyn6DSQka/rdQtkku+0Kzr6sB1zsoX9sOzb6Gbl7vIoV/DwfGdmGfvvrd1nOtfe57BfpVKbsEeK6DfdYEZgPfLJStAcwEvtusa/HQsc1D0rqkX/hXVMoiYg5wJSn1YF/1JeCmSHmJKy4j5f/dspP9Xo2IWysFEXEv6Zd2K/896r3ePin/b7S7+uR3W+e19jkRUa23/QDw4Q522ya/X1U4zkukkbmmfacOtFa2fn5/vFT+GPAhSVXX8ixYSdLrkmZJekDSLj3fxLqsT+maIuIFUg9v/ap71Ngve6yT/Zqt3uutuDDf/3tZ0umSlm5EI5usr36386Ovf6+DgSc72L4+MCEippTKm/qd+h6tla2c398qlb9Z2D6xxr5Pk4YnHwCWJyWb/5ukXSPiqhr79JaV+eA1QbqulauUd2W/dXugXY1S7/XOAM4m5TGeDGwF/JR0j3ennm1i0/XV77Yeff57zRMyv0rHebXr/d99QznQLgQkrQis3lm9iKj2677LIuKvpfOOAO4EfkVhKMdaV0S8DBxUKBor6VXgD5L+OyIealLTbD709e81zxq+BLgmIi5qamPq4KHjhcPupKGTzl4wt+e6YukYK5e2dyrSTISrgI268hhNg73JB68J0nV1dE317tdsPdnu4fn90/PVotbTV7/bntInvldJHwJuJKXv3KeT6i35nTrQLgQi4o8Roc5euXqlV1u+n7E+8EZE1Bo2rnn6/Gq2xyldk6S1gGWofp+u5n5Zrft7raLe660mSu8Lir763faUlv9eJS0DXAcsQZoxPa2TXR4H1pK0bKm8qd+pA63NIyKeJU022L1SJmmR/PnG7hxLkoBdgYci4r2ebGcdbgS2lbR8oWxPYDpwSyf7rSZpi0qBpDbSPbxu/T16Wb3XW81u+X1cTzSshfTV77antPT3Kmkx0tMOHwW2i4jXurDbqPy+c+E4HyY9F92877TZz0r51XovYC/gPdLCBp8DLiL9B3rDQp0tSc+rbVkouwX4IWmK/c7ADaRFAL7SAte0MvAyMBr4AjAEmAIcX6r3NHB+qewm4FlgF9JkjCeA25p9TY24XtKz0Kfla/0CcGz+7v/W7Gvq5HqXIQWO3YC7gH8XPi+zgH233b7Wvvi9kjLyRP5vymal15K5zhhgTGm/c0kLXHyDtHjO3cBTwFJNu5Zm/zH9as0XsH/+P+sM4H5g69L2rfL/CbYqlJ2f/6M1HZgK3AZ8qdnXUmjfx4Gbc/teBo4DFi3VGQ9cVCpbCbiQNJtxMmlSxgcepm+1Vz3XC3wNaCet8DUz/2/g2Mp/2Fr1xdwFGaq9Bi1I320919oXv9d8DZ1d51hKC3YASwKnk56OmEr6wb9OM6/FafLMzMwayPdozczMGsiB1szMrIEcaM3MzBrIgdbMzKyBHGjNzMwayIHWzMysgRxozczMGsiB1szMrIH+PwatoWw11Nx0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cf_ind_adult.draw(sample_ind)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ASV\n"
     ]
    },
    {
     "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=\"1534pt\" height=\"131pt\"\n",
       " viewBox=\"0.00 0.00 1534.44 131.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 127)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-127 1530.44,-127 1530.44,4 -4,4\"/>\n",
       "<!-- Education&#45;Num noise -->\n",
       "<g id=\"node1\" class=\"node\"><title>Education&#45;Num noise</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"88.3922\" cy=\"-105\" rx=\"88.2844\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"88.3922\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">Education&#45;Num noise</text>\n",
       "</g>\n",
       "<!-- predicted income &gt;50k -->\n",
       "<g id=\"node2\" class=\"node\"><title>predicted income &gt;50k</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"148.392\" cy=\"-18\" rx=\"92.8835\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"148.392\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">predicted income &gt;50k</text>\n",
       "</g>\n",
       "<!-- Education&#45;Num noise&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge1\" class=\"edge\"><title>Education&#45;Num noise&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"#ff0051\" stroke-width=\"5\" d=\"M100.534,-86.799C109.129,-74.6221 120.8,-58.0889 130.464,-44.3977\"/>\n",
       "<polygon fill=\"#ff0051\" stroke=\"#ff0051\" stroke-width=\"5\" points=\"134.28,-46.5781 136.473,-35.8854 127.132,-41.532 134.28,-46.5781\"/>\n",
       "<text text-anchor=\"middle\" x=\"134.892\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#ff0051\">2.08</text>\n",
       "</g>\n",
       "<!-- Hours per week noise -->\n",
       "<g id=\"node3\" class=\"node\"><title>Hours per week noise</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"282.392\" cy=\"-105\" rx=\"87.9851\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"282.392\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">Hours per week noise</text>\n",
       "</g>\n",
       "<!-- Hours per week noise&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge2\" class=\"edge\"><title>Hours per week noise&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"#ff0051\" stroke-width=\"1.26733\" d=\"M256.556,-87.6111C235.515,-74.2646 205.632,-55.309 182.606,-40.7028\"/>\n",
       "<polygon fill=\"#ff0051\" stroke=\"#ff0051\" stroke-width=\"1.26733\" points=\"184.469,-37.7396 174.149,-35.3385 180.719,-43.6506 184.469,-37.7396\"/>\n",
       "<text text-anchor=\"middle\" x=\"236.892\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#ff0051\">0.51</text>\n",
       "</g>\n",
       "<!-- Sex -->\n",
       "<g id=\"node4\" class=\"node\"><title>Sex</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"415.392\" cy=\"-105\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"415.392\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">Sex</text>\n",
       "</g>\n",
       "<!-- Sex&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge3\" class=\"edge\"><title>Sex&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"#ff0051\" stroke-width=\"4.0258\" d=\"M394.308,-93.377C389.469,-91.1263 384.301,-88.8688 379.392,-87 324.313,-66.0331 259.847,-47.5376 213.033,-35.1497\"/>\n",
       "<polygon fill=\"#ff0051\" stroke=\"#ff0051\" stroke-width=\"4.0258\" points=\"213.888,-31.7323 203.322,-32.6002 212.1,-38.5465 213.888,-31.7323\"/>\n",
       "<text text-anchor=\"middle\" x=\"335.892\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#ff0051\">1.67</text>\n",
       "</g>\n",
       "<!-- Relationship noise -->\n",
       "<g id=\"node5\" class=\"node\"><title>Relationship noise</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"536.392\" cy=\"-105\" rx=\"76.0865\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"536.392\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">Relationship noise</text>\n",
       "</g>\n",
       "<!-- Relationship noise&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge4\" class=\"edge\"><title>Relationship noise&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"#ff0051\" stroke-width=\"3.49475\" d=\"M491.321,-90.4596C454.252,-79.6365 400.229,-64.5779 352.392,-54 312.463,-45.1707 267.721,-37.2749 230.336,-31.2274\"/>\n",
       "<polygon fill=\"#ff0051\" stroke=\"#ff0051\" stroke-width=\"3.49475\" points=\"230.424,-27.697 219.996,-29.5715 229.317,-34.6089 230.424,-27.697\"/>\n",
       "<text text-anchor=\"middle\" x=\"426.892\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#ff0051\">1.45</text>\n",
       "</g>\n",
       "<!-- Occupation noise -->\n",
       "<g id=\"node6\" class=\"node\"><title>Occupation noise</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"702.392\" cy=\"-105\" rx=\"72.2875\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"702.392\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">Occupation noise</text>\n",
       "</g>\n",
       "<!-- Occupation noise&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge5\" class=\"edge\"><title>Occupation noise&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"#ff0051\" stroke-width=\"0.39214\" d=\"M649.272,-92.6646C639.987,-90.7366 630.414,-88.7827 621.392,-87 542.459,-71.4024 522.939,-66.0795 443.392,-54 374.946,-43.6062 297.171,-34.4894 238.669,-28.168\"/>\n",
       "<polygon fill=\"#ff0051\" stroke=\"#ff0051\" stroke-width=\"0.39214\" points=\"239.031,-24.6869 228.715,-27.0993 238.284,-31.6469 239.031,-24.6869\"/>\n",
       "<text text-anchor=\"middle\" x=\"542.892\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#ff0051\">0.14</text>\n",
       "</g>\n",
       "<!-- Age -->\n",
       "<g id=\"node7\" class=\"node\"><title>Age</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"819.392\" cy=\"-105\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"819.392\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">Age</text>\n",
       "</g>\n",
       "<!-- Age&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge6\" class=\"edge\"><title>Age&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"#ff0051\" stroke-width=\"1.93522\" d=\"M798.866,-93.0367C793.91,-90.7488 788.552,-88.5609 783.392,-87 735.232,-72.4328 609.371,-59.804 559.392,-54 452.012,-41.53 328.362,-31.5966 244.966,-25.5442\"/>\n",
       "<polygon fill=\"#ff0051\" stroke=\"#ff0051\" stroke-width=\"1.93522\" points=\"244.971,-22.0356 234.746,-24.8072 244.468,-29.0174 244.971,-22.0356\"/>\n",
       "<text text-anchor=\"middle\" x=\"695.892\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#ff0051\">0.79</text>\n",
       "</g>\n",
       "<!-- Marital Status noise -->\n",
       "<g id=\"node8\" class=\"node\"><title>Marital Status noise</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"946.392\" cy=\"-105\" rx=\"81.7856\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"946.392\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">Marital Status noise</text>\n",
       "</g>\n",
       "<!-- Marital Status noise&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge7\" class=\"edge\"><title>Marital Status noise&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"#008bfb\" stroke-width=\"1.38788\" d=\"M895.195,-90.7639C848.045,-79.1242 775.98,-62.7174 712.392,-54 552,-32.0114 364.035,-23.8277 250.691,-20.788\"/>\n",
       "<polygon fill=\"#008bfb\" stroke=\"#008bfb\" stroke-width=\"1.38788\" points=\"250.722,-17.2878 240.635,-20.5266 250.541,-24.2854 250.722,-17.2878\"/>\n",
       "<text text-anchor=\"middle\" x=\"811.892\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#008bfb\">&#45;0.56</text>\n",
       "</g>\n",
       "<!-- Workclass noise -->\n",
       "<g id=\"node9\" class=\"node\"><title>Workclass noise</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"1115.39\" cy=\"-105\" rx=\"69.5877\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"1115.39\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">Workclass noise</text>\n",
       "</g>\n",
       "<!-- Workclass noise&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge8\" class=\"edge\"><title>Workclass noise&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"#ff0051\" stroke-width=\"0.073609\" d=\"M1064.71,-92.5185C1055.63,-90.5864 1046.24,-88.669 1037.39,-87 945.843,-69.7404 923.044,-63.7346 830.392,-54 625.411,-32.4634 384.32,-24.0367 250.691,-20.8461\"/>\n",
       "<polygon fill=\"#ff0051\" stroke=\"#ff0051\" stroke-width=\"0.073609\" points=\"250.548,-17.3419 240.469,-20.6073 250.384,-24.34 250.548,-17.3419\"/>\n",
       "<text text-anchor=\"middle\" x=\"953.892\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#ff0051\">0.01</text>\n",
       "</g>\n",
       "<!-- Capital Gain noise -->\n",
       "<g id=\"node10\" class=\"node\"><title>Capital Gain noise</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"1279.39\" cy=\"-105\" rx=\"76.8869\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"1279.39\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">Capital Gain noise</text>\n",
       "</g>\n",
       "<!-- Capital Gain noise&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge9\" class=\"edge\"><title>Capital Gain noise&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"#008bfb\" stroke-width=\"0.893073\" d=\"M1223.69,-92.4353C1213.94,-90.5328 1203.88,-88.6469 1194.39,-87 1095.24,-69.7935 1070.57,-63.4865 970.392,-54 712.276,-29.5582 406.938,-22.181 251.107,-19.9573\"/>\n",
       "<polygon fill=\"#008bfb\" stroke=\"#008bfb\" stroke-width=\"0.893073\" points=\"251.109,-16.457 241.062,-19.8181 251.012,-23.4564 251.109,-16.457\"/>\n",
       "<text text-anchor=\"middle\" x=\"1104.89\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#008bfb\">&#45;0.35</text>\n",
       "</g>\n",
       "<!-- Capital Loss noise -->\n",
       "<g id=\"node11\" class=\"node\"><title>Capital Loss noise</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"1450.39\" cy=\"-105\" rx=\"76.0865\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"1450.39\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">Capital Loss noise</text>\n",
       "</g>\n",
       "<!-- Capital Loss noise&#45;&gt;predicted income &gt;50k -->\n",
       "<g id=\"edge10\" class=\"edge\"><title>Capital Loss noise&#45;&gt;predicted income &gt;50k</title>\n",
       "<path fill=\"none\" stroke=\"#008bfb\" stroke-width=\"0.128419\" d=\"M1395.16,-92.4381C1385.27,-90.5145 1375.04,-88.6208 1365.39,-87 1258.34,-69.0138 1231.54,-63.292 1123.39,-54 806.374,-26.7633 429.61,-20.653 251.486,-19.3299\"/>\n",
       "<polygon fill=\"#008bfb\" stroke=\"#008bfb\" stroke-width=\"0.128419\" points=\"251.363,-15.8291 241.339,-19.2582 251.314,-22.8289 251.363,-15.8291\"/>\n",
       "<text text-anchor=\"middle\" x=\"1267.89\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#008bfb\">&#45;0.03</text>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.files.Source at 0x7f72a948b090>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# asv\n",
    "print('ASV')\n",
    "cf_flow_adult.fold_noise = False\n",
    "g = cf_flow_adult.draw_asv(sample_ind, max_display=10, show_fg_val=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate difference with the ground truth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "gt = get_direct_effect_ground_truth(cf_flow_adult.graph)\n",
    "theta_ind = get_effect_ind(cf_ind_adult)\n",
    "theta_manifold = get_effect_manifold(cf_on_manifold_adult)\n",
    "theta_asv = get_effect_asv(cf_flow_adult)\n",
    "theta_flow = get_direct_effect_flow(cf_flow_adult)\n",
    "\n",
    "gt_indirect = get_indirect_effect_ground_truth(cf_flow_adult.graph)\n",
    "theta_flow_indirect = get_indirect_effect_flow(cf_flow_adult)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXxU5fn//9fFJpsRBQIRRCwu4FapQbQIdcMNBHdAK+KGG7VqrVqrHxGlVdHSWpcWXNDfTwUKUhRqhYpUVFQWVwStVlQom4ICBlDg+v5xn8TJMFkmmWROkvfz8cgjc+6zzHXOLNec+9znvs3dERERiZt62Q5AREQkFSUoERGJJSUoERGJJSUoERGJJSUoERGJpQbZDiDTWrVq5R07dsx2GCIiUk4LFiz40t1bJ5fXugTVsWNH5s+fn+0wRESknMzss1TlquITEZFYUoISEZFYUoISEZFYUoISEZFYUoISEZFYUoISEZFYqtZm5mb2KNAXWO3uB0ZluwETgI7AUuBsd19nZgb8CTgZKACGuPvC6oxXRKrP+vXrWb16Nd9//322Q5EMadiwIbm5ueTk5FRo/eq+D2occD/wRELZjcCL7n6nmd0YTd8AnATsE/11Bx6K/otILbN+/XpWrVpFu3btaNKkCeH3qdRk7s6mTZtYvnw5QIWSVLVW8bn7y8DapOL+wOPR48eBUxPKn/DgdaCFmeVVT6QiUp1Wr15Nu3btaNq0qZJTLWFmNG3alHbt2rF69eoKbSMOPUm0cfcV0eOVQJvocTvgi4TllkVlK0hiZkOBoQAdOnSoukhFqtBz19xW6vxTRt9aTZFUv++//54mTZpkOwypAk2aNKlwtW2sGkl4GN437SF+3X2Mu+e7e37r1jt05yQiNYDOnGqnyryucUhQqwqr7qL/heeCy4E9EpZrH5WJiEgdEIcE9SxwfvT4fGBqQvlgCw4HvkmoChQRkVquupuZPw0cBbQys2XArcCdwEQzuwj4DDg7WvwfhCbmHxOamV9QnbGKSPaNnvlRVp73mt77pr3OkCFDeP/99ys9msK0adM45ZRT+PTTT6nrQwdVa4Jy90ElzDo2xbIOXFm1EYmIZMYtt9zCpk2bsh1GrRKHVnwiIjVep06dsh1CrROHa1AiIjXekCFDyM/PB2DcuHGYGe+99x69e/emWbNmdO7cmWeeeabYOu7O8OHDyc3NZeedd2bw4MGsX78+G+HHkhKUiEgVOeecc+jXrx9Tpkxhn332YeDAgSxbtqxo/n333ceIESMYOnQokyZNokmTJlx//fVZjDheVMUnIlJFrrnmGi688EIADj30UNq0acO0adO47LLL2LZtG3fddReXXnopd9xxBwAnnHACvXv3LuoeqK7TGZSISBU5/vjjix63bNmS3NzcojOoL774ghUrVtC/f/9i65x++unVGmOcKUGJiFSRFi1aFJtu1KgRmzdvBmDlypUA5ObmFlsmebouU4ISEcmCtm3bAuzQkWpFO1atjZSgRESyYI899qBt27ZMnTq1WHlyS7+6TI0kRESyoH79+lx//fVcd911tGrVip49ezJ58mQWL16c7dBiQwlKRGKrIl0O1SRXX301a9eu5S9/+Qt//OMf6devH3fffTfnnntutkOLBQs9CtUe+fn5Xtm+sESyoS6PB7V48WK6dOmS7TCkipT1+prZAnfPTy7XNSgREYklJSgREYklJSgREYklJSgREYklJSgREYklJSgREYklJSgREYklJSgREYklJSgREYklJSgRkVpg6dKlmBnTpk0rKvv2228ZOHAgLVu2xMwYN25cubY1e/ZszIz333+/1OWuu+46OnbsWImoS6e++EQkvl76fXae9+jfZOd5KyEvL4+5c+fSuXPnorKHHnqI5557jieeeIJ27drRqVOnLEaYPiUoEZFaYKedduLwww8vVrZkyRL2228/zjjjjCxFVTmq4hMRyYAhQ4aQn5/P9OnT2X///WnatCl9+vRh7dq1fPzxxxx99NE0a9aM/Px83n333aL17r33Xrp168Yuu+xCmzZtOOWUU/j444+Lbfuoo47izDPP5KmnnmLvvfcmJyeHk046qWj4eNixiq9jx4488sgjvPXWW5gZZla07KxZs+jevTuNGzemTZs2XHHFFWzcuLHU/fv6668555xzaN68OXl5eYwcOTITh61USlAiIhny+eef83//93/ccccdjBkzhtdee42hQ4cycOBABg4cyKRJk9i6dSsDBw6kcCSJZcuWMWzYMKZOncrYsWPZtm0bP/3pT/nmm2+KbfuNN97g/vvv595772XMmDEsXLiQoUOHlhjLlClTOPnkk+ncuTNz585l7ty5ACxatIgTTzyRVq1aMXnyZG677TaeeuopzjzzzFL37YILLuD5559n9OjRjBkzhhkzZjB+/PhKHrHSqYpPRCRD1q5dy9y5c4uu9bz77ruMGjWKxx9/nMGDBwPg7vTp04clS5bQpUsXRo8eXbT+tm3b6N27N7m5uUydOrVoHYD169czffp0dt11VwBWrlzJNddcw6ZNm2jSpMkOsXTt2pXWrVuzatWqYlV/t99+O3vuuSfPPvss9evXB2C33XZjwIABzJ07lyOOOGKHbS1atIi///3vjB8/ngEDBgBw9NFH06FDB3Jycip72EqkMygRkQzp2LFjsYYIe++9NwDHHHPMDmXLly8H4PXXX6d37960bNmSBg0a0LRpUzZu3MhHH31UbNvdunUrSk4A+++/f7HtlNebb77JaaedVpScAM444wwaNGjAK6+8knKdefPmAdC/f/+isubNm9O7d++0njtdSlAiIhnSokWLYtONGjXaobywbPPmzXz++eccf/zxuDt//etfefXVV5k3bx65ubls3ry5XNtOXq4sK1asoE2bNsXK6tevT8uWLVm7dm3KdVauXMnOO+9M48aNi5Xn5uam9dzpUhWfiEiW/POf/6SgoICpU6fSrFkzALZu3VpiosiEvLw8Vq9eXaxs27ZtfPXVV+y2224p12nbti0bNmxg8+bNxZJU8nYyTWdQIiJZsmnTJurVq0eDBj+cK0ycOJGtW7dW2XN2796dKVOmsG3btqKyZ555hq1bt3LkkUemXKdbt24ATJ06tahs48aNzJw5s8riBJ1BiYhkzTHHHMO2bdu44IILuOiii1i0aBH33HPPDtV5mXTzzTfTtWtXTj31VC6//HKWLVvGDTfcwAknnJCygQTAAQccQL9+/bj88stZv349eXl5jBo1iqZNm1ZZnKAEJSJxVgN7dEjHQQcdxLhx4xg+fDhTpkzhxz/+MX/729+KWspVhQMOOIDnn3+em266idNPP52cnBwGDRrE3XffXep648aN4/LLL+fqq6+mefPmXHnllXTr1o1JkyZVWaxW2Ba/tsjPz/f58+dnOwyRtD13zW2lzj9l9K3VFEn1W7x4MV26dMl2GFJFynp9zWyBu+cnl+salIiIxJISlIiIxFJsEpSZXWNmi8zsfTN72swam9leZvaGmX1sZhPMrFG24xQRkeoRiwRlZu2Aq4B8dz8QqA8MBO4CRrv73sA64KLsRSkiItUpFgkq0gBoYmYNgKbACuAYoLCJyOPAqVmKTUREqlksEpS7LwfuAT4nJKZvgAXA1+5eeMfaMqBdqvXNbKiZzTez+WvWrKmOkEVEpIrFIkGZ2a5Af2AvYHegGXBiedd39zHunu/u+a1bt66iKEVEpDrFIkEBxwGfuvsad/8eeAboAbSIqvwA2gPpddsrIiI1VlwS1OfA4WbW1MKwj8cCHwAvAYWjaJ0PTC1hfRERqWVikaDc/Q1CY4iFwHuEuMYANwDXmtnHQEvgkawFKSIi1So2ffG5+61Acl8u/wUOy0I4IhIDD779YFae94pDrsjK80pxsTiDEhERSaYEJSIisaQEJSKSIXPnzqVfv37k5eXRrFkzDjnkEJ588smi+V9//TUXX3wxu+++O40bN6ZDhw5ccsklAMyePRszY9GiRcW2uW7dOho1asTDDz9crfsSB7G5BiUiUtN99tln9OjRg8suu4zGjRvz6quvcsEFF1CvXj0GDRrEtddey2uvvcbo0aNp27YtX3zxBS+//DIAvXr1Ii8vj4kTJ3LbbT8MvTJlyhQAzjjjjKzsUzYpQYmIZMjAgQOLHrs7vXr1YtmyZYwdO5ZBgwbx5ptvcuWVVxYbkPDnP/85APXq1eOss85iwoQJxRLUhAkTOP7449l1112rb0diQlV8IiIZsm7dOq666ir23HNPGjZsSMOGDRkzZgwfffQRAIcccgijRo3iwQcfLCpLNGDAAD788EPeeecdAL788ktmzZpVpSPsxpkSlIhIhgwZMoQJEybw61//mhkzZjBv3jwuvPBCNm/eDMD999/PqaeeyogRI9hvv/3YZ599GD9+fNH6RxxxBB06dGDChAkATJ48mQYNGnDqqXWzn2wlKBGRDNi8eTPTpk3jtttuY9iwYRxzzDHk5+ezffv2omVatGjBfffdx8qVK3nnnXfo3r075557Lh988AEAZsbZZ5/NxIkTgVC9d9JJJ7HzzjtnZZ+yTQlKRCQDtmzZwvbt29lpp52KyjZs2MCzzz6bcvmDDz6YUaNGsX37dpYsWVJUPnDgQD755BOmTZvGv//972LXteoaNZIQEcmAXXbZhW7dujFixAhycnKoV68ed955J7vssgvr168H4Mgjj+S0007jwAMPxMwYO3YszZo147DDfugw59BDD2Xvvfdm6NChNGnShL59+2Zrl7JOCUpEYqumdTn01FNPcemllzJ48GBatmzJsGHDKCgo4P777wfCNaZx48axdOlS6tevT9euXXn++edp3759se0MGDCAkSNHMnDgQJo2bZqNXYkFc/dsx5BR+fn5Pn/+/GyHIZK25665rdT5p4xO7qqy9li8eDFdunTJdhhSRcp6fc1sgbvnJ5frGpSIiMSSEpSIiMSSEpSIiMSSEpSIiMSSEpSIiMRS2gnKzA42swlm9omZbTGzn0TlI83spMyHKCIidVFaCSpKQAuAtsATQMOE2VuAX2QuNBERqcvSPYP6PTDO3X8GjEya9zZwSEaiEhGROi/dBNUZmBA9Tr7Ddz2wW6UjEhERIf0EtRr4UQnzDgA+r1w4IiI114gRI2jXrh316tWjY8eOmBnvv/9+tsOqsdLti288MMLMPgDmRmVuZvsCNwCPZDI4Eanb1vz5/qw8b+tfDEt7nfnz53Prrbfyu9/9jqOOOooPPviAiy++uAqiqzvSTVC3APsD/wZWRmVTCY0mZgC/y1xoIiI1R+GQGVdeeSU5OTls2bIlyxHVfGlV8bn7FnfvCxwPPA48DDwF9HH3vu7+fRXEKCISa0OGDOG8884DwrAbZpZyuYKCAq666iratm1L48aN6datGzNmzCia/9hjj9GsWTO+//6Hr9Ldd9+dli1bUtix9/bt22nRogVjx46twj2KhwrdqOvuL7r7Te4+1N1vdPeZmQ5MRKSmuOWWW7j55psBmDVrFnPnzi0aAyrRJZdcwmOPPcZvf/tbpkyZwh577EGfPn145ZVXAOjZsycFBQUsXLgQgP/85z+sXr2aDRs2FI26+8477/DNN9/Qs2fPatq77Emris/MOpQyezuw3t13fFVERGqxTp060alTJwC6detG8+bNmT17drFlFi9ezNNPP81jjz3G+eefD8AJJ5zAwQcfzO23384LL7zA3nvvTV5eHnPmzKF79+7MmTOHH//4xzRq1Ig5c+ZwwAEHMGfOHFq3bk3nzp2rezerXbpnUEuBT0v4+wxYZ2afmtk1mQxSRKSmmzdvHu7OWWedVVRWr149zjrrrKIzKAhnUXPmzAHg5ZdfplevXvTq1atY2ZFHHlm9wWdJugnqHGAZ8E/gCuCs6P8LwHLgSuBF4G4lKRGRH6xYsYLmzZvvMEJumzZtKCgoKGpU0bNnT1555RXcnTlz5tCzZ89iSeuVV16pE9V7kH4rvuOAZ909uUujv5rZn4GfuvtgM9sIXAaMzkSQIiI1XV5eHhs3bqSgoKBYklq1ahVNmzZlp512AkKCWrt2LTNnzuTTTz+lZ8+eNGjQgOXLlzNjxgxWrVqlBFWCs4AzSpj3LDApevw8IUGJSDm1/2ZBtkOQKtStWzfMjEmTJjF48GAA3J1JkyYVq7I76KCDaNGiBSNHjqRz5860bt0agAMPPJCRI0fSvHlzunbtmpV9qG7pJqjNQA/gXynm9YjmAxjwbSXiEhGpVbp06cKgQYMYNmwYGzZsoFOnTowdO5YlS5bw0EMPFS1Xr149evTowfTp07n00kuLynv27MkDDzxA7969qV+/fjZ2odqlm6DGALeYWUvgOWAN0BroTzhjKrxR96fAO5kKUkTqpor06BBnY8eO5YYbbmDEiBF8/fXXHHTQQUybNm2HRg89e/Zk+vTp9OrVq1jZAw88UGcaSABY4c1f5V4hNH74NaH3CCecLa0ERrn76GiZA4Bv3X1pRqMth/z8fJ8/f351P61Ipb11Yb9S53d99NlqiqT6LV68mC5dumQ7DKkiZb2+ZrbA3fOTy9M9g8LdR5vZn4A9CElqJfCFu29PWGZRutsVERFJVNGeJLa7+2fu/kb0f3vZa5XOzFqY2SQzW2Jmi83sCDPbzcxmmtl/ov+7VvZ5RESkZkj7DMrMdiZcc9oXaJw8392vr2AsfwL+6e5nmlkjoClwE/Ciu99pZjcCNxJ6TRcRkVou3a6OOgGvAU2AZoRGErtF21kHfAOknaDMbBegFzAEwN2/A74zs/7AUdFijwOzUYISEakT0q3iGw3MA9oQGkecTEhWPwc2AgMqGMdehGT3mJm9ZWYPm1kzoI27r4iWWRk97w7MbKiZzTez+WvWrKlgCCKSTek22JKaoTKva7oJ6jDgL0DhQCeN3H2buz8F3EuopquIBsBPgIfcvSvhHqobExfwsJcp99Tdx7h7vrvnF97UJiI1R8OGDdm0aVO2w5AqsGnTJho2bFihddNNUI0JPZZvB9YCuyfMex/4cYWiCP37LXP3N6LpSYSEtcrM8gCi/6sruH0RibHc3FyWL19OQUGBzqRqCXenoKCA5cuXk5ubW6FtpNtI4iNgz+jxW8BlZvYPYBtwEfC/igTh7ivN7Asz28/dPwSOBT6I/s4H7oz+T63I9kUk3nJycgD43//+V2ywPqnZGjZsSJs2bYpe33Slm6DGA4cA/x9h+PcXgPWEsaAaEJJIRf0CeDJqwfdf4ALCGd5EM7uIMJzH2ZXYvojEWE5OToW/yKR2SitBufsfEh6/bmYHAicSGkrMcvf3KxqIu78N7HAnMeFsSkRE6ph0m5n3Aha6+0YAd/8CGBvNa25mvdz95cyHKSIidU26jSReAvYvYd5+0XwREZFKSzdBWSnzmgMFlYhFRESkSJlVfFG13lEJRReb2YlJizUG+gDvZS40ERGpy8pzDao7oYUdhBtlzwK2Ji3zHbCEMAyHiIhIpZWZoNx9FDAKwMw+BU51dw1GKCIiVSrdZuZ7VVUgIiIiiSoy3EZjQs/j7dlxuA1394cyEZiIiNRt6d4HdSQwGSipR1YHlKBERKTS0m1mfh+hG6KuwE7uXi/pr37mQxQRkboo3Sq+/YDT1UhCRESqWrpnUO8CbasiEBERkUTpJqjLgWvM7GdVEYyIiEihdKv4ZgJNgVlm9h2wIXkBd6/YyFQiIiIJ0k1QD1DCsOsiIiKZlO6NusOrKA4REZFi0r5RF8DMdgUOBPYAnnf3ddENvN+5+/ZMBigiInVTWo0kzKyBmd0NLAP+TRj6vbD7o8nArZkNT0RE6qp0W/GNBC4BhgE/ovj4UFOBUzIUl4iI1HHpVvENBm5098fMLLnXiE8ISUtERKTS0k1QLQiJKJVGgLo6EqmgD21tqfO7VlMcInGRbhXf+0D/EuadBCysXDgiIiJBumdQdwCTzawJ8DfCPVGHmNlpwKVAvwzHJyIidVRaZ1DuPhU4BzgOeJ7QSOJhYAhwnru/kOkARUSkbkr7Pih3nwhMNLN9gVbAWuBDd1cPEyIikjEVulEXwN0/Aj7KYCwiIiJF0r1R91EzG1/CvKfNbGxmwhIRkbou3VZ8vQk9RqQyGTihcuGIiIgE6Sao1oRrTqmsAzTUhoiIZES6CeozoFcJ83oR+ugTERGptHQT1DjgBjO70syaA5hZczO7Arie0ORcRESk0tJtxXcX0An4M3CfmX0LNCPcDzUmmi8iIlJp6Q5YuB242MxGAUcDLYGvgFlRs3MREZGMKHeCigYk/AYY4O5/Bz6ssqhERKTOK/c1KHffDKwGtlZdOCIiIkG6jST+ClxlZg2rIhgzq29mb5nZtGh6LzN7w8w+NrMJZtaoKp5XRETipyLjQR0ILDWzF4FVhB7NC7m731CJeH4JLAZyoum7gNHuPt7M/gJcBDxUie2LiEgNkW6COgPYEj3umWK+AxVKUGbWHuhDGFb+WjMz4BhC7+kAjwPDUYISEakT0m3Ft1dVBQL8kXAv1c7RdEvga3cvvOa1DGiXakUzGwoMBejQoUMVhigiItUl3WtQVcLM+gKr3X1BRdZ39zHunu/u+a1bt85wdCIikg1pJygzOzhqsPCJmW0xs59E5SPN7KQKxtED6GdmS4HxhKq9PwEtzKzwLK89sLyC2xcRkRom3eE2TgIWAG2BJ4DE1nxbgF9UJAh3/427t3f3jsBAwo2/5wIvAWdGi50PTK3I9kVEpOZJ9wzq98A4d/8ZoTFDoreBQzIS1Q9uIDSY+JhwTeqRDG9fRERiKt1WfJ2B66LHyUO8rwd2q2xA7j4bmB09/i9wWGW3KSIiNU+6Z1CrgR+VMO8A4PPKhSMiIhKkm6DGAyPM7MiEMjezfQnVcU9mLDIREanT0q3iuwXYH3gZWBGVTSU0mpgB/C5zoYmISF1WrgRlZk2Ak4GOwNPAU4Quj1oRhoB/0d1nVlGMIiJSB5WZoMzsR8C/CMmp0HrCsBsvVFFcIiJSx5XnGtTdwHZC33tNCY0h3kJ94omISBUqT4I6ArjZ3V91983uvhi4FNjTzPKqNjwREamrypOg8oD/JpV9AhihcYSIiEjGlbeZefJNuSIiIlWqvM3MXzCzVEO9v5hc7u65lQ9LRETquvIkqNuqPAoREZEkZSYod1eCEhGRaheLAQtFRESSKUGJiEgsKUGJiEgsKUGJiEgsKUGJiEgsKUGJiEgsKUGJiEgsKUGJiEgsKUGJiEgsKUGJiEgsKUGJiEgsKUGJiEgsKUGJiEgsKUGJiEgsKUGJiEgsKUGJiEgsKUGJiEgsKUGJiEgsKUGJiEgsKUGJiEgsKUGJiEgsKUGJiEgsKUGJiEgsxSJBmdkeZvaSmX1gZovM7JdR+W5mNtPM/hP93zXbsYqISPWIRYICtgK/cvf9gcOBK81sf+BG4EV33wd4MZoWEZE6IBYJyt1XuPvC6PEGYDHQDugPPB4t9jhwanYiFBGR6haLBJXIzDoCXYE3gDbuviKatRJoU8I6Q81svpnNX7NmTbXEKSIiVStWCcrMmgOTgavdfX3iPHd3wFOt5+5j3D3f3fNbt25dDZGKiEhVi02CMrOGhOT0pLs/ExWvMrO8aH4esDpb8YmISPWKRYIyMwMeARa7+x8SZj0LnB89Ph+YWt2xiYhIdjTIdgCRHsB5wHtm9nZUdhNwJzDRzC4CPgPOzlJ8IiJSzWKRoNz9FcBKmH1sdcYiIiLxEIsqPhERkWRKUCIiEktKUCIiEktKUCIiEktKUCIiEktKUCIiEktKUCIiEkuxuA9KpNZ76ffZjkCkxtEZlIiIxJISlIiIxJISlIiIxJKuQYnUEA++/WCp86845IpqikSkeugMSkREYkkJSkREYkkJSkREYkkJSkREYkkJSkREYkkJSkREYkkJSkREYkkJSkREYkk36orUFJ/OKX2+btSVWkZnUCIiEktKUCIiEktKUCIiEktKUCIiEktKUCIiEktqxSdSS2g4DqltlKBEMqDM5FBNcYjUJkpQUuvF4cziwa/fLXOZ3ao8CpGaRQlKpI4oK1GDqgElXtRIQkREYkkJSkREYklVfFLjlafqqiY8R6WV1VffXj2rJ44YiMN1R6k8JSiRTCgrOYhI2pSgpM57cMqgsheqQ2cfpSnzzGTdN2Vv5OjfZCgaqe1in6DM7ETgT0B94GF3vzPLIUlNo7OboBzHoQZUZEodEusEZWb1gQeA3sAyYJ6ZPevuH2Q3MpFaqrLJvMXBZS5S1hnrFac9XbkYpNaIdYICDgM+dvf/ApjZeKA/oASVAZW98F+eC801onFBeegsrFzKc0NymduIQQLTPWPxYO6e7RhKZGZnAie6+8XR9HlAd3cflrTcUGBoNLkf8GG1BrqjVsCXWY6hJtBxKpuOUfnoOJUtzsdoT3dvnVwY9zOocnH3McCYbMdRyMzmu3t+tuOIOx2nsukYlY+OU9lq4jGK+426y4E9EqbbR2UiIlLLxT1BzQP2MbO9zKwRMBB4NssxiYhINYh1FZ+7bzWzYcALhGbmj7r7oiyHVR6xqW6MOR2nsukYlY+OU9lq3DGKdSMJERGpu+JexSciInWUEpSIiMSSElQVM7NfmZmbWatsxxJHZjbKzJaY2btmNsXMWmQ7prgwsxPN7EMz+9jMbsx2PHFjZnuY2Utm9oGZLTKzX2Y7pjgzs/pm9paZTct2LOWlBFWFzGwP4Hjg82zHEmMzgQPd/WDgI0A9iVKsm6+TgP2BQWa2f3ajip2twK/cfX/gcOBKHaNS/RJYnO0g0qEEVbVGA9cDaolSAnef4e5bo8nXCfe6SUI3X+7+HVDYzZdE3H2Fuy+MHm8gfPm2y25U8WRm7YE+wMPZjiUdSlBVxMz6A8vd/Z1sx1KDXAg8n+0gYqId8EXC9DL05VsiM+sIdAXeyG4ksfVHwo/l7dkOJB2xvg8q7szsX0DbFLN+C9xEqN6r80o7Tu4+NVrmt4QqmyerMzap+cysOTAZuNrd12c7nrgxs77AandfYGZHZTuedChBVYK7H5eq3MwOAvYC3jEzCNVWC83sMHdfWY0hxkJJx6mQmQ0B+gLHum7MK6RuvsrBzBoSktOT7jWiNtkAAAyPSURBVP5MtuOJqR5APzM7GWgM5JjZ/+/uP89yXGXSjbrVwMyWAvnuHteehLMmGpDyD8DP3H1NtuOJCzNrQGg0ciwhMc0DzqkhPalUCwu//h4H1rr71dmOpyaIzqCuc/e+2Y6lPHQNSrLtfmBnYKaZvW1mf8l2QHEQNRwp7OZrMTBRyWkHPYDzgGOi987b0VmC1BI6gxIRkVjSGZSIiMSSEpSIiMSSEpSIiMSSEpSIiMSSEpSIiMSSElQlmdnwqLdyN7PtZrbOzOaZ2Ugza5u0bMdoub4JZc3MbLyZfRXNGxKVX2Jmn5rZVjObXb17VTOY2Tgzm58wfZiZDU+x3HAzq7P3oJnZfDMbV8Yy12e6lwEz2zc69i2SyodE7/XmGXqee6J7Datk++V4/qFmdmoGt1et3wlmdnbhc8SNElRmfAMcAfwUGAg8Q7g/4z0zOzRhuRXRcq8klF0OnAIMjeZNjxLbQ8BU4GfAFVW9AzXU7cCQhOnDgFuzE0qNdz1wVIa3uS/h9ajuIVSmEz5LBdX0fEOBjCUoqv874WyKf45iQ10dZcZWd389YfoFM3sIeBkYb2ad3X2bu28h9NidqDPwobtPLiwwsyOB+sCj7v5uZQIzsybuvqky24grd/8k2zGkEvVwsJO7b852LJlWE95PUY8kpfZKEvP9qNLvhBrF3fVXiT9gOPBlCfNOJAy1cWI03TGa7htNL42mE/+GpygbEi3fGLib0Mv1FuAd4OSk51wK3AvcQugB+/uovB5wI/BxtO5HwPlJ684GJgHnRMutJ/Qu3j5puSZRHJ9F2/oU+H3SMhcDi6L5nwHXl3EcbwM+SphuBnwPLEwoa0Xojbl3ND0OmB89HpLiuM1OfI0IvV2/Tvhl/RbQsxyvb4foGGyK9nNIdIxmJ78HgCMJXRJtAc6L5p0NvBeVfQGMBBqU9f6J4h+W9LreA1wTva7rCENwtEha70DgVWAzoQeKfsB8YFwp+7g0xbE7KiGOawm9Ya8hDAFSFE/Sdgpfg+aEs7HkbS5NWu4gwnhg3wJLgNPL8Xq0AJ4CNhJqJH4bHZelqeJI+tydCzwBfA38K5q3GzAGWBUds9eA7knPWZ8wTtlH0eu4rPB4Ej4zKT+vJcRf6me4hNdieEnPUdb2ErZ7CeF9uDna10nALoTP0A7Pl43v0pTHK9sB1PQ/Sk9QjQlfssOj6cIPSmGC6kqojlhMGHDtcEKnoFdEy50TlbWOlp8GrCZUARxPGNtlK3BIwnMujT64/yJ8OZ0elT8QfaivB44D7gK2FcYSLTM7eqO/Rhh7aGD0Zv5HwjIGzAA2RNs6FhgMjE1Y5tfRfo8EehMS4xYSvnBTHKvjon1uE033JiSFbUBOVHZatL+FXzzj+CFBtSZ8UXnCsdw/4TUqAN4FLiAMAvg64Qu3aSkxGfB2dEwHRc+/MDpGs5PeAwXAJ8ClwNHAftFr5IT+4k6MjtcW4C9lvX9InaA+j94DJxOqfzYCDyYs04TQb987wOmEL+T/Ru+ZcaXsZ1fCl/bDCccuJyGOFcCEaB9OTointASVA/wqmj4t2mbXpOXeA34RHafngO9I+jGUItYphOR8CaEa7N+EhLE0VRxJn7sVhM9Bb+AYYKfo9fwv4T18IqEKbQPQNmF7D0ex3RGtO4DQ9RSEwSQXEz7HhceudSnxl/oZpmq+E24m/LC7P9rH04FHCMO3dAJmRceh6Pmy/b1aFHu2A6jpf5SSoKL5K4CHoseFH5TEpDCO6Es2oeyoaLkDE8qOjcp+lrTsy8DfEqaXRs/ZOKFs7+gNen7Suk8A8xKmZxOup+2aUHZ19LxNoukToul+JexvDuGL89ak8hHASqB+CesVnjGdmbD8ZOB//HAG+oekeIsdO0LfdV7Ca+TAMQllh5BwdltCTH2iZbollLWL4pydYvv9k9Z/HXgpqex6QtJtX9r7h9QJ6hOKn339EViZMH1FFFv7hLIe0bbGlfE+/pIUv5yjdRemKF9KKQkqmu4bTXcsYbkLE8paEr5YLyslxgOi9QYklDUH1lK+BDUlaXsXERLPPgllDaLjPCqa7hyte1UpcZV6hpqwXHk/w8Xe11HZUVTgO4FwxlkA/KGUuIrVCMTpT40kqp5laDvHEb7gXzWzBoV/wItAftKyL3rx6x/HEhLUlBTrHhINL15onruvS5j+IPpfOFjeMYTeo58tIc4jCMnmb0nPNQtoQwkj5rr7t4RfcT2jol6ED9qcpLI5JTxvWb4jJOBChftV2gi+3QgJYF5CnMuBBSmWdRIGW4yO6U+AvyUtN4FQ3XpEeQNP8JL/MPowhH3IjYacgNBIZIG7L0uI91XCL+zK+Ecl1y/JjMIH7v4VIc6yXg8IZzmF620kVBOWx/Sk6eMIr+WnCe9TCGdlhZ+po6P/48r5HKVJ5zOcqe0dQTizfqzy4Vc/NZKoQmbWmPDLcFUGNteKMOjf9ynmbUuaTn6+VoR69G9K2HYeoZoEQlVPou+i/42j/y0JZ2ilxQnh+lMqexCuSaUyh9AzdSOgO+HaxzbgbDPbmXDWM7KU5y7NBncvGk3U3b+LxupqXPIqtCX1xfY1hB7YE63zMDR7oVZAQ3Z8LQqndytP0ElSvTZGqKr6Poo3VTKqbILKxPs3lVT7U9brscF3bHxS3v1L9bk4nNSfqcIGOC2Bbz0zAyGm8xnO1PZaRv9L+8zGlhJU1TqacIznZmBbawnXF8rTnNVTrLuVUN2TasjndL7AviIktJKsjf73JfUX24elrDuH0AjgWMKX1duED9o9hGNZn+JN9KvaSsK1rWStCRebEyUf8y8JXxy5SeVtov+Fx2kz0ChxATPbNe1Ig5WEKqlkyTGkK3nfIEXcQEXjLq+VwM5m1jgpSZV3/1J9LuYTrt8k2xL9/wpoZmY5GUhS6XyGM7W9r6L/eYT3ZI2iBFVFopsT7yK0hvtXBjb5IuGi80Z3X5LmurMIX+67uHt5q0NKi+N6M+vr7tNSzJ9LaNywu7snV6mUZQ7hjOBG4FV3325m70Xb+xWwxEsf1PA7CGeuKX5lV8Q84NZoJOQ3o223Aw4ltJQrkbtvM7MFwFmE+1cKnU34kVD4o2UZ4Uu3XVR9COFid0XjPdfM2hdW85lZD8r3BV7W2UuyZUCXpLLkuJPPviursKq1P6GqtHC4996EFqfpepEQ8+fuXtKPtFnR/8GERgaplPfYVeYzXNHtFX4ezweuK2GZdF/7aqMElRkNzOzw6PHOhC+wy4GmhIvwFTl9TzaTMHjdTDO7i1CFlkOo9mrs7r8paUV3/zAaCHC8md1N+NXYmHDReV93v7gCcTxlZiMI143ygF7ufqm7fx315vAnM9uTcB2pHuGmzaPd/bRS4lxrZh8QrjX9JirbbmavEhosjC0jtsIP6S/NbBaw3t1LO2Mryz8ILeImmtlvCB/0WwlnhqnORJPdSrgn7jFCk/CDCDcXj024TvTPaLuPmtm9wF7AZRWM9zFCi63p0WvQJHq+8vxyXgL0MbN/Ehq5fOjuG0pZfgrwZzO7iZA4ziC8nxIVHvtLzWw8UODu75V3Z5K5+yIzexZ4yMxyCNVWv6biN+Q+QTjWs83sHkJrvpaEa3kr3X109NkZA9xrZrmE93MLQmOegdF2lgAnmNkJhDOWT6Nraskq/BkuQZnbiz6PtwMjo6rzfxCqhPsAt0U/ipYA/aPeMJYB/3P3/6UZS9XIdiuNmv5H8XsUthPq1ecTrpW0TVq2IxVsxReV70S4X+hjwq+elYQvuD4JyywlqXVVVG6EFnmF9yatIVwMHpywzGxgUlmxEL747iG8mQvvgxqZtN7PCRegNxGaBb8BXFuO4/lQ9Hw/TSi7ISobnLRssWMX7ePdhJZ/20m6DyrFcxVrKVdCPHtGx3gz4drZUMLF/b8nvQdKutVgAKE59XfR8Sp2H1S0zEnR61JAOIvskhxbqteVpNZqUdnBhNsEthASxKmUo5UZ4UfV64R7kpzi90HtcIwI19f+EL0H1wF/io5Ncjy/io7bVna8D6p50jZ32McUz7srIdl/S/ih8H+U/z6ovim2t0sU+xcJr9EzQI+EZeoDNxESWOEyjybM/xGhluQbyr4Pqjyf4XFk8DshWu5SQqOaLdEyE/nhVoJWhB8ca4nZfVAaUVckDWa2C+GL6n53vzXb8YjUZqriEymFmV1GOBv7D6FxxLWEX62PZjMukbpACUqkdJsJVYx7Eqo/3gSOc/eSmsqLSIaoik9ERGJJPUmIiEgsKUGJiEgsKUGJiEgsKUGJiEgsKUGJiEgs/T/jGQ+tLBFVwAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "       ind: 0.0 ($\\pm$ 0.0)\n",
      "  manifold: 0.2 ($\\pm$ 0.3)\n",
      "       asv: 0.4 ($\\pm$ 0.6)\n",
      "      flow: 0.0 ($\\pm$ 0.0)\n"
     ]
    }
   ],
   "source": [
    "errors = {\n",
    "    'ind': get_error(gt, theta_ind),\n",
    "    'manifold': get_error(gt, theta_manifold),\n",
    "    'asv': get_error(gt, theta_asv),\n",
    "    'flow': get_error(gt, theta_flow)\n",
    "}\n",
    "\n",
    "bins=np.histogram(np.hstack(list(errors.values())), bins=40)[1] # get the bin edges\n",
    "for k, v in errors.items():\n",
    "    plt.hist(v, bins, label=k, alpha=0.5, weights=np.zeros_like(v) + 100. / v.size)\n",
    "\n",
    "plt.legend(fontsize=15)\n",
    "plt.xlabel('Difference with ground truth direct effect', fontsize=15)\n",
    "plt.ylabel('Percentage', fontsize=15)\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "for k, v in errors.items():\n",
    "    print(f\"{k:>10s}: {np.abs(v).mean().round(1)} ($\\pm$ {np.abs(v).std().round(1)})\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAEYCAYAAAAaryJBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deZgU1dn38e/NJrvsMCKIQSO4QhxERQiiuIu7oiaIS3B9jBqjZjEqhjwqMSRGYwIu6PuogCDBoEZQRFHRgLgiYDSiQFhUUMABFLjfP07N2NP0THfP9Ez1zPw+1zXXdJ1a+u5a+u46deqUuTsiIiLVrV7cAYiISN2kBCQiIrFQAhIRkVgoAYmISCyUgEREJBYN4g4g19q1a+fdunWLOwwRkTrnjTfe+Nzd22c6fa1LQN26dWP+/PlxhyEiUueY2SfZTK8qOBERiYUSkIiIxEIJSEREYqEEJCIisVACEhGRWCgBiYhILKq1GbaZPQCcAKxx932jsjbARKAbsBQ4093XmZkBfwKOA4qA4e6+oDrjFZGKWb9+PWvWrOHbb7+NOxTJkYYNG9KhQwdatmyZs2VW931A44G7gYcTym4Annf328zshmj4euBYYM/ory9wb/RfRPLY+vXrWb16NZ07d6ZJkyaE35JSk7k7mzZtYsWKFQA5S0LVWgXn7i8Ba5OKTwIeil4/BJycUP6wB68BrcysoHoiFZGKWrNmDZ07d6Zp06ZKPrWEmdG0aVM6d+7MmjVrcrbcfOgJoaO7r4xerwI6Rq87A8sSplsela0kiZmNAEYAdO3ateoilby3YdYL5Y5vMejwaoqk7vr2229p0qRJ3GFIFWjSpElOq1XzqhGCh8ezZv2IVncf6+6F7l7Yvn3G3RCJSBXRmU/tlOvtmg8JaHVx1Vr0v/j8bgXQJWG6XaMyERGpBfIhAT0JnBe9Pg+YllA+zIKDga8SqupERKSGq+5m2I8BA4F2ZrYcuAm4DZhkZhcCnwBnRpM/TWiC/SGhGfb51RmriOTWc++vjuV9j9y7Y/qJkgwfPpz33nuv0j3rT58+nRNPPJGPP/4YPSZmR9WagNz97DJGHZFiWgcur9qIRER2dOONN7Jp06a4w6j18qEVnIhIXunevXvcIdQJ+XANSEQkrwwfPpzCwkIAxo8fj5nx7rvvMnjwYJo1a0aPHj144oknSs3j7tx888106NCBFi1aMGzYMNavXx9H+DWGEpCISAbOOecchgwZwtSpU9lzzz0ZOnQoy5cvLxl/1113MXLkSEaMGMHkyZNp0qQJ1113XYwR5z9VwYmIZODqq6/mggsuAODAAw+kY8eOTJ8+nUsuuYRt27Zx++23c/HFF/Pb3/4WgKOPPprBgweXdF8jO9IZkIhIBo466qiS123btqVDhw4lZ0DLli1j5cqVnHTSSaXmOfXUU6s1xppGCUhEJAOtWrUqNdyoUSM2b94MwKpVqwDo0KFDqWmSh6U0JSARkUrq1KkTwA4ddeay487aSAlIRKSSunTpQqdOnZg2bVqp8uSWclKaGiGIiFRS/fr1ue6667j22mtp164d/fv3Z8qUKSxatCju0PKaEpCIVIuKdIlTk1x11VWsXbuWv/71r/zxj39kyJAh3HHHHZx77rlxh5a3LPR4U3sUFhZ6ZftvkppLzwOK36JFi+jZs2fcYUgVKW/7mtkb7l6Y6bJ0DUhERGKhBCQiIrFQAhIRkVgoAYmISCyUgEREJBZKQCIiEgslIBERiYUSkIiIxEIJSEREYqEEJCKSZ5YuXYqZMX369JKyr7/+mqFDh9K2bVvMjPHjx2e0rNmzZ2NmvPfee+VOd+2119KtW7dKRJ099QUnItVjyTPxvO9ex8bzvpVQUFDA3Llz6dGjR0nZvffeyz/+8Q8efvhhOnfuTPfu3WOMMDeUgERE8sxOO+3EwQcfXKps8eLF7LXXXpx22mkxRZV7qoITEUkyfPhwCgsLeeqpp9h7771p2rQpxx9/PGvXruXDDz/k8MMPp1mzZhQWFvLOO++UzHfnnXfSp08fdt55Zzp27MiJJ57Ihx9+WGrZAwcO5PTTT+fRRx9ljz32oGXLlhx77LElj/eGHavgunXrxv3338+bb76JmWFmJdPOmjWLvn370rhxYzp27Mhll13Gxo0by/18X375Jeeccw7NmzenoKCAUaNG5WK1ZU0JSEQkhU8//ZTf/OY3/Pa3v2Xs2LG8+uqrjBgxgqFDhzJ06FAmT57M1q1bGTp0KMVPFVi+fDlXXHEF06ZNY9y4cWzbto1DDz2Ur776qtSyX3/9de6++27uvPNOxo4dy4IFCxgxYkSZsUydOpXjjjuOHj16MHfuXObOnQvAwoULOeaYY2jXrh1Tpkzhlltu4dFHH+X0008v97Odf/75PPPMM4wZM4axY8cyY8YMJkyYUMk1lj1VwYmIpLB27Vrmzp1bcq3lnXfeYfTo0Tz00EMMGzYMAHfn+OOPZ/HixfTs2ZMxY8aUzL9t2zYGDx5Mhw4dmDZtWsk8AOvXr+epp56idevWAKxatYqrr76aTZs20aRJkx1i6d27N+3bt2f16tWlquZuvfVWdtttN5588knq168PQJs2bTjrrLOYO3cuhxxyyA7LWrhwIX//+9+ZMGECZ511FgCHH344Xbt2pWXLlpVdbVnRGZCISArdunUrdaF/jz32AGDQoEE7lK1YsQKA1157jcGDB9O2bVsaNGhA06ZN2bhxIx988EGpZffp06ck+QDsvffepZaTqX/961+ccsopJckH4LTTTqNBgwa8/PLLKeeZN28eACeddFJJWfPmzRk8eHBW750LSkAiIim0atWq1HCjRo12KC8u27x5M59++ilHHXUU7s7f/vY3XnnlFebNm0eHDh3YvHlzRstOni6dlStX0rFj6SfN1q9fn7Zt27J27dqU86xatYoWLVrQuHHjUuUdOnTI6r1zQVVwIiI58M9//pOioiKmTZtGs2bNANi6dWuZiSAXCgoKWLNmTamybdu28cUXX9CmTZuU83Tq1IkNGzawefPmUkkoeTnVQWdAIiI5sGnTJurVq0eDBt/9rp80aRJbt26tsvfs27cvU6dOZdu2bSVlTzzxBFu3buWwww5LOU+fPn0AmDZtWknZxo0bmTlzZpXFWRadAYmI5MCgQYPYtm0b559/PhdeeCELFy7k97///Q7Vbbn061//mt69e3PyySdz6aWXsnz5cq6//nqOPvrolA0QAPbZZx+GDBnCpZdeyvr16ykoKGD06NE0bdq0yuIsixKQiFSPGtgjQTb2228/xo8fz80338zUqVM54IADePzxx0tamlWFffbZh2eeeYZf/vKXnHrqqbRs2ZKzzz6bO+64o9z5xo8fz6WXXspVV11F8+bNufzyy+nTpw+TJ0+uslhTseL267VFYWGhz58/P+4wJCavTrqr3PGHnnllNUVSdy1atIiePXvGHYZUkfK2r5m94e6FmS5L14BERCQWSkAiIhKLvElAZna1mS00s/fM7DEza2xmu5vZ62b2oZlNNLNGcccpIiK5kRcJyMw6A1cChe6+L1AfGArcDoxx9z2AdcCF8UUpIiK5lBcJKNIAaGJmDYCmwEpgEFDcLOMh4OSYYhMRkRzLiwTk7iuA3wOfEhLPV8AbwJfuXnwX13Kgc6r5zWyEmc03s/mfffZZdYQsIiKVlBcJyMxaAycBuwO7AM2AYzKd393Hunuhuxe2b9++iqIUEZFcyosEBBwJfOzun7n7t8ATQD+gVVQlB7ArkF1XsSIikrfyJQF9ChxsZk0tPOrvCOB94AWg+MlK5wHTyphfRERqmLxIQO7+OqGxwQLgXUJcY4HrgWvM7EOgLXB/bEGKiEhO5U1fcO5+E3BTUvF/gINiCEdEcmz2stmxvO/ALgNjeV9JLy/OgEREpO5RAhIRkVgoAYmIJJk7dy5DhgyhoKCAZs2a0atXLx555JGS8V9++SUXXXQRu+yyC40bN6Zr16785Cc/AWD27NmYGQsXLiy1zHXr1tGoUSPuu+++av0s+SxvrgGJiOSLTz75hH79+nHJJZfQuHFjXnnlFc4//3zq1avH2WefzTXXXMOrr77KmDFj6NSpE8uWLeOll14CYMCAARQUFDBp0iRuueWWkmVOnToVgNNOOy2Wz5SPlIBERJIMHTq05LW7M2DAAJYvX864ceM4++yz+de//sXll19e6mFzP/rRjwCoV68eZ5xxBhMnTiyVgCZOnMhRRx1F69atq++D5DlVwYmIJFm3bh1XXnklu+22Gw0bNqRhw4aMHTuWDz74AIBevXoxevRo/vKXv5SUJTrrrLNYsmQJb7/9NgCff/45s2bNqtKno9ZESkAiIkmGDx/OxIkT+fnPf86MGTOYN28eF1xwAZs3bwbg7rvv5uSTT2bkyJHstdde7LnnnkyYMKFk/kMOOYSuXbsyceJEAKZMmUKDBg04+WT1p5xICUhEJMHmzZuZPn06t9xyC1dccQWDBg2isLCQ7du3l0zTqlUr7rrrLlatWsXbb79N3759Offcc3n//fcBMDPOPPNMJk2aBITqt2OPPZYWLVrE8pnylRKQiEiCLVu2sH37dnbaaaeSsg0bNvDkk0+mnH7//fdn9OjRbN++ncWLF5eUDx06lI8++ojp06fz4osvlrquJIEaIYiIJNh5553p06cPI0eOpGXLltSrV4/bbruNnXfemfXr1wNw2GGHccopp7DvvvtiZowbN45mzZpx0EHfddxy4IEHssceezBixAiaNGnCCSecENdHyltKQCJSLWpSlziPPvooF198McOGDaNt27ZcccUVFBUVcffddwPhGs/48eNZunQp9evXp3fv3jzzzDPsuuuupZZz1llnMWrUKIYOHUrTpk3j+Ch5zdw97hhyqrCw0OfPnx93GBKTVyfdVe74Q8+8spoiqbsWLVpEz5494w5Dqkh529fM3nD3wkyXpWtAIiISCyUgERGJhRKQiIjEQglIRERioQQkIiKxyDoBmdn+ZjbRzD4ysy1m9oOofJSZHZv7EEVEpDbKKgFFCeYNoBPwMNAwYfQW4H9yF5qIiNRm2Z4B/S8w3t1/CIxKGvcW0CsnUYmISK2XbQLqAUyMXiffwboeaFPpiEREpE7INgGtAb5Xxrh9gE8rF46ISH4YOXIknTt3pl69enTr1g0z47333os7rFol277gJgAjzex9YG5U5mb2feB64P5cBicitceGWS/E8r4tBh2e9Tzz58/npptu4ne/+x0DBw7k/fff56KLLqqC6Oq2bBPQjcDewIvAqqhsGqFRwgzgd7kLTUQkHsWPVbj88stp2bIlW7ZsiTmi2imrKjh33+LuJwBHAQ8B9wGPAse7+wnu/m0VxCgiUm2GDx/Oj3/8YyA8msHMUk5XVFTElVdeSadOnWjcuDF9+vRhxowZJeMffPBBmjVrxrfffve1uMsuu9C2bVuKO4Hevn07rVq1Yty4cVX4ifJXhW5Edffn3f2X7j7C3W9w95m5DkxEJA433ngjv/71rwGYNWsWc+fOLXkOUKKf/OQnPPjgg/zqV79i6tSpdOnSheOPP56XX34ZgP79+1NUVMSCBQsA+Pe//82aNWvYsGFDyZNT3377bb766iv69+9fTZ8uv2RVBWdmXcsZvR1Y7+47bikRkRqie/fudO/eHYA+ffrQvHlzZs+eXWqaRYsW8dhjj/Hggw9y3nnnAXD00Uez//77c+utt/Lss8+yxx57UFBQwJw5c+jbty9z5szhgAMOoFGjRsyZM4d99tmHOXPm0L59e3r06FHdHzMvZHsGtBT4uIy/T4B1ZvaxmV2dyyBFRPLJvHnzcHfOOOOMkrJ69epxxhlnlJwBQTgLmjNnDgAvvfQSAwYMYMCAAaXKDjvssOoNPo9km4DOAZYD/wQuA86I/j8LrAAuB54H7lASEpHaauXKlTRv3nyHp5x27NiRoqKikkYL/fv35+WXX8bdmTNnDv379y+VlF5++eU6W/0G2beCOxJ40t2Tu9z5m5n9GTjU3YeZ2UbgEmBMLoIUEcknBQUFbNy4kaKiolJJaPXq1TRt2pSddtoJCAlo7dq1zJw5k48//pj+/fvToEEDVqxYwYwZM1i9enWdTkDZngGdQWh2ncqTwEnR62eA3SoalIhIPuvTpw9mxuTJk0vK3J3JkyeXqlLbb7/9aNWqFaNGjaJHjx60b9+e1q1bs++++zJq1CiaN29O79694/gIeSHbM6DNQD/guRTj+kXjAQz4uhJxiYjkrZ49e3L22WdzxRVXsGHDBrp37864ceNYvHgx9957b8l09erVo1+/fjz11FNcfPHFJeX9+/fnnnvuYfDgwdSvXz+Oj5AXsk1AY4Ebzawt8A/gM6A94cznEr67EfVQ4O1cBSkiNV9FeiTIZ+PGjeP6669n5MiRfPnll+y3335Mnz59h0YF/fv356mnnmLAgAGlyu6555463QABwIpviMp4htC44OeE3g+ccLazChjt7mOiafYBvnb3pTmNNgOFhYU+f/786n5byROvTrqr3PGHnnllNUVSdy1atIiePXvGHYZUkfK2r5m94e6FmS4r2zMg3H2Mmf0J6EJIQquAZe6+PWGahdkuV0RE6paK9oSw3d0/cffXo//b089VPjNrZWaTzWyxmS0ys0PMrI2ZzTSzf0f/W1f2fUREJD9kfQZkZi0I13y+DzROHu/u11Uwlj8B/3T3082sEdAU+CXwvLvfZmY3ADcQet0WSenjz8tv+3JoNcUhIull2xVPd+BVoAnQjNAIoU20nHXAV0DWCcjMdgYGAMMB3P0b4BszOwkYGE32EDAbJSARkVoh2yq4McA8oCOh8cFxhGT0I2AjcFYF49idkMweNLM3zew+M2sGdHT3ldE0q6L33YGZjTCz+WY2/7PPPqtgCCKSK9k2bpKaIdfbNdsEdBDwV6D44RiN3H2buz8K3EmoRquIBsAPgHvdvTfhHqIbEifw8MlTfnp3H+vuhe5e2L59+wqGICK50LBhQzZt2hR3GFIFNm3aRMOGDXO2vGwTUGNCj9fbgbXALgnj3gMOqGAcy4Hl7v56NDyZkJBWm1kBQPR/TQWXLyLVpEOHDqxYsYKioiKdCdUS7k5RURErVqygQ4cOOVtuto0QPuC7LnbeBC4xs6eBbcCFwH8rEoS7rzKzZWa2l7svAY4A3o/+zgNui/6X1Q2QiOSJli1bAvDf//631MPYpGZr2LAhHTt2LNm+uZBtApoA9AL+H+Hx3M8C6wnPAmpASBIV9T/AI1ELuP8A5xPO0CaZ2YWExz2cWYnli0g1admyZU6/qKR2yioBufsfEl6/Zmb7AscQGiLMcvf3KhqIu78FpLqD9oiKLlNERPJXts2wBwAL3H0jgLsvA8ZF45qb2QB3fyn3YYqISG2TbSOEF4C9yxi3VzReREQkrWwTkJUzrjlQVIlYRESkDklbBRdVuw1MKLrIzI5JmqwxcDzwbu5CExGR2iyTa0B9CS3UINwIegawNWmab4DFhMc0iIiIpJU2Abn7aGA0gJl9DJzs7nrYnIiIVEq2zbB3r6pARESkbqnI4xgaE3qu3pUdH8fg7n7vjnOJiIiUlu19QIcBU4Cyevx0QAlIRETSyrYZ9l2EbnJ6Azu5e72kv/q5D1FERGqjbKvg9gJOVSMEERGprGzPgN4BOlVFICIiUrdkm4AuBa42sx9WRTAiIlJ3ZFsFNxNoCswys2+ADckTuHvunlYkIiK1VrYJ6B7KeCy2iIhINrK9EfXmKopDRETqmKxvRAUws9bAvkAX4Bl3XxfdoPqNu2/PZYAiIlI7ZdUIwcwamNkdwHLgRcKjuYu755kC3JTb8EREpLbKthXcKOAnwBXA9yj9fKBpwIk5iktERGq5bKvghgE3uPuDZpbc68FHhKQkIiKSVrZnQK0IiSaVRoC64hERkYxkm4DeA04qY9yxwILKhSMiInVFtlVwvwWmmFkT4HHCPUG9zOwU4GJgSI7jExGRWiqrMyB3nwacAxwJPENohHAfMBz4sbs/m+sARUSkdsr6PiB3nwRMMrPvA+2AtcASd1cPCSIikrEK3YgK4O4fAB/kMBYREalDsr0R9QEzm1DGuMfMbFxuwhIRkdou21Zwgwk9HqQyBTi6cuGIiEhdkW0Cak+45pPKOkCPYhARkYxkm4A+AQaUMW4AoY84ERGRtLJNQOOB683scjNrDmBmzc3sMuA6QpNsERGRtLJtBXc70B34M3CXmX0NNCPcDzQ2Gi8iIpJWtg+k2w5cZGajgcOBtsAXwKyoWbaIiEhGMk5A0QPnvgLOcve/A0uqLCoREan1Mr4G5O6bgTXA1qoLR0RE6opsGyH8DbjSzBpWRTBmVt/M3jSz6dHw7mb2upl9aGYTzaxRVbyviIhUv2wbIbQC9gWWmtnzwGpCj9jF3N2vr0Q8PwUWAS2j4duBMe4+wcz+ClwI3FuJ5YuISJ7INgGdBmyJXvdPMd6BCiUgM9sVOJ7w2O9rzMyAQYTetwEeAm5GCUhEpFbIthXc7lUVCPBHwr1ELaLhtsCX7l58zWk50DnVjGY2AhgB0LVr1yoMUUREciXba0BVwsxOANa4+xsVmd/dx7p7obsXtm/fPsfRiYhIVcg6AZnZ/lGDgI/MbIuZ/SAqH2Vmx1Ywjn7AEDNbCkwgVL39CWhlZsVnabsCKyq4fBERyTPZPo7hWOANoBPwMJDYGm4L8D8VCcLdf+Huu7p7N2Ao4cbWc4EXgNOjyc4DplVk+SIikn+yPQP6X2C8u/+Q0Fgg0VtAr5xE9Z3rCQ0SPiRcE7o/x8sXEZGYZNsKrgdwbfQ6+RHc64E2lQ3I3WcDs6PX/wEOquwyRUQk/2R7BrQG+F4Z4/YBPq1cOCIiUldkm4AmACPN7LCEMjez7xOqyx7JWWQiIlKrZVsFdyOwN/ASsDIqm0ZolDAD+F3uQhMRkdosowRkZk2A44BuwGPAo4QuedoRHtH9vLvPrKIYRUSkFkqbgMzse8BzhORTbD3hsQzPVlFcIiJSy2VyDegOYDuh77emhMYGb6I+2UREpBIySUCHAL9291fcfbO7LwIuBnYzs4KqDU9ERGqrTBJQAfCfpLKPACM0PhAREclaps2wk286FRERqZRMm2E/a2apHsX9fHK5u3eofFgiIlLbZZKAbqnyKEREpM5Jm4DcXQlIRERyLi8eSCciInWPEpCIiMQi277gROK15Jm4IxCRHNEZkIiIxEIJSEREYqEEJCIisVACEhGRWCgBiYhILJSAREQkFkpAIiISCyUgERGJhRKQiIjEQglIRERioQQkIiKxUAISEZFYKAGJiEgslIBERCQWSkAiIhILJSAREYmFEpCIiMRCCUhERGKhBCQiIrFQAhIRkVjkRQIysy5m9oKZvW9mC83sp1F5GzObaWb/jv63jjtWERHJjbxIQMBW4GfuvjdwMHC5me0N3AA87+57As9HwyIiUgvkRQJy95XuviB6vQFYBHQGTgIeiiZ7CDg5nghFRCTX8iIBJTKzbkBv4HWgo7uvjEatAjqWMc8IM5tvZvM/++yzaolTREQqJ68SkJk1B6YAV7n7+sRx7u6Ap5rP3ce6e6G7F7Zv374aIhURkcrKmwRkZg0JyecRd38iKl5tZgXR+AJgTVzxiYhIbuVFAjIzA+4HFrn7HxJGPQmcF70+D5hW3bGJiEjVaBB3AJF+wI+Bd83srajsl8BtwCQzuxD4BDgzpvikhmi6YWncIYhIhvIiAbn7y4CVMfqI6oxFRESqR15UwYmISN2jBCQiIrFQAhIRkVgoAYmISCyUgEREJBZKQCIiEgslIBERiYUSkIiIxEIJSEREYqEEJCIisVACEhGRWORFX3AiJZY8E3cEIlJNdAYkIiKxUAISEZFYqApOapS3ln0ZdwgikiM6AxIRkVgoAYmISCyUgEREJBZKQCIiEgslIBERiYUSkIiIxEIJSEREYqEEJCIisVACEhGRWCgBiYhILNQVj1Qv9XYtIhGdAYmISCyUgEREJBaqgpO8ot6uReoOJSCpW9Jdg9rr2OqJowzPvb+63PFH1l9Q/gJijl8kG0pAUq10hiMixZSAJLfSnGEs2PTvcsf/oMmeuYwme3l+hpRWTY9f6hQlIMkpneFUTrv/zip/gi6tqieQCpq9bHa54wd2GVgtcUjNoAQkpdXyX9DpEmSvqv6C131QIiWUgCQ7VfwFWtVVdOmW34s+5S8gzxNI2gS7VzUFIpKBvE9AZnYM8CegPnCfu98Wc0g1WyW/QB/4YF6546v6Gk66BNK6ksuvbBViujOodMuvdIIUqUHyOgGZWX3gHmAwsByYZ2ZPuvv78UZWe9X1azjpEkBay6o2AVd2+8x+7c5yxw88+GeVWr5INvI6AQEHAR+6+38AzGwCcBKQtwko3UXYVktWljv+y70Kyh2f9iJumjOcyv4CT6fSX+A1XNxViGltKn/0wEruP6TZf9N567nHyh2f7vNfcOJvKvX+Ur3yPQF1BpYlDC8H+iZPZGYjgBHR4EYzW1INsWWrHfB53EFkqSbGDDUz7ryI+cLsJs+LmBNdyE3pJsm7mDNQk2LeLZuJ8z0BZcTdxwJj446jPGY2390L444jGzUxZqiZcSvm6qGY80u+d0a6AuiSMLxrVCYiIjVcviegecCeZra7mTUChgJPxhyTiIjkQF5Xwbn7VjO7AniW0Az7AXdfGHNYFZXXVYRlqIkxQ82MWzFXD8WcR8zd445BRETqoHyvghMRkVpKCUhERGKhBFSNzGy0mS02s3fMbKqZ5XfXxoCZnWFmC81su5nldVNQMzvGzJaY2YdmdkPc8WTCzB4wszVm9l7csWTKzLqY2Qtm9n60b/w07pjSMbPGZvYvM3s7ivmWuGPKlJnVN7M3zWx63LHkmhJQ9ZoJ7Ovu+wMfAL+IOZ5MvAecCrwUdyDlSei26Vhgb+BsM9s73qgyMh44Ju4gsrQV+Jm77w0cDFxeA9b1FmCQux8A9AKOMbODY44pUz8FFsUdRFVQAqpG7j7D3bdGg68R7mvKa+6+yN3zsWeJZCXdNrn7N0Bxt015zd1fAtbGHUc23H2luy+IXm8gfDl2jjeq8nmwMRpsGC+s9HoAAA78SURBVP3lfQssM9sVOB64L+5YqoISUHwuAPK7b/+aJVW3TXn9pVgbmFk3oDfweryRpBdVZb0FrAFmunvexwz8EbgO2B53IFUhr+8DqonM7DmgU4pRv3L3adE0vyJUYzxSnbGVJZOYRZKZWXNgCnCVu6+PO5503H0b0Cu69jrVzPZ197y99mZmJwBr3P0NMxsYdzxVQQkox9z9yPLGm9lw4ATgCM+Tm7DSxVxDqNumamRmDQnJ5xF3fyLueLLh7l+a2QuEa295m4CAfsAQMzsOaAy0NLP/c/cfxRxXzqgKrhpFD9e7Dhji7kVxx1PLqNumamJmBtwPLHL3P8QdTybMrH1xq1Mza0J4xtjieKMqn7v/wt13dfduhP15Vm1KPqAEVN3uBloAM83sLTP7a9wBpWNmp5jZcuAQ4CkzezbumFKJGncUd9u0CJhUE7ptMrPHgLnAXma23MyyfCJCLPoBPwYGRfvxW9Gv9HxWALxgZu8QfqzMdPda16y5plFXPCIiEgudAYmISCyUgEREJBZKQCIiEgslIBERiYUSkIiIxKLGJyAzu9nMPPrbbmbrzGyemY0ys05J03aLpjshoayZmU0wsy+iccOj8p+Y2cdmttXMZlfvp6oZzGy8mc1PGD7IzG5OMd3NZvZ5tQaXR8xsvpmNTzPNdbm+293Mvh+t+1ZJ5cOjfb15BZbp0VOKcxHf5MRjq7r3k1yvczMrMLOnzeyraD0NjMp/Y2Yrou+n8Tl8vxFmdnKulpfhew6MPtu+uVhejU9Aka8I96kcSrhh6wnCfQrvmtmBCdOtjKZ7OaHsUuBEYATf3evSCbgXmAb8ELisqj9ADXUrMDxh+CDgpnhCqfGuAwbmeJnfJ2yPXD724xDg8RwuL9F9wNFVtOxUcr3OfwUcAJxNWE8LokeY3EK4B7Af4ZjJlRFAtSagXKstXfFsdffXEoafNbN7CY8QmGBmPdx9m7tvIfRCnagHsMTdpxQXmNlhQH3gAXd/pzKBmVkTd99UmWXkK3f/KO4YUonu1N/J3TfHHUuuxb0/JR1nuV72ckInsmWK+/On0QN43d2fLi4wsx7Ry3tqQn951c7da/QfcDPweRnjjiF0uX5MNNwtGj4hGl4aDSf+3ZyibHg0fWPgDkKvy1uAt4Hjkt5zKXAncCPhYPo2Kq8H3AB8GM37AXBe0ryzgcnAOdF06wk9Zu+aNF2TKI5PomV9DPxv0jQXAQuj8Z8A16VZj7cAHyQMNwO+BRYklLUj9Mo7OBoeD8yPXg9Psd5mJ24jQq/JrwFFwJtA/wy2b9doHWyKPufwaB3NTt4HgMMId7lvAX4cjTsTeDcqWwaMAhqk23+i+K9I2q6/B66Otus6wiMfWiXNty/wCrCZ0CPDEGA+ML6cz7g0xbobmBDHNYRekT8jPHKiJJ6k5RRvg+aEX/bJy1yaNN1+hGdUfU3olubUDLZH8nqZTWb7bBfg6Wg7LiXsnym3Y8Jw8Wc4mtCt0kbg/oT9YgLhURZFhB4w9sr0OClvnZfxudsAY4HV0bZ9FeibtF5KrWvC8VHWdi13edE09QnPDPsgin850X4UrfeU31NlxN8OeAj4Ilpfs4HCFPthuft4wjbZNxqelLgNk7blaqBhuftTuh0u3/+Sd9qkcY0JX6I3R8PdKJ2AegNPEb4oDo7+diVUuTnhoDoYaB9NP53QlfulwFGEKoOtQK+kjbgSeI7w5XNqVH4P4QC6DjgSuB3YVhxLwk61LNoZTyJUJ64Gnk6YxoAZwIZoWUcAw4BxCdP8PPrcowh9Xt0Q7cBXlLMej4w+c8doeDDhy2Ib0DIqOyX6vM2j4fF8l4DaRzuvJ6zLvRO2URHwDnA+4aFxrxG+UJuWE5MBb0Xr9Ozo/RdE6yj5i6sI+Ai4GDgc2CvaRk448Ir74dsC/DXd/kPqBPRptA8cR6j+2Aj8JWGaJoQOUN8mPMTvXOA/0T4zvpzP2Rv4krA/Fa+7lglxrAQmRp/huMQvi6TlDOe7BNQS+Fk0fEq0zN5J070L/E+0nv4BfENS4shgvcwms312QbT+zonWzbvRusokAS0nVF0NIlSzt4mW9SbhB8YJhGr1ZUCTTI6T8tZ5is+8UxT/f6JlHEOont8AdIqmOTiaZlbxuga6R3E7YZ88ONouaZcXLfO+aJv8lnA8nkXoYgrCQxcXEb6/iuNvX852exlYRTj+TiTUDm0A9shyHy/eJsUJ6GjCj9Ldk7b3x8Cdab+/qzI5VMdf8k6bYvxK4N7odTcSElBUNp7oS7SslRyVHRGV/TBp2peAx5M24kqgcULZHtFGOi9p3oeBeUkH81dA64Syq6L3LT6wjo6Gh5TxeVtGO81NSeUjox2wfhnzFZ/xnJ4w/RTgv3x3BvmHpHhLrTtCX2xexjZywhMpi8t6kXB2WkZMx0fT9Eko6xzFmfzF5cBJSfO/BryQVHYdIanuWt7+Q+oE9BGlz57+CKxKGL4sim3XhLJ+0bLGp9mPPyf6oZQijgUpypdSTgKKhk+IhruVMd0FCWVtCT8uLkkTZ6oElG6fPS4aTjxj2C16v+TtmCoBjUmK4VbCL/k2CWWtozguz+Q4KW+dp5juQkIi2DOhrEG0P4xOWheTy9smmS6PUJ3nwJXlxFXumXXCdMU1QT9MKGtG+AH4tyz38eJtUpyA6hGS1i0J0wwi6fuzrL/a0gihPJaj5RxJ+AJ/xcwaFP8BzwOFSdM+76WvPxxBSEBTU8zbK3qcdLF57r4uYfj96H/xw9UGAWvdvayeng8h7FyPJ73XLKAjZTyF1d2/Jvwq6x8VDSAk1zlJZXPKeN90viEcoMWKP1d5T4XtQ9j55yXEuQJ4I8W0TsID/qJ1+gN2vGA+kXDQHJJp4Ale8O+eaAvhM3SIHk0AoRHGGx6uZRTH+wrhDKgynk4/SYXMKH7h7l8Q4qzIU3rT7bMHAas94QFw7v4JqbdjKk8lDR9JqDpcn7B/b4iWV3wspjtOsnFktOyPE94P4EV2PPZztbzDo//jKxZyKQcRniv0YnFBdLxPJ1RbJ0q3j5fi7tuBB4Fh0bVXCEl3vmfwrKXa0gghJTNrTPhltzoHi2tHeGjbtynGbUsaTn6/doT63K/KWHYB3118/TJp3DfR/8bR/7aEM6zy4oRw/SeVLoQ68VTmEHo4bgT0JVx72AacaWYtCGcto8p57/JsiHZWANz9m2h/bVz2LHQi/EpL9hmhV/FE6zw8irtYO8Jjl5O3RfFwm0yCTpJq2xihSuXbKN5UyaayCSgX+28qqT5Pedsjm+WQsKzy1kvydkwl1fF0MKFKKtnz0f90x0k2it8v1bFfkYY4mSyvLfC156bhQgGp1/9qdjwO0u3jqTxIuOZ9uJnNA04Drs0ksFqdgAi/IhoQuruvrLWEOutMmj16inm3EqpjUj1aN5svqC8IO1RZ1kb/TyD1F9eScuadQ7gAeQRhx3uLkIB+T1iX9SndhL2qrSJcW0rWnnDhNlHyOv+ccMB0SCrvGP0vXk+bgUaJE5hZ66wjDVYRqk6SJceQreTPBiniJlRD5aNVpF4HHQjXGdNJdTw9SeomzRui/+mOk2ysJVR3XZpi3JYqWt4XQDMza5mDJLSS1Ou/I98dBxXm7kstPFV5OLA7oYbhsUzmrbVVcNHNd7cTWuY8l4NFPk/4JbfR3ecn/6WZdxbhy3vnVPMm/XLPJI42iTfTJplLOKh3KeO9NpQxH4QEZIRGC69EZyzvRsv7GbDY3VOdkRT7BkrOPHNhHtDJzA4qLjCzzsCBZc8SeHj88hvAGUmjziT8CCj+UbIcaBEtt9hRlYj3QDMrqcYys35kloCyPftYDvRMKkuOO/lMJC7zgI5m1re4wMy6EqpIK+J5YB9gYYr9e0nCNOUdJ5D5On+ecB330xTv924F40+3vFnR/2E5iP91QjXagOICM2tKuMaaqx+U9xPOfC4D/u7uyWdSKdWWM6AGZnZw9LoF4QvqUqAp4SJ3chVZRcwkNPWcaWa3E6q4WhKqpRq7+y/KmtHdl0QPn5tgZncQfv00JhxE33f3iyoQx6NmNpJw3aYAGODuF3t43PDNwJ/MbDfCdZx6hJsSD3f3U8qJc62ZvU+41vOLqGy7mb1C2FnHpYmt+AmTPzWzWcD6hC+Einia0KJskpn9gpAIbyKc2aU6k0x2E+GesAcJzUn3I/xqHpdwneaf0XIfMLM7Cb/gLqlgvA8CvybczHwzoVXcrYSzsXQWA8eb2T8JjUiWpPmxMBX4s5n9kvAFfxphf0pUvO4vNrMJQFEFvzArq3g7Pm5m1xN+5d9Cxasm/wD8CJhlZn8m1Ex0JNw0/rK7P0aa4yRaTqbr/GHCPjHbzH5PaL3WlnBtZZW7j8ky/rTLi74zxgJ3mlkHwnHcitBIaGhC/Eeb2dGEM6aPo2t5pbj7s2b2KjDRzG6Ipr2WsH+OzjL2svwd+AvhR0WZ34U7SNdKId//KH3fznZCHeZ8wrWKTknTdqOCreCi8p0IB86HhF8fqwhfYMcnTLOUpNZJUbkRWgcV35vzGeGi47CEaWazYyuaHWIh7Di/J/wKLr6/YVTSfD8inAFsIrTnfx24JoP1eW/0focmlF0flQ1LmrbUuos+4x2ElnPbSboPKMV7lWpRVUY8u0XreDPh2tUIwsXzvyftA2U1xT+LcBb3TbS+St0HFE1zbLRdighngT2TY0u1XUndwml/QpPkLYQEcDIZtFYi/Gh6jXBPjlP6PqAd1hHh+tYfon1wHfCnaN0kx/OzaL1tZcf7gJonLXOHz5hum2Wxz3blu2T/CaG5fKb3Ae3QmgrYhZDwV0freinwf8A+mR4nZa3zMj73ztE6XpawLz0B9EuzLspa15ksrz7wS0KCKp7mgYTx3yPU7nxF+vuA2hMS37poG7xIQuvSTPfxNNvk/wgt4uql+54p/tMTUaVGMbOdCQfk3e5+U9zxiAhELfk+ISTIGzOdr7ZUwUktZWaXEM6m/k34FXcN4Uz0gTjjEhGIWsweQLjBuC3wt2zmVwKSfLeZUAW4G+HU/1/AkR7uIxGReO1COCbXABd7wj1wmVAVnIiIxKLWNsMWEZH8pgQkIiKxUAISEZFYKAGJiEgslIBERCQW/x+VE6IirP8C2gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "       ind: 0.2 ($\\pm$ 0.5)\n",
      "  manifold: 0.4 ($\\pm$ 0.5)\n",
      "       asv: 0.2 ($\\pm$ 0.5)\n",
      "      flow: 0.0 ($\\pm$ 0.0)\n"
     ]
    }
   ],
   "source": [
    "errors = {\n",
    "    'ind': get_error(gt_indirect, theta_ind),\n",
    "    'manifold': get_error(gt_indirect, theta_manifold),\n",
    "    'asv': get_error(gt_indirect, theta_asv),\n",
    "    'flow': get_error(gt_indirect, theta_flow_indirect)\n",
    "}\n",
    "\n",
    "bins=np.histogram(np.hstack(list(errors.values())), bins=40)[1] #get the bin edges\n",
    "for k, v in errors.items():\n",
    "    plt.hist(v, bins, label=k, alpha=0.3, weights=np.zeros_like(v) + 100. / v.size)\n",
    "\n",
    "plt.legend(fontsize=15)\n",
    "plt.xlabel('Difference with ground truth indirect effect only', fontsize=15)\n",
    "plt.ylabel('Percentage', fontsize=15)\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "for k, v in errors.items():\n",
    "    print(f\"{k:>10s}: {np.abs(v).mean().round(1)} ($\\pm$ {np.abs(v).std().round(1)})\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Only Shapley Flow makes no mistake on both the direct and indirect effects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}