examples/LongitudinalLiNGAM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Longitudinal LiNGAM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import and settings\n",
    "In this example, we need to import `numpy`, `pandas`, and `graphviz` in addition to `lingam`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['1.24.4', '2.0.3', '0.20.1', '1.8.3']\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import graphviz\n",
    "import lingam\n",
    "from lingam.utils import print_causal_directions, print_dagc, make_dot\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "print([np.__version__, pd.__version__, graphviz.__version__, lingam.__version__])\n",
    "\n",
    "np.set_printoptions(precision=3, suppress=True)\n",
    "np.random.seed(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test data\n",
    "We create test data consisting of 5 variables. The causal model at each timepoint is as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# setting\n",
    "n_features = 5\n",
    "n_samples = 200\n",
    "n_lags = 1\n",
    "n_timepoints = 3\n",
    "\n",
    "causal_orders = []\n",
    "B_t_true = np.empty((n_timepoints, n_features, n_features))\n",
    "B_tau_true = np.empty((n_timepoints, n_lags, n_features, n_features))\n",
    "X_t = np.empty((n_timepoints, n_samples, n_features))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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.43.0 (0)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"156pt\" height=\"392pt\"\n",
       " viewBox=\"0.00 0.00 156.45 392.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 388)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-388 152.45,-388 152.45,4 -4,4\"/>\n",
       "<!-- x0(0) -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>x0(0)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"112.45\" cy=\"-105\" rx=\"34.39\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"112.45\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">x0(0)</text>\n",
       "</g>\n",
       "<!-- x4(0) -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>x4(0)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"39.45\" cy=\"-18\" rx=\"34.39\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"39.45\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x4(0)</text>\n",
       "</g>\n",
       "<!-- x0(0)&#45;&gt;x4(0) -->\n",
       "<g id=\"edge6\" class=\"edge\">\n",
       "<title>x0(0)&#45;&gt;x4(0)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M99.06,-88.41C88.01,-75.55 72.17,-57.11 59.6,-42.46\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"62.17,-40.09 53,-34.78 56.86,-44.65 62.17,-40.09\"/>\n",
       "<text text-anchor=\"middle\" x=\"95.45\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\">0.10</text>\n",
       "</g>\n",
       "<!-- x1(0) -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>x1(0)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"34.45\" cy=\"-192\" rx=\"34.39\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"34.45\" y=\"-188.3\" font-family=\"Times,serif\" font-size=\"14.00\">x1(0)</text>\n",
       "</g>\n",
       "<!-- x1(0)&#45;&gt;x0(0) -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>x1(0)&#45;&gt;x0(0)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M43.21,-174.21C49.01,-163.98 57.23,-151 66.45,-141 71.5,-135.52 77.49,-130.28 83.45,-125.59\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"85.57,-128.38 91.47,-119.58 81.37,-122.78 85.57,-128.38\"/>\n",
       "<text text-anchor=\"middle\" x=\"82.45\" y=\"-144.8\" font-family=\"Times,serif\" font-size=\"14.00\">0.50</text>\n",
       "</g>\n",
       "<!-- x1(0)&#45;&gt;x4(0) -->\n",
       "<g id=\"edge7\" class=\"edge\">\n",
       "<title>x1(0)&#45;&gt;x4(0)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M33.37,-173.69C32.26,-153.04 30.88,-117.48 32.45,-87 33.14,-73.52 34.6,-58.61 36,-46.18\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"39.51,-46.35 37.21,-36.01 32.56,-45.53 39.51,-46.35\"/>\n",
       "<text text-anchor=\"middle\" x=\"50.95\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">&#45;0.70</text>\n",
       "</g>\n",
       "<!-- x2(0) -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>x2(0)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"111.45\" cy=\"-279\" rx=\"34.39\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"111.45\" y=\"-275.3\" font-family=\"Times,serif\" font-size=\"14.00\">x2(0)</text>\n",
       "</g>\n",
       "<!-- x2(0)&#45;&gt;x0(0) -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>x2(0)&#45;&gt;x0(0)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M111.55,-260.88C111.72,-231 112.08,-169.11 112.29,-133.27\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"115.79,-133.07 112.35,-123.05 108.79,-133.03 115.79,-133.07\"/>\n",
       "<text text-anchor=\"middle\" x=\"129.95\" y=\"-188.3\" font-family=\"Times,serif\" font-size=\"14.00\">&#45;0.30</text>\n",
       "</g>\n",
       "<!-- x2(0)&#45;&gt;x1(0) -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>x2(0)&#45;&gt;x1(0)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M88.77,-265.33C79.52,-259.39 69.2,-251.7 61.45,-243 55.22,-236.01 49.9,-227.33 45.67,-219.18\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"48.73,-217.46 41.22,-209.98 42.43,-220.51 48.73,-217.46\"/>\n",
       "<text text-anchor=\"middle\" x=\"79.95\" y=\"-231.8\" font-family=\"Times,serif\" font-size=\"14.00\">&#45;0.30</text>\n",
       "</g>\n",
       "<!-- x3(0) -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>x3(0)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"43.45\" cy=\"-366\" rx=\"34.39\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"43.45\" y=\"-362.3\" font-family=\"Times,serif\" font-size=\"14.00\">x3(0)</text>\n",
       "</g>\n",
       "<!-- x3(0)&#45;&gt;x1(0) -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>x3(0)&#45;&gt;x1(0)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M41.3,-347.99C39.66,-334.25 37.51,-314.41 36.45,-297 34.87,-271.06 34.43,-241.42 34.35,-220.43\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"37.85,-220.31 34.34,-210.32 30.85,-220.32 37.85,-220.31\"/>\n",
       "<text text-anchor=\"middle\" x=\"52.45\" y=\"-275.3\" font-family=\"Times,serif\" font-size=\"14.00\">0.40</text>\n",
       "</g>\n",
       "<!-- x3(0)&#45;&gt;x2(0) -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>x3(0)&#45;&gt;x2(0)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M56.24,-349.01C66.4,-336.31 80.77,-318.34 92.31,-303.93\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"95.2,-305.91 98.71,-295.92 89.73,-301.54 95.2,-305.91\"/>\n",
       "<text text-anchor=\"middle\" x=\"97.45\" y=\"-318.8\" font-family=\"Times,serif\" font-size=\"14.00\">0.30</text>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.graphs.Digraph at 0x7f341f313a00>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# B(0,0)\n",
    "B_t_true[0] = np.array([[0.0, 0.5,-0.3, 0.0, 0.0],\n",
    "                        [0.0, 0.0,-0.3, 0.4, 0.0],\n",
    "                        [0.0, 0.0, 0.0, 0.3, 0.0],\n",
    "                        [0.0, 0.0, 0.0, 0.0, 0.0],\n",
    "                        [0.1,-0.7, 0.0, 0.0, 0.0]])\n",
    "causal_orders.append([3, 2, 1, 0, 4])\n",
    "make_dot(B_t_true[0], labels=[f'x{i}(0)' for i in range(5)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.43.0 (0)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"209pt\" height=\"305pt\"\n",
       " viewBox=\"0.00 0.00 209.31 305.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 301)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-301 205.31,-301 205.31,4 -4,4\"/>\n",
       "<!-- x0(1) -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>x0(1)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"79.86\" cy=\"-18\" rx=\"34.39\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"79.86\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x0(1)</text>\n",
       "</g>\n",
       "<!-- x1(1) -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>x1(1)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"79.86\" cy=\"-192\" rx=\"34.39\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"79.86\" y=\"-188.3\" font-family=\"Times,serif\" font-size=\"14.00\">x1(1)</text>\n",
       "</g>\n",
       "<!-- x1(1)&#45;&gt;x0(1) -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>x1(1)&#45;&gt;x0(1)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M57.43,-178.08C39.46,-166.31 15.54,-147.1 4.86,-123 -1.62,-108.37 -1.62,-101.63 4.86,-87 13.96,-66.47 32.65,-49.5 49.1,-37.65\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"51.17,-40.47 57.43,-31.92 47.2,-34.7 51.17,-40.47\"/>\n",
       "<text text-anchor=\"middle\" x=\"20.86\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">0.20</text>\n",
       "</g>\n",
       "<!-- x2(1) -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>x2(1)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"79.86\" cy=\"-105\" rx=\"34.39\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"79.86\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">x2(1)</text>\n",
       "</g>\n",
       "<!-- x1(1)&#45;&gt;x2(1) -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>x1(1)&#45;&gt;x2(1)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M79.86,-173.8C79.86,-162.16 79.86,-146.55 79.86,-133.24\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"83.36,-133.18 79.86,-123.18 76.36,-133.18 83.36,-133.18\"/>\n",
       "<text text-anchor=\"middle\" x=\"95.86\" y=\"-144.8\" font-family=\"Times,serif\" font-size=\"14.00\">0.30</text>\n",
       "</g>\n",
       "<!-- x4(1) -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>x4(1)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"166.86\" cy=\"-105\" rx=\"34.39\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"166.86\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">x4(1)</text>\n",
       "</g>\n",
       "<!-- x1(1)&#45;&gt;x4(1) -->\n",
       "<g id=\"edge6\" class=\"edge\">\n",
       "<title>x1(1)&#45;&gt;x4(1)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M95.65,-175.68C101.94,-169.56 109.25,-162.45 115.86,-156 125.06,-147.02 135.17,-137.12 143.94,-128.51\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"146.57,-130.84 151.25,-121.34 141.67,-125.84 146.57,-130.84\"/>\n",
       "<text text-anchor=\"middle\" x=\"149.36\" y=\"-144.8\" font-family=\"Times,serif\" font-size=\"14.00\">&#45;0.40</text>\n",
       "</g>\n",
       "<!-- x2(1)&#45;&gt;x0(1) -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>x2(1)&#45;&gt;x0(1)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M79.86,-86.8C79.86,-75.16 79.86,-59.55 79.86,-46.24\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"83.36,-46.18 79.86,-36.18 76.36,-46.18 83.36,-46.18\"/>\n",
       "<text text-anchor=\"middle\" x=\"98.36\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\">&#45;0.10</text>\n",
       "</g>\n",
       "<!-- x3(1) -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>x3(1)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"79.86\" cy=\"-279\" rx=\"34.39\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"79.86\" y=\"-275.3\" font-family=\"Times,serif\" font-size=\"14.00\">x3(1)</text>\n",
       "</g>\n",
       "<!-- x3(1)&#45;&gt;x1(1) -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>x3(1)&#45;&gt;x1(1)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M79.86,-260.8C79.86,-249.16 79.86,-233.55 79.86,-220.24\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"83.36,-220.18 79.86,-210.18 76.36,-220.18 83.36,-220.18\"/>\n",
       "<text text-anchor=\"middle\" x=\"95.86\" y=\"-231.8\" font-family=\"Times,serif\" font-size=\"14.00\">0.40</text>\n",
       "</g>\n",
       "<!-- x4(1)&#45;&gt;x0(1) -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>x4(1)&#45;&gt;x0(1)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M153.23,-88.23C144.31,-78.15 132.24,-64.97 120.86,-54 116.09,-49.4 110.81,-44.66 105.68,-40.22\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"107.71,-37.35 97.83,-33.54 103.18,-42.68 107.71,-37.35\"/>\n",
       "<text text-anchor=\"middle\" x=\"154.36\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\">&#45;0.50</text>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.graphs.Digraph at 0x7f3478092ca0>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# B(1,1)\n",
    "B_t_true[1] = np.array([[0.0, 0.2,-0.1, 0.0,-0.5],\n",
    "                        [0.0, 0.0, 0.0, 0.4, 0.0],\n",
    "                        [0.0, 0.3, 0.0, 0.0, 0.0],\n",
    "                        [0.0, 0.0, 0.0, 0.0, 0.0],\n",
    "                        [0.0,-0.4, 0.0, 0.0, 0.0]])\n",
    "causal_orders.append([3, 1, 2, 4, 0])\n",
    "make_dot(B_t_true[1], labels=[f'x{i}(1)' for i in range(5)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.43.0 (0)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"197pt\" height=\"218pt\"\n",
       " viewBox=\"0.00 0.00 196.89 218.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 214)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-214 192.89,-214 192.89,4 -4,4\"/>\n",
       "<!-- x0(2) -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>x0(2)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"110.45\" cy=\"-192\" rx=\"34.39\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"110.45\" y=\"-188.3\" font-family=\"Times,serif\" font-size=\"14.00\">x0(2)</text>\n",
       "</g>\n",
       "<!-- x2(2) -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>x2(2)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"67.45\" cy=\"-105\" rx=\"34.39\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"67.45\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">x2(2)</text>\n",
       "</g>\n",
       "<!-- x0(2)&#45;&gt;x2(2) -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>x0(2)&#45;&gt;x2(2)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M101.8,-174.29C98.82,-168.5 95.48,-161.97 92.45,-156 88.45,-148.12 84.14,-139.53 80.24,-131.72\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"83.33,-130.08 75.73,-122.69 77.06,-133.2 83.33,-130.08\"/>\n",
       "<text text-anchor=\"middle\" x=\"108.45\" y=\"-144.8\" font-family=\"Times,serif\" font-size=\"14.00\">0.20</text>\n",
       "</g>\n",
       "<!-- x4(2) -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>x4(2)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"154.45\" cy=\"-105\" rx=\"34.39\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"154.45\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">x4(2)</text>\n",
       "</g>\n",
       "<!-- x0(2)&#45;&gt;x4(2) -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>x0(2)&#45;&gt;x4(2)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M119.14,-174.21C125.45,-162.02 134.09,-145.32 141.25,-131.49\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"144.43,-132.96 145.92,-122.47 138.22,-129.74 144.43,-132.96\"/>\n",
       "<text text-anchor=\"middle\" x=\"150.45\" y=\"-144.8\" font-family=\"Times,serif\" font-size=\"14.00\">0.30</text>\n",
       "</g>\n",
       "<!-- x1(2) -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>x1(2)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"137.45\" cy=\"-18\" rx=\"34.39\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"137.45\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x1(2)</text>\n",
       "</g>\n",
       "<!-- x2(2)&#45;&gt;x1(2) -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>x2(2)&#45;&gt;x1(2)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M78.22,-87.79C85.09,-77.78 94.38,-64.82 103.45,-54 106.9,-49.88 110.72,-45.66 114.52,-41.64\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"117.26,-43.84 121.71,-34.22 112.23,-38.97 117.26,-43.84\"/>\n",
       "<text text-anchor=\"middle\" x=\"121.95\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\">&#45;0.70</text>\n",
       "</g>\n",
       "<!-- x3(2) -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>x3(2)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"34.45\" cy=\"-18\" rx=\"34.39\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"34.45\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x3(2)</text>\n",
       "</g>\n",
       "<!-- x2(2)&#45;&gt;x3(2) -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>x2(2)&#45;&gt;x3(2)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M49.69,-89.37C44.05,-83.66 38.5,-76.64 35.45,-69 32.64,-61.99 31.63,-53.96 31.51,-46.43\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"35.02,-46.37 31.85,-36.26 28.02,-46.14 35.02,-46.37\"/>\n",
       "<text text-anchor=\"middle\" x=\"53.95\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\">&#45;0.40</text>\n",
       "</g>\n",
       "<!-- x4(2)&#45;&gt;x1(2) -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>x4(2)&#45;&gt;x1(2)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M151.01,-86.8C148.68,-75.16 145.56,-59.55 142.89,-46.24\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"146.28,-45.29 140.88,-36.18 139.41,-46.67 146.28,-45.29\"/>\n",
       "<text text-anchor=\"middle\" x=\"162.45\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\">0.50</text>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.graphs.Digraph at 0x7f341f347400>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# B(2,2)\n",
    "B_t_true[2] = np.array([[0.0, 0.0, 0.0, 0.0, 0.0],\n",
    "                        [0.0, 0.0,-0.7, 0.0, 0.5],\n",
    "                        [0.2, 0.0, 0.0, 0.0, 0.0],\n",
    "                        [0.0, 0.0,-0.4, 0.0, 0.0],\n",
    "                        [0.3, 0.0, 0.0, 0.0, 0.0]])\n",
    "causal_orders.append([0, 2, 4, 3, 1])\n",
    "make_dot(B_t_true[2], labels=[f'x{i}(2)' for i in range(5)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create B(t,t-τ) and X\n",
    "for t in range(n_timepoints):\n",
    "    # external influence\n",
    "    expon = 0.1\n",
    "    ext = np.empty((n_features, n_samples))\n",
    "    for i in range(n_features):\n",
    "        ext[i, :] = np.random.normal(size=(1, n_samples));\n",
    "        ext[i, :] = np.multiply(np.sign(ext[i, :]), abs(ext[i, :]) ** expon);\n",
    "        ext[i, :] = ext[i, :] - np.mean(ext[i, :]);\n",
    "        ext[i, :] = ext[i, :] / np.std(ext[i, :]);\n",
    "\n",
    "    # create B(t,t-τ)\n",
    "    for tau in range(n_lags):\n",
    "        value = np.random.uniform(low=0.01, high=0.5, size=(n_features, n_features))\n",
    "        sign = np.random.choice([-1, 1], size=(n_features, n_features))\n",
    "        B_tau_true[t, tau] = np.multiply(value, sign)\n",
    "\n",
    "    # create X(t)\n",
    "    X = np.zeros((n_features, n_samples))\n",
    "    for co in causal_orders[t]:\n",
    "        X[co] = np.dot(B_t_true[t][co, :], X) + ext[co]\n",
    "        if t > 0:\n",
    "            for tau in range(n_lags):\n",
    "                X[co] = X[co] + np.dot(B_tau_true[t, tau][co, :], X_t[t-(tau+1)].T)\n",
    "    \n",
    "    X_t[t] = X.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Causal Discovery\n",
    "To run causal discovery, we create a `LongitudinalLiNGAM` object by specifying the `n_lags` parameter. Then, we call the `fit` method. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = lingam.LongitudinalLiNGAM(n_lags=n_lags)\n",
    "model = model.fit(X_t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the `causal_orders_` property, we can see the causal ordering in time-points as a result of the causal discovery. \n",
    "All elements are nan because the causal order of B(t,t) at t=0 is not calculated.\n",
    "So access to the time points above t=1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[nan, nan, nan, nan, nan]\n",
      "[3, 1, 2, 4, 0]\n",
      "[0, 4, 2, 3, 1]\n"
     ]
    }
   ],
   "source": [
    "print(model.causal_orders_[0]) # nan at t=0\n",
    "print(model.causal_orders_[1])\n",
    "print(model.causal_orders_[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also, using the `adjacency_matrices_` property, we can see the adjacency matrix as a result of the causal discovery.\n",
    "As with the causal order, all elements are nan because the B(t,t) and B(t,t-τ) at t=0 is not calculated. \n",
    "So access to the time points above t=1.\n",
    "Also, if we run causal discovery with n_lags=2, B(t,t-τ) at t=1 is also not computed, so all the elements are nan."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "B(0,0):\n",
      "[[nan nan nan nan nan]\n",
      " [nan nan nan nan nan]\n",
      " [nan nan nan nan nan]\n",
      " [nan nan nan nan nan]\n",
      " [nan nan nan nan nan]]\n",
      "B(0,-1):\n",
      "[[nan nan nan nan nan]\n",
      " [nan nan nan nan nan]\n",
      " [nan nan nan nan nan]\n",
      " [nan nan nan nan nan]\n",
      " [nan nan nan nan nan]]\n",
      "B(1,1):\n",
      "[[ 0.     0.     0.     0.    -0.611]\n",
      " [ 0.     0.     0.     0.398  0.   ]\n",
      " [ 0.     0.328  0.     0.     0.   ]\n",
      " [ 0.     0.     0.     0.     0.   ]\n",
      " [ 0.    -0.338  0.     0.     0.   ]]\n",
      "B(1,0):\n",
      "[[ 0.029  0.064 -0.27   0.065 -0.18 ]\n",
      " [ 0.139 -0.211 -0.43   0.558  0.051]\n",
      " [-0.181  0.178  0.466  0.214  0.079]\n",
      " [ 0.384 -0.083 -0.495 -0.072 -0.323]\n",
      " [-0.174 -0.383 -0.274 -0.275  0.457]]\n",
      "B(2,2):\n",
      "[[ 0.     0.     0.     0.     0.   ]\n",
      " [ 0.     0.    -0.696  0.     0.487]\n",
      " [ 0.231  0.     0.     0.     0.   ]\n",
      " [ 0.     0.    -0.409  0.     0.   ]\n",
      " [ 0.25   0.     0.     0.     0.   ]]\n",
      "B(2,1):\n",
      "[[ 0.194  0.2    0.031 -0.473 -0.002]\n",
      " [-0.376 -0.038  0.16   0.261  0.102]\n",
      " [ 0.117  0.266 -0.05   0.523 -0.019]\n",
      " [ 0.249 -0.448  0.473 -0.001 -0.177]\n",
      " [-0.177  0.309 -0.112  0.295 -0.273]]\n"
     ]
    }
   ],
   "source": [
    "t = 0 # nan at t=0\n",
    "print('B(0,0):')\n",
    "print(model.adjacency_matrices_[t, 0])\n",
    "print('B(0,-1):')\n",
    "print(model.adjacency_matrices_[t, 1])\n",
    "\n",
    "t = 1\n",
    "print('B(1,1):')\n",
    "print(model.adjacency_matrices_[t, 0])\n",
    "print('B(1,0):')\n",
    "print(model.adjacency_matrices_[t, 1])\n",
    "\n",
    "t = 2\n",
    "print('B(2,2):')\n",
    "print(model.adjacency_matrices_[t, 0])\n",
    "print('B(2,1):')\n",
    "print(model.adjacency_matrices_[t, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAADQCAYAAADI+yJFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8WgzjOAAAACXBIWXMAAAsTAAALEwEAmpwYAAA36klEQVR4nO3dd3gUVRfA4d+FEDqh9yogRaVIaIYmAQREqgUCiHSxgGLDgooIIoL6KSCgFAmIJKL0CCqE0ExoARRRug1CT0IIhCT3++MuMYFAQrK7s7s57/Ps4+7O7MwZ5HDuzNy5V2mtEUIIIYRrymV1AEIIIYS4OSnUQgghhAuTQi2EEEK4MCnUQgghhAuTQi2EEEK4MCnUQgghhAuTQi1SKKU6K6Wa297fqZQao5Q6qpR68Ra/yaOUaq+UmqOUildKlUy1rKdSqoETQhdC3EJWcjuD7UluO5EU6hxEKXVMKRWplNqslNqhlFqmlKpmW1YBeBmIsK3+NJAI/J7BZqsA/YEfgXzXLfsB+EwpVcBuByGEuIGDchul1GCl1G7bK8S2LZDcdiop1DnPEK11C621LxAGzLB9PxaYqbVOAtBaj9JaTwFO3mpjWutDWuvHtdaL01kWC6wEnrHrEQgh0mPX3FZKPQi8ArTTWjcEVgFf27Yhue1EUqhzKKVUHqA8kEsplRvojkk8e/sO6OOA7Qoh0mHH3B4KfKK1Pmv7PBOooZSqa/ssue0kUqhznllKqW3AeaA+MBKoDpzVWsfZe2da69+AO5VS3vbethAiDXvndi3g12sfbGfkv9u+l9x2IinUOc9wrXVzoAzmPtPDQElMcjvKBaCEA7cvhLB/bl8Frp8M4vrPF5DcdjgvqwMQ1tBaxymlpgAJmPtZBR24u0JAjAO3L4SwsWNu/w7UBUIBlFIKqE3aTmiS204gZ9Q520OYFvFeoKLtfpZdKaXKADGOuKwuhLgpe+T2bGCkUqq47fNA4IjWej9IbjuTnFHnPF8opS4C+YFYoLPWOloptR1oBmzJaANKqcXARq31zEzsrwOwIjsBCyEyxa65rbX+QSk1DdiklIoDzgG9U60uue0kSuajFgBKKT/gWa1170ysexDoqLU+nMF6CnPprZ/W+rh9IhVC3A7Jbfcnl74FAFrrLcBhW1LflFKqHvBHRols0xf4UhJZCOtIbrs/OaMWKWytZHQGfymUUnm01lczsb1cWutke8UnhMgayW335vQzaqXUw0qpIKXUnzdZ/qhSKkIptVMpNdXZ8eVk2iYT62WYyLb1JJFzEMlt1yW57d6suPR9GngKuOEheaVUFWA80B7wxfRW7OXc8IQQWSS5LYQDOL3Xt9Z6I4DtSsz1OgJLtdbRtnVmYR4JWHr9ikqpYcAwgIIFCzaqXbu2o0IWwiliY+HwYUhK2nlGa13K6nhul+S2EOn76y84dQoga7ntao9nlSDtQPEngNLprai1no15zg9fX1+9Y8cOx0cnhIPMmQNPPgk1a8KBA8oTO+hIboscJyYGeveGnTth5Ej45JOs5bar9fqOIm3ylrV9J4RHSkqCF16AIUPA3x9+/tnqiBxGclvkKEeOQPPm8MMPMHMm/O9/Wd+WqxXqNUAPpVRh2+dBwHIL4xHCYWJioGtX+PBD09petQp8fKyOymEkt0WOsWkTNG0KJ07A2rUwfHj2tucShVop9bVSqoHW+gQwEQhTSoUDUVrrG+5hCeHujh6F++4zSfzZZ6a17eVqN6LsQHJb5DTz5pmrY8WLQ3g4tG2b/W1a9k+D1rpsqve9U71fBCyyJCghnGDzZujRw1z2XrvWJLUnkdwWOVFSEowZA1OmQLt2EBQExYrZZ9sucUYtRE4xf75pYRcvbu5He1qRFiInio01je8pU+Dpp2HNGvsVaZBCLYRTJCXByy/DwIHQqpUp0nfeaXVUQojsOnYM/PxMcZ42zbzy5LHvPjzwrpgQriU2Fvr1gxUrYMQIcz/a3okshHC+LVvMmXRCAoSEQPv2jtmPnFEL4UDHj0OLFrB6tWlpz5ghRVoITxAYaG5j+fiYK2SOKtIghVoIh9m2DZo0McV6zRpz70oI4d6Sk+G11+Dxx80l7/BwcPTgeVKohXCAhQuhTRsoXNi0tjt0sDoiIUR2XbwIvXrBe+/BsGHmqY3ixR2/XynUQtjRtdZ2//7mOWlntLaFEI7355/mNtaKFaafycyZzruNJZ3JhLCTuDhToL/7DoYONfekvW+YR0oI4W5+/hm6d4f4eNPfpGNH5+5fzqiFsIO//jKt7eXL4eOPYdYsKdJCeILFi81trIIFTb8TZxdpkEItRLaFh5tOY4cPm/G6R42C9Gd6FEK4i+RkGDsWAgLMuN3h4VC3rjWxSKEWIhsWL4bWrSF/ftPa7tTJ6oiEENkVFwePPgrvvguDBpkZsEqWtC4eKdRCZEFyMrz1lmltN2kCERFw111WRyWEyK5//jGjB377LUydCl98Yf1tLOlMJsRtunQJnngCgoPNkKAzZ1qfyEKI7Nuxw0w9Gxtrend36WJ1RIacUQtxG661tr/5xgzAP2eOFGkhPEFQELRsCXnzmttYrlKkQQq1EJm2Ywc0bgy//25a2y+8IJ3GhHB3WsO4cfDYY9Cokek0dvfdVkeVlhRqITIhONicSXt7w9atrtXaFkJkTXw89OkDb79thgT96ScoXdrqqG4khVqIW9Aaxo83PUAbNjSdxu65x+qohBDZdeKEeWIjKAjef9/MFZ83r9VRpU86kwlxE/Hx5tGMr782re3Zs103kYUQmbdrl+k0duGCGUmwWzerI7o1OaMWIh3XWttLlsCkSa7d2hZCZN7SpWYUwVy5zHzSrl6kQQq1EDfYvds8G71/v3mW8pVXpNOYEO5Oa5gwAR5+GOrXN7ex6te3OqrMkUItRCrffWda20rB5s1mIH4hhHu7fNlMmPPGG9C3L2zYAGXLWh1V5kmhFgLT2p44EXr2NJ3FIiKgQQOroxJCZNfJk3D//bBokTmjDgyEfPmsjur2SGcykeNdvgxDhphEDggwg5i4WyILIW60Zw889BCcPWvuTffsaXVEWSNn1CJHi4qCtm1NkX73XVi4UIq0EJ5g+XLw8zPj8m/a5L5FGqRQixxs717TaSwy0gwJ+vrr0mlMCFeTmJjIhx9+SP369alevTrPPPMMJ0+evOn6Wpvnonv0MNNSbt8O997rxIAdQAq1yJFWrID77oOkJNNprFcvqyMSQqRn6NChrFq1iunTp7Ny5Uq8vb1p0aIF0dHRN6x75YqZMGfMGDNI0caNUK6c82O2N7lHLXIUreGDD0wiN2pkLo+VL291VEKI9Bw8eJA1a9Zw9OhRChQoAMCHH37IP//8w5dffsnIkSNT1j11ylze3rLFjN09dqznXCGTQi1yjCtX4MknzeAljz4K8+aBLfeFEC5o165dtGzZMqVIX9OxY0c2bNiQ8nnfPtNp7NQpMyToI484O1LHkkvfIkc4fRratTNF+u23zbCgUqSFcG1VqlRh3759aK3TfL93716qVq0KwKpV5jbW1asQFuZ5RRosKNRKqUeVUhFKqZ1KqanpLA+97tXE2TEKz/LLL6bT2I4dpkC/9ZbnXBJzJZLbwt6aNm1KiRIlePHFF4mNjSU5OZng4GAWLVrE4MFDmDrVjNl9551m7ANfX6sjdgynXvpWSlUBxgNNgBjga6VUL6310lSr5dVaN3dmXMJzrV4NvXtD4cKmtd24sdUReSbJbeEISimWL1/OU089Rfny5fHy8uKOO+4gOHg577xTmblzzZCgX37p2VfInH1G3RFYqrWO1uZaxiyg+7WFSikvoKhSKkgpFaaUGq+Uyu3kGIUH0Bo+/NDct7rzTvOIhhRph5LcFna1du1aWrduTd26dTlx4gRz5szhwIEDfPDBPJ57ri5z55oOY0uWeHaRBucX6hJA6gfgTgCpp+kuBIQCw4A2QDlgSHobUkoNU0rtUErtOH36tEOCFe4pIQGGDoUXXjC9QMPCoEIFq6PyeJLbwm5Wr17NwIEDefbZZ/nll1948cUXGTlyJA0b9uWBB4qyZ08+ihR5kvr1l5IrB/S0cvYhRpE2ecvavgNAa31Baz3C9t9k4FvMpbQbaK1na619tda+pUqVcmjQwn2cOQPt25thQN94w/QALVjQ6qhyBMltYTfvvPMOs2bN4uGHH6ZMmTJ06dKFhAR/Tp78lhIlKrFtW15++mkII0aM4MCBA1aH63DOLtRrgB5KqcK2z4OA5dcWKqXKKqVeUyqlq09HYJeTYxRuav9+aNoUwsPNkKDjx5MjWtsuQnJb2E1kZCT+/v6AuY317LOHOH9+AVofITxc07Qp+Pr6MnjwYObPn29tsE7g1H/GtNYngIlAmFIqHIjSWi+19QC91gIvBOxSSm0CFDDbmTEK9/T999C8OcTFmdGIAgKsjihnkdwW9lSjRg127tzJ1atm7IOZM2tTsuQ2qlbtT5Uq/5WtKlWqcObMGQsjdQ6nD3iitV4ELLruuzapPr5mewlxA601gYGBLFiwgIsXL/LAAx0pUOBlXnutAPfcY4YGrVzZ6ihzJsltYS8vvfQSgwe/hI/PD+zYUZjevY8QFHQ/L788MWWd5ORkFi9ezLBhwyyM1DlkZDLhVkaPHk1YWBhjx46lSJESPPnkVQ4eLECXLoksXuxFoUJWRyiEyI7Q0FCCg/dx7NjXXL3qjbf3EMLCQmjXri3z5s2jVKlSFC1alNmzzQWZRzxxhJPrSKEWbuP48eMsWLCAw4cPk5xclEcegYMHoXr1YDp3PkehQsOtDlEIkQ2ff/45r78eysWLcylQAPz8pnDgwAa2bdtNqVKlWLlyJQsWLCAuLo4HH3yQwYMH4+3tbXXYDieFWriN8PBw2rRpQ1RUUR56CI4fhwULICEhhrCwrYwYIYVaCHcVHx/P888f5PLlhdStq1i5EqpUeZ0hQ44yffp0xo0bR9euXenatavVoTqd9IkVbqNcuXLs2lWCZs00Fy7A+vXQv7+ZYaecJ8xlJ0QOlZgI/fvHEBc3mc6dFVu2QJUqZlnv3r3ZuHGjtQFaTM6ohUv5/fffWbhwIbGxsXTu3Jn27dtz7YmevXv9OHasOWXKnGbTpmLUrJmHDRs2MGfOHLZu3Wpx5EKIrDh/3sxm9+OPZShQYBpBQcPJly9PyvJjx46R05+nlzNq4TIWLFhAy5YtuXLlCuXKlWP06NH07duXhIRknnkGnnkmF/7+V6lZ8wn8/CpQvXp1Bg8ezIIFC6hZs6bV4QshbtPBg9CsmXmkcu5c8PNbwVtvvcHVq1cBOHToEBMmTMgRPbtvRc6ohUuIjo5m1KhRbN26lTp16gAwatQomjTpQJMmZ9izpzQvvQTvvZef3LnX8Oeff3Lx4kVq1apF7twyZLQQ7mb9ejOhRu7c8NNP0LIlPPjgQgICAqhcuTIVK1bkyJEjvPPOO7Rv397qcC0lhVq4hPXr19OsWbOUIg3w11/5OHVqOadPF2bOHBg06L/1K8vD0kK4rVmz4JlnoFYtWLkSqlUz35cuXZoff/yRQ4cOcerUKerXr09BGQNYLn0L15A3b17i4uJSPm/YYIYDjY3NS5cu/0tTpIUQ7ikxEUaNMqONdegAW7f+V6RTq1GjBvfdd58UaRsp1MIl+Pv788cff/DDDz8we7ZJ4lKlEila9AFefDHduRuEEG4kOhq6dIFPPoHnnzejCBYpYnVU7kEufQuXkDdvXr76KoguXX4jPh7Kl9/D6dPdeOGFobRs2dLq8IQQ2XD4sJkb/uBB+PxzGJLuBKfiZqRQC5cQHQ0ffNCK+PhWdOr0Oz167KRjx01UqlTJ6tCEENmwcaOZFx7ghx+gTRtLw3FLUqiF5VK3tmfNgmHDagG1rA5LCJFNc+bAk09qatSAlSsVNWpYHZF7umWhVkotBnR6y7TWMpGgyLawMNPa1hrWrYP777c6opxBcls4UlISjBp1henT85I7908cPNibUaOaMnXqVGrXrm11eG4no85kM4FZmLlk99je/wkccHBcIgeYOxfatYNSpSA8XIq0k0luC4eIiYGuXTXTp+elbt2fOH78HmJi/qRDhw60bds2R8wfbW+3LNRa641a641Ada31+7bPY4A2TolOeKSkJHjxRRg82BTnbduQS2JOJrktrklMTCQsLIwNGzZw5cqVbG3r6FG47z5YuxZKlXqTffvup0KFMhQoUIBRo0bxwAMPMH/+fPsEnoNk9h51BaVUQa11nFKqCFDUgTEJDxYTAwEBsHq1GfDgo4/AS3pKWElyOwfbtGkTffr0oWzZsnh5eXHs2DHmzZtHp06dbntbmzdDjx7mWekXXviBf/45Qq5cac8FmzVrxo4dO+wVfo6R2X8ipwCRSqlfgbuAVxwXkvBUx46ZTmO//QYzZsCIEVZHJJDczrGio6Pp2bMngYGBdOzYEYAtW7bQrVs39u3bd1sz0s2fD8OGQdWqsGoVxMaW4OGHt5CYmIhXqpb4xo0badq0qZ2PxPNlasATrfVioAnwLtBYa/2tQ6MSHmfLFmjSBP7+G77/Xoq0q5Dczrm+/fZbWrVqlVKkAfz8/OjevTuLFy/O1DaSkuCVV2DgQGjVCrZuTeLff0P5+++/qVGjBn379uWPP/4gKiqKcePGsWXLFgYMGOCoQ/JYmTqjVkoVA94EKgBTlFKFtNbrHRqZ8Bhffmla21WqmHF9a8mTVy5DcjvnOn/+POXLl7/h+/Lly3P+/PkMfx8bC/36mRHGRoyAoUN/oVmz7vj4+FC6dGkiIiLw9fWldevWXLp0iQcffJCNGzdStGhRBxyNh9NaZ/gClgNdgVDAG9iYmd8569WoUSMtXE9SktavvKI1aN22rdZnz1odkXsAdmgn5Y7kds4VGRmpK1SooGNiYlK+i4+P17Vq1dKhoaG3/O2xY1rXq6d1rlxaf/qp1omJibpGjRp6/vz5KescP35cV6pUSW/atMlhx+BusprbmR3ru4jWegWQrLVO4CbPXwpxzcWL5vno9983A/B//z0UL251VCIdkts5VP369enWrRt+fn7MmTOHL7/8klatWnHvvffSqlWrm/5u2zZzG+v4cQgJMZ1Ct27dSsGCBdNc1q5cuTLPPfccX375pTMOx6NltjNZnFKqH5BbKeUPRDswJuHm/vwTunaFffvg00/h6adBKaujEjchuZ2DTZs2jWXLlhEUFERSUhKjR4/mkUceQd0kYRcuNI9VVqpkhga9NnZJTEwMJUuWvGH9UqVKER4e7shDyBEyW6gHY3qHFgcGAkMdFpFwa9u2mUc04uNhzRp44AGrIxIZkNzOwZRS9OjRgx49etyw7M8//2TChAmEhoZSrFgJSpb8jNWr69OmDXzzDZQo8d+6LVq0oF+/fhw+fJjq1asDkJyczLx58+jfv7+TjsZzZbZQ19dap/xpK6VaAqccE5JwV4sWmdZ2hQpmPuk6dayOSGSC5La4wcmTJ2nUqBFFixYlf/6S/P33x4SH16dBg+2sXdsYb+//1j137hyfffYZFStWpF69egwcOJBGjRqxYMEClFIEBMiItNmV2XvUY669UUp5AR84JhzhjpKT4Y03TA/QZs0gIkKKtBuR3M5hkpKSCAwMpGvXrnTp0oW5c+eSmJiYZp1+/fpx9epVnn56EvHx6/j3X1+KFXuHI0faExt7NmW906dP07RpUw4ePMiECRMYMmQI8+bNY86cOQwcOJDvv/+evHnzOvsQPU5Gk3K0AAKBskqpI4ACkgF5fEMAEBcHjz8O335r5pidPp00rW3hmiS3c57Lly/j5eXFoEGD+OOPP3j++efJlSsXn3zyCWvWrCE4OBilFNHR0WzcuJFXX/2O99/vQlycGUkwKOgo69f7sG/fPtrY5qr86KOPaNeuHZ999hkAXbt2pX///nTr1o0+ffqQJ08eC4/Yg2SmazjwZFa6lDvrJY9wWOOvv7Ru2NA8ovHhh1onJ1sdkWfAuY9nSW57uG3btmk/Pz+dJ08eXaBAAV2kSBF95syZlOVXrlzRtWvX1uvXr9daa71x40ZdpMgw7eWVoKtV0/qXX8x63333nfb29taHDh1K+W2TJk3SffyqVq1aeu/evY49MDeU1dzO7MhkM5VSDZVSrZRSrW29REUOFhEBjRvDoUNmwIPnn5ee3e5IctuzHTp0iK5du/LUU09x6dIlRo8eTdmyZXnyySdT1vH29qZXr15s3LiR5GRYsqQuMTGzUGoHU6aEUbeuJjo6msmTJ1O6dOmUzmIAxYsX559//kmzzytXrnDmzBmKy/OYdpOpQq2UWgg8A8wHhgB3OjAm4SSJiYlMnTqVe+65h6pVqzJ8+PAbki4kJASlVJpX8eIjaN0a8uc3vbwffDBz+5oyZUrKvp588kn+/fdfBx2ZyCzJbc82ffp0hg0bRkBAAF5eXlSrVo06deoQFhbG4cOHU9Y7fvw4RYqUpXdvmDGjJGXLhvDggx/zyiuDKVeuHBUrViQyMpI5c+ak2f7QoUN5++23U/7dSExMZOzYsTRt2pQKFSo49Vg9WWY7k1XRWg8GjmvTQzRbXYWUUo8qpSKUUjuVUlPTWT7StjxSKfVidvYlbu7JJ59k1apVzJw5k3Xr1lGsWDFatGjBuXPnUtbp3Llzyvt8+QoAb3H+/GfkzbuX8HC4667M7Wvo0KGEhIQwa9Ys1q5dS5EiRWjRogUXLlyw70GJ22W33Ja8dj2HDh2icePGKZ8ffvhhtm3bRtmyZTl06BAAa9euZfXqSBYsGMw338CUKRAZ2YjLl2M4f/48Pj4+5M+fn1mzZtGhQ4c02+/Zsyf9+/fnrrvuolWrVlStWpXt27czb948px6nx8vM9XEgDKgHrAHqA/uycp3dtq0qwO+AD6YDyxKgV6rlfsA2zHCG3sBmwPdW25T7WLfv8OHDumTJkvrixYtpvu/fv7+ePHmy1lrrQoUKaUCPGDFCx8Vp/eijWoPWMFeDd6b3dfDgQV26dGkdFxeX5vuAgAA9derU7B+Mh8G596jtktuOyGstuZ1tr776qh41alSa73788UedK1cuXbNmTX333XfrsmW76JIlL+tChbReuTLt748fP6537dqlL1++fMv9nDt3Tv/00096//79dj4Cz5LV3M7sGfUgIA/m0Y35wPTbag2k1RFYqrWOtgU+C+ieankXYJ7WOkGbIQ3nAt2ysT+RjsjISPz8/ChYsGCa7zt06MDu3bsBuHjxIgBvvDGD1q0hOBgmTwbz1yEh0/vavXs3LVq0oECBAjfdl7CMvXJb8toFPfXUUyxZsoQpU6YQFRXFrl27eO+99xg0aBBfffUVAQHfEh29goIF87J1K3Tpkvb3lStXpmHDhhk+YlWsWDHatm1LHXku0yEy25nskNZ6p9Z6g9a6odZ6Zjb2WQI4merzCaD0bSwHQCk1TCm1Qym14/Tp09kIJ2eqWrUqe/bsISkpKc33kZGRVK1aFcA26fu9NGyYyIEDsGwZvPTS7e+rWrVq7Nmzh+Tk5JvuS1jDjrltl7wGyW17qlixIhs2bCAiIoI6derw6KOP4u/vz4wZnxES4strr9WkYUNFRATcc4/V0YqbyWxnspZKqWVKqfXXXtnYZxRpE7Ss7bvMLgdAaz1ba+2rtfYtVapUNsLJmRo2bEjVqlUZNWoU0dHRJCUlERwczIIFCxg61Iwi+cUXF4BNnDr1Dw0bPkOpUttSxgC+nanqGjVqRPny5Xn++edT9rVkyRIWLVrEkCFD7H9wItPsmNt2yWuQ3La32rVrExQUxLlz5zh06BDPPfcqjz2WyJtvwr33/sInn/xC6XSbTMJlZOb6OBAJ3AdUv/bKynV227bKAb8ChW2fA0l7L8sXc98sD5AbM/2e3KN2gLNnz+qAgABdqFAhXaRIEd2oUSO9efNmnZys9fjx5n504cJ7NZTWmFmVUl6368yZM7p3794p+2rcuLHeunWrA47K/eHce9R2yW1H5LWW3La7f//V+s47z2tI0g0bLtbPPjtSly9fXo8ePVrHx8frefPm6T59+ujhw4frbdu2WR2ux8lqbmc2Cb8HcmVlBzfZXl9gNxAOTLF9FwqUtb1/0bZ8O/BCRtuTZM6eixcvpgyAcOmS1n36mL8Z/fppHR9v1hkxYoRu0qSJXfcl0ufkQm233LZ3XmvJbbvatUvrChWSNFzUkyYdSPn+/PnzumbNmrpBgwa6TZs2et68eXry5Mm6QoUKetq0aRZG7HmymtvK/PbWlFI9gXdtre9rZ+IuM9K6r6+v3rFjh9VhuL2TJ6F7dwgPh4kTYcwYGcTECkqpnVprXyftS3I7B/juOzMWf4EC8VSrNoqIiNlplgcEBBAaGsrff/9t65sCR48epWHDhhw9epRixYpZEbbHyWpuZ3b2rHHAm4D07PBQkZFmDumzZ8243enMeic8k+S2B9MaJk2C116Dpk1hxIgwFi48esN6v//+OzVr1kwp0mA6gTZt2pTNmzfz0EMPOTNscZ3MFuo/gO+01kkZrinczrJl0LcvFC8OmzdDw4ZWRyScSHLbQ12+DEOHwsKFEBAAc+ZAYqIfo0cHsHPnTho1agRAbGwsR44coVevXjds49SpU/j4+Dg7dHGdzBZqL2C3UuqXa1+40uUxkTVaw/vvw6uvQpMmpmCXK2d1VMLJJLc9UFSUuSq2bRu8+645oza3sQoxZ84cOnToQI8ePShevDhLlizB39+fkJAQjh07lvLI5IIFC4iNjcXPz8/KQxFkvlB/6NAohNNduWJa24GB0Ls3zJ1rxu4WOY7ktofZuxceeghOn4ZvvoHrT5S7d+9OkyZNWLJkCRcvXmTp0qX4+voyffp0GjZsSJMmTTh16hSxsbEsW7aM3LlzW3MgIkWmOpO5OulwcntOnTKt7a1bYfx4eP116TTmSpzZmczVSW6nT2vNkiVLmDVrFqdPn6Z169aMGTOGtWvz8uyzxSlWTLFqVW7q109iw4YNREVF0aJFC6pUqXLL7Z4/f55Nmzbh4+NDixYtpEjbmUM6kymlvtdad1RKncA8PwtmHF+ttS6fhTiFxVK3toOD4eGHrY5IWEFy27299957LFq0iAkTJlClShUCAxdyxx0zSUwcT6FCfxAX14OPP27Mzz//jI+PD3fccQejRo3iiSee4IMPPkgZuOh6xYoVo2vXrk4+GpGRjC599wLQWsudSw+wcqXpVFKkCISFga+cs+VkkttuSGvN+fPn+eCDD9i3bx8VK1bkyhUICblCYmIzatWKZNeuBly5spUqVarQtm1bli1bBsCFCxdo06YNQUFBPPbYY9YeiLgttxxCVGsdB6CU+ir190qpWY4MStiX1mbqum7doHZtiIiQIp3TSW67l+TkZN5//30qVqxIiRIlSEhIYN++fZw+Df7+mgMHmtGr1z7KlRtNgQKmt3b+/Pk5cOBAyjaKFi3KmDFjCAwMtPBIRFZkdOn7bkzLu5lS6k3b13mBdo4OTNjHlSswYgTMmwePPALz58N1k1iJHEhy272MGzeOH374gbVr11KwYEEaNGhAv36T8Pa+nwsX8uLl1Q8/v0bs3l0RgLi4OHx8fDh79mya7RQtWpS4uDgrDkFkQ0aXvv8EjgFXgOO27zQgNzHcwOnTpsfnpk3w1lvw5puQK7MTmwpPJ7ntJuLj45k2bRq7d++mcuXKANSp8yI7dozGyyuOTZvyMWDALsaPX8OaNWsAqFevHlFRUSnPSoO5bD579mw6d+5syXGIrLtlodZaxwBfKqWWaK0vAyiliti+Fy7s119Np7F//4XFi80jWEJcI7ntPk6ePEmhQoWoXLkyWsNHH0FExBsULnyY2Ni29OoFCQkJJCYmsnbtWs6dO8dPP/1Erly52L17Ny+//DLVq1cnKCiIS5cu8dRTT1l9SOI2Zfb8apJSqqJSqhewVCm1yJFBiewJCYHmzSE+3nQakyItbkFy20X9888/hISEEB0dTVxcHAcOHGHYMHjhBejZU/HOO6Hcf39N/ve//xEcHEyvXr1YunQpb775Jrlz52bPnj1s374dLy8vIiIiGDBgAKGhoRQsWNDqQxO3KbMDntTTWv+tlHpVa91eKbXZoVGJLNEaPv4YXnwR6teHFSugYkWroxIuTnLbxSQnJ/Pcc8+xcOFCfH192b9/Pz4+d+Dre5a4uDt46aUrVKs2n1GjnqVQoUI8++yznDhxgs6dOzNgwADmzZtHVFQUlSpVQinFxIkTrT4kkU2ZPaPOrZR6H4hQShUFLjkuJJEVCQkwfDiMHm1mwNq0SYq0yBTJbRczY8YMdu7cydGjR1m3bh0hIcc4c2YVly7Vo1y5F5k2rRivv/4qPXr0IDIykri4ONatW8fevXtp2LAh27dvZ+fOnaxbt87qQxF2ktlC/RiwUWv9pe03ox0XkrhdZ8/CAw/A55+bUcaCg0GubolMktx2MXPnzmXixIn4+Pjw/ffQooUX+fKVJF++Thw8OI5ff/2VPHnyEBgYyMaNG2nfvj3+/v689tprzJo1i/z58zNgwABWrVpl9aEIO7lloVZKPQ6gtT4JbLW9PwdIt0EXceCAmb5u2zYzS86770rPbpExyW3Xdf78ecqUKcunn8KDD0K1arBlSyK5ckXw1FNP8dZbb1GyZEm8vb3Jnz8/MTGm/1/16tWJiooCICYmhvwyeL/HyOif9CdSvf821fuO9g9FXHP69GlmzJjBpEmT2LVr103XW7cOmjWD2FjYsMFMVSlEJj2R6r3ktgvx9+9IQEA0I0eaJzc2bLhK7973kZyczL333kvZsmX57bffmDJlCp06dWLnzp1s2LCBoKAgWrVqxbFjx/jiiy/oK/8geIyMOpOpTLwXdrR27VoCAgLo3LkzpUqVonv37nTp0oXp06enjM+rNUybBs89B3ffbTqNZTDWvhDXk9x2QefOwW+/fcTu3fmoV28NffvGM3BgIHv27CEkJIR27cx4NHnz5uWVV14hf/78jBkzhs6dO6OUwt/fnwYNGjBhwgTq169v8dEIe8moUOtMvBd2cvnyZR5//HGWL19OixYtADMiUfPmzVm1ahUPPfQQV6/CqFHw2WfQtSssWgSFClkcuHBHktsu5vffzRn08eP5+PTTGM6f38lXX+1i//79jB8/PqVIA4wfP57ly5ezePFicuXKxciRI6lduzbe3t7MmTOH0qVLW3gkwt4yKtQtlFL/YlrZxVO9L+bwyHKgsLAw7rzzzpQiDVC4cGGefvppvvnmG1q0eIhHHoGffoJXXoGJE+V+tMgyyW0X8uOP8Mgjmty5k/nxR0XLlkWAsQAMHDgw3Wefc+fOzcSJE2nVqpWToxXOltGkHN5a6/Ja63Ja67yp3udzVoA5idY63ennlFJER5ehWTMzgMn8+TBpkhRpkXWS267jww8v06FDEjEx+9G6CYMH12blypUpy/v168fHH3+c0lEMIDg4mAsXLuDn52dFyMLJMjvgiXCC1q1b079/f7Zu3cp9990HwMWLF5k8eSenTk0nf35Yvx5SnXALIdxUYqLpZzJ9ej4qVoxky5YqVKq0g9DQUPr06cPq1atp1KgR/v7+PPHEE9SpU4f27dtz8uRJDh8+zIoVK8idO7fVhyGcQGnt/rekfH199Y4dO6wOwy5CQkLo168fXbp0oXTp0syd68358+9w1125WLFCUa2a1REKR1NK7dRay0SkeFZup3bhAjz6KPzwA+TPP43Tp4dQsOB/FzOmTp3Kr7/+yty5c1O+++eff9iwYQNFixalQ4cOeHt7WxC5yI6s5racUbuYTp068dtvv/HVV0EsWdKUc+ca8+CD8NVXUKSI1dEJITLj6tWrfPvtt2zatIkyZcowYMCAlJmvDh2CLl3gyBEYM+YgoaGLKFjwmTS/r1evHqtXr07zXYUKFejXr5/TjkG4DinULsjbuzQhIc/w889mAP733we5wiWEe7h06RIPPPAAAL169eLIkSPce++9LF68GC+v9vTqZfqX/Pgj3HNPSWbOPMDJkycpW7ZsyjbWrFlD48aNrToE4WKkULuYQ4fMIxqHDsEXX8DgwVZHJIS4HTNmzKB48eJ899135LL1+OzWrRuPPvojMTHtuPNOxcqVcMcdAMUYOXIknTp14r333qNKlSp8/fXXBAUFERERYelxCNchhdqFhIZCr16glGltt25tdURCiNu1atUqXn311ZQinZgIK1f6c+6cP35+MaxeXQQfn//Wf/vtt6lWrRrjxo3jzJkztG7dms2bN1OhQgWLjkC4GinULuLzz+Gpp6BmTVi5EqpXtzoiIURW5MuXj4sXLwIQHQ19+pg54gsV+oJZs1ri45O2s4lSiieeeIInnnjCgmiFO5AncS2WlATPPw/DhoG/v5lcQ4q0EO4rICCASZMmsXdvHM2bw7p1ydSoMRkvr5dYv34dcXFxVoco3IxTC7Uy3lNKhSulIpVSN4war5TyUkqdUUqFpnp55HMIMTHmfvTHH5thQVetIs0lMSHcgeR1Wv369aNChT40aHCZP/6IQamOnD07ibFjx/LDDz/Qpk2blDNuITLD2WfUAUBNoBnQCnhdKVXuunUqAeu01m1SvRKcHKfDHTkCzZub5yhnzjTF2ktuRAj3JHmdyvz5uVizZjRVqhTE27sln3zSg6ioKEaPHs3y5cupVKkSX3zxhdVhCjfi7ELdBZitjRjgG26c/7YqUFopFaKU2qSU6u3kGB1u0yYzh/SJE7B2LQwfbnVEQmSL5DXmNtZLL5knNe6/Hz79dAf16xdgxIgR5MmTBzD3ox9//HHWrVtncbTCnTjkHE4p1RZ4M51FCcDJVJ9PANdP83IJCAXeAwoB65VSe7XW+6/bxzBgGJAykIA7mDfPFOZq1cyl7po1rY5IiMxxRl7b9uN2uR0TY+aDX7UKnnkGPvoI9u8vQlRUFMnJySk9wAFOnjxJ0aJFrQtWuB2HnFFrrddfd4mrjda6DRBF2gQua/su9W/Dtdbvaq2TtNbRwE9Ao3T2MVtr7au19i1VqpQjDsOurrW2Bw0yj139/LMUaeFenJHXtnXdKrePHQM/P9Oze/p0+PRTcxvrnnvuoXjx4nz00UdcG6r5r7/+YvLkyQwaNMjaoIVbcfal7+XAYAClVAGgJxCSegWllN+1y2JKqbxAG2C3c8O0r9hY6NEDpkyBp5+GNWugmEwmKDxHjsxrgC1boEkT+Ptv+P5784jlNUopgoKCCAwMpE6dOrRr14569erx9NNPp5lbWoiMOLv70lKguVJqB2aC+kla6xNKqQbAGK11b+A3YJRS6gUgEXPv6xcnx2k3x45B166wfz9Mm2YKtRAeJsflNcCCBTB0KFSpYsY+qFXrxnXuuOMOdu/ezfbt2zl79ixNmzalePHizg9WuDWZPcuBtmwxZ9IJCRAcDO3bWx2RcAcye9Z/XDG3k5PhtdfMGPxt25rcltorMiOruS0DnjhIYKBJYh8fcz+6fXszo86ECROoUaMGpUqVom/fvhw+fNjqUIUQmXTxIvTsaYr08OHmcrcUaeFoUqjt7Fpr+/HHTQeT8HCoXdssGzp0KGFhYQQHBxMZGUndunVp3bo1p06dsjZoIUSG/vwTWrQwl7k/+QQ++wxsT10J4VAyxIYdXbwI/fvDsmVmSNBp0/5L5CNHjrB69WqOHz9OgQIFAHj99dc5fvw4s2bNYuzYsdYFLoS4pZ9/hu7dIT7edAa1zWIphFPIGbWdXGttr1gB//ufGW0sdWv7l19+oWnTpilF+hp/f3/27t3r5GiFEJn11VfQpg0ULGgKthRp4WxSqO3g55/NIxpHj8Lq1TBypJmqMrXq1asTGRlJYmJimu+3b99OjRo1nBitECIzkpPhjTfMQCbNmkFEBNSpY3VUIieSQp1Nixf/19retg06dkx/vbvuuov69eszZMgQTp06xdWrVwkMDGTBggUMlzFEhXApcXHwyCMwYQIMGQLr1kGJElZHJXIqKdRZlJwMY8dCQIAZtzs8HOrWvfVvvv76a/Lly0f16tUpXLgwn3/+OatXr6Zq1apOiVkIkbG//4ZWrUxfkw8/hNmzwdsj5/kS7kI6k2VBXBwMGABLl5ohQT/7LHOJXLhwYWbOnMm0adNISEi44X61EMJaERGm09jFi6a/yYMPWh2REHJGfdv++ce0tr/9FqZOhS++uP3WtpeXlxRpIVzMkiVmHP68eWHrVinSwnVIob4NO3ZA48bwxx+mtT169I2dxoQQ7kVrePtt6N0bfH3NWfXdd1sdlRD/kUKdSUFB0LKlaW1v2wZdulgdkRAiu+LjTYEeNw6eeAJ+/BHcYMIukcNIoc6A1iaJH3sMGjUyncaktS2E+/v3X3MbKzgYJk+GuXNNQ1wIVyOdyW4hPh4GDjT3rh5/3PT+lEQWwv3t3GlmtYuONr27u3a1OiIhbk7OqG/ixAnTsSQoyAzAP3++FGkhPME335jbWF5eptOYFGnh6qRQp2PXLtNpbP9++O47ePll6TQmhLvTGt591wxk0qCB6TRWr57VUQmRMSnU11m61IzZnSuXmU+6WzerIxJCZFd8PPTrZwYp6tcP1q+HMmWsjkqIzJFCbaO1GS7w4Yehfn3T2q5f3+qohBDZdfIk3H+/mVxj4kRYsADy5bM6KiEyTzqTAZcvm/F8Fy0yA/B/8YUkshCeIDLS3IM+e9YMUtSjh9URCXH7cvwZ9bXW9qJF5ow6MFCKtBCeYNky8PMzV8s2b5YiLdxXji7Ue/aY6Sn37jX3pl97TTqNCeHutIZJk6BnTzPmQUQENGxodVRCZF2OLdTLl5vWdnIybNpkkloI4d6uXDET5rz6qhmkKDQUypWzOiohsifHFWqtzXPRPXqYaSm3b4d777U6KiFEdp06BW3bmttX77xjOo/lz291VEJkX47qTHblCgwbZnp9PvYYzJsniSyEJ9i3Dx56yBTr4GDz9IYQniLHnFGfOgX+/qZIjxsHixdLkRbCE6xcCffdB1evQliYFGnheXJEod63z3Qa27nTjNv95pvSaUwId6c1TJliBiWqVct0GvP1tToqIezP4wv1qlWmtZ2QYFrbjz5qdURCiOxKSIDBg+Gll8wZdFgYVKhgdVRCOIbHFmqtYepUM9jBnXeaTmONG1sdlRAiu86cgXbtTB+Tt96Cr7+GAgWsjkoIx/HIzmQJCTBihJlf9uGH4csvJZGF8AS//mo6jf37r+ln0ru31REJ4Xged0Z95gy0b2+K9Nix5p60FGkh3F9ICDRvbibY2LhRirTIOTyqUO/fD02bQni4GRL0nXfMLFhCCPelNXz8MXTpAtWrm05jTZtaHZUQzuPUMqaUyquUGqmUClNKLb7JOkop9Z5SKlwpFamU6puZbV9rbcfFmdZ2QIB9YxdCpM+Rea01DB8Ozz8P3bubMbsrVbJr+EK4PGefbyYCB4D3gJs9IBUA1ASaAa2A15VStxwE8NQp09quVs10GpPWthBO5ZC8BvjjD/j8czMOf3AwFCxor5CFcB9OLdRa6ySt9Tog/hardQFmayMG+AbofKvt/vWX6d0trW0hnM9ReQ3mCllgoJnZTm5jiZzKIb2+lVJtgTfTWdRba30yg5+XAFKvcwIonc4+hgHDbB+vLFumfilcOCvRWqYkcMbqIG6TxOwctawOID3OyGvbftLkdv/+6pf+/W83Wku54985idk5spTbDinUWuv1wPos/jyKtAlcFjiezj5mA7MBlFI7tNZuNSaRxOwc7hqz1TGkxxl5bduP5LaTSczOkdXcdsWLScuBwQBKqQJATyDE0oiEENkleS1EFrlEoVZKlVVKhdo+LgX+tbU8NgKTtNYnLAtOCJElktdC2IclI5NprUOB0FSfTwJtbO818MJtbnK2nUJzJonZOSRmJ3FAXoN7/llIzM6RY2JWJn+EEEII4Ypc4tK3EEIIIdInhVoIIYRwYW5ZqB05ZKGjZCYepZSXUuqMUio01cvbglgfVUpFKKV2KqWmprN8pG15pFLqRWfHl55MxBx63auJFXFeF9PDSqkgpdSfN1l+y2PyNJLXTolXctvBHJLXWmu3ewG5gQ5AJ+Drm6zTFzP6kQKKAPuBchbGnGE8QDXgK4v/bKsAvwM+tliXAL1SLfcDtgHettdmwNeVY7ats83KGG8Sd2vMoA0ns3JMnvaSvHZ4rJLbzonZ7nntlmfU2oFDFjpQZuKpCpRWSoUopTYppayYyK8jsFRrHa3N36xZQPdUy7sA87TWCVrrBGAu0M35YaZxy5iVUl5AUVsrN0wpNV4plduiWFNorTdqrW82slJG/x88juS1w0luO4Ej8tqSx7MySzlpyEJ7ukXMCZmI5xLm8Zb3gELAeqXUXq31fgeEejMZ/bmVwLS6Uy+3ehqUjGIuhPlzfRWIwTwiMQSTJK7Kkr+/ziB5bUleg+S2K8jS31+XLtTaSUMW2tPNYlZKBWYUj9Y6HAi3fYxWSv0ENMJcTnOWKMylumvK2r5Lvfz640i93Aq3jFlrfQEYce2zUupboBeum8yQ8f8HtyV5bUleg+S2K8hSXrvlpe9McrUhCzOMRynld+2ymFIqL2awiN3ODZM1QA+l1LUpTgZhYr9mOfC4UiqP7RLTAGCFk2O83i1jVmaErNeUUtemYOwI7HJyjLcro/8POZXkddZJblsvS3ntUYVaufaQhenGo5RqoJT62rbOb0BPpdR2zOWc2VrrX5wZpO3PaCIQppQKB6K01kttvSnLaq13YJI3AvgZWGn7zjIZxYxpsRYCdimlNmE6cbjkqEZKqa+VUg1udkwWh2cJyWv7kNy2TnbzWkYmE0IIIVyYR51RCyGEEJ5GCrUQQgjhwqRQCyGEEC5MCrUQQgjhwqRQCyGEEC7MpQc8EY6nlCoErLJ9rIoZaelf2+dOWutbDecohHBRktueQx7PEimUUm9jBpKfaXUsQgj7kdx2b3JGLdKllKoKzACOYQZsKIEt0W0D4R/SWldVShXEDN5fDtNiH6a1PmJN1EKIjEhuux+5Ry1uxRf4TGv96S3WeRX4VWvdCngR+NApkQkhskNy243IGbW4lb+11vsyWKcBUMY2uxBAXseGJISwA8ltNyKFWtxKQqr30UAZ2/vuwLXODXuBo1rrz5VSuTCTzwshXJvkthuRS98isxYDDyil1gO1gSu27ycCbZRSYcAmoLJF8QkhskZy28VJr28hhBDChckZtRBCCOHCpFALIYQQLkwKtRBCCOHCpFALIYQQLkwKtRBCCOHCpFALIYQQLkwKtRBCCOHC/g8ce+/j+HFclQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 504x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAADQCAYAAADI+yJFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8WgzjOAAAACXBIWXMAAAsTAAALEwEAmpwYAAA28ElEQVR4nO3deZxN9f/A8dfHMJYhZpAtjbQQSjL2fSt7dpLoG1lK9ZXSIkVkC31LSrJki0gypkiy1zD2REokZV8aBmNmzLx/f3wuv6ExM2bm3nPvzPv5eNyHe88598z7MG/vzznncz4fIyIopZRSyjtlczoApZRSSt2YFmqllFLKi2mhVkoppbyYFmqllFLKi2mhVkoppbyYFmqllFLKi2mhVlcZY5obY2pk4P7aGWMeyKj9KaXSJi25bYzJYYxpYoyZZoyJNsYUSrROc9uDtFBnIcaYg8aYHcaYDcaYLcaYr4wxd7jWlQAGARGuz3cbY5a7tt1qjHnbGJPk74sxpocxJsK1bbgxpqFr1XfAR8aYPJ44PqWyKjfldjDwOLASyHXdOs1tD9JCnfX0EpHaIhICrAM+dC0fAkwWkXhjjB8QBkwQkdpAdaAW0P36nRlj6gEvAg+5tu0JfGGMySciUcBSoL/bj0oplaG5LSK/i0h3EZmXxDrNbQ/SQp1FGWNyAMWBbK7kbYNNvCueFJEVACISB+wG7kliV78AnUQk0vX5N+zvVTHX58XAoxkdv1IqaRmY2ynR3PYQo0OIZh3GmIPAaSAWuA8Ix7aIBVgsIuVv8L3OwLtAfRH5LZn9BwCTgSARaZFo+QUgUERiM+hQlFKJeCC3BSgsIqeuW6657QF6Rp319BGRGkAR7H2mDkAh4J/rNzTGBBlj5gN9gDopJHIdYDNwENuCTywSKJgBsSulbswtuZ2CSDS33S670wEoZ4jIBWPMOGwLfB0QkHi9MSYYWAKMF5HZye3LGNMD6Ae0uUHC5wXOZUjgSqlkZWRup4Lmtgdooc7aWmFbxD8Btxlj/FwdTrJj7z89JyLrktuB65GPZ4F6InIhifVFgHNJrVNKuU26czslmtueo4U665lqjDkP5AaigOYictYYsxnbA/QHoAa2c8lIY0zi724WkQHGmHnAWhGZDPQCgoBvr9v2LVeHlYeAUHcflFIqw3M7JZrbHqKdyRQAxphawLMi0iUV2+4DmorI/hS2M9hLb91E5M+MiVQpdTM0t32fdiZTAIjID8B+V1LfkDHmfuC3lBLZ5TFgpiayUs7R3PZ9ekatrnK1kpEUfimMMTlcz1+mtL9sIpKQUfEppdJGc9u3efyM2hjTwRizwBhz6AbrO7mGo9xqjBnv6fiyMnFJxXYpJrJrO03kLERz23tpbvs2Jy59nwSeBvyvX+F6bGA40AQIwfZWbO/Z8JRSaaS5rZQbeLzXt4isBbiux+EVTYFFInLWtc3HwH+ARddvaIzpDfQGCAgIqFy2bFl3hayUR0RFwf79EB+/9ZSIFHY6npulua1U0v76C06cAEhbbnvb41kFgWOJPh8Fbk1qQxGZAkwBCAkJkS1btrg/OqXcZNo06NsX7r4b9u41mbGDjua2ynLOnYMuXWDrVnjuOXj//bTltrf1+j7Otclb1LVMqUwpPh4GDoRevaBRI9i40emI3EZzW2UpBw5AjRrw3XcweTK8917a9+VthfoboK0xJp/r85PYoe6UynTOnYPWrWHCBNvaDguD/PmdjsptNLdVlrF+PVSrBkePwrffQp8+6dufVxRqY8x8Y8wDInIUGAmsM8ZsAo6LyL/uYSnl6/74A2rWtEn80Ue2tZ3d225EZQDNbZXVzJhhr44FBcGmTdCwYfr36dh/DSJSNNH7LonezwXmOhKUUh6wYQO0bWsve3/7rU3qzERzW2VF8fHwyiswbhw0bgwLFkBgYMbs2yvOqJXKKj791Lawg4Ls/ejMVqSVyoqiomzje9w4eOYZ+OabjCvSoIVaKY+Ij4dBg+A//4G6dW2Rvucep6NSSqXXwYNQq5Ytzh98YF85cmTsz8iEd8WU8i5RUdCtG4SGQr9+9n50RieyUsrzfvjBnknHxsKyZdCkiXt+jp5RK+VGf/4JtWvD11/blvaHH2qRViozmD3b3sbKn99eIXNXkQYt1Eq5TXg4VK1qi/U339h7V0op35aQAK+9Bt2720vemzaBuwfP00KtlBvMmQP160O+fLa1/dBDTkeklEqv8+ehfXsYNQp697ZPbQQFuf/naqFWKgNdaW0//rh9TtoTrW2llPsdOmRvY4WG2n4mkyd77jaWdiZTKoNcuGAL9OLF8NRT9p60/7/mkVJK+ZqNG6FNG4iOtv1Nmjb17M/XM2qlMsBff9nW9pIl8L//wccfa5FWKjOYN8/exgoIsP1OPF2kQQu1Uum2aZPtNLZ/vx2v+/nnIemZHpVSviIhAYYMga5d7bjdmzZBuXLOxKKFWql0mDcP6tWD3Llta7tZM6cjUkql14UL0KkTjBgBTz5pZ8AqVMi5eLRQK5UGCQnw5pu2tV21KkREQPnyTkellEqvw4ft6IFffgnjx8PUqc7fxtLOZErdpIsX4YknYOFCOyTo5MnOJ7JSKv22bLFTz0ZF2d7dLVs6HZGlZ9RK3YQrre0vvrAD8E+bpkVaqcxgwQKoUwdy5rS3sbylSIMWaqVSbcsWqFIFfv3VtrYHDtROY0r5OhEYNgw6d4bKlW2nsQoVnI7qWlqolUqFhQvtmbS/P/z4o3e1tpVSaRMdDY8+CkOH2iFBv/8ebr3V6aj+TQu1UskQgeHDbQ/QSpVsp7H77nM6KqVUeh09ap/YWLAAxoyxc8XnzOl0VEnTzmRK3UB0tH00Y/5829qeMsV7E1kplXrbttlOY5GRdiTBRx5xOqLk6Rm1Ukm40tr+/HMYPdq7W9tKqdRbtMiOIpgtm51P2tuLNGihVupftm+3z0bv2WOfpXz5Ze00ppSvE4G334YOHaBiRXsbq2JFp6NKHS3USiWyeLFtbRsDGzbYgfiVUr7t0iU7Yc7rr8Njj8Hq1VC0qNNRpZ4WaqWwre2RI6FdO9tZLCICHnjA6aiUUul17Bg0aABz59oz6tmzIVcup6O6OdqZTGV5ly5Br142kbt2tYOY+FoiK6X+bedOaNUKTp+296bbtXM6orTRM2qVpR0/Dg0b2iI9YgTMmaNFWqnMYMkSqFXLjsu/fr3vFmnQQq2ysJ9+sp3GduywQ4IOHqydxpTydSL2uei2be20lJs3w4MPOh1V+mihVllSaCjUrAnx8bbTWPv2TkeklEqvmBg7Yc4rr9hBitauhWLFnI4q/bRQqyxFBMaOtb25773Xdhrz9da2UpnBuXPn2LlzJ2fOnEnT90+cgEaNYNYsO3b3vHl2nvjMQAu1yjJiYuxIYy+/DB072tZ28eJOR6VU1iYiDBkyhODgYLp168add95J//79iYuLS/U+du2yt7G2bbNDgr7xRua6jaWFWmUJJ09C48Z2hLGhQ+2woHnyOB2VUmrSpEl8++237Nmzh127drF//35+//13hg4dmqrvh4XZ21hxcbBunW2EZzYefzzLGNMJeBHwA9aIyMDr1q+57iuDRCTCQ+GpTOL48eMsXLiQCxcucNddbXjxxTIcO2YLdOfOTkeXOWluq7SYNGkS06dPp5jrZnJQUBAfffQRISEhDB8+nGzZkj6fFIEJE+Cll+yEOaGhUKKEJyP3HI8WamNMMDAcqAqcA+YbY9qLyKJEm+UUkRqejEtlLqGhoTzxxBO0atWKs2dr8+qrJciT5yxr195C1aqZ6HqYF9HcVimJjo5m6dKlnD59mnr16lGuXDkAjh07xt13333NtsHBwURFRREbG0uuJJ6XjI2Ffv1g+nQ7JOjMmZn7CpmnL303BRaJyFkREeBjoM2VlcaY7EABY8wCY8w6Y8xwY4yfh2NUPuzChQv85z//Yfnyb6lYcSahoU9x//25uPXWlpw7973T4WVmmtvqhrZt28add97J9OnT2bZtG40aNaJ///6ICLVq1WLRokXXbB8WFsZ9992XZJE+dQqaNLFFesgQO3FOZi7S4PlL3wWBY4k+HwUST9OdF1gDvIptlU8BemGT/hrGmN5Ab4Dbb7/dPdEqn7Ny5UoqVarGlClVmDbNPnY1c2Z2pkxpx8KFC2ncuLHTIWZWmtsqSQkJCXTp0oV3332Xzq77TlFRUdStW5eFCxcybNgwmjZtyunTp2nYsCERERG8/fbbzJkz51/72rPHjjR2+DB89hk8+qinj8YZnj6jPs61yVvUtQwAEYkUkX6uPxOAL7GX0v5FRKaISIiIhBQuXNitQSvfcfZsDrZtG820aXYA/gULICAATGbqAuqdNLdVkrZt24afnx+dOnW6uixfvny88MILzJ8/n8qVK7N27VoOHDjAc889x6ZNmwgLC6NJkybX7GfZMqhRAy5csE9sZJUiDZ4/o/4GWGmMGSMiUcCTwFdXVhpjirqWjXJdPmsKbPNwjMpH7dkDQ4c25Z9/Yhkx4gCDB5cG7POZkydPZuLEiQ5HmKlpbqskxcXFkStXrn81lnPnzk1sbCwA5cqVY+rUqUl+XwTefx9eeMFOmLN0KZQs6fawvYpHz6hF5CgwElhnjNkEHBeRRcaYNa5EPo69RLbNGLMeMNhLZEola/ly29q+eDEbY8ZsYsKEKjz55JO8+OKLlC9fnocfflgve7uR5ra6kZCQEI4fP8769euvLrt8+TIffvghjzzySLLfjYuDvn3hv/+F1q3tKIJZrUgDGNu49W0hISGyZcsWp8NQDhCBiRNhwADb2g4Nhdtvtz1JFyxYwPnz52natCkP+tDwY8aYrSIS4nQc3kBzO3NYvnw53bp1o1OnTtx+++0sWLCAYsWKsXjxYvz9/ZP8zpkztkf36tXw2mswfDjc4Ektn5HW3NZpLpXPiouDZ5+Fjz+GRx6xM1/lzWvXFS1alOeee87ZAJVSADRt2pQdO3YwZ84cTp06xfDhw2natCl+fkl3/N+713YaO3TIDgn6+OMeDtjLaKFWPunMGTsC0apVdgD+t9/2/da2UpnZbbfdxiuvvJLidt99Z3Pb39+eTdes6YHgvJz+16Z8zq+/QvXq9n7VrFkwapQWaaUyg0mToFkze/tq82Yt0lfof2/Kp3z3nS3SkZH2bDqrXxJTKjO4fBmeeQb694fmzeGHHyA42OmovIcWauUzPvzQtrZvu81OT1mrltMRKaXS659/bF5/+KEdt3vxYsiXz+movIveo1Ze7/Jl+3jGpEnQsqUdkUgTWSnfdOnSJY4ePUrRokX5++/ctGwJf/xhhwT9z3+cjs47aaFWXi0yEjp1spe8X3rJ3o++QUdRpZQXExFGjhzJhAkTyJs3L2fOPMDly/PImzc3339vqFPH6Qi9l176Vl5r3z57P3rNGpg2DcaO1SKtlK+aNGkSixcvZsOGDdx11zucP/8Fly4dIGfOumTPHu50eF5NC7XySqtXQ7VqdqaclSvhySedjkgplR7vvfce77//IY888gerVnWiYcN41q69THT0Hlq3bs2+ffucDtFraaFWXmfKFHjoIShWzHYaq1vX6YiUUknZtWsXw4YNY/jw4ezZs+eG212+fJkDB07TpEkM+/Y1p0KF75g58x9q1bqPf/75hz59+jBp0iQPRu5btFArr3H5MrRrd4g+faBgwa307j2D4sUvOR2WUuo6IsITTzxBtWrVWLx4MT/99BP169dnwoQJSW7fpctg4Eeio2twxx0jeeih5dSvX4fPP/+ckJAQqlWrpmfUydBCrTxGRNiyZQurV6/m4sWL16w7exbKl9/P4sW306DBDiZM+J2vv55H06ZNr86wo5TyDl27dmXOnDm88sordO/enR07dtC+fXtGjRrF/v37r9l27ty/+fLLlwkIKE3evO04fvxtOnfuTFBQEH379mX48OGsXr2a++67z6Gj8QEi4vOvypUri/Jue/bskfLly0uZMmWkZs2aUrBgQZk1a5aIiPz+u8jdd8cJxMrYsf9c/U58fLzUr1//6nZZBbBFvCCvvOGlue19IiIipECBAvL0009fXXbmzBkpVqyYdOrUScaPH391+dSpIn5+8ZIv31+yb59IeHi4lCpVSvLkySNlypSRevXqyfjx46VIkSJy6NAhJw7Ho9Ka28k+nmWMmQckOb2WiHR1Q7tBZULx8fG0bt2aQYMG0atXL4wx7N69m0aNGhEXV4NBg+4iJiaemjWH89JLI65+L1u2bDz22GOsXLmSx3UIsgylua3SatmyZVSsWJFsicbtDQwMpH379kRERFCrVi3i42HQIJgwAapVi+LEidaULr2Fu+6qzr59+xg3bhyjR4/mzz//pFChQqxevZqSWXH+ylRK6Tnqya4/2wDHgI1AM+Dijb6g1PXWrVvHLbfcwlNPPXV1Wfny5alefQpPPVWKe+6BV1/dysSJ3wIjrvnu4cOHCQwM9HDEWYLmtkqV33//ndDQUIwxnD59milTpnDu3Dk2btxIu3btaNCgAQAHDx7k559/pkmT9rRuDd98Y2e3Gz/+Fho0yMPAgQMZNmwYAQEB3Hnnnfj7+7Nz506CdazQlKXmtBsIve7zqrScvrvrpZfHvNsXX3whrVq1uvr58mWRgQNFQKRYsZ/kn39ELl++LKVLl5Zp06ZJQkKCiIjs3r1bihQpIjt27HAocmfgwUvfmtsqOe+8844UKlRI+vbtK2XLlpXs2bPLf//7XylQoIDUr19fjDHSuXNnadasmRhj5I03Zkj58iJ+fiIfffT/+zl58qR06tRJ8ubNK/nz55cHH3xQfvjhB+cOzCFpze3UJvNWIMD1/hZgW1p+mLtemsze7dixY1KgQAE5cuSInD0r0qKF/c277bYvZerUT69ut3v3bilbtqzce++9Urt2bQkKCpKZM2c6GLkzPFyoNbdVknbt2iVFixaVw4cPy19//SWBgYGyZcsWKVCggHzyyScSGBgoxYoVk4IFC0pAQIAMHrxMChUSKVBAZOXKpPd57tw5OXHixNXGeFaT1txO7RCi44AdxpjdQHng5fSfy6usokiRIrz00ktUr96Fy5e/5NixQMqV+5CCBRfQrdt3V7crV64cu3fvZvPmzZw/f57q1asTEBDgYORZgua2ukZMTAyjRo3i/fffJyYmhiFDhtCgQQMqV65M5cqVadGiBfHx8fz555+8/fbbrFq1ip491/PsszkpVQrCwuCee5Led758+cinA/XftFQVahGZZ4xZDtwJ/C4ikW6NSmU69eq9xtixL3LxYhxVq75J9+7FeeKJb8mZM+c122XLlo1q1ao5FGXWo7mtrte1a1fi4uLo0KEDfn5+BAYG8uqrrxITE0NsbCx+fn6ICPny5cPPz5+EhFH07ZuTRo1g4ULQLiUZL1WF2hgTCLwBlADGGWPyisgqt0amMo2ZM6F3bwgO9mfpUn/KlBnudEjKRXNbJbZr1y42bdrEgQMH+OWXX2jevDnbt29n//797Nmzh06dOrFmzRpGjBjBzJlfMm5cDWJjG9GvH7z3HuTI4fQRZE6pHfDkU2A1cCuwA3jTTfGoTCQhAV55BZ54AmrXho0boUwZp6NS1/kUzW3lsnPnTurWrYu/vz8VK1akb9++VKxYkTNnzhAXF8fy5cuJiYnh7rsb06dPOWJjH6J8+cmUK/cBly9HOx1+ppXaQn2LiIQCCSISyw2ev1TqivPnoV07GDMG+vaF5cshKMjpqFQSNLfVVXfeeSdbt24lISEBgCFDhrBy5UoiIyMpVaoUe/bsYcWKKLJn30psbHFq1x5JmTLfMWfOHJo0acKlSzrkrzuktlBfMMZ0A/yMMY2As26MSfm4Q4fsGfTSpTBxInz4oV4S82Ka2+qq6tWrU7hwYZ599lnOnDlDbGws4eHh/PXXX3z66af8+GNpGjf24+LF4+TL14R77/2L+++/n6ioKPbt28ecOXOcPoRMKbWFuifwMBAE/Ad4KvnNVVYVHg5Vq8Iff9gBD/r3B2OcjkolQ3NbXWWMITQ0lKioKEqWLEn+/PmZO3cu33yznIkTi/L443DXXSfJkaMOYWHjmDJlCm+++Sbbt28nX758TJ8+3elDyJRS+3hWRRG5OoajMaYOcMI9ISlfNXcu9OwJJUrY+aTvvdfpiFQqaG6rawQFBTFr1iymTp1KXFwcEMDjj8PixfDUU3DvvQsZPDiS2rVrX/2Ov78/lSpVYtu2bc4Fnoml9oz6lStvjDHZgXfcE47yRQkJ8Prr0K0bVK9u55DWIu0zNLdVkvz9/TlzJoDatWHJEnj3Xfj4Y6hSpSLx8fFMnjz5yqA57Nq1ixUrVlC6dGmHo86cUpqUozYwGyhqjDkAGCAB0Mc3FAAXLkD37vDll9CrF0yaBP7+TkelUqK5rVISEQGPPCJERSUwefJRnnrqNgBq1KhB4cKFGTNmDBMmTKBgwYLs27ePAgUK0Lt3b4ejzpySPaMWkQ0icgcwQERKi8gdInKniOh9LMXff0OdOvDVV3aWnClTtEj7Cs1tlZz586FOnXhOnjxEkSJtGTq0GjVq1ODAgQP4+fnxxRdfcPHiRYoXL07RokXJmTMnLVq0oEOHDk6HnimldmSyycaYSkA+bMu7pIho974szLa27Rl1aCi0aOF0RCotNLdVYgkJMGwYvPUW5MgRwdKl2WnWLJT4+Hjee+89WrVqxa5du6hevToHDhzgq6++4vTp0wwfPpz77rvP6fAzrVTdozbGzAH6YwdH6AXcYCRXldmEhYVhjLnmVaBAH+rVg9y5bS9vLdK+S3M7c9i7dy/9+vWjXr169OnThz179lxdt3PnTh599FHKli1Ls2bNWLFiRZL7uHgRunSxRbpChc288MIymjWrAoCfnx8DBgwgR44crF+/HoC8efPSrVs3nn/+eS3SbpbazmTBItIT+NPVQzRdXYWMMZ2MMRHGmK3GmPFJrH/OtX6HMebF9PwslT6tWrW6+t7fPxfwJmfPfoy//042bYLy5Z2LTWWIDMttzWtnREREUKdOHYoXL87QoUMpWbIk9erVY+PGjWzfvp3GjRtTrVo1vvzyS7p160bPnj2ZN2/eNfs4fBjq1oUvvhD69dtPjhx9OXPmGNHR/z/amDGG0qVLc/z4cU8fokrNFFvAOuB+4BugIrArLVN1ufYVDPwK5MdeavscaJ9ofS0gHPB3vTYAIcntU6fCc488efIIIL169ZILF0Q6dbLTU8J0AX+nw8u08Ow0lxmS2+7Ia9HcTpVGjRrJ9OnTr1k2c+ZMqV+/vrRt21Y++OCDa9Zt2LBBSpcuLfHx8SIisnmzSPHiInnzJkhIyFApW7asNGnSRAoVKiQlS5aUn3/+WURETp8+LYGBgXLo0CHPHFgmlNbcTm0S3gVUBhoA24G+aflhrn31AUYm+twQmJ3o8yigd6LPTwLDk9unJrN7YIeTlMOHRUJCRIwRGTv2/5cr9/Bwoc6Q3HZHXovmdqpkz55d9u7dK0eOHLm6LDo6Wvz8/KREiRJy8ODBa7ZPSEiQwMBAOXnypCxYIJI7t0hwsMiAAdOlWbNmEhcXJ9u2bZNbbrlFADHGSIMGDaRcuXIyaNAgDx9d5pLW3E7VpW8R+V1EtorIahGpJCKTU/O9GygIHEv0+Sh2QoDUrgfAGNPbGLPFGLPl5MmT6QhH3YgxBniQSpXi2LvX9u5+6SWno1IZKQNzO0PyGjS3b8bmzZsxxlCzZk0qVKhAvXr12L9/PwcPHqRw4cKULFmS3bt3X/Odv/76CxH44IMCdOoElSrZzqEbNnzEoEGDiIyMpEWLFgwePJihQ4fi5+fH3r17ERFGjRrl0JFmbantTFbHGPOVMWbVlVc6fuZxrk3Qoq5lqV0PgIhMEZEQEQkpXLhwOsJRNzJ1aiSwnhMnjnDffX0JCtrgKt7o5O+ZRAbmdobkNWhup9bJkydp2bIlbdq0oWLFivzyyy+0adOGJk2a0L9/f/r27cvzzz/PgAED+OWXXwA4duwYTz75DMWKrWLYsOx07w6rVsGtt0JCQgLZs2fn008/5eGHH2bQoEG8/vrrBAYG8uOPPxIfH09ERITDR501pXYI0YnA09wgsW7SN8BKY8wYEYnCXgL7KtH6JcAEY8xM7AAMPQDteOJBIvD22zBkyC0EBPzEhQtNCA8/QZ06H1/d5ty5cw5GqDJQRuW25rWHzZkzh+bNm/PJJ58wYMAAypQpw913382hQ4coU6YMgwcPJnv27Jw6dYoGDRrg7+/P2bN5yJPnW44fv53Ro2HQoP8fi79du3ZMmDCBoKAgqlevDsDcuXMJDg6mVKlSVK1alV9//fXqOuVBqbk+DiwHsqXl2voN9vcY9n7YJmCca9kaoKjr/Yuu9ZuBgSntT+9jZZyLF0UefVQERLp1E4mOtst79eollSpVcja4LALP3qPOsNzO6LwWze1kvfDCCzJ27Nirn0+cOCHh4eHSvXv3f3Ugi4mJkbCww1KiRLzkySOyePG/93fhwgWpX7++3HbbbVK+fHnp2LGjFC1aVLZv3y5xcXESHBwsW7dudfNRZW5pze3UJmA7YA/w2ZVXWn6Yu16azBnj6FGRatXsb8XIkSIJCU5HlDV5uFBrbvuozz//XGrXri0JiRL10qVLEhwcLFu2bLm67MSJE9K79zLJkSNWihaNke3bb7zPy5cvy/z58yV//vzSokUL+fXXX2Xv3r3SsWNHad68uRuPJmtIa26n9jnqYcAbwMeJXioT2bHDTk+5a5cdt/vVV3V6yixCc9tHtW3blri4OLp27Up4eDirVq2iZcuWVK1alcqVKwMQGrqU4ODJTJnSlKCgwyQkVOGzzwZdaaT9i5+fH507d2bnzp0UKFCAqlWr0qBBA0qVKsXChQs9eXgqkdTeo/4NWCwi8e4MRjnjq6/gsccgKAg2bLC9QFWWobnto3LkyMF3333H+PHjefrpp/H396dLly70798fgFOnztOxYzSxsUPo2hWmTStFdPQaqlSpQrNmzWjQoMEN9x0cHMycOTqSrLdIbaHODmw3xvx8ZYGIdHVPSMpTRGDMGHv2XLWqLdjFijkdlfIwzW0fli9fPoYOHcrQoUOvWX78ONSrF0NsbCdGjIDXXrNXyHLlCqRv374sWLAg2UKtvEtqC/UEt0ahPC4mxk4CP3u2Hd93+nQ7drfKcjS3M5mffoJWreDYsfxUrDicwYOHXLPez8/vhpe+lXdK7exZa90diPKcEyegbVv48UcYPhwGD9b70VmV5rbviIqKYt68eezdu5dy5crRpUsX8ubNyx9//MGsWbM4ffo0AQGPMnFidQoUMKxcGcMjj7zLjh2teOCBBwD7WOXHH3/MxIkTnT0YdVOS7UxmjFnu+vOoMeaI63XUGHPEM+GpjPbTT1ClCmzfDgsXwuuva5HOijS3fcvBgwe5//77Wb58OcWLFycsLIyKFSsybdo0qlSpwj//RHLgQAdGj66Gn98+NmyIpU6dAD755BMaNmzIU089xaBBg6hQoQJNmjShcePGTh+SugkmuUsgxpgAEbngwXjSJCQkRLZs2eJ0GF5v6VLo2hVuuQWWLIGQEKcjUkkxxmwVEbf+62hu+5aOHTtSsWJFXn/99avLXnvtNd59911WrFjL9OlV+fRT6NgxgRMnWtK1axt69+4NwOHDh/n888+5cOECzZo1I0QT3zFpzu3UPMPFdc9WAh+n5Vkwd730WcvkJSSIvPOOnVQjJETk77+djkglB88+R6257eUSEhLE399fIiMjr1m+dOlSMeZWqV1bBESGDrW5vmDBAmnZsqVD0arkpDW3k71HbYypALQHqhtj3nAtzgnodRMfERMD/frBjBnQsSN8+inkyeN0VMppmtu+JUeOHFy6dIn8+fNfXfbnn/kQ2ciWLcL8+YbOne3y8+fPkytXLociVe6Q0oAnh4CDQAzwp+v1K9DavWGpjHDyJDRpYov0m2/C/PlapNVVmts+whjDo48+yrBhw6721g4LEwYMqIoxuRg8eMXVIh0ZGcn48ePp2lWfsMtMkj2jFpFzwExjzOcicgnAGHOLa7nyYrt320c0jhyBefPsI1hKXaG57VvGjBlD8+bNqVjxAfLmHUJ4eDsCAg4xb14kzz/fk1Wr7uG2225j+fLldO/enTZt2jgdsspAqX2OerQxZhxQDehrjDkhIo+5MS6VDsuWQefOEBAA69bZwUyUugHNbR8QGBjIuHHv8+yz2QgPD6Fu3RN8/fU95M1raNlyP8uWLeP06dO88cYb3HXXXU6HqzJYagv1/SLytzHmVRFpYozZ4NaoVJqIwP/+By++CBUrQmgo3Hab01EpL6e57eXOnDlDy5Y92L79dS5dCiEgYAIlSmwjZ84ZQA5y5sypZ9CZXGon5fAzxowBIowxBYCL7gtJpUVsLPTpAy+8AG3awPr1WqRVqmhue7kePUaze/d0RKoydy6cOvU0Z86cYty4cU6HpjwktYW6M7BWRGa6vvOC+0JSN+v0aXj4YfjkEzvK2MKF9rK3Uqmgue3FliyJISxsMLlyFWTNGkPXrpArVy7eeustZs2a5XR4ykNSGpmsO4CIHAN+dL0/AzR3f2gqNfbuhWrVIDwc5syBESMgW2qbXyrL0tz2biIwcSK0a+dPtmx/Eh6eQPXq/7++YMGCREVFOReg8qiU/kt/ItH7LxO9b5rxoaibtWIFVK8OUVGwerWdqlKpVHoi0XvNbS8SF2fHPnjuOWjVyhASMoAff5x/zTZTp06leXNtU2UVKXUmM6l4rzxMBD74AP77X6hQwXYaCw52OirlYzS3vdCZM3ZgolWr4JVX4O23YevW0bRo0YKIiAgefPBBVqxYQXh4OOvXr3c6XOUhKZ1RSyreKw+Ki4NnnrGt7ZYt4YcftEirNNHc9jK//mqvkG3YADNnwqhR9jZWlSpV2LZtG0FBQaxatYoqVaqwfft2btPeollGSmfUtV2z6RggKNH7QLdHpv7ln39sa/v77+Hll2HkSL0frdJMc9uLrFxpcztHDns2nTv3Nl56aR6xsbG0bt2ahg0bMnToUKfDVA5J9r95EfEXkeIiUkxEciZ6rwPJethvv9nW9rp1drzu0aO1SKu009z2Hh9+CE2b2scpIyJg06YJtGrVijx58lC8eHGefvpp+vXrd3X4UJX1pHbAE+VhJ0+e5Pz58wQHB7N6dTY6dIDs2W1ru3Ztp6NTSqXX5cu2n8mkSfY21mefwdmzfzNixAh27dpFiRIlAHjmmWeoVKkS69ato169es4GrRyh52Re5sSJE7Ru3Zp77rmH2rVrU7ToUB56KOFqa1uLtFK+LzISmje3RfrFF+GrryBfPli2bBktWrS4WqQB8ubNS48ePQgLC3MsXuUsLdReRERo27Yt9957L4cOHaFt2785efItsmdfyYwZv3HHHU5HqJRKr99/t7ex1qyBadPgnXfAz8+uy507N+fPn//Xd6KiosidO7dnA1VeQy99e5Ft27Zx/PhxXn55FB06ZGPFChg4EPz91zNv3reEhIx3OkSlVDqsXg3t29v+JStXQt26dvnBgwf5+++/qVOnDs8//zzh4eHUqFEDgP379zNjxgzWrFnjXODKUVqovcjhw4cpWbIBtWpl4/ffYepU6NkTPvvsXpYsWeJ0eEqpdJgyxT5aec89sHQplC4N586do0ePHvzwww/ceeed/PLLL7Rs2ZIWLVpQs2ZNcufOzXfffcfYsWMpV66c04egHKKF2ovExNRg7draBAYmsHJlNq70G1m6dCnVqlVzNjilVJpcvgz9+0fz8ce5qVnzHF9/nY8CBey4Ms888wxBQUEcOnSIXLlyceTIEZo1a8bo0aPJnz8/sbGxTJo0iVtvvdXho1BO0kLtJT75BJ5+ujCBgccoUaI9Fy48xZYtt/Lpp5+ydetWPvroI6dDVErdpLNnoWrV3/ntt7soWfILjhx5lWbNCrFo0SICAgIIDQ29WqQBihcvzsiRIxk9erSOPKau0s5kDouPhwEDoHdvaNQI9u8vwvPPt2TkyJE8+eST5M6dmw0bNlCgQAGnQ1VK3YT9+6FChSh++y2YMWMiOXSoA/v3/8pDDz1Et27dOHv2LAEBAdxyyy3XfC84OJiTJ086FLXyRh49ozbGGGAk0BDICbwjInOv2yY7cAz4OdHih0Qk1mOBesi5c9ClCyxbBs8/D+PGQfbshp49e9KzZ0+nw1MqVTSv/23dOmjXzub4a6+tZdCgxgBky5aN119/nZIlSxIXF0fu3LlZv349da/0KgPmz59P/fr1HYpceSNPX/ruCtwNVAfyARuNMatE5GiibUoCK0Skq4dj86gDB6BVKzvi2OTJ0KeP0xEplWaa14lMnw59+9rOYqVK9aZRo6euWZ8jRw4KFSrE2bNnmTBhAp06dWLQoEFUqFCBr7/+mkWLFvHDDz84FL3yRp6+9N0SmCLWOeAL/j3/bSngVmPMMmPMemNMFw/H6Hbr19s5pI8ehW+/1SKtfJ7mNfY21ksv2Sc1GjSAjRuhVauyzJ49+5rttm3bxunTp6lQoQKPPPIIYWFh7Nq1i7Fjx5InTx4iIiII1pl2VCJuOaM2xjQE3khiVSz28tcVR4HruzNeBNYAo4C8wCpjzE8isue6n9Eb6A1w++23Z0zgHjBjhi3Md9wBYWFw991OR6RU6ngir10/x+dy+9w5Ox98WBj07w/vvmuH/H3uueeoXbs2nTt3pmPHjuzfv593332X//3vf/j7+wMQEhLCjBkzHD4C5c3cUqhFZBWw6vrlxpjZXJvARYE/r/vuJmCT6+NZY8z3QGVgz3XbTQGmAISEhHj9aPXx8XZ+2XHjoHFjWLAAAnWeIuVDPJHXrm19KrcPHrS3sX75xQ4J+vTT/78uMDCQH3/8kalTpzJnzhyKFCnC119/TeXKlR2LV/keT9+jXgL0BFYaY/IA7YCHEm9gjKkFlBSR+caYnEB9YJaH48xQUVG2tb10qR3w4N137XR2SmUSWTKvwc4H37atnSd++XLbCL9e/vz5GThwIAMHDvR8gCpT8PQ96kXAEWPMFmAtMFpEjhpjHjDGzHdt8wvQzhizGXupbIqI/Jz07rzfwYNQqxZ88w188IF9aZFWmUyWy2uAWbOgYUMoUMDej06qSCuVETx6Ri12QtV/NStFZAfQxfX+DNDJk3G5y5XWdmysfQSrSROnI1Iq42W1vE5IgNdegzFjbKFeuBCCgpyOSmVmOuCJm8yebZM4f37b2tYirZTvO3/ePh89ZoztFLp8uRZp5X5aqDPYldZ29+72kvemTVC2rNNRKaXS69AhOx/80qXw/vvw0Ud6G0t5ho71nYHOn4fHH7eTwPfurfejlcosNm6ENm0gOtr2N3n4YacjUlmJnlFnkCut7dBQeO89O9qYFmmlfN9nn0H9+hAQYAu2FmnlaVqoM8DGjVC1KvzxB3z9NTz3HBjjdFRKqfRISIDXX7ePVlavDhERcO+9TkelsiIt1Ok0b97/t7bDw6FpU6cjUkql14UL0LEjvP029OoFK1ZAwYJOR6WyKr1HnUp///038+bNIzIykiZNmlCnTj2GDjWMGAF168KiRVCokNNRKqXS6++/4ZFHYMcOmDAB/vtfvUKmnKVn1KkQGhrKAw88wP79+8mePTu9ew/gjjs2M2IEPPkkfPedFmmlMoOICHsba98+299kwAAt0sp5ekadgujoaHr27Mk333xD1apVOXwYli59k337oHv3HUyd+oAmslKZwOefwxNPQNGi9lJ3hQpOR6SUpWfUKdiwYQNly5alatWqbNkCVarAvn3ZePbZlcTFjdUirZSPE4GhQ6FLFwgJsWfVWqSVN9FCnQI/Pz/i4uJYsADq1IGcOW2nsXLlDpA9u16QUMqXRUfbAj1smD2bXrkSChd2OiqlrqWFOgW1atVmz54OdO4MlSvbkcZKljzL+++/T+fOnZ0OTymVRkeO2I6gCxfC2LEwfbptiCvlbbRQJyM6Gnr08Ccq6kVy5vycggU7MXTo05QtW5aHH36Y5s2bOx2iUioNtm61t7F++cWOJPjSS9ppTHkvvXZ7A0eP2kc0tmy5MgB/M7766hKRkZF8//33lCtXzukQlVJp8MUXdiz+woXhxx/h/vudjkip5GmhTsK2bdC6NURGwuLFtmDDLfTo0cPhyJRSaSViBzAZMgRq1LC5XaSI01EplTK99H2dRYvsmN3Zstn5pG2RVkr5suho6NbNFulu3WDVKi3SyndooXa50tru0AEqVrSPaFSs6HRUSqn0OnYMGjSwk2uMHAmzZkGuXE5HpVTq6aVv4NIlO57v3Ll2AP6pUzWRlcoMduywt7FOn4Yvv4S2bZ2OSKmbl+XPqK+0tufOtWfUs2drkVYqM/jqK6hVy14t27BBi7TyXVm6UO/cacf1/ekne2/6tdf0EQ2lfJ0IjB4N7drZEcYiIqBSJaejUirtsmyhXrLEtrYTEmD9epvUSinfFhMDPXrAq69C586wZg0UK+Z0VEqlT5Yr1CL2uei2baFcOdi8GR580OmolFLpdeIENGxob1+99ZbtPJY7t9NRKZV+WaozWUwM9O5te3127gwzZmgiK5UZ7NoFrVrZYr1woX16Q6nMIsucUZ84AY0a2SI9bBjMm6dFWqnMYOlSqFkT4uJg3Tot0irzyRKFetcu22ls61Y75+wbb2inMaV8nQiMG2cHJSpTxnYaCwlxOiqlMl6mL9RhYba1HRtrW9udOjkdkVIqvWJjoWdPO5lGhw42t0uUcDoqpdwj0xZqERg/3g52cM89ttNYlSpOR6WUSq9Tp6BxY9vH5M03Yf58yJPH6aiUcp9M2ZksNhb69bPzy3boADNnaiIrlRns3m07jR05YvuZdOnidERKuV+mO6M+dQqaNLFFesgQe09ai7RSvm/ZMjvrVXQ0rF2rRVplHZmqUO/ZA9WqwaZNdkjQt96ys2AppXyXCPzvf9CyJdx5p+00Vq2a01Ep5TkeLWPGmJzGmOeMMeuMMfNusI0xxowyxmwyxuwwxjyWmn1faW1fuGBb2127ZmzsSqmkuTOvRaBPHxgwANq0sWN2lyyZoeEr5fU8fb55GdgLjAJu9IBUV+BuoDpQFxhsjEl2EMATJ2xr+447bKcxbW0r5VFuyWuA336DTz6x4/AvXAgBARkVslK+w6OFWkTiRWQFEJ3MZi2BKWKdA74Amie337/+sr27tbWtlOe5K6/BXiGbPdvObKe3sVRW5ZZe38aYhsAbSazqIiLHUvh6QSDxNkeBW5P4Gb2B3q6PMV99ZX7Oly8t0TqmEHDK6SBuksbsGWWcDiApnshr18+5Jrcff9z8/PjjNxuto3zxd05j9ow05bZbCrWIrAJWpfHrx7k2gYsCfybxM6YAUwCMMVtExKfGJNKYPcNXY3Y6hqR4Iq9dP0dz28M0Zs9Ia25748WkJUBPAGNMHqAdsMzRiJRS6aV5rVQaeUWhNsYUNcascX1cBBxxtTzWAqNF5KhjwSml0kTzWqmM4cjIZCKyBliT6PMxoL7rvQADb3KXUzIoNE/SmD1DY/YQN+Q1+ObfhcbsGVkmZmPzRymllFLeyCsufSullFIqaVqolVJKKS/mk4XanUMWuktq4jHGZDfGnDLGrEn08ncg1k7GmAhjzFZjzPgk1j/nWr/DGPOip+NLSipiXnPdq6oTcV4XUwdjzAJjzKEbrE/2mDIbzWuPxKu57WZuyWsR8bkX4Ac8BDQD5t9gm8ewox8Z4BZgD1DMwZhTjAe4A/jM4b/bYOBXIL8r1s+B9onW1wLCAX/XawMQ4s0xu7YJdzLGG8RdDztow7G0HFNme2leuz1WzW3PxJzhee2TZ9TixiEL3Sg18ZQCbjXGLDPGrDfGODGRX1NgkYicFfub9THQJtH6lsAMEYkVkVhgOvCI58O8RrIxG2OyAwVcrdx1xpjhxhg/h2K9SkTWisiNRlZK6d8h09G8djvNbQ9wR1478nhWahkPDVmYkZKJOTYV8VzEPt4yCsgLrDLG/CQie9wQ6o2k9PdWENvqTrze6WlQUoo5L/bv9VXgHPYRiV7YJPFWjvz+eoLmtSN5DZrb3iBNv79eXajFQ0MWZqQbxWyMmZ1SPCKyCdjk+njWGPM9UBl7Oc1TjmMv1V1R1LUs8frrjyPxeickG7OIRAL9rnw2xnwJtMd7kxlS/nfwWZrXjuQ1aG57gzTltU9e+k4lbxuyMMV4jDG1rlwWM8bkxA4Wsd2zYfIN0NYYc2WKkyexsV+xBOhujMnhusTUAwj1cIzXSzZmY0fIes0Yc2UKxqbANg/HeLNS+nfIqjSv005z23lpyutMVaiNdw9ZmGQ8xpgHjDHzXdv8ArQzxmzGXs6ZIiI/ezJI19/RSGCdMWYTcFxEFrl6UxYVkS3Y5I0ANgJLXcsck1LM2BZrXmCbMWY9thOHV45qZIyZb4x54EbH5HB4jtC8zhia285Jb17ryGRKKaWUF8tUZ9RKKaVUZqOFWimllPJiWqiVUkopL6aFWimllPJiWqiVUkopL+bVA54o9zPG5AXCXB9LYUdaOuL63ExEkhvOUSnlpTS3Mw99PEtdZYwZih1IfrLTsSilMo7mtm/TM2qVJGNMKeBD4CB2wIaCuBLdNRD+7yJSyhgTgB28vxi2xd5bRA44E7VSKiWa275H71Gr5IQAH4nIxGS2eRXYLSJ1gReBCR6JTCmVHprbPkTPqFVy/haRXSls8wBQxDW7EEBO94aklMoAmts+RAu1Sk5sovdngSKu922AK50bfgL+EJFPjDHZsJPPK6W8m+a2D9FL3yq15gEPG2NWAWWBGNfykUB9Y8w6YD1wu0PxKaXSRnPby2mvb6WUUsqL6Rm1Ukop5cW0UCullFJeTAu1Ukop5cW0UCullFJeTAu1Ukop5cW0UCullFJeTAu1Ukop5cX+D8clfq2aABUaAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 504x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for t in range(1, n_timepoints):\n",
    "    B_t, B_tau = model.adjacency_matrices_[t]\n",
    "    plt.figure(figsize=(7, 3))\n",
    "\n",
    "    plt.subplot(1,2,1)\n",
    "    plt.plot([-1, 1],[-1, 1], marker=\"\", color=\"blue\", label=\"support\")\n",
    "    plt.scatter(B_t_true[t], B_t, facecolors='none', edgecolors='black')\n",
    "    plt.xlim(-1, 1)\n",
    "    plt.ylim(-1, 1)\n",
    "    plt.xlabel('True')\n",
    "    plt.ylabel('Estimated')\n",
    "    plt.title(f'B({t},{t})')\n",
    "\n",
    "    plt.subplot(1,2,2)\n",
    "    plt.plot([-1, 1],[-1, 1], marker=\"\", color=\"blue\", label=\"support\")\n",
    "    plt.scatter(B_tau_true[t], B_tau, facecolors='none', edgecolors='black')\n",
    "    plt.xlim(-1, 1)\n",
    "    plt.ylim(-1, 1)\n",
    "    plt.xlabel('True')\n",
    "    plt.ylabel('Estimated')\n",
    "    plt.title(f'B({t},{t-1})')\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Independence between error variables\n",
    "To check if the LiNGAM assumption is broken, we can get p-values of independence between error variables. The value in the i-th row and j-th column of the obtained matrix shows the p-value of the independence of the error variables $e_i$ and $e_j$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_values_list = model.get_error_independence_p_values()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.    0.026 0.064 0.289 0.051]\n",
      " [0.026 0.    0.363 0.821 0.581]\n",
      " [0.064 0.363 0.    0.067 0.098]\n",
      " [0.289 0.821 0.067 0.    0.059]\n",
      " [0.051 0.581 0.098 0.059 0.   ]]\n"
     ]
    }
   ],
   "source": [
    "t = 1\n",
    "print(p_values_list[t])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.    0.715 0.719 0.593 0.564]\n",
      " [0.715 0.    0.78  0.7   0.504]\n",
      " [0.719 0.78  0.    0.532 0.591]\n",
      " [0.593 0.7   0.532 0.    0.401]\n",
      " [0.564 0.504 0.591 0.401 0.   ]]\n"
     ]
    }
   ],
   "source": [
    "t = 2\n",
    "print(p_values_list[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bootstrapping\n",
    "We call `bootstrap()` method instead of `fit()`. Here, the second argument specifies the number of bootstrap sampling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = lingam.LongitudinalLiNGAM()\n",
    "result = model.bootstrap(X_t, n_sampling=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Causal Directions\n",
    "Since `LongitudinalBootstrapResult` object is returned, we can get the ranking of the causal directions extracted by `get_causal_direction_counts()` method. In the following sample code, `n_directions` option is limited to the causal directions of the top 8 rankings, and `min_causal_effect` option is limited to causal directions with a coefficient of 0.01 or more."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "cdc_list = result.get_causal_direction_counts(n_directions=12, min_causal_effect=0.01, split_by_causal_effect_sign=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x2(1) <--- x0(0) (b<0) (100.0%)\n",
      "x4(1) <--- x1(0) (b<0) (100.0%)\n",
      "x4(1) <--- x1(1) (b<0) (100.0%)\n",
      "x3(1) <--- x4(0) (b<0) (100.0%)\n",
      "x3(1) <--- x2(0) (b<0) (100.0%)\n",
      "x3(1) <--- x0(0) (b>0) (100.0%)\n",
      "x2(1) <--- x2(0) (b>0) (100.0%)\n",
      "x1(1) <--- x3(0) (b>0) (100.0%)\n",
      "x1(1) <--- x2(0) (b<0) (100.0%)\n",
      "x1(1) <--- x3(1) (b>0) (100.0%)\n",
      "x4(1) <--- x4(0) (b>0) (100.0%)\n",
      "x0(1) <--- x4(1) (b<0) (100.0%)\n"
     ]
    }
   ],
   "source": [
    "t = 1\n",
    "labels = [f'x{i}({u})' for u in [t, t-1] for i in range(5)]\n",
    "print_causal_directions(cdc_list[t], 100, labels=labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x0(2) <--- x0(1) (b>0) (100.0%)\n",
      "x4(2) <--- x1(1) (b>0) (100.0%)\n",
      "x4(2) <--- x0(1) (b<0) (100.0%)\n",
      "x3(2) <--- x2(1) (b>0) (100.0%)\n",
      "x3(2) <--- x1(1) (b<0) (100.0%)\n",
      "x3(2) <--- x0(1) (b>0) (100.0%)\n",
      "x3(2) <--- x2(2) (b<0) (100.0%)\n",
      "x2(2) <--- x3(1) (b>0) (100.0%)\n",
      "x2(2) <--- x1(1) (b>0) (100.0%)\n",
      "x4(2) <--- x3(1) (b>0) (100.0%)\n",
      "x1(2) <--- x3(1) (b>0) (100.0%)\n",
      "x1(2) <--- x2(1) (b>0) (100.0%)\n"
     ]
    }
   ],
   "source": [
    "t = 2\n",
    "labels = [f'x{i}({u})' for u in [t, t-1] for i in range(5)]\n",
    "print_causal_directions(cdc_list[t], 100, labels=labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Directed Acyclic Graphs\n",
    "Also, using the `get_directed_acyclic_graph_counts()` method, we can get the ranking of the DAGs extracted. In the following sample code, `n_dags` option is limited to the dags of the top 3 rankings, and `min_causal_effect` option is limited to causal directions with a coefficient of 0.01 or more."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "dagc_list = result.get_directed_acyclic_graph_counts(n_dags=3, min_causal_effect=0.01, split_by_causal_effect_sign=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DAG[0]: 2.0%\n",
      "\tx0(1) <--- x1(1) (b>0)\n",
      "\tx0(1) <--- x3(1) (b>0)\n",
      "\tx0(1) <--- x4(1) (b<0)\n",
      "\tx0(1) <--- x0(0) (b<0)\n",
      "\tx0(1) <--- x1(0) (b>0)\n",
      "\tx0(1) <--- x2(0) (b<0)\n",
      "\tx0(1) <--- x3(0) (b<0)\n",
      "\tx0(1) <--- x4(0) (b<0)\n",
      "\tx1(1) <--- x3(1) (b>0)\n",
      "\tx1(1) <--- x0(0) (b>0)\n",
      "\tx1(1) <--- x1(0) (b<0)\n",
      "\tx1(1) <--- x2(0) (b<0)\n",
      "\tx1(1) <--- x3(0) (b>0)\n",
      "\tx1(1) <--- x4(0) (b>0)\n",
      "\tx2(1) <--- x1(1) (b>0)\n",
      "\tx2(1) <--- x0(0) (b<0)\n",
      "\tx2(1) <--- x1(0) (b>0)\n",
      "\tx2(1) <--- x2(0) (b>0)\n",
      "\tx2(1) <--- x3(0) (b>0)\n",
      "\tx2(1) <--- x4(0) (b>0)\n",
      "\tx3(1) <--- x0(0) (b>0)\n",
      "\tx3(1) <--- x1(0) (b<0)\n",
      "\tx3(1) <--- x2(0) (b<0)\n",
      "\tx3(1) <--- x3(0) (b<0)\n",
      "\tx3(1) <--- x4(0) (b<0)\n",
      "\tx4(1) <--- x1(1) (b<0)\n",
      "\tx4(1) <--- x0(0) (b<0)\n",
      "\tx4(1) <--- x1(0) (b<0)\n",
      "\tx4(1) <--- x2(0) (b<0)\n",
      "\tx4(1) <--- x3(0) (b<0)\n",
      "\tx4(1) <--- x4(0) (b>0)\n",
      "DAG[1]: 2.0%\n",
      "\tx0(1) <--- x1(1) (b>0)\n",
      "\tx0(1) <--- x4(1) (b<0)\n",
      "\tx0(1) <--- x0(0) (b<0)\n",
      "\tx0(1) <--- x1(0) (b>0)\n",
      "\tx0(1) <--- x2(0) (b<0)\n",
      "\tx0(1) <--- x3(0) (b<0)\n",
      "\tx0(1) <--- x4(0) (b<0)\n",
      "\tx1(1) <--- x3(1) (b>0)\n",
      "\tx1(1) <--- x0(0) (b>0)\n",
      "\tx1(1) <--- x1(0) (b<0)\n",
      "\tx1(1) <--- x2(0) (b<0)\n",
      "\tx1(1) <--- x3(0) (b>0)\n",
      "\tx1(1) <--- x4(0) (b>0)\n",
      "\tx2(1) <--- x1(1) (b>0)\n",
      "\tx2(1) <--- x0(0) (b<0)\n",
      "\tx2(1) <--- x1(0) (b>0)\n",
      "\tx2(1) <--- x2(0) (b>0)\n",
      "\tx2(1) <--- x3(0) (b>0)\n",
      "\tx2(1) <--- x4(0) (b>0)\n",
      "\tx3(1) <--- x0(0) (b>0)\n",
      "\tx3(1) <--- x1(0) (b<0)\n",
      "\tx3(1) <--- x2(0) (b<0)\n",
      "\tx3(1) <--- x3(0) (b<0)\n",
      "\tx3(1) <--- x4(0) (b<0)\n",
      "\tx4(1) <--- x1(1) (b<0)\n",
      "\tx4(1) <--- x0(0) (b<0)\n",
      "\tx4(1) <--- x1(0) (b<0)\n",
      "\tx4(1) <--- x2(0) (b<0)\n",
      "\tx4(1) <--- x3(0) (b<0)\n",
      "\tx4(1) <--- x4(0) (b>0)\n",
      "DAG[2]: 2.0%\n",
      "\tx0(1) <--- x4(1) (b<0)\n",
      "\tx0(1) <--- x0(0) (b<0)\n",
      "\tx0(1) <--- x1(0) (b>0)\n",
      "\tx0(1) <--- x2(0) (b<0)\n",
      "\tx0(1) <--- x3(0) (b>0)\n",
      "\tx0(1) <--- x4(0) (b<0)\n",
      "\tx1(1) <--- x3(1) (b>0)\n",
      "\tx1(1) <--- x0(0) (b>0)\n",
      "\tx1(1) <--- x1(0) (b<0)\n",
      "\tx1(1) <--- x2(0) (b<0)\n",
      "\tx1(1) <--- x3(0) (b>0)\n",
      "\tx1(1) <--- x4(0) (b>0)\n",
      "\tx2(1) <--- x1(1) (b>0)\n",
      "\tx2(1) <--- x3(1) (b<0)\n",
      "\tx2(1) <--- x4(1) (b<0)\n",
      "\tx2(1) <--- x0(0) (b<0)\n",
      "\tx2(1) <--- x1(0) (b>0)\n",
      "\tx2(1) <--- x2(0) (b>0)\n",
      "\tx2(1) <--- x3(0) (b>0)\n",
      "\tx2(1) <--- x4(0) (b>0)\n",
      "\tx3(1) <--- x0(0) (b>0)\n",
      "\tx3(1) <--- x1(0) (b>0)\n",
      "\tx3(1) <--- x2(0) (b<0)\n",
      "\tx3(1) <--- x3(0) (b<0)\n",
      "\tx3(1) <--- x4(0) (b<0)\n",
      "\tx4(1) <--- x1(1) (b<0)\n",
      "\tx4(1) <--- x3(1) (b>0)\n",
      "\tx4(1) <--- x0(0) (b<0)\n",
      "\tx4(1) <--- x1(0) (b<0)\n",
      "\tx4(1) <--- x2(0) (b<0)\n",
      "\tx4(1) <--- x3(0) (b<0)\n",
      "\tx4(1) <--- x4(0) (b>0)\n"
     ]
    }
   ],
   "source": [
    "t = 1\n",
    "labels = [f'x{i}({u})' for u in [t, t-1] for i in range(5)]\n",
    "print_dagc(dagc_list[t], 100, labels=labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DAG[0]: 5.0%\n",
      "\tx0(2) <--- x0(1) (b>0)\n",
      "\tx0(2) <--- x1(1) (b>0)\n",
      "\tx0(2) <--- x2(1) (b>0)\n",
      "\tx0(2) <--- x3(1) (b<0)\n",
      "\tx0(2) <--- x4(1) (b>0)\n",
      "\tx1(2) <--- x2(2) (b<0)\n",
      "\tx1(2) <--- x4(2) (b>0)\n",
      "\tx1(2) <--- x0(1) (b<0)\n",
      "\tx1(2) <--- x1(1) (b<0)\n",
      "\tx1(2) <--- x2(1) (b>0)\n",
      "\tx1(2) <--- x3(1) (b>0)\n",
      "\tx1(2) <--- x4(1) (b>0)\n",
      "\tx2(2) <--- x0(2) (b>0)\n",
      "\tx2(2) <--- x0(1) (b>0)\n",
      "\tx2(2) <--- x1(1) (b>0)\n",
      "\tx2(2) <--- x2(1) (b<0)\n",
      "\tx2(2) <--- x3(1) (b>0)\n",
      "\tx2(2) <--- x4(1) (b<0)\n",
      "\tx3(2) <--- x2(2) (b<0)\n",
      "\tx3(2) <--- x0(1) (b>0)\n",
      "\tx3(2) <--- x1(1) (b<0)\n",
      "\tx3(2) <--- x2(1) (b>0)\n",
      "\tx3(2) <--- x3(1) (b>0)\n",
      "\tx3(2) <--- x4(1) (b<0)\n",
      "\tx4(2) <--- x0(2) (b>0)\n",
      "\tx4(2) <--- x0(1) (b<0)\n",
      "\tx4(2) <--- x1(1) (b>0)\n",
      "\tx4(2) <--- x2(1) (b<0)\n",
      "\tx4(2) <--- x3(1) (b>0)\n",
      "\tx4(2) <--- x4(1) (b<0)\n",
      "DAG[1]: 2.0%\n",
      "\tx0(2) <--- x0(1) (b>0)\n",
      "\tx0(2) <--- x1(1) (b>0)\n",
      "\tx0(2) <--- x2(1) (b>0)\n",
      "\tx0(2) <--- x3(1) (b<0)\n",
      "\tx0(2) <--- x4(1) (b<0)\n",
      "\tx1(2) <--- x2(2) (b<0)\n",
      "\tx1(2) <--- x4(2) (b>0)\n",
      "\tx1(2) <--- x0(1) (b<0)\n",
      "\tx1(2) <--- x1(1) (b<0)\n",
      "\tx1(2) <--- x2(1) (b>0)\n",
      "\tx1(2) <--- x3(1) (b>0)\n",
      "\tx1(2) <--- x4(1) (b>0)\n",
      "\tx2(2) <--- x0(2) (b>0)\n",
      "\tx2(2) <--- x0(1) (b>0)\n",
      "\tx2(2) <--- x1(1) (b>0)\n",
      "\tx2(2) <--- x3(1) (b>0)\n",
      "\tx2(2) <--- x4(1) (b<0)\n",
      "\tx3(2) <--- x2(2) (b<0)\n",
      "\tx3(2) <--- x0(1) (b>0)\n",
      "\tx3(2) <--- x1(1) (b<0)\n",
      "\tx3(2) <--- x2(1) (b>0)\n",
      "\tx3(2) <--- x3(1) (b>0)\n",
      "\tx3(2) <--- x4(1) (b<0)\n",
      "\tx4(2) <--- x0(2) (b>0)\n",
      "\tx4(2) <--- x0(1) (b<0)\n",
      "\tx4(2) <--- x1(1) (b>0)\n",
      "\tx4(2) <--- x2(1) (b<0)\n",
      "\tx4(2) <--- x3(1) (b>0)\n",
      "\tx4(2) <--- x4(1) (b<0)\n",
      "DAG[2]: 2.0%\n",
      "\tx0(2) <--- x0(1) (b>0)\n",
      "\tx0(2) <--- x1(1) (b>0)\n",
      "\tx0(2) <--- x2(1) (b<0)\n",
      "\tx0(2) <--- x3(1) (b<0)\n",
      "\tx0(2) <--- x4(1) (b<0)\n",
      "\tx1(2) <--- x2(2) (b<0)\n",
      "\tx1(2) <--- x4(2) (b>0)\n",
      "\tx1(2) <--- x0(1) (b<0)\n",
      "\tx1(2) <--- x1(1) (b<0)\n",
      "\tx1(2) <--- x2(1) (b>0)\n",
      "\tx1(2) <--- x3(1) (b>0)\n",
      "\tx1(2) <--- x4(1) (b>0)\n",
      "\tx2(2) <--- x0(1) (b>0)\n",
      "\tx2(2) <--- x1(1) (b>0)\n",
      "\tx2(2) <--- x2(1) (b<0)\n",
      "\tx2(2) <--- x3(1) (b>0)\n",
      "\tx2(2) <--- x4(1) (b<0)\n",
      "\tx3(2) <--- x2(2) (b<0)\n",
      "\tx3(2) <--- x0(1) (b>0)\n",
      "\tx3(2) <--- x1(1) (b<0)\n",
      "\tx3(2) <--- x2(1) (b>0)\n",
      "\tx3(2) <--- x3(1) (b>0)\n",
      "\tx3(2) <--- x4(1) (b<0)\n",
      "\tx4(2) <--- x0(2) (b>0)\n",
      "\tx4(2) <--- x0(1) (b<0)\n",
      "\tx4(2) <--- x1(1) (b>0)\n",
      "\tx4(2) <--- x2(1) (b<0)\n",
      "\tx4(2) <--- x3(1) (b>0)\n",
      "\tx4(2) <--- x4(1) (b<0)\n"
     ]
    }
   ],
   "source": [
    "t = 2\n",
    "labels = [f'x{i}({u})' for u in [t, t-1] for i in range(5)]\n",
    "print_dagc(dagc_list[t], 100, labels=labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Probability\n",
    "Using the `get_probabilities()` method, we can get the probability of bootstrapping."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[0.   0.37 0.1  0.12 1.  ]\n",
      "  [0.   0.   0.   1.   0.  ]\n",
      "  [0.   0.98 0.   0.5  0.24]\n",
      "  [0.   0.   0.   0.   0.  ]\n",
      "  [0.   1.   0.11 0.28 0.  ]]\n",
      "\n",
      " [[0.91 0.93 1.   0.94 0.97]\n",
      "  [0.99 0.99 1.   1.   0.94]\n",
      "  [1.   1.   1.   0.99 0.84]\n",
      "  [1.   0.98 1.   0.92 1.  ]\n",
      "  [0.98 1.   1.   1.   1.  ]]]\n"
     ]
    }
   ],
   "source": [
    "probs = result.get_probabilities(min_causal_effect=0.01)\n",
    "print(probs[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "B(1,1):\n",
      "[[0.   0.37 0.1  0.12 1.  ]\n",
      " [0.   0.   0.   1.   0.  ]\n",
      " [0.   0.98 0.   0.5  0.24]\n",
      " [0.   0.   0.   0.   0.  ]\n",
      " [0.   1.   0.11 0.28 0.  ]]\n",
      "B(1,0):\n",
      "[[0.91 0.93 1.   0.94 0.97]\n",
      " [0.99 0.99 1.   1.   0.94]\n",
      " [1.   1.   1.   0.99 0.84]\n",
      " [1.   0.98 1.   0.92 1.  ]\n",
      " [0.98 1.   1.   1.   1.  ]]\n",
      "B(2,2):\n",
      "[[0.   0.   0.   0.   0.  ]\n",
      " [0.06 0.   1.   0.04 1.  ]\n",
      " [0.8  0.   0.   0.   0.07]\n",
      " [0.03 0.02 1.   0.   0.1 ]\n",
      " [0.91 0.   0.   0.01 0.  ]]\n",
      "B(2,1):\n",
      "[[1.   1.   0.91 1.   0.92]\n",
      " [1.   0.86 1.   1.   0.96]\n",
      " [0.95 1.   0.95 1.   0.82]\n",
      " [1.   1.   1.   0.92 0.99]\n",
      " [1.   1.   0.97 1.   1.  ]]\n"
     ]
    }
   ],
   "source": [
    "t = 1\n",
    "print('B(1,1):')\n",
    "print(probs[t, 0])\n",
    "print('B(1,0):')\n",
    "print(probs[t, 1])\n",
    "\n",
    "t = 2\n",
    "print('B(2,2):')\n",
    "print(probs[t, 0])\n",
    "print('B(2,1):')\n",
    "print(probs[t, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Total Causal Effects\n",
    "Using the `get_total_causal_effects()` method, we can get the list of total causal effect. The total causal effects we can get are dictionary type variable.\n",
    "We can display the list nicely by assigning it to pandas.DataFrame. Also, we have replaced the variable index with a label below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>from</th>\n",
       "      <th>to</th>\n",
       "      <th>effect</th>\n",
       "      <th>probability</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>x1(1)</td>\n",
       "      <td>x0(1)</td>\n",
       "      <td>0.257084</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>x4(1)</td>\n",
       "      <td>x4(2)</td>\n",
       "      <td>-0.278507</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>x3(1)</td>\n",
       "      <td>x4(2)</td>\n",
       "      <td>0.185780</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>x1(1)</td>\n",
       "      <td>x4(2)</td>\n",
       "      <td>0.351397</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>x2(2)</td>\n",
       "      <td>x3(2)</td>\n",
       "      <td>-0.428210</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>x3(1)</td>\n",
       "      <td>x3(2)</td>\n",
       "      <td>-0.161284</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>x2(1)</td>\n",
       "      <td>x3(2)</td>\n",
       "      <td>0.495256</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>x1(1)</td>\n",
       "      <td>x3(2)</td>\n",
       "      <td>-0.579338</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>x0(1)</td>\n",
       "      <td>x3(2)</td>\n",
       "      <td>0.186140</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>x3(1)</td>\n",
       "      <td>x2(2)</td>\n",
       "      <td>0.400577</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>x1(1)</td>\n",
       "      <td>x2(2)</td>\n",
       "      <td>0.326661</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>x0(1)</td>\n",
       "      <td>x2(2)</td>\n",
       "      <td>0.161875</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>x2(2)</td>\n",
       "      <td>x1(2)</td>\n",
       "      <td>-0.692908</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>x0(1)</td>\n",
       "      <td>x1(2)</td>\n",
       "      <td>-0.563879</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>x4(2)</td>\n",
       "      <td>x1(2)</td>\n",
       "      <td>0.476373</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>x3(1)</td>\n",
       "      <td>x0(2)</td>\n",
       "      <td>-0.495518</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>x4(1)</td>\n",
       "      <td>x0(1)</td>\n",
       "      <td>-0.586968</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>x3(1)</td>\n",
       "      <td>x1(1)</td>\n",
       "      <td>0.388875</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>x0(1)</td>\n",
       "      <td>x0(2)</td>\n",
       "      <td>0.202197</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>x1(1)</td>\n",
       "      <td>x0(2)</td>\n",
       "      <td>0.191862</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>x1(1)</td>\n",
       "      <td>x4(1)</td>\n",
       "      <td>-0.356674</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>x1(1)</td>\n",
       "      <td>x2(1)</td>\n",
       "      <td>0.357268</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>x1(1)</td>\n",
       "      <td>x1(2)</td>\n",
       "      <td>-0.100172</td>\n",
       "      <td>0.99</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>x2(1)</td>\n",
       "      <td>x1(2)</td>\n",
       "      <td>0.169769</td>\n",
       "      <td>0.99</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>x3(1)</td>\n",
       "      <td>x4(1)</td>\n",
       "      <td>-0.108293</td>\n",
       "      <td>0.98</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>x4(1)</td>\n",
       "      <td>x3(2)</td>\n",
       "      <td>-0.158863</td>\n",
       "      <td>0.98</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26</th>\n",
       "      <td>x2(1)</td>\n",
       "      <td>x2(2)</td>\n",
       "      <td>-0.064596</td>\n",
       "      <td>0.97</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>27</th>\n",
       "      <td>x0(1)</td>\n",
       "      <td>x4(2)</td>\n",
       "      <td>-0.146124</td>\n",
       "      <td>0.97</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28</th>\n",
       "      <td>x3(1)</td>\n",
       "      <td>x0(1)</td>\n",
       "      <td>0.080405</td>\n",
       "      <td>0.97</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29</th>\n",
       "      <td>x3(1)</td>\n",
       "      <td>x2(1)</td>\n",
       "      <td>0.032170</td>\n",
       "      <td>0.94</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>30</th>\n",
       "      <td>x2(1)</td>\n",
       "      <td>x4(2)</td>\n",
       "      <td>-0.099157</td>\n",
       "      <td>0.94</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>31</th>\n",
       "      <td>x3(1)</td>\n",
       "      <td>x1(2)</td>\n",
       "      <td>0.079244</td>\n",
       "      <td>0.93</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32</th>\n",
       "      <td>x4(1)</td>\n",
       "      <td>x0(2)</td>\n",
       "      <td>-0.005440</td>\n",
       "      <td>0.92</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>33</th>\n",
       "      <td>x0(2)</td>\n",
       "      <td>x4(2)</td>\n",
       "      <td>0.261939</td>\n",
       "      <td>0.91</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>34</th>\n",
       "      <td>x2(1)</td>\n",
       "      <td>x0(2)</td>\n",
       "      <td>0.019144</td>\n",
       "      <td>0.91</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35</th>\n",
       "      <td>x0(2)</td>\n",
       "      <td>x1(2)</td>\n",
       "      <td>-0.029275</td>\n",
       "      <td>0.90</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>36</th>\n",
       "      <td>x4(1)</td>\n",
       "      <td>x1(2)</td>\n",
       "      <td>-0.014277</td>\n",
       "      <td>0.90</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>37</th>\n",
       "      <td>x4(1)</td>\n",
       "      <td>x2(2)</td>\n",
       "      <td>-0.019646</td>\n",
       "      <td>0.85</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>38</th>\n",
       "      <td>x0(2)</td>\n",
       "      <td>x3(2)</td>\n",
       "      <td>-0.106739</td>\n",
       "      <td>0.84</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>39</th>\n",
       "      <td>x0(2)</td>\n",
       "      <td>x2(2)</td>\n",
       "      <td>0.250640</td>\n",
       "      <td>0.80</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>40</th>\n",
       "      <td>x4(1)</td>\n",
       "      <td>x2(1)</td>\n",
       "      <td>-0.169832</td>\n",
       "      <td>0.24</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>41</th>\n",
       "      <td>x2(1)</td>\n",
       "      <td>x0(1)</td>\n",
       "      <td>0.015604</td>\n",
       "      <td>0.20</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>42</th>\n",
       "      <td>x4(2)</td>\n",
       "      <td>x3(2)</td>\n",
       "      <td>-0.147539</td>\n",
       "      <td>0.18</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>43</th>\n",
       "      <td>x2(1)</td>\n",
       "      <td>x4(1)</td>\n",
       "      <td>-0.171814</td>\n",
       "      <td>0.11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>44</th>\n",
       "      <td>x4(2)</td>\n",
       "      <td>x2(2)</td>\n",
       "      <td>0.155502</td>\n",
       "      <td>0.07</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>45</th>\n",
       "      <td>x3(2)</td>\n",
       "      <td>x1(2)</td>\n",
       "      <td>-0.155433</td>\n",
       "      <td>0.05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>46</th>\n",
       "      <td>x1(2)</td>\n",
       "      <td>x3(2)</td>\n",
       "      <td>-0.174134</td>\n",
       "      <td>0.02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>47</th>\n",
       "      <td>x2(2)</td>\n",
       "      <td>x4(2)</td>\n",
       "      <td>0.045734</td>\n",
       "      <td>0.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>48</th>\n",
       "      <td>x3(2)</td>\n",
       "      <td>x4(2)</td>\n",
       "      <td>-0.146344</td>\n",
       "      <td>0.01</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     from     to    effect  probability\n",
       "0   x1(1)  x0(1)  0.257084         1.00\n",
       "1   x4(1)  x4(2) -0.278507         1.00\n",
       "2   x3(1)  x4(2)  0.185780         1.00\n",
       "3   x1(1)  x4(2)  0.351397         1.00\n",
       "4   x2(2)  x3(2) -0.428210         1.00\n",
       "5   x3(1)  x3(2) -0.161284         1.00\n",
       "6   x2(1)  x3(2)  0.495256         1.00\n",
       "7   x1(1)  x3(2) -0.579338         1.00\n",
       "8   x0(1)  x3(2)  0.186140         1.00\n",
       "9   x3(1)  x2(2)  0.400577         1.00\n",
       "10  x1(1)  x2(2)  0.326661         1.00\n",
       "11  x0(1)  x2(2)  0.161875         1.00\n",
       "12  x2(2)  x1(2) -0.692908         1.00\n",
       "13  x0(1)  x1(2) -0.563879         1.00\n",
       "14  x4(2)  x1(2)  0.476373         1.00\n",
       "15  x3(1)  x0(2) -0.495518         1.00\n",
       "16  x4(1)  x0(1) -0.586968         1.00\n",
       "17  x3(1)  x1(1)  0.388875         1.00\n",
       "18  x0(1)  x0(2)  0.202197         1.00\n",
       "19  x1(1)  x0(2)  0.191862         1.00\n",
       "20  x1(1)  x4(1) -0.356674         1.00\n",
       "21  x1(1)  x2(1)  0.357268         1.00\n",
       "22  x1(1)  x1(2) -0.100172         0.99\n",
       "23  x2(1)  x1(2)  0.169769         0.99\n",
       "24  x3(1)  x4(1) -0.108293         0.98\n",
       "25  x4(1)  x3(2) -0.158863         0.98\n",
       "26  x2(1)  x2(2) -0.064596         0.97\n",
       "27  x0(1)  x4(2) -0.146124         0.97\n",
       "28  x3(1)  x0(1)  0.080405         0.97\n",
       "29  x3(1)  x2(1)  0.032170         0.94\n",
       "30  x2(1)  x4(2) -0.099157         0.94\n",
       "31  x3(1)  x1(2)  0.079244         0.93\n",
       "32  x4(1)  x0(2) -0.005440         0.92\n",
       "33  x0(2)  x4(2)  0.261939         0.91\n",
       "34  x2(1)  x0(2)  0.019144         0.91\n",
       "35  x0(2)  x1(2) -0.029275         0.90\n",
       "36  x4(1)  x1(2) -0.014277         0.90\n",
       "37  x4(1)  x2(2) -0.019646         0.85\n",
       "38  x0(2)  x3(2) -0.106739         0.84\n",
       "39  x0(2)  x2(2)  0.250640         0.80\n",
       "40  x4(1)  x2(1) -0.169832         0.24\n",
       "41  x2(1)  x0(1)  0.015604         0.20\n",
       "42  x4(2)  x3(2) -0.147539         0.18\n",
       "43  x2(1)  x4(1) -0.171814         0.11\n",
       "44  x4(2)  x2(2)  0.155502         0.07\n",
       "45  x3(2)  x1(2) -0.155433         0.05\n",
       "46  x1(2)  x3(2) -0.174134         0.02\n",
       "47  x2(2)  x4(2)  0.045734         0.01\n",
       "48  x3(2)  x4(2) -0.146344         0.01"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "causal_effects = result.get_total_causal_effects(min_causal_effect=0.01)\n",
    "\n",
    "df = pd.DataFrame(causal_effects)\n",
    "\n",
    "labels = [f'x{i}({t})' for t in range(3) for i in range(5)]\n",
    "df['from'] = df['from'].apply(lambda x : labels[x])\n",
    "df['to'] = df['to'].apply(lambda x : labels[x])\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can easily perform sorting operations with pandas.DataFrame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>from</th>\n",
       "      <th>to</th>\n",
       "      <th>effect</th>\n",
       "      <th>probability</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>x2(1)</td>\n",
       "      <td>x3(2)</td>\n",
       "      <td>0.495256</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>x4(2)</td>\n",
       "      <td>x1(2)</td>\n",
       "      <td>0.476373</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>x3(1)</td>\n",
       "      <td>x2(2)</td>\n",
       "      <td>0.400577</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>x3(1)</td>\n",
       "      <td>x1(1)</td>\n",
       "      <td>0.388875</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>x1(1)</td>\n",
       "      <td>x2(1)</td>\n",
       "      <td>0.357268</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     from     to    effect  probability\n",
       "6   x2(1)  x3(2)  0.495256          1.0\n",
       "14  x4(2)  x1(2)  0.476373          1.0\n",
       "9   x3(1)  x2(2)  0.400577          1.0\n",
       "17  x3(1)  x1(1)  0.388875          1.0\n",
       "21  x1(1)  x2(1)  0.357268          1.0"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.sort_values('effect', ascending=False).head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And with pandas.DataFrame, we can easily filter by keywords. The following code extracts the causal direction towards x0(2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>from</th>\n",
       "      <th>to</th>\n",
       "      <th>effect</th>\n",
       "      <th>probability</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>x3(1)</td>\n",
       "      <td>x0(2)</td>\n",
       "      <td>-0.495518</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>x0(1)</td>\n",
       "      <td>x0(2)</td>\n",
       "      <td>0.202197</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>x1(1)</td>\n",
       "      <td>x0(2)</td>\n",
       "      <td>0.191862</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32</th>\n",
       "      <td>x4(1)</td>\n",
       "      <td>x0(2)</td>\n",
       "      <td>-0.005440</td>\n",
       "      <td>0.92</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>34</th>\n",
       "      <td>x2(1)</td>\n",
       "      <td>x0(2)</td>\n",
       "      <td>0.019144</td>\n",
       "      <td>0.91</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     from     to    effect  probability\n",
       "15  x3(1)  x0(2) -0.495518         1.00\n",
       "18  x0(1)  x0(2)  0.202197         1.00\n",
       "19  x1(1)  x0(2)  0.191862         1.00\n",
       "32  x4(1)  x0(2) -0.005440         0.92\n",
       "34  x2(1)  x0(2)  0.019144         0.91"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df[df['to']=='x0(2)'].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because it holds the raw data of the total causal effect (the original data for calculating the median), it is possible to draw a histogram of the values of the causal effect, as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 6.,  8., 13., 12., 19., 16.,  9., 10.,  6.,  1.]),\n",
       " array([0.026, 0.057, 0.088, 0.119, 0.15 , 0.181, 0.212, 0.243, 0.274,\n",
       "        0.304, 0.335]),\n",
       " <BarContainer object of 10 artists>)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAD7CAYAAAB0d9PAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8WgzjOAAAACXBIWXMAAAsTAAALEwEAmpwYAAATaElEQVR4nO3de4xedZ3H8XeptKAtus6ON25dpf2GKG4FiXGjazZBUbKlKLJiFHBxBcXYZKVZo2aRuIEQS7OIYCAa5bKGNeimpUjwkoihXrISbQTNfq2itLYQx1HT1oVeZmb/mDP4MDvz6zNzznObvl9J0/Oc83vO7/ebM+f5zHMuv7NoYmICSZJmc1SvGyBJ6m8GhSSpyKCQJBUZFJKkIoNCklT0rF43oA1LgTOBx4GxHrdFkgbFYuDFwA+B/XVWNAhBcSbwYK8bIUkD6vXA1jorGISgeBzgD3/4E+Pj/XfPx9DQMkZH9/W6GbXYh/6wEPoAC6MfC6UPlcfrrmsQgmIMYHx8oi+DAujbds2FfegPC6EPsDD6sRD6UKl9yN6T2ZKkIoNCklRkUEiSigwKSVKRQSFJKjIoJElFBoUkqWgQ7qOQDuvAwTGGh5d3vd6n9h9i754nu16v1E0GhRaEJUcvZs2Vm7te75aNa9nb9Vql7vLQkySpyKCQJBUZFJKkIoNCklRkUEiSigwKSVKRQSFJKjIoJElFBoUkqcigkCQVGRSSpKK2xnqKiOuB84EVwGmZ+UhErAA2tRR7HnBcZj5/hvdfDVwB7K5mfTczPzjfRkuSuqfdQQE3AZ8GHpyakZm/BlZPvY6IGw6zvjsyc/1cGyhJ6q22giIztwJExIzLI2IJ8C7g7MZaJknqC00NM34usCszf1Qoc2FEvAl4AvhEZn6/obolSR3UVFBcCnyhsPwW4JrMPBgRbwQ2R8SpmTnabgVDQ8vqtrFjevHAnKYthD70SpM/u4WyHRZCPxZCH5pSOygi4njgDcBFs5XJzCdapr8ZETuBVwDfabee0dF9jI9P1GlqRwwPL2dkZLAfXbNQ+tArTf3sFsJ2gIXRj4XSh6Y0cXnsJcDXSt8OqjCZml7N5NVT2UDdkqQOa/fy2BuBtwEvAr4VEaOZ+fJq8XuAdTO85z7gqsx8CLg2Is4AxoADwEWt3zIkSf2r3aue1jFDGFTLVs0y/5yW6Uvm1TpJUs95Z7YkqcigkCQVGRSSpCKDQpJUZFBIkooMCklSkUEhSSoyKCRJRQaFJKnIoJAkFRkUkqSipp5HIR2RDhwc69nzKJ7af4i9e55srG5pNgaFVMOSoxez5srNPal7y8a1DPYTEzQoPPQkSSoyKCRJRQaFJKnIoJAkFRkUkqSidp+ZfT1wPrACOC0zH6nm/xp4qvoH8JHM/PoM73828EXgDOAQsD4z763beElS57V7eewm4NPAgzMse/tUcBSsB/Zk5ikRsRJ4MCJOycx97TdVktQLbR16ysytmbmzRj3vAG6t1rUdeAh4S431SZK6pIkb7r4UEYuArcDHMvOPM5Q5CXis5fUO4MQG6pYkdVjdoHh9Zu6MiKXADcBNwLtrt2oGQ0PLOrHaRjQ5hEOvLIQ+HIn6dbv1a7vmYiH0oSm1gmLqcFRm7o+IzwL3zFJ0B3AyMFK9Pgn49lzqGh3dx/j4xHyb2jHDw8sZGRnsgRQWSh+ORP243RbK79NC6ENT5n15bEQ8JyKeW00vAi4Ets1S/G7g8qrsSuBM4P751i1J6p52L4+9EXgb8CLgWxExCqwBvhoRi4HFwM+AK1resw04JzN3AxuA2yLiF8AYcFlmDnZcS9IRoq2gyMx1wLoZFr2q8J7VLdN/Ai6Ya+MkSb3nndmSpCKDQpJUZFBIkooMCklSkUEhSSoyKCRJRQaFJKnIoJAkFRkUkqQig0KSVGRQSJKKDApJUpFBIUkqMigkSUUGhSSpyKCQJBUZFJKkoraecKfBsvy4Yzlm6dw2bVMPYn9q/yH27nmykXVJ6g/tPjP7euB8YAVwWmY+EhFDwJ3Ay4ADwHbg8swcmeH9twFnAb+rZt2dmdfUbr1mdMzSZ7Hmys09qXvLxrX4MHRpYWn3z85NwKeBB1vmTQCfyswHACJiA3Ad8N5Z1nFdZt40v2ZKknqlraDIzK0AEdE67/fAAy3FfgB8oMG2SZL6QCMnsyPiKCZD4p5CsQ9HxMMRsSkiTm2iXklS5zV1MvszwD5gtkNLHwcez8zxiLgYuD8iXpqZY+1WMDS0rIFmdkZTJ4IXCn8e3dOvP+t+bddcLIQ+NKV2UFQnulcCazJzfKYymbmrZfqOiPh34ATgsXbrGR3dx/j4RN3mNm54eDkjI/11+rbXv+C9+Hn0us+90m+/e9Cf+8RcLZQ+NKXWoaeIuBY4AzgvM/cXyh3fMn02MAbsmq28JKl/tHt57I3A24AXAd+KiFHgH4CPAj8Hvled6P5VZr61es824JzM3A3cHhEvBMaBPcC5mXmo4b5Ikjqg3aue1gHrZli0qPCe1S3TZ825ZZKkvuAQHpKkIofwUKMOHBw7Yk8sSwuVQaFGLTl6cU+GD9mycW3X65SOFB56kiQVGRSSpCKDQpJUZFBIkooMCklSkUEhSSoyKCRJRQaFJKnIoJAkFRkUkqQig0KSVGRQSJKKDApJUpFBIUkqMigkSUUGhSSp6LAPLoqI64HzgRXAaZn5SDV/FXA7MASMAhdn5vYZ3r8YuBF4MzABXJeZn2+qA5KkzmrnG8Um4G+Bx6bNvwW4OTNXATcDt87y/ncBpwArgdcCV0fEivk0VpLUfYcNiszcmpk7W+dFxAuA04G7qll3AadHxPAMq3gH8LnMHM/MESaD54JarZYkdc18n5l9IrArM8cAMnMsInZX80emlT2JZ34b2VGVm5OhoWXzbGrnDQ8v73UTdITqxe/egYNjLDl6cbFMp9rVTt1Ncb/+s/kGRdeNju5jfHyi1834f4aHlzMysrfXzXgGf8GPHL343RseXs6aKzd3vV6ALRvXdqXP/bhfz1WTnwPzveppJ3B8daJ66oT1S6r50+0ATm55fdIs5SRJfWheQZGZvwW2Ae+sZr0T+HF1DmK6u4H3RcRR1TmM84CvzKdeSVL3HTYoIuLGiPgNcALwrYj4abXo/cCHIuLnwIeq11PvuS8iXl29vBN4FNgO/AD4ZGb+qsE+SJI66LDnKDJzHbBuhvn/A7xmlvec0zI9BnygRhslST3kndmSpCKDQpJUZFBIkooMCklSkUEhSSoamDuzJT3TgYNj3oWvrjAopAG15OjFPRlKY8vGtV2vU73loSdJUpFBIUkqMigkSUUGhSSpyKCQJBUZFJKkIoNCklRkUEiSigwKSVKRd2Z30PLjjuWYpf6IJQ02P8U66Jilz3KIBUkDr1ZQRMQKYFPLrOcBx2Xm86eVuxq4AthdzfpuZn6wTt2SpO6oFRSZ+Wtg9dTriLihsM47MnN9nfokSd3X2KGniFgCvAs4u6l1SpJ6r8lzFOcCuzLzR7MsvzAi3gQ8AXwiM7/fYN2SpA5pMiguBb4wy7JbgGsy82BEvBHYHBGnZuZouysfGlrWRBs7wofHSN3Trf3N/frPGgmKiDgeeANw0UzLM/OJlulvRsRO4BXAd9qtY3R0H+PjE3Wb2rjh4eWMjOyddZmkZs22vzWptF8PiiY/f5q64e4S4GuzfUOogmRqejWwAsiG6pYkdVBTh57eA6xrnRER9wFXZeZDwLURcQYwBhwALmr9liFJ6l+NBEVmrpph3jkt05c0UY8kqfsc60mSVGRQSJKKDApJUpFBIUkqMigkSUUGhSSpyKCQJBUZFJKkIoNCklRkUEiSigwKSVJRk8+j6EvLjzuWY5Z2tpsOJy51x4GDYz15HsVT+w+xd8+TXam3Hy34oDhm6bNYc+XmntS9ZePantQrLVRLjl7ck/15y8a1DPbTKerx0JMkqcigkCQVGRSSpCKDQpJUZFBIkopqX/UUEb8Gnqr+AXwkM78+rcyzgS8CZwCHgPWZeW/duiVJndfU5bFvz8xHCsvXA3sy85SIWAk8GBGnZOa+huqXJHVItw49vQO4FSAztwMPAW/pUt2SpBqa+kbxpYhYBGwFPpaZf5y2/CTgsZbXO4ATG6pbktRBTQTF6zNzZ0QsBW4AbgLe3cB6n2FoaFnTq5Skth3JQ/XUDorM3Fn9vz8iPgvcM0OxHcDJwEj1+iTg23OpZ3R0H+PjE3Nu35G8cSU1Z2RksAbxaPKzr9Y5ioh4TkQ8t5peBFwIbJuh6N3A5VW5lcCZwP116pYkdUfdk9kvBB6IiJ8AjwCrgCsAImJbRLykKrcBeF5E/AK4F7gsMwcrniXpCFXr0FNmPgq8apZlq1um/wRcUKcuSVJveGe2JKnIoJAkFRkUkqQig0KSVGRQSJKKDApJUpFBIUkqMigkSUUGhSSpyKCQJBUZFJKkIoNCklRkUEiSigwKSVKRQSFJKjIoJElFBoUkqcigkCQV1XoUakQMAXcCLwMOANuByzNzZFq524CzgN9Vs+7OzGvq1C1J6o5aQQFMAJ/KzAcAImIDcB3w3hnKXpeZN9WsT5LUZbWCIjN/DzzQMusHwAfqrFOS1F8aO0cREUcxGRL3zFLkwxHxcERsiohTm6pXktRZdQ89tfoMsA+Y6fDSx4HHM3M8Ii4G7o+Il2bmWLsrHxpa1lAzJWnuhoeX97oJPdNIUETE9cBKYE1mjk9fnpm7WqbviIh/B04AHmu3jtHRfYyPT8y5bUfyxpXUnJGRvb1uwpw0+dlX+9BTRFwLnAGcl5n7ZylzfMv02cAYsGumspKk/lL38tiXAx8Ffg58LyIAfpWZb42IbcA5mbkbuD0iXgiMA3uAczPzUK2WS5K6ou5VTz8FFs2ybHXL9Fl16pEk9Y53ZkuSigwKSVKRQSFJKjIoJElFBoUkqcigkCQVGRSSpCKDQpJUZFBIkooMCklSkUEhSSpq8nkUkrQgHTg41rNHFjy1/xB79zzZk7qnGBSSdBhLjl7Mmis396TuLRvX0usnYXjoSZJUZFBIkooMCklSkUEhSSoyKCRJRQaFJKmo9uWxEbEKuB0YAkaBizNz+7Qyi4EbgTcDE8B1mfn5unVLkjqviW8UtwA3Z+Yq4Gbg1hnKvAs4BVgJvBa4OiJWNFC3JKnDan2jiIgXAKcDb6xm3QXcFBHDmTnSUvQdwOcycxwYiYhNwAXAhjaqWQxw1FGL5t3OF/zFsfN+b129qts+L/x6e1m3fe6uOp9/VJ+hdSyamJiY95sj4gzgjsx8ecu8nwHvzswftcx7GLg0M39Yvf4X4ITMXNdGNa8DHpx3IyXpyPZ6YGudFQzCEB4/ZLKjjwNjPW6LJA2KxcCLmfwMraVuUOwEjo+IxZk5Vp20fkk1v9UO4GT+3OCTgMfarGM/NdNQko5Qv2xiJbVOZmfmb4FtwDurWe8Efjzt/ATA3cD7IuKoiBgGzgO+UqduSVJ3NHHV0/uBD0XEz4EPVa+JiPsi4tVVmTuBR4HtwA+AT2bmrxqoW5LUYbVOZkuSFj7vzJYkFRkUkqQig0KSVGRQSJKKBuGGu56oO9hhRFwNXAHsrop/NzM/2J3WP92+dvrwJuBa4DTgM5m5vmVZzwdzbKAPVzMY2+FfgQuZvKn0IPCxzPx6tezZwBeBM4BDwPrMvLd7PWikD7cBZwG/q4rfnZnXdKf1T7evnT78I/DPwDiTN6x9LjNvrJb1fH+o2lG3H1czx33CbxSza2Kwwzsyc3X1r6sfTpV2+vAo8E/MPO5WPwzmWLcPMBjb4b+BMzPzlcClwJcjYmpwofXAnsw8BVgDfD4ilnWh3a3q9gEmP1intkNXQ6LSTh++Cvx1Zq4G/ga4MiJeWS3rh/0B6vcD5rhPGBQzaBns8K5q1l3A6dXNgq2eHuywuslwE5ODHfZcu33IzF9k5jYm/1Kdrqf9a6gPPTWHPnw9M/+3evkTYBGTfzHC5Ha4tSq3HXgIeEuHm/60hvrQU3Pow57MnLpn4NnA0Ux+e4A+2N8b6secGRQzOxHYlZljANX/u6v5raYPRbJjWpkLI+InEfGNiHhtJxs8g3b7UHK4/nVaE32AwdsOFwO/zMzfVK8HcTtM7wPAhyPi4YjYFBGndq65M2q7DxFxbkT8lMmf+YbMfLha1OvtAM30A+a4TxgUnXML8FfV1/ANwOaI6Iu/ro4wA7UdIuINwL/x52FxBs4sffg4cEpmngb8F3B/dcy/72TmPdWI2KuAiyIiet2m+Sj0Y877hEExs6cHO4SnT2KVBjucctJUmcx8IjMPVtPfrOa/osPtbtVuH0pm7V+X1O7DIG2H6i+7/wDOy8xsWTQw22G2PmTmrup5NGTmHcAy4IQutH3KnH+XMnMHk+dd/r6a1evtAA30Yz77hEExgyYGO4yI46cKRcRqYAWQdMkc+lDS08Ecm+jDoGyHiDgT+DLw9tZnuVTuBi6vyq0EzgTu72Czn6GJPkzbDmczeWXUrg42+xnm0IdTW6b/Evg7YOqQTc8HN22iH/PZJ7w8dnbvB26PiKuAPzB5zJWIuA+4KjMfYnKww9cwOdghPHOww2urBzuNAQeAizLziW52gDb6EBGvA/4TOA5YFBEXAu+tLmss9W9Q+jAQ2wH4LHAscGvLkY6LquPKG4DbIuIXTPbjsszcO2B9uD0iXsjk5Zp7gHMzs9sXH7TTh8uqy60PMnky/qbM/Eb1/n7YH6B+P+a8TzgooCSpyENPkqQig0KSVGRQSJKKDApJUpFBIUkqMigkSUUGhSSpyKCQJBX9H4SZUTzWovMhAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set()\n",
    "%matplotlib inline\n",
    "\n",
    "from_index = 5 # index of x0(1). (index:0)+(n_features:5)*(timepoint:1) = 5\n",
    "to_index = 12 # index of x2(2). (index:2)+(n_features:5)*(timepoint:2) = 12\n",
    "plt.hist(result.total_effects_[:, to_index, from_index])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}