IDC_MarkovianInput.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyPAH7jrHC03YDlAPhT2cZRK",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/aliash98/BPM_CausalReasoning/blob/main/IDC_MarkovianInput.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "!pip install ylearn"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "AbUoElPOIrcT",
        "outputId": "9096e6b5-0e72-4970-fe67-0442a502eebe"
      },
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/\n",
            "Requirement already satisfied: ylearn in /usr/local/lib/python3.9/dist-packages (0.2.0)\n",
            "Requirement already satisfied: networkx in /usr/local/lib/python3.9/dist-packages (from ylearn) (2.8)\n",
            "Requirement already satisfied: pandas>=0.25.3 in /usr/local/lib/python3.9/dist-packages (from ylearn) (1.4.4)\n",
            "Requirement already satisfied: scikit-learn>=0.22.1 in /usr/local/lib/python3.9/dist-packages (from ylearn) (1.2.2)\n",
            "Requirement already satisfied: torch in /usr/local/lib/python3.9/dist-packages (from ylearn) (1.13.1+cu116)\n",
            "Requirement already satisfied: numpy>=1.16.5 in /usr/local/lib/python3.9/dist-packages (from ylearn) (1.22.4)\n",
            "Requirement already satisfied: tqdm in /usr/local/lib/python3.9/dist-packages (from ylearn) (4.65.0)\n",
            "Requirement already satisfied: joblib in /usr/local/lib/python3.9/dist-packages (from ylearn) (1.1.1)\n",
            "Requirement already satisfied: scipy in /usr/local/lib/python3.9/dist-packages (from ylearn) (1.10.1)\n",
            "Requirement already satisfied: pydot in /usr/local/lib/python3.9/dist-packages (from ylearn) (1.3.0)\n",
            "Requirement already satisfied: ipython in /usr/local/lib/python3.9/dist-packages (from ylearn) (7.9.0)\n",
            "Requirement already satisfied: python-dateutil>=2.8.1 in /usr/local/lib/python3.9/dist-packages (from pandas>=0.25.3->ylearn) (2.8.2)\n",
            "Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.9/dist-packages (from pandas>=0.25.3->ylearn) (2022.7.1)\n",
            "Requirement already satisfied: threadpoolctl>=2.0.0 in /usr/local/lib/python3.9/dist-packages (from scikit-learn>=0.22.1->ylearn) (3.1.0)\n",
            "Requirement already satisfied: traitlets>=4.2 in /usr/local/lib/python3.9/dist-packages (from ipython->ylearn) (5.7.1)\n",
            "Requirement already satisfied: prompt-toolkit<2.1.0,>=2.0.0 in /usr/local/lib/python3.9/dist-packages (from ipython->ylearn) (2.0.10)\n",
            "Requirement already satisfied: pickleshare in /usr/local/lib/python3.9/dist-packages (from ipython->ylearn) (0.7.5)\n",
            "Requirement already satisfied: decorator in /usr/local/lib/python3.9/dist-packages (from ipython->ylearn) (4.4.2)\n",
            "Requirement already satisfied: pygments in /usr/local/lib/python3.9/dist-packages (from ipython->ylearn) (2.6.1)\n",
            "Requirement already satisfied: pexpect in /usr/local/lib/python3.9/dist-packages (from ipython->ylearn) (4.8.0)\n",
            "Requirement already satisfied: backcall in /usr/local/lib/python3.9/dist-packages (from ipython->ylearn) (0.2.0)\n",
            "Requirement already satisfied: jedi>=0.10 in /usr/local/lib/python3.9/dist-packages (from ipython->ylearn) (0.18.2)\n",
            "Requirement already satisfied: setuptools>=18.5 in /usr/local/lib/python3.9/dist-packages (from ipython->ylearn) (63.4.3)\n",
            "Requirement already satisfied: pyparsing>=2.1.4 in /usr/local/lib/python3.9/dist-packages (from pydot->ylearn) (3.0.9)\n",
            "Requirement already satisfied: typing-extensions in /usr/local/lib/python3.9/dist-packages (from torch->ylearn) (4.5.0)\n",
            "Requirement already satisfied: parso<0.9.0,>=0.8.0 in /usr/local/lib/python3.9/dist-packages (from jedi>=0.10->ipython->ylearn) (0.8.3)\n",
            "Requirement already satisfied: six>=1.9.0 in /usr/local/lib/python3.9/dist-packages (from prompt-toolkit<2.1.0,>=2.0.0->ipython->ylearn) (1.15.0)\n",
            "Requirement already satisfied: wcwidth in /usr/local/lib/python3.9/dist-packages (from prompt-toolkit<2.1.0,>=2.0.0->ipython->ylearn) (0.2.6)\n",
            "Requirement already satisfied: ptyprocess>=0.5 in /usr/local/lib/python3.9/dist-packages (from pexpect->ipython->ylearn) (0.7.0)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "!pip install networkx==2.8"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "PVg8YdFKgBqj",
        "outputId": "06bfcd2d-cc87-4451-ce6e-8bd0ba6eb829"
      },
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/\n",
            "Requirement already satisfied: networkx==2.8 in /usr/local/lib/python3.9/dist-packages (2.8)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import networkx as nx\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "from ylearn.causal_model.graph import CausalGraph\n",
        "from ylearn.causal_model import CausalModel\n",
        "import sys\n",
        "from itertools import chain, combinations"
      ],
      "metadata": {
        "id": "-K2RFYL6Isk7"
      },
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# This function is modified for the IDC algorithm\n",
        "\n",
        "def markov_to_semimarkov(G):\n",
        "\n",
        "\n",
        "  G_sm = nx.DiGraph()\n",
        "\n",
        "  # Rule 1 ------- Adding all observed nodes to the new graph\n",
        "  for node in G.nodes(data=True):\n",
        "    if node[1]['color'] == 'blue':\n",
        "      G_sm.add_node(node[0])\n",
        "  \n",
        "  # Rule 2 ------- Importing edges between two blue nodes\n",
        "  for edge in G.edges(data=True):\n",
        "    from_node = edge[0]\n",
        "    to_node = edge[1]\n",
        "    from_node_color = G.nodes[from_node]['color']\n",
        "    to_node_color = G.nodes[to_node]['color']\n",
        "    if from_node_color =='blue' and to_node_color =='blue':\n",
        "      G_sm.add_edge(edge[0], edge[1])\n",
        "\n",
        "  # Rule 3 ------- Adding edges for directed paths through grey nodes between two blue nodes\n",
        "  grey_nodes = [n for n,v in G.nodes(data=True) if v['color'] == 'grey']  \n",
        "  blue_nodes = [n for n,v in G.nodes(data=True) if v['color'] == 'blue']\n",
        "\n",
        "  for source in blue_nodes:\n",
        "    for sink in blue_nodes:\n",
        "      if source == sink or G.has_edge(source, sink):\n",
        "        continue\n",
        "      for path in nx.all_simple_paths(G, source=source, target=sink):\n",
        "        counter = 0\n",
        "        for between_node in path:\n",
        "          if between_node in blue_nodes:\n",
        "            counter += 1\n",
        "        if counter == 2:\n",
        "          G_sm.add_edge(source, sink)\n",
        "          break\n",
        "    \n",
        "  # Rule 4 ------- Extracting the hidden confounders\n",
        "  confounding_list = []\n",
        "\n",
        "  for source in grey_nodes:\n",
        "    potential_confounding = []\n",
        "    for sink in blue_nodes:\n",
        "      for path in nx.all_simple_paths(G, source=source, target=sink):\n",
        "        red_flag = 0\n",
        "        for between_node in path:\n",
        "          if between_node in blue_nodes and between_node != sink:\n",
        "            red_flag = 1\n",
        "            break\n",
        "        if red_flag == 0:\n",
        "          potential_confounding.append(sink)\n",
        "          break\n",
        "    if len(potential_confounding) > 1:\n",
        "      confounding_list.append(potential_confounding)\n",
        "    \n",
        "  # we construct a graph for confounding arcs\n",
        "  G_confounded = nx.Graph()\n",
        "  number = 1\n",
        "  hidden_vars = []\n",
        "  for item in confounding_list:\n",
        "    name = \"U\" + str(number)\n",
        "    G_sm.add_node(name)\n",
        "    for node in item:\n",
        "      G_sm.add_edge(name, node)\n",
        "    hidden_vars.append(name)\n",
        "    number += 1\n",
        "\n",
        "  return G_sm, hidden_vars\n",
        "  # returns graph + list of hidden variables"
      ],
      "metadata": {
        "id": "G_269SkFJG-o"
      },
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# function to compute the power set \n",
        "def powerset(iterable):\n",
        "    s = list(iterable)\n",
        "    return list(chain.from_iterable(combinations(s, r) for r in range(1, len(s)+1)))"
      ],
      "metadata": {
        "id": "orvE2q_XJPBm"
      },
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "def remove_redundant_sets(s):\n",
        "  final_set = []\n",
        "  for a in s:\n",
        "    counter = 0\n",
        "    for b in s:\n",
        "      if set(b).issubset(a):\n",
        "        counter += 1\n",
        "    if counter == 1 :\n",
        "      final_set.append(a)\n",
        "  return final_set"
      ],
      "metadata": {
        "id": "YoGnH1KOE8zc"
      },
      "execution_count": 6,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "id": "6xM4oJIlfFDj"
      },
      "outputs": [],
      "source": [
        "# IDC - ID implementation by @mattez90\n",
        "\n",
        "def condset(var,orderlist):\n",
        "    pos = int(np.where(orderlist==var)[0])\n",
        "    if pos==0:\n",
        "        condset = np.array([])\n",
        "    else:\n",
        "        # plus one in order to get the node itself\n",
        "        condset = orderlist[0:pos]\n",
        "    return condset   \n",
        "\n",
        "\n",
        "def toposort(varsarray, order):\n",
        "    return np.array([orderi for orderi in order if orderi in varsarray])\n",
        "\n",
        "\n",
        "def removeedge(x,G,hiddvar):\n",
        "    G_del = G.copy()\n",
        "    for edge in x:\n",
        "        neighbors = list(G.predecessors(edge)) + list(G.successors(edge))\n",
        "         #Get the hidden edge\n",
        "        delete = [u for u in hiddvar.tolist() if u in neighbors] + [edge]\n",
        "        #Remove node and hidden edge connected to the removed edge\n",
        "        G_del.remove_nodes_from(delete)\n",
        " \n",
        "    return G_del\n",
        "\n",
        "def incomingedgedel(nodes,G):\n",
        "  \n",
        "    #Retreive the adjency matrix\n",
        "    adjmtx = nx.to_numpy_matrix(G)\n",
        "    for nodei in nodes.tolist():\n",
        "        adjmtx[:,list(nx.nodes(G)).index(nodei)]=0\n",
        "    #Convert again to a graph object\n",
        "    G_bar_x = nx.from_numpy_matrix(adjmtx,create_using=nx.DiGraph)\n",
        "    #Map edge according to the list in the previous graph\n",
        "    G_ = nx.relabel_nodes(G_bar_x, {i:j for i,j in enumerate(list(nx.nodes(G)))})\n",
        "      \n",
        "    return G_\n",
        "\n",
        "def outcomeedgedel(nodes,G):\n",
        "\n",
        "    #Retreive the adjency matrix\n",
        "    adjmtx = nx.to_numpy_matrix(G)\n",
        "    for nodei in nodes.tolist():\n",
        "        adjmtx[list(nx.nodes(G)).index(nodei),:]=0\n",
        "    #Convert again to a graph object\n",
        "    G_bar_x = nx.from_numpy_matrix(adjmtx,create_using=nx.DiGraph)\n",
        "    #Map edge according to the list in the previous graph\n",
        "    G_ = nx.relabel_nodes(G_bar_x, {i:j for i,j in enumerate(list(nx.nodes(G)))})\n",
        "      \n",
        "    return G_\n",
        "\n",
        "\n",
        "def ispresent(x,graph):\n",
        "    return np.array([node for node in x.tolist() if node  in list(nx.nodes(graph))])\n",
        "\n",
        "\n",
        "def C_components(graph,hiddvar,orderlist):\n",
        " \n",
        "    adj_mtx  = np.zeros((len(list(nx.nodes(graph))),len(list(nx.nodes(graph)))))\n",
        "    listnode = list(nx.nodes(graph))\n",
        "    #Trace all the bi-derected path starting from the empty matrix\n",
        "    for u in hiddvar:\n",
        "        if u in listnode:\n",
        "            neighbors = list(graph.neighbors(u))\n",
        "            adj_mtx[listnode.index(neighbors[0]),listnode.index(neighbors[1])] = adj_mtx[listnode.index(neighbors[0]),listnode.index(neighbors[1])] + 1\n",
        "            adj_mtx[listnode.index(neighbors[1]),listnode.index(neighbors[0])] = adj_mtx[listnode.index(neighbors[1]),listnode.index(neighbors[0])] + 1\n",
        " \n",
        "    # Delete nn-observable node from list of nodes and row and column from the\n",
        "    # adjacency matrix\n",
        "    adj_mtx2 = np.delete(np.delete(adj_mtx,[listnode.index(u) for u in hiddvar if u in listnode],axis=1),[listnode.index(u) for u in hiddvar if u in listnode],axis=0)\n",
        " \n",
        "    # Delete non observable component\n",
        "    for u in hiddvar:\n",
        "         if u in listnode:\n",
        "            listnode.remove(u)\n",
        "    Gsub = nx.relabel_nodes(nx.from_numpy_matrix(adj_mtx2), {i:j for i, j in enumerate(listnode)})\n",
        "    \n",
        "    # I have to guarantee the topological order inside each C-component and \n",
        "    # among different C-component\n",
        "    \n",
        "    #Get C-component\n",
        "    C_comp = [np.array(list(component_i)) for component_i in nx.connected_components(Gsub)]\n",
        "    \n",
        "    #Order the c-component\n",
        "    C_topo = []\n",
        "    posmax = []\n",
        "    for ii in range(0,len(C_comp)):\n",
        "        posmax = [orderlist.tolist().index(Sub_C_i) for Sub_C_i in C_comp[ii]]\n",
        "        C_topo.append(max(posmax))\n",
        "    \n",
        "    C_topo_sort = C_topo.copy()\n",
        "    C_topo_sort.sort()\n",
        "    \n",
        "    newind = [C_topo_sort.index(ii) for ii in C_topo]\n",
        "    C_comp = [C_comp[idx] for idx in newind]\n",
        "    #Order element inside C-component\n",
        "    C_comp = [toposort(C_comp[ii], orderlist) for ii in range(0,len(C_comp))]\n",
        "    return C_comp\n",
        "\n",
        "def probabilitystr():\n",
        "      return{\"01_var\" : np.array([]), \"02_cond\" : np.array([]), \"03_sumvar\" : np.array([]),\"05_fraction\" : False, \"07_sum\" : False, \"08_prod\" : False, \"09_branch\" : [] ,\"04_numerator\" : [], \"06_denominator\" : []}\n",
        "\n",
        "def an(x,G,orderlist):\n",
        "\n",
        "    anc = []\n",
        "    for xx in x.tolist():\n",
        "        anc.extend(list(nx.ancestors(G,xx)))\n",
        "        anc.extend(list(x))\n",
        "    anc = np.array([ordi for ordi in orderlist if ordi in anc])\n",
        "    return np.array(anc)\n",
        "\n",
        "\n",
        "def ID(y, x, P, G, orderlist, hiddvar,**kwargs):\n",
        "\n",
        "    P_out  = P.copy()\n",
        "    if (\"debug\" in kwargs):\n",
        "        debug  = kwargs[\"debug\"]\n",
        "    else :\n",
        "        debug = False\n",
        "    if (\"debug\" in kwargs):\n",
        "        debug2  = kwargs[\"debug2\"]\n",
        "    else :\n",
        "        debug2 = False\n",
        "\n",
        "    # Line 0 (preliminary): i) Get observed (current) graph\n",
        "    #                       ii) Get the list of (currente) node\n",
        "    if len(ispresent(hiddvar,G)) > 0:\n",
        "        G_obs = removeedge(ispresent(hiddvar,G),G,ispresent(hiddvar,G))\n",
        "        v     = np.array(list(nx.nodes(G_obs)))\n",
        "        v     = np.array([ordi for ordi in orderlist if ordi in v])\n",
        "    else:\n",
        "        G_obs = G.copy()\n",
        "        v     = np.array(list(nx.nodes(G_obs)))\n",
        "        v     = np.array([ordi for ordi in orderlist if ordi in v])\n",
        "        \n",
        "    if debug2:\n",
        "        print(\" \")\n",
        "        print(\"--------- Graph information ------------------\")\n",
        "        print(\"> The edge of the current graph are:\")\n",
        "        print(\">>>\" + str(list(nx.edges(G))) + \";\")\n",
        "        print(\" \")\n",
        "        print(\"> The edge of the current observed graph are:\")\n",
        "        print(\">>>\" + str(list(nx.edges(G_obs))) + \";\")\n",
        "        print(\"---------------------------------------------\")\n",
        "        print(\" \")\n",
        "    \n",
        "    #Line 1 of the algorihtm\n",
        "    if len(x)==0:\n",
        "        if debug:\n",
        "            print(\">> Enter in Line 1 of the algorithm\")\n",
        "        #Correct cases in Line 6 and 7\n",
        "        if P_out[\"08_prod\"]==True or P_out[\"05_fraction\"]==True:\n",
        "            P_out[\"03_sumvar\"] = toposort(np.union1d(np.setdiff1d(v,y),P_out[\"03_sumvar\"]), orderlist)  #list((set(v)-set(y)) | set(P_out[\"03_sumvar\"]))\n",
        "            P_out[\"07_sum\"]    = True\n",
        "        else:\n",
        "            if len(v)==1: #==0 Means just P(var)\n",
        "                 P_out[\"07_sum\"]    = True\n",
        "                 P_out[\"03_sumvar\"] = toposort(np.union1d(np.setdiff1d(v,y),P_out[\"03_sumvar\"]), orderlist)   #list((set(v)-set(y)) | set(P_out[\"03_sumvar\"]))\n",
        "                 \n",
        "            P_out[\"01_var\"]    = y\n",
        "\n",
        "        return P_out\n",
        "\n",
        "\n",
        "    anc = an(y,G_obs,orderlist)\n",
        "    if len(np.setdiff1d(v,anc)) != 0: #len(set(v) - set(anc)) != 0:\n",
        "        if debug:\n",
        "            print(\">> Enter in Line 2 of the algorithm\")\n",
        "        #Remove non ancestors(Y) node form the graph\n",
        "        #G_an    = G.copy()\n",
        "        edgedel = np.array([u for u in list(nx.nodes(G_obs)) if u not in anc])\n",
        "        G_an = removeedge(edgedel,G,ispresent(hiddvar,G))\n",
        "        \n",
        "        #G_an.remove_nodes_from(edgedel)\n",
        "        \n",
        "        if P_out[\"08_prod\"]==True or P_out[\"05_fraction\"]==True:\n",
        "            P_out[\"03_sumvar\"] = toposort(np.union1d(np.setdiff1d(v,anc),P_out[\"03_sumvar\"]), orderlist)  #list((set(v)-set(anc)) | set(P_out[\"03_sumvar\"]))\n",
        "            P_out[\"07_sum\"]    = True\n",
        "        else:\n",
        "            P_out[\"01_var\"]    = anc    \n",
        "        P_iter = ID(y,np.intersect1d(x,anc),P_out,G_an,orderlist,hiddvar,**kwargs)\n",
        "        return P_iter\n",
        "    G_bar_x  = incomingedgedel(x,G_obs)\n",
        "    an_bar_x = an(y,G_bar_x,orderlist)\n",
        "    W        = np.setdiff1d(np.setdiff1d(v,x),an_bar_x) # list((set(v)-set(x)) - set(an_bar_x))\n",
        "    if len(W) > 0:\n",
        "        if debug:\n",
        "            print(\">> Enter in Line 3 of the algorithm\")\n",
        "        #x      = x + W\n",
        "        P_iter = ID(y,np.union1d(x,W),P_out,G,orderlist,hiddvar,**kwargs)\n",
        "        return P_iter\n",
        "    G_less_x   = removeedge(x,G,hiddvar)\n",
        "    C_G_less_x = C_components(G_less_x,hiddvar,orderlist)\n",
        "\n",
        "    #-> Line 4\n",
        "    if len(C_G_less_x) > 1:\n",
        "        if debug:\n",
        "            print(\">> Enter in Line 4 of the algorithm\")\n",
        "            print(\">>> Condition ' if C(G\\X)={S1, ..., Sk} is valid';\")\n",
        "\n",
        "        P_iter = []\n",
        "        for p_iter in range(0,len(C_G_less_x)):\n",
        "            P_iter_i = ID(C_G_less_x[p_iter],np.setdiff1d(v,C_G_less_x[p_iter]),P_out,G,orderlist,hiddvar,**kwargs)\n",
        "            P_iter.append(P_iter_i)  \n",
        "                        \n",
        "        if len(np.setdiff1d(v,np.union1d(y,x))) != 0:\n",
        "            P_out[\"07_sum\"]    = True\n",
        "        P_out[\"08_prod\"]   = True\n",
        "                \n",
        "        P_out[\"03_sumvar\"] = toposort(np.setdiff1d(v,np.union1d(y,x)), orderlist)\n",
        "        P_out[\"09_branch\"] = P_iter\n",
        "        \n",
        "        return P_out\n",
        "\n",
        "    if len(C_G_less_x) == 1:\n",
        "        if debug:\n",
        "            print(\">> Enter in Line 4 of the algorithm\")\n",
        "            print(\">>> Condition ' if C(G\\X)={S} is valid';\")\n",
        "            \n",
        "        C_G = C_components(G,hiddvar,orderlist)\n",
        "        # Line 5 \n",
        "        if len(C_G)==1: \n",
        "            if len(np.setdiff1d(v,C_G))==0:  \n",
        "                raise ValueError('Graph form an Hedge: causal effect is not identifiable')\n",
        "\n",
        "        if any([C_G_less_x[0].tolist() == C_Gi.tolist() for C_Gi in C_G]):\n",
        "            if debug:\n",
        "                print(\">> Enter in Line 6 of the algorithm\")\n",
        "            P_prod_b = []\n",
        "            \n",
        "            for ii in range(0,len(C_G_less_x[0])):\n",
        "                if P_out[\"08_prod\"]==True:\n",
        "                    #Compute condition set \n",
        "                    cond_list = condset(C_G_less_x[0][ii],v)\n",
        "                    P_prod_i  = complex_prob(P_out,C_G_less_x[0][ii],cond_list,v,orderlist)                 \n",
        "                else:\n",
        "                    P_prod_i         = P_out.copy()\n",
        "                    P_prod_i[\"01_var\"]  = np.array(C_G_less_x[0][ii])\n",
        "                    cond_list        = condset(C_G_less_x[0][ii],v)\n",
        "                    P_prod_i[\"02_cond\"] = cond_list\n",
        "                P_prod_b.append(P_prod_i)                    \n",
        "                    \n",
        "            #If the C-component have more than one variable\n",
        "            if len(C_G_less_x[0]) > 1:\n",
        "                P_6 = probabilitystr()\n",
        "                P_6[\"03_sumvar\"]  = toposort(np.setdiff1d(v,np.union1d(y,x)), orderlist)\n",
        "                P_6[\"09_branch\"]  = P_prod_b.copy()   \n",
        "                P_6[\"08_prod\"]    = True\n",
        "            else:\n",
        "                P_6 = P_prod_b[0].copy()\n",
        "                if P_6[\"08_prod\"] or P_6[\"05_fraction\"]:\n",
        "                    P_6[\"03_sumvar\"]  = np.union1d(P_6[\"03_sumvar\"],np.setdiff1d(C_G_less_x[0],y))\n",
        "                else:\n",
        "                    P_6[\"01_var\"] =  np.setdiff1d(P_6[\"01_var\"], np.union1d(P_6[\"03_sumvar\"],np.setdiff1d(C_G_less_x[0],y)))\n",
        "            \n",
        "            return P_6        \n",
        "                    \n",
        "        #-> Line 7\n",
        "        #For testing if C_G_less_x is a subset of C_Gi I use the difference \n",
        "        #function with array: if length is 0 the is a subset\n",
        "        if any([len(np.setdiff1d(C_G_less_x,C_Gi))==0 for C_Gi in C_G]): \n",
        "            if debug:\n",
        "                print(\">> Enter in Line 7 of the algorithm\")\n",
        "\n",
        "            S            = [C_Gi for C_Gi in C_G if len(np.setdiff1d(C_G_less_x,C_Gi))==0][0]\n",
        "            P_           = probabilitystr()\n",
        "            P_prod_b     = []\n",
        "\n",
        "            for iii in range(0,len(S)):\n",
        "                \n",
        "                if P_out[\"08_prod\"]==True:\n",
        "                    cond_list = condset(S[iii],v)\n",
        "                    P_prod_i  = complex_prob(P_out,S[iii],cond_list,v,orderlist)\n",
        "                    \n",
        "                else:\n",
        "                    P_prod_i         = P_out.copy()\n",
        "                    P_prod_i[\"01_var\"]  = np.array(S[iii])\n",
        "                    cond_list        = condset(S[iii],v)\n",
        "                    P_prod_i[\"02_cond\"] = cond_list\n",
        "                    \n",
        "                P_prod_b.append(P_prod_i) \n",
        "                     \n",
        "            if len(S) > 1:\n",
        "                P_[\"08_prod\"]   = True                        \n",
        "                P_[\"09_branch\"]  = P_prod_b.copy()\n",
        "            else:\n",
        "                P_ = P_prod_b[0].copy()\n",
        "                \n",
        "            edgedel = [u for u in list(nx.nodes(G_obs)) if u not in S]\n",
        "            G_s     = removeedge(edgedel,G,hiddvar)\n",
        "            P_iter  = ID(y,np.intersect1d(x,S),P_,G_s,orderlist,hiddvar,**kwargs)\n",
        "    \n",
        "            return P_iter  \n",
        "        \n",
        "def complex_prob(P, var, condset, v, orderlist):\n",
        "    \n",
        "    if len(condset) == 0:\n",
        "        P_complex = P.copy()    \n",
        "        if P_complex[\"07_sum\"]==True:\n",
        "            P_complex[\"03_sumvar\"] = toposort(np.union1d(P_complex[\"03_sumvar\"],np.setdiff1d(v,np.union1d(var,condset))), orderlist)\n",
        "        else:\n",
        "            P_complex[\"07_sum\"] = True\n",
        "            P_complex[\"03_sumvar\"] = toposort(np.setdiff1d(v,np.union1d(var,condset)), orderlist)\n",
        "            \n",
        "   #If we have the condition set express the conditional probability in terms of \n",
        "   # conjunction probability require a denominator\n",
        "   \n",
        "    else:\n",
        "        P_complex = probabilitystr()    \n",
        "        P_complex[\"05_fraction\"]    = True\n",
        "        P_complex[\"04_numerator\"]   = P.copy()\n",
        "        P_complex[\"06_denominator\"] = P.copy()\n",
        "        if P[\"07_sum\"]==True:\n",
        "            P_complex[\"04_numerator\"][\"03_sumvar\"]   = toposort(np.union1d(P[\"03_sumvar\"],np.setdiff1d(v,np.union1d(var,condset))), orderlist)\n",
        "            P_complex[\"06_denominator\"][\"03_sumvar\"] = toposort(np.union1d(P[\"03_sumvar\"],np.setdiff1d(v,condset)), orderlist)\n",
        "        else:\n",
        "            P_complex[\"04_numerator\"][\"03_sumvar\"]   = toposort(np.setdiff1d(v,np.union1d(var,condset)), orderlist)\n",
        "            P_complex[\"06_denominator\"][\"03_sumvar\"] = toposort(np.setdiff1d(v,condset), orderlist)\n",
        "           \n",
        "       \n",
        "    return P_complex\n",
        "      \n",
        "    \n",
        "def probabilitystr():\n",
        "      return{  \"01_var\" : np.array([])\n",
        "             , \"02_cond\" : np.array([])\n",
        "             , \"03_sumvar\" : np.array([])\n",
        "             , \"04_numerator\" : []\n",
        "             , \"05_fraction\" : False\n",
        "             , \"06_denominator\" : []\n",
        "             , \"07_sum\" : False\n",
        "             , \"08_prod\" : False\n",
        "             , \"09_branch\" : [] }\n",
        "    \n",
        "    \n",
        "    \n",
        "def IDC(y, x, z, P, G, orderlist, hiddvar,**kwargs):\n",
        "    \n",
        "    debug  = kwargs[\"debug\"]\n",
        "    \n",
        "    #Transform array y, z and z in array\n",
        "    y_set = set(list(y))\n",
        "    \n",
        "    if all(x) != \"\":\n",
        "        x_set = set(list(x))\n",
        "    z_set = set(list(z))\n",
        "    #Get graph without incoming edge on x and outgoing edge to z\n",
        "    G_under_x_bar_z = incomingedgedel(x,outcomeedgedel(z,G))\n",
        "    for zi in z:\n",
        "        c_set = x_set.copy()\n",
        "        c_set | (z_set - {zi})\n",
        "        if nx.d_separated(G_under_x_bar_z,y_set,{zi},c_set):\n",
        "            if debug:\n",
        "                print(\"> Iterate over IDC: Z variable is \" +str(zi))\n",
        "            P_numerator = IDC(y, np.union1d(x,np.array([zi])),np.setdiff1d(z,np.array([zi])), P, G, orderlist, hiddvar,**kwargs)\n",
        "            return P_numerator      \n",
        "    if debug:\n",
        "        print(\"> Iterate over ID\")\n",
        "    P_numerator = ID(np.union1d(y,z), x, P, G, orderlist, hiddvar,**kwargs)\n",
        "\n",
        "    return P_numerator        \n",
        "\n",
        "        \n",
        "def ce_identification(y,x,z,P,G,orderlist,hiddvar,**kwargs):\n",
        " \n",
        "   debug  = kwargs[\"debug\"]\n",
        "\n",
        "   if z.tolist():\n",
        "       \n",
        "       if debug:\n",
        "           print(\"@@ Condition set z is not empty: IDC version is used\")\n",
        "           \n",
        "       P_numerator = IDC(y, x, z, P, G, orderlist, hiddvar,**kwargs)\n",
        "       \n",
        "       P_prime = probabilitystr()\n",
        "       \n",
        "       #Construct the denominator\n",
        "       P_prime_denominator = P_numerator.copy()\n",
        "       \n",
        "       if P_prime_denominator[\"07_sum\"]:\n",
        "           P_prime_denominator[\"03_sumvar\"] = np.union1d(y,P_prime_denominator[\"03_sumvar\"])\n",
        "       else:\n",
        "           P_prime_denominator[\"07_sum\"]    = True\n",
        "           P_prime_denominator[\"03_sumvar\"] = y\n",
        "           \n",
        "       #Add the denominator to whole distribution\n",
        "       P_prime[\"05_fraction\"]    = True\n",
        "       P_prime[\"06_denominator\"] = P_prime_denominator\n",
        "       P_prime[\"04_numerator\"]   = P_numerator\n",
        "\n",
        "   else:\n",
        "       if debug:\n",
        "           print(\"@@Condition set z is empty: ID version is used\")\n",
        "           \n",
        "       P_prime = ID(y,x,P,G,orderlist,hiddvar,**kwargs)\n",
        "       \n",
        "       \n",
        "   return P_prime\n",
        "      "
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# function\n",
        "\n",
        "# INPUT: Markovian graph (including observed(blue)/unobserved(grey)), Targer, Intervention (node names), Context variables (as a set)\n",
        "# OUTPUT: True/False (for the identifiability), [] a set of all possible adjustment sets\n",
        "\n",
        "def IDC_with_markov(G, T, I, C):\n",
        "\n",
        "  grey_nodes = [n for n,v in G.nodes(data=True) if v['color'] == 'grey']\n",
        "\n",
        "  # TODO: encapsulate all these pre-condtions in a auxiliary function\n",
        "  if T in grey_nodes:\n",
        "    print('The target variable (T) cannot be unobserved.')\n",
        "    return False\n",
        "  if I in grey_nodes:\n",
        "    print('The intervention variable (I) cannot be unobserved.')\n",
        "    return False\n",
        "  for c in C:\n",
        "    if c in grey_nodes:\n",
        "      print('The contextual variables (C) cannot be unobserved.')\n",
        "      return False\n",
        "\n",
        "  final_graph, hidden_confounders = markov_to_semimarkov(G)\n",
        "\n",
        "  DG_1 = nx.DiGraph()\n",
        "  DG_1.add_edges_from(list(final_graph.edges()))\n",
        "    \n",
        "  hiddvar          = np.array(hidden_confounders)                # declare an array containign the hidden variables\n",
        "  Gobs1            = removeedge(hiddvar,DG_1,hiddvar)            # gets the observed graph (so without hidden variable)\n",
        "  order            = np.array(list(nx.topological_sort(Gobs1)))  # retreive the topological ordering of the graph\n",
        "  kwargs           = dict()                                      # specify others input\n",
        "  kwargs[\"debug\"]  = False                                        # diplay the enter interrogated at each step\n",
        "  kwargs[\"debug2\"] = False                                       # if true, edges of graphs at each steps are diplayed\n",
        "  P                = probabilitystr()                            # initalize probability structure\n",
        "\n",
        "  try:\n",
        "    Pi1 = ce_identification(np.array([T]),np.array([I]),np.array(C), P,DG_1,order,hiddvar,**kwargs)\n",
        "  except:\n",
        "    print('Based on the current setting the query is not identifiable.')\n",
        "\n",
        "    power_set = powerset(grey_nodes)\n",
        "    adj_dict = {}\n",
        "    adj_sets = []\n",
        "\n",
        "    for grey_set in power_set:\n",
        "      H = G.copy()\n",
        "\n",
        "      grey_list = list(grey_set)\n",
        "      for grey_node in grey_list:\n",
        "        nx.set_node_attributes(H, {grey_node: \"blue\"}, name=\"color\")\n",
        "\n",
        "      final_graph, hidden_confounders = markov_to_semimarkov(H)\n",
        "\n",
        "      DG_1 = nx.DiGraph()\n",
        "      DG_1.add_edges_from(list(final_graph.edges()))\n",
        "      hiddvar          = np.array(hidden_confounders)                       \n",
        "      Gobs1            = removeedge(hiddvar,DG_1,hiddvar)            \n",
        "      order            = np.array(list(nx.topological_sort(Gobs1)))                                       \n",
        "      P                = probabilitystr()                            \n",
        "\n",
        "      try:\n",
        "        Pi1 = ce_identification(np.array([T]),np.array([I]),np.array(C), P,DG_1,order,hiddvar,**kwargs)\n",
        "      except:\n",
        "        print(\"The query is still not identifiable with \" + str(grey_list) + \" being observed.\")\n",
        "      else:\n",
        "        print(\"The query becomes identifiable if \" + str(grey_list) + \" are observed.\")\n",
        "        adj_dict[str(grey_list)] =  Pi1\n",
        "        adj_sets.append(grey_list)\n",
        "\n",
        "    return False, remove_redundant_sets(adj_sets)       # for probability expressions we can use the dict\n",
        " \n",
        "  else:\n",
        "    print(\"Based on the current setting the query is already identifiable.\")\n",
        "    return True, Pi1"
      ],
      "metadata": {
        "id": "05OpjO09IUKx"
      },
      "execution_count": 14,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Example\n",
        "G = nx.DiGraph()\n",
        "G.add_nodes_from([\n",
        "    ('V5', {\"color\": \"grey\"}),\n",
        "    ('V4', {\"color\": \"blue\"}),\n",
        "    ('V3', {\"color\": \"blue\"}),\n",
        "    ('V2', {\"color\": \"blue\"}),\n",
        "    ('V1', {\"color\": \"grey\"}),\n",
        "])\n",
        "G.add_edge('V5', 'V2')\n",
        "G.add_edge('V5', 'V3')\n",
        "G.add_edge('V4', 'V3')\n",
        "G.add_edge('V3', 'V2')\n",
        "G.add_edge('V1', 'V5')\n",
        "G.add_edge('V1', 'V2')\n",
        "G.add_edge('V4', 'V5')\n",
        "\n",
        "\n",
        "nx.draw(G, with_labels=True, font_weight='bold', node_size=1000, node_color = [x[1] for x in G.nodes('color')])\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 319
        },
        "id": "w2OChWAdKvPk",
        "outputId": "1fd6a91b-397c-4158-9b2e-b107af79ec2b"
      },
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAABBwklEQVR4nO3deViU9f7/8efAAIoL7qamuOGGiRsgoJK54L4vKZqWlmnZBn0759Q51TllHTXz16JWauVyzHJX3DXFBUFEFBRQXHFX3FgFZub3x8DkuLHNzD3MvB/XxYUOw32/uUpfft73Z1HpdDodQgghhJ1wULoAIYQQwpIk+IQQQtgVCT4hhBB2RYJPCCGEXZHgE0IIYVck+IQQQtgVCT4hhBB2RYJPCCGEXZHgE0IIYVck+IQQQtgVCT4hhBB2RYJPCCGEXZHgE0IIYVck+IQQQtgVtdIFCFFW5ebmcu3aNa5cuUJmZiYajQa1Wk3lypWpW7cuNWrUwMFB/m0phLWR4BOiGHJzczl+/DgHDhwgNTUVJycnNBoNeXl5hvc4OTmhUqnQaDTUr1+fgIAAmjRpgkqlUrByIUQBlRxEK0ThNBoN4eHhREREoFKpyMnJKfL3Ojs74+zsTFBQEJ6enhKAQihMgk+IQly9epU//viDtLQ0cnNzS3wdJycnGjRowODBg6lYsaIJKxRCFIcEnxBPcfToUTZu3GjUyiwNBwcHnJ2deemll6hTp45JrimEKB4JPiGeICYmhs2bN5ss9B7k7OzM+PHjqVu3rsmvLYR4OplyJsRjnDp1ymyhB5CTk8PixYu5c+eOWa4vhHgyCT4hHpKVlcXq1avNFnoFcnJyWLVqFdJ0EcKyJPiEeMiGDRtKNYmlqHQ6HdeuXePw4cNmv5cQ4i8SfEI84PLlyyQnJ6PRaCxyv9zcXLZv326RoBVC6EnwCfGAiIgIs7c4H+f48eMWv6cQ9kqCT4h8WVlZJCYmWvyZW05ODvv27bPoPYWwZ7JlmRD5Tp8+/dS9NRctWsSFCxcIDAykW7duhtcvXbrETz/9hIODA0OHDiUyMpIrV66Ql5eHu7s7L7/8cqH3vnv3Lvfu3aNy5com+VmEEE8mIz4h8l28ePGpW5G1adMGgPj4eKPX4+LiAPDw8CAzM5Pc3Fxq1apVrHs7Ojpy+fLlYlYshCgJCT4h8p0/f/6pX/f09MTR0ZHU1FRDSGm1WkMQtmnTBh8fH15//XWee+65Yt07JyeHS5culaxwIUSxSPAJka+wxeTly5fHw8MD+GvUd+7cOdLT03FxcaF58+YlvnfB0gYhhPlJ8AmRT6vVFvqeB9udOp3O0Ob09PRErS7dI3NZ0iCEZUjwCZGvKMcFNWvWjHLlynHv3j3OnDlDQkIC8FcgloYcWiuEZcifNCHyOTk5FfoetVpNq1atAAgLCyM7Oxs3Nzfc3d1Lff/y5cuX+hpCiMJJ8AmR75lnninS+wpGd7du3QLgueeeM4wWz58/z5o1awwt0Js3b7JmzRr27t371Guq1Wrq169f0tKFEMUgwSdEPnd3dxwdHYv0Pjc3N8PvH2xz3rp1i6NHjxpmfWZkZHD06FGSk5Ofes2cnBzef/99pkyZwvLly7l48WIJfwohRGHkPD4h8p0/f57//e9/T13LZy6Ojo706tWL/fv3Ex4ezt69e6lcuTJdu3Y1fDRp0qRIzyGFEE8nwSdEPp1Ox+zZs0lPT7fofVUqFa1atWL48OGG17RaLYmJiYSHhxs+tFqtIQS7dOmCp6enTIgRogQk+IR4wIEDB9i9e7dFlxY4OTkxYcKEp57GrtPpOHv2LHv37jUE4a1bt+jcubMhDNu1a1fqJRVC2AMJPiEekJWVxZw5cyzW7lSpVNSuXZvJkycX+3svX75sFITnz5+nU6dOhiD08fGhXLlyZqhaiLJNgk+Ihxw9epSwsDCLjPrUajWTJ0+mRo0apb5Wamqq4RlheHg4J06coH379obWqL+/P5UqVTJB1UKUbRJ8QjxEp9OxZMkSzp8/X6TdXErKycmJrl270rlzZ7NcPy0tjYiICMOo8PDhw7Rs2dIwIuzcuTPVq1c3y72FsGYSfEI8RkZGBnPnziUzM9Ms11er1TRs2JDRo0dbbIJKdnY2hw4dMowIIyIiaNCggdGEmXr16lmkFiGUJMEnxGOcOHGCQYMGMWHCBLRarUlHfk5OTtSrV4/g4GBFJ6Pk5eURGxtrCMK9e/dStWpVunTpYgjDxo0byxIKYXMk+IR4yIkTJ+jRowczZsxg8ODBLF++nBs3bpjkmZ9araZ169b079+/SIvlLUmr1XLixAlDa3TPnj2oVCqjEWGrVq1kCYUo8yT4hHjA8ePH6dmzJzNmzGDs2LGA/pnfoUOH2LFjBxqNpkSjP7VajbOzM8OGDaNx48amLtssdDodZ86cMRoR3r5922hE2LZtW1lCIcocCT4h8hWE3syZMwkODn7k63fv3uXAgQPExsYCFGnJg7OzM87Ozvj5+dGhQwdcXFxMXbZFXbp0yWgJRUpKCn5+foYw9Pb2liUUwupJ8AnBX6E3a9YsxowZ89T35ubmkpCQwOnTp7l48SJ37txBpVKhUqko+ONUo0YN3N3dadasmU0/J0tNTWXfvn2GIExISKBDhw6GEaGfnx8VK1ZUukwhjEjwCbsXHx9Pr169ihR6j6PRaMjJyUGj0eDo6Ei5cuVsNugKk5aWxoEDBwyjwpiYGFq1amW0hKJatWpKlynsnASfsGvx8fH07NmT2bNnM3r0aKXLsTnZ2dlERUUZRoQHDx6kYcOGhtZoly5dnrpVmxDmIMEn7FZcXBy9evWS0LOg3NzcR5ZQVK9e3WjmaKNGjex2xCwsQ4JP2KWC0Pv666958cUXlS7Hbmm1Wo4fP240YcbBwcHoOKaWLVtKEAqTkuATdufYsWMEBQUxZ84cRo0apXQ54gE6nY7Tp08bjQjv3r1rtITCy8tLllCIUpHgE3ZFQq/suXjxotGI8OLFi/j7+xtao97e3mV+mYiwLAk+YTeOHj1KUFAQ33zzDSNHjlS6HFFCN2/eNFpCkZiYSMeOHY2WUFSoUEHpMoUVk+ATdkFCz3bdu3fPaAnFkSNH8PT0NFpCUbVqVaXLFFZEgk/YvILQ+/bbbxkxYoTS5Qgzy8rKemQJRePGjQ2t0S5dulCnTh2lyxQKkuATNi02NpbevXtL6Nmx3Nxcjhw5YgjCffv2UaNGDaOZo+7u7jJz1I5I8AmbVRB63333HcOHD1e6HGEltFot8fHxhlmj4eHhqNVqoyBs0aKFBKENk+ATNunIkSP06dOH77//nmHDhildjrBiOp2O5ORkw4gwPDyc9PT0R5ZQWNsxUqLkJPiEzZHQE6WVkpJitITi8uXLRksoOnbsKEsoyjAJPmFTYmJi6NOnD/PmzWPo0KFKlyNsxI0bN4yWUCQlJeHt7W0YEXbq1MlullDk5sKJE3D4MJw5A+np4OwMNWtCu3bQvj1Y+z7kEnzCZkjoCUspOJux4DlhbGwsrVu3NgRhQECATS2h0Gph506YMQP27IFy5fSvZWT89R4nJ3B1hawsqFUL3n4bXnnFOkNQgk/YhILQmz9/PkOGDFG6HGFnsrKyiIyMNIwIIyMjadKkidEpFM8884zSZZbIunUwdSrcu6cf3RWVq6s+HMeNg6++gkqVzFdjcUnwiTLv8OHD9O3bV0JPWI3c3FxiYmKMllDUqlXLEIJlYQlFaiq8+ips3QqZmSW/Trly+tBbvhy6dzddfaUhwSfKtOjoaPr168cPP/zA4MGDlS5HiMfSaDTEx8cbTZhxdnY2WkLRvHlzqwnChATo2hXS0uD+fdNc09UV/vY3+OgjUPrHlOATZVZB6P34448MGjRI6XKEKDKdTsepU6eM1hJmZGQYLaFo06aNIksoTpwAf399a9PU6eDqqn/2N326aa9bXBJ8okyS0BO25sKFC0YjwitXrhAQEGAIw44dO+Ls7GzWGm7cgJYt4dYt04deAVdXmD0bJk82z/WLQoJPlDmHDh2if//+/PTTTwwcOFDpcoQwi+vXrxstoTh58iQ+Pj6G54TmWEIxYABs2wY5OSa97CNcXSE+Hho1Mu99nkSCT5Qphw4dol+/fixYsEBCT9iVB5dQhIeHExsbS5s2bYyWUFSpUqXE11+5EiZMMF6iYC6OjtChA0REgIOD+e/3MAk+UWZERUXRv39/Fi5cyIABA5QuRwhFZWZmGpZQ7N2717CEoiAIu3TpQu3atR/5vjVr1lC3bl18fX0Nr2k0ULcuXL9uuforVoRly0CJf79K8IkyQUJPiKfLyckxWkKxf/9+ateu/cgSCnd3d65evcrChQsZO3YsABs3wpgx+lmclhQQAPv2WfaeIMEnyoDIyEgGDBjAokWL6N+/v9LlCFEmFCyheHDzbScnJ65cuYJWq6VcuXK8/fbbfPHFFwQGqti71/I1lisHcXHQtKll7yvBJ6yahJ4QpqHT6Zg/fz7vvfce2dnZhtf9/XsQFbWdvLwnfWdXYC/wMfDJA68fAnwAJ+BDIAxIBnKAFsC/gKf3MV1c9Esb3nuvJD9RySnwWFGIojl48CADBgzg559/ltATopRUKhXnz58nJycHV1dXOnfuTEhICCNGfIar69O+c2z+5+UPvf6//M99gD+BO+iD7jngMDAMOPLUmu7f1+/9aWky4hNW6eDBgwwcOJBffvmFvn37Kl2OEDbh8uXL3L59m5YtW+KQP53yu+/g//5Pv7n0490G6gD30Qdae0ALPAtcAf4AGue/DqABmgFngK+Apw/n6taFS5dK8UOVgIz4hNWJiIiQ0BPCDOrWrYunp6ch9EC/PdmTQw+gKlDw57Bg1Pcn+tBzAwbwV+gVKFgI+GyhNV29WuhbTE6CT1iViIgIBg0aJKEnhIUUbd1eQbtzBaDjrwAcATx8IO97wEXAHyj8eDCdjqc8XzQPtWVvJ8STHThwgMGDB/Prr7/Sp08fpcsRwi6oi5QC/YAqQAqwE1iV//rYB96jASYDC4GOwEaKEjE6nX5BuyXJiE9YhQMHDjBo0CAWL14soSeEmeh0Oh6e1lGrVlFOS3BBP7oDmIJ+Ios7+hmfANnoJ7MsBHqib4UW7SDe8uUtf1qDBJ9QXMFIb8mSJfTu3VvpcoSwWSEhIVSvXp1Bgwbxww8/kJiYSLt2uiIeElswukvO/zwGKEisicA6oBzgAXwEvMNfMz+frEWLotdvKtLqFIrav38/Q4YMYcmSJQQFBSldjhA2rXHjxqSnp7N+/XrCwsLQaDR4eQ0iN3dtEb67C/pR3vn83z/Y5iyYlpkNzH3g9fHoA/LxVCro0qXI5ZuMBJ9QzL59+xg6dKiEnhBmcv/+fWJiYjhw4AD79+9nz5495OXPJHFwcKBFixZs3jyfVq0Km9kJ+tHduSd8bXeJ6qtYUX/graVJq1MoYt++fQwZMoSlS5dK6AlhIjdu3GDdunV88MEHdO7cmerVq/PGG29w7tw5RowYwaFDh3B2dsbV1ZWXXnqJI0eOUKfOM7z2mn4XFUtTqUCJvSlkAbuwuILQW7ZsGb169VK6HCHKJK1WS0JCAgcOHDCM6K5fv06nTp0ICAjA398fX19fKlasaPR9EydOpEOHDkydOtXwWkoKeHjod1KxFGdn/WnsM2ZY7p4FJPiERe3du5ehQ4fyv//9j549eypdjhBlRkZGBlFRUYaQO3jwINWqVTOEXEBAAK1atTJanF4cgwfD5s3mP4S2gKsrJCZC/fqWud+DJPiExYSHhzNs2DAJPSGK4OLFi+zfv98QdAkJCXh5eRmCzt/f/7Hn7ZXUtWv6UZ8ljiaqUAH+8x94913z3+txJPiERRSE3vLly+nRo4fS5QhhVfLy8jh27JhR0GVnZxtGcv7+/nTo0IFy5cqZtY7ffoOJEyEz03z3cHSEtm0hKkqZ09dBgk9YgISeEMbu3LnDwYMHDUF36NAhGjRoYBR0TZs2RWXhld06Hbz+Oixdap7wU6mgRg04fFiZFmcBWc4gzGrPnj0MHz6c3377je7duytdjhAWp9PpOH36tNFo7vz583h7exMQEEBoaCidOnWiatWi7XRiTrdupVKx4pd06jSNgwcbmDT8HB2hWjXYv1/Z0AMZ8Qkz2r17NyNGjJDQE3YlOzubmJgYQ9AdOHAAZ2dnAgICDKM5Ly8v1EXbJNPsMjMz2bBhA/Pnz2f//v3k5uby1VezuXfvXWbMgOxs/UiwNCpU0J+yvn49NGhgmrpLQ4JPmEVB6K1YsYIXXnhB6XKEMJvr168bRnIHDhzg6NGjtGjRwqhtWV/pIc5TvPLKKyxevBiNRgOAk5MTt27domLFisTGwogRcOVKUU9xMObkpN8E+9NP9aesW3oz6ieR4BMmJ6EnbJVWq+XEiRNGQXfz5k38/PwMQeft7f3I2jlrdvPmTdq3b09KSgoAXbp0ITw83PD1nBxYtky/3i4lRT8CzM/IJ6pYUT9KfPllfeA1amTOn6D4JPiESf3555+MHDmS33//nW7duildjhClkp6eTlRUlCHkDh48SI0aNYzWzj14mnlZFBMTQ1BQEBUqVODixYvMmzePV1999bHvPXRI364MD4ejRyE9XT+K02r1E1caNwZ/f+jeHYYM0a/Vs0YSfMJkJPREWXfhwgWj0VxiYiJt27Y1WjtXq1Ytpcs0mfj4eHr06MG8efPo0aMHoaGhTJ8+nerVqxfp+7Oz9Xt8OjnpjxeyllZmYST4hEns2rWLUaNG8ccff/D8888rXY4QhcrNzeXo0aNGQZeTk2MUch06dMBFiU0sLSApKYkXXniBWbNmMXr0aKXLsSgJPlFqO3fu5MUXX2TlypUEBgYqXY4Qj3X79m0iIiIMQRcdHU3Dhg2N2paNGze2+No5JZw9e5bAwEA+/fRTXn75ZaXLsTgJPlEqEnrCGul0OpKTk43WzqWkpBjWzvn7+9OpUyeqVKmidKkWl5KSQmBgIKGhoUYbVdsTCT5RYhJ6wlpkZ2cTHR1t1LZ0dXU1WlLQpk0bq1k7p5SrV6/StWtXJk+eTEhIiNLlKEaCT5TIjh07GD16NKtWraKrEidJCrt27do1o9HcsWPHaNWqlVHQPfvss0qXaVVu3rzJ888/z4svvshHH32kdDmKkuATxbZ9+3bGjBkjoScsQqPRcOLECaOgu337Nn5+foaQ8/b2pkKFCkqXarVu375N9+7d6d27N59//rldPMd8Ggk+USwFobd69Wq6dOmidDnCBqWlpT2ydq527dpGo7kWLVqU6bVzlpSWlkbPnj3p1KkTX3/9td2HHkjwiWLYtm0bwcHBEnrCZHQ6HRcuXDDa1/LkyZO0a9fOEHR+fn7UrFlT6VLLpMzMTPr06UPLli2ZN2+ehF4+CT5RJNu2bWPs2LGsXr2azp07K12OKKNyc3OJjY01altqNBqjJQXt2rWz2bVzlpSdnc3AgQOpU6cOP//8s4yQHyDBJwpVMNJbs2aNhJ4ollu3bhEREWEIusOHD9O4cWOjtmWjRo1kJGJiOTk5DBs2DFdXV5YtW2b3s1kfJsEnnmrr1q2MGzeONWvWEBAQoHQ5worpdDpOnjxptKTg4sWL+Pr6GoLO19cXNzc3pUu1aXl5eYwePZqcnBxWrlyJk5OT0iVZHQk+8URbtmzhpZdektATj5WVlfXI2rmKFSsatS1bt24tow0L0mq1jB8/nuvXr7Nu3TrKlSundElWSYJPPFZB6K1duxZ/f3+lyxFW4MqVK0YhFxcXh6enp9HelvXq1VO6TLul0+mYPHkyJ0+eZNOmTbha69EIVkCCTzxCQk9oNBri4+ONgu7u3buGgCtYOyd/uVoHnU7HO++8Q1RUFNu2baNSpUpKl2TVJPiEkc2bNzN+/HjWrVuHn5+f0uUIC0lLS+PgwYOGoIuMjKROnTpGbctmzZrJzEArpNPp+Pvf/8727dvZuXOnXe4/WlwSfMJg06ZNTJgwQULPxul0Os6fP2+0pCA5OZn27dsbgs7Pz48aNWooXaoogv/85z/8/vvv7N69u8jn6Nk7CT4B/BV669evp1OnTkqXI0woJyfnkbVzOp2OgIAAQ9C1a9cOZ2dnpUsVxTRr1ix++ukn9uzZwzPPPKN0OWWGBJ8gLCyMl19+WULPRqSmphp2QSlYO9e0aVOjtXMNGzaUtXNl3Pfff89XX31FeHi4bMhdTBJ8dm7jxo288sorbNiwAV9fX6XLEcWk0+lISkoyGs1duXIFX19fw3ZfnTp1onLlykqXKkxo0aJFfPLJJ+zZs4dGjRopXU6ZI8FnxyT0yp7MzEyio6ON9rZ0c3MzGs21bt0aR0dHpUsVZvK///2P999/nz///JNmzZopXU6ZJMFnpyT0yobLly8bLSmIj4/nueeeM9rAuW7dukqXKSxk9erVTJ06lR07dtC6dWulyymzJPjs0IYNG5g4cSIbN27Ex8dH6XJEPo1GQ1xcnFHQpaWlGdbNBQQE0LFjR8qXL690qUIBmzZt4uWXX2bLli20a9dO6XLKNAk+O7NhwwYmTZrExo0b8fb2Vrocu3b37l0iIyMNIRcVFUW9evWM2pbNmjWTSSiCnTt3Mnr0aJmAZiISfHZk/fr1TJo0ibCwMAk9C9PpdJw9e9ZoNHf69Gk6duxoCLpOnTrJOizxiH379jFkyBBWrVpF165dlS7HJkjw2Yl169bx2muvyUjPQnJycoiJiTEKOgcHB6OdULy8vGTtnHiqqKgo+vfvz7Jly+jZs6fS5dgMCT47UBB6YWFhdOzYUelybNLNmzcNsyz379/PkSNH8PDwMAq6Bg0aSNtSFFlsbCxBQUEsWLCAAQMGKF2OTZHgs3ESeqan1WpJTEw0Gs1dvXqVTp06GYLO19dXNgoWJXbixAm6d+/Ot99+y/Dhw5Uux+ZI8NmwtWvXMnnyZDZt2kSHDh2ULqfMyszMJCoqyhB0ERERVK1a1WgSiqenp6ydEyaRnJzM888/z5dffsnYsWOVLscmSfDZqDVr1vD6669L6JXApUuXjHZCOXHiBG3atDE6d072RRTmcP78ebp27cpHH33Eq6++qnQ5NkuCzwYVhN7mzZtp37690uVYtby8POLi4oyCLjMz02g017FjRznJWpjdpUuXCAwM5K233uKtt95SuhybJsFnY1avXs2UKVMk9J7g7t27HDx40GjtXP369Y2CzsPDQyahCIu6fv06gYGBTJgwgQ8++EDpcmyeBJ8NWbVqFVOnTpXQy6fT6Thz5ozRJJSzZ88+snauWrVqSpcq7NitW7d4/vnnGTJkCJ9++qnS5dgFCT4bsWrVKt544w02b95st9sZ3b9/n5iYGKMNnNVqtdG5c15eXjg5OSldqhCAvgPRo0cPunXrxn//+1/pNFiIBJ8NsNfQu379OhEREYagi42NpXnz5kZtywYNGihdphCPlZ6eTlBQEO3bt+ebb76R0LMgCb4ybuXKlbz55pts2bKFtm3bKl2O2Wi1WhISEozaltevX8fPz88QdD4+PlSsWFHpUoUoVFZWFv369aNx48b8+OOPODg4KF2SXZHgK8NsOfQyMjIeWTtXvXp1o51QWrVqJX9hiDLn/v37DB48mGrVqrF48WJZ/6kACT4LysrK4sqVK1y/fp2cnBx0Oh3Ozs7UrFmTOnXqUKFChad+/+3btxkwYAA///wzsbGxTJs2ja1bt+Ll5WWhn8B8UlJSjEZziYmJeHl5GdbN+fv7U7t2baXLFKJUcnNzGTlyJA4ODqxYsQK1Wq10SXZJgs/M0tLSOHToEEeOHCErKwu1Wk1eXh4ajQYAR0dHw2suLi54eXnh4+NDlSpVHrnWokWLeP3113FxccHFxYWdO3eWydDLy8vj6NGjRkGXnZ1t9GyuQ4cOsnZO2BSNRsPYsWO5d+8ea9askQ3KFSTBZyb37t0jLCyM06dPAxiCrjCOjo6oVCrq169P//79jabad+7cmf379wPg5uZGdHQ0TZs2NX3xJnbnzh0iIiIMQXfo0CHc3d2N2pZNmjSRh/vCZmm1WiZOnEhKSgobN26Uf9QpTILPxHQ6HbGxsWzZsoW8vDy0Wm2JrqNSqXB0dKR79+74+vpy+/ZtateuTV5eHg4ODjg6OjJy5EiWLl1q4p+gdHQ6HadPnzbaCeX8+fN4e3sbgq5Tp05UrVpV6VKFsAidTscbb7xBXFwcW7ZsKfSRhjA/aTCbkEajYeXKlZw+fZrc3NxSXUun05GXl8euXbtISkri5MmT5OXlUblyZcaMGUNwcDD+/v4mqrzksrOzOXz4sGHd3IEDB3BxcTGM5F577TW8vLzkWYawSzqdjtDQUKKjo9mxY4eEnpWQEZ+JaDQali9fzoULF0odeg9Tq9VUqlSJevXqMWTIEEVnMl67ds3o3LmjR4/SsmVLo+dz9evXV6w+IazJP//5TzZs2MCuXbtkhyArIv8MN5H169ebJfRAPxkkLS2NrKwsiz4H02q1nDhxwqhtmZqaalg79/nnn+Pt7S1r54R4jOnTp7Nq1Sr27NkjoWdlZMRnAklJSaxatcosofcgJycn+vTpY7bdWdLT04mKijIE3cGDB6lZs6bRaK5ly5aydk6IQsyZM4fvv/+e8PBw6tSpo3Q54iESfKWUlZXFN998Q3Z2tkXu5+TkxJtvvknlypVLfa0LFy4YLSlISkqibdu2hqDz8/OjVq1aJqhaCPvxww8/8MUXXxAeHi5b5lkpaXWWUnh4uNlHeg/Ky8tj+/btDBs2rFjfl5ub+8jauZycHMNIbsyYMbRv3x4XFxczVS6E7fv111/57LPP2L17t4SeFZPgK4W8vDxiYmKKvEbPFHQ6HYmJiWRlZVG+fPknvu/27dtGa+eio6Np2LAhAQEB9OvXj+nTp9O4cWNZOyeEifz+++/8/e9/Z9euXTRp0kTpcsRTSPCVwvHjxxW7d2xsLH5+foA+DE+dOmU0mktJScHHxwd/f3/+9re/4evr+9jdYIQQpbdu3Treeusttm3bRosWLZQuRxRCnvGVwuLFizl79uxjv7Zo0SIuXLhAYGAg3bp1M7x+6dIlfvrpJxwcHOjYsSNRUVGPfO/7779f6HqfSpUqkZeXZwg6V1dXo51QnnvuOVk7J4QFbN26lXHjxrFp0yY6duyodDmiCORvxlK4evXqE7/Wpk0bLly4QHx8vFHwxcXFAeDh4WFoVbZs2dJoskpRDkq9e/cud+/eJTg4mO+//55nn322pD+GEKKEdu/ezbhx41i7dq2EXhkiwVdCGRkZ5OTkPPHrnp6ebN68mdTUVC5fvkzdunXRarXEx8cD+mC8fv06AD4+PjRq1KhY93dxceHDDz+kZs2aJf8hhBAlFhERwciRI1mxYoVV7KIkik6Cr4Ru3bqFWq1+4sSW8uXL4+HhQWJiIvHx8dStW5dz586Rnp6Oi4sLzZs3NwTfb7/9hkajMZw316ZNm0Lv7+DgQGpqqgSfEAo4fPgwgwYNYvHixUYdHVE2yErkEirKEoaCAIuPj0en0xnanJ6enqjVahwcHHB3d6d169a4u7tz7do1Vq9eTXJycqHXLtjLUwhhWXFxcfTr148ff/yR3r17K12OKAEZ8ZVQUXYvadasGeXKlePevXucOXOGhIQE4K9A7Nq1K4GBgYb3r1y5kvj4eBISEop03JDsoCKEZSUlJREUFMScOXMYPHiw0uWIEpK/OUvIxcWFwibEqtVqWrVqBUBYWBjZ2dm4ubnh7u4O6Nulj1OUtXUqlUoWmwthQWfOnKFHjx5Mnz6dF198UelyRCnIiK+EatasWaRWY5s2bYiJiTGE3HPPPWcItiVLllCxYkVq1arF3bt3OX36NCqVitatWxd63by8PJ555pnS/RBCiCJJSUmhe/fu/OMf/2DChAlKlyNKSUZ8JaRWq4u0X6a7uztubm6G3z84caV9+/bk5uYSHx/PxYsXqV+/PqNHj6Zhw4aFXtfJyUnO9hLCAq5cuUL37t2ZNm0aU6ZMUbocYQKygL0U1q1bx9GjRwtteZpD06ZNCQ4Otvh9hbAnN27c4Pnnn2fMmDF8+OGHSpcjTERGfKXQsWNHRXZHcXJywtfX1+L3FcKe3L59m169ejF48GAJPRsjwVcK9erVM8nxQMV17949pk6dyrZt2xQZbQph6+7du0fv3r15/vnn+eyzz5QuR5iYBF8pBQYGFmmLMVNxcnJiwIABBAcHExISgpeXF4sXL37qLjJCiKLLyMigf//+tG/fntmzZ8sJJjZIgq+UWrduTd26dS2ypk6lUlG1alX8/PwYP348x44dY+bMmSxZsoTGjRszY8YM7ty5Y/Y6hLBV2dnZDB48mMaNG/P9999L6NkoCb5SUqlUDB06FEdHR7Pfy9HRkREjRhhCVqVSERQUxPbt29m4cSNxcXE0adKEkJAQLly4YPZ6hLAlOTk5DB8+nOrVq7Nw4ULZIMKGyX9ZE6hcuTLDhg0z60QXtVpNv379qFGjxmO/3rZtW5YsWUJsbCwODg60a9eO4OBgYmJizFaTELYiLy+PMWPGoFarWbJkiUX+ISuUI8FnIs2bN2fgwIFmCT+1Wk3Pnj1p27Ztoe+tX78+M2fO5MyZM7Rv355BgwbRvXt3Nm/eLBNhhHgMjUbDhAkTSE9PZ8WKFRZ9Zi+UIev4TOzUqVOsWrWK3NxctFptqa7l4OCAo6MjAwcOLNJuLo+Tm5vL77//zsyZM8nLyyMkJIQxY8bIdmdCoN/s/bXXXiM5OZmwsDBcXV2VLklYgASfGWRkZLBu3TrOnTtXpFMcHsfJyYk6deowbNgwkyyZ0Ol07Ny5k1mzZnHs2DGmTZvG66+/TtWqVUt9bSHKIp1Ox9tvv010dDRbt26lUqVKSpckLESCz4ySk5PZu3cvly9fRqfTPfHsvgIODg44ODhQu3ZtAgICaNGihVlmlR07dozZs2ezfv16xo0bxzvvvFPsg3CFKMt0Oh1/+9vf2LFjBzt37qRKlSpKlyQsSILPjFavXs2IESM4deoUKSkpnDt3jqtXr5KTk2N4eK7ValGr1dSuXRt3d3fatGnzxAkspnbp0iW+/fZbFixYQPfu3QkNDcXb29si9xZCSf/+97/5448/2L17N9WrV1e6HGFhEnxmMm/ePN5++23y8vI4fvw4LVu2NHwtOzvbsODcycmJ8uXLK1UmAGlpaSxcuJCvv/6ahg0b8v7779O3b1+Zzi1s0syZM1m4cCF79uyhdu3aSpcjFCDBZ2JarZb333+f+fPnk5mZSYUKFdi6dSsBAQFKl1aovLw8Vq5cycyZM8nMzCQkJISxY8dSrlw5pUsTwiS+++47vv76a/bs2cOzzz6rdDlCIfJPehNLTk7mm2++MUxqcXR0fOKBs9ZGrVbz4osvEh0dzdy5c1mzZg0NGzbks88+IzU1VenyhCiVhQsXMnPmTHbu3CmhZ+ck+EysWbNmXLp0CU9PTypUqEB6enqZCb4CKpWKbt26ERYWxs6dOzl79iweHh68+eabnD59WunyhCi2ZcuW8fHHH7Njx44inXcpbJsEnxk4Ojpy7tw5zpw5w+bNmwkKClK6pBLz9PRk4cKFHD9+HDc3Nzp16sTw4cM5ePCg0qUJUSSrVq0iNDSUbdu24eHhoXQ5wgrIMz4zmDVrFvHx8fzyyy9Kl2Jy6enp/Pzzz8yePZt69eoRGhrKwIEDZSKMsEphYWG88sorbNmyhXbt2ildjrASEnwmptVqadasGcuWLbPpw2Lz8vJYs2YNM2fO5M6dO4SEhPDSSy8pPkNViAI7duxgzJgxbNiwwab/LIrik3+mm9j27dtxc3PDx8dH6VLMSq1WM2LECCIjI1m4cCFhYWE0bNiQTz/9lBs3bihdnrBze/fuZfTo0axatUpCTzxCgs/E5s6dy9SpU+3mHC+VSkWXLl1Yv349e/bs4dKlSzRr1owpU6Zw8uRJpcsTdigqKophw4axfPlyunTponQ5wgpJ8JnQhQsX2LdvHy+++KLSpSiiRYsW/PjjjyQmJlKzZk06d+7MkCFD2L9/v5wMISwiNjaWAQMGsGjRInr06KF0OcJKyTM+E/roo49IT09nzpw5SpdiFTIzM/nll1+YPXs2NWvWJDQ0lMGDB8tZZ8IsTpw4Qffu3fnuu+8YNmyY0uUIKybBZyI5OTk0aNCAPXv20Lx5c6XLsSoajYZ169Yxc+ZMrl+/znvvvceECROoUKGC0qUJG3Hq1Cm6devGl19+ydixY5UuR1g5aXWayOrVq/H09JTQewxHR0eGDh1KREQES5YsYefOnTRq1Ih//vOfXLt2TenyRBl37tw5evTowSeffCKhJ4pEgs9E5s2bx9SpU5Uuw+r5+/uzevVq9u/fT2pqKi1btuS1114jMTFR6dJEGXTp0iXDySKTJk1SuhxRRkjwmUB8fDzJyckMHDhQ6VLKDA8PD+bOnUtSUhLPPvssgYGBDBw4kPDwcJkII4rk2rVrdO/encmTJzNt2jSlyxFliDzjM4E33niDmjVr8sknnyhdSpmVlZXF4sWL+eqrr6hSpQqhoaEMHToUtVqtdGnCCqWmptKtWzeGDh0qf+5EsUnwlVJaWhru7u7ExcVRr149pcsp87RaLRs2bGDWrFlcvHiRd999l1deeYWKFSsqXZqwEnfu3KFHjx50796dL7/80m7WzArTkVZnKS1btoxu3bpJ6JmIg4MDgwYNYu/evSxfvpy9e/fSsGFD/vGPf3DlyhWlyxMKS09Pp1+/fvj7+0voiRKT4CsFnU5n2KlFmF6nTp34448/iIyMJC0tDU9PT1555RWOHz+udGlCAVlZWQwYMIBWrVoxZ84cCT1RYhJ8pXDgwAFycnJ44YUXlC7FpjVp0oRvv/2WU6dO0aRJE3r06EHfvn3ZtWuXTISxE/fv32fIkCHUq1eP+fPny2kgolTkGV8pBAcH4+3tzTvvvKN0KXYlOzubZcuWMWvWLFxdXQkNDWX48OE4OTkpXZowg9zcXEaMGIFarea3336TCU+i1CT4Suj69es0b96cM2fOULVqVaXLsUtarZZNmzYxa9Yszp49yzvvvMOkSZOoVKmS0qUJE9FoNAQHB5Oens7q1atxdnZWuiRhA6RfUEKLFi1i2LBhEnoKcnBwoH///uzevZuVK1cSGRlJo0aN+OCDD7h06ZLS5YlS0mq1TJo0idTUVFauXCmhJ0xGgq8ENBoN8+fPZ8qUKUqXIvJ5e3vz22+/ER0dzf3793nuuecYP348x44dU7o0UQI6nY4333yT06dPs3btWsqVK6d0ScKGSPCVwJYtW6hduzYdOnRQuhTxkIYNGzJnzhxOnz5Ny5Yt6d27N0FBQWzfvl0mwpQROp2OkJAQDh8+zMaNG2Uzc2Fy8oyvBPr168fIkSMZP3680qWIQty/f5/ly5cza9Ys1Go1oaGhjBo1SibCWLF//vOfbNy4kV27dsmjBGEWEnzFdPbsWXx8fLhw4QLly5dXuhxRRDqdjq1btzJr1iySkpJ4++23efXVV3Fzc1O6NPGA6dOns2zZMnbv3k3NmjWVLkfYKGl1FtMPP/zA+PHjJfTKGJVKRe/evdmxYwfr168nNjaWxo0bExoaSkpKitLlCeDrr7/m559/ZseOHRJ6wqwk+IohOzubRYsW8frrrytdiiiFdu3asXTpUo4cOYJOp8PLy4uxY8dy5MgRpUuzW/Pnz+ebb75h586d1KlTR+lyhI2T4CuGlStX0q5dO5o2bap0KcIEGjRowFdffcWZM2fw8vJiwIAB9OjRgy1btshEGAv69ddf+fzzz9mxYwcNGjRQuhxhB+QZXzEEBATwf//3fwwaNEjpUoQZ5OTksGLFCmbNmoVWqyUkJITRo0fj4uKidGk2a8WKFbz77rvs2rWLFi1aKF2OsBMSfEUUGxvLwIEDOXPmjGyZZON0Oh07duxg1qxZxMfHM23aNCZPniwzDE1s3bp1TJ48me3bt/Pcc88pXY6wI9LqLKJ58+bx2muvSejZAZVKRc+ePdm6dSubNm0iISGBJk2a8M4773Du3Dmly7MJW7du5bXXXiMsLExCT1icBF8R3L17lz/++INJkyYpXYqwMC8vL3799VeOHTuGs7MzHTp0YPTo0URHRytdWpm1e/duxo0bx5o1a2QTCKEICb4iWLJkCb169eKZZ55RuhShkGeffZYZM2Zw9uxZvL29GTp0KN26dSMsLAytVqt0eWVGREQEI0eO5Pfff8ff31/pcoSdkmd8hdDpdHh6ejJv3jwCAwOVLkdYidzcXP744w9mzpzJ/fv3CQkJITg4WPaUfIrDhw/Tt29fFi9eTFBQkNLlCDsmI75ChIeHo1Kp6Nq1q9KlCCvi5OTEmDFjiImJ4bvvvmPVqlU0atSIzz//nNTUVKXLszpxcXH069ePH3/8UUJPKE6CrxBz585l6tSpqFQqpUsRVkilUvHCCy+wadMmtm/fzunTp/Hw8GDatGmcOXNG6fKsQlJSEkFBQfy///f/ZCmQsAoSfE9x5coVtm3bxrhx45QuRZQBrVu3ZtGiRcTHx1OpUiV8fHwYOXIkkZGRSpemmDNnztCjRw+++OILRo0apXQ5QgASfE+1cOFCRo0aReXKlZUuRZQhdevWZfr06Zw9e5aAgABGjRpF165dWb9+vV1NhElJSaF79+58+OGHcpKJsCoyueUJ8vLyaNSoERs3bsTLy0vpckQZlpeXx6pVq5g5cyZpaWmEhIQwbtw4m97o/MqVKwQGBjJlyhTeffddpcsRwoiM+J5g48aNNGjQQEJPlJparWbUqFEcOnSIH3/8kQ0bNtCoUSP+/e9/c/PmTaXLM7kbN27Qo0cPJkyYIKEnrJIE3xPMmzePqVOnKl2GsCEqlYrAwEA2bNjAn3/+SUpKCh4eHkydOpVTp04pXZ5J3L59m169ejFkyBD+8Y9/KF2OEI8lrc7HOHXqFJ07d+bChQuyQbEwq6tXr/L9998zf/58unTpQmhoqEUXdp86BVu3Qng4REbCjRuQmwuOjuDmBu3aQWCg/sPPD542ufnevXv07NmTzp07M2vWLJkJLayWBN9jhISE4OTkxJdffql0KcJOZGRk8MsvvzB79mxq165NaGgogwYNwtHR0eT30mhg/Xr473/h2DH9a1lZT36/kxO4uED16vD++/DSS1Cp0qP19+nTB09PT+bOnSuhJ6yaBN9DsrKyaNCgAVFRUTRq1EjpcoSd0Wg0rF27lpkzZ5Kamsp7773H+PHjcXV1Ncn1ExJg5Eg4dw7S04v//RUq6ENw6VLo00f/WnZ2NgMGDODZZ59l4cKFODjIExRh3eT/0IesWLECX19fCT2hCEdHR4YNG0ZERAQ///wzW7dupWHDhnz88cdcv369xNfV6fQjvA4d4PjxkoUeQEYG3LoFw4fD6NFw+3YOw4cPp3r16ixYsEBCT5QJ8n/pQ+bOncuUKVOULkPYOZVKRefOnVm7di179+7l2rVrNG/enMmTJ5OUlFSsa2m1MGkS/Oc/+pamKXo8mZmwdi00aXIFna4yS5YsMUtbVghzkOB7QHR0NDdu3KB3795KlyKEQfPmzZk/fz5JSUnUqVOHLl26MGjQIPbu3UvBk4p79+4xffp0NBqN0ffqdDB1Kvz2m360ZkrZ2ZCe/iyXLy8lN9fJtBcXwozkGd8DJk6cSLNmzfjggw+ULkWIJ8rMzGTx4sV89dVXVKtWjdDQUJKTk/nwww955513mD17tuG9v/4Kb7xh+tB7ULlyMGIELF5svnsIYUoSfPlu375N48aNOXnyJDVr1lS6HCEKpdFo2LBhAzNmzCAyMhKtVkv58uX56aefCA4O5vJlaN685M/zisPVFdasgV69zH8vIUpLrXQB1uLXX3+lX79+EnqizHB0dGTw4MGkp6cTGxtLVlYWWVlZjB8/ntq1a/Pf//YgO9sytWRmQnAwnD+vD0EhrJk84wO0Wq3h+CEhypq1a9eSnZ2Ns7Mzrq6uqNVq/vWv5ezfD3l5lqsjK0v/LFEIayfBB+zatYvy5cvj5+endClCFNsff/yBRqPh/v37ZGRkkJ2dTZs2C8nNtWwdGRnw5ZemmTUqhDlJ8PHXvpyy24Qoi1QqldH/uzk5sGTJ00Z7XQEV8MlDrx/Kf90Z+BRoCpQHKgMdgd8LreXyZThypHj1C2Fpdh98Fy9eZPfu3QQHBytdihAmERcH6qc+vR+b/3n5Q6//L/9zH+AO0Bp4GWgHHAZGA8lPvbdGA/v2Fa9eISzN7oPvp59+YvTo0VSsWFHpUoQwicOHC3u2NwJwAU4CMfmvaYEV+b8eB3wNrAXmArsBt/z3XHjqvbOz9RteC2HN7HpWZ25uLgsWLGDbtm1KlyKEyezfr59l+WRVgb7AGvSjvvbAn8AV9AE3IP99W4Aw4BhwF+ic//F0hw6VsHAhLMSuR3zr1q3Dw8MDT09PpUsRwmSuXi3KuwranSsAHX+1PQtGgwAHge+AcKAc+rAsfIeWu3eLXqsQSrDr4JN9OYUtyskpyrv6AVWAFGAnsCr/9bEPvOcTIBc4CtQC/sFfzwGfzJJLKIQoCbsNvoSEBBISEhgyZIjSpQhhUuXLF+VdLuhHdwBT0E9mcUc/4xMgLf+zGmgDtMj//clCr+zsXLQ6hVCK3Qbf/PnzmThxIs7yp1SUYZ999hnTpk0jLCyMe/fuAfptyoq2MqdgdFcwU3MM+uUMAHWAgehDsQewDf1fFz0LvWrdukWtXghl2OXkloyMDJYuXcoRWXAkyrijR4+ycuVKfv31V7KysqhSpQq9ey+iYsUBpKUV9t1d0I/yzuf//sE2Zw/06/q2ol/HFwiEUpTJLZ0Lf4sQirLL4Fu+fDldunShQYMGSpciRLFlZ2dz5MgRIiMjOXv2LCqVirT8lEtPT6dr1wqsXl2UK6mAc0/42toS1VaxIvj7l+hbhbAYuzudQafT0b59e7788kuCgoKULkeIp9JqtZw6dYqoqCgiIyOJjIzkxIkTtGjRAh8fH2rXrs2MGTMAaNmyJZs2baJWrdrUrw+XLlm+XhcXOHNG2p3Cutld8B08eJCxY8dy8uRJHBzs9hGnsFI3btwwCrlDhw5RuXJlfH19DR/t27enfP4MltzcXNzc3AgODmbu3Lk4OemXG3z9NXz0UWHr+UxLpYI+fSAszHL3FKIk7C74xo8fT5s2bQgJCVG6FGHnHmxZFnzcunULb29vQ8gVjOqe5tatW1SrVs3otTt3oE4dLHYsEUCFCrBhA3TrZrl7ClESdhV8N2/exMPDg+TkZKpXr650OcKOFLQsCwIuKirK0LJ8MOSaN29usk7Ev/4FX31lmVGfWg2dOum3K5O93oW1s6vgmzlzJsePH+eXX35RuhRh427cuGEIuIKWpZubGz4+Po9tWZpDbi60bg2nTpn/qKAKFSAhAerXN+99hDAFuwk+rVaLh4cHy5cvx8fHR+lyhA0xVcvSHI4dAz8/8476XF3h22/hlVfMdw8hTMlugm/Lli18+OGHREdHy7l7osQebllGRkaSkJBg1LL09fWlWbNmVjN5KiwMRozQn5BuahUqwFtvwfTppr+2EOZiN8E3cOBABg0axMSJE5UuRZQhBS3LB2dZVqlSxWgkZ+6WpSls2vRX+JnqT7yrK4SGwiefyHM9UbbYRfCdP3+eDh06cOHCBVxdXZUuR1iprKwsQ8uy4Nnc7du38fb2NjybU6plaQrHj+vD7/z50rU+XVz0obd4MfTvb7r6hLAUuwi+Dz/8kIyMDObMmaN0KcJKlMWWpSnk5cEXX+g/HBwgI6Po3+vioh/ZDR4Mc+dC1apmK1MIs7L54MvJyaFBgwbs2bOH5s2bK12OUEhhLUtfX1/atWtn9S1LU0lLg6VLYcYM/fl9zs76ENRo/nqPSqXfgkyrBScneOMNmDIF6tVTrm4hTMHmg++3335jwYIF7NixQ+lShIU82LIs+Lhz584jsyxr1aqldKlW4fJlOHwYoqLg5En9c0AXF33A+fpChw7QtKl+hCiELbD54OvatStvv/02w4YNU7oUYQZarZaTJ08aLQxPSEigZcuWhoCzxZalEKLkbDr44uLi6N27N+fOnTPsYSjKtuvXrz+yMNyeW5ZCiOKz6eB74403qFWrFh9//LHSpYgSkJalEMIcbDb40tLScHd3Jy4ujnryNN7qPdyyjIyMJDEx0dCyLAg5aVkKIUrLZg+iXbp0KS+88IKEnpUqaFkWtC0PHTpE1apVDSE3duxYaVkKIczCJkd8Op0OLy8v5syZwwsvvKB0OXYvKyuLmJgYo2dzd+7cwcfHx2hhuLQshRCWYJPBt2/fPiZNmkRCQoLsy2lhWq2WpKQko8NUH25Z+vr64uHhIS1LIYQibLLVOXfuXKZMmSKhZwEPtiwLZllWq1ZNWpZCCKtlcyO+69ev07x5c86ePUuVKlWULsemPNiyLPi4e/euUbtSWpZCCGtnc8H3xRdfcPr0aRYsWKB0KWVaQcvywedyiYmJtGrVyugwVWlZCiHKGpsKPo1GQ+PGjVmzZg3t27dXupwy5dq1a0bP5R5uWRYsDC9XrpzSpQohRKnY1DO+zZs3U6dOHQm9QhTWsnznnXfw8fGhZs2aSpcqhBAmZ1Mjvr59+zJq1CjGjx+vdClW48GWZcFHUlISrVq1MtrLUlqWQgh7YfXBp9PB6dP63ePPnYPsbFCroXp1aNcOvLygXDk4c+YMvr6+XLhwwa5nEF67ds0o5KKjo6levbpRyEnLUghhz6yy1anTwYEDMHMmbNumPw7FwUF/XEpenv7X5crpAzAzEzw8wN09huDgV+0q9DIzM4mJiTF6Nvdgy/Ldd9+VlqUQQjzE6kZ8W7boD7y8dk0fakWtztk5B0dHNdOmOfDpp/pgtCWFtSwfnGUp6xeFEOLJrCb47t7Vn+68bp0+8EqqfHmoVQt+/x18fExXn6U9rWVZ8NG2bVtpWQohRDFZRfCdOgVdusCdO3D/vmmuWb68vlX6xhumuZ45FbQsH9y0+d69e4/sZSktSyGEKD3Fgy8pCfz89KFn6kpcXeHf/4aQENNetzS0Wi2JiYlGz+WkZSmEEJajaPDdvAmtWuk/m6sKV1dYtAhGjTLP9QsjLUshhLAuigbfoEH6ySw5Oea9T8WKcPIk1Klj3vs83LKMjIwkLS3NaIsvb29vaVkKIYSCFAu+1ath3LjSTWQpKicnCAzUL40wVfewoGX5YMidPHkST0/PRxaGS8tSCCGshyLBp9HoR183bljunhUq6GeMdu9esu+/evWq0eST6OhoatSoYTSak5alEEJYP0WCb8MGCA6GtDTL3jcoSN9aLUxRWpY+Pj7UqFHD/EULIYQwKUWCr3Nn2L/f0nfVL2pPToZ69f56rbCWZcFH06ZNpWUphBA2wOLBl54O1apBbu6T3tEV2At8DHzywOuHAB/ACbgCVAeSgbZABuAFxD713q6u8K9/3aVFi91Gsyxr1qxpNJKTlqUQQtguiwffvn3Qrx/cu/ekd/wITAaaAUkPvP4uMAcYCKwD8oDOwOH8XxcefAAuLit5/vkF0rIUQgg7ZfFNqqOjC9udZQTwFnASiAHaA1pgRf7Xx+V//hw4BoQA/y3y/evVG86WLcOLV7QQQgibYfED2OLjCwu+qkDf/F8vz//8J/r2phswAH3b8zNgFtCiWPe/eLFYbxdCCGFjLB58RZvJOTb/8wpAx18BOALQ5H+9JzC12Pd/8rNFIYQQ9sDirU51ke7YD6gCpAA7gVX5r49FP9o7CVQD+gMFQ7iz+b/f+NQryyHjQghh3yweA0XbrcsF/egOYApwB3BHP+OzYC7OQSAMOJr/+3v5v386OzqnVgghxGNYPPg6dtTvnVm4gnZncv7nMYAKeB59+BV8/Jz/dS/+CsUna9my6LUKIYSwPRYPvg4divrOLuhHeQXGPumNRaZSQdeupb6MEEKIMszi6/g0GnBzg4wMS95Vr1IlWLwYBg+2/L2FEEJYB4uP+Bwd4aWX9CcmWJpKBX37Fv4+IYQQtkuROY5vv60PQEtydobXX9d/FkIIYb8UCb7mzcHPr6hLG0xDrYY337Tc/YQQQlgnxVa1/fILuLhY5l4VKsC//w3161vmfkIIIayXYsHXoAHMmqUPJXNydIRmzeDdd817HyGEEGWDIufxFdDpYORI2LQJMjNNf32VSn8E0uHD4O5e+PuFEELYPkU38FKpYPly6N5df1aeKTk46ENv3z4JPSGEEH9RfOdKtRrWrIGJE00Xfq6u+vbmoUPQoniHNwghhLBxigcf6J/DffMNbNkCdeqUPAAdHfV7cX7wAcTFQaNGpq1TCCFE2afoM77HycyEH36Ar76Cu3f1O7wUVmGFCqDVwosv6kOveXPL1CqEEKLssbrgK6DTwa5dsHYt7N0LiYn61x0d9V/LzdWPDn18ICgIRo2CypUVLVkIIUQZYLXB9zCtFm7dgqws/XZnbm5yxJAQQojiKzPBJ4QQQpiCVUxuEUIIISxFgk8IIYRdkeATQghhVyT4hBBC2BUJPiGEEHZFgk8IIYRdkeATQghhVyT4hBBC2BUJPiGEEHZFgk8IIYRdkeATQghhVyT4hBBC2BUJPiGEEHZFgk8IIYRd+f+/xFA2SWjxPgAAAABJRU5ErkJggg==\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "G_folan, l = markov_to_semimarkov(G)\n",
        "nx.draw(G_folan, with_labels=True, font_weight='bold', node_size=1000, node_color = 'red')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 319
        },
        "id": "QA_T-9j1PnK8",
        "outputId": "722c3073-21e8-43f1-a902-8d108958ad2f"
      },
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAABJcElEQVR4nO3deVxU1RvH8Q+IQqipaZp7miSGASqIa5hLe6lZPzfcQiA1JcTdNjXNMjW11BSETC0t0xaXMk0tU1lEwAUJd1zIDUmRZeD+/rhgTqKyzMyd5Xm/XrzQmWHOYxpfnnPPOddOURQFIYQQwkbYa12AEEIIYUoSfEIIIWyKBJ8QQgibIsEnhBDCpkjwCSGEsCkSfEIIIWyKBJ8QQgibIsEnhBDCpkjwCSGEsCkSfEIIIWyKBJ8QQgibIsEnhBDCpkjwCSGEsCkSfEIIIWyKBJ8QQgib4qB1AUIIK5WeDrt3Q3Q07NgBZ89CTg5UqAB16oCvL3h7Q9u2ULWq1tUKG2InN6IVQhhUTAzMng3r16shl5kJOt3tr3NwAGdnNQx79IDQUPDyMnW1wgZJ8AkhDCM1Ffz81A4vKwvy84v/tfb24OSkdoArVkC9esarU9g8ucYnhCgbRYFly8DVFXbtUju8koQeqK/PzFS/3tVVfT/5mVwYiXR8QojSUxQICYGlS9XgMhRnZwgIgLlzwc7OcO8rBNLxCSFKqzD0wsIMG3qgvl9YGIwebdj3FQIJPiFEaUVEqJ3e9evGef/r12HJEnXaUwgDkqlOIUTJpaaq1+KMFXq3qlQJkpKgbl3jjyVsgnR8QoiS8/OD7GzTjJWVBf37m2YsYRMk+IQQJRMdrX4UtTfPGHQ6dbyYGNOMJ6yeBJ8QomRmz1a7MFPKyoI5c0w7prBaEnxCiOJLT4fvvy9yn94TgB3w3n8ejy54vAIwBWgC3AfcD3gBa4ozbn4+rFunji9EGUnwCSGKb/du9RiyIvgVfP7qP4+vKvj8LJAONAeGAC2AWKAvkFKcsR0dYc+ektUrRBEk+IQQxRcdfcc9e68CjkAysK/gsXxgdcGvBwBzgfXAQmA7UKXgNaeKM/b16+r4QpSRBJ8Qovh27LjjopZqwHMFvy7s+n4DzqEG3IsFj20GRgKdgKtAh4KPe9LpYPv2UhQthD4JPiFE8Z09e9enC6c7VwMK/wZgYTcIsAf4FNgJOKGGZfnijn/uXPFrFeIOJPiEEMWXk3PXp58HqgKnga3A2oLH/W55zXtALhAP1AQm8e91wHsy1d5BYdUk+IQQt/H39+fJJ59k0aJFnDp1yxW4OyxsKeSI2t0BDENdzNIQdcUnwD8Fnx0Ad8C14PfJxS3M0fHerxHiHiT4hBC3uXbtGtu3b2f06NE0adKE+++/n7lz56p3Tr+Hwu6ucKVmP9TtDAC1gZdQQ7Er8AvqN6FuxS2sdu3ivlKIO3LQugAhhPbOnz9PfHz8zY8dO3YAkJWVhb29PQ4ODnh5ecE//8DOnXc9taUjapd3suD3t05zdkXd1/cz6j4+X2AMxVzc4uAAnTqV8E8mxO0k+ISwIbm5uSQlJemFXHx8PDqdDg8PDzw8POjWrRudOnUiODgYOzs7nn/+eSIiIqhYsSJcu6beKy8j445j2AEn7vDc+rIUX7Gieod2IcpI7s4ghJW6dOnSbQF35MgRGjRocDPkCj/q1q2L3S03fE1PT6dx48bMmjULf39/bnlCnW409ZFlAE5O6qrOqlVNP7awKhJ8Qli4vLw8/vrrr9tC7p9//sHd3V0v4Jo3b46zs3Ox3ldRFL0wvKlPH/jmmyKPLTMae3vo3RtWFXv9pxB3JMEnhAW5evUqCQkJegF38OBBHnroodu6uIcffrjo4CqrmBjw9TX8XdfvxtlZvbbYqpXpxhRWS4JPCDOUn5/P8ePHb+vi/v77b5o3b64XcO7u7tx///2mLbBTJ9i1yzS3JnJwgPbt5dQWYTASfEJo7Pr16yQmJuoFXGJiIlWqVLmti2vSpAnlypXTumS5A7uwaBJ8QpiIoiikpqbqBdz+/ftJTU3F1dX1ti6uevXqWpd8d8uWwciRxp3ydHaGBQvgtdeMN4awORJ8QhhBVlYWhw4dum2qskKFCrd1cU2bNqV8+WKfVmk+FAVCQiAszDidX8WKEBAAc+ca/r2FTZPgE6KM/rv5Oz4+nqNHj9KkSZPbQq5WrVpal2tYioJu1CjyFi/G0ZDX+5ydITBQveu6MRboCJsmwSdEMRVn83fhx2OPPYajjZwrOW7sWOr9+isjU1Kwy8oq04KXvHLlyLWzw2HhQhwCAgxYpRD/kuAToghl2fxtS3bs2EHfvn2Jj4/nwexs8PNTbxablVWyfX729uQ6OLArJ4f+QF6tWixZsoQXX3zRZv/bCuOR4BM27U6bvzMyMorc/F2xYkWtSzYbV69excPDg4ULF/Lcc8/9+0RMjDpFuW6dejeF69eL7gIdHNTreNnZ0LMn0R070nr48JtPOzs74+HhwZo1a6hXr54J/kTCVkjwCZtxp83ftWrVKnLzt7293LzkbgYMGEDlypVZuHBh0S9IT4c9e9QOcPt29bix7Gw1DGvXVvcCentDmzZQtSpxcXG0b9+eGzdu3HwLZ2dnfv/9d1q2bGmKP5KwERJ8wuqY/eZvK7B69Wreeecd4uLiin0E2r2cOHECNzc3Mgu2R9jZ2bF582aeeuopg7y/EIUk+IRFs8jN3xYuNTWVli1bsmHDBrwNeLeE9PR0qlWrhpOTEwDDhw8nPj6eLVu2yHU+YVASfMIiWN3mbwuVn5/PU089RadOnXjrrbcM/v5bt26lY8eOLFiwgC1btnD58mX8/f0JCgoy+FjCdknwCbNjE5u/LdTcuXP55ptv2LlzJw4OxrudZ05ODu7u7gQHB/POO+8QExNDw4YNjTaesC0SfEJTNr3528IcOHCAJ598kj179vDII48YfbxNmzYRHBzMoEGD2L59O7/88otMeQqDkOATJiGbvy1bdnY2rVu3Jjg4mNdMeG7mCy+8QIcOHfjuu+8YOnQogYGBJhtbWC8JPmFwd9r8Xb9+fb2A8/T0tOnN35Zk7NixpKSk8N1335n07ys5OZl27dqxdu1aXnnlFZnyFAYhwSdKTTZ/24bt27fTr18/9XSWBx80+fhjx47l0qVLuLi4sG3bNpnyFGUmwSeKJSMj47bN3wcOHJDN31YuPT0dDw8PFi9ezLPPPqtJDRkZGbi6urJ27VpGjRpFQECATHmKMpHgE3ry8/M5ceLEbV1cWlqabP62QX5+flSpUoXPPvtM0zoiIiJYsmQJS5YsoXPnzjLlKcpEgs+GXb9+nQMHDugFXEJCgmz+FgB8/fXXvPfee+zbt89gp7OUVn5+Pm3atGHkyJGkpqby22+/8fPPP8uUpygV8w++9HTYvVs972/HDjh7FnJyoEIFqFMHfH3V8/7atoWqVbWu1iwVtfk7Pj6eU6dO3bb528PDQzZ/C06fPk2rVq3YuHEjXl5eWpcDwJ49e+jVqxcHDhzgqaeeIjAwkAC5dZEoBfMNvpgYmD0b1q9XQy4z884nvDs7q2HYoweEhoKZ/I+qhezsbA4ePHhbF+fg4HBbwLm6usrmb3Gb/Px8unXrRufOnZk8ebLW5egZOHAgdevWpX///jz55JPExsbSoEEDrcsSFsb8gi81tUz39MLJSe0AV6wAK7+VSVpamt7xXYWbvx955JHbQu6hhx7SulxhIebMmcPatWvZsWOHUU9nKY2zZ8/y+OOPExUVxerVq9m+fbtMeYoSM5/gUxSIiIBRo9Rbl5ThLs44OKi3Ppk/H4YMAQv/nyI3N5cjR47cNlWZk5NT5ObvwkN+hSipxMREOnfuzN69e2ncuLHW5RTpgw8+YO/evXz77be0adOGoKAgmfIUJWIewacoEBICS5eqU5qG4uwMAQEwd67FhN/ly5dvC7ikpKTbNn97eHhQr149+UlXGExWVhatW7cmJCSEIUOGaF3OHWVlZeHm5sbixYupXbu2THmKEtM++ApDLyxMvVOzoVWs+G/4mZG8vDxSUlJuC7mrV6/y+OOP651uIpu/hSmMGTOGY8eOsXbtWrP/ger7779n4sSJxMfHM2vWLJnyFCWiffAtWwYjRxq20/svZ2dYsABMeMbgrYra/H3w4EEefPDB27q4Ro0ayeZvYXK//fYbfn5+xMfHU6NGDa3LuSdFUXj66ad5/vnnGTFiBG3atOH1119n6NChWpcmLIC2wZeaCq6uxun0/qtSJUhKgrp1jTaEoiicOHHi5kKTWzd/u7m53bb5u0qVKkarRYjiKjyd5fPPP+eZZ57RupxiO3ToEL6+vhw6dIi0tDSZ8hTFpm3wdeoEu3aVbSFLcTk4QPv2sH27Qd4uMzOzyM3flStXvq2Lc3Fxkc3fwmz179+fatWq8emnn2pdSokFBweTnZ3N4sWLmT59Ojt37mTz5s0y5SnuSrvgi45Wg8+YU5z/5eysboIvwT4/RVE4c+ZMkZu/mzZtelsXZwnTREIU+uqrr5g6dSqxsbGan85SGleuXMHV1ZWff/6Z5s2by5SnKBbtgq9PH/jmm5Lt0ysre3vo3RtWrSry6ezs7CLv/F2uXLmbC01k87ewFoWns2zatIlWrVppXU6pff7556xcuZIdO3Zw4MABOnfuLFOe4q60Cb70dKhdW92gXoSHgZPAOqBHwWNvAvOAQUAkMAT4reB1FPy6U3HGdnKCc+dIy86+LeBSUlJo3LhxkZu/ZepEWJP8/Hy6du1K165dmTRpktbllEleXh5eXl5MmDCB3r17M336dH7//Xc2bdok/9+KImlzLMPu3eoxZHcIvmK9BfA4cA7IKcHXXdPpCGjUiM1wM9g6d+5MSEgIbm5usvlb2IS5c+eSk5PD+PHjtS6lzMqVK8f8+fPp378/L774IuPGjeO7775j2bJl+Pv7a12eMEPaBF90dJmv7SUVfK5KyYLPWVFY5O9PlVmz5KdBYZMSEhKYOXMmUVFRVrPoqmPHjrRr144PP/yQKVOmEBkZSefOnXnqqaeoX7++1uUJM6PNhrEdO0yzkrMI9nl5VI2Lk9ATNikrKws/Pz9mzZpFo0aNtC7HoD766CM+/fRTTp48yeOPP05wcDABAQFovVVZmB9tgu/sWU2GvencOW3HF0IjkydPxsXFhUGDBmldisE1aNCAUaNGMXbsWADGjx/PhQsXWLZsmcaVCXOjTfDl3H1y8r7Cl93yWHbBZ4MsuM7OvvdrhLAy27ZtY/Xq1Xz++edWO+MxduxYoqKi2L59O+XLlycyMpIJEyZw+vRprUsTZkSb4KtQ4a5PP1LweXfB5zxgb8GvmxhifEdHQ7yLEBbjypUrDB48mPDwcKvea+rs7MysWbMIDg5Gp9PJlKcokjbBV6fOXZ8uvMHIJ0BrwBWIA+4H+hY8NwYYDBQukZlZ8PvCRS93Vbt28WsVwgqMGDGC7t278/TTT2tditG98sorVKtWjaVLlwL/TnlGRERoXJkwF9oEn6+veoTYHXQHVgAtUYPsKtAN2AoURta3wBdAbsHvfy74/fl7je3goJ4YI4SNWLVqFXFxcXz44Ydal2ISdnZ2zJ8/n/fee4/Lly9Tvnx5IiIiGD9+vEx5CkCrDeybNqknt2RkmHxoqlSBr78GCzqMV4jSOnXqFF5eXmzevJmWLVtqXY5JDRs2jPLlyzN//nwApk2bxp9//snGjRut9hqnKB6zPLnFqApObqFqVdOPLYQJ5efn06VLF5566ikmTpyodTkmd/HiRR577DG2bdtG8+bNyc3NxcfHhzfeeIPXNLpFmTAP2kx1Vq0K3burZ2eakr099OwpoSdswpw5c9DpdIwbN07rUjRRo0YN3n77bd58800URbm5ynP8+PGkpqZqXZ7QkHaHVMfEqNf6TH13hp07wYIP5BWiOBISEujSpQtRUVFWt1G9JHQ6HZ6enkybNo2ePXsCMuUptOr4QL01kLf3XRe5GJSDgzqehJ6wcllZWfTv35+PP/7YpkMPwMHBgXnz5hEaGkpWwaWVCRMmkJaWJqs8bZjcgV0IKzN69GhOnTrFN998Ix1NgZdffhkvL6+bd6Io7Ijj4uKoV6+extUJU9M2+ACWLYORI4075ensDAsWgFzQFlZu69atDBo0iPj4eKpXr651OWbj2LFjtG7dmvj4eOoW/PA7depUdu/eLVOeNki7qc5CQ4ZAQABUrGic969YEQIDJfSE1bty5QpDhgxh2bJlEnr/0bhxY4KCgvRuwzRx4kTOnz9PZGSkdoUJTWjf8QEoCoSEwNKlhu38nJ3V0JszB+QnOmHFFEWhb9++1KxZ8+a+NaHv2rVruLq6smbNGtq1awdAfHw8Xbt2lSlPG6N9xwdqKM2dq05HVqpU9gUvDg7q+yxYoL6vhJ6wcqtWrSIhIcFmTmcpjUqVKvHhhx8yatQo8vPzAfVm1CNHjiQwMFDO8rQh5hF8oIbTa6/B4cPQvr3arZV0n5+9vfp17durC1lkelPYgFOnThESEsLKlSu577777v0FNqxfv344OjrqTW9OnDiRc+fOyZSnDTGPqc6ixMSoU5Tr1ql3U7h+veib1zo4oDg7k5WRQf5LL1Hx7bfVrRJC2IC8vDy6dOnCM888w4QJE7QuxyLExMTw4osvkpSURJUqVQCZ8rQ15ht8hdLTYc8eiI6G7dvV48ays9UwrF1bPXDa25vRa9ZQ/ZFHmDx5ssYFC2E6s2bN4scff+S3336jXLlyWpdjMfz9/alWrRoff/zxzcemTJnC3r172bBhg6zytHLmH3zFtGfPHgYOHMiRI0fkH62wCYVdSnR0NA8//LDW5ViUtLQ0mjdvzh9//EHTpk0ByM3NpXXr1gQHBzN48GBtCxRGZT7X+MrIx8cHOzs79uzZo3UpQhhd4ekss2fPltArhVq1ajFhwgRCQkJuPlZ4+6Jx48Zx5swZDasTxmY1wWdnZ8fgwYPlArWwCRMnTqRZs2YMGDBA61Is1siRIzl69CgbN268+ZinpycjRoyQVZ5WzmqmOgFSU1Nxd3fnzJkzsrpNWK1ff/2VwYMHy+ksBrBx40ZCQkJITEykQoUKAOTk5NC6dWvefPNNmfK0UlbT8QHUq1ePVq1a8cMPP2hdihBGcfnyZYYMGUJERISEngE899xzNGnSRG/Tf4UKFYiMjJQpTytmVR0fwMqVK1mxYgWbNm3SuhQhDEpRFPr06cNDDz3EvHnztC7HaiQnJ9OuXTsOHDjAQw89dPPxKVOmEBUVxU8//SQL5qyM1QVfZmYmdevW5eDBg9SpU0frcoQwmBUrVvDBBx8QExMjU/kGNnbsWC5dusSyZctuPlY45RkSEsKgQYM0rE4YmtUFH6h7dFxdXRk7dqzWpQhhECdPnsTb25tffvkFT09PrcuxOhkZGTRt2pQffvgBb2/vm4/v37+fp556iri4uJt3dRCWz6qu8RUqXN1phZkubFBeXh6DBg0iNDRUQs9I7r//fqZPn05wcLDe9w1PT0+GDx9OUFCQfD+xIlYZfB06dCArK4vY2FitSxGizGbPno2iKIwZM0brUqza4MGDyc3NZeXKlXqPT5o0idTUVJYvX65RZcLQrHKqE9QL0xcvXmTBggValyJEqRVOtUVHR9OwYUOty7F6u3fv5tVXXyUpKYlKlSrdfFymPK2L1Qbf8ePHad26NampqTg6OmpdjhAlduPGDby9vRk/frxsVDehAQMGUL9+fWbMmKH3+HvvvUdMTAw//vijrPK0cFY51QnQqFEj3Nzc2LBhg9alCFEqEydO5LHHHsPPz0/rUmzKzJkz+fzzzzl27Jje44VTnl9++aVGlQlDsdqODyAiIoL169fz/fffa12KECWyZcsWXnvtNeLj43nggQe0LsfmzJgxg+joaNatW6f3eOGU5/79+2W7lAWz6uD7559/qF+/PsnJydSsWVPrcoQolsuXL+Ph4UFERARdu3bVuhyblJWVhZubG59//vltfwfvvvsusbGxMuVpwax2qhOgcuXKvPTSS6xatUrrUoQoFkVReP311+nVq5eEnoacnJyYPXs2wcHB5Obm6j03efJkTp8+LVOeFsyqgw9g0KBBfPHFF1qXIUSxrFixgkOHDvHBBx9oXYrN6969O3Xq1GHRokV6jxee5TlmzBjOnj2rUXWiLKx6qhMgPz+fhx9+mB9//BEPDw+tyxHijk6ePImXlxe//vqr/Fs1EwcPHqRTp04cPnyYGjVq/PtEejpfDh+OXWws/evVw+7sWcjJgQoVoE4d8PUFb29o2xaqVtWsflE0qw8+gLfeeovMzEzmzJmjdSlCFCkvL4/OnTvz/PPPM27cOK3LEbcYNWoUubm5aucXEwOzZ8P69SgVKqDLyKB8UV/k4ADOzmoY9ugBoaHg5WXiysWd2ETwJScn07FjR1JTUylfvsh/pkJo6sMPP2TTpk1s3bqVcuXKaV2OuMWVK1d40sWF3U2acF9iImRlQX5+8d/A3h6cnNQOcMUKqFfPeMWKYrH6a3wAjz76KE2aNGHz5s1alyLEbeLi4pg9ezZffPGFhJ65URSqrVvHvhs3cIqNhczMkoUeqK/PzIRdu8DVFZYtA+vvN8yaTQQfyCIXYZ5u3LiBn58fc+bMkSPJzI2iQEgIjByJfWYmdjpd2d5Pp4Pr12HkSPV9Jfw0YxNTnQDp6ek0bNiQY8eOyZ2rhdkIDg4mLS2Nr776SvaEmZPC0AsLU8PK0CpWhIAAmDvX8O8t7slmOr6qVavy3HPP8fXXX2tdihCAejrLd999x6JFiyT0zE1EBCxdapzQA/V9lyxRpz2FydlMxwewefNm3nnnHaKiorQuRdi4y5cv4+7uzhdffEGXLl20LkfcKjVVvRZnrNC7VaVKkJQEcscHk7KZjg+gW7dupKamcujQIa1LETZMURSCgoL43//+J6Fnjvz8IDvbNGNlZUH//qYZS9xkU8FXrlw5BgwYIItchKa+/PJLkpKSbrvtjTAD0dHqR1kXshSXTqeOFxNjmvEEYGNTnQCHDh2iW7dunDp1SpaOC5M7ceIE3t7ecjqLuerTB775puRbFsrC3h569wY5U9hkbKrjA3jssceoU6cOv/76q9alCBuTl5fHwIEDGTdunISeOUpPh++/LzL0ngDsgPf+83h0weMVgEsFj6UAlQoe9yzOuPn5sG6dOr4wCZsLPoDBgwcTGRmpdRnCxsyaNYty5coxevRorUsRRdm9Wz1rswiFtwL+6j+PF/ZozwLVAV3Ba0t8hdDREfbsKelXiVKyyeDr06cPmzZt4urVq1qXImxEXFwcc+bMkdNZzFl0tHrCShFeBRyBZGBfwWP5wOqCXw8o+DwdSABCSzr29evq+MIkbDL4qlevTpcuXVizZo3WpQgbcOPGDfr378/cuXNp0KCB1uWIO9mx446LWqoBzxX8urDr+w04B1QBXkSd9nwf+BhwLenYOh1s317SrxKlZJPBBzLdKUxn/PjxeHh40K9fP61LEXdzj3vrFU53rgYU/g3AV4G8gue7AcNLO/65c6X9SlFCNht8zzzzDCkpKfz1119alyKs2C+//ML69etZuHChnM5i7nJy7vr080BV4DSwFVhb8LgfareXDFwBXgA+KXjueMHvi8VUeweF7QZf+fLl6devH8uXL9e6FGGlLl26xGuvvUZERATVqlXTuhxxL3dY2FLIEbW7AxgGpAMNUVd8Fu4J2wNsAOILfp9R8PticXQsdqmibGw2+ECd7ly+fDn5ptyzI2xC4eksvXv3ltNZLEWdOvd8SeF0Z0rB536o2xY6oYZf4UdEwfMe/BuK91S7dnFfKcrIpoPPw8ODatWqsV0uKgsDW758OcnJyUyfPl3rUkRx+fqqd06/i46oXV4hvzu9sKQcHKBTJ0O9m7gHmw4+kPv0CcM7fvw4Y8aMYcWKFTg5OWldjigub29wdr7rS+yAE/zb2T12h9cNLnh+f3HHrlhRHV+YhM0dWfZff//9N48++iinT5+mcuXKWpcjLFxeXh6dOnWiR48ehIaWeDeX0FJ6ujrdmJVl+rGdnNRVnVWrmn5sG2TzHV/NmjV54oknWLt27b1fLMQ9fPTRR5QvX56QkBCtSxElVbUqdO+unp1pSvb20LOnhJ4J2XzwgUx3CsPYt28fc+fO5YsvvsDe1N88hWGMGaN2X6bk5AQyO2BS8n8n8MILL5CYmMjx48e1LkVYqMzMTPr3788nn3xC/fr1tS5HlNJVFxcuNW58z0UuhqKzs0Px9oZWrUwynlBJ8AGOjo706dOHL7/8UutShIUaP348LVq0kNNZLNgPP/yAm5sbizt0MNmeumw7O75/9dV7v1AYlM0vbikUHR1Nnz59SElJkRM2RIn8/PPPBAQEEB8fLxvVLVBaWhqjRo1i3759hIWF4evrC8uWwciRdzy02iCcnTkeGkqbzz8nPj6ehx56yHhjCT3S8RXw8vLCycmJP/74Q+tShAW5dOkS/v7+REZGSuhZGEVRWL58Oe7u7jRq1IiEhAQ19ACGDIGAAHWbgTFUrAiBgTSaOhV/f3+GDRuG9CCmIx3fLT766COSk5MJCwvTuhRhARRF4ZVXXqFRo0Z8/PHHWpcjSuDkyZMEBQWRlpZGeHg4LVu2vP1FigIhIbB0qWE7P2dnCAyEOXPAzo7s7GxatmzJW2+9Rd++fQ03jrgzRdx05swZpWrVqsr169e1LkVYgIiICOXxxx9XsrKytC5FFJNOp1Pmz5+vVK9eXZkxY4aSk5Nz9y/Iz1eU8HBFqVRJURwcFEWNw9J9ODio7xMeftswUVFRSs2aNZXz588b6U8ubiUd3388++yz+Pn50b9/f61LEWbs+PHjtG7dmm3btvH4449rXY4ohkOHDjF06FDKlStHWFgYTZs2Lf4Xp6aCn596s9isLCjJ+b729uqWBW9vWLkS6tYt8mWTJk0iKSmJtWvXyjoDI5NrfP8he/rEveTl5TFgwAAmTpwooWcBcnJymDZtGr6+vgwYMIAdO3aULPQA6tVTbxS7Ywf07q0GWZUqd9724OCgPu/kpL5+xw716+8QegDvvvsuR44cYfXq1Xd8jTAM6fj+48aNG9StW5eEhATq1aundTnCDM2YMYOtW7eyZcsW2ahu5qKjo/H396dBgwYsWrTIcHss09Nhzx61A9y+XT1uLDtb3QZRu7Z64LS3N7RpU6ITWaKjo3nhhRdISEigVq1ahqlV3EaCrwhBQUE8/PDDTJw4UetShJmJjY3l2WefJTY2Vjaqm7HMzEzeeecdVqxYwZw5c+jbt6/FTB9OnDiRI0eOyJSnEcmPq0UonO6UnwnErTIzM/Hz82P+/PkSemas8LrruXPnSExMpF+/fhYVIO+9955MeRqZdHxFUBSFpk2b8uWXX+Lj46N1OcJMvPHGG1y5coWVK1dqXYooQnp6OmPGjOGXX35h0aJFPP/881qXVGpRUVG8+OKLMuVpJNLxFcHOzo5BgwYRGRmpdSnCTGzevJkff/yRzz77TOtSRBHWrVuHm5sbjo6OHDhwwKJDD6B169YMGTKE4cOHy8yTEUjHdwenTp2iRYsWnDlzRm4mauMuXryIh4cHK1as4Mknn9S6HHGL8+fPM3LkSBISEggLC6Njx45al2QwWVlZtGzZknfffZfevXtrXY5VkY7vDho0aICnpyc//vij1qUIDSmKQlBQEP369ZPQMyOKohAZGYm7uzsuLi7Ex8dbVegBODk5ERkZSXBwMGlpaVqXY1Wk47uLL7/8kq+//poNGzZoXYrQSEREBJ988glRUVE4mujEfnF3x48fJygoiEuXLhEeHo6np6fWJRnVhAkT+Ouvv/j2228tapGOOZOO7y5efvll/vzzT86fP691KUIDx44dY9y4caxYsUJCzwzk5eXxySef4O3tTdeuXdm7d6/Vhx6oqzwPHz7MmjVrtC7FakjHdw9DhgyhefPmhModkm2KTqejU6dO9OrVi5CQEK3LsXkHDx7E398fJycnli5diouLi9YlmdTevXvp3r07CQkJ1KxZU+tyLJ50fPcwePBgIiMjZWWVjfnwww9xcnIiODhY61JsWk5ODlOmTKFTp04MGTKEbdu22VzoAfj4+DBo0CBZ5WkgEnz30LFjR65du0ZcXJzWpQgTiYmJYd68eURGRsqRZBrau3cvLVu2JDY2lri4OIKCgmz672PKlCkcPHiQb775RutSLJ7t/isqJnt7ewYOHCgHV9uIwtNZFixYIGe1auT69euEhITQo0cP3n77bb7//nv5u+DfVZ6jRo3i77//1rociybX+Irh6NGjtGnThjNnzlChQgWtyxFGNGLECK5evcqKFSu0LsUm/frrrwQGBtKhQwfmzp1L9erVtS7J7IwfP55jx45J51cG0vEVwyOPPEKzZs3YuHGj1qUII9q0aRM//fQTn376qdal2JwrV67w2muv4e/vz8KFC1m+fLmE3h1MmTKFAwcOyCrPMpDgKya5T591u3jxIkOHDuWLL76gagluIyPKbu3atbi5uVGxYkUOHDjAM888o3VJZs3JyYmIiAiZ8iwDmeospoyMDBo0aMBff/3Fgw8+qHU5woAUReHll1/GxcWFjz76SOtybMa5c+cYMWIEhw8fJjw8nHbt2mldkkUZN24cJ06ckM6vFKTjK6b777+fF154ga+++krrUoSBRUREcPz4caZNm6Z1KTZBURTCw8Px8PDAzc2NuLg4Cb1SmDp1KomJiXKtrxSk4yuBLVu2MGHCBGJjY7UuRRjIsWPH8PHx4bfffqN58+Zal2P1jh49SmBgIBkZGYSHh+Pu7q51SRZtz5499OjRg8TERJmJKgHp+Eqgc+fO/P333yQmJmpdijAAnU7HgAEDmDx5soSekeXl5TF79mx8fHx47rnn2L17t4SeAbRp04aBAwcyYsQIrUuxKBJ8JVCuXDkGDBggi1ysxMyZM7nvvvsYNWqU1qVYtcTERNq2bcuGDRvYu3cvoaGhODg4aF2W1ZgyZQoJCQky5VkCMtVZQkeOHKFTp06cPn1a/ue1YDExMTz//PPExsbK5mgjyc7OZvr06SxevJgZM2bg7+8vdxcwkj179tCzZ08SEhJkyrMYpOMroaZNm/Lwww/z888/a12KKCU5ncX4/vzzT1q0aEFiYiL79+9n6NChEnpG1KZNG/z8/HjjjTe0LsUiSPCVguzps2xjx47F29ub//3vf1qXYnWuXbvGqFGjeOWVV5g6dSrfffcdderU0bosmzB16lTi4+P59ttvtS7F7MlUZylcuXKFRo0acfz4capVq6Z1OaIENm7cyPDhw4mPj6dKlSpal2NVfv75Z4KCgnjyySeZPXs2DzzwgNYl2Zzdu3fz8ssvy5TnPUjwlVLv3r3p1KkTw4YN07oUUUwXLlzAw8ODr776Cl9fX63LsRqXLl1i9OjR7Ny5k88//5ynnnpK65Js2tixYzl16hSrV6/WuhSzJVOdpSTTnZZFURQCAwMZMGCAhJ6BKIrCmjVraN68OdWqVSMxMVFCzwzIlOe9ScdXSjqdjvr16/Pbb7/h6uqqdTniHsLDw1mwYAF79+7F0dFR63Is3pkzZxg+fDgpKSmEhYXRtm1brUsSt5Apz7uTjq+UHBwc8PPzk67PAhw9epQJEyawcuVKCb0yys/PZ8mSJXh6etKiRQv27dsnoWeG2rZtS//+/Rk5cqTWpZgl6fjKoPAk+ZMnT1KuXDmtyxFF0Ol0PPHEE/Tu3Zvg4GCty7FoKSkpBAQEkJmZSXh4uJx2Y+Zu3LiBp6cnM2bMoFevXlqXY1ak4yuD5s2b89BDD7F161atSxF38MEHH1CxYkX5ybcMdDods2bNok2bNrz00kv8+eefEnoW4L777iMiIoI33niDixcval2OWZGOr4wWLFjAnj17WLlypdaliP+Ijo7mhRdeYN++fdStW1frcixSfHw8/v7+VK1alSVLltC4cWOtSxIlNGbMGFJTU/n666+1LsVsSMdXRn379mXDhg1kZGRoXYq4xfXr12+eziKhV3JZWVlMnjyZbt26MWLECLZs2SKhZ6GmTZtGXFwca9eu1boUsyHBV0Y1atSgc+fOckCsmRk7diw+Pj5yOksp/PHHH3h6enLkyBHi4+MZMmSIHDdmwWTK83Yy1WkA33//PR9//DG///671qUIYMOGDYwYMUJOZymhf/75h4kTJ7Ju3ToWLFjAyy+/rHVJwoBCQ0M5e/as3Ewb6fgM4tlnn+XIkSMcPXpU61Js3oULFwgICGD58uUSeiWwceNGmjdvzo0bNzhw4ICEnhV6//332bdvH999953WpWhOOj4DefPNN6lSpQpTpkzRuhSbpSgKPXv2xNXVlZkzZ2pdjkW4ePEib775Jn/++SdLliyha9euWpckjGjXrl288sorJCYmUqNGDa3L0Yx0fAZSeIRZfn6+1qXYrPDwcE6ePMnUqVO1LsXsKYrCV199RfPmzalVqxaJiYkSejagffv29OvXz+a390jHZyCKouDh4cH8+fPp1KmT1uXYnJSUFNq2bcv27dtxc3PTuhyzlpqayrBhwzhx4gTh4eG0bt1a65KECWVmZuLp6cnMmTNtdkpbOj4DsbOzY/DgwXKEmQZ0Oh0DBgzgrbfektC7i/z8fBYvXkyLFi3w9vYmNjZWQs8GOTs72/wqT+n4DCgtLQ1XV1dOnz5NpUqVtC7HZkybNo3ff/+dzZs3Y28vP8sVJTk5mYCAAHJycggLC5MfEASjR4/m/PnzrFq1SutSTE6+SxhQrVq1aN++vayaMqGoqCg+/fRTIiIiJPSKkJuby8yZM2nXrh0vv/wyf/zxh4SeANRVnjExMaxbt07rUkxOvlMYmEx3mk7h6SyffvqpnM5ShLi4OHx8fNi2bRvR0dEEBwfLYeripsIpzxEjRnDp0iWtyzEpmeo0sOzsbOrWrUtsbCwNGzbUuhyrNmzYMDIzM+UHjf+4ceMGU6dOJTw8nFmzZjFw4EA5eUXckS1OeUrHZ2COjo7873//48svv9S6FKv2008/sXnzZubPn691KWZl586deHh4cPToURISEhg0aJCEnrgrW5zylI7PCKKioujfvz/JycnyTccI/v77bzw9PVm9ejUdO3bUuhyzkJGRwfjx4/nxxx/59NNP6dGjh9YlCQvyxx9/8L///Y/ExESqV6+udTlGJx2fEXh7e+Pg4MCff/6pdSlWR1EUAgICGDRokIRegZ9++onmzZuj0+k4cOCAhJ4osQ4dOtC7d29GjRqldSkmIR2fkcycOZNjx46xZMkSrUuxKkuXLmXRokXs2bOHChUqaF2Opi5cuEBwcDB79+5l6dKldO7cWeuShAXLzMzEw8ODWbNmWf0PT9LxGcmAAQP49ttvuXHjhtalWI2UlBQmTZrEihUrbDr0FEVh5cqVNG/enLp165KYmCihJ8qscJXn8OHDrX6Vp3R8RvT0008zePBg+vbtq3UpFk+n09GhQwf69etnM9MxRTl16hSvv/46qamphIeH4+3trXVJwsqEhITw999/s3LlSq1LMRrp+Iyo8OBqUXYzZsygSpUqvPHGG1qXoon8/Hw+++wzWrVqRbt27YiJiZHQE0Yxffp0oqKiWL9+vdalGI10fEaUmZlJvXr1SExMlA3WZbB3715eeukl4uLiqFOnjtblmFxSUhIBAQHk5+cTFhZGs2bNtC5JWLnff/+d3r17F73KMz0ddu+G6GjYsQPOnoWcHKhQAerUAV9f8PaGtm2halUtyr8nCT4jCwgIoEmTJowfP17rUizS9evXadGiBTNmzOCVV17RuhyTys3NZdasWcyZM4f33nuP4cOHy7FswmTefPNNLl68yIoVK9QHYmJg9mxYv14NucxM0Olu/0IHB3B2VsOwRw8IDQUvL1OWfk8SfEb2xx9/EBgYyMGDB2VPXym8/vrrZGVlERkZqXUpJhUbG4u/vz+1a9dm8eLFcgqQMLnCVZ6R779P+0WL1A4vKwtKcs9Re3twclI7wBUroF494xVcAvLjo5G1b9+enJwcYmJitC7F4vz444/8/PPPNnU6S2ZmJuPGjeO5554jNDSUjRs3SugJTTjfdx+7AwJo+9prsGuX2uGV9Ebb+fnq1+3aBa6usGwZmEGvJcFnZHZ2dgwaNMjmOpay+vvvvwkMDGT58uXcf//9WpdjEtu3b8fDw4NTp06RmJjIgAEDZJZAaENRICSEGlOmYH+nKc2S0Ong+nUYORJCQjQPP5nqNIETJ07g5eXFmTNncHR01Locs6coCt27d8fNzY0PPvhA63KM7urVq4wbN46NGzfy2Wef8dJLL2ldkrBlBaFHWJgaVoZWsSIEBMDcuYZ/72KSjs8EHn74Ydzd3fnpp5+0LsUihIWFkZqaypQpU7Quxeh++OEH3NzcsLOz48CBAxJ6QnsREbB0qXFCD9T3XbJEnfbUiHR8JvLFF1/w7bff8uOPP2pdiln766+/aNeuHTt37rTqZftpaWmMGjWKffv2ERYWhq+vr9YlCQGpqeq1OGOF3q0qVYKkJNBgq5d0fCbSq1cvfv/9d9LS0rQuxWzpdDoGDBjAO++8Y7WhpygKy5cvx93dnYcffpiEhAQJPWE+/PwgO9s0Y2VlQf/+phnrP6TjM6HBgwfj4eFBSEiI1qWYpSlTprB79242btxolfvVTp48SVBQEOfPnyc8PJxWrVppXZIQ/4qOhk6d1FWYpuLsrG6CN/E+P+v77mLGZHXnne3du5dFixaxbNkyqwu9vLw8FixYQKtWrfD19SU6OlpCT5if2bPVLsyUsrJgzhzTjokEn0n5+vpy9epV9u/fr3UpZuXatWv4+fnx2WefWd2RZIcOHaJjx46sWbOGXbt2MXHiRMqXL691WULoS0+H77+/4z69hwE7YP0tj71Z8Nhg4DjwSsHrnIC6QCBw5V7j5ufDunXq+CYkwWdC9vb2DBw4UA6u/o/Q0FA6dOhAr169tC7FYHJycpg2bRq+vr74+fmxY8cOmjZtqnVZQhRt9271GLJSOgn8BDQD+gEZwFJgaHG+2NER9uwp9dil4WDS0QQDBw6kffv2fPTRR/KTP+rpLFu2bLGqLjg6Ohp/f3/q169PbGwsDRo00LokIe4uOrpM1/YeRe36ahf8/glgCLClOF98/bo6/jPPlHr8kpKOz8SaNGmCi4sLmzZt0roUzaWlpVnV6SyZmZmMGTOGF198kQkTJvDTTz9J6AnLsGNHmU5nqcO/oQeQU/C5WCdz6nSwfXupxy4NCT4NDB482OanOxVFYejQoQwZMoQOHTpoXU6Zbdu2jccff5xz586RmJhIv3795LgxYTnOnjXYWx0BJqOGy0fF/aJz5ww2fnHIVKcGXn31VcaMGcOlS5duv9eVjVi6dClnz55l7dq1WpdSJunp6YwZM4ZffvmFRYsW8fzzz2tdkhAll5Nz16fvK3zZLY8V7vZzvuWxaOA51EUty4AXiju+qfYOFpCOTwNVqlThueee46uvvtK6FE0kJyczefJkVqxYQYUyXFDX2rp163Bzc6NChQocOHBAQk9Yrnv8f/hIwefdBZ/zgL0Fv25S8HkL0Bm4BqwFBpVkfBOfYSzBpxFbne7Mzc1lwIABvPvuuxZ7Osv58+d59dVXmTBhAl999RULFy60imuUwjZdu3aN3AcfvOtrAgo+fwK0BlyBOOB+oC9wELW7uwa0BH5D3e7wJnC5OEXUrn3v1xiQBJ9GunTpwrlz5zh48KDWpZjU9OnTeeCBBxgxYoTWpZSYoihERkbi7u6Oi4sL8fHxPPHEE1qXJUSx5eXlkZCQwNKlSwkICMDd3Z1atWqxw85OvXP6HXQHVqCGWhJwFegGbEVd1HKBf6dB/wTm3fKRca+iHBzUE2NMSI4s09CECRPIz8/no4+KfQnYou3Zs4cePXoQFxdHbRP/hFdWx48fJygoiIsXLxIeHk6LFi20LkmIu1IUhdTUVKKioti7dy979+5l37591K1bFx8fH1q3bo2Pjw/u7u5U2LoV+vSBjHvGlOFVqQJff23S7QwSfBo6fPgwXbp04dSpUzjc5acta3Dt2jVatGjBhx9+yMsvv6x1OcVWeNzY+++/z9ixYwkNDbX6vythmTIyMoiJiWHv3r03w06n0+Hj43Pzw8vLi2rVqt3+xenp6nSjqY8sA3ByUld1Vq1qsiEl+DTm4+PDe++9x7PPPqt1KUYVGBhIbm4uERERWpdSbAcPHsTf3x9HR0eWLl3Ko48+qnVJQgDqnUwSExP1Qu7kyZN4enrqdXMNGzYs/raaPn3gm2/ueGyZUdjbQ+/esGqV6cZEgk9zCxcuZOfOnXz99ddal2I0P/zwA2+++Sb79++3iEUgOTk5fPDBB3z66ae8//77BAQEWN3B2cJyKIrCyZMn9aYs9+/fT8OGDW8GnI+PD82bNy/baVAxMeDra/q7M+zcCSY+tF2CT2OXL1+mcePGnDhxgqombPVNJS0tDU9PT7799lvat2+vdTn3tHfvXvz9/WnUqBGLFi2iXr1inT0hhMGkp6cTHR2t183Z29vrTVm2atWKKlWqGH7wTp1g164yneJSbA4O0L69yU9tAQk+s/Dqq6/StWtXgoKCtC7FoBRF4cUXX8TDw4Pp06drXc5dXb9+nbfeeouvv/6auXPn0rt3bzl5RRhdTk4OCQkJeiF35swZWrZsqTdlWa9ePdP8e7SRO7DLVXozMGjQIGbMmGF1wbdkyRLOnz/Pd999p3Upd/Xrr78SGBhI+/btSUxMpEaNGlqXJKyQoigcO3ZMb8oyISGBRx55hNatW9OhQwdCQ0N57LHHtFtAVa8ezJ8PI0cad8rT2RnmzdMk9EA6PrOQm5tL/fr12blzp9UsoEhOTqZ9+/b8/vvvuLq6al1Oka5cuUJoaChbt25l8eLFVr/ASJjW5cuXb4ZcVFQUUVFRODo63pyubN26Na1ataJy5cpal6pPUSAkBMLCjNP5VawIAQEwd67h37uYJPjMRGhoKE5OTmY/JVgcubm5tG/fnkGDBpntRvW1a9cycuRIevXqxYwZM8zvm4+wKNnZ2ezfv1+vm0tLS8PLy+vmdGXr1q2pq1GHU2KF4bd0qWE7P2dnCAxU77qu4aUECT4zkZCQwAsvvMCJEycsfgXhu+++S1RUFBs3bjS762Tnzp1jxIgRHD58mLCwMItYcCPMi6IopKSk3Ay4qKgoDhw4gIuLi14316xZM8qVK6d1uaWnKBARAcHB6v6+six4cXBQ9+vNmwevvWa4GktJgs+MtGzZklmzZtGlSxetSym13bt307NnT7M7nUVRFCIiIpgwYQKBgYG89dZbODk5aV2WsAAXLly4OVVZGHSVK1fWC7mWLVtSsWJFrUs1jtRU8PNTbxablVWyfX729mrgeXvDypWaXdP7Lwk+MzJv3jxiYmL48ssvtS6lVK5du4anpyezZs2iZ8+eWpdz07FjxwgMDCQ9PZ3w8HA8PDy0LkmYqRs3brB//369bu7SpUt4e3vrTVk+9NBDWpdqejEx6hTlunXq3RSuXy+6C3RwUK/jZWdDz54wejR4eZm+3ruQ4DMjFy5cwMXFhdOnT1vkNaeAgADy8vJYtmyZ1qUA6nFj8+bNY8aMGUyYMIE333xTjhsTN+Xn55OcnKy3leDw4cO4urrqdXNNmza1+MsPBpWeDnv2qB3g9u3qcWPZ2WoY1q6t7gX09oY2bUx6DFlJSPCZmR49etC9e3eGDBmidSkl8v333zN69Gj2799vFqGdmJjI0KFDcXZ2ZunSpTRp0uTeXySsWlpamt7ik+joaB544AG9/XItWrTgvvvuu/ebCYsmwWdm1q1bxyeffMKOHTu0LqXYzp8/j6enJ9999x3t2rXTtJbs7GxmzJjBwoULmTFjBkOHDjW7BTbC+DIzM9m3b5/elOXVq1f1jvjy9vamZs2aWpcqNCDBZ2ZycnKoW7cuUVFRNGrUSOty7klRFF544QVatGjB+++/r2ktu3fvxt/fn0cffZTPPvvMcpaOizLJy8sjKSlJb8oyOTkZNzc3vSlLFxcX+SFIABJ8ZmnUqFFUr16dd999V+tS7mnx4sWEh4fz559/lu2A3DK4du0akydPZs2aNcyfP59XXnlFvsFZsbNnz+pNWcbGxlKzZk29KUtPT08cHR21LlWYKQk+MxQbG8urr75KSkqKWV9UP3LkCB06dOCPP/6gadOmmtTwyy+/EBQUhK+vL7Nnz6Z69eqa1CGM49q1a8TGxup1czdu3LhtylL+3kVJSPCZIUVRePzxx1m4cCFPPPGE1uUUKTc3l3bt2jFkyBCGDx9u8vEvX77M6NGj2bFjB4sXL+bpp582eQ3CsPLy8jh48KBeyB09ehR3d3e9bq5x48bS0YsykbXdZsjOzo7BgwfzxRdfmG3wTZs2jQcffJBhw4aZdFxFUfj2228JDg7m1VdfJTExkUqVKpm0BlF2iqKQmpqqN2W5b98+6tatezPkgoKCcHd3p0KFClqXK6yMdHxm6ty5czz22GOkpqaa3YkQhaez7N+/36Qbec+ePcvw4cNJTk4mPDyctm3bmmxsUTYZGRnExMTodXM6nU7vHnNeXl5Uq1ZN61KFDZDgM2PPPfcc/fr1w8/PT+tSbvrnn3/w9PRk9uzZ9OjRwyRjKopCWFgYkyZNYvjw4UyaNEkWLpgxnU5HYmKiXjd38uRJPD099aYsGzZsKFOWQhMSfGZszZo1LF26lC1btmhdyk1Dhw5FURTCw8NNMl5KSgqBgYFcu3aN8PBwHn/8cZOMK4pHURROnTqlt18uLi6Ohg0b6i1Aad68uWarfoX4Lwk+M5aVlUXdunXZv38/9evX17oc1q9fT2hoqElOZ9HpdHzyySfMnDmTSZMmERwcbNkn3VuJ9PR0oqOj9Q5strOz09sv5+XlRZUqVbQuVYg7kuAzc6+//joNGjRg0qRJmtZReDrLunXrjH5tLT4+Hn9/f6pUqcLSpUtp3LixUccTRcvJySEhIUFvyvLMmTO0bNlSb8qyXr16MmUpLIoEn5nbs2cPgwYNIikpSbNvLoqi8Pzzz9OqVSumTZtmtHGysrJ4//33WbJkCTNnzmTIkCHyDdVEFEXh+PHjelOW8fHxPPLII3pTlo899pgc9C0snvwLNnM+Pj6AGoBarWJcvHgxFy5c4J133jHaGLt27WLo0KE89thjxMfHm9W9/KzR5cuXb7vHnKOj482AmzFjBq1atTKLA8eFMDTp+CzAjBkzOHXqFIsXLzb52MY+neWff/5h0qRJrF27lgULFtCrVy+Dj2HrsrOziY+P1+vmzp8/j5eXl9495uRsU2ErJPgswOnTp/H09OTMmTMmvWt44eksr732mlE2qm/atInXX3+drl278vHHH8seLgNQFIWUlBS9/XIHDhzAxcVFbwFKs2bNZLGQsFkSfBaiW7duDB06lN69e5tszLfffpt9+/bx008/GfRa28WLFwkJCWHXrl0sWbKErl27Guy9bc3Fixf1Fp9ERUVRuXJlvcUnLVu2NLtDEITQkgSfhVixYgWrVq1i48aNJhnvzz//pFevXsTFxRnsdBZFUVi9ejUhISH07duXadOmyTfkEsjKyiIuLk4v5C5evIi3t7deN2fK03SEsEQSfBbi+vXr1KtXj0OHDhl94Ufh6Sxz5syhe/fuBnnP1NRUhg0bxvHjxwkPD7+5aEcULT8/n+TkZL0py8OHD+Pq6qoXck2bNjXrO3gIYY4k+CyIv78/rq6ujB071ujj2NnZERYWVub3ys/PZ8mSJbz99tuMHDmSCRMmyKHDRUhLS9ObsoyOjuaBBx7Qm7Js0aIF9913n9alCmHxJPgsyM6dOxk+fDiJiYlG29+2fv16xowZw/79+8t814Pk5GQCAgLIyckhLCwMNzc3A1Vp2TIzM9m3b5/elOXVq1dvu8dczZo1tS5VCKskwWdB8vPzcXFxYc2aNbRq1crg73/u3DlatGhR5tNZdDods2fPZtasWbz99tu88cYbNruCMC8vj6SkJL0py+TkZNzc3PSmLF1cXGSzvhAmIhvYLYi9vT0DBw4kMjLS4MGnKAr+/v4EBgaWKfTi4uLw9/enRo0aREdH06hRIwNWaf7Onj2rN2UZGxtLzZo1bwbckCFD8PT0lLtLCKEh6fgsTEpKCq1ataJDhw5ER0dz7tw5g3RTCxcuJDIykl27dpXqFP0bN24wdepUwsPDmTVrFgMHDrT6DubatWvExsbqdXM3bty4bcqyevXqWpcqhLiFBJ8FCQ8PZ8yYMfzzzz/k5eXh5OTEjRs3yvy+SUlJdOzYkV27dvHoo4+W+Ot///13hg4dioeHB/Pnz7fK5fR5eXkcPHhQr5s7evQo7u7uegtQGjdubPWBL4Slk6lOC1KrVi2ysrLIy8sD4MEHHyz1e504cQJfX18WLFjA1KlTmTZtWolDLyMjgwkTJvDDDz/w6aefmuzGtMamKApnzpzRW3yyb98+6tSpczPggoKCcHd3lxWqQlgg6fgszK+//kr37t3JzMykZcuWxMbGlup91q9fT9++fdHpdDz00EMkJSWVaDP5hg0bGDZsGE8//TSzZs2iatWqparDHGRkZBATE6PXzel0Or3FJ97e3nKkmhBWQjo+C9O1a1e2bdtG+/btKcvPLEeOHEGn06HT6fj7779xc3PjyJEj91x0ceHCBYKDg9m7dy+RkZF07ty51DVoQafTkZiYqBdyJ0+exNPTEx8fH/r06cPcuXNp2LChTFkKYaUk+CyQj48PW7du5Z9//oH0dNi9G6KjYccOOHsWcnKgQgWoUwd8fcHbG9q2hVu6spiYGHQ6HeXLl6dChQqMHDnyrtN2iqKwatUqQkNDGTBgAImJiTg7Oxv/D1sGiqJw6tQpvSnLuLg4GjZseHPKcuTIkTRv3rxUC3qEEJZJpjotVUwMzJ4N69erIZeZCTrd7a9zcABnZzUMe/SA0FDw8uLBBx/k0qVLhISE8Pbbb991qvL06dO8/vrrnD59mvDwcLy9vY31pyqT9PR0oqOj9e4xZ2dnpzdl6eXlRZUqVbQuVQihIQk+S5OaCn5+aoeXlQX5+cX/Wnt7cHICb2/COnXiiX79blvQkpGRQWhoKPPmzcPJyYnFixfz7rvvEhwczLhx48xmMUdOTg6JiYl63VxqaiotW7bU205Qr149mbIUQuiRqU5LoSgQEQGjRkF2dtHd3b3k56ud4a5dDI2JgQYNwMUFbgmGSZMmsWzZMrKzszl+/Dj5+fns3LmTZs2aGfAPUzKKonD8+HG9/XLx8fE0btwYHx8fOnTowOjRo3Fzc8PBQf5JCyHuTjo+S6AoEBICS5eqwWUozs4QEABz54KdHYcPH6Zly5ZkZWUBMHbsWGbOnGny0/8vX75MVFSU3pSlo6Oj3n65Vq1aUblyZZPWJYSwDhJ85q4w9MLC4Pp1w79/xYoQEIAyZw4eHh4kJibefOqRRx4hKSnJqF1UdnY28fHxelOW58+fx8vL62bItW7dmrp16xqtBiGEbZF5IXMXEWH4Tu9W16/DkiXE5eaSmJiIvb09NWrUoHbt2jRr1gydTmew4FMUhZSUFL0pywMHDuDi4oKPjw+dO3dmwoQJNGvWzGYPtRZCGJ90fOYsNRVcXY3T6f2HUqkSF3bu5EFPT4MtBrl48aLefrmoqCgqV66sN2XZsmVLuQu7EMKkJPjMWadOsGtX6RaylJSDA7RvD9u3l+rLs7KyiIuL0wu5ixcv4u3trbedwBrP8RRCWBYJPnMVHa0Gn7GmOIvi7KxugvfyAuDSpUv07NmTTp06MXXq1Jsvy8/PJzk5WW/K8vDhw7i6uuqFXNOmTU2+MEYIIe5FrvGZq9mz1X16ppSVBXPmwKpV7N+/n2eeeYaLFy9y+fJlvL29b3Zz0dHRPPDAAzcDrn///rRo0YL77rvPtPUKIUQpSMdnjtLToXbtIoPvCeB34F3gvVsejwZaA+WBycAGIAXIAVyBd4CXijO2kxNzx41j/AcfkJube/Phbt260aZNm5v3mKtZs2Yp/mBCCKE9CT5ztGkT9OkDGRm3PbUECAIeBY7c8ngI8AlquF0FzgLtCl6zB7W1jwJa3GvsKlXwc3DgqytXcHR0JCsri/LlyxMdHY27u3vZ/lxCCGEG5AKMOYqOvuO1vVcBRyAZ2FfwWD6wuuDXA4A5Bc9HAn8AjQEd8Ftxxr5+nRXBwVy+fJnVq1fz+uuv06hRI65du1bKP4wQQpgX6fjMUZcusG3bHZ9+GVgHjAFmAVuBrkAVIA01GAvlAQ8Dqajh+L/ijN+5M2zdWorChRDC/EnHZ47Onr3r034Fn1cDCvBVwe8Lu8FbjUYNvXaogVks584V95VCCGFxZFWnOcrJuevTzwNVgdOo3d7agsf9bnlNHuq1wHDAC/iJEvxlZ2cXu1QhhLA00vGZo3vc+scRtbsDGAakAw1RV3wCZAG9UEOvG+q1vWolGf8ed2EXQghLJsFnjurUuedLCru7lILP/YDCg8b8ge8BJ8AFeAt4E1hV3PFr1y7uK4UQwuJI8JkjX1/1CLG76Ija5RW6dZrzTMHnLGAhMK/g45fijO3goJ4YI4QQVkqu8Zkjb2/1+LAi9vEVsgNO3OG57WUZu2JFdXwhhLBSsp3BHN3l5Bajc3JSV3VWrWr6sYUQwgRkqtMcVa0K3buDqQ94treHnj0l9IQQVk2Cz1yNGaN2X6bk5AShoaYdUwghTEyCz1x5eanX2gx09/N7cnBQx2vVyjTjCSGERuQanzkz4R3YqVQJkpKgbl3jjyWEEBqSjs+c1asH8+erKzyNydkZ5s2T0BNC2AQJPnM3ZAgEBKjbDIyhYkUIDITXXjPO+wshhJmR4DN3dnYwdy4MHWr4zs/ZWQ3VOXMM+75CCGHG5BqfpVAUiIiA4GB1f59OV/r3cnBQV3DOmyednhDC5kjwWZrUVPDzU29Wm5UF+fnF/1p7ezXwvL1h5Uq5pieEsEky1Wlp6tWD7dthxw7o3VsNsipV7rztwcFBfd7JSX39jh3q10voCSFslHR8li49HfbsUTvA7dvV48ays9VbC9WurR447e0NbdrIiSxCCIEEnxBCCBsjU51CCCFsigSfEEIImyLBJ4QQwqZI8AkhhLApEnxCCCFsigSfEEIImyLBJ4QQwqZI8AkhhLApEnxCCCFsigSfEEIImyLBJ4QQwqZI8AkhhLApEnxCCCFsigSfEEIIm/J/gMC7WImQV+UAAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "print(IDC_with_markov(G, 'V2', 'V3', ['V4']))"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "exBrZ7nJQaER",
        "outputId": "a9c522bc-dec2-4a87-e86e-fe68898f4de3"
      },
      "execution_count": 15,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Based on the current setting the query is not identifiable.\n",
            "The query becomes identifiable if ['V5'] are observed.\n",
            "The query is still not identifiable with ['V1'] being observed.\n",
            "The query becomes identifiable if ['V5', 'V1'] are observed.\n",
            "(False, [['V5']])\n"
          ]
        }
      ]
    }
  ]
}