Report Lab Rotation 2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Report"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scientific background "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The topic of this lab rotation is very close to that of the lab rotation 1, where I was testing a Matlab implementation of conditional spectral Granger causality [1]. Granger-causality [2] is a common metric in Neuroscience that is sometimes interpreted as directionality (i.e. feedforward vs. feedback, [3], also see [4]) in propagating electrophysiological signals. If the activation is periodic, it can be calculated in the spectral domain [5]. The spectral Granger-causality is interpreted as the influence that rhythms around one frequency in one brain area have on rhythms around some other frequency in another brain area. A famous example is the work by Fries et al. around differential directionality of the influence by $\\alpha$- and $\\gamma$-frequency oscillations in the visual cortical areas of mammals [6]. Based on these metric, more complicated models of functional connectiviy, such as dynamic causal models can be built [7].\n",
    "\n",
    "With more and more experimental datasets obtained from MEG and fMRI studies, the insights gained from causal analysis may also increase. Therefore it is crucial to have easy tools for assessment and interpretation of directed influences available whenever performing experimental studies that include mutliple brain regions that are functionally connected with oscillations. In a previous Lab Rotation I had a deeper look at the implementation for computing Granger causality as provided by the Matlab-package Fieldtrip [8]. While Matlab is still popular at some universities, Python has a considerable bigger user base, has better interfacing with other programming languages and is open source. Therefore, it would be of value to have a Python implementation of conditional Granger causality, which I was not able to find. Therefore, I suggested to build a package myself that could be used in Neuroscience studies.\n",
    "\n",
    "The basic problem description is the same as in the Matlab lab rotation: \n",
    "Given data that we can assume to represent three distinct recording channels ($c_1, c_2, c_3$), we want to be able to differentiate which of the two following influences is greater: \n",
    "\n",
    "1. $c_1 \\rightarrow c_2 \\rightarrow c_3$\n",
    "\n",
    "2. $c_1 \\rightarrow c_2, c_1 \\rightarrow c_3$\n",
    "\n",
    "The first case is sometimes called “sequential drive”, while the second is called “differentially delayed drive”. \n",
    "By influence we again mean Granger causality, denoted $F_{X\\rightarrow Y}$. It is opposed to conditional Granger causality, denoted $F_{X\\rightarrow Y \\vert Z}$.\n",
    "\n",
    "The goal here was to construct a pair of conditions, which shows the shortcomings of regular Granger-causality on the first hand and the advantage behind using conditional causality on the second hand. In the second part of the report, we will see that regular Granger-causality can not be used to distinguish cases 1 and 2, whereas conditional causality does. This main goal was reached. There are a few minor things on the side, i.e. code quality and full-grained testing, that the time did not allow for. However, I think that for \"just\" a lab rotation it should be sufficient."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## My tasks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As laid out in the proposal, the tasks were the following:\n",
    "\n",
    "1. **Formalise implementation:** The base structure of the pycondgranger-package re-uses an \"Analyzer\"-class that was introduced in another Python package popular in Neuroscience: the nitime package. I chose this package, because it is still actively maintained. At a later stage I would like to add my code to this package, such that more people will find and use it. The base nitime-Analyzer only supports regular Granger causality in the frequency domain, but it also features multi-tapered spectral analysis. This is an advanced method for computing the spectra between which we want to assess the influence. I implemented the logic of conditional Granger causality from scratch, which took a lot of time. Copying the Matlab-formulas would have been an option, but that way we would no longer be able to cross-validate both implementations with each other.\n",
    "2. **Test functionality:** Where a mathematical formula can be disproven, using a software implementation for that formula can lead to incorrect results for a variety of reasons. Examples range from typing errors over computational constraints to a bad structure in your code. Therefore a test functionality is absolutely necessary. I provided this in the 'test' subdirectory of the pycondgranger code repository. The basic approach is to call all functions to check if it runs successfully.\n",
    "3. **Provide test script:** Some more complicated software errors do let the program at questions complete without formal error, but the computed results do not make sense. To avoid this kind of error, a sanity check of the results can be useful. To avoid this kind of non-sense errors, I concluded with the test script attached below this part of the report. It will show that the package is functioning as intended."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time investment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Formalise implementation\n",
    "    - implement the 3-VAR model - _6 hrs_\n",
    "    - adapt the processing pipeline of the base class - _2 hrs_\n",
    "    - adapt the normalisation matrix $P$ - _4 hrs_\n",
    "    - debug variables handover between functions  - _3 hrs_\n",
    "2. Test functionality\n",
    "    - find adequate test case in literature - _3 hrs_\n",
    "    - implement the test case in VAR-system - _2 hrs_\n",
    "    - create full output - _2 hrs_\n",
    "3. Provide test script\n",
    "    - re-write test script for interpretability - _2 hrs_\n",
    "    - create visualisations - _1 hrs_\n",
    "    - polish up for readability - _1 hrs_\n",
    "4. Write report - _3 hrs_\n",
    "\n",
    "Sum: _29 hrs_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reflecting thoughts\n",
    "\n",
    "As a beginning, I would like to thank Paul Tiesinga for working remotely on this project. Working remotely on a project is not always that useful, but in projects that mostly do coding it is (hopefully) okay. Asking Paul for the lab rotation made sense, since he is my second supervisor for the thesis and actually taught me about Granger causality in his lectures. Since more and more people start using Python in his group, I was hoping to provide some additional tool that might be useful in the future.\n",
    "The concept of a functional hierarchy of cortical areas knwon from Fries lab [6] has already been replicated in other parts of the brain than low-level sensory areas [9]. These hierarchies often share a moderator that sets the fraction of \"bottom-up vs. top-down processing\" inside any area of the hierarchy. For the visual hierarchy these moderators are thalamic nuclei. With conditional Granger causality, we can now formulate hypotheses like the following: \n",
    "- around $\\alpha / \\beta$-frequencies hierarchy-area 1 receives feedback from hierarchy-area 1 given the activation of the moderator\n",
    "- without moderation hierarchy-area 1 sends feedforward influence to hierarchy-area 2 around $\\gamma$-frequencies\n",
    "\n",
    "Unfortunately the moderator area in our case are sub-cortical can not easily be recorded simultaneously with the hierarchy-areas. However, with a targeted setup and combined methods like MEG/EG or optogenetic manipulation, this might become possible in the nearer future."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Literature"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- [1] Chen, Y., Bressler, S. L., & Ding, M. (2006). Frequency decomposition of conditional Granger causality and application to multivariate neural field potential data. *Journal of neuroscience methods, 150(2), 228-237.* <a id='ref1'></a>\n",
    "- [2] C. W. J. Granger (1969). Investigating causal relations by econometric models and cross-spectral methods, *Econometrica, 37 , 424-438.* <a id='ref2'></a>\n",
    "- [3] Blinowska, K. J., Kuś, R., & Kamiński, M. (2004). Granger causality and information flow in multivariate processes. *Physical Review E, 70(5), 050902.* <a id='ref3'></a>\n",
    "- [4] Measures of Coupling between Neural Populations Based on Granger Causality Principle, Kaminski M Brzezicka A Kaminski J Blinowska K *Frontiers in Computational Neuroscience, 2016 vol: 10 pp: 114* <a id='ref4'></a>\n",
    "- [5] Chicharro, D. (2011). On the spectral formulation of Granger causality. *Biological cybernetics, 105(5-6), 331-347.* <a id='ref5'></a>\n",
    "- [6] Bastos, A. M., Vezoli, J., Bosman, C. A., Schoffelen, J. M., Oostenveld, R., Dowdall, J. R., ... & Fries, P. (2015). Visual areas exert feedforward and feedback influences through distinct frequency channels. *Neuron, 85(2), 390-401.* <a id='ref6'></a>\n",
    "- [7] Bastos, A. M., Litvak, V., Moran, R., Bosman, C. A., Fries, P., & Friston, K. J. (2015). A DCM study of spectral asymmetries in feedforward and feedback connections between visual areas V1 and V4 in the monkey. _Neuroimage, 108, 460-475._ <a id='ref7'></a>\n",
    "- [8] Oostenveld, R., Fries, P., Maris, E., Schoffelen, JM (2011). FieldTrip: Open Source Software for Advanced Analysis of MEG, EEG, and Invasive Electrophysiological Data. *Computational Intelligence and Neuroscience, Volume 2011 (2011)* <a id='ref8'></a>\n",
    "- [9] van Pelt, S., Heil, L., Kwisthout, J., Ondobaka, S., van Rooij, I., & Bekkering, H. (2016). Beta-and gamma-band activity reflect predictive coding in the processing of causal events. *Social cognitive and affective neuroscience, 11(6), 973-980.* <a id='ref9'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "# Example analysis with the pycondgranger package"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup of the package"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import warnings\n",
    "import sys\n",
    "from itertools import permutations\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set_context('notebook')\n",
    "\n",
    "sys.path.append(\"/media/dan/usbdata1/labrots/labrot_pycondgranger/pycondgranger\")\n",
    "import granger_class as gc\n",
    "import granger_utils as utils\n",
    "from timeseries import TimeSeries\n",
    "\n",
    "def float_formatter(x):\n",
    "    return \"%.3f\" % x\n",
    "\n",
    "\n",
    "np.set_printoptions(\n",
    "    precision=3,\n",
    "    formatter={'float_kind': float_formatter}\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulating two cases of three (in-)dependent channels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sample process a: differentially delayed driving\n",
    "\n",
    "mu = 0.5\n",
    "a1 = np.array([\n",
    "    [0.0, 0.0, 0.0],\n",
    "    [1.0, 0.0, 0.0],\n",
    "    [0.0, 0.0, mu]]\n",
    ")\n",
    "a2 = np.array([\n",
    "    [0.0, 0.0, 0.0],\n",
    "    [0.0, 0.0, 0.0],\n",
    "    [1.0, 0.0, 0.0]]\n",
    ")\n",
    "a = np.array([a1, a2])\n",
    "\n",
    "# Sample process b: sequential driving\n",
    "\n",
    "b = np.array([\n",
    "    [0.0, 0.0, 0.0],\n",
    "    [1.0, 0.0, 0.0],\n",
    "    [0.0, 1.0, mu]]\n",
    ")\n",
    "b = np.array([b])\n",
    "\n",
    "cov = np.array(\n",
    "    [[1.00, 0.00, 0.00],\n",
    "     [0.00, 0.04, 0.00],\n",
    "     [0.00, 0.00, 0.09]]\n",
    ")\n",
    "\n",
    "N = 5000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Process $a$ resembles the case $c_1 \\rightarrow c_2, c_1 \\rightarrow c_3$, where the signal detected at channel 1 drives channel 2 with at lag 1 and channel 3 at lag 2.\n",
    "For process $b$, we have $c_1 \\rightarrow c_2 \\rightarrow c_3$, where both influences happen at lag 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "# inputs a, b should have the shape (order, n_channels, n_channels)S\n",
    "signal_1_, noise_1 = utils.generate_mar(a, cov, N)\n",
    "signal_2_, noise_2 = utils.generate_mar(b, cov, N)\n",
    "\n",
    "signal_1 = TimeSeries(\n",
    "    np.nan_to_num(signal_1_), # + noise_1,\n",
    "    sampling_rate=1000,\n",
    "    time_unit='s'\n",
    ")\n",
    "signal_2 = TimeSeries(\n",
    "    np.nan_to_num(signal_2_), # + noise_2,\n",
    "    sampling_rate=1000,\n",
    "    time_unit='s'\n",
    ")\n",
    "# including the noise will make you see the jitter on the spectral plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Granger causality metrics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing regular Granger causality"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "with warnings.catch_warnings():\n",
    "    warnings.simplefilter('ignore')\n",
    "    \n",
    "    G1 = gc.GrangerAnalyzer(signal_1, max_order=5, n_freqs=2048)\n",
    "    granger1 = G1.causality_xy\n",
    "\n",
    "    G2 = gc.GrangerAnalyzer(signal_2, max_order=5, n_freqs=2048)\n",
    "    granger2 = G2.causality_xy\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing conditional granger causality"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "with warnings.catch_warnings():\n",
    "    warnings.simplefilter('ignore')\n",
    "\n",
    "    cond_granger_full_1 = G1.conditional_causality_xyz\n",
    "    cond_granger_1 = np.nanmean(cond_granger_full_1, axis=-2)\n",
    "\n",
    "    cond_granger_full_2 = G2.conditional_causality_xyz\n",
    "    cond_granger_2 = np.nanmean(cond_granger_full_2, axis=-2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the case of only three channels we can neglect the third index, i.e. $Z$ in $F_{X\\rightarrow Y \\vert Z}$, because there is no freedom for choosing it. \n",
    "Now, the conditional Granger-causality of $i$ on $j$ via $k$ can be read from entry $(j, i)$ in the matrix.\n",
    "Let' have a look at the triplets of channels that we calculated conditional influence for:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Validate results with original paper\n",
    "\n",
    "We are aiming here to reproduce the Figures 1b, 2b, 3 from [_Y. Chen et al. / Journal of Neuroscience Methods 150 (2006) 228–237_](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.64.1269&rep=rep1&type=pdf)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAJRCAYAAAANj2DmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdfZBldX3v+/eHBh8OYOaIA1emgdE44zUiISOEWMngyYmoGePRpCSEY9SblBrIyVCTlH94TcUh5JCyIlVnYhjPcIgJCmaiRIIPmTIn5FyTwUBFEwZ0Ehye5wHJTKFD0CRoer73j16jm+3u6d3dq3t1736/qlbtvX/7t37ruxR+fHrt9ZCqQpIkSQvruK4LkCRJWo4MYZIkSR0whEmSJHXAECZJktQBQ5gkSVIHDGGSJEkdMIRJkiR1wBCmOUny3CR/muSbSR5J8l8XYJu/kuSLSZ5KcsN8b0/SaOpiLuliztTidXzXBWjJ2wp8CzgNOBf4syR3V9Xuedzmo8B/B14DPHsetyNptHUxl3QxZ2qRinfM12wlORH4OnB2Ve1p2m4EDlTVu2cxXmoG/0Am+e/AeFX9PzPdliQd1cZcMsz81facqaXPnyM1F2uBiaOTSeNu4KUzHSjJemBHkme1VZwkLYQZzF+tzZkaDYYwzcVJwBN9bU8AJ89irM8DB4FPGcQkLTHDzl9tzpkaAYYwzcU3gOf0tT0HeHKqFZK8Nkn1L8AE8FbgIuDyeatYkmaphflrxnOmRpsn5msu9gDHJ1lTVfc1bT8ITHmCaVV9Fkh/e5LjgD8Eng9sm4daJWlOWpi/ZjxnarR5JEyzVlXfBG4BrkpyYpIfBd4A3DiL4X6UyauF3lBV/3qsjkmObw75jwFjSZ6VxD8oJM1Ii3PJUPNXy3OmRoBXR2pOkjwX+AMmD8M/Dry7qv5olmMNdXVkkiuBzX3Nv1lVV85mu5KWpzbnkhnMX63NmVr6DGGSJEkd8OdISZKkDsxLCEuyublq5Oz5GF+ShpXkmiQPHWtOSjKWZGuSB5Lcn+TtC12npOWn9RCWZB3wI8DetseWpFm4FbgQeOQYfd4MvAhYA7wCuDLJ6nmvTNKy1moIS/JMJp+L9cuAJ5tJ6lxV3V5V+6bpdglwfVUdqapDTAa3i+e/OknLWduX9V8F3FRVDyXfcysVAJKsAFb0NT8DeCFwH5M3vZM0+saYvK/SF6rqqY5rOZOnHynbC5zR38n5S1KjlfmrtRCW5BXA+cB0DyHdxPdeEixp+VoP3N51EUNy/pLUa07zV5tHwl4J/N/A0aNg48CfJ/mFqvrfPf22ADf0rXsW8LmdO3cyPj7eYkmSFqv9+/ezfv16gK92XQuTR77OAr7QfO4/MnaU85ek1uav1kJYVb0PeN/Rz0keBn6qqr7c1+8wcLi37ehPl+Pj46xevbqtkiQtDYvhJ7ybgXckuQU4BXgjkyfzP43zl6Q+c5q/vE+YpJGW5ANJ9jN5dP62JLub9h1Jzmu63Qg8yOR5XXcCV1XVg50ULGnZmLfn7VXV6vkaW5KGVVVXAFcMaN/Q834CuHwh65Ikj4RJkiR1wBAmSZLUAUOYJElSBwxhkiRJHTCESZIkdcAQJkmS1AFDmCRJUgcMYZIkSR0whEmSJHXAECZJktQBQ5gkSVIHDGGSJEkdMIRJkiR1wBAmSZLUAUOYJElSBwxhkiRJHTi+6wJ67d4NX//6wmwrmZ++baw7aJ1hxzlWv9mMO+w6s2mbaa297XPZ3mzqmss4/f2m+jxMDccaa7pxpqtnuvEH9TnW/xaSpGNbVCHsox+Fk0/uugpJczVVUOtt/+d/XtiaJGmxWVQhTNJoqBr8vtfExMLUIkmL1aIKYT/wA3DKKfO/nan+o7CQ40y17rBjHqvfoO+mG3fY72czThvrDvP9dNub6vMw4w7bd5h+x9qnqb6bbp+mG3O67cxk3GECliRpeosqhP38z8Pq1V1XIWk2jhUgB3338MOwffuClCZJi9KiCmHcdRc89ljXVSxdy/GwxHLc57mYx/+9MsX7qZzw1a/OVymStCQsrhB2882emS8tF08+2XUFktQp7xMmaaQlWZvkjiR7mtc1A/qcmuTPktyT5N4kH0yyuP5IlTRyFtck80M/BM97XtdVaKlZSjepWkq1tq1/3w8dWqiTwrYBW6vqpiQ/D1wH/Oe+Pu8B/rGqXpfkBOB24GeAjy9EgZKWp8UVwi6+2DPzpeXi4YfhXe+a100kORVYB1zUNG0Hrk2ysqoO9XQt4OQkxwHPBJ4BHJjX4iQte4srhElSu84ADlTVBEBVTSR5tGnvDWG/BXwC+CpwInBtVX2+f7AkK4AVfc3j81G4pNHnOWGSBBcD9wDPB1YBFyZ504B+m4CH+padC1WkpNFiCJM0yvYBq5KMATSvpzftvTYCH62qI1X1BPBJ4McHjLcFeEHfsn6eapc04gxhkkZWVR0EdgGXNk2XAnf1nQ8Gk0e0XguQ5BnAq4AvDxjvcFU93LsA++erfkmjzRAmadRdBmxMsofJI16XASTZkeS8ps8mYH2SLzEZ2vYA13dRrKTlo9UT85PcyuTh+SPAN4CNVbWrzW1I0kxU1b3ABQPaN/S8f4DvXkEpSQui7asj39acT0GSNwB/wOTl4ZIkSerR6s+RRwNY4/uYPCImSZKkPq3fJyzJ7wOvZvIZvq8d8L332ZEkScte6yGsqt4OkOQtwPuBDX1dNgGb296uJEnSUjJvV0dW1Y3Ajyc5pe8r77MjSZKWvdaOhCU5CfiPVbWv+fx64GvN8h1VdRg43LduW2VIkiQtCW3+HHkicHOSE4EJJsPX66uqWtyGJEnSSGgthFXVPwE/0tZ4kiRJo8w75kuSJHXAECZJktQBQ5gkSVIHDGGSJEkdMIRJkiR1wBAmSZLUAUOYJElSBwxhkiRJHTCESZIkdcAQJkmS1AFDmCRJUgcMYZIkSR0whEmSJHXAECZJktQBQ5gkSVIHDGGSRlqStUnuSLKneV0zRb+fTfKlJF9uXk9b6FolLS+GMEmjbhuwtarWAluB6/o7JDkPuBK4qKrOBn4MeGIhi5S0/BjCJI2sJKcC64DtTdN2YF2SlX1dfxW4pqoeA6iqJ6rq3xauUknL0fFdFyBJ8+gM4EBVTQBU1USSR5v2Qz39fgB4KMlfAycBtwBXV1X1DpZkBbCibxvj81W8pNFmCJOkybnwHOAi4BnAZ4G9wEf6+m0CNi9saZJGlT9HShpl+4BVScYAmtfTm/ZejwB/UlVPVdWTwCeBHx4w3hbgBX3L+nmqXdKIM4RJGllVdRDYBVzaNF0K3FVVh/q6/hHw6kw6AfgJ4O4B4x2uqod7F2D//O2BpFFmCJM06i4DNibZA2xsPpNkR3NVJMAfAweBf2AytO0GPtRBrZKWEc8JkzTSqupe4IIB7Rt63h8Bfq1ZJGlBeCRMkiSpA4YwSZKkDhjCJEmSOmAIkyRJ6oAhTJIkqQOGMEmSpA4YwiRJkjrQaghLckpzA8SvJLknyS1JVra5DUmSpFHQ9pGwAn6nql5cVecADwDva3kbkiRJS16rIayqvlZVn+tpuhM4q81tSJIkjYJ5e2xRkuOAy4FP9bWvAFb0dR+frzokSZIWo/l8duTvAd8Aru1r3wRsnsftSpIkLXrzEsKSXAOsAV7fPBi31xbghr62cWDnfNQiSZK0GLUewpJcDbwceF1VPdX/fVUdBg73rdN2GZIkSYtaqyEsyUuB9wB7gL9pwtVDVfXTbW5HkiRpqWs1hFXVbsDDWpIkSdPwjvmSJEkdMIRJkiR1wBAmSZLUAUOYJElSBwxhkiRJHTCESZIkdcAQJkmS1AFDmKSRlmRtkjuS7Gle1xyj74uT/Evz6DVJmleGMEmjbhuwtarWAluB6wZ1SjLWfHfrAtYmaRkzhEkaWUlOBdYB25um7cC6JCsHdH838BkmH7smSfOu9Qd4S9IicgZwoKomAKpqIsmjTfuho52SnAO8Bvhx4DemGizJCmBFX/N420VLWh4MYZKWtSQnANcDv9CEtGN13wRsXpDCJI08Q5ikUbYPWJVkrAlYY8DpTftRzwe+H9jRBLAVQJI8p6re2TfeFuCGvrZxYOd8FC9ptBnCJI2sqjqYZBdwKXBT83pXVR3q6bMXeN7Rz0muBE6qqncNGO8wcLi3bZojZ5I0JU/MlzTqLgM2JtkDbGw+k2RHkvM6rUzSsuaRMEkjraruBS4Y0L5hiv5XzndNkgQeCZMkSeqEIUySJKkDhjBJkqQOGMIkSZI6YAiTJEnqgCFMkiSpA4YwSZKkDhjCJEmSOmAIkyRJ6oAhTJIkqQOGMEmSpA4YwiRJkjpgCJMkSeqAIUySJKkDhjBJkqQOtBrCklyT5KEkleTsNseWJEkaJW0fCbsVuBB4pOVxJUmSRsrxbQ5WVbcDJGlzWEmSpJHTaggbRpIVwIq+5vGFrkOSJKlLCx7CgE3A5g62K0mStGh0EcK2ADf0tY0DOxe+FEmSpG4seAirqsPA4d42zyGTJEnLTdu3qPhAkv1MHtm6LcnuNseXJEkaFa2GsKq6oqrGq+r4qvq/quqlbY4vSTOVZG2SO5LsaV7XDOjzG0l2J7k7yd8leU0XtUpaXrxjvqRRtw3YWlVrga3AdQP6/C1wflX9IPCLwMeSPHsBa5S0DBnCJI2sJKcC64DtTdN2YF2Slb39qurPq+pfmo/3AAFOWbBCJS1LXVwdKUkL5QzgQFVNAFTVRJJHm/ZDU6zzVuCBqtrf/4X3OZTUJkOYJDWSvBL4LeCiKbp4n0NJrTGESRpl+4BVScaao2BjwOlN+9MkeQVwE/CGqvrKFON5n0NJrTGESRpZVXUwyS7gUiYD1qXAXVX1tJ8ik5wPfAx4U1X9/THG8z6HklrjifmSRt1lwMYke4CNzWeS7EhyXtPng8CzgeuS7GqWl3VTrqTlwiNhkkZaVd0LXDCgfUPP+/MXtChJwiNhkiRJnTCESZIkdcAQJkmS1AFDmCRJUgcMYZIkSR0whEmSJHXAECZJktQBQ5gkSVIHDGGSJEkdMIRJkiR1wBAmSZLUAUOYJElSBwxhkiRJHTCESZIkdcAQJkmS1AFDmCRJUgcMYZIkSR0whEmSJHXAECZJktQBQ5gkSVIHDGGSJEkdMIRJkiR1wBAmSZLUgVZDWJK1Se5Isqd5XdPm+JI0U8PMS0nGkmxN8kCS+5O8vYtaJS0vbR8J2wZsraq1wFbgupbHl6SZGmZeejPwImAN8ArgyiSrF6pASctTayEsyanAOmB707QdWJdkZVvbkKSZmMG8dAlwfVUdqapDwK3AxQtXqaTl6PgWxzoDOFBVEwBVNZHk0ab90NFOSVYAK/rWPQtg//79LZYjaTHr+fd9bB43M9S8BJwJPNLzeW/T52mONX99+sIrOPWEk1osXdJidfDb3zj6dk7zV5shbFibgM2Dvli/fv0ClyJpEVgDPNB1EUOacv66Yt+nF7gUSYvAnOavNkPYPmBVkrHmr80x4PSmvdcW4Ia+thcCfwm8ksm/QJeycWAnsB5Y6of2RmVfRmU/YLT25Uzgr4AH53Ebw85Le5k8ovWFntoe4Xs5fy0d7sviMyr7AS3NX62FsKo6mGQXcClwU/N6V3N+RW+/w8Dh3rYkR9/uraqH26qpCz37st99WRxGZT9gZPflW/O1jWHnJeBm4B1JbgFOAd4IXDhgPOevJcJ9WXxGZT+gvfmr7asjLwM2JtkDbGw+S1KXBs5LSXYkOa/pcyOTf9HeB9wJXFVV83mETpLaPSesqu4FLmhzTEmai6nmpara0PN+Arh8IeuSJO+YL0mS1IHFEsIOA79J37kWS5T7sviMyn6A+7IYjcp+gPuyWI3KvozKfkBL+5KqaqccSZIkDW2xHAmTJElaVgxhkiRJHTCESZIkdcAQJkmS1AFDmCRJUgcMYZIkSR0whEmSJHXAECZJktQBQ5gkSVIHDGGakyTPTfKnSb6Z5JEk/3Wet/fMJB9qtvVkkruS/OR8blPS6OlqLlnoOVOL2/FdF6AlbyvwLeA04Fzgz5LcXVW752l7xwP7gFcCe4ENwMeTvKyqHp6nbUoaPV3NJQs9Z2oR89mRmrUkJwJfB86uqj1N243Agap69yzGS83iH8gk9wC/WVWfmOm6knTUXOaSYeavtudMLX3+HKm5WAtMHJ1MGncDL53pQEnWAzuSPGuG653W1OFfkZJmbS5zyQzmr9bmTI0GQ5jm4iTgib62J4CTZzHW54GDwKeGDWJJTgA+Cny4qu6dxTYlqY25ZNj5q805UyPAEKa5+AbwnL625wBPTrVCktcmqf4FmADeClwEXD7dhpMcB9zI5LkVvzLbHZC0vM1kLmlh/prxnKnR5on5mos9wPFJ1lTVfU3bD3KMw/lV9Vkg/e3NRPiHwPOBbcfaaJIAH2LyxNYNVfXt2ZUvaTmb6VzSwvw14zlTo80jYZq1qvomcAtwVZITk/wo8AYm/6qcqR9lciJ8Q1X96zR9/yfwEuD1Q/SVpKm0NZcMNX+1PGdqBHh1pOYkyXOBP2DyMPzjwLur6o9mOdYwVxedBTwMPAX8e89Xv1RVH53NdiUtP23PJcNe3d3mnKmlzxAmSZLUAX+OlCRJ6sBQISzJrUnubh7rsDPJuQP6jCXZmuSBJPcneXv75UrSzCS5JslDzZVsZ0/Rx/lL0oIb9urIt1XVEwBJ3sDk79nr+vq8GXgRsAY4BbgryW0+SkZSx24FfhfYeYw+zl+SFtxQR8KOBrDG9wFHBnS7BLi+qo5U1SEmJ76L516iJM1eVd1eVfum6eb8JWnBDX2fsCS/D7yayXukvHZAlzOBR3o+7wXOGDDOCmBFX/MzgBcC9zF50ztJo2+MyfsqfaGqnuq4FucvSTPRyvw1dAirqrcDJHkL8H4mnzg/G5uAzbNcV9LoWQ/c3nURQ3L+ktRrTvPXjO+YX1U3JvlfSU6pqsd7vtoLnAV8ofnc/5flUVuAG/razgI+t3PnTsbHx2dakqQlaP/+/axfvx7gq13XgvOXpBloa/6aNoQlOQn4j0fPqUjyeuBrzdLrZuAdSW5h8sTWNwIX9o9XVYeBw33bAGB8fJzVq1fPeCckLWmL4Sc85y9JszGn+WuYE/NPBG5O8qUku4BfZfIRD5VkR5Lzmn43Ag8yeV7EncBVVfXgXIqTpLlK8oEk+4Fx4LYku5t25y9JnZr2SFhV/RPwI1N8t6Hn/QTHfnq8JC24qroCuGJAu/OXpE55x3xJkqQOGMIkSZI6YAiTJEnqgCFMkiSpA4YwSZKkDhjCJEmSOmAIkyRJ6oAhTJIkqQOGMEmSpA4YwiRJkjpgCJMkSeqAIUySJKkDhjBJkqQOGMIkSZI6YAiTJEnqgCFMkiSpA4YwSZKkDhjCJEmSOjBtCEtySpIdSb6S5J4ktyRZOaDfDUn2J9nVLL8+PyVLkiQtfcMcCSvgd6rqxVV1DvAA8L4p+r6vqs5tlqtbq1KSJGnEHD9dh6r6GvC5nqY7gctnu8EkK4AVfc3jsx1PkiRpKZrROWFJjmMygH1qii6/luRLSW5N8pIp+mwCHupbds6kDkmSpKVu2iNhfX4P+AZw7YDvfh34alUdSfJW4LNJXlhVE339tgA39LWNYxCTJEnLyNBHwpJcA6wBLqmqI/3fV9WBo+1V9RHgJAb8zFhVh6vq4d4F2D/bHZCkY0myNskdSfY0r2sG9Dk1yZ81Fx/dm+SDSWb6R6okzchQISzJ1cDLgTdW1VNT9FnV8/41wARwoI0iJWkOtgFbq2otsBW4bkCf9wD/2Fx89DIm57ufWbgSJS1H0/6ll+SlTE5Qe4C/SQLwUFX9dJJdwIaqehT4cJLTgCPAPwP/par+ff5Kl6RjS3IqsA64qGnaDlybZGVVHerpWsDJzXmvzwSegX9ESppnw1wduRvIFN+d2/P+VS3WJUltOAM4cPTc1KqaSPJo094bwn4L+ATwVeBE4Nqq+nz/YF7dLalN3jFfkuBi4B7g+cAq4MIkbxrQz6u7JbXGE0+17FR97+vRZdDnI0ee3r//c2//qcYZtO4w2xn03VRjTLXMdaz+/sNud7oxDh4c7v+vOdoHrEoy1hwFGwNOb9p7bQR+sbm46IkknwR+HPiTvn5e3S2pNYsqhH384/C8533389H/aExn2H7D9u/9/lh9B3033br9bb3/oZxunEHrHmusY43R/9107ccaY9h1pgo/bbTNZD0tDk8+Of/bqKqDzbmrlwI3Na939Z0PBpNHtF4L/G2SZwCvAm4ZMN5h4HBvW3OerCTN2KIKYbt2wcknd12FpBFzGZMXDr0X+DrwVoAkO4D3VtUXmfyZcVuSLwFjwP8HXN9RvZKWiUUVwqSFdPQARvLdpbf9uOO+t89U/Y87buoxp1q3dztTjdNbw1T9plp61xm0L4OW6bZ39PvZLP1j/9M/wfbtw/1/NRdVdS9wwYD2DT3vH+C7V1BK0oJYVCHs4ovh9NNnt+6xfhGYya8Fg/q20TZd/2HeT/XdoD7TjTFo3anWm+p1pmNNNf5xx31v26Bxj9U20z7q3sMPd12BJHVrUYWwH/ohWL266yokSZLmn7eokCRJ6oAhTJIkqQOGMEmSpA4sqnPC2LYNnvvcrquYH4vtbPDFVs9MLYb651LDQte/UNsbdjsJPP74/NYiSYvc4gphe/fC17/edRWSFsJC3K1VkhYxf46UJEnqwOI6EvbOd8L4eNdVtG85PitnKezzXGrsav8Waruz2c5M1zlwYGHu1ipJi9TiCmGrV3ujMGm5OPHEriuQpE75c6QkSVIHDGGSJEkdmDaEJTklyY4kX0lyT5Jbkqwc0O8/JPlYkvuT3Jvkp+anZEmSpKVvmCNhBfxOVb24qs4BHgDeN6Dfu4Anq+pFwOuB309yUnulSpIkjY5pQ1hVfa2qPtfTdCdw1oCulwDbmnXuA74I/GQLNUqSJI2cGV0dmeQ44HLgUwO+PhN4pOfzXuCMAWOsAFb0NY/gfSkkSZKmNtNbVPwe8A3g2jlscxOweQ7rS5IkLXlDXx2Z5BpgDXBJVR0Z0GUvT/+Z8kxg34B+W4AX9C3rh61DkiRpFAx1JCzJ1cDLgddV1VNTdLsZ+CXgi0nWAOcDl/Z3qqrDwOG+8WdSsyRJ0pI3zC0qXgq8Bzgd+Jsku5L8afPdriSnN13fD6xIcj/wGeCdVeUTeiVJkgaY9khYVe0GBh6qqqpze95/E7i4vdIkSZJGl3fMlyRJ6oAhTJIkqQOGMEkjLcnaJHck2dO8rpmi388m+VKSLzevpy10rZKWF0OYpFG3DdhaVWuBrcB1/R2SnAdcCVxUVWcDPwY8sZBFSlp+DGGSRlaSU4F1wPamaTuwLsnKvq6/ClxTVY8BVNUTVfVvC1eppOVopnfMl6Sl5AzgQFVNAFTVRJJHm/ZDPf1+AHgoyV8DJwG3AFdXVfUO5mPXJLXJECZJk3PhOcBFwDOAzzL5FJCP9PXzsWuSWuPPkZJG2T5gVZIxgOb1dL73kWqPAH9SVU81N5n+JPDDA8bzsWuSWmMIkzSyquogsIvvPkLtUuCuqjrU1/WPgFdn0gnATwB3DxjvcFU93LsA++dvDySNMkOYpFF3GbAxyR5gY/OZJDuaqyIB/hg4CPwDk6FtN/ChDmqVtIx4TpikkVZV9wIXDGjf0PP+CPBrzSJJC8IjYZIkSR0whEmSJHXAECZJktQBQ5gkSVIHDGGSJEkdMIRJkiR1wBAmSZLUgaFCWJJrkjyUpJKcPUWfK5McTLKrWba2W6okSdLoGPZmrbcCvwvsnKbfR6rqXXMrSZIkafQNFcKq6naAJPNbjSRJ0jLR9mOLfi7Jq4HHgM1VdUd/hyQrgBV9zeMt1yFJkrSotRnCtgFXV9W3k1wEfDLJS6rq8b5+m4DNLW5XkiRpyWnt6siqeqyqvt28/wtgHzDoJP4twAv6lvVt1SFJkrQUtHYkLMmqqjrQvD8XWA18pb9fVR0GDvet21YZkiRJS8Kwt6j4QJL9TJ67dVuS3U37jiTnNd1+O8mXk9wNXA+8paoem5eqJUmSlrhhr468ArhiQPuGnvdva7EuSZKkkeYd8yVJkjpgCJMkSeqAIUySJKkDhjBJkqQOGMIkSZI6YAiTJEnqgCFMkiSpA4YwSSMtydokdyTZ07yuOUbfFyf5lyTXLGSNkpYnQ5ikUbcN2FpVa4GtwHWDOiUZa767dQFrk7SMtfbsSElabJKcCqwDLmqatgPXJllZVYf6ur8b+AxwUrMMGm8FsKKveby9iiUtJx4JkzTKzgAOVNUEQPP6aNP+HUnOAV4D/I9pxtsEPNS37Gy5ZknLhEfCJC1rSU4Argd+oaomkhyr+xbghr62cQxikmbBECZplO0DViUZawLWGHB6037U84HvB3Y0AWwFkCTPqap39g5WVYeBw71t04Q2SZqSIUzSyKqqg0l2AZcCNzWvd/WeD1ZVe4HnHf2c5ErgpKp61wKXK2mZ8ZwwSaPuMmBjkj3AxuYzSXYkOa/TyiQtax4JkzTSqupe4IIB7Rum6H/lfNckSeCRMEmSpE4YwiRJkjowbQhLck2Sh5JUkrOn6DOWZGuSB5Lcn+Tt7ZcqSZI0OoY5EnYrcCHwyDH6vBl4EbAGeAVwZZLVcy1OkiRpVE17Yn5V3Q7T3gvnEuD6qjoCHEpyK3Ax8P7+jj72Q5Ikqb2rI8/k6UfK9tL3WJAem4DNLW1XkiRpSeriFhU+9kOSJC17bYWwvcBZwBeaz/1Hxr7Dx35IkiS1d4uKm4F3JDkuyUrgjcAnWhpbkiRp5Axzi4oPJNnP5E+GtyXZ3bT3PvLjRuBB4D7gTuCqqnpwnmqWJEla8oa5OvIK4IoB7Rt63k8Al7dbmiRJ0ujyjvmSJEkdMIRJkiR1wBAmSZLUAUOYJElSBwxhkiRJHTCESZIkdcAQJkmS1AFDmCRJUgcMYZIkSR0whEmSJHXAECZJktQBQ5gkSVIHDGGSRlqStUnuSLKneV0zoM9vJNmd5O4kf5fkNV3UKml5MYRJGnXbgK1VtRbYClw3oM/fAudX1Q8Cvwh8LMmzF7BGScvQ8V0XIEnzJcmpwDrgoqZpO//N2gMAABsaSURBVHBtkpVVdehov6r6857V7gECnALs7xtvBbCibzPjbdctaXkwhEkaZWcAB6pqAqCqJpI82rQfmmKdtwIPVNX+Ad9tAjbPS6WSlh1DmCQ1krwS+C2+e+Ss3xbghr62cWDnPJYlaUQNFcKSrAU+zOTh+ceBt1bVfX19rgR+GXi0afp8Vf239kqVpBnbB6xKMtYcBRsDTm/anybJK4CbgDdU1VcGDVZVh4HDfeu1X7WkZWHYE/OHObEV4CNVdW6zGMAkdaqqDgK7gEubpkuBu3rPBwNIcj7wMeBNVfX3C1ulpOVq2hDWc2Lr9qZpO7Auycr5LEySWnIZsDHJHmBj85kkO5Kc1/T5IPBs4Loku5rlZd2UK2m5GObnyJmc2PpzSV4NPAZsrqo7+gfz6iJJC6mq7gUuGNC+oef9+QtalCTR7on524Crq+rbSS4CPpnkJVX1eF8/ry6SJEnL3jDnhH3nxFaAqU5srarHqurbzfu/aL4/e8B4W4AX9C3rZ7sDkiRJS9G0R8Kq6mCSoye23sTUJ7auqqoDzftzgdXA91xh5NVFkiRJw/8ceRnw4STvBb7O5M0MSbIDeG9VfRH47SQvByaAbwFvqarH5qFmSZKkJW+oEDbkia1va7EuSZKkkeYDvCVJkjpgCJMkSeqAIUySJKkDhjBJkqQOGMIkSZI6YAiTJEnqgCFMkiSpA4YwSZKkDhjCJEmSOmAIkyRJ6oAhTJIkqQOGMEmSpA4YwiRJkjpgCJMkSeqAIUySJKkDhjBJkqQOGMIkSZI6YAiTJEnqwFAhLMnaJHck2dO8rhnQZyzJ1iQPJLk/ydvbL1eSZsb5S9JiNeyRsG3A1qpaC2wFrhvQ583Ai4A1wCuAK5OsbqFGSZoL5y9Ji9Lx03VIciqwDrioadoOXJtkZVUd6ul6CXB9VR0BDiW5FbgYeH/feCuAFX2bOQtg//79s9oJSUtPz7/vY/O1jYWcvz594RWcesJJ87AXkhabg9/+xtG3c5q/pg1hwBnAgaqaAKiqiSSPNu29k9iZwCM9n/c2ffptAjYP2tD69euHqVnSaFkDPDBPYy/Y/HXFvk+3UrCkJWVO89cwIaxtW4Ab+tpeCPwl8EomJ7+lbBzYCawHlvqhvVHZl1HZDxitfTkT+Cvgwa4LmQHnr6XDfVl8RmU/oKX5a5gQtg9YlWSs+StyDDi9ae+1l8nD8l/oKfCRvj5U1WHgcG9bku+MUVUPD139ItSzL/vdl8VhVPYDRnZfvjWPm3H+moER/efLfVkkRmU/oL35a9oT86vqILALuLRpuhS4q+98CoCbgXckOS7JSuCNwCfmUpwkzYXzl6TFbNirIy8DNibZA2xsPpNkR5Lzmj43MnlY7j7gTuCqqlpKPzNIGk3OX5IWpaHOCauqe4ELBrRv6Hk/AVzeXmmSNHfOX5IWq8Vyx/zDwG/Sd67FEuW+LD6jsh/gvixGo7If4L4sVqOyL6OyH9DSvqSq2ilHkiRJQ1ssR8IkSZKWFUOYJElSBwxhkiRJHTCESZIkdcAQJkmS1AFDmCRJUgcMYZIkSR0whEmSJHXAECZJktQBQ5gkSVIHDGGakyS/kuSLSZ5KcsMCbO+ZST6U5JEkTya5K8lPzvd2JY2WruaShZ4ztbgd33UBWvIeBf478Brg2QuwveOBfcArgb3ABuDjSV5WVQ8vwPYljYau5pKFnjO1iBnCNCdVdQtAkvOA8bmMlSQ1zRPlq+qbwJU9TZ9J8hDwcuDhuWxf0vLR9lwyzPzVbLe1OVNLnz9HalFIsh7YkeRZM1zvNGAtsHteCpO0LMxlLpnt/CUZwrRYfB44CHxq2IksyQnAR4EPV9W981mcpNHVwlwy4/lLAkOYFliS1yap/gWYAN4KXARcPsQ4xwE3At8CfmVei5Y0smYyl7Q1f0lHeU6YFlRVfRZIf3szEf4h8Hxg27HGSBLgQ8BpwIaq+vY8lCppxM10Lmlj/pJ6GcI0J0mOZ/KfozFgrDkU/+9V9e8zHOpHmZwI31BV/zpN3/8JvAR41RB9JWkqbc0lQ89fLc6ZGgEZ4mIOaUpJrgQ29zX/ZlVdOYuxpr26KMlZTF659BTQO2n9UlV9dKbblLQ8tT2XDHt1ZJtzppY+Q5gkSVIHPDFfkiSpA0OFsCS3Jrm7eazDziTnDugzlmRrkgeS3J/k7e2XK0kzk+SaJA81V7KdPUUf5y9JC27YE/PfVlVPACR5A/AHwLq+Pm8GXgSsAU4B7kpym4+SkdSxW4HfBXYeo4/zl6QFN1QIOxrAGt8HHBnQ7RLg+qo6AhxKcitwMfD+3k5JVgAr+tZ9BvBC4D4m77ciafSNMXlJ/xeq6qn52khV3Q4weTeCKTl/SZqJVuavoW9RkeT3gVczeY+U1w7ocibwSM/nvcAZA/pt4nuvDJG0fK0Hbu+4BucvSbMxp/lr6BBWVW8HSPIWJv863DDLbW4BbuhrOwv43M6dOxkf93mm0nKwf/9+1q9fD/DVrmuZAecvSa3NXzO+WWtV3ZjkfyU5paoe7/lqL5OT0Reaz/1/WR5d/zBwuLft6M8E4+PjrF69eqYlSVraFsNPeM5fkmZjTvPXtFdHJjkpyRk9n18PfK1Zet0MvCPJcUlWAm8EPjGX4iRpgTh/SVpww9yi4kTg5iRfSrIL+FXg9VVVSXYkOa/pdyPwIJMnp94JXFVVD85L1ZI0pCQfSLIfGAduS7K7aXf+ktSpaX+OrKp/An5kiu829LyfwKfHS1pkquoK4IoB7c5fkjrlHfMlSZI6YAiTJEnqgCFMkiSpA4YwSZKkDhjCJEmSOmAIkyRJ6oAhTJIkqQOGMEmSpA4YwiRJkjpgCJMkSeqAIUySJKkDhjBJkqQOGMIkSZI6YAiTJEnqgCFMkiSpA4YwSZKkDhjCJEmSOmAIkyRJ6sC0ISzJKUl2JPlKknuS3JJk5YB+NyTZn2RXs/z6/JQsSZK09A1zJKyA36mqF1fVOcADwPum6Pu+qjq3Wa5urUpJkqQRc/x0Harqa8DnepruBC6f7QaTrABW9DWPz3Y8SZKkpWjaENYryXFMBrBPTdHl15L8EpNHy/7fqvrHAX02AZsHrv1//g+cdtpMSppaVTvjLIYajjXOXLYxk3XbqGG6foO+n806x+o31/HaGKf3u5nu82zWnW6duaw7l3Uef3zqWiRpGZhRCAN+D/gGcO2A734d+GpVHUnyVuCzSV5YVRN9/bYAN/S1jQM7ue02OPnkGZYkaUl68smuK5CkTg19dWSSa4A1wCVVdaT/+6o6cLS9qj4CnMSAnxmr6nBVPdy7APtnuwOSdCxJ1ia5I8me5nXNgD6nJvmz5uKje5N8MMlM/0iVpBkZapJJcjXwcuB1VfXUFH1WVdWB5v1rgAngwIyq+U//qb2fIycLaW+s+RivrbGnW3emY8+mlkHrDDvObOqfyz63Nd5c2gZ9P9ca2q7xWG1T1TqTdQ4cgO3bB4/Trm3A1qq6KcnPA9cB/7mvz3uAf6yq1yU5Abgd+Bng4wtRoKTladoQluSlTE5Qe4C/yeRE+lBV/XSSXcCGqnoU+HCS04AjwD8D/6Wq/n1G1bz61bB69cz2QNLS9H3fN++bSHIqsA64qGnaDlybZGVVHerpWsDJzXmvzwSewYA/Ir2wSFKbhrk6cjcw8M/eqjq35/2rWqxLktpwBnDg6LmpVTWR5NGmvTeE/RbwCeCrwInAtVX1+QHjTX1hkSTNkHfMlyS4GLgHeD6wCrgwyZsG9NsCvKBvWb9QRUoaLZ54KmmU7QNWJRlrjoKNAac37b02Ar/YXFz0RJJPAj8O/Elvp6o6DBzubct8nisqaaR5JEzSyKqqg8Au4NKm6VLgrr7zwQAeAl4LkOQZwKuALy9UnZKWJ0OYpFF3GbAxyR4mj3hdBtA8E/e8ps8mYH2SLzEZ2vYA13dRrKTlw58jJY20qroXuGBA+4ae9w/w3SsoJWlBeCRMkiSpA4YwSZKkDhjCJEmSOmAIkyRJ6oAhTJIkqQOGMEmSpA4sqltU/Nu/wb/+a9dVtGsUb6a9GPZpLjUsVP1tb2fY8bra7kKNI0mjYlGFsKuugpNP7roKSQvhySe7rkCSuuXPkZIkSR1YVEfCnvlMeNazuq5i6avquoLpzaXGhd6/+dzesGO3XUNb4y2Ff9YkabFaVCFs82ZYvbrrKiQthIcfhu3bu65Ckrrjz5GSJEkdmDaEJTklyY4kX0lyT5Jbkqwc0O8/JPlYkvuT3Jvkp+anZEmSpKVvmCNhBfxOVb24qs4BHgDeN6Dfu4Anq+pFwOuB309yUnulSpIkjY5pQ1hVfa2qPtfTdCdw1oCulwDbmnXuA74I/GQLNUqSJI2cGZ2Yn+Q44HLgUwO+PhN4pOfzXuCMAWOsAFb0NY/PpA5JkqSlbqZXR/4e8A3g2jlscxOweQ7rS5IkLXlDXx2Z5BpgDXBJVR0Z0GUvT/+Z8kxg34B+W4AX9C3rh61DkiRpFAx1JCzJ1cDLgddV1VNTdLsZ+CXgi0nWAOcDl/Z3qqrDwOG+8WdSsyRJ0pI3zC0qXgq8Bzgd+Jsku5L8afPdriSnN13fD6xIcj/wGeCdVeXT4SRJkgaY9khYVe0GBh6qqqpze95/E7i4vdIkSZJGl3fMlyRJ6oAhTJIkqQOGMEkjLcnaJHck2dO8rpmi388m+VKSLzevpy10rZKWF0OYpFG3DdhaVWuBrcB1/R2SnAdcCVxUVWcDPwY8sZBFSlp+ZnqzVklaMpKcCqwDLmqatgPXJllZVYd6uv4qcE1VPQZQVQMDmE/8kNQmQ5ikUXYGcKCqJgCqaiLJo017bwj7AeChJH8NnATcAlxdVdU3nk/8kNQaQ5gkTc6F5zB5xOwZwGeZfArIR/r6bQFu6GsbB3bOc32SRpAhTNIo2wesSjLWHAUbY/LG0/2PVHsE+JPmiSBPJfkk8MP0hTCf+CGpTZ6YL2lkVdVBYBfffYTapcBdfeeDAfwR8OpMOgH4CeDuhatU0nJkCJM06i4DNibZA2xsPpNkR3NVJMAfAweBf2AytO0GPtRBrZKWEX+OlDTSqupe4IIB7Rt63h8Bfq1ZJGlBeCRMkiSpA4YwSZKkDhjCJEmSOmAIkyRJ6oAhTJIkqQOGMEmSpA4YwiRJkjowVAhLck2Sh5JUkrOn6HNlkoNJdjXL1nZLlSRJGh3D3qz1VuB3mf4htR+pqnfNrSRJkqTRN1QIq6rboZ0H1SZZAazoax6f88CSJElLSNvnhP1cknuS/O8kr5iizybgob5luiNskiRJI6XNELYNeEFVnQO8H/hkklMG9NsCvKBvWd9iHZIkSYteaw/wrqrHet7/RZJ9wNnAX/X1Owwc7m1r42dOSZKkpaS1I2FJVvW8PxdYDXylrfElSZJGybC3qPhAkv1MnkB/W5LdTfuOJOc13X47yZeT3A1cD7yl9+iYJEmSvmvYqyOvAK4Y0L6h5/3bWqxLkiRppHnHfEmSpA4YwiRJkjpgCJMkSeqAIUySJKkDhjBJkqQOGMIkSZI6YAiTJEnqgCFM0khLsjbJHUn2NK9rjtH3xUn+Jck1C1mjpOXJECZp1G0DtlbVWmArcN2gTknGmu9uXcDaJC1jrT3AW5IWmySnAuuAi5qm7cC1SVZW1aG+7u8GPgOc1CyDxlsBrOhrHm+vYknLiUfCJI2yM4ADVTUB0Lw+2rR/R5JzgNcA/2Oa8TYBD/UtO1uuWdIy4ZEwSctakhOA64FfqKqJJMfqvgW4oa9tHIOYpFkwhEkaZfuAVUnGmoA1BpzetB/1fOD7gR1NAFsBJMlzquqdvYNV1WHgcG/bNKFNkqZkCJM0sqrqYJJdwKXATc3rXb3ng1XVXuB5Rz8nuRI4qaretcDlSlpmPCdM0qi7DNiYZA+wsflMkh1Jzuu0MknLmkfCJI20qroXuGBA+4Yp+l853zVJEngkTJIkqROGMEmSpA5MG8KSXJPkoSSV5Owp+owl2ZrkgST3J3l7+6VKkiSNjmGOhN0KXAg8cow+bwZeBKwBXgFcmWT1XIuTJEkaVdOemF9Vt8O098K5BLi+qo4Ah5LcClwMvL+/o4/9kCRJau/qyDN5+pGyvfQ9FqTHJmBzS9uVJElakrq4RYWP/ZAkScteWyFsL3AW8IXmc/+Rse/wsR+SJEnt3aLiZuAdSY5LshJ4I/CJlsaWJEkaOcPcouIDSfYz+ZPhbUl2N+29j/y4EXgQuA+4E7iqqh6cp5olSZKWvGGujrwCuGJA+4ae9xPA5e2WJkmSNLq8Y74kSVIHDGGSJEkdMIRJkiR1wBAmSZLUAUOYJElSBwxhkiRJHTCESZIkdcAQJkmS1AFDmCRJUgcMYZIkSR0whEmSJHXAECZJktQBQ5ikkZZkbZI7kuxpXtcM6PMbSXYnuTvJ3yV5TRe1SlpeDGGSRt02YGtVrQW2AtcN6PO3wPlV9YPALwIfS/LsBaxR0jJ0fNcFSNJ8SXIqsA64qGnaDlybZGVVHTrar6r+vGe1e4AApwD7+8ZbAazo28x423VLWh4MYZJG2RnAgaqaAKiqiSSPNu2HpljnrcADVbV/wHebgM3zUqmkZccQJkmNJK8EfovvHjnrtwW4oa9tHNg5j2VJGlFDhbAka4EPM3l4/nHgrVV1X1+fK4FfBh5tmj5fVf+tvVIlacb2AauSjDVHwcaA05v2p0nyCuAm4A1V9ZVBg1XVYeBw33rtVy1pWRj2xPxhTmwF+EhVndssBjBJnaqqg8Au4NKm6VLgrt7zwQCSnA98DHhTVf39wlYpabmaNoT1nNi6vWnaDqxLsnI+C5OkllwGbEyyB9jYfCbJjiTnNX0+CDwbuC7JrmZ5WTflSlouhvk5ciYntv5cklcDjwGbq+qO/sG8ukjSQqqqe4ELBrRv6Hl//oIWJUm0e2L+NuDqqvp2kouATyZ5SVU93tfPq4skSdKyN8w5Yd85sRVgqhNbq+qxqvp28/4vmu/PHjDeFuAFfcv62e6AJEnSUjTtkbCqOpjk6ImtNzH1ia2rqupA8/5cYDXwPVcYeXWRJEnS8D9HXgZ8OMl7ga8zeTNDkuwA3ltVXwR+O8nLgQngW8BbquqxeahZkiRpyRsqhA15YuvbWqxLkiRppPkAb0mSpA4YwiRJkjpgCJMkSeqAIUySJKkDhjBJkqQOGMIkSZI6YAiTJEnqgCFMkiSpA4YwSZKkDhjCJEmSOmAIkyRJ6oAhTJIkqQOGMEmSpA4YwiRJkjpgCJMkSeqAIUySJKkDhjBJkqQOGMIkSZI6MFQIS7I2yR1J9jSvawb0GUuyNckDSe5P8vb2y5WkmXH+krRYDXskbBuwtarWAluB6wb0eTPwImAN8ArgyiSrW6hRkubC+UvSonT8dB2SnAqsAy5qmrYD1yZZWVWHerpeAlxfVUeAQ0luBS4G3t833gpgRd9mzgLYv3//rHZC0tLT8+/72HxtYyHnr09feAWnnnDSPOyFpMXm4Le/cfTtnOavaUMYcAZwoKomAKpqIsmjTXvvJHYm8EjP571Nn36bgM2DNrR+/fphapY0WtYAD8zT2As2f12x79OtFCxpSZnT/DVMCGvbFuCGvrYXAn8JvJLJyW8pGwd2AuuBpX5ob1T2ZVT2A0ZrX84E/gp4sOtCZsD5a+lwXxafUdkPaGn+GiaE7QNWJRlr/oocA07//9u7m9C4qjCM4/+nLajQVW1VimBFpRWDiFb8KupGRFEMCkqx4ErQRdciCMWCIFgQVBQE3SiIivVjofhB7UY0+JHWKthUaCpZiFQNWPED9XVxzqRjNJOJmeacM3l+cEnuTUjeJznzcu7MvXPy8W7fkJ6W/7irwCOzvoeImAamu49JmvkZETHZd/UV6soy5Sx1GJYcMLRZfj+Bv8b9awGGdHw5SyWGJQcMrn/Ne2F+RHwH7AO25kNbgfFZ11MAvAzcLWmFpHXAKPDKYoozM1sM9y8zq1m/d0feA2yXNAFsz/tIelPS5vw9z5GeljsEfATsjIiWXmYws+Hk/mVmVerrmrCI+Aq47D+O39j1+Z/AvYMrzcxs8dy/zKxWtbxj/jTwILOutWiUs9RnWHKAs9RoWHKAs9RqWLIMSw4YUBZFxGDKMTMzM7O+1fJMmJmZmdmy4kmYmZmZWQHFJ2H9LK5bK0mn5jusDkr6XNLufHs7ki6XtD/neicvn1I9STskhaSRvN9cDkknS3pK0iFJByQ9nY83N9Yk3SRpXNK+PMZuzcerziJpl6TD3WMpH5+z7tozzaXhut2/KuT+VYcl62ERUXQD9gDb8ufbgD2la1pA7WuAa7v2HwGeAQR8DWzJxx8Ani1dbx95LgbeIr1J5UjDOR4DHuX4NY+n549NjbX89/8RGMn7FwI/kU6eqs4CbCEt+zPZqX++/0HtmXpkbbVu968KN/evOral6mGlQ55GurNgZd5fmffXlf4H/M88twHvAZcCX3QdXwscK13fPLWfBHwInN0ZdI3mWJ3H0OpZx5sba7mJfQ9clfevBiZaytLdwHrV3VKm1sdVjyzuX+VzuH9Vtp3oHlb65ch/La4LdBbXbYqkFaT3GXqDWUueRMRRYIWkNYXK68dO4PmIONx1rMUc55Ae+DskfSJpr6TOGU1TYy3SI/l24HVJR4DXgLtoMEvWq+5hzNQM969quH/VbeA9rPQkbJg8DhwDnihdyEJJuoJ01vhk6VoGYBVpQeXxiNgM3AfsJp1hNkXSKuB+4JaIOAu4GXiRBrNY9dy/6uD+tcyUnoTNLK4LoLkX162apF3AecAdEfEXxxcD7nx9LenE4IdCJc7nGmATcFjSJGml+7eBc2krB6Qz3z+AFwAiYgw4CvxCe2PtImB9RHwAkD/+DPxKe1mg9+O91V7Qat0z3L+q4v5Vt4H3sKKTsOh/cd1qSXoIuAQYjYjf8uFPgVPy08iQ1qp7qUR9/YiIhyNifURsiIgNwBRwPelC3WZywMxLDu8D10G6W4X0Wv0E7Y21KeBMSRsBJJ0PnEFa37C1LD0f7632glbr7nD/qov7V91OSA+r4KK3TcAYaZCNARtL17SA2i8AAjiY//j7gFfz164EDpAG3LvkO1xa2PjnhYjN5SA9nb831/0ZcEOrYw24M+fYn7fRFrKQ7vCaIp3Vfwt8OV/dtWfqkbXVut2/Ktzcv+rYlqqHedkiMzMzswJKXxNmZmZmtix5EmZmZmZWgCdhZmZmZgV4EmZmZmZWgCdhZmZmZgV4EmZmZmZWgCdhZmZmZgV4EmZmZmZWwN++aulDyA8WJgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x720 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(nrows=3, ncols=2, sharex=True, figsize=(10, 10))\n",
    "\n",
    "for i, (x, y) in enumerate([(0,1), (1,0), (0,2), (2,0), (1,2), (2,1)]):\n",
    "        x_range = np.arange(len(granger1[x, y, :]))\n",
    "        ax.flatten()[i].plot(granger1[x, y, :], color='b', lw=3, alpha=0.6)\n",
    "        ax.flatten()[i].plot(granger2[x, y, :], color='r', lw=3, alpha=0.6)\n",
    "        ax.flatten()[i].set_title(r\"{} $\\rightarrow$ {}\".format(x, y))\n",
    "        ax.flatten()[i].set_xlim((0, 100))\n",
    "        if i == 0:\n",
    "            ax.flatten()[i].set_ylim((0, 4))\n",
    "        else:\n",
    "            ax.flatten()[i].set_ylim((0, 3))\n",
    "        if i % 2 == 1:\n",
    "            ax.flatten()[i].set_ylim((0, 1))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Replication of Figure 1b, 2b from Chen et al.**\n",
    "Granger causality spectra by pairwise analysis. The blue lines mark the influence for the case of differentially delayed drive, the red lines for sequential drive. From- and to-channel are indicated above each figure. The influences mimic very well the results from Chen et al., where a slight deviation on the $y$-scale may result from a different resolution when transitioning to the spectral domain. The value chosen here was 1024 Hz, whereas Chen et al. do not document a resolution. Note that in the case of sequential drive, we would prefer not to see a direct influence $0 \\rightarrow 2$, whereas for differentially delayed drive there should not be any direct influence $1 \\rightarrow 2$. At least the two different cases of drive should be easier to distinguish."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/media/dan/usbdata1/conda_envs/py3/lib/python3.7/site-packages/ipykernel_launcher.py:3: RuntimeWarning: Mean of empty slice\n",
      "  This is separate from the ipykernel package so we can avoid doing imports until\n",
      "/media/dan/usbdata1/conda_envs/py3/lib/python3.7/site-packages/ipykernel_launcher.py:4: RuntimeWarning: Mean of empty slice\n",
      "  after removing the cwd from sys.path.\n"
     ]
    }
   ],
   "source": [
    "# we have nan on the main diagonal, because we do not calculate relective influence of a channel onto iteself\n",
    "# this causes the warnings, which can be ignored.\n",
    "cond_granger_mean1 = np.nanmean(cond_granger_full_1, axis=-2).transpose(1, 0 ,2)\n",
    "cond_granger_mean2 = np.nanmean(cond_granger_full_2, axis=-2).transpose(1, 0 ,2)\n",
    "# the transposition is necessary, because internally, we compute the influence of y on x via z, not x on y via z.\n",
    "# This is not (yet) consistent with the original GrangerAnalyzer syntax."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAJTCAYAAABAR8HtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdf7Ddd13v++erO/y6NJxca8ul3W3Dj4SDQGVCa2UgcBgtSIQLR0WsaM/oBW0dw0THO5ebAzexSofBnjk9SrS1By1SjAinFi5E8OI51RTLEa9JS6slpTRN04KJeNOTohZM3/eP9Y2urq6d7CRrr8/aaz0fM99Za332Z3++7y80n3ntz/r+SFUhSZKk8TqtdQGSJEmzyBAmSZLUgCFMkiSpAUOYJElSA4YwSZKkBgxhkiRJDRjCJEmSGjCESZIkNWAI07KVZGuSra3rkKQT5fwlMIRpAUm+LckfJPlGkvuT/OgY9vmzSf4iyaNJbjjFsZ6S5ANd7YeT7EryuhGVKmmCjXIuOYF9jmzOdP6aHStaF6CJtQ34JvBM4CXAp5LcXlV3LeE+HwJ+GXgt8LRTHGsF8ADwKmAfsAH4/SQvrqq9pzi2pMk2yrlksUY5Zzp/zQhXwvQESZ4O/CDw7qp6pKpuBT4B/PhJjpfF9Kuqm6rqZuDrJ7OfgbG+UVVbq2pvVT1WVZ8E7gNeeqpjS5pso5xLFjN/jXrOdP6aHYYwDbMWOFJVe/rabgdeeKIDJVkP7Ejy1FEVdzKSPJPecS3lSp6kKXIC89fI5swF6nD+mlKGMA1zOvDwQNvDwMqTGOtzwAHgE62CWJInAR8GPlhVd7eoQdKytNj5a5Rz5uM4f003Q5iGeQR4xkDbM4DDC/1Cku9LUoMbcAS4DLgEuGLJKl64rtOAD9E7V+Nnx71/SZNvBPPXCc+Zi6zL+WvKeWK+htkDrEiypqru6dq+k2MshVfVp4EnnDvRTSK/DTwLuHYJal1Qdy7HB+idKLuhqr41zv1LWh5GMH+d8Jx5PM5fs8GVMD1BVX0DuAm4MsnTk7wceCO9v8hO1MvpTSJvrKp/OFbHJCu6Jf85YC7JU5Ocyh8KvwG8AHjD8fYtaXqMcC5Z1Pw14jnzKOevGWAI00J+ht6l3QeA7cAVJ3OpdVXtBF63yEnkXcA/AO8Efqx7/64T3SdAkvOBn6Z3qfjXkjzSbW89mfEkLSsjmUtOcP4ayZwJzl+zJFXVugbppBy923RVbW1biSSdGOcvgSthkiRJTRjCtJzd0m3SgpJcneS+7oq3Fy3QZy7JtiT3JvlykreNu07NnFtw/pp5fh0paaoleQVwP7ATeH1V3Tmkz2XAW4HXAWcAu4BX+IgYSUvJlTBJU62qbq2qB47T7S3A9d0jYg4CNwNvXvrqJM2yibhPWJKnABcBX6V3czxJ02+O3v2XvlBVjzau5Tx6q2VH7QPOHeyUZBWwaqD5ycBzgHtw/pJmxUjmr4kIYfQC2M7WRUhqYj1wa+siFmkTsKV1EZImxinNX5MSwr4KsHPnTubn51vXImkM9u/fz/r166H799/YPuB84Avd58GVsaOuAW4YaDsfuMX5S5odo5q/JiWEHQGYn59n9erVjUuRNGaT8BXeR4G3J7mJ3on5bwJeOdipqg4Bh/rbek+Xcf6SZtQpzV+emC9pqiX51ST7gXngs0nu6tp3JLmw6/Yh4Cv0zuv6PHBlVX2lScGSZsakrIRJ0pKoqncA7xjSvqHv/RHginHWJUmuhEmSJDUwWSthmzfDypWtq5A0DocPt65AkppyJUySJKkBQ5gkSVIDk/V15FVXgZd4S7Nh717Yvr11FZLUjCthkiRJDRjCJEmSGjCESZIkNWAIkyRJasAQJkmS1IAhTJIkqQFDmCRJUgOGMEmSpAYMYZIkSQ0YwiRJkhowhEmSJDVgCJMkSWrAECZJktSAIUySJKkBQ5gkSVIDhjBJkqQGDGGSJEkNGMIkSZIaMIRJmmpJ1ia5Lcme7nXNkD5nJflUkjuS3J3k15OsaFGvpNlhCJM07a4FtlXVWmAbcN2QPpuBv66qC4AXAy8FfmB8JUqaRYYwSVMryVnAOmB717QdWJfkzIGuBaxMchrwFODJwINjK1TSTFpUCEtyc5Lbk+xKsjPJS4b0mUuyLcm9Sb6c5G2jL1eSTsi5wINVdQSge32oa+/3S8Ba4KvA14DPVNXnBgdLsirJ6v4NmF/C+iVNscWe8/DvquphgCRvBH6L3l+X/d4KPA9YA5wB7Ery2araO6JaJWmpvBm4A/geYCXwh0l+qKo+NtBvE7Bl3MVJmk6LWgk7GsA6/wp4bEi3twDXV9VjVXUQuJnexCZJrTwAnJNkDnor9sDZXXu/jcCHu/nrYeDjwKuHjHcN8OyBbf0S1S5pyi366p8k/xl4DRDg+4Z0OQ+4v+/zPp645E+SVcCqgWaX8yWNXFUdSLIbuBS4sXvd1f2h2O8+evPanyd5MvC9wE1DxjsEHOpvS7IUpUuaAYs+Mb+q3lZV59G7iuhXTmGfm+hNeP3bzlMYT5KO5XJgY5I99Fa8LgdIsiPJhV2fTcD6JF8EdgN7gOtbFCtpdpzwfXCq6kNJfjPJGVX19b4f7QPOB77QfR5cGTvqGuCGgbZ5DGKSlkBV3Q1cPKR9Q9/7e4FLxlmXJB03hCU5Hfifq+qB7vMbgL/rtn4fBd6e5CZ6J+a/CXjl4Hgu50uSJC1uJezpwEeTPB04Qi98vaGqKskO4P+qqr8APkTvr817ut+7sqq+shRFS5IkLXfHDWFV9TfAdy/ws/7l/CPAFaMrTZIkaXp5x3xJkqQGDGGSJEkNGMIkSZIaMIRJkiQ1YAiTJElqwBAmSZLUgCFMkiSpAUOYJElSA4YwSZKkBgxhkiRJDRjCJEmSGjCESZIkNWAIkyRJasAQJkmS1IAhTJIkqQFDmCRJUgOGMEmSpAYMYZIkSQ0YwiRJkhowhEmSJDVgCJM01ZKsTXJbkj3d65oF+v1wki8mubN7fea4a5U0WwxhkqbdtcC2qloLbAOuG+yQ5EJgK3BJVb0IeAXw8DiLlDR7DGGSplaSs4B1wPauaTuwLsmZA11/Dri6qr4GUFUPV9U/jq9SSbNoResCJGkJnQs8WFVHAKrqSJKHuvaDff2+A7gvyZ8CpwM3Ae+pquofLMkqYNXAPuaXqnhJ080QJkm9ufAC4BLgycCngX3A7wz02wRsGW9pkqaVX0dKmmYPAOckmQPoXs/u2vvdD3ysqh6tqsPAx4HvGjLeNcCzB7b1S1S7pCl33BCW5IwkO5J8KckdSW4acj4FSW5Isj/J7m7790tTsiQtTlUdAHYDl3ZNlwK7qurgQNffBV6TnicB3wPcPmS8Q1W1t38D9i/dEUiaZotZCSvgfVX1/Kq6ALgXeO8Cfd9bVS/ptveMrEpJOnmXAxuT7AE2dp/p/ri8sOvze8AB4K/ohba7gA80qFXSDDnuOWFV9XfALX1NnweuWKqCJGmUqupu4OIh7Rv63j8G/Hy3SdJYnNA5YUlOoxfAPrFAl5/vbnJ4c5IXLDDGqiSr+ze8ukiSJM2YE7068teAR4D3D/nZvwe+WlWPJbkM+HSS5xy9NLyPVxdJkqSZt+iVsCRXA2uAt3RL949TVQ8eba+q36F3r51hK1xeXSRJkmbeolbCkrwHeCnw/VX16AJ9zqmqB7v3rwWOAA8O9quqQ8Chgd89wbIlSZKWt+OGsCQvBDYDe4A/6wLTfVX1b5PsBjZU1UPAB7sH3j4G/A/gf62qf1q60iVJkpavxVwdeRcwdKmqql7S9/57R1iXJEnSVPOO+ZIkSQ0YwiRJkhowhEmSJDVgCJMkSWrAECZJktSAIUySJKkBQ5gkSVIDhjBJkqQGDGGSJEkNGMIkSZIaMIRJkiQ1YAiTJElqwBAmSZLUgCFMkiSpAUOYJElSA4YwSZKkBgxhkiRJDRjCJE21JGuT3JZkT/e65hh9n5/k75NcPc4aJc0mQ5ikaXctsK2q1gLbgOuGdUoy1/3s5jHWJmmGGcIkTa0kZwHrgO1d03ZgXZIzh3R/J/BJYM+YypM041a0LkCSltC5wINVdQSgqo4keahrP3i0U5ILgNcCrwbevdBgSVYBqwaa50ddtKTZYAiTNNOSPAm4HviJLqQdq/smYMtYCpM09QxhkqbZA8A5Sea6gDUHnN21H/Us4LnAji6ArQKS5BlV9VMD410D3DDQNg/sXIriJU03Q5ikqVVVB5LsBi4Fbuxed1XVwb4++4BvP/o5yVbg9Kr6hSHjHQIO9bcdZ+VMkhbkifmSpt3lwMYke4CN3WeS7EhyYdPKJM00V8IkTbWquhu4eEj7hgX6b13qmiQJXAmTJElq4rghLMkZ3bL9l5LckeSmYffYSfI/JflIki8nuTvJ65emZEmSpOVvMSthBbyvqp5fVRcA9wLvHdLvF4DDVfU84A3Af05y+uhKlSRJmh7HDWFV9XdVdUtf0+eB84d0fQu9x4NQVfcAfwG8bgQ1SpIkTZ0TOjE/yWnAFcAnhvz4POD+vs/76N2VenAM7zgtSZJm3oleHflrwCPA+09hn95xWpIkzbxFXx2Z5GpgDfCWqnpsSJd9PP5ryvN4/F2pj7oGePbAtn6xdUiSJE2DRa2EJXkP8FLg+6vq0QW6fRT4aeAvkqwBLqJ3d+rH8Y7TkiRJi7tFxQuBzfSet/ZnSXYn+YPuZ7uTnN11/RVgVZIvA58EfqqqDi9R3ZIkScvacVfCquouYOhSVVW9pO/9N4A3j640SZKk6eUd8yVJkhowhEmSJDXgA7wlaQQ2b4aVK1tXIWkcDo/ojHdXwiRJkhowhEmSJDXg15GSNAJXXQWrV7euQtI47N0L27ef+jiuhEmSJDVgCJMkSWrAECZJktSAIUySJKkBQ5gkSVIDhjBJkqQGDGGSJEkNeJ8wSVMtyVrgg8AZwNeBy6rqnoE+7wZ+BPinbttcVZ85oR353CJpdozouUWuhEmadtcC26pqLbANuG5Inz8HLqqq7wR+EvhIkqeNsUZJM8gQJmlqJTkLWAccvbf1dmBdkjP7+1XVZ6rq77uPdwCht3ImSUvGryMlTbNzgQer6ghAVR1J8lDXfnCB37kMuLeq9g/+IMkqYNVA8zzgc4ukWTKi5xYZwiSpk+RVwC8BlyzQZROwZXwVSZpmhjBJ0+wB4Jwkc90q2Bxwdtf+OEleBtwIvLGqvrTAeNcANwy0zQM7R1eypFlhCJM0tarqQJLdwKX0AtalwK6qetxXkUkuAj4C/FBV/eUxxjsEHBr43ZHXLWk2eGK+pGl3ObAxyR5gY/eZJDuSXNj1+XXgacB1SXZ324vblCtpVrgSJmmqVdXdwMVD2jf0vb9orEVJEq6ESZIkNWEIkyRJasAQJkmS1IAhTJIkqQFDmCRJUgOLCmFJrk5yX5JK8qIF+mxNcqDv8u5toy1VkiRpeiz2FhU3A/+J498V+neq6hdOrSRJkqTpt6gQVlW3gneGliRJGpVR36z1R5K8BvgasKWqbhvskGQVsGqgeX7EdUiSJE20UYawa4H3VNW3klwCfDzJC6rq6wP9NgFbRrhfSZKkZWdkV0dW1deq6lvd+/8HeAAYdhL/NcCzB7b1o6pDkiRpORjZSliSc6rqwe79S4DVwJcG+1XVIeDQwO+OqgxJkqRlYbG3qPjVJPvpnbv12SR3de07klzYdbsqyZ1JbgeuB368qr62JFVLkiQtc4u9OvIdwDuGtG/oe//vRliXJEnSVPOO+ZIkSQ0YwiRJkhowhEmSJDVgCJMkSWpg1HfMPyWbN8PKla2rkDQOhw+3rkCS2nIlTJIkqQFDmCRJUgMT9XXkVVfB6tWtq5A0Dnv3wvbtrauQpHZcCZMkSWrAECZpqiVZm+S2JHu61zVD+swl2Zbk3iRfTvK2FrVKmi2GMEnT7lpgW1WtBbYB1w3p81bgecAa4GXA1iSrx1WgpNlkCJM0tZKcBawDjp59th1Yl+TMga5vAa6vqseq6iBwM/Dm8VUqaRZNyon5cwD79+9vXYekMen79z63hLs5F3iwqo4AVNWRJA917Qf7+p0H3N/3eV/X53GSrAJWDTSfD85f0iwZ1fw1KSFsDcD69etb1yFp/NYA97YuYpE2AVuG/cD5S5pJpzR/TUoI+0r3+ip6f4EuZ/PATmA9sNz/NJ6WY5mW44DpOpbzgD/hX/79L4UHgHOSzHWrYHPA2V17v330VrS+0Ffb/TzRNcANA23PAf4Y569J47FMnmk5DhjR/DUpIeyb3eu+qtrbspBTleTo2/0ey2SYluOAqT2Wbx6r36moqgNJdgOXAjd2r7u68776fRR4e5KbgDOANwGvHDLeIeBQf1vfcTh/TRCPZfJMy3HA6OYvT8yXNO0uBzYm2QNs7D6TZEeSC7s+H6L3F+09wOeBK6tqKVfoJGliVsIkaUlU1d3AxUPaN/S9PwJcMc66JMmVMEmSpAYmJYQdAn6RgXMtlimPZfJMy3GAxzKJpuU4wGOZVNNyLNNyHDCiY0lVjaYcSZIkLdqkrIRJkiTNFEOYJElSA4YwSZKkBgxhkiRJDRjCJEmSGjCESZIkNWAIkyRJasAQJkmS1IAhTMtWkq1JtrauQ5JOlPOXwBCmBST5tiR/kOQbSe5P8qNLvL+nJPlAt6/DSXYled0pjvmzSf4iyaNJbhhRqZIm2FLMJYvc70jnTOev2bCidQGaWNuAbwLPBF4CfCrJ7VV11xLtbwXwAPAqYB+wAfj9JC+uqr0nOeZDwC8DrwWeNooiJU28pZhLFmPUc6bz1wxwJUxPkOTpwA8C766qR6rqVuATwI+f5Hg5Xp+q+kZVba2qvVX1WFV9ErgPeOnJ7LMb86aquhn4+smOIWl5GfVcspj5a9RzJjh/zQpDmIZZCxypqj19bbcDLzzRgZKsB3YkeeoJ/t4zuzqWauVN0gw4lbnkBOavkc2Zmi2GMA1zOvDwQNvDwMqTGOtzwAHgE4sNYkmeBHwY+GBV3X0S+5SkUcwli52/RjlnaoYYwjTMI8AzBtqeARxe6BeSfF+SGtyAI8BlwCXAFcfbcZLTgA/RO7fiZ0/2ACTNthOZS0Ywf53wnCmBJ+ZruD3AiiRrquqeru07OcZyflV9GnjCuRPdRPjbwLOAa4+10+7ciw/QO7F1Q1V96+TKlzTLTnQuGcH8dcJzpgSuhGmIqvoGcBNwZZKnJ3k58EZ6f1WeqJfTmwjfWFX/cJy+vwG8AHjDIvoeV5IV3VcIc8Bckqcm8Q8PafqNai5Z1Pw14jkTcP6aFYYwLeRn6F0WfQDYDlxxMpdaV9VO4HXHmwiTnA/8NL1Lu7+W5JFue+uJl/7P3gX8A/BO4Me69+86hfEkTbhRziWLnb86I5kz+zh/zYBUVesapJNy9G7TVbW1bSWSdGKcvwSuhEmSJDVhCNNydku3SQtKcnWS+7or3l60QJ+5JNuS3Jvky0neNu46NXNuwflr5vl1pKSpluQVwP3ATuD1VXXnkD6XAW8FXgecAewCXrHEj7mRNONcCZM01arq1qp64Djd3gJc3z3m5iBwM/Dmpa9O0iybiMtdkzwFuAj4Kr2b40mafnP07r/0hap6tHEt59FbLTtqH3DuYKckq4BVA81PBp4D3IPzlzQrRjJ/TUQIoxfAdrYuQlIT64FbWxexSJuALa2LkDQxTmn+mpQQ9lWAnTt3Mj8/37oWSWOwf/9+1q9fD92//8b2AecDX+g+D66MHXUNcMNA2/nALc5f0uwY1fw1KSHsCMD8/DyrV69uXIqkMZuEr/A+Crw9yU30Tsx/E/DKwU5VdQg41N/We0KO85c0o05p/vLEfElTLcmvJtkPzAOfTXJX174jyYVdtw8BX6F3XtfngSur6itNCpY0MyZlJUySlkRVvQN4x5D2DX3vjwBXjLMuSXIlTJIkqYGJWgnbvBlWrmxdhaRxOHy4dQWS1JYrYZIkSQ0YwiRJkhqYqK8jr7oKvMJbmg1798L27a2rkKR2XAmTJElqwBAmSZLUgCFMkiSpAUOYJElSA4YwSZKkBgxhkiRJDRjCJEmSGjCESZIkNWAIkyRJasAQJkmS1IAhTJIkqQFDmCRJUgOGMEmSpAYMYZIkSQ2MNIQleWqS30hyT5IvJvnNUY4vSZI0LVaMeLz3Af8IrK2qSvLMEY8vSZI0FUYWwpKcDlwGzFdVAVTV34xqfEmSpGkyypWw5wJfB7YkeTXwCPCuqrq1v1OSVcCqgd+dH2EdkiRJE2+U54StAJ4D7KqqC4H/A7gpyTMG+m0C7hvYdo6wDkn6Z0nWJrktyZ7udc2QPmcl+VSSO5LcneTXk4z6dA1JepxRhrD7gX8CtgNU1X8H/hZYO9DvGuDZA9v6EdYhSf2uBbZV1VpgG3DdkD6bgb+uqguAFwMvBX5gfCVKmkUj+0uvqv42yX8DLgH+KMla4CzgywP9DgGH+tuSjKoMSfpnSc4C1tGbl6D3R+L7k5xZVQf7uhawMslpwFOAJwMPjrVYSTNn1MvtlwO/leQ/AN8CfrwLXZLUwrnAg1V1BKCqjiR5qGvvD2G/BPwX4KvA04H3V9XnBgfznFZJozTSEFZVXwH+zSjHlKQxeDNwB/A9wErgD5P8UFV9bKDfJmDLuIuTNJ28Y76kafYAcE6SOYDu9eyuvd9G4MNV9VhVPQx8HHj1kPE8p1XSyBjCJE2tqjoA7AYu7ZoupXcF98GBrvcB3weQ5MnA9wJ3DhnvUFXt7d+A/UtVv6TpZgiTNO0uBzYm2UNvxetygCQ7klzY9dkErE/yRXqhbQ9wfYtiJc0O74MjaapV1d3AxUPaN/S9v5d/uYJSksbClTBJkqQGDGGSJEkNGMIkSZIaMIRJkiQ1YAiTJElqwBAmSZLUgCFMkiSpAUOYJElSA4YwSZKkBgxhkiRJDRjCJEmSGjCESZIkNWAIkyRJasAQJkmS1IAhTJIkqQFDmCRJUgOGMEmSpAYMYZIkSQ0YwiRJkhowhEmSJDVgCJMkSWrAECZJktSAIUySJKkBQ5ikqZZkbZLbkuzpXtcs0O+Hk3wxyZ3d6zPHXauk2WIIkzTtrgW2VdVaYBtw3WCHJBcCW4FLqupFwCuAh8dZpKTZYwiTNLWSnAWsA7Z3TduBdUnOHOj6c8DVVfU1gKp6uKr+cXyVSppFK1oXIElL6Fzgwao6AlBVR5I81LUf7Ov3HcB9Sf4UOB24CXhPVVX/YElWAasG9jG/VMVLmm6GMEnqzYUXAJcATwY+DewDfmeg3yZgy3hLkzSt/DpS0jR7ADgnyRxA93p2197vfuBjVfVoVR0GPg5815DxrgGePbCtX6LaJU05Q5ikqVVVB4DdwKVd06XArqo6OND1d4HXpOdJwPcAtw8Z71BV7e3fgP1LdwSSppkhTNK0uxzYmGQPsLH7TJId3VWRAL8HHAD+il5ouwv4QINaJc2QJTknLMkWepd7v7iq7lyKfUjSYlTV3cDFQ9o39L1/DPj5bpOksRj5SliSdcB30zupVZIkSUOMNIQleQq9myH+DFDH6S5JkjSzRv115JXAjVV1X5KhHbzPjiRJ0ghDWJKXARcB7zxOV++zI0mSZt4ov458FfCv6d11ei+91a3PJHnNQD/vsyNJkmbeyFbCquq9wHuPfu6C2OsHr46sqkPAof62hb66lCRJmlbeJ0ySJKmBJXt2ZFWtXqqxJUmSljtXwiRJkhowhEmSJDVgCJMkSWrAECZJktSAIUySJKkBQ5gkSVIDhjBJkqQGDGGSJEkNGMIkSZIaMIRJkiQ1YAiTJElqwBAmSZLUgCFMkiSpAUOYJElSA4YwSVMtydoktyXZ072uOUbf5yf5+yRXj7NGSbPJECZp2l0LbKuqtcA24LphnZLMdT+7eYy1SZphhjBJUyvJWcA6YHvXtB1Yl+TMId3fCXwS2DOm8iTNuBWtC5CkJXQu8GBVHQGoqiNJHuraDx7tlOQC4LXAq4F3LzRYklXAqoHm+VEXLWk2GMIkzbQkTwKuB36iC2nH6r4J2DKWwiRNPUOYpGn2AHBOkrkuYM0BZ3ftRz0LeC6wowtgq4AkeUZV/dTAeNcANwy0zQM7l6J4SdPNECZpalXVgSS7gUuBG7vXXVV1sK/PPuDbj35OshU4vap+Ych4h4BD/W3HWTmTpAV5Yr6kaXc5sDHJHmBj95kkO5Jc2LQySTPNlTBJU62q7gYuHtK+YYH+W5e6JkkCV8IkSZKaMIRJkiQ1YAiTJElqwBAmSZLUgCFMkiSpAUOYJElSA4YwSZKkBgxhkiRJDXizVkkahc2bYeXK1lVIGofDh0cyjCthkiRJDRjCJEmSGvDrSEkahauugtWrW1chaRz27oXt2095GFfCJEmSGhhpCEtyRpIdSb6U5I4kNyU5c5T7kCRJmgajXgkr4H1V9fyqugC4F3jviPchSZK07I00hFXV31XVLX1NnwfOH+U+JEmSpsGSnZif5DTgCuATA+2rgFUD3eeXqg5JkqRJtJRXR/4a8Ajw/oH2TcCWJdyvJEnSxFuSEJbkamAN8Iaqemzgx9cANwy0zQM7l6IWSZKkSTTyEJbkPcBLge+vqkcHf15Vh4BDA78z6jIkSZIm2khDWJIXApuBPcCfdeHqvqr6t6PcjyRJ0nI30hBWVXcBLmtJkiQdh3fMlyRJasBnR0qaaknWAh8EzgC+DlxWVfcM9Hk38CPAP3Xb5qr6zInsZ/NmWLlyNDVLmmyHD49mHFfCJE27a4FtVbUW2AZcN6TPnwMXVdV3Aj8JfCTJ08ZYo6QZZAiTNLWSnAWsA7Z3TduBdYPPtK2qz1TV33cf76B3busZYytU0kzy60hJ0+xc4MGqOgJQVUeSPNS1H1zgdy4D7q2q/YM/ONYTP666ClavHlXZkibZ3r2wfftxux2XIUySOkleBfwScMkCXXzih6SRMYRJmmYPAOckmetWweaAs7v2x0nyMuBG4I1V9aUFxvOJH5JGxhAmaWpV1YEku4FL6QWsS4FdVfW4ryKTXAR8BPihqvrLY4znEz8kjYwn5kuadpcDG5PsATZ2n0myI8mFXZ9fB54GXJdkd7e9uE25kmaFK2GSplpV3Q1cPKR9Q9/7iwjwPMwAABeDSURBVMZalCThSpgkSVIThjBJkqQGDGGSJEkNGMIkSZIaMIRJkiQ1YAiTJElqwBAmSZLUgCFMkiSpAUOYJElSA4YwSZKkBgxhkiRJDRjCJEmSGjCESZIkNWAIkyRJasAQJkmS1IAhTJIkqQFDmCRJUgOGMEmSpAYMYZIkSQ2saF3A42zeDCtXtq5C0jgcPty6AklqypUwSZKkBgxhkiRJDUzW15FXXQWrV7euQtI47N0L27e3rkKSmnElTNJUS7I2yW1J9nSva4b0mUuyLcm9Sb6c5G0tapU0WwxhkqbdtcC2qloLbAOuG9LnrcDzgDXAy4CtSVaPq0BJs8kQJmlqJTkLWAcc/d5zO7AuyZkDXd8CXF9Vj1XVQeBm4M3jq1TSLJqUc8LmAPbv39+6Dklj0vfvfW4Jd3Mu8GBVHQGoqiNJHuraD/b1Ow+4v+/zvq7P4yRZBawaaD4fnL+kWTKq+WtSQtgagPXr17euQ9L4rQHubV3EIm0Ctgz7gfOXNJNOaf6alBD2le71VfT+Al3O5oGdwHpguf9pPC3HMi3HAdN1LOcBf8K//PtfCg8A5ySZ61bB5oCzu/Z+++itaH2hr7b7eaJrgBsG2p4D/DHOX5PGY5k803IcMKL5a1JC2De7131VtbdlIacqydG3+z2WyTAtxwFTeyzfPFa/U1FVB5LsBi4Fbuxed3XnffX7KPD2JDcBZwBvAl45ZLxDwKH+tr7jcP6aIB7L5JmW44DRzV+emC9p2l0ObEyyB9jYfSbJjiQXdn0+RO8v2nuAzwNXVtVSrtBJ0sSshEnSkqiqu4GLh7Rv6Ht/BLhinHVJkithkiRJDUxKCDsE/CID51osUx7L5JmW4wCPZRJNy3GAxzKppuVYpuU4YETHkqoaTTmSJElatElZCZMkSZophjBJkqQGDGGSJEkNGMIkSZIaMIRJkiQ1YAiTJElqwBAmSZLUgCFMkiSpAUOYJElSA4YwLVtJtibZ2roOSTpRzl8CQ5gWkORnk/xFkkeT3DCG/T0lyQeS3J/kcJJdSV53imN+W5I/SPKNbtwfHVW9kibTUswli9zvSOdM56/ZsKJ1AZpYDwG/DLwWeNoY9rcCeAB4FbAP2AD8fpIXV9XekxxzG/BN4JnAS4BPJbm9qu4aQb2SJtNSzCWLMeo50/lrBrgSpqGq6qaquhn4+qmOlSSL2N83qmprVe2tqseq6pPAfcBLT3KfTwd+EHh3VT1SVbcCnwB+/GTGk7Q8LMFcctz5q9vvKOdM568ZYQjTkkqyHtiR5Kkn+HvPBNYCJ/tX31rgSFXt6Wu7HXjhSY4naRk6lbnkZOevEXD+mhGGMC21zwEHgE8sdiJL8iTgw8AHq+ruk9zv6cDDA20PAytPcjxJy8wI5pITnr9GxPlrRhjCNBJJvi9JDW7AEeAy4BLgikWMcxrwIXrnQvzsKZT0CPCMgbZnAIdPYUxJy8SJzCWjmr9GyPlrRnhivkaiqj4NPOHciW4i/G3gWcC1xxqjO/fiA/RORN1QVd86hZL2ACuSrKmqe7q27+Tkv96UtEyc6FwyivlrxJy/ZoQrYRoqyYpu+X0OmEvy1CQnE9pfTm8ifGNV/cNx+v4G8ALgDYvoe0xV9Q3gJuDKJE9P8nLgjfT+MpY03UY1lyx6/hrhnOn8NUMMYVrIu4B/AN4J/Fj3/l0nOkhV7QRet4gJ7Hzgp+ldiv21JI9021tPuPJ/8TP0LhU/AGwHrvDybmm6jXIuWez81RnJnNnH+WsGpKpa1yCdlKN3m66qrW0rkaQT4/wlcCVM0pRLcnWS+7qTrV+0QJ+5JNuS3Jvky0neNu46Jc0eT8zXcnZL6wK0LNwM/Cdg5zH6vBV4HrAGOAPYleSzS3yHdc22W1oXoPb8OlLSTEiyF3h9Vd055GefAn67qj7WfX4/cH9V/cp4q5Q0SyZiJSzJU4CLgK/Suy+LpOk3R+/S/y9U1aONazkPuL/v8z7g3MFOSVYBqwaanww8B7gH5y9pVoxk/pqIEEYvgB3rqwJJ02s9cGvrIhZpE7CldRGSJsYpzV+TEsK+CrBz507m5+db1yJpDPbv38/69euh+/ff2D7gfOAL3efBlbGjrgFuGGg7H7jF+UuaHaOavyYlhB0BmJ+fZ/Xq1Y1LkTRmk/AV3keBtye5id6J+W8CXjnYqaoOAYf623o3Z3f+kmbUKc1f3qJC0lRL8qtJ9gPzwGeT3NW170hyYdftQ8BX6J3X9Xngyqr6SpOCJc2MSVkJk6QlUVXvAN4xpH1D3/sjjPcBzZLkSpgkSVILhjBJkqQGJuvryM2bYeXK1lVIGofDh1tXIElNuRImSZLUwAmHsGEPw01yRnel0ZeS3JHkpiRnjr5cSZKk6XAyX0cOexhuAe+rqlsAkvwK8F7gfzuhka+6CrzPjjQb9u6F7dtbVyFJzZzwSlhV3VpVDwy0/d3RANb5PL27SEuSJGmIkZ+Yn+Q0evfb+cQCPx/2AFyf9SFJkmbKUlwd+WvAI8D7F/i5D8CVJEkzb6QhLMnVwBrgDVX12ALdhj0Ad57Hn2MmSZI01UYWwpK8B3gp8P1V9ehC/Y71AFxJkqRZcTK3qHjCw3CTvBDYDJwN/FmS3Un+YMS1SpIkTY0TXglb6GG4gMtZkiRJi+Qd8yVJkhowhEmSJDVgCJMkSWrAECZJktSAIUySJKkBQ5gkSVIDhjBJkqQGDGGSJEkNGMIkSZIaMIRJkiQ1YAiTNNWSrE1yW5I93euaIX3OSvKpJHckuTvJryc54ce6SdKJMIRJmnbXAtuqai2wDbhuSJ/NwF9X1QXAi4GXAj8wvhIlzSJDmKSpleQsYB2wvWvaDqxLcuZA1wJWJjkNeArwZODBsRUqaSa53C5pmp0LPFhVRwCq6kiSh7r2g339fgn4L8BXgacD76+qzw0OlmQVsGqgeX4pCpc0/VwJkyR4M3AH8CzgHOCVSX5oSL9NwH0D285xFSlpuhjCJE2zB4BzkswBdK9nd+39NgIfrqrHquph4OPAq4eMdw3w7IFt/RLVLmnKGcIkTa2qOgDsBi7tmi4FdlXVwYGu9wHfB5DkycD3AncOGe9QVe3t34D9S1W/pOlmCJM07S4HNibZQ2/F63KAJDuSXNj12QSsT/JFeqFtD3B9i2IlzQ5PzJc01arqbuDiIe0b+t7fC1wyzrokyZUwSZKkBgxhkiRJDRjCJEmSGjjhEJbk6iT3JakkL+prP+7z2SRJktRzMithNwOvBO4faF/M89kkSZLESYSwqrq1qh53o8MTeD6bJEmSGN0tKhb7fDafvSZJkkSb+4RtArY02K8kSdLEGFUI++fns3WrYAs9nw16z167YaBtHh+CK0mSZshIQlhVHUhy9PlsN7Lw89moqkPAof62JKMoQ5Ikadk4mVtU/GqS/fRWrz6b5K7uR0OfzyZJkqQnOuGVsKp6B/COIe1Dn88mSZKkJ/KO+ZIkSQ0YwiRJkhowhEmSJDVgCJMkSWrAECZJktSAIUySJKkBQ5gkSVIDhjBJkqQGDGGSJEkNGMIkSZIaMIRJkiQ1YAiTNNWSrE1yW5I93euaBfr9cJIvJrmze33muGuVNFsMYZKm3bXAtqpaC2wDrhvskORCYCtwSVW9CHgF8PA4i5Q0ewxhkqZWkrOAdcD2rmk7sC7JmQNdfw64uqq+BlBVD1fVP46vUkmzaEXrAiRpCZ0LPFhVRwCq6kiSh7r2g339vgO4L8mfAqcDNwHvqarqHyzJKmDVwD7ml6p4SdPNECZJvbnwAuAS4MnAp4F9wO8M9NsEbBlvaZKmlV9HSppmDwDnJJkD6F7P7tr73Q98rKoerarDwMeB7xoy3jXAswe29UtUu6QpZwiTNLWq6gCwG7i0a7oU2FVVBwe6/i7wmvQ8Cfge4PYh4x2qqr39G7B/6Y5A0jQzhEmadpcDG5PsATZ2n0myo7sqEuD3gAPAX9ELbXcBH2hQq6QZ4jlhkqZaVd0NXDykfUPf+8eAn+82SRoLV8IkSZIaMIRJkiQ1MNIQluT1SXYl2Z3kjiQ/MMrxJUmSpsXIzglLEuBDwPqqujPJBcDnktzcnW8hSZKkzqi/jnwM+Ffd+1XAVw1gkiRJTzSylbCqqiQ/DHw8yTeAlcD3D/bzsR+SJEmj/TpyBfB/Am+sqs8leTnwkSTfUVWP9HX1sR+SJGnmjfLryJcAZ1fV5wC6128ALxjo52M/JEnSzBvlzVr3A/NJnl9VX0ryAuB/Ae7t71RVh4BD/W29c/olSZJmxyjPCftakiuAjyU5ejL+T1TV341qH5IkSdNipI8tqqoPAx8e5ZiSJEnTyDvmS5IkNWAIkyRJasAQJkmS1IAhTJIkqQFDmCRJUgOGMEmSpAYMYZIkSQ0YwiRJkhowhEmSJDVgCJMkSWrAECZJktSAIUySJKkBQ5ikqZZkbZLbkuzpXtcco+/zk/x9kqvHWaOk2WQIkzTtrgW2VdVaYBtw3bBOSea6n908xtokzTBDmKSpleQsYB2wvWvaDqxLcuaQ7u8EPgnsGVN5kmbcitYFSNISOhd4sKqOAFTVkSQPde0Hj3ZKcgHwWuDVwLsXGizJKmDVQPP8qIuWNBsMYZJmWpInAdcDP9GFtGN13wRsGUthkqaeIUzSNHsAOCfJXBew5oCzu/ajngU8F9jRBbBVQJI8o6p+amC8a4AbBtrmgZ1LUbyk6WYIkzS1qupAkt3ApcCN3euuqjrY12cf8O1HPyfZCpxeVb8wZLxDwKH+tuOsnEnSgjwxX9K0uxzYmGQPsLH7TJIdSS5sWpmkmeZKmKSpVlV3AxcPad+wQP+tS12TJMGIV8KSPDXJbyS5J8kXk/zmKMeXJEmaFqNeCXsf8I/A2qqqJM8c8fiSJElTYWQhLMnpwGXAfFUVQFX9zajGlyRJmiajXAl7LvB1YEuSVwOPAO+qqlv7O3mzQ0mSpNGGsBXAc+hd/v2/J7kY+L+TPK+q/kdfP292KEmSZt4oT8y/H/gnume0VdV/B/4WWDvQ7xrg2QPb+hHWIUmSNPFGthJWVX+b5L8BlwB/lGQtcBbw5YF+3uxQkiTNvFFfHXk58FtJ/gPwLeDHu9AlSZKkPiMNYVX1FeDfjHJMSZKkaeRjiyRJkhowhEmSJDVgCJMkSWrAECZJktSAIUySJKkBQ5gkSVIDhjBJkqQGDGGSJEkNGMIkSZIaMIRJkiQ1YAiTJElqwBAmSZLUgCFMkiSpAUOYJElSA4YwSVMtydoktyXZ072uGdLn3UnuSnJ7kv83yWtb1CppthjCJE27a4FtVbUW2AZcN6TPnwMXVdV3Aj8JfCTJ08ZYo6QZZAiTNLWSnAWsA7Z3TduBdUnO7O9XVZ+pqr/vPt4BBDhjbIVKmkkrWhcgSUvoXODBqjoCUFVHkjzUtR9c4HcuA+6tqv2DP0iyClg10Dw/wnolzRBDmCR1krwK+CXgkgW6bAK2jK8iSdPMECZpmj0AnJNkrlsFmwPO7tofJ8nLgBuBN1bVlxYY7xrghoG2eWDn6EqWNCsMYZKmVlUdSLIbuJRewLoU2FVVj/sqMslFwEeAH6qqvzzGeIeAQwO/O/K6Jc0GT8yXNO0uBzYm2QNs7D6TZEeSC7s+vw48Dbguye5ue3GbciXNiiVZCUuyBdgKvLiq7lyKfUjSYlTV3cDFQ9o39L2/aKxFSRJLsBKWZB3w3cC+UY8tSZI0LUYawpI8hd7NEH8GqFGOLUmSNE1G/XXklcCNVXXfQierep8dSZKkEYaw7vLui4B3Hqer99mRJEkzb5RfR74K+NfAfUn20lvd+kyS1wz0uwZ49sC2foR1SJIkTbyRrYRV1XuB9x793AWx1w9eHel9diRJkrxPmCRJUhNLdsf8qlq9VGNLkiQtd66ESZIkNWAIkyRJasAQJkmS1IAhTJIkqQFDmCRJUgOGMEmSpAYMYZIkSQ0YwiRJkhowhEmSJDVgCJMkSWrAECZJktTAkj078mRs3gwrV7auQtI4HD7cugJJasuVMEmSpAYMYZIkSQ1M1NeRV10Fq1e3rkLSOOzdC9u3t65CktpxJUzSVEuyNsltSfZ0r2uG9JlLsi3JvUm+nORtLWqVNFsMYZKm3bXAtqpaC2wDrhvS563A84A1wMuArUlWj6tASbPJECZpaiU5C1gHHP3iczuwLsmZA13fAlxfVY9V1UHgZuDN46tU0iyalHPC5gD279/fug5JY9L3731uCXdzLvBgVR0BqKojSR7q2g/29TsPuL/v876uz+MkWQWsGmg+H5y/pFkyqvlrUkLYGoD169e3rkPS+K0B7m1dxCJtArYM+4HzlzSTTmn+mpQQ9pXu9VX0/gJdzuaBncB6YLn/aTwtxzItxwHTdSznAX/Cv/z7XwoPAOckmetWweaAs7v2fvvorWh9oa+2+3mia4AbBtqeA/wxzl+TxmOZPNNyHDCi+WtSQtg3u9d9VbW3ZSGnKsnRt/s9lskwLccBU3ss3zxWv1NRVQeS7AYuBW7sXnd15331+yjw9iQ3AWcAbwJeOWS8Q8Ch/ra+43D+miAey+SZluOA0c1fnpgvadpdDmxMsgfY2H0myY4kF3Z9PkTvL9p7gM8DV1bVUq7QSdLErIRJ0pKoqruBi4e0b+h7fwS4Ypx1SZIrYZIkSQ1MSgg7BPwiA+daLFMey+SZluMAj2USTctxgMcyqablWKblOGBEx5KqGk05kiRJWrRJWQmTJEmaKYYwSZKkBpqHsCRrk9yWZE/3uqZ1TYuV5IzuMvcvJbkjyU1Hn0mX5LuT3N4d1x91z7CbeEm2JKkkL+o+L7vjSPLUJL+R5J4kX0zym137svtvLcnrk+xKsrv7b+wHuvaJPpYkVye5r/+/pa59wbon/ZgWsozrdv6aQM5fk2Fsc1hVNd2A/wr8WPf+x4D/2rqmE6j924B/0/f5V4APAAG+DLyia38X8Fut613E8awD/pDencJftIyP41eB/8i/nPP4zO51Wf231v3v//8BL+o+XwAcpvfH00QfC/AKes9e3Hu0/uP9fzDpx3SMY12udTt/TeDm/DUZ27jmsNYHeRa9Kwvmus9z3eczW/8fcJLH84PAZ4GLgDv72r8deKR1fcep/SnAbcCzj/5Ht0yP4/Tuv6HTB9qX3X9r3ST2deDl3edXAnuW07H0T2DHqns5HdNy/+/qGMfi/NX+OJy/Jmxb6jms9deR5wIPVu9GiXSvD3Xty0qS0+jd7PETDDx3rqr+Fjgtybc1Km8xrgRurKr7+tqW43E89/9v795ZswiiAAy/B9PYW0gUkkJJxEbQRhsrEQsxlRIs/BliZSVY2Cn+BdHCSykWprEIeIkEC228pRAVK0UF8VjMJGjzxUB0dpL3geWD3eacj9nDmd1ZhnLjn4+IhxExFxHLM5quxlqWO/kkcCciXgO3gTN0mEs1Ku6NmFM3rF+DYf0atnWvYa2bsI3kMvAZuNI6kLWKiIOUWePV1rGsgzHKhspPMvMAcBa4SZlhdiUixoBzwInMnACOA9fpMBcNnvVrGKxfm0zrJuwtsCMitgDU3/F6vhsRcQnYDZzKzJ/AG2Dit+vbKBODT41CXM1hYBp4GRGvKDvd3wV20VceUGa+P4BrAJk5D3wEvtLfWNsHjGfmA4D6+wX4Rn+5wOj7vdda0GvcK6xfg2L9GrZ1r2FNm7DMfA8sALP11CxlBvChXVRrExEXgP3ATGZ+r6cfAVvrY2QoGwbfaBHf38jMi5k5npmTmTkJLAFHKQt1u8kDVl453AeOQPlahfKu/gX9jbUlYGdETAFExB5gO2WT6d5yGXm/91oLeo17mfVrWKxfw/ZPatgAFr1NA/OUQTYPTLWOaQ2x7wUSeF7//AXgVr12CFikDLh71C9cejj4cyFid3lQHufP1bgfA8d6HWvA6ZrH03rM9JAL5QuvJcqs/h3wbLW4h57TiFx7jdv6NcDD+jWM43/VMLctkiRJaqD1mjBJkqRNySZMkiSpAZswSZKkBmzCJEmSGrAJkyRJasAmTJIkqQGbMEmSpAZswiRJkhr4BXfwy3leyAjSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x720 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(nrows=3, ncols=2, sharex=True, figsize=(10, 10))\n",
    "\n",
    "for i, (x, y) in enumerate([(0,1), (1,0), (0,2), (2,0), (1,2), (2,1)]):\n",
    "        x_range = np.arange(len(granger1[x, y, :]))\n",
    "        ax.flatten()[i].plot(cond_granger_mean1[x, y, :], color='b', lw=3, alpha=0.6)\n",
    "        ax.flatten()[i].plot(cond_granger_mean2[x, y, :], color='r', lw=3, alpha=0.6)\n",
    "        ax.flatten()[i].set_title(r\"{} $\\rightarrow$ {} $\\vert$ {}\".format(x, y, 3-x-y))\n",
    "        ax.flatten()[i].set_xlim((0, 100))\n",
    "        if i % 2 == 1:\n",
    "            ax.flatten()[i].set_ylim((0, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember, we would like $F_{0 \\rightarrow 2 \\vert 1}^{seq} << F_{0 \\rightarrow 2 \\vert 1}^{diff}$ and $F_{1 \\rightarrow 2\\vert 0}^{diff} << F_{1 \\rightarrow 2\\vert 0}^{seq}$. Sequential drive is blue and differential drive is red. With our abbreviated notation, the two plots are located in panels $(2, 1)$ and $(3, 1)$. While we do not have, $F_{0 \\rightarrow 2 \\vert 1} = F_{1 \\rightarrow 2\\vert 0} = 0$, the cases of differentially delayed drive and sequential drive and be distinguished easily. The difference to the paper by Chen et al. (who actually achieve conditional causality close to zero) may have several reasons, such as the repeated computation over 500 samples or a different approach for computing the spectral density of the two processes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "normed_cond_granger1 = cond_granger_mean1 / np.nanmax([cond_granger_mean1, cond_granger_mean2])\n",
    "normed_cond_granger2 = cond_granger_mean2 / np.nanmax([cond_granger_mean1, cond_granger_mean2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 4)"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARIAAAELCAYAAAAV7glwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deZxU5Z3v8c+3G7OBogKaALJEdgm0u3jDEF43OJG4MJqoiZNAEhPNdXnlZphrxpVxybjlopkgJjEOkQQXRgWNTBJyR0YiLqiQaKKISSObKCAoOi5I/+4f53R7KKq7q/tU09XN9/161avrnOepU79TfepXz9meRxGBmVkeVe0dgJl1fE4kZpabE4mZ5eZEYma5OZGYWW5OJGaWmxNJAUm3SLq0veNoCUkhaVAbLHeKpN+Xe7klvvciSWelz8+U9Nvd9L4D0s+zS4n1Z0m6qonysZJWlC/CyrRHJRJJqyS9LWmbpK2Slkg6R1LD5xAR50TElbshjs+25Xt0JhHxy4g4rpS6kqZJ+kVbx1SqiFgcEUPbO462tkclktSJEbE30B+4BrgQ+Fn7hmSdUamtms5gT0wkAETE6xFxP3A6MFnSSNi5qSqpp6Rfpa2X1yQtrm+9SOot6R5JGyXVSrqgftnpr+Ldkm5PWz9/knREWjYb6Ac8IOlNSf8nnT9X0gZJr0t6WNIhmeXNkjRD0oPp8h6XdHCx9ZLUQ9IDkt6QtFTSVdndE0nDJC1M12eFpNMKXnt/+tongKLvkdat3wX4mqQ1krakrbsjJf0x/cx+VPCar0t6Lq37G0n9M2UTJD2frv+PAGXKphSsw03pe74h6SlJY9P5nwMuAk5PP9s/pPO7S/qZpJclrUs/k+q0rFrSDZI2Sfor8PnG1jmtf6ikp9P/w13ARzJln5G0VtKFkjYA/1Y/Ly3/nqR/L1jeTZJ+2FycFS8i9pgHsAr4bJH5q4Fvp89nAVelz/8FuAXYK32MJdnAq4CngMuADwGfBP4K/G36umnAO8BEoDpdzmNNxQF8Hdgb+DBwI7A8UzYLeA04CugC/BK4M1MewKD0+Z3p42PACGAN8Pu0rGs6/bV0OYcBm4BDMq+9O603ElhX/9oin9mA9H1vIfkyHZeu8zzgAKAP8CowLq0/CXgRGJ6+9yXAkrSsJ/AG8IX0c/7fwPvAWWn5lGwcwN8DPdLl/AOwAfhI5rP/RUGs84Afp+t1APAEcHZadg7wPHAQsD/wULpeXYqs84eAl9L49krj3c4H28tn0rivTf+PH03nrU3L+wP/DeyTTlcDLwPHNBdnpT/aPYDdurKNJ5LHgIvT57MyG8YVwPz6L2mm/tHA6oJ5/wT8W2Zj/l2mbATwdnNxZMr3TTfm7pmYbs2UTwSez0wHMCjdMLcDQzNlV/FBIjkdWFzwXj8GLs+8dlim7Ps0n0j6ZOZtBk7PTN8DfCd9/h/ANzJlVemXqj/wVXZOtALW0kgiKRLLFmB05rP/RabsQOBd4KOZeV8CHkqf/ydwTqbsOBpPJH8DrAeUmbeEnRPJe6RJLTNvbWb698BX0+cTgL+UEmelP/aYfbhm9CH5xS90PcmG+VtJAD+JiGtINv7ekrZm6lYDizPTGzLP/xv4iKQuEfF+4ZukzdergS8CvYC6tKgn8Hojy+tWJN5eJL/SazLzss/7A0cXxN0FmN3Ia18q8h6FXsk8f7vIdH2c/YGbJP0gUy6Sz7539n0jIiRl49iJpH8AzkpfF8A+JJ9VMf1JWg8vp/9DSJJY/fJ3em+aXufewLpIv+WN1N8YEe80sYw5JAniduDL6XQpcVa0PT6RSDqSZGPe5TRnRGwjaTr/Q3rM4iFJS0n+ubURMbiVb1t4y/WXgZOBz5K0VrqT/MqKltlI0rTuC7yQzjsoU74G+K+ImFD4wjSZvZ/Wfz6d3a+F79+UNcDVEfHLIu89OBunkm/SQYX10rKxJAfI/yfwp4iok5T9rAo/2zUkv/Q9iyVxkl2L7Hs1tc4vA30kKZNM+gF/ydRp7nb6ucAPJPUF/g4YU2KcFW2PPdgqaR9JJ5AcF/hFRDxTpM4JkgalG/YbwI708QTwRnpQ7aPpAbuRaVIqxSskx1Xq7U2yEW0mObbx/dasU0TsAO4Fpkn6mKRhJLsN9X4FDJH0FUl7pY8jJQ0v8toRwOTWxNGIW4B/qj+InB5Y/GJa9iBwiKRTlJzpuAD4eCPL2Zsk4W0Euki6jKRFUu8VYIDSg+IR8TLwW5Iv7z6SqiQdLGlcWv9u4AJJfSXtB3yviXV4NH3vCyR1kXQKyXGrkkXERmAR8G8kP0bPlRhnRdsTE8kDkraR/AJcDPxfkoOPxQwGfge8SbIR3RwRi9Iv3YlADVBLcsDyVpKWRCn+BbgkPbMxlaSZ+xLJwc0/kxyzaa3z0jg2kOyy3EGSpOpbWMcBZ5Ds62/ggwOD9a/tls6fRbKxl0VE3Je+152S3gCeBY5PyzaR7NZdQ5JMBwOPNLKo35Acb3mB5DN7h52b/3PTv5slPZ0+/yrJgdI/k7T0/h34RFr203SZfwCeJkmmja3De8ApJMdstpAcc2q0fhPmkLQ+5xTMbyrOiqadd/ess5F0LfDxiChn68JsJ3tii6RTU3KdyCgljgK+AdzX3nFZ59aiRCLpciUXIY0sUvYxSXdJejG9sOiE8oVpLbA3SXP7LZL9/x+QnMI2azMln7WRdBhwDMnFW8VMBbZFxKD0KPxiSYMi4s0yxGklioilJNeUmO02JbVIJH0YmAH8Lxo/vXU6yZF5ImIl8CTpwTQz69xKbZFcQXKKtDZzsUyhfux8cc5qilwLIGlfkis3s+ovM19JcnrVzHaPapIzQ0sj4t3WLqTZRCJpDHAkTZ9fb4nvkFySbWaVYyxFLsosVSktknHAMKC+NdIX+I2kr0VEtrOZ1SSX+W5Mp/uR3ABV6EaSaxSy+gOLFi9eTN++fUuP3sxyWbt2LWPHjoXkqt1WazaRpPeWXFM/LWkVcEJEPFtQdS5wNvBkerD1SJJ7CgqXtxXI3utB/e5S3759GTBgQItWwMzKItchhVzXkUhaLql3Onk9sK+kF0kuxf5WeiWlmXVyLb5pLyIGZJ7XZJ6/RXKZs5ntYfb4u39LUVdXx9q1a3nrrbfaOxSzVuvatSt9+/alqqr8F7Q7kZRg06ZNSGLo0KFt8k8wa2t1dXWsW7eOTZs2ccABB5R9+f5WlGDr1q0ceOCBTiLWYVVVVXHggQfy+uuvN1+5Nctvk6V2Mjt27GCvvfZq7zDMctlrr714//226TPJiaRETVzRa9YhtOU27ERiZrk5kZhZbk4kHdSAAQMYNmwYo0ePZtCgQZx88sksWbKkobympoa3334bgHnz5jF8+HAOPfRQVqxYscv07jBt2jTee++9ovE1RRJvvtl8TxTN1Sv1/Ur1wgsvMGbMGIYMGcKYMWNYuXJlq5e1efNmJk6cyNChQxk1ahSnnHIKGzdubCifMmUKixYtavT1U6dOZeDAgUji2WcLLzjfTdp7PIzIjJFSW1sblejPf/5ze4ewi/79+8czzzzTMH3PPfdE9+7d47HHHtul7uc+97m4++67G50u1fbt21sXbEQAsW3btjZ7XWP18sTclPHjx8fs2bMjImL27Nkxfvz4Vi9r8+bN8dBDDzVMT506Nb7+9a83TE+ePHmn8kKLFy+O1atX77JNFFO4LdfW1gZJ1yADwuPa7F5nn9327/HjH7es/imnnMITTzzBDTfcwNy5c5HEtm3buPTSS1m8eDErVqzg5ptvpqamZqfphx56iMcff5zvfe97vPHGGwBcccUVfP7zyciVkrjuuut48MEHGTt2LCeccEKTda+++mruu+8+Nm/ezPXXX8+pp57KueeeC8Cxxx5LVVUVixYtYr/99mPbtm1069aNM888kxUrVvDuu+8yaNAgbrvtNvbbb78m1/fee+/loosuYv/992fixIk7lRXGfOWVVzZ8HtOnT+e1115j+vTpQNIaGDJkCKtXr6Zr165Nfhb1Xn31VZ5++mkWLlwIwJe+9CXOO+88Nm7cSK9evRqNOSKYOnUql1xyyU7rt//++/OZz3ymYfqYY45h5syZTa5/1qc//emS67YVJ5JO5Oijj+b+++/fad706dNZtmwZU6dO5YQTkt4vs9Nbt27lnHPOYcGCBXziE5/g5Zdf5sgjj+TZZ59l332TbmPq6upYtGgRW7duZfz48U3W3WeffVi6dCmPPPIIp512GqeeeiozZszg5ptvZsmSJXTrtuu4XjfddBM9eybjW11yySVce+21XHPNNbvUq/fqq6/yzW9+kyVLljB06FCuu+66XerUx1xo8uTJHH300Vx//fV06dKFOXPmcPLJJ9O1a9eSPguANWvW0KdPH6qrk2F5q6ur6d27N2vWrGkykUji4IMPZsKECSxcuLBosqyrq2PmzJmcdNJJjS6nEjmRdCLRihEBlixZQm1tLccf/0FndpJ48cUXOeKII4Dky1dq3TPOOANIflXXr1/PO++8w0c+0jDOdlG33347v/zlL3nvvfd46623GDJkSJP1H3vsMQ477DCGDh0KwLe+9S0uvPDCnerUx1yoX79+jBgxggULFnDSSScxa9YsbrzxxpLXr1R33nkn5513XtGyLVu2cOGFF/KTn/xkl7Lzzz+fbt26NfraSuVE0got3e3YXZYuXcrIkbv0y92kiGDUqFE8/PDDjdapb0WUUrc+adT/Wjd3AdTixYuZOXMmS5YsoVevXsyZM6foF6ww5uYUa/nUmzJlCj//+c/55Cc/yeuvv17fH0dJ6wdw0EEHsW7dOnbs2EF1dTU7duxg/fr1HHTQBx0CnnHGGQ1JNeuuu+7iiiuuYNq0abuUTZ06lZUrV/LAAw90uKuoO1a01qj58+czc+ZMvvvd77bodcceeywrV67koYc+6INq6dKlRb+sLalbaO+99y56efbWrVvp3r07PXr04N133+W2225rdlljxoxh2bJlDWdKbr311mZfk3Xqqafy8MMPc8MNNzBlypSG+aWu3wEHHEBNTQ133HEHAHfccQeHHnpok7s1kCSqX//61yxcuJDevXvvVHbxxRfz1FNPMW/ePD784Q83soQKludIbbke+KxNi/Xv3z+GDh0ao0aNioMPPjhOPPHEeOSRRxrKyZzFGDduXDzwwAMNZYXTTzzxRIwbNy5GjRoVw4YNi4kTJ8aOHTt2WU5L62anp02bFkOGDInRo0fHli1bGsq2b98ep512WgwePDjGjx8f//iP/xjjxo0ruoyse+65J4YOHRpjxoyJ6dOn71Sv2GsK533jG9+IqqqqeOmll3aq19T6ZT333HNx1FFHxeDBg+Ooo46K559/fpc6pXr22WcDaPh8Ro8eHZMmTWoob+6szfnnnx99+vSJ6urqOPDAA2PEiBGN1m2rszYVMdKepAFAbW1tbUX2kPbcc88xfPjw9g7D9lBTpkxhypQpO53Zaa3CbXnVqlUMHDgQYGBErGrtcr1rY2a5OZGYVbhJkyZVZEs9q6SzNpLmAQOBOuBN4PyIWF5QZxrJAFrr01mPRMS55QvVbM80adKk9g6hWaWe/p0cEa8DSDoZuA04rEi92yNiarmCM7OOoaRdm/okkupO0jIxMwNaNoj4rcBxgIDPNVLtDEnHARuAyyPi0SLLKTZkp0fFMuvASk4kEXEWgKSvkIxhM7Ggyi3A1RGxXdIEYL6k4RGxuaCeh+w062RafNYmImYD4yX1KJi/ISK2p88XAmuAYtdr30hy4Db7GNvSOMyscpQyiHg3YL+IWJNOnwi8lj6y9fpExLr0eQ3J1aq79JoTTQzZaWYdUym7Nl2BuZK6kowP+hpwYkSEpAXAZRHxJPB9SYendd4DvhIRG9oqcDOrHKUMIv4KcEwjZRMzz4vft21mnZ6vbO2g5s6dy6GHHkpNTQ3Dhg3jy1/+cnuHtIu27qe1ubrl7qe1nH2j5u2ntZx9xpZFnjv+yvXAd/+2yPr166Nnz56xevXqiIioq6uLZcuWtXNUu6KN+2ltqm5b9NXakr5Rm5O3n9bW9hnrPlsrSTt32rphwwb22msvevRITpxJoqamBqDJPkcL+zm99NJL2bZtG5s2beKII45g06ZNQHJHaHa6uT5dd1c/rcXWIauwr9arrrqqbP20Quv6Ro026Ke1tX3GtiXv2nRAo0eP5qijjqJfv3584Qtf4MYbb2Tz5s0NfY7OmTOHp556il/96lecffbZbN26taGf0/nz57NkyRI+9KEPlfReTS2zXn0/rbNnz+aCCy4AYMaMGUDSfeHy5ct36vMUkn5an3zySZ555hkOOeQQrr322mZjKWUd6vtqvfLKKxvmTZ48mTvvvLOht7Zi/bQ2tX55ZPtp3bJlS9E6Le2ntak+Y9uLE0kHVFVVxbx581i0aBHjx4/nwQcfZNSoUSxYsKChz9GamhqOP/74hj5Hi/VzWopsP6aFy6xXrJ/W5tx+++0cfvjhfOpTn2LOnDksX7682deUsg7F+mrN9tMKMGvWLL72ta+VvH6luvPOO+nZs+cuj8suu4xly5bt0q9svY7aT2uWd21ao0I6bR05ciQjR47k3HPPZcSIEU32OTp//vxGl9OlSxfq6j64fSqbCJpaZr3d0U9rfSzNaayv1rz9tJZid/XTWkqfsbubWyQd0Lp163j00Q9uY1q7di0bN25kxIgRjfY52lQ/px//+MfZvn17w6/wnDlzGsoqpZ9WyNdXa95+Wlsr2qCf1tb2Gdum8hypLdcDn7VpkVWrVsWECRMa+vj81Kc+FbfccktENN3naFP9nP7sZz+LAQMGxLhx42LatGnRo0ePhverlH5am1uHpmKJyN9Pa0v6Rm1O3n5aW9tnrPtsbUedtc/W+tHnmhq6wdpHOftpzXKfrWZWsXywdQ9WCa1RK64j9NOa5URiVoE6Qj+tWd61MbPcnEhK5N0A6+jacht2IilBdXU127dvb+8wzHLZvn07Xbq0zdEMJ5IS7Lvvvrzyyis7Xf1p1pHU1dXxyiuv0L179zZZvg+2lqBnz56sXbuWFSt26TnSrMPo2rUrPXv2bJNlO5GUoKqqin79+rV3GGYVy7s2ZpZbSYlE0jxJf5C0TNLitJf4wjrVkmZI+oukFyWdVf5wzawSlXPs3zOBQcBgoAewTNLv8ly/b2YdQ0mJJEob+/d04KcRUQdslDQP+CLJqHwNPGSnWedTzrF/+wEvZaZXA8V6WvGQnWadTMkHWyPirIjoB1xEQSujhTxkp1kn0+LTvxExW9JPJPWInQcIXw30B5am04UtlPrXe8hOs06m2RaJpG6SDspMFx37F5gLfFNSlaRewCTgnnIGa2aVqZxj/84Gjgbqh/y6IiL+2hZBm1llKefYvzuAb5cvNDPrKHxlq5nl5kRiZrk5kZhZbk4kZpabE4mZ5eZEYma5OZGYWW5OJGaWmxOJmeXmRGJmuTmRmFluTiRmlpsTiZnl5kRiZrk5kZhZbk4kZpabE4mZ5eZEYma5ldL5cw9JCyStkPRHSfemnTsX1pslaa2k5enj4rYJ2cwqTSktkgCui4ihETEK+AtwTSN1r4mImvRxddmiNLOKVkrnz68BizKzHiNHJ88estOs82nRMRJJVSRJ5P5GqnxX0jOS5kka3kid7wC1BY/FLYnDzCpLS0fa+1fgTeBHRcouBl6OiDpJXwV+LemT6TAVWTcCswrm9cXJxKzDaskg4jcAg0kGx6orLI+IdZnnt0uaTpIgXiqo5yE7zTqZknZtJF0NHA5Mioh3G6nTJ/P8b0lG5VtXrK6ZdS7NtkgkHQJcBLwALElbD7UR8XeSlgMTI2I98HNJBwJ1wBvASRHxftuFbmaVopSzNn8Ciu57RERN5vlnyxiXmXUgvrLVzHJzIjGz3JxIzCw3JxIzy82JxMxycyIxs9ycSMwsNycSM8vNicTMcnMiMbPcnEjMLDcnEjPLzYnEzHJzIjGz3JxIzCw3JxIzy82JxMxycyIxs9zKOWTnxyTdJelFSc9LOqFtQjazSlPOITunAtsiYhBwInCrpG7lC9XMKlWziSQiXouIRZlZjwH9i1Q9Hbglfc1K4Eng+DLEaGYVrkUj7TUzZGc/dh4MazVwUJFleOxfs06mnEN2luo7wOU5Xm9mFabkszaZITtPLzZkJ0kLJLvL0w9YU6TejcDAgsfYUuMws8pTUoskM2Tn5xsbshOYC5wNPClpMHAk8KXCSh7716zzKeX0b/2Qnb1JhuxcLum+tGy5pN5p1euBfSW9CPwK+FZEbGujuM2sgpRzyM63gC+WLzQz6yh8ZauZ5eZEYma5OZGYWW5OJGaWmxOJmeXmRGJmuTmRmFluTiRmlpsTiZnl5kRiZrk5kZhZbk4kZpabE4mZ5eZEYma5OZGYWW5OJGaWmxOJmeXmRGJmuZWUSCTdIKlWUkga2UidaZJeTftxXS5pRnlDNbNKVeq4NvOAm4DFzdS7PSKm5gvJzDqakhJJRPwePGyEmRXX0pH2mnOGpOOADcDlEfFoYQUP2WnW+ZQzkdwCXB0R2yVNAOZLGh4RmwvqechOs06mbGdtImJDRGxPny8kGa6z2IFZD9lp1smUrUUiqU9ErEuf1wADgBWF9Txkp1nnU+rp3x9KWktyLON3kv6Uzl8g6Yi02vclPSvpD8BPga9ExIY2idrMKkqpZ20uAC4oMn9i5vnkMsZlZh2Ir2w1s9ycSMwsNycSM8vNicTMcnMiMbPcnEjMLDcnEjPLzYnEzHJzIjGz3JxIzCw3JxIzy82JxMxycyIxs9ycSMwsNycSM8vNicTMcnMiMbPcnEjMLLdmE0mJw3VWS5oh6S+SXpR0VvlDNbNKVUqLZB7wN8BLTdQ5ExgEDAbGANMkDcgbnJl1DM0mkoj4fUSsaaba6cBPI6IuIjaSJJ8vliNAM6t85RrXph87t1hWAwcVq+ghO806n3KP/VsKD9lp1smU66zNaqB/ZrofyZCdxXjITrNOplwtkrnANyXdC/QAJpEcoN2Fh+w063xKOf1bynCds4G/AiuBx4ArIuKvbRSzmVWYZlskJQ7XuQP4dnlDM7OOwle2mlluTiRmlpsTiZnl5kRiZrk5kZhZbk4kZpabE4mZ5eZEYma5OZGYWW5OJGaWmxOJmeXmRGJmuTmRmFluTiRmlpsTiZnl5kRiZrk5kZhZbk4kZpZbSYlE0hBJj0p6If07uEidaZJelbQ8fcwof7hmVolKbZHcAsyIiCHADODHjdS7PSJq0se5ZYnQzCpeKb3IHwAcBtyRzroDOExSr7YMzMw6jlJaJAcB69Ke4ut7jF9P8SE5z5D0R0m/lTSm2MIk7StpQPaBh+w069DKOWTnLcDVEbFd0gRgvqThEbG5oJ6H7DTrZEppkawB+kiqBkj/9qZgSM6I2BAR29PnC9PykUWW5yE7zTqZUgbIelXScuBLwC/Sv8siYmO2nqQ+EbEufV4DDABWFFmeh+w062RK3bU5B/i5pMuALcBXIRm2E7gsIp4Evi/pcGAH8B7wlYjY0AYxm1mFKSmRRMTzwNFF5meH7ZxcxrjMrAPxla1mlpsTiZnl5kRiZrk5kZhZbk4kZpabE4mZ5eZEYma5OZGYWW7lvGkvt4sugr33bu8ozPYc27aVZzlukZhZbk4kZpabIqK9YyDt3Ki2traWAQMGtG8wZnuQVatWMXDgQICBEbGqtctxi8TMcnMiMbPcnEjMLDcnEjPLzYnEzHJzIjGz3JxIzCy3ki6RlzQE+DnQA9gMfDUiVhbUqQZ+CHwOCOCaiLi1RdH4Gnmz3atM18iXc+zfM4FBwGBgDDAtvdDMzDq5ZlskmbF/J6Sz7gB+JKlXwdg2pwM/jYg6YKOkecAXgesLlrcvsG/B2/QHWPvmm61aCTNrncx3rjrPckrZtdll7F9J9WP/ZhNJP+ClzPRqio8P3OiQnWMfeKCUmM2s/AYDf2nti9ujG4EbgVkF8z4J/D9gHEkC6gj6AotJhhtd286xlKKjxQuOeXfoB/wX8Nc8CyklkTSM/Zu2RoqO/UuSAPoDSzMBvlRQp7khO1fnuXFod8rEvLYjxNzR4gXHvDtk4n0vz3KaPdgaEa8C9WP/QiNj/wJzgW9KqpLUC5gE3JMnODPrGEo9a3MOcL6kF4Dz02kkLZB0RFpnNknzaCXwGHBFRORqLplZx1DOsX93AN8uX2hm1lFUypWtW4F/puDYSYXraDF3tHjBMe8OZYm3InpIM7OOrVJaJGbWgTmRmFlu7Z5IJA2R9KikF9K/g9s7pkKSeqRnqFZI+qOke9NT3Eg6RtIf0vh/m95SUDEkXS4pJI1Mpys2XkkfkTRT0kpJz0j6STq/IrcRSSdIWiZpebpdnJLOr5h4Jd0gqTa7DTQXY6vij4h2fQD/Cfx9+vzvgf9s75iKxLg/8JnM9PXAzwABLwKfTudfAtzW3vFm4jwM+A+SCwNHdoB4fwhM54NjdwdW6jaSfpZbgJHp9ChgG8mPc8XEC3ya5FaVVfWxNveZtib+9v5nHEBytLg6na5Op3u194bSTNynAr8DjgSezczvCbzZ3vGlsXwYeBQYWL8RVXi83dL/fbeOsI2kiWQz8D/S6b8BXqjgeBsSSVMxtjb+9t612eWGQKD+hsCKJKmK5HqZ+ym4DSAiNgFVkvZvp/CyrgB+ERG1mXmVHO/BJF/MyyU9KWmRpPpf04rbRiL5lp0GzJf0EjAPmEyFxlugqRhbFX97J5KO6F+BN4EftXcgjZE0hqT1cXN7x9ICXUhu3lwWEUcAFwL3krRUKo6kLsA/ASdHRH/gROAuKjTettbeiaThhkBo6GWt2A2BFUHSDSS3W58eSb8r9Tcq1pf3JPmxeq2dQqw3DhgG1EpaRXJH6m9IOp6qxHghaSm9T9LfDRHxOLAJeJvK3EZqgN4R8QhA+vct4B0qM96spr53rfpOtmsiidJvCGx3kq4GDgcmRcS76eyngI+mTXBI7kG6uz3iy4qIayKid0QMiIgBJLez/y3JQeKKixcadrMeIu1AK+3e8wCS4w6VuI2sBfpKGgogaTjwcZJ7zSox3gZNfe9a/Z1szwNA6cGcYcDjJBvM48DQ9jeHPpcAAACASURBVI6pSIyHkPRDuyL9kJcD96VlxwLPkGxAC0nPNFTSg50PtFVsvCS7NovS+J4Gjq/kbYSke9FngD+kj0mVFi/JmbC1JK29DcCfmouxNfH7Enkzy629j5GYWSfgRGJmuTmRmFluTiRmlpsTiZnl5kRiZrk5kZhZbk4kZpbb/wd6r6F6jCDbOQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(4, 4))\n",
    "\n",
    "ax.plot(normed_cond_granger1[1, 2, :], color='b', lw=3, alpha=0.6, label=r\"Differential drive $0\\rightarrow 2\\vert$ 1\")\n",
    "ax.plot(normed_cond_granger2[0, 2, :], color='r', lw=3, alpha=0.6, label=r\"Sequential drive $1\\rightarrow 2\\vert$ 0\")\n",
    "ax.set_title(\"Disentanlged mediated drive\")\n",
    "ax.legend()\n",
    "ax.set_xlim((0, 100))\n",
    "ax.set_ylim((0, 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Replication of Figure 3 from Chen et al** When scaling the causality the the same order of magnitude as in the regular causality plots, we actually get values that are visually indistinguishable from Chen et al.'s results. Such a normalisation has not been mentioned in the article, but makes sense to be able to compare direct and indirect Granger-influence better. Not quite there, but close."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:py3]",
   "language": "python",
   "name": "conda-env-py3-py"
  },
  "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.3"
  },
  "toc-autonumbering": true,
  "toc-showcode": false,
  "toc-showmarkdowntxt": false,
  "toc-showtags": false
 },
 "nbformat": 4,
 "nbformat_minor": 4
}