PCalgorithm_run_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## notebook detailing what the PC algorithm wrapper does: PCalgWrapper.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the PC algorithm is part of the Java Tetrad causal search suite,\n",
    "# repository here: https://github.com/cmu-phil/tetrad\n",
    "# and project website here: https://sites.google.com/view/tetradcausal\n",
    "# Specifically we use the Tetrad command line wrapper: https://github.com/bd2kccd/causal-cmd\n",
    "\n",
    "# We include in our repository the Tetrad command line version we used in Sanchez-Romero et al.(2023)\n",
    "# Importantly, it includes the PC modification where we remove some of the orientation rules\n",
    "# that assume acyclicity in the true model (see Sanchez-Romero et al.(2023) for details).\n",
    "\n",
    "# we provide here an example of a detailed PC algorithm run\n",
    "\n",
    "# Ruben Sanchez-Romero 2023"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import os\n",
    "import numpy as np\n",
    "import h5py\n",
    "import subprocess\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### single PC run\n",
    "The idea of the PCalgWrapper is to create a string containing PC run commands and call it from system command line. Then transform the output (list of edges) into matrix form."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### build the string with PC commands"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "java -Xmx20G -jar DirectedActflow_release/tetrad_causal-cmd-1.1.3/causal-cmd-1.1.3-jar-noMeekRules234.jar --thread 360 --algorithm pc-all --stableFAS --concurrentFAS --test fisher-z-test --alpha 0.01 --depth -1 --colliderDiscoveryRule 3 --conflictRule 3 --dataset DirectedActflow_release/temp/example_temp.txt --no-header --data-type continuous --delimiter comma --out DirectedActflow_release/temp --prefix PCgraph_temp --skip-latest  \n"
     ]
    }
   ],
   "source": [
    "# when using Tetrad is recommended to use full paths\n",
    "# adjust paths as necessary!!\n",
    "\n",
    "current_dir = os.getcwd()\n",
    "\n",
    "# 1. PC requires the dataset [timepoints x variables] to be passed as a comma-separated txt file\n",
    "\n",
    "# download an example hcp bold data subject\n",
    "file_name = f'{current_dir}/hcpData/100206_hcpData.h5'\n",
    "f = h5py.File(file_name, 'r')\n",
    "bold_data = np.asarray(f['rest_data']).T # transpose to get [timepoints x variables]\n",
    "# convert the data into a comma-separated txt file and save temporarly\n",
    "input_data = f'{current_dir}/temp/example_temp.txt'\n",
    "np.savetxt(input_data, bold_data, fmt='%0.6f', delimiter=',')        \n",
    "\n",
    "# 2. Define output file dir and name\n",
    "\n",
    "out_dir = f'{current_dir}/temp'\n",
    "out_name = 'PCgraph_temp'  # this is used by Tetrad and needs to be only the file name, no path, no extension!\n",
    "\n",
    "# 3. Define the path where the command line Tetrad Java file is located\n",
    "\n",
    "tetrad_path = f'{current_dir}/tetrad_causal-cmd-1.1.3/causal-cmd-1.1.3-jar-noMeekRules234.jar'\n",
    "        \n",
    "# 4. Define the PC algorithm parameters. See Sanchez-Romero et al.(2023) for details.\n",
    "\n",
    "algorithm = 'pc-all'\n",
    "alpha = 0.01\n",
    "test = 'fisher-z-test'  #z-test conditional independence\n",
    "depth = -1  # maximum size of the conditioning set, -1 is all the variables\n",
    "# Use collider MaxP heuristics from Ramsey (2016) https://arxiv.org/abs/1610.00378\n",
    "# Collider discovery: 1 = Lookup from adjacency sepsets, 2 = Conservative (CPC), 3 = Max-P\n",
    "collider_discovery_rule = 3                \n",
    "# Collider conflicts: 1 = Overwrite, 2 = Orient bidirected, 3 = Prioritize existing colliders\n",
    "conflict_rule = 3                     \n",
    "# input-output parameters\n",
    "data_type = 'continuous'\n",
    "delimiter = 'comma'\n",
    "# others\n",
    "thread = '360'  #for parallelization, usually equal to the number of variables\n",
    "\n",
    "# 5. Form the string to call it from system\n",
    "\n",
    "call_PC = (\n",
    "          f'java -Xmx20G -jar {tetrad_path} ' #Xmx**G is used to allocate memory, increase or decrease accordingly\n",
    "          f'--thread {thread} '\n",
    "          f'--algorithm {algorithm} ' \n",
    "          f'--stableFAS '  # This runs the order independent PC version from Colombo and Maathuis (2014) https://jmlr.org/papers/v15/colombo14a.html\n",
    "          f'--concurrentFAS '  # to run in parallel when possible\n",
    "          f'--test {test} '\n",
    "          f'--alpha {alpha} '\n",
    "          f'--depth {depth} '\n",
    "          f'--colliderDiscoveryRule {collider_discovery_rule} '\n",
    "          f'--conflictRule {conflict_rule} '\n",
    "          f'--dataset {input_data} '\n",
    "          f'--no-header ' # do not use if your data contain headers, ie. 1st row with variable names\n",
    "          f'--data-type {data_type} ' \n",
    "          f'--delimiter {delimiter} '\n",
    "          f'--out {out_dir} '\n",
    "          f'--prefix {out_name} '\n",
    "          f'--skip-latest  '  # to avoid checking if there is a new Tetrad version\n",
    "          )\n",
    "# print to check it looks correct\n",
    "print(call_PC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### run call_PC from system command line"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PC graph for DirectedActflow_release/temp/example_temp.txt in DirectedActflow_release/temp/PCgraph_temp.txt\n"
     ]
    }
   ],
   "source": [
    "# run as is to avoid having temporary output printed in the notebook\n",
    "p = subprocess.call(call_PC,\n",
    "                    shell=True,\n",
    "                    stdout=subprocess.DEVNULL,\n",
    "                    stderr=subprocess.STDOUT)\n",
    "# if want to print the temporary output\n",
    "#p = subprocess.call(call_PC,shell=True)\n",
    "\n",
    "# print if there was an error in the PC run\n",
    "if p != 0:\n",
    "    print(f'{p}, error in the PC run')\n",
    "    print(f' copy call from above and run it from a terminal window')\n",
    "else:\n",
    "    print(f'PC graph for {input_data} in {out_dir}/{out_name}.txt')\n",
    "# TODO print the Tetrad error here. As of today, in case of error I copy the string and call it directly from\n",
    "# a terminal window (outside the notebook) to see the error."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### open the PC output and transform it into a connectivity matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(360, 360)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f78c1a134a8>"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQEAAAD8CAYAAAB3lxGOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABMTklEQVR4nO29d3gb55Xv/z0zg0KAHSwSi0hRpKguSpRYZCdxiyU7zpVbYsmJ4ziOldjO3mT37m7K3r2b/e3u3Xa3ZrOO4zhxiePYcZy4xL0ksaPerS6KpERSYu8gQWBm3t8fAxYQbQYYVM7nefhwMHjnnXeAwZm3nPM9xBiDgYHBwoVLdAMMDAwSi2EEDAwWOIYRMDBY4BhGwMBggWMYAQODBY5hBAwMFjgxMwJEtI2IzhJRCxF9K1bnMTAwiA6KhZ8AEfEAzgH4JIBOAAcA7GSMndL9ZAYGBlERq55AA4AWxlgrY8wN4OcAtsfoXAYGBlEgxKjeUgAdc153AmgMVthMFmaFPUZNMZgLy7KBc7rAZDniOojj4MnPgNDvBIjAsjJAoxOhD7JnABwBY2HKGcSMMQz1M8YK5++PlREICxHtArALAEyZedhS+RCkk2f9yvU+vAXFPzyI3i/Xo/CRPfFuZkIYuL8ZRb88BWl4RPe6+ZJlkNs6wDzuyCthAAaAlv9oQu5pQvHv+yGdPh/6mOnfPkV+WoPoeIe9cDHQ/ljNCTQD+C5jbKv39bcBgDH294HKZ+aVs6sLd2J8TSGmsnnkPrUwfuyxZPgLzch//ghklytkOc5qDVsmGEJZKTxLCsAfPhtxHQbx4x32wiHG2Kb5+2M1J3AAQA0RLSUiM4AdAF4OVpgbdkJu70Dm787DMip5d/Lw3OjXXgOV5P/qRNgfJlksOPPv6yI+h9jZBf7wWfTv3BBxHQaJJyZGgDEmAvgagDcBnAbwPGPsZMhjRBHSyCjsraPeHTIs3eOxaN4MQlUlLv/ZFlx8fi3YVXUxPVesEJZWAOTfx5bHxsIey9xurHgkfLlQyC4X8n+yB5zNBmHxoqjqMkgMMfMTYIy9xhhbzhhbxhj7O1UHyRLk42emK5jdBgCOV/50ZGRDMSaLZVT95SToD0d1rTtejG5YBFCEXyNjkI+d1qUdXE42zn99KfjcHF3qS2ZIEAIa3lQlJnMCWsmmfNZI14csw66qA0kysPf4zD7avBbOMhtsv9oX8bn53BxIo+OALEVch4ECn5sDaWQUSIJ7Kpa4Pt0A+9kBSOcuJLopmoj3nIBm+OXLMHlrQ8D3SBDQvy7DxwAAAOecgmXIM/vabge3boUmKy0NjxgGQCek4ZGwBmDsriYIFeVxalFssL6yP+UMQCiSxghI51uR8ev9fvuFxYswuW1jwOVB6dQ58L89DD43B+OfaQSVFKP9tnyQ2RyHFvsjlJZg6lObE3LuaJm4rRF8od8Ssu5kPbcX4sWO8AUN4kZSDAeycsvYVWX3Qjp1zu89EgRwuTnwrCgH9+HRwBVwPPj8XEj9A7FtaBjIZAaXaYc0NJTQdkQC78iHPDIKJoqJbkpQ+OqloCkPxI7ORDclJUnq4QA3PgXWegkggrCoeGa/sHgR+MICSIPDcOeFeLrLUsINAAAwjzuhBoDPy8PQF5sjOlYaGNRkAKY+tRn8quVhy5Eg6NbDIElO+/mGRJAURoDZLKCyxeCzsnDqb5YAUG6eoU9U4vQ/lIB4Hu7MpGhqUiMNDSHvyb3qDyACbV6rzKNoxPLawYA9t/lweXlwNi3VXH8gxLaLEDu7dKnLYJbk+GWNTUBqaVN8yznF0jNRRNbP92Lld7rBPG5kPzvn5m6K3MEl7VH5pBQql2D8zgZwZy8CbRH8sFSeR+rrg/WV/cZ3lsQkhRFgOTbwq5ZDHnei5se+XVKx6/LMdt9XmwGOh2zS119gISK2X0LmL/ZBGh1V5VgULbKJR9+DzeCLi2JzAo4H78iPTd1pTlIYARqZgHT6PJgohnTaKfrRAUCWwH1wJH6NS0HG7mpKdBP84D44gkUvXsDo1foMDebD52Rj9NqamNStF+5tm+OyAqOVpDACAFR1L5Nl5ppfWQPXpwP7NCQDeXu1de85u13VJF+0SD29sP9yHzr+YgtctzTo6gEqDQ3B/kLkTmPxwHasA3IMIkOjJXmMgAb4vLyZbbJYwNlscT2/dPq8Ms5NUuQsu/LUmfM5hYJMAqRsq99+tcdrZelPOzBRuPCGdOKV7uhCuMNAm9ZA/oT2YK7UMwJE6N8+O5tNtUshbvR9irGr6iCUl8W7ZUmDfPIsMvaew9BNtQAUQzm5vQHOOwLrukjDI37emAAw8OkVip+8zogXO5D/kz2Gp6bOsMOnwX3g/z2GIymchdTEDmhBWFQMeWwcstOpW53Jgnx1HYhBW8ATx0OoKAPjOWUVJgAjn2+C4912iFe69WmoSnhHPijTnl5ehETg1tb6BsAlAUntLKQX/bsURxmxu0cXA0AmMzjrnG4yEbisrKjrjQbuw6PaIx5lCWLbxaAGQCgtAT/F4m4AAIAy7Wj5cpnv55zqEIepRZmJboVq0soIFP30GDibDe5t+vjv08oqeJpWzbzmbDYM3rpGl7qTCfHyFWS9fDQx577YgWX/fBJUkUbDN1mC6a2DiW6FapLSCPQ+vAX88mUgk7pAIH5lDUAEeWICzO1GxgXFhZhfXQvPDfVgzesjaod8/Az43x6efe10IvfpNJQ+Ywxsaiphp5fGxoD+wYSdXy/cWzeBq1sVvmCSkZRGoOj7u+Fakgsym1SVn1iaC+KV2WYmipDOt4Kz2+GsykHfRgsGVyurB9K1G4O6yPLFRUm5vr4gYAzSQOobAfObByEfTb3UGklpBADA9M4h1eN6y2sH/HwI+u5eB+sr+1HyT7uReVl5T9h9EuxUS8A6pN4+5LyifWbVIDZM3ZxcIdmseX3svB0TTNIagWhxPLZHmRWvqkTGJcUtlk1NBXc4YgzyhKGJnyxMf2exoPdrWzQfY+oZAZuYRP+u5vSaxEQSG4GxHU2YunkzuPUrI66DTALGVxVBPnEGZLHo2LpI2mIGiDB0b7PquY6FTCy/s6L/2q35GLG1HfLYGAp+uMdPxVmoqsTUTcnVc9FC0hqB3NdOIeN3p8DOtEZcB5uagvXV/eCLi9D2lxsTKoJ58dubwK2tRcFLp2LqNTYNWSyQrt0Y8/PEkqG7Yv+d6SFrL13qRMbvU28uYJqojAARtRPRR0R0lIgOevflE9HbRHTe+z8i31PKyVa67ypmrXlHPgbvCy6mQRyHitcmYpLRRy1Ln+sF2rri1waZwTSQ2sOb3Kf2xPzzMvdNRl0HE8WUdkzToydwLWOsbo4n0rcAvMsYqwHwrvd1SFiWDfzyZT77Tv9pGfhpHXuOR9v/bQbVrw54vDQwqLihBoPnIZu9vupEPq6wsXCLDYR0tiUuIbvTMI9bncda07qUXNbSC3YkZDqMBUFUbsNE1A5gE2Osf86+swCuYYxdIaLFAH7LGKsNVU8OX8Aazdt8nvp8Xp4iX+31LxfKyyAPDEY9eTd6dxN4N4P9hX0QykoxvKUcmc9rUONJMzibDZBlI43YAiBWbsMMwFtEdMibYBQAihljV7zb3QCKAx86pxJZ9jEAvQ9tUbT65gSYiB2d2g0AEYTSEp9dmR1TyP5g1n1WMsNH13ChIU9MJL0BIIslKePw04VojcDVjLGNAG4C8DARfXzum0zpZgTsahDRLiI6SEQHPfAd9xfplH2Ys1gw0uyrcW9u7YFrvaJjOPixcuQ+dxhid48u5zOIDXyBA64NFTOvJ25rXNCGOxgjn4/M2U23KEIi+i6AcQAPQONwQO8oQrUIZaUY/Fg5si5OgnYfi/v5DSJDWFSMCw8uw7KfdEJsv5To5iQNwqLikA803YcDRGQnoqzpbQA3AjgBJfvwvd5i9wJ4KZL6h+5t1k0zLthTQ+zsQvbP96Gvzg62JbL4AoP4I3b3YNlPkiv3AGe1YuCByOTe9SLSHm3EPQEiqgLwK+9LAcDPGGN/R0QOAM8DWALgIoDPMsZCOoYH7AlwvG6iE8P3NIcO/JmbtiwJ9BUMUhQd71ldIQIY078nwBhrZYyt9/6tns48zBgbYIxdzxirYYzdEM4ABEXHDzP36T2gDauDe6AxBjCG0R2NSa0dGE+ExYsgVC7RflzlkgWRopxfXeufci4ZDQCA3odC91CS1mMwUgbuD+zbLWcIoDCJSrOf3QvTWHAx096vbUHfg4nt8sUNswnM5O9DQRYLer+2BZzdHvAwZjYBcfK98IHjY6aJGAjp5FlYfnMgbueLhqLvh3aTTjsjUPDkgYBLXrT7mKqlsLn6AfMpfvQg8s9MRfSEVAO7qg7C0orwBeOAeLED0nl/l22uogy5591BPeSkcxcSkiuQz8lGyzdXxL0XkkzfWaSklRGgDat1kSU/94MGDDzQDNqwGr0PbQFfWw0+Nwf9X6yHpX0Acl9s8h4K57og9/aHL5hApHMXYH5TX9Ucbv3KqD03paEh1PxHG+ShYc3HkiBE3MNLhe8sHGllBNyFGbrUU/ULCY7H94MdOYn8M1OQ8myAxYLCn5/Ala0lYG6PT/nhe5r9uqJ8bg7A8TO6h2qQ+vp8nrBC5RJM3BZYIVgrnM0Gzm4P2o1PJJ78DICivxXFK93Aiiqc++8GTbH/TBRRGKFvyvzvLBVJS7XhWOK6pQH2cwOQzl2Y2UcWi1+g0+B9zch/Yi/IbI5cuosIJJh0iToUr6tHx41mlHwgpsxYNiI4HsLiYpz/owrUfO+iTxq7hc6CUBuOB9ZX94NZTRj4cjP47GwA8PmRkyDg3CMNyDszEb12H2NhDUDfg83ga6vDPuGF9w6h+vEeTQaALBbwtdWqyycFsgSx6zJqvncRl3ZWYuCB5rgnp0k1jJ6ACpx3NqJ/LY+Kv1JmWadu3gzLa8nzNGVb1kM4fxlSX5+u9XJZWfBsqgH/fvDJUoPUwegJRAHJwNIXh2ZeBzIA3PqV4Ktjk2wzHLT7mO4GAADksbGUNwBTN21OuKpUspOABV39ECqXQO7tj7k2oO3FfZC9265PN8B+YRjwiJDbLs2sRnC9Q2Du2CsGRcrUTZuR0TWWdFlxYo2tdQiSR0Tv17aAdzEUPX8S0uhoopuVVKS0EXBXOGByTgJxEAjlV9bAVZoN29vHIblccG/bDGvH5RkjECp7z+B9zch/cn9CPcqs7xyDLGk8PxE4i2XGv4IsFjCPmLSecYGQzirq0tO6gqnT8viR0sMB7ndHVHeDPTduisolWD7fBuuxSxi4S8n6an4jsFNSIApeOBH2hyNeXx/Tbqv72nXgV2qb5OPz8zB0R93M68lPrgdfFRtHKYPEkdJGYJq+rzaH/QFZLo8jozPy9VwmipD6+pD35Ox6Ml9cBK5uFc4/uTGky6oaWTFL9zig9UmtAdNbByGdPKvpGGlgEDnPzKouWV/djyvbFuvdtLjT+7B2yfFY4ryjEXxNVcLOb6wORABfvRSTyxwwD02BH3KCXe7B4O3r0jNFWYSMf6YROe+3QOqPjXdltIztaELuG6cTKj4bb4KtDqT0nEA08LXVILcHYttFzcfKF7tgvdIL2ekEVVVi9KY1sIwao8255Lx5GtJ48nrS5b5xWtGwNEi94cDo3U1+uoGRQBMuYDIybT3mmQ2gEVvbYX/xIGyXnBAql2D8s+mTz1CoKPfVWtCANDqa1BOI0vAIwBi6v7FlQYQ+hyLljED2z/aqdgUd3dkU1FtM7OjUrMQSNMhFlsCOnETH7WXgPIkfXgVDa5DOyOaSmUSv6UrZc62YXJtGadEjIOWMgBby3j4PeTL65BLT9N+7OaQLavlz7XAu4uKqjCtfXaco2qig9wFtqbLsL+zTJSozmRGvdMP0lr5RkalGWhsBqX9AV7kwx+N7ZhyThEXFfj8+sesyCh/ZExPvPT+I0PdgMzhRBpgcvjygOVKOL3Bg6N7UFlEZ+LL62AGhtASjd6fPcE4taW0EYsnI1ZXgMmYVjNhVdZCu3eibmJLIL+ORbjCGwh/sBfYej5kuotQ/gLynUjsxi+PxvTOGmywWjN0V/Ecudl2GySkvuBwHC3Z1IFrsL8y6EpPJjO4GG0qfOA0zz814pY3ubITJKaPzk4Tl3zikf9c6Hsu7SbCEHBVz2s/cbuT9rg2hvoXM3W2Qh4ZClEg/jJ6ADjCPG4v+bTekoSGfdfHsn+1Fxkv7UfM19WPruQo3WoQxQkK04J5ugRi8rwmyM7SLudTXByaKIEEAX+CIU8sUyfLp0PR4syCMgFBeBvG6ep99tGkN+NUhc6IkhKLHZiMUB7YuC1EyNJzNhonbFVUizmLB0A2R15Uu5D+xX3VSWC4nG2Mfj5+WAlWUwbMhMd9RWCNARD8mol4iOjFnX8D046Twn0TUQkTHiWhjLBuvFrmvH5aPfDPVcBe6gI4rQY7wZeI2fdw6adOasGXm9hhyn4rcA1F2TSF7f6d324XsZ/eCr146YxgWJBr8FqSBQdhe3BfDxsw739kWcL87ErfzzUVNT+AJANvm7QuWfvwmADXev10AHtGnmdEhu1x+M/bS0BCk0VElS3FxUchsR7Zf7QuovKsVT7YlrCw2X1OlT+4DWYLY2eWzS2ppi+uNbZAahDUCjLHfA5ifQGQ7gCe9208CuHXO/qeYwl4Aud58hElLzvMHIfX04vQ/V4FM5pieS3jvEE7/UzU4mw3OOwM/kaWWNmS8ntpCHgapRaRzAsHSj5cC6JhTrtO7z49QWYm1wq1bEbF763T3e9Xf9usi6BmOVf+3F/LkJHIOB/FWZCzhDjp8bTU8N8zOoXB1qxLYGoNYE/XEYKj042GO+yFjbBNjbJMJ0cXRO6uyo3ZvFVvbozpe9XnaL4HLzITY2h5wNpgEIeHCmGK+HdYjbTOvXf/shHh9fYgjFg6DXwqc4SreDDzQHPGDbz6RGoGe6W6+93+vd38XgPI55cq8+2KKK5fXRbc+HozubMTQ9tUAgIFbV8/cUBO3NWLi9kbwS8owec1qXc85dfNmTTeubOV9bjDr7YMw/f4jXduUqhT87AikjbUQygJ2cDUxeWvkcz+FTx/RzYcj0l9OsPTjLwP4gneVoAnAyJxhg+4IVZXg7HbkPbFH9668sLQiJok68g72IecZZXIu96k9M+pEWacGkHVyAKx/EPaTQaTKOD6iTDm2c/2ahhiW1j6wOWHA3V9YG5ehUiogu1wwXeqHrEMYcubpyHL1TrdDL9QsET4LYA+AWiLqJKL7AfwDgE8S0XkAN3hfA8BrAFoBtAB4DMBDurU0AO6SXFBGbLpmnpI8UAy65dK5C34WnK+txmRlLqSzLaD8XIiLg6wgyFJEmXKkljbVRoAvLMRQc6nPTRYuoWUqQiYzBu+LLC5C7OxS7W8Qimn9w0RjKAvFgambNsPy5uGg69QkCADPK4lKmtYp8QCJgghcRkbMFZyTAc5ujzqFGNWvBt87kpAkrGoRyssgFeXg7YN/beQdABRPOqHcP36cXx47by3bpdGQkX5MFAGZJd4AAABjC8IAANAlhyDfOwKmQ68glrBxJ/je4DJqC88IZGfBXeXvRz++0qHbbOt8pJNnQ07i8LXVcF8bOwMw8vkmzfMbMYl8TEPEjs6k1ymUhoZC9lQWnBEQu3sCumdmvLRf82wrt25FwF6FVuQL7bB8eMpvv3urX88tIizDsuZrY5sM3wCtcGtWKJJsKcaCMwLTSNdsBBrWRlUHNzoBpoNyERPFgF1wa48+3XJb+6j27EiJHpakINyYEyxMlGIykvRGQCgtweR2deupQ19UP9vL//YwsD+6tW+x/VJMJbXlo6f8u+UNa1UFIvnUc+KM6tWBWLtOpzPixY6Z+yGVPsekNwJSTy8y31eXP88yok5mK5UY2rnZJ7iJO3EBdCJ2S0t999Wr1iw0CAyZzOj/Yup4WCa9EWCiqDqBpL0z9bpioRCWVsA0KUMamHUqkScmwCQ5ZjLZBT/c47OU6byzEcLSipicK11hHjccj836c+gxbxRLktIIaOnWz4UdiKJ7n4RPP7HtIuwvzAn99baRy8mCc2N8JqDsL+wDy7AoqsYG2uF4DG0xjIBmnItjs1QXir5dDUlpCKYRyssw9hlFxFTqH4DlNwfCHBEaLT9q1noJpuPR6yksSGQJWc8lt1hrUhqBsr/X5qYqLCqOeo2/8Ad7Ep4xh8/LCxroI3Z0arqZuKwsDH4peI+K5Nklw74HQ0ekyS4XpHGnfpqHBklFUhoBLYzubMLgtUtBgilu5xSWVvhpFgKI2hC51y8FlYdIscbxGNvRFPJcnhvqIZSXQR4bQ/4T+/2On2nq7mMz24WPhveR4Ow2jH5saegLMEhJUt4I5P/+ErKf3Rv7KLemdTOb8pUe9K+z+OkB9AaL8GtYq8pA8L89HFTGjM/OBremBnkfKFqJozsaA+ZktB5th9TrlVKb07Nxb9uMvl8HEc4M0QPiV9eCz82BPDbmOz9hEBFCWSmEyiWJboYPKW8E1OYljJSBB5pBFgtky+xTVHa5sOg/9gBm395H0X8HHsbIFh1ccDkCM/Ez15v97F6IpQ7IH9vgU0zqH1ACkeZhfuMAHP+iPSqSmXggzfMRxhWBBxNCf57xlodPeSMQawp+cgBsasrP1Xjss43o+x/LVdXBfRC9AIQ0PAJ26KTvzgMn4CyxqJ7Q1Kpmy+fmoO22XFASKOmkC2L7JUgtbSHL9N8cP6lzwDACYQnmaZe3u9N/zK0S552NEKoqw5aj+lmFod6HtgBEkK+uA7uqzts4hry9XX7deW7dCt90aBEijztR+fII5IFB8Lk5ISca0x2131kwOr+zBcNf8P38+OKigLEGeU9GLjXPWa3g1qzQdEzS6gkM39MMx6tnIHlTQnFWKxhjAbu66YrnhnqY3jmkrnDDWowttUOykF++As5mQ/+O9cj/ceQ3VyD4mio4awtgfTUyY7jQESrKceH+MlT+3WHd7mvOboe4oQbch0f93nuHvZBaegJ5Pz8EuaoEfK3SNZLrlgNrawKWnby1IaLQ12RPxOFnADg+aJs5t4SMPg/ynvU3GvLEBBxPH4LzDn2vV2ppQ8YbqSePrvfnECnixQ5dDQCgaCQEMgChSFojwDxucJd6gR7vTPfe42AHTwQsm3W8F0zSvsafdbw3fKEYEGnXkjhCX13gr0w+egrCe4eUz239Sr/3mceN7GM6X2+M5dEHv9QMPjdH93p1/xyiYNoA9D68JWFtSNrhgB4M3dsMxwvHdVGQSSTyxzaAn3D7TwwGwfXpBtjePg7Z5QK/fBkmluXD8np0HoYGqU/KDQf0wPHckZgbAGHxotmJuhjB7z0JHDururz9wjDcVymTilJLO6zvGdoABsFJayMQSpaZLy7SJ/W0tyclXl+vOc5f9Sk8bnA1gb31+MJCDN2rzDqTyQy+pgrSqXMQ3j0E16cb4PpUPSZvXK+UrakKOQ/Crww852IQXy7/+RZd8hqoJdKsxN8loi4iOur9u3nOe9/2ZiU+S0RbY9VwLQhVlZj6lO+SGSvMB/JzQx5HghBSHMK9bTMm1ygRYsK7h4LOWcyFz8vD6N2K6y+XlYWRzzeFPQYAnNWBZcilvr6ZJSUym+BaOqs9YH1l/8wfAEjnW4MnJCWCc1noZKkG8WHJzy7i/ENLYpL3IhBh5wSI6OMAxqEkGl3j3fddAOOMsf83r+wqAM8CaABQAuAdAMsZYyFn7WIuOc7x4KwWzSq63PqV8ORngH8/8Az41Kc2I6PLCfnoKXBrVoBzTkJsuxi6UiJwNpsyTJm7HQNo0xrwVwZj7lVpoD+c3a7crzrO2UU8JxAkK3EwtgP4OWNsijHWBiUJiQ55ttXB11QFfkOWNBkAMpkhVC6BfOw0+PcPQygt8bPKA19uhu2Ds5CPKgKh3PAY2Nh4+MoZm/3Rz92OkL4Hm4Muj/L9o2Dz6heWVqS1knDQeyDFkJ1OgDFwdnvAGBE9iWZO4GtEdNw7XJjuR+qelXjwPuUmVyM0MrY2vM81WXyTn5LJ7Bfcw2VY4VxZPPPaXVUMbt5SleNHe3wUj8TOrpjqDQaj8JE9AZfpyGJRXFTnyWE7VxaBy8pKeNLTWKHmHoglkWY1CgaXm4PzX6sImLxWt3NEeNwjAJYBqANwBcC/aK1AbVbiwl8rIpkFL4XXGQw63p1D3xc3+rxu+249+GrfSTdpdNRnSY374EjKdakHPrcx4A/d8toBTG2oQst31/u+QZQWmYfV3AOxpPDX6vQw1SJ2XUb1v5zD1ObYTdpGZAQYYz2MMYkxJkPJOTjd5dc9K/G02/D0f6G0BGSxRKzvXvCor+vssp8NgnUFSQCapPQ92OwXNDTf2ST/x3uCDoGE9w6h6s/3QCgr9RkamAb1S3K5UJm+T3Wtc3AY5oFJcFZrTLQlIzIC02nJvdwGYHpa/GUAO4jIQkRLAdQA0NWxfGJtKbisTIxsLtFFDkw6eTbl0m4VPuKvglT0/d3gc3NmVh6CMXZX00zXcryuFFymd66DMbAj6pyRVNOw1icIyiBCZEmZfM7LxeRq/ZcOw84QebMSXwOggIg6AfwVgGuIqA4AA9AO4CsAwBg7SUTPAzgFQATwcLiVAa2Y3zgACYD9hfiPv5Menocrn0Oo0WPuO+cgedOOW1/dDwmA/IkNmsOM1cCduABAuUkMoke80g3TlW7w1UtBHhHixY7wB6kgrBFgjO0MsPvxEOX/DsDfRdOoSOALHJDHxsGmpsA78sGcE7rmcE8FpIFBFP1XaH1GedyJc/9dj9pHnTNPfm5Shf8/x0MoLoR4Rf3QKdV6WKkCeURA1O/Zmjoeg0QBl7b41bVgzevhqqsEX1gAAHCvqwS3KHpRTLWOPKmCUFEO9yfWYvlDh3y7/ioyMXF2G0abjfwDyYB4sWN2olqHIXHKLBjLV63Hhc9aUPN1X1FM1tYBQRBAo6OYfp7x7x+GHrFtjnfbdaknWZC6e2Edd0KKQFVZHhtL+My7wTw4Hn27GhSl7ChI6yjCZGDayUitU5CwqBhid08sm6QJPi8P8rgz9kKuBjFnQUYRJgNseSXYSvVS3f2f1N/jjSwWOO+MTEhjoqka/Jyh1ejO9BoipQvubZsjzguRFkaANq2Ja/Yg3pEPfvkyVWXZkZMBA4vmhh/zBY4Zd9fcp/cocx9zJM6jhbndyNkT2Uyy5fUDEDs6Z17nf9gZonQKwfExi/pMBBlHLkIeHI7o2JQ0ApzV6uMNJ9njl3gEACAIYNbw5+Rzc4IaJ/rD0dltQQCzzEYrMlEE9s5qAEjXbIz6hmXzVkr6vjrr3iqUlmDsLnVP+LkGIdWJ+30TQ6Se3oiHbCk5J0AbVoNZeJ8fSjIycXsjsv7QBqknSjkrjgeYHHFEGWe1YuiOOuQ8M5vGjExmn5uGBCGmUmFzGf9MIzJ/YUwyTiNdsxGWlh6InVE514YlreYE2JGTEK4MoeMvE6fLpgbbi/uiNwCA4h0YhbGWXS4fAwDA76mh1gDooYWXezC13LRjjeVMF+RBfd2NaYN6T82UNAIAIA8Oo/S3kwCAgfubQ4p/BKN/V+hEnLGAz8vz8z/o/4q6yLORzzfh8p/F3/D175ptX9H3tSWLDURYzYUFhtjdo7tj1VRRBga/1OwXNRuIlBwOzIcslohkmzmrVZNXoeuWBmS8fSxqiej57VXbDrJYQDzvc8NMbm9AxiuHYppRWevnZJAczL/P0mo4MJ/pC+Ud+Zq02QLd2HxhIYTFiyBUVfoJidjP9IF5gnebObs9pJQ4X+BQhEnnGZHuL28McoQvbGrK74mReapfmS+IIXoZALJYVK+qxLVuIvCra/VtUBKg9mGVukbAK83lsysrE3JBlDr1uVlg+TnwLMoBzatfamkL+cSlDCs8i3ODv5+dBdnh/344f/9QSOdbdZWgUsOMceR4cBryFHIWC9yl+ucRAAB+cTHOPRChqzhxmFgSO9GOZCDUkDM5hgP2EnZVwV1hl5/4vDy41y8F/9vD4AscGL6hBlk/3xvyGAP9Gbq3GY5fHAOVFGNieQEsrwXIacDxmNq20ee9qZs3w/LG4dgMXYjAZWSkddASX1gIcXmpz/KyWji7HW+NP5m8wwESJbCxsbDl2OQkzJ1ekZH+gZQ3AHx2NgbuT70kn3lPegVLeD6gAeCXLwOYDFu7r7SZrX0kdkMXxtLaAAAAm5iA6XJkqwih3NaTwggwtwfS8AhGPtcELisraDnZ5Qqb1jmVkEZH4Xhc3ySh8eTKDUWYvLXBL9vu6W8pkpPSqXNKJmMicFlZkM/Ff+iSTshOZ0xWVpLCCEyT9/JJyOMqFHsN9IcI7q1+PUUASggyt84/3XXxo/uR+bvzsIz6du9XfuPczI+98IWTkK+uw5l/Xgm5MX3cdOMNt2YFhMolsak7JrVGiDw2BjCmi0NKIhM8piSMwdoduMvInBPghuYYZ46HULkETBQhjYzC3jrqU36uCrM0OgpT7xhW/MmJiMaygRj6YjM6v7Owvl9uZBzMOanpmP6vNKuauE2OicEgfgLzXVsXAkJVJSaWF8L8RuQJRAe/1AzHUwfUeQFOxzZomKzjrFY4t65Dxku+8pFkMmNoRz1yn07+IY6ae2tyewOyDnamjNI0V7cKks0E2n0s4Psp5ydAJjPa//emiFWFUxWp4zKsf4hOtjr33CSmrq9TVVa+ah2ofpWm+mWXy88AAIorsuNVfSW3Y8XQznrweaHTrmW+f0Yft+94caYV/NHzmg9Lmp7AltpdGF/lQMavZ28uoaIcrfeVo+KlEeCjs3ELcDFQR+/DW3RxIwbHQyhZFPMAmnghlJYoUX1Jdr8mfU9AOnfBxwAAipZa1aOt4P51CNyyysQ0LJ2JUoOh6L8Dd/sDaUGGSn3GZVgxurksqrbowfA9+sSSjG0qiyjD0+jOJk3OV3qhJitxORG9T0SniOgkEX3duz+fiN4movPe/3ne/URE/+nNTHyciNT5xAZBvNIN7BDR+rnCmM2OLlSGP9cQXXr2IL3I3gd8M0CDCL1f3hywLKAsfdl+lfjQYsdbF3RZwsx4ab/P5Kha8t65ADnKuJRIUNMTEAH8L8bYKgBNAB72Zh/+FoB3GWM1AN71vgaAm6AkHakBsAtKyrKokHp6sezHXWi7pyzuUX+pSO9D6mbOc5/eEzZ/ovPORghLtakML/614ssxky2HsajFMONBIsb/V/5kCwYeUPwspL4+VUaIX75M18SrarISX2GMHfZujwE4DSXJ6HYAT3qLPQngVu/2dihpzBljbC+A3HkZiyJCbL+E8r/dA9enNkNYVBz+gAVM0Q/0e6raf7lfs4PKwHWVIEHA4DWVyg7DcAdl8b/vg+NH6j1fOasVu37zJk5/M1+3z1XTnAARVQLYAGAfgGLG2BXvW90Apn+ZqjITq81K7ANjyDx4ERceXGYMDUKhp29+BN3jnGf2gokisp9Vbu7RnY2Y3B63DPWphUbBGNnlwg9v2YYVf3xWN+9L1UaAiDIB/BLANxhjPgMepiwxaGqR2qzE8xG7e7DsJ50AlzRzmgZhyP7ZXlgGPYluRlLT92D4GBKhohzOOxvhWZQFeWWlz3t8Xp7P5KtQUR7SBX8uqn5JRGSCYgCeYYy96N3dM93N9/6fHlDpnpl4PmL7JYit7fDcuMkYGqQI3Ae+uQ6ddzaGXDFYaBQ9Ft45TLzYAfuLB8H9/ihwwFfB2nl1zUwGLgBou6ccYl21qnOrWR0gKLkHTzPG/nXOWy8DuNe7fS+Al+bs/4J3laAJwMicYYOuZBzvgDQ0HIuqUx4ymcHVaXMC0gPxunpVysjZx/pAa5ZHdS7OZksbMRAmiqANq8Mbxunhw7yhgPWV/T55Isv/dref4Q2Gmp7AVQDuAXAdER31/t0M4B8AfJKIzgO4wfsaAF4D0AqgBcBjAB5S1ZIIELt7wKam0L+rOeLEC2kLRxCz1A+zQla1fiXE6+tVlRXeOxQwz8JcaMNqTFQ7IGaG14XkbLagOnkkCBCz47+uHgv6vtoMKdsMkO9PkgQhaLeez87WJd9G0ngMRqMxyBcXoeWPl6Hmexchdl0Gu6oOpkv9aaWRn1A4HsSRfh5wGuIV5E9sgKl/AtLJs/qcO0kJppMpVC6Bc1VxQN0G5x2NyP59q7K0qIKk9xgMBb+yJqQnldTTi5rvXUTHjkoAgOnCFd0lnFMNLisLA19WL1gy8vmm2XX9+ciSJgPA2e2h17FlSfUKBve7IylpAMTr6kGb1/rt54uLIJSW+O0Ppgcotl8KrNwEwP7LfaoNQChSwgh4CjJB5tBdR7HrMhb/y26ACEPXLAWVRe2a4Adnt2Po3tRQApLHxuD4kXoHnZyf7vUZU0YDWS0QC9XNTE8zdfNmXR1gEo3w3iGwAwFSvmfZwbLt/vsTSEoYAe6DI6rcMCe3N4DLzMRwDQcMjoQtrxXZ6YTjOXWTLakE78iHdG1U3t0+SAODQcNZ58KuqoOwqBiczQbGE6QL6Z+PQGppg3Rae6RfIMTrwkdCqiEljIBaMs8MQh4bQ/nf7AY8bvC11Rh4oBlcVlbAkOTeh7do9rrSU3+fBCFmEtxaYM4JWNr6o6pDi4gL78iHsHgRTJf6IY+Ng7ndsJ/pB2QJnN0+40ZrEBprW78uuopJaQTm6/2rRTrbMrNNWVloubcQRS9fAJdph6dEedrNlckq+u894FRkaIkVZDbDVRm9JQ/E0L3qVGUAxbCJ7ZeiOp+mkOK8HMiFuRA7OiE7nYpC0flW8LXVcF29Eo7Hkj/OIBkQ2y5GnQgHSNLVgYEHmnW5EfjiIvTftAx5Tyh1kcUCSLOTXEJpCUaay5H9kRJEM9eIpDqc3R5SYVZvpm7eHHQCSy0kCCBBMLIdxYiUWh3Q60kg9fQi74k9Smag0hJwuTngMmd7GWLXZdhf2AcaGQONhJc8TyXiYQCmP1cAyLgU/efHRBGy26MM4dbM9tg4m01TZikDbSSlEQCAsR1Nukx6AACXm4Pzf1SBybolQEG+3/tidw/E7h5dzrWQ4PJy4alUnLTkE2eU76ywMCoHFjIJsF+RwM63zTgJcdlZmKrWxz1cvK4e/CoNnoocn/buzUk5HAAUbyjJqz6sB3xhIS58vRrVj3VCvNgR/gAd4bKyIK5fBu7Do3E9b7zhs7Phrq+GMDoFduhkVHX1f6UZwiTgePkUpGH9Vno4mw3MI86IjHpu3ATTWwcDF25aB9FmAskM/G8P69aGRJFSwwHAK1uto4GS+vpQ/VgnWu8rD+isEVM8HpgG0zs7DqB8Z/z7h4MagLnh37wjH4P3BV8FKHh0D3Kf2gNXQ43iU2+x4NJfbQG7qg6jO5uCHhcOeWLCR2XY3BdcxpsfmoD5wDnEywBodfDS7bxxP2MkcLyi/xYl4sUOVD3WjnN/VBHXWAPZ5YJ06lzczpcoyGTGyOe9P1COR9s/NPt4zY1sWjwzVJAGBpH/k/BzP6a3DoIdOQkigmuJG84SC7J/rp9oCjsSvMcinW1RcmHECdUOXk3rdA0OSw0jIEtwvK7Pj0jsuoyaf7uA1oeXLTg581jDPG7kv+F1hJElZFwhjNTMTsTaX9gXkeCJdM1GyC4XVn6zDZm/2JfUqcyka/RzugoGd7wFONOqX3261RRjwmnhaaqrpxfLHu/EhS+VG3oEOtH+N81w3tno8z2VvXgJ+b9XP/8ydwVg7jZJDHxxEaT+gaRfJeAnYi+eIk9M6LqMmjJGQG/Eix1Y9sgFODcYMmV6UP14FwSnb8ZhsaNTUy6BwY+VB9w2t/bAtV75ngY+UZ7cs/X7A8QL6BDuG0sWrBEAlKVByxsHwbasT3RTUp6ppQWw7b8QVR3TmoQAkHVp9kkndl2emcEfXcph5LObUuc743j07UpufcW0NQKc1Qo+N0d1+d6Ht6T80IAvLkqYsi/nliGPjqN/lz6z267CwFGj5X+7Gzm/OIyWnYpLNF/gAJnCi5MkDFlKern1tDUCVLYY7g0qgnMYA+0+hpJfXsBYkzZ9/WRj6LoqkGBKyLnpD0fBPG4U/Mg/R6G2ihQjNpUVvAvNPG7UfF3pGYxftQw9DwROqW6gjrQ1AlJLG/j356zvcnxI7TuxuwcZv94PobQkPqsGDf6CE9GS/ezeyLM4E+nTpijlzgfubwKXlYXcp/fM+AcEO4+wqBjuLA4lL8c/BJkvcKSN/kHaGoH5EEeQMlQ8JQUerffFfvJJsiXmiR2KRLSJLBYfDT3Hj/bMrM3LGQIoxPBG7O5RxFDmTD6y5vjMFZAgAObk+w4jIWndhhOJmtz1sYBdVQfT5SHNGX9SGX5lDTxFmeB+F51YC5+bg4mrlsPy+mG/3oh8dR3MnYMRhUuPf6ZR8U1IA1LObTiR9N5fDzKZwRcWxtWhSDjXBbk3OnGPVEM6fR7CsCvqnlfvZ1bBfqoXtH6F33um85ch94XwMyEK6oGXeyAmavm6QhtWRzUhHE1W4u8SUdc8GfLpY77tzUp8loi2Rty6BFH4gz1gHjfIZgXL1J5iOhhksYQWTO3ri6sGgFY4uz0mwyRPfoaf1LZWclrdEC92wu2Y/XzJZFYChsadYO7ATjwkCOBzsuF2ZAR8P1qxlWghQQgrsuPJj052Xc03Op2V+DARZQE4RERve9/7N8bY/5tb2JuxeAeA1QBKALxDRMsZYzomyIsP86MNXbc0IOPtYxGrudCKKohZlpSNJpz8+CrYz/ZBbG3XpT75ExswVGOF40d7tOWwC4Dw7iEAgOmdQzP7uOoKTJVkg/EEa/sQpHP+fgz8omKMNpSBZKbqxxBv+PJSOFcWBRVs4VfXgutxQo5iWB9NVuJgbAfwc8bYFGOsDUoSkuT2llBJ5qleMI927f2+B5sVxZxjp1PWAACA5fUDuhiAacEQ84VeFB4Yjro+QJmtn1aC5lfXAhwPxvMQbTysnaMBDQAAiJ1dsL9yBD31yenVd/lTpaEVm7r7QN3RyY5Hk5UYAL5GRMeJ6MdENK0Aoiorsd5wVmvM3TPF1nZwGaG7XiQIGLjf12Gm8JE9+iXuSAMmy5XVALGzC/Kx07rUKfUPIO9JxSlnsiwLxBFad+QBDHCVKOfjCwsDhi8zjxuVf5mcDj1F/xVau1EaGPSJ1xj8UrPmIVs0WYkfAbAMQB2AKwD+RcuJI0pNHoKJT66DsEQ/WyNfXecXbszZ7Wj9zjoI5WVBj2OiiMJnw8ttL2Qsr0enRRgO85sHwUQRy/7pJDLePjYzRBjYVo3C16JzbU52Cn5+TPMDJ+KsxIyxHsaYxBiToeQcnO7yq8pKrCY1ubCoGLzDXw4sENZX9us6iWPuHAQb952kk51OVP+gA2AsZDSbHjLQ6Qy/siYm9Y5/tgnC0jlen6XF4CrKZnqIuU/vgdTTG+ToyBBKS5ScgCrha9VlCo6USO69iLMST6cl93IbgOkslC8D2EFEFiJaCqAGQES+pHJxPpCfG8mhUSO2Xwo4Uz+d3/D8Q0tU539PWYhUy5ZrwVkdG5n1zOf3+vhYjC/Pw0R1PsgkBFfsIfLxShz/TCP44iJl2KBieCmVOEB5KmNUiOBcru6hphdksYS9jrDOQkR0NYAPAHwEYDpW9DsAdkIZCjAA7QC+Mp2CnIj+AsCXoKwsfIMx9nqocySbs5AauKwsYFk5uIHRtE18yjvyMbhtOXKe2Ru+cAJwb90E81uHVImM8NnZM1ms+OIiiNUlEO0CMk53Y3RzKWwvKtNcnN0OedIFzm4LqSoklJVCKswNqUyUDLhuaYD9TB+klragzkLJ6zFIBGFphW7LUbFgZm5AliF2XU5sY5IE8bp6mMbcgfPw6Qy3bgXk42c0H8dnZwPFBWAWM+QT2o8HAD4vD5STlXA/Ai2knMcg8TxG1yd3aK/Y0QlPmQNt91UmuikxQ2uYbtBEnNP1CYJu4c6RGAAAIEcexlcVRGwAAEAaGoqLAQglxhoKfnUt5E9sUFU2aY0AE0XYfpX8PtvcgZOo+I+PwK+sSTrNQm79yqg0EnhHPobv0lczb/zW+oRLhDGzCR67cus772iMSURnpMz/zgp/HZmhklvaIRxQl9I9aY1AqsBEEfLYGGjCBeaKfqlTT8jpAvNErnknDQwi56f6zgfYX9gXcg5l4MvNMZ9wlc62IPtnynVl9LvBjykqRslgxMnpgrzIMfNaGhqKqB42NaV6pcAwAkGY+tRmTU8s8WIHBrYug/PORtXLmrFGammDNDAYupBe8tVE0WUe8jq4zA0ljgejFdaZVOGn/7oo4fqFUksbRmrVLznqgWEEgmDb3QKpx9cdkyyWkNp2jtdbkP3+eUhD+mXMCYd8dV1UPz7uowvAufao2yFULsHErdoUfrj1K8E78sFnZ2Pgns1RnX/6cxDKSjWle3e8PpuEduX/NwDU+UchxhN+dS1yf6tNTlz+2IbYRhEuVKShIX9NAZmBcwePgyKbFSPXLwdkaSYUOdZwogwwOXzBIMhOpy7OTWLbxZllNrWQRwJkBml0VFUiklBwooyxz2yGWJIPEtXHqkl9s4aeCTxa7spKbG9AlABptv1ksYAvcIQ4AOA80cXmJe8SYapBpKxXe//3PrQF2R0irK9EqblnoJ7p7yDY2xYLRm/dgKzngs9zTG5vQObuNh/jkEiEinJM1hT5REeGJcjnkHJLhKkEmczof8Cbfsv74Zf8qtUwAPEmzAONud3I+11byDIZL+1PGgMAKHNNWgwAX+DA8D3acjUaRkAHmMeNgh/6dmfFK90z2/1f0R7ZBUBzvsR45ldMBfoenLfGzljiU9A3rFXG8EHgc3OCi6uqQOofQO5T2oZWhhGIA4WPH8DU9XUQFi/SdNzAVvUTXAAwdH1VzDT4nXc2xlXff2xHE8hkhvPOxojrKHosttGKEXHgRFBNCc5mw4U/WwVOx2hYNRhGQCeEyiVBHXOYKMJ2vBOdd2mTqNZq0bN/FoXkeBhy9nWBifrl2eMd+Rj6YnBvuLw9yvly9qlPYzafcCG14vX1oPrVIcuEkqmPCMaCDltk1xSqH7mEKzcqDwuhqjLk5PLI55p0SZhjTAzqBL+yBjTh8pMkm8+V/7UFnAco/s/QYhHpAu/ID++rEGXdiqAMp3qVQ0ubxOvqYTnSCkgS5ImJgIaFr6mCc0WB7nNA/Opa0PBYxHEpA/c3o+Cnh2fk8IyJwRgjnT4f1gAAQNlPW2AdiHxJL9Xov6U2ojXs8c+EHwb036LIiGFFFaT16jUK+j9dq7qs8N4hDNyyAuPXrQBfujhgGamlDbZ3T2DitsiHLgHrPXk2qsC0gqcOqdLDNHoCBkmJUFWZVBGkzjsbkXO4J3ibiCBULklczggi9D7UjKLvB+9hLuieAAkCOJs66XAtKjHR0L+rWRHtSHdhkgjpviHwUzcR8NnZyHr9BK5sXQw+OxvDXwgwl8FYQAMwurMpYqcxTfciYz4GQEti2AVhBPglZZi8NvQEEACACH13qiinA0VPHQGXmYmh/6Gcb+rmzTMqPmQyw3VLWgg0R0zRU9FlJNKT/ttXY3zrGix6rw9MkmAZVTz05n5nwch58Qikfu0JZTirFQO3qr8XpWs3+ngWavn8jOFAEsDn5sCzvgrcB8eVFFocD35ZBaTz2nzIg8Lx4GurZgJlUoYkbDcJArjKciXhbfVSSK2XQiZhHb6nGY63LuiubTgfoaIcck8fZJcSEcnVrYKr2Abzmwdnyizo4UCyQ1lZuHSjdVbOXJb0MwBQBFomK3K1HVO/WgnKSSDEkeZ2h0O8vl7JSxAhTBQVA1BbjcllDkCWQmYIioW4aSDEix0zBgAA5KOnfAxAKAwjkCQsecMVsxRkzOOG+Q2NjjPHzoLfdyom7VELE0Xt7Q7C1KeUKEXThycgn41edly+0A7LbxUFpb6710VdXyIxjIAXYVExhu9pVkJb4xD9ByhdS6GqEmxsDOZ2xV+dy8pKuPIOoPwAE5GZOVbY2hWRUTY1pUsiGM5mA5+vqCY7HkuexCV9X23W7NlpGAEvYncPcp/eAxTkg5UUxOWcZDbDuaIQ0vDIjNoOl5MNT0V8jFCyQiaz7tmkpJPqpLYCEWjyb2pTDYZ/YgNns0GoqsTkrQ0Y/JI6mfJYMp1MVwtq8g5YiWg/ER3zZiX+a+/+pUS0z5t9+DkiMnv3W7yvW7zvV0ZyMfPhly8DX71Uj6pCIp1t0S01VjjkiQm/PHNiZxfoD0fjcv5kZeq6deBXBk7SweflgTUrwi6eGzfpJloair7PbfB7upp+fwy5X3RCnpyEdKkT9rdPImNAArcmNolVYomansAUgOsYY+uh5BnYRkRNAP4RSlbiagBDAO73lr8fwJB3/795y0UNOSdB4wsrs0/vw1sS3YSEYH7zYNAnN3O7IQwqcyeWntimcR/b0QShohyOx/2frkwUlYhExhSdSacTvEsGuNmflE82JB2RrtkINOk3D6FpiZCIbAA+BPAggN8AWMQYE4moGcB3GWNbiehN7/YeIhIAdAMoZCFOtNCXCENBgrAgkpmSyZz0cxDB2jh/PwkChnZuhmlShv2F2ClmO+9oRPbvW1XrH0S1REhEPBEdBdAL4G0AFwAMM8am7865mYdnshJ73x8BEFofySAoQzs3+wmXStdujEs3OBhkMoeMidcKn52Nwc/V61ZfrOi9P3Ab5+9nogjHa2djagAAIPvdM7oEZ6kyAt7Eo3VQkos2AIhajVHvrMTpSu7Te/y+aGEswU9MJoMf168N0ugo8p5Inhn2YBT+wJv6fHuDj5jp9P65xCpy0uccwyMhHZXUoml1gDE2DOB9AM0Acr3dfcA38/BMVmLv+zkABjAPNVmJDQLA8WAHT6jKv6f3eadhohg+B1+CZ8ljScZL+yG3JihQSCX8quVo+/tmkMWC4XtCxxGoWR0oJKJc73YGgE8COA3FGNzpLXYvgJe82y97X8P7/nuh5gMMwsDxYFfVzWz3P6BPTIFWb8CB+xs0yV6N3L1Zlb9For0SI6X3gegk0oHYXjtrvYTqx7vB3G443grtHKWmJ7AYwPtEdBzAAQBvM8ZeBfBNAH9CRC1QxvyPe8s/DsDh3f8nAL4V4XUYTCN7bagsoeDRPSBBiNqhiWRtdtnx2B5VsenT5Px0r6oJK63t0Au+uEhTb4V35M8YQT4vD0U/Phz2mJ7/uQX8quVB9SWnr31u3Xohu1yQWi+BLyoM67ZsBBClIHxtNbhHnWA75MQLZ6YC3lh/19ICCO8dwujdTeCnGLLfPj2Trjwc7q2bwLtlcFMSpAwB1tb+sNoBQlkpzv3PJaj5x7PKHAERRnc0IvtZX8lz99ZNsJ3pUSVKowUuKwtjW1fNTFCmXmryNIBbvxJouRRdTADHA5tWAftnM/2SIIBWVOPKtfko/t7CkCmLBtctDcjocuLKJ3JQ+uMToKzMiBR7OJsNZLGozg8oXl8Pc5/TJ3uyUFoCsdQBHIjdvA6/uhbo6lYmDudgRBEmAGbifZxHIkU2+3ZbmShCPnEGxd/bjb4Hm3URm4w1nM0WMtouFvCrazFxeyMminiwIyex6N93Qxod9TEAJAjgc3N8jwsy1JInJjQlCJWsHJjZdyggdl2GbBEAUu6LcNmFIoGZeIBX7pnJ7Q3ga0IL3BpGIIawgyeiT64pS0ElqgFg0YsXMHJ1ZXTniAOsdink1drUlqNFOnUOmW+fgmU0uKYjX1wE59VzQouJMLAtsMuyViy/OaCs5MxjYpEZxCl+HoNba3SXcpePnppZosx45RCkltAJVxamESACbdBPQYhbs0L3iZ1gyFfXza4WAJB6emPulBIN058zO3LSZ0gTFxiD7JwIKVsudl2G9dX9PsfkPh3eZyEal+6c4wPo/NMGuD7dAFuPJ7aekrIUdtiRkkaAz8uLrgLi4MkLLQulBTEvA8THZ12c+/Bo8AAjorhpJAZDumajj5a/J9+aUO9GyNJMhKaehBL0DAe73IPFH07A+sp+nxRjXFZWyFwMgFd3UGcfjNQzAkTo3x6lw6IsQXhPQ4LHMHAfHNEls2+08FlZaPnOagilJQlrg+D0gJucTVIivHso4JNo8taGxGb/TQCeG+rBFzggj40FHOLJY2PIfyb0fTmybRX4vJyQZaaZ3N6gaqiRekaAsZRwMU0E0ugoar53EeA4CBXlCWnD+BIbphaFV1DmPOFnxsliQf9XfJ+MvQ9v0aVnwa/UfywejoyzPZDHxv32O+9ohOvTDZi6aXPYoUHm83tVuyRnfdSrKmvUglsi5LKyop+sS3KE8jKc/UYZ7Jc4lL7eg4ll+bC8rn9evv5dzX6JWNXi3roJlnePho2QJIsFkJnu42bp2o0w7T8bM0k3+WMbwE+4wQ6Fca8OAAkCSBB8NAP1wFgi9NJ3l8655ZIQsaMTy/+pFSYng9TSDut7x2NynqKfHov4WPObB1WFSItbVoNWKqsK3JoVPoE7kSIsXgTOLcfMAAAAv/ckcCwyNaOerzRg4pPrMHVT9K7JalhwRsDxo4UxlJB6epVrlSVN7r5aCDUPMvyFZvCra6NeB+ffPzzjbEN9g8DgcFT1AYA8OgZTp19Mmw9ksUSlZMU87qBGji8sxNC9wScAF7/XB+tvDsF2IfaRiEASG4HRu5uiXwUwmEFtBia96pzJqOzw/Q7JYol4TC/19ELqD/3jDQeZzJAnZxPHBrsGLsOKqQpFx4HPy8Po3U2q6h/4cvjMP1JfH/KeDP4wkk6fV2Tnz0WviqyGpDUCuS99BGl4ONHNSAs4mw2t31kPobxM13p7P78+5PvSybOQzrb47Bu/pS4haspTN20Gmcy49OebfPQLOx/y9buYRhoeUVY2AEjDw8h9SZ2PQ+Fz/s5Bc3FvC93Fp01rfFZ3hEXFQMNaVeeOlKQ1ArLTGf+Y+TRFnphA9aMdaP/8El3rjWRS0P7LfTFZt5/PxO2NPu6yGV1jYJKEyl8P+OgXlj1yDELPSKAqZmFM9fxBuEln6+Ux9D3YHHR5lO8fBZtzLmlxAXo3xzZf5YJbHQCU2dfBz20O2SVLZshiidk4X1tDCMTzPmPfmbZxPIa+0KBtOTdAfTNvpfh3lgwYqwNzYKKIgpcSm10nGvru3RiV1xi/arm/HwERxOu16fzx1UsxuW3jbBUmM/rv8b6WJRS8dCbIkUHqW1EN1ycDaxdKTWtQ+O4lTfUZqGNBGgEAfmGWqUTBD/dEpS1Hzkkwl39PwjSobV1aOt/q43fPPG6f1Ze5EXdCZfihiHT6fFB/BmHYBVnDpODEbY1ho+eihbPbMfCA70SgsLQCzjsbffbJV9fh/H81hjTc3LoVcVsSnE/KDgeG72lG7jP7dRFaXKg472xE9gdtUSfMDCqLPmcVYPzOBohWTpfvbOL2RmR9oF5qOykgCjrHpVZWXigvw9jGEmS8tD9s2UCk3XDA8fo5wwBESfYHbei+LTrnG85qRf8XAz/B5KvWA41rAcaQ+Yt9un1nWe+fi4uar64EMQB8bg4GPxe6B8Ca14OzWsEsJrjy9A9US9megIF+kMkMLj83Lim0AaUbTSYBEATIY+NgU1PgHflgzgndXWVTiYEHmlH4zDEfJyyhrBSe8gKcf8CE2u9Nhld5DkHa9QTmwmdnw3lHY/iCBgHh8nNx/o+XafPumzu+DeX8E2AcTBWlkFZUwLVxKa48VA/ObkfLn9aCW1SkodUpBJEqBynHj/b6eWEOfrwctOcYancdBTsam8ns9OgJcDx4R35qjRGTDL7AAXl4RNXYlLPb0bdjHRyPK5OA7q2bYO12Bkzk2vfVZhQ+ujdod1hYvAhiTx/4QkfceiLxRr66DiQz0O7IYy3E6+thGnTFpCeQHkbAQDf6vtqMwh+m3oRr/1eaUfj4gRkjNnlrAzJP9kM63xq7kxKBLyhImYdPxMOBEKnJnyCiNiI66v2r8+4nIvpPb2ry40S0MeQJDJKKRBsA5x2NEcU5FDy236cXk/HSgdgaAABcRgaGr5szsRpLBaUo6h7dGTruIZrU5ADwZ4yxOu/fUe++mwDUeP92AXgkgnYnJWQy66pNGC9IEHwkv0LiNQDC4kUJESbJ2dsBOYAPQ1jmGy6de7jOOxohVFXOvOYLHBjftha5xwfAFziUVZIH1AUZaYUsFvTvCl03v7oWXFZg9+L8D0O7aYc1AkxhWg7F5P0L9QlvB/CU97i9UHIWLg53nmiZuM33S4oFvffXw+0IrE3IO/IhVJRjbEcT+neFjyRTAwmCPpqBxEGyaVTRMZvALPoq7/TvCu4zP43YdVm3nghZLEF/GD7lBAH9u5oVFaMgT1z7L/dBbG33OcaVw0G2mkEmE2SXK2KBlXCwqSkUPBq6bjnDBOID/5zDxWqomhMgIh7AIQDVAL7PGPsmET0BJTHpFIB3AXyLMTZFRK8C+AfG2IfeY98F8E3G2MF5de6C0lOAFbb6q+nmsO0I2UZBAJPCK6tGfY5Afu0mM4bv2oicZ/YpgqPE6aKEI5SXwbluMSy/0V8VKBGQyRxbZd158LXVcJdkg38/fMqwaakxPdpHgoDx7fUAIagStHh9PawnO3XLIOW5cROsxy6FnFyNaolwfmpyIloD4NtQUpRvBpAPJTehavTOSsxEMeZRh8FmzpnHjZyfKjPgTBR1u9HFjs60MQCAPj8wLUhnW1QZAMArAhKmfaN3NykahxwfUnKcSRJyDl1BzqHuoGWsp7og6+i6nnHysk99WoatkaYm38YYu+Lt8k8B+AmA6XS5M6nJvcxNW26QQvgMRZJAzjyRTNzWCMlMyD8zBeL50JLjjEFsv6TkKuR4n/iC6WxH4pVuVY5R87MjBUPsuuwTWTpVlKHqOCDy1ORnpsf5REQAbgUwrabwMoAveFcJmgCMMMauqG6RQcKRrt0IvrAQ/betnhnD81lZGLplVYJbljjsrx5B3pN7Ibx7SFuPRpZQ8JS3N0KEvtuUz1C6ZqOSGTkMfbdH9pmb3zwYvpCXaFKTP0NEHwH4CEABgL/1ln8NQCuAFgCPAXhIfdOTB6GqMu6585IFS2sf2NgY8p7cMzMEkkZHkf2zvWGOjA/D9/jmX5Q/sQFC5RJdREiDwTxucOtXwr1VGVL3fk19BqKZJzRjyP+JMsFnudALFkB+fD75P9ZnspFbFzxXh+EsFAS2ZT2E85eTwhFEvL4ewpg7/mm8UoT+Xc0YumoKbELA8gcji7BLd1y3NODDV/48eT0GiagPgBNAf6LbEkcKYFxvupNs11zBGPNLuZwURgAAiOhgICuVrhjXm/6kyjWnRRShgYFB5BhGwMBggZNMRuCHiW5AnDGuN/1JiWtOmjkBAwODxJBMPQEDA4MEkHAjQETbiOisV3/gW4luj14Q0Y+JqJeITszZl09EbxPRee//PO/+lNdgIKJyInqfiE55dSe+7t2fltccQmdjKRHt817Xc0Rk9u63eF+3eN+vTOgFzIUxlrA/ADyACwCqAJgBHAOwKpFt0vHaPg5gI4ATc/b9E5RoSwD4FoB/9G7fDOB1AASgCcC+RLc/gutdDGCjdzsLwDkAq9L1mr3tzvRumwDs817H8wB2ePf/AMCD3u2HAPzAu70DwHOJvoaZa0nwB9kM4M05r78N4NuJ/lB0vL7KeUbgLIDF3u3FAM56tx8FsDNQuVT9A/ASlDiTtL9mADYAhwE0QnEOErz7Z+5vAG8CaPZuC95ylOi2M8YSPhwoBdAx53Wnd1+6Usxmg6m6AUw7wKfV5+Dt6m6A8nRM22smIp6IjgLoBfA2lF7tMGNsOuZ87jXNXK/3/REAGuSdY0eijcCChSmPhLRbmiGiTAC/BPANxtjo3PfS7ZrZPJ0NKPoaKUeijcBC0x7omROCvRjKEwRIk8+BiExQDMAzjLEXvbvT+poBH52NZihyetMaanOvaeZ6ve/nAFCfXDGGJNoIHABQ451RNUOZMHk5wW2KJS8DuNe7fS+UcfP0/pTWYPDqSjwO4DRj7F/nvJWW1xxEZ+M0FGNwp7fY/Oud/hzuBPCet2eUeBI9KQFllvgclPHUXyS6PTpe17MArgDwQBkb3g9lDPgugPMA3gGQ7y1LAL7v/Qw+ArAp0e2P4HqvhtLVPw7gqPfv5nS9ZgDrABzxXu8JAP/Hu78KwH4oehq/AGDx7rd6X7d4369K9DVM/xkegwYGC5xEDwcMDAwSjGEEDAwWOIYRMDBY4BhGwMBggWMYAQODBY5hBAwMFjiGETAwWOAYRsDAYIHz/wMMWfwpUgPzwgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# the PCalgWrapper contains a function to open the PC graph and transfrom it into a connectivity matrix\n",
    "# tetrad2matrix encodes directed edges X->Y as 2 (M(y,x) = 2, direction column to row) \n",
    "# and undirected X--Y as 1 (for M(x,y) = M(y,x) = 1)\n",
    "\n",
    "from PCalgWrapper import tetrad2matrix\n",
    "\n",
    "output_graph_file = f'{out_dir}/{out_name}.txt'\n",
    "graph_matrix = tetrad2matrix(output_graph_file)\n",
    "print(graph_matrix.shape)\n",
    "plt.imshow(graph_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}