diff --git a/LICENSE b/LICENSE index 8d27c7e..470e2af 100644 --- a/LICENSE +++ b/LICENSE @@ -1,6 +1,6 @@ MIT License -Copyright (c) 2019 T.Ikeuchi, G.Haraoka, S.Shimizu +Copyright (c) 2019 T.Ikeuchi, G.Haraoka, W.Kurebayashi, S.Shimizu Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal diff --git a/docs/conf.py b/docs/conf.py index aa31d70..cbf9d5c 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -21,8 +21,8 @@ # -- Project information ----------------------------------------------------- project = 'LiNGAM' -copyright = '{}, T.Ikeuchi, G.Haraoka, S.Shimizu'.format(datetime.datetime.now().year) -author = 'T.Ikeuchi, G.Haraoka, S.Shimizu' +author = 'T.Ikeuchi, G.Haraoka, W.Kurebayashi, S.Shimizu' +copyright = '{}, {}'.format(datetime.datetime.now().year, author) import lingam version = lingam.__version__ diff --git a/docs/reference/index.rst b/docs/reference/index.rst index 57a0ffd..689dd9b 100644 --- a/docs/reference/index.rst +++ b/docs/reference/index.rst @@ -11,6 +11,8 @@ API Reference multi_group_direct_lingam var_lingam varma_lingam + longitudinal_lingam bootstrap + longitudinal_bootstrap causal_effect utils diff --git a/docs/reference/longitudinal_bootstrap.rst b/docs/reference/longitudinal_bootstrap.rst new file mode 100644 index 0000000..b741d7c --- /dev/null +++ b/docs/reference/longitudinal_bootstrap.rst @@ -0,0 +1,8 @@ +.. module:: lingam + +LongitudinalBootstrapResult +=========================== + +.. autoclass:: LongitudinalBootstrapResult + :members: + :inherited-members: diff --git a/docs/reference/longitudinal_lingam.rst b/docs/reference/longitudinal_lingam.rst new file mode 100644 index 0000000..8a3cf01 --- /dev/null +++ b/docs/reference/longitudinal_lingam.rst @@ -0,0 +1,8 @@ +.. module:: lingam + +LongitudinalLiNGAM +================== + +.. autoclass:: LongitudinalLiNGAM + :members: + :inherited-members: diff --git a/docs/tutorial.rst b/docs/tutorial.rst index c791eaf..d7132bf 100644 --- a/docs/tutorial.rst +++ b/docs/tutorial.rst @@ -84,7 +84,7 @@ Then, we call :func:`~lingam.DirectLiNGAM.bootstrap` method instead of :func:`~l Causal Directions ^^^^^^^^^^^^^^^^^ -Since :class:`~lingam.BootstrapResult` object is returned, we can get the ranking of the causal directions extracted by :func:`~lingam.BootstrapResult.get_causal_direction_counts` method. +Since :class:`~lingam.BootstrapResult` object is returned, we can get the ranking of the causal directions extracted by :func:`~lingam.BootstrapResult.get_causal_direction_counts` method. .. code-block:: python @@ -305,7 +305,7 @@ we use lingam package and :func:`~lingam.utils.make_prior_knowledge`: .. code-block:: python import lingam - form lingam.utils import make_prior_knowledge + from lingam.utils import make_prior_knowledge First, we create a prior knowledge matrix: @@ -328,6 +328,12 @@ First, we create a prior knowledge matrix: [-1 0 1 -1 0 -1] [-1 0 -1 -1 -1 0]] +The values of the prior knowledge matrix elements are represented as follows: + +* ``0`` : :math:`x_i` does not have a directed path to :math:`x_j` +* ``1`` : :math:`x_i` has a directed path to :math:`x_j` +* ``-1`` : No prior knowledge is available to know if either of the two cases above (0 or 1) is true. + Then, if we use a prior knowledge, we set prior knowledge matrix to :class:`~lingam.DirectLiNGAM` object: .. code-block:: python @@ -394,7 +400,7 @@ Using the :attr:`~lingam.MultiGroupDirectLiNGAM.causal_order_` property, we can print(model.causal_order_) -Also, using the :attr:`~lingam.MultiGroupDirectLiNGAM.adjacency_matrices_` property, we can see the adjacency matrix as a result of the causal discovery. +Also, using the :attr:`~lingam.MultiGroupDirectLiNGAM.adjacency_matrices_` property, we can see the adjacency matrix as a result of the causal discovery. Since :attr:`~lingam.MultiGroupDirectLiNGAM.adjacency_matrices_` property returns a list, we can access the first matrix by indexing as follows: .. code-block:: python @@ -495,6 +501,23 @@ In the following example, we estimate the intervention value at variable index 1 Optimal intervention: 7.871 +Use a known causal model +^^^^^^^^^^^^^^^^^^^^^^^^ + +When using a known causal model, we can specify the adjacency matrix when we create :class:`~lingam.CausalEffect` object. + +.. code-block:: python + + m = np.array([[0.0, 0.0, 0.0, 3.0, 0.0, 0.0], + [3.0, 0.0, 2.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 6.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [8.0, 0.0,-1.0, 0.0, 0.0, 0.0], + [4.0, 0.0, 0.0, 0.0, 0.0, 0.0]]) + + ce = lingam.CausalEffect(causal_model=m) + effects = ce.estimate_effects_on_prediction(X, target, reg) + For details, see also https://github.com/cdt15/lingam/blob/master/examples/CausalEffect.ipynb https://github.com/cdt15/lingam/blob/master/examples/CausalEffect(LassoCV).ipynb https://github.com/cdt15/lingam/blob/master/examples/CausalEffect(LogisticRegression).ipynb @@ -601,7 +624,7 @@ The output of the :attr:`~lingam.VARMALiNGAM.adjacency_matrices_` property is as [-0.392, 0. , 0.182, 0. , 0. ], [ 0. , 0. , 0. , 0. , 0. ], [ 0.523, -0.149, 0. , 0. , 0. ], - [ 0. , 0. , 0. , 0. , 0. ]], + [ 0. , 0. , 0. , 0. , 0. ]], [[-0.145, -0.288, -0.418, 0.041, 0.592], [-0.324, 0.027, 0.024, 0.231, 0.379], [-0.249, 0.191, -0.01 , 0.136, 0.261], @@ -620,3 +643,43 @@ For example, we can draw a causal graph by using graphviz as follows: For details, see also https://github.com/cdt15/lingam/blob/master/examples/VARMALiNGAM.ipynb +LongitudinalLiNGAM +------------------ + +We use lingam package: + +.. code-block:: python + + import lingam + +First, if we use datasets from several time-points, we create a list like this: + +.. code-block:: python + + X_list = [X1, X2, X3] + +Then, if we want to run Longitudinal-LiNGAM algorithm, we create a :class:`~lingam.LongitudinalLiNGAM` object and call the :func:`~lingam.LongitudinalLiNGAM.fit` method: + +.. code-block:: python + + model = lingam.LongitudinalLiNGAM() + model.fit(X_list) + +Using the :attr:`~lingam.LongitudinalLiNGAM.causal_orders_` property, we can see the causal ordering in time-points as a result of the causal discovery. + +.. code-block:: python + + print(model.causal_orders_) + +Also, using the :attr:`~lingam.LongitudinalLiNGAM.adjacency_matrices_` property, we can see the adjacency matrix as a result of the causal discovery. + +.. code-block:: python + + t = 1 + print('B(1,1):\n', model.adjacency_matrices_[t, 0]) + print('B(1,0):\n', model.adjacency_matrices_[t, 1]) + + t = 2 + print('B(2,2):\n', model.adjacency_matrices_[t, 0]) + print('B(2,1):\n', model.adjacency_matrices_[t, 1]) + diff --git a/examples/DirectLiNGAM(Kernel).ipynb b/examples/DirectLiNGAM(Kernel).ipynb new file mode 100644 index 0000000..a9a70be --- /dev/null +++ b/examples/DirectLiNGAM(Kernel).ipynb @@ -0,0 +1,542 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# DirectLiNGAM by Kernel Method" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Import and settings\n", + "In this example, we need to import `numpy`, `pandas`, and `graphviz` in addition to `lingam`." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "ExecuteTime": { + "end_time": "2019-09-09T02:01:39.097825Z", + "start_time": "2019-09-09T02:01:33.841227Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "['1.16.2', '0.24.2', '0.11.1', '1.3']\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "import graphviz\n", + "import lingam\n", + "from lingam.utils import make_dot\n", + "\n", + "print([np.__version__, pd.__version__, graphviz.__version__, lingam.__version__])\n", + "\n", + "np.set_printoptions(precision=3, suppress=True)\n", + "np.random.seed(0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Test data\n", + "We create test data consisting of 6 variables." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "ExecuteTime": { + "end_time": "2019-09-09T02:01:39.159633Z", + "start_time": "2019-09-09T02:01:39.110757Z" + } + }, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
x0x1x2x3x4
00.2052601.0773220.2363190.102727-0.150883
1-1.121347-1.001230-3.7368390.562784-0.494015
20.8570910.1222820.0194670.230076-0.997795
31.237355-0.4922370.5687120.094054-0.133962
40.090289-0.558081-2.480684-0.165689-2.380579
\n", + "
" + ], + "text/plain": [ + " x0 x1 x2 x3 x4\n", + "0 0.205260 1.077322 0.236319 0.102727 -0.150883\n", + "1 -1.121347 -1.001230 -3.736839 0.562784 -0.494015\n", + "2 0.857091 0.122282 0.019467 0.230076 -0.997795\n", + "3 1.237355 -0.492237 0.568712 0.094054 -0.133962\n", + "4 0.090289 -0.558081 -2.480684 -0.165689 -2.380579" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "n = 1000\n", + "e = lambda n: np.random.laplace(0, 1, n)\n", + "x3 = e(n)\n", + "x2 = 0.3*x3 + e(n)\n", + "x1 = 0.3*x3 + 0.3*x2 + e(n)\n", + "x0 = 0.3*x2 + 0.3*x1 + e(n)\n", + "x4 = 0.3*x1 + 0.3*x0 + e(n)\n", + "X = pd.DataFrame(np.array([x0, x1, x2, x3, x4]).T ,columns=['x0', 'x1', 'x2', 'x3', 'x4'])\n", + "X.head()" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "ExecuteTime": { + "end_time": "2019-09-09T02:01:39.420171Z", + "start_time": "2019-09-09T02:01:39.165610Z" + } + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\r\n", + "\r\n", + "\r\n", + "\r\n", + "\r\n", + "\r\n", + "%3\r\n", + "\r\n", + "\r\n", + "x0\r\n", + "\r\n", + "x0\r\n", + "\r\n", + "\r\n", + "x4\r\n", + "\r\n", + "x4\r\n", + "\r\n", + "\r\n", + "x0->x4\r\n", + "\r\n", + "\r\n", + "0.30\r\n", + "\r\n", + "\r\n", + "x1\r\n", + "\r\n", + "x1\r\n", + "\r\n", + "\r\n", + "x1->x0\r\n", + "\r\n", + "\r\n", + "0.30\r\n", + "\r\n", + "\r\n", + "x1->x4\r\n", + "\r\n", + "\r\n", + "0.30\r\n", + "\r\n", + "\r\n", + "x2\r\n", + "\r\n", + "x2\r\n", + "\r\n", + "\r\n", + "x2->x0\r\n", + "\r\n", + "\r\n", + "0.30\r\n", + "\r\n", + "\r\n", + "x2->x1\r\n", + "\r\n", + "\r\n", + "0.30\r\n", + "\r\n", + "\r\n", + "x3\r\n", + "\r\n", + "x3\r\n", + "\r\n", + "\r\n", + "x3->x1\r\n", + "\r\n", + "\r\n", + "0.30\r\n", + "\r\n", + "\r\n", + "x3->x2\r\n", + "\r\n", + "\r\n", + "0.30\r\n", + "\r\n", + "\r\n", + "\r\n" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m = np.array([[0.0, 0.3, 0.3, 0.0, 0.0],\n", + " [0.0, 0.0, 0.3, 0.3, 0.0],\n", + " [0.0, 0.0, 0.0, 0.3, 0.0],\n", + " [0.0, 0.0, 0.0, 0.0, 0.0],\n", + " [0.3, 0.3, 0.0, 0.0, 0.0]])\n", + "\n", + "make_dot(m)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Causal Discovery\n", + "To run causal discovery, we create a `DirectLiNGAM` object by specifying 'kernel' in the `measure` parameter. Then, we call the `fit` method. " + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "ExecuteTime": { + "end_time": "2019-09-09T02:01:39.557802Z", + "start_time": "2019-09-09T02:01:39.423164Z" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "model = lingam.DirectLiNGAM(measure='kernel')\n", + "model.fit(X)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Using the `causal_order_` properties, we can see the causal ordering as a result of the causal discovery." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "ExecuteTime": { + "end_time": "2019-09-09T02:01:39.568772Z", + "start_time": "2019-09-09T02:01:39.560796Z" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[3, 2, 1, 0, 4]" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "model.causal_order_" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "ExecuteTime": { + "end_time": "2019-09-09T01:24:30.429100Z", + "start_time": "2019-09-09T01:24:30.422118Z" + } + }, + "source": [ + "Also, using the `adjacency_matrix_` properties, we can see the adjacency matrix as a result of the causal discovery." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "ExecuteTime": { + "end_time": "2019-09-09T02:01:39.583732Z", + "start_time": "2019-09-09T02:01:39.574757Z" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[0. , 0.34 , 0.273, 0. , 0. ],\n", + " [0. , 0. , 0.304, 0.275, 0. ],\n", + " [0. , 0. , 0. , 0.261, 0. ],\n", + " [0. , 0. , 0. , 0. , 0. ],\n", + " [0.26 , 0.239, 0. , 0. , 0. ]])" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "model.adjacency_matrix_" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can draw a causal graph by utility funciton." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "ExecuteTime": { + "end_time": "2019-09-09T02:01:39.863981Z", + "start_time": "2019-09-09T02:01:39.589716Z" + } + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\r\n", + "\r\n", + "\r\n", + "\r\n", + "\r\n", + "\r\n", + "%3\r\n", + "\r\n", + "\r\n", + "x0\r\n", + "\r\n", + "x0\r\n", + "\r\n", + "\r\n", + "x4\r\n", + "\r\n", + "x4\r\n", + "\r\n", + "\r\n", + "x0->x4\r\n", + "\r\n", + "\r\n", + "0.26\r\n", + "\r\n", + "\r\n", + "x1\r\n", + "\r\n", + "x1\r\n", + "\r\n", + "\r\n", + "x1->x0\r\n", + "\r\n", + "\r\n", + "0.34\r\n", + "\r\n", + "\r\n", + "x1->x4\r\n", + "\r\n", + "\r\n", + "0.24\r\n", + "\r\n", + "\r\n", + "x2\r\n", + "\r\n", + "x2\r\n", + "\r\n", + "\r\n", + "x2->x0\r\n", + "\r\n", + "\r\n", + "0.27\r\n", + "\r\n", + "\r\n", + "x2->x1\r\n", + "\r\n", + "\r\n", + "0.30\r\n", + "\r\n", + "\r\n", + "x3\r\n", + "\r\n", + "x3\r\n", + "\r\n", + "\r\n", + "x3->x1\r\n", + "\r\n", + "\r\n", + "0.28\r\n", + "\r\n", + "\r\n", + "x3->x2\r\n", + "\r\n", + "\r\n", + "0.26\r\n", + "\r\n", + "\r\n", + "\r\n" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "make_dot(model.adjacency_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.3" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": {}, + "toc_section_display": true, + "toc_window_display": false + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/LongitudinalLiNGAM.ipynb b/examples/LongitudinalLiNGAM.ipynb new file mode 100644 index 0000000..69bc193 --- /dev/null +++ b/examples/LongitudinalLiNGAM.ipynb @@ -0,0 +1,595 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Longitudinal LiNGAM" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Import and settings\n", + "In this example, we need to import `numpy`, `pandas`, and `graphviz` in addition to `lingam`." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "['1.16.2', '0.24.2', '0.11.1', '1.3']\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "import graphviz\n", + "import lingam\n", + "from lingam.utils import make_dot\n", + "\n", + "print([np.__version__, pd.__version__, graphviz.__version__, lingam.__version__])\n", + "\n", + "np.set_printoptions(precision=3, suppress=True)\n", + "np.random.seed(0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Test data\n", + "We create test data consisting of 5 variables. The causal model at each timepoint is as follows." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# setting\n", + "n_features = 5\n", + "n_samples = 200\n", + "n_lags = 1\n", + "n_timepoints = 3\n", + "\n", + "causal_orders = []\n", + "B_t_true = np.empty((n_timepoints, n_features, n_features))\n", + "B_tau_true = np.empty((n_timepoints, n_lags, n_features, n_features))\n", + "X_t = np.empty((n_timepoints, n_samples, n_features))" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\r\n", + "\r\n", + "\r\n", + "\r\n", + "\r\n", + "\r\n", + "%3\r\n", + "\r\n", + "\r\n", + "x0\r\n", + "\r\n", + "x0\r\n", + "\r\n", + "\r\n", + "x4\r\n", + "\r\n", + "x4\r\n", + "\r\n", + "\r\n", + "x0->x4\r\n", + "\r\n", + "\r\n", + "0.10\r\n", + "\r\n", + "\r\n", + "x1\r\n", + "\r\n", + "x1\r\n", + "\r\n", + "\r\n", + "x1->x0\r\n", + "\r\n", + "\r\n", + "0.50\r\n", + "\r\n", + "\r\n", + "x1->x4\r\n", + "\r\n", + "\r\n", + "-0.70\r\n", + "\r\n", + "\r\n", + "x2\r\n", + "\r\n", + "x2\r\n", + "\r\n", + "\r\n", + "x2->x0\r\n", + "\r\n", + "\r\n", + "-0.30\r\n", + "\r\n", + "\r\n", + "x2->x1\r\n", + "\r\n", + "\r\n", + "-0.30\r\n", + "\r\n", + "\r\n", + "x3\r\n", + "\r\n", + "x3\r\n", + "\r\n", + "\r\n", + "x3->x1\r\n", + "\r\n", + "\r\n", + "0.40\r\n", + "\r\n", + "\r\n", + "x3->x2\r\n", + "\r\n", + "\r\n", + "0.30\r\n", + "\r\n", + "\r\n", + "\r\n" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# B(0,0)\n", + "B_t_true[0] = np.array([[0.0, 0.5,-0.3, 0.0, 0.0],\n", + " [0.0, 0.0,-0.3, 0.4, 0.0],\n", + " [0.0, 0.0, 0.0, 0.3, 0.0],\n", + " [0.0, 0.0, 0.0, 0.0, 0.0],\n", + " [0.1,-0.7, 0.0, 0.0, 0.0]])\n", + "causal_orders.append([3, 2, 1, 0, 4])\n", + "make_dot(B_t_true[0])" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\r\n", + "\r\n", + "\r\n", + "\r\n", + "\r\n", + "\r\n", + "%3\r\n", + "\r\n", + "\r\n", + "x0\r\n", + "\r\n", + "x0\r\n", + "\r\n", + "\r\n", + "x1\r\n", + "\r\n", + "x1\r\n", + "\r\n", + "\r\n", + "x1->x0\r\n", + "\r\n", + "\r\n", + "0.20\r\n", + "\r\n", + "\r\n", + "x2\r\n", + "\r\n", + "x2\r\n", + "\r\n", + "\r\n", + "x1->x2\r\n", + "\r\n", + "\r\n", + "0.30\r\n", + "\r\n", + "\r\n", + "x4\r\n", + "\r\n", + "x4\r\n", + "\r\n", + "\r\n", + "x1->x4\r\n", + "\r\n", + "\r\n", + "-0.40\r\n", + "\r\n", + "\r\n", + "x2->x0\r\n", + "\r\n", + "\r\n", + "-0.10\r\n", + "\r\n", + "\r\n", + "x3\r\n", + "\r\n", + "x3\r\n", + "\r\n", + "\r\n", + "x3->x1\r\n", + "\r\n", + "\r\n", + "0.40\r\n", + "\r\n", + "\r\n", + "x4->x0\r\n", + "\r\n", + "\r\n", + "-0.50\r\n", + "\r\n", + "\r\n", + "\r\n" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# B(1,1)\n", + "B_t_true[1] = np.array([[0.0, 0.2,-0.1, 0.0,-0.5],\n", + " [0.0, 0.0, 0.0, 0.4, 0.0],\n", + " [0.0, 0.3, 0.0, 0.0, 0.0],\n", + " [0.0, 0.0, 0.0, 0.0, 0.0],\n", + " [0.0,-0.4, 0.0, 0.0, 0.0]])\n", + "causal_orders.append([3, 1, 2, 4, 0])\n", + "make_dot(B_t_true[1])" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\r\n", + "\r\n", + "\r\n", + "\r\n", + "\r\n", + "\r\n", + "%3\r\n", + "\r\n", + "\r\n", + "x0\r\n", + "\r\n", + "x0\r\n", + "\r\n", + "\r\n", + "x2\r\n", + "\r\n", + "x2\r\n", + "\r\n", + "\r\n", + "x0->x2\r\n", + "\r\n", + "\r\n", + "0.20\r\n", + "\r\n", + "\r\n", + "x4\r\n", + "\r\n", + "x4\r\n", + "\r\n", + "\r\n", + "x0->x4\r\n", + "\r\n", + "\r\n", + "0.30\r\n", + "\r\n", + "\r\n", + "x1\r\n", + "\r\n", + "x1\r\n", + "\r\n", + "\r\n", + "x2->x1\r\n", + "\r\n", + "\r\n", + "-0.70\r\n", + "\r\n", + "\r\n", + "x3\r\n", + "\r\n", + "x3\r\n", + "\r\n", + "\r\n", + "x2->x3\r\n", + "\r\n", + "\r\n", + "-0.40\r\n", + "\r\n", + "\r\n", + "x4->x1\r\n", + "\r\n", + "\r\n", + "0.50\r\n", + "\r\n", + "\r\n", + "\r\n" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# B(2,2)\n", + "B_t_true[2] = np.array([[0.0, 0.0, 0.0, 0.0, 0.0],\n", + " [0.0, 0.0,-0.7, 0.0, 0.5],\n", + " [0.2, 0.0, 0.0, 0.0, 0.0],\n", + " [0.0, 0.0,-0.4, 0.0, 0.0],\n", + " [0.3, 0.0, 0.0, 0.0, 0.0]])\n", + "causal_orders.append([0, 2, 4, 3, 1])\n", + "make_dot(B_t_true[2])" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "# create B(t,t-τ) and X\n", + "for t in range(n_timepoints):\n", + " # external influence\n", + " expon = 0.1\n", + " ext = np.empty((n_features, n_samples))\n", + " for i in range(n_features):\n", + " ext[i, :] = np.random.normal(size=(1, n_samples));\n", + " ext[i, :] = np.multiply(np.sign(ext[i, :]), abs(ext[i, :]) ** expon);\n", + " ext[i, :] = ext[i, :] - np.mean(ext[i, :]);\n", + " ext[i, :] = ext[i, :] / np.std(ext[i, :]);\n", + "\n", + " # create B(t,t-τ)\n", + " for tau in range(n_lags):\n", + " value = np.random.uniform(low=0.01, high=0.5, size=(n_features, n_features))\n", + " sign = np.random.choice([-1, 1], size=(n_features, n_features))\n", + " B_tau_true[t, tau] = np.multiply(value, sign)\n", + "\n", + " # create X(t)\n", + " X = np.zeros((n_features, n_samples))\n", + " for co in causal_orders[t]:\n", + " X[co] = np.dot(B_t_true[t][co, :], X) + ext[co]\n", + " if t > 0:\n", + " for tau in range(n_lags):\n", + " X[co] = X[co] + np.dot(B_tau_true[t, tau][co, :], X_t[t-(tau+1)].T)\n", + " \n", + " X_t[t] = X.T" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Causal Discovery\n", + "To run causal discovery, we create a `LongitudinalLiNGAM` object by specifying the `n_lags` parameter. Then, we call the `fit` method. " + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "model = lingam.LongitudinalLiNGAM(n_lags=n_lags)\n", + "model = model.fit(X_t)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Using the `causal_orders_` property, we can see the causal ordering in time-points as a result of the causal discovery." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[3, 1, 2, 4, 0], [0, 4, 2, 3, 1]]\n" + ] + } + ], + "source": [ + "print(model.causal_orders_)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Also, using the `adjacency_matrices_` property, we can see the adjacency matrix as a result of the causal discovery." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "B(1,1):\n", + " [[ 0. 0.099 0. 0. -0.52 ]\n", + " [ 0. 0. 0. 0.398 0. ]\n", + " [ 0. 0.384 0. -0.162 0. ]\n", + " [ 0. 0. 0. 0. 0. ]\n", + " [ 0. -0.249 -0.074 0. 0. ]]\n", + "B(1,0):\n", + " [[ 0.025 0.116 -0.202 0.054 -0.216]\n", + " [ 0.139 -0.211 -0.43 0.558 0.051]\n", + " [-0.135 0.178 0.421 0.173 0.031]\n", + " [ 0.384 -0.083 -0.495 -0.072 -0.323]\n", + " [-0.206 -0.354 -0.199 -0.293 0.468]]\n", + "B(2,2):\n", + " [[ 0. 0. 0. 0. 0. ]\n", + " [ 0. 0. -0.67 0. 0.46 ]\n", + " [ 0.187 0. 0. 0. 0. ]\n", + " [ 0. 0. -0.341 0. 0. ]\n", + " [ 0.25 0. 0. 0. 0. ]]\n", + "B(2,1):\n", + " [[ 0.194 0.2 0.031 -0.473 -0.002]\n", + " [-0.384 -0.037 0.158 0.255 0.095]\n", + " [ 0.126 0.275 -0.048 0.502 -0.019]\n", + " [ 0.238 -0.469 0.475 -0.029 -0.176]\n", + " [-0.177 0.309 -0.112 0.295 -0.273]]\n" + ] + } + ], + "source": [ + "t = 1\n", + "print('B(1,1):\\n', model.adjacency_matrices_[t, 0])\n", + "print('B(1,0):\\n', model.adjacency_matrices_[t, 1])\n", + "\n", + "t = 2\n", + "print('B(2,2):\\n', model.adjacency_matrices_[t, 0])\n", + "print('B(2,1):\\n', model.adjacency_matrices_[t, 1])" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAADQCAYAAADI+yJFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmcTfX/wPHXZwxjZ8aWlCSiUsTYRVkK2UbZvmQn35Jf31LpS4VvC/oiLfZdQtn3ZQzZzfAlRCRbdmEw1lnevz8+1zRTw4yZuXPunXk/H4/7cO+55577Prd59z7nfD7n8zEiglJKKaU8k4/TASillFLqzrRQK6WUUh5MC7VSSinlwbRQK6WUUh5MC7VSSinlwbRQK6WUUh5MC7WKZYx53hgz3w3b9TPG/GKMKZja21ZKJc6NuV3IGLPPGOOX2ttWf9JCnYEYY44YY64bYyKMMReNMUuMMQ/GWeVTYFCc9f9jjNltjIkyxvRPZNvPGWPWGGMuGWOOxH1PRG4CE4H3Um9vlFK3uTm3jTFmsDHmvOsxxBhjAETkDLAG6O6G3VIuWqgznsYikhMoDJwBvgIwxlQE8ojIljjrHgTeBZYkYbtXscX4nTu8/x3QQY+8lXIbd+V2d6AZUBZ4CmgEvBrn/el/ea1SmRbqDEpEbgCzgcddixoAP/5lnSkisgy4koTthYrINODQHd4/DlwEqqQkbqXU3aV2bgMdgKEiclxETgBDgY5x3t8KFDfGPJTS2FXCtFBnUMaY7EAr4PZR9pPAfjd/7T7sUblSyk3ckNtPAD/Fef2TaxkAIhKFPUPX3HYTX6cDUGluvjEmCsgJnAVecC3PS9KOrlPiiut7lFKpz125nRO4FOf1JSCnMcbIn5NFaG67kZ5RZzzNRCQv4Af0BH40xtyHvSydy83fnQsId/N3KJVRuSu3I4DccV7nBiIk/oxOmttupIU6gxKRaBGZC0QDNYBdwKNu/trHiH8JTSmVytyQ2z8T/7J2WdcyAIwxvkAJNLfdRgt1BuW65aIp4I9tO14K1PrLOpmNMVmxfye+xpisxphMrveKGWPEGFPM9drHtW5m1+azGmOyxNlWESCAP9vNlFJukNq5DUwF3jLGFDHG3A+8DUyOs7lKwBEROerG3crQtI0641lkjIkGBDgKdBCRnwFc90BXFpGtrnXHYXt83tYX6IRN0gddnz/heq8m9n7K265je5o+63r9D2CK655qpVTqc1dujwGKA7tdr8e7lt3WFhid6nujYpn4zQwqIzPGPA+8JiLNkrBuP+CciIxJwrp+2MtiNUXkbMojVUrdCzfmdkHsAfnTrtvClBtooVZKKaU8mKNt1MaYicaYs8aYPXd43xhjvjTGHDTG7DLGlE/rGJVS90bzWqnU5XRnsslA/bu83wAo6Xp0B0alQUxKqZSZjOa1UqnG0UItIuuAC3dZpSkwVawtQF5jTOG0iU4plRya10qlLk/v9V0E+D3O6+OuZafirmSM6Y5r9pYcOXJUKF26dJoFqFRqEIHff4dz5wC2/yEiBZyOyY2SlNegua28340bcPAg3LyZ/Lz29EJtElj2t95vIjIWGAsQGBgo27Ztc3dcSqWaCxegRQv43//g3XdhyBCT3u9HTVJeg+a28m6rVtnczp0bzp1Lfl473UadmOPYe/puewA46VAsSqW6AwegShVYvx4mT4bBg52OKE1oXqt0b+RIaNAAihaFsLCUbcvTC/VCoL2rl2gV4JKI/O3ymFLeaPVqqFwZLl6EkBDo0CHxz6QTmtcq3YqKgp494fXXbaHeuBEeSuEEoI5e+jbGzMCOXJXfGHMc+Ag7BCUiMho79F1D7BRq17Aj5yjl9UaPtslcujQsWgQPP+x0RKlH81plVOHh0LKlveT9zjvw2WeQKVPKt+tooRaRNom8L8DraRSOUm4XFQVvvQVffQUNG8KMGbb9Kj3RvFYZ0cGD0KgRHDoEEyZA586pt21P70ymVLoRHg6tWsHKlbZYDxmSOkfbSilnrVkDL70EPj4QHAw1a6bu9j29jVqpdOHgQaha1bZFjxsHQ4dqkVYqPRg3Dp5/HgoXhtDQ1C/SoGfUSrndjz9C8+b2+apV8OyzjoajlEoF0dHQuzd88QXUrw8zZ0KePO75Lj2jVsqNJkyAunWhYEHYulWLtFLpweXL0LixLdJvvmk7hLqrSIMWaqXcIjoa3n4bunaFOnVg82YoUcLpqJRSKXXokG3GWrXK3r0xfDj4uvnatF76ViqVXb4MbdrA0qXQq5dtj3Z3Iiul3G/9etuMFR0NK1ZA7dpp8716Rq1UKjp8GKpVs0k8ahSMGKFFWqn0YNIke3UsXz7bjJVWRRq0UCuVajZsgEqV4MQJW6h79HA6IqVUSkVH28FLOne2fUy2bIGSJdM2Bi3USqWCyZPtEba/vz3arlPH6YiUUil15QoEBcF//2uHBF26FPLmTfs4tFArlQLR0fDee9Cpk71/cutWePRRp6NSSqXU0aNQvbotzl9/bR9ONWNp65lSyRQRAW3bwsKF9jL3l19C5sxOR6WUSqlNm+yZ9M2bsGwZ1KvnbDx6Rq1UMhw7Zo+2Fy+243aPHKlFWqn0YNo0eO45Owb/li3OF2nQQq3UPdu8GSpWhCNH7GWxnj3BGKejUkqlREwM/Pvf0L69PQjfutXObucJtFArdQ+mT7dH2zlz2qPtF15wOiKlVEpdvQovv2ynpezWzd61ERDgdFR/0kKtVBLExEC/ftCuHVSpYgfff+wxp6NSSqXU779DjRqwYIEdEnTMGM9rxtLOZEol4upVezls7lzo0sW2R2fJ4nRUSqmUCg2Fpk1tji9eDA0aOB1RwvSMWqm7OH7c3nY1fz4MG2antNMirZT3mzkTatWCbNlsvxNPLdKghVqpOwoLsyON/fqrvQXrX//STmNKebuYGPjoIzsef8WK9qz6iSecjurutFArlYBZs+yZtJ+fvafyxRedjkgplVLXrkHr1jBwoB2kKDgY8ud3OqrEaaFWKg4R6N/fJnNgoD3aLlPG6aiUUil14oQ9+J492w4JOmGC9zRjaWcypVyuX4eOHeH776FDB9v708/P6aiUUim1fTs0aWKnoF24EBo1cjqie6Nn1EoBJ0/ajiU//ACDB9sp7bRIK+X9Zs+GZ56xt1xt2uR9RRq0UCvF//5nO43t3Qvz5sG772qnMaW8nQh8/DG0aAFPP22bsZ580umokkcvfasMbc4ceOUV26Fk40YoW9bpiJRSfxUVFcWSJUvYtGkThQsXpm3bthQoUOCO61+/bsc8mDHD5vfYsZA1axoGnMr0jFplSCLwySd22MCyZe3RthZppTzPtWvXqFu3Lp9++il58uRhx44dPPHEE2zZsiXB9U+ftsP8zphhhwSdMsW7izToGbXKgG7cgK5d7bjdbdvC+PHen8hKpVdfffUVefPmJSQkBB8fe245Z84cunTpwp49ezBx2ql27rSdxs6ftyMJBgU5FXXq0jNqlaGcOWOPtqdPt+1X06ZpkVbKk82fP59evXrFFmmA5s2bc/nyZQ4ePBhnPTvrlQhs2JB+ijRooVYZyE8/2ZGIdu2yPUH79tVOY0p5Ol9fX27duhVvWUxMDFFRUfj6+iJi79Ro3tyOeRAaajuPpSeOFmpjTH1jzH5jzEFjTJ8E3u9ojDlnjNnpenR1Ik7l/RYutEfbMTGwfj289JLTEaVvmtsqtbRu3ZrBgwdz8+bN2GXjx4+naNGi3H//w3TsCH36QKtWsHYtFC7sWKhu41gbtTEmE/ANUA84DoQZYxaKyN6/rDpLRHqmeYAqXRCBzz+3iRwYaKeyS4+J7Ek0t1VqevXVV1m/fj2lS5emYcOG7N+/n19//ZXp01dRu7a9N3rgQDsNbXq9QubkGXUl4KCIHBKRW8BMoKmD8ah05uZNO57ve+9By5bw449apNOI5rZKNb6+vsycOZPvv/+enDlzcv36dWJiHqdevdxs2xbF99/DBx+k3yINzhbqIsDvcV4fdy37q5eMMbuMMbONMQ8mtCFjTHdjzDZjzLZz5865I1blZc6dgzp17K0Z/fvbWzWyZXM6qgxDc1ulOj8/PyZOnMjjj7/LhQtLyJYtD/nzv8TlyxOcDs3tnCzUCR3/yF9eLwKKichTQDAwJaENichYEQkUkcC73QSvMoY9e+xIY9u32zlnP/oofR9teyDNbZXqPvnkU6pVm8OECU157DEfdu/OxuLF/fnoo4+IiopyOjy3cvI+6uNA3KPoB4CTcVcQkfNxXo4DBqdBXMqLLVliZ77KlQvWrbO9vFWa09xWqSIiIoKlS5dy5cpNli4NIiKiJi+/bK+UZc8ORYo8TUxMDKdOneLBBxO8KJMuOHlGHQaUNMY8bIzJArQGFsZdwRgTt0WxCbAvDeNTXkQEhg2Dxo2hZEl7i4YWacdobqsUW716NQ8//DBjx87l/fcrEBHRinLlFjBrli3SAKdOneL69evky5fP2WDdzLFCLSJRQE9gBTZJvxeRn40xA40xTVyr9TLG/GyM+QnoBXR0JlrlSa5du8Y333xDs2bNaN++PStXrqV7d3j7bTvIwfr18MADTkeZcWluq5S6evUqrVu35r//XcrhwzO5fPlxunVbw+7dLzFr1gzAFukuXbrQqVMnst+u3OmUEflr05F3CwwMlG3btjkdhkqmmJgY5syZww8//ADASy+9RIsWLWJHJbp+/Tp16tTB39+fDh06cPjwZT76qAw3b1ahb197m4aPlw/jY4zZLiKBTsfhaTS3M4aoqCh69+7N1KlniYgYT86cmVi2zI/KlaFJkyasWbOGXLlycf36dTp16sSgQYPIkiWL02EnKiV5rWN9K4/SrVs3du7cyRtvvIExhs8//5xly5YxadIkjDFMmzaNPHnysHjxYn75xfDvfwMI2bN35+23B+Pj4+/0LiilkikmJoYWLVqyYcPTXLw4lCJFLuDr25ylS2tTufIAatasyYMPPkifPn3Ily9fuj+Tvk0LtfIYoaGhhISE8PPPP8cmYIsWLShTpgyhoaFUrlyZ4OBg2rVrx8qVhpYt7Tjda9YY+vc/ysaNG2nkjbPCK6UAWLhwGevXt+H8+RZkzryYFSuepFCheZQuXZo2bdowadIkhg0blq47jiXEyy8SqvRk9erVvPTSS/GOkrNnz87LL7/MqlWrAPD3D+CHHwrRsCEUK2Y7jVWpIpw4cYKAgACHIldKpdSFC/DGGyU5f74F778Pw4Ydo27dKowYMYKiRYtSu3ZtKlWqxPPPP+90qGlOC7XyGP7+/pw8efJvy0+ePElAQACRkXD+/AAWLKjLs89GsHEjFC0qjBo1ChGhSpUqDkStlEqp/fuhcmU4dao4DRrM4NNPoWfP11ixYgWRkZH88ccfdO3alYkTJ8ab1jKj0EKtPEbLli1ZuXIlwcHBsctCQkJYvnw5zz/figYNYM6cQtSv/xP/+99DNGjwDKVLl2bUqFHMnz8/3jR4SinvsGqVLdKXLsG0aScJC+vF7t27AXjqqaeoXr06t27d4v3338+QRRq0jVp5kICAAH744QfatWtH4cKFMcZw4sQJhg5dyIsv5uPwYZg8GTp0KMuVK0cIDQ0lT548VKhQIcMmsFLebORI6NULHnsMFi2CYsWKAl/y7LPP8tRTT3H16lVOnTrF/PnzyZaBxwDW27OUx4mMjGTz5s2ICDduVKd1a198fWHePKhRw+no3E9vz0qY5nb6ERUFb74J33wDjRrBd9/Z0QRvi4iIYP369WTNmpVnnnkGX1/vP6fU27NUupI5c2Zq1qzJ6NHQsyeULm2Pth9+2OnIlFIpFR5uZ7NbtQreeQc++wwyZYq/Ts6cOWnQoIEzAXogLdTK40RFwVtvwVdfQcOGduar3LmdjkoplVIHD9oz6EOHYMIE6NzZ6Yi8w10LtTHmrbu9LyLDUjccldGFh0OrVrBypS3WQ4b8/WhbpZzmtkpra9bASy/ZkQODg6FmTacj8h6JnVHfbjUoBVTkz4H1GwPr3BWUypgOHrSTahw8COPGQdeuTkeUrmluqzQzbhy89ho8+qhtxipe3OmIvMtdC7WIDAAwxqwEyovIFdfr/sAPbo9OZRg//gjNm9vnq1bBs886Gk66p7mt0kJ0NPTuDV98AfXr2/nh8+RxOirvk9QbT4sCt+K8vgUUS/VoVIY0YQLUrQsFC8LWrVqk05jmtnKLy5ftFbIvvrA9vBct0iKdXEntTDYNCDXGzAMECAKmui0qlSFER8O779p5pF94wR5t583rdFQZjua2SnWHDtkifeAAjB4Nr77qdETeLUmFWkQ+McYsA55xLeokIjvcF5ZK7y5fhjZtYOlSO+DB0KGQDm6V9Dqa2wrgypUrzJw5k0OHDlGuXDmCgoKSPXXk+vW2GSs6GlasgNq14dixY8TExPDQQw/p4ETJcC9jLmYHLovICOC4MUbvalXJcvgwVKtmk3jUKBgxQou0wzS3M7D9+/fz+OOPs3z5cnLmzMmoUaOoXLkyFy5cuOdtTZoEdepAQIBtxrrvvr1UqVKFwMBAqlSpQoUKFdixQ48D75mIJPoAPgIWAQdcr+8HNibls2n9qFChgijPtX69SP78InnzigQHOx2NZwK2SRrli+a2qlu3rowYMSL2dUxMjHTr1k3efPPNJG8jKkqkd28REKlbV+TCBZGrV6/KAw88IGPGjJGoqCiJjo6WadOmSaFChSQ8PNwdu+LRUpLXST2jDgKaAFddxf0kf97eoVSSTJlij7b9/e3Rdp06Tkek0NzO0K5cucKmTZt4NU4jsjGGf/3rX8ybNy+J24CgIPjvf+0tWIMG7eKTT3rTqFEjChQoQOfOncmUKRM+Pj60a9eOWrVqMXPmTHftUrqU1EJ9y3VEIADGmBzuC0mlN9HR8N570LEjPPOMLdKPPup0VMpFczsDuz3jXFRUVLzlkZGRSRpf++hRqF7d9jX5+msIDJzEiy8+T44cOcidOzcnT54kKCgo3vZLly7NiRMnUndH0rmkFurvjTFjgLzGmG5AMDDefWGp9CIiwnYsGTIEevSAZcvsGbXyGJrbGViOHDmoW7cuQ4cOjV0WHR3NJ598QqtWre762U2boFIlOHYMRo48zPLljenSpQvZs2cnICCAN998kzx58nDu3Dm+//772G0vXrxY546/R0nt9f1fY0w94DJ2JKMPRWSVWyNTXu/YMXuLxp49dtzu11+HtWvXMHHiRMLDw6lXrx5dunQhRw49iXOK5rb6+uuvqVevHsHBwZQvX57g4GAKFizIxIkT7/iZadPsyIFFi8I33xzjlVeq0KRJE5555hkGDRrE22+/zfHjx3n00Uc5cuQIEyZMoEiRIgwfPpw8efLwwgsvpOEepgNJacgGBidlmSc8tMOJZ9i0SaRgQZHcuUWWL7fLhg8fLg899JB8/fXXMnv2bGncuLFUqlRJrl696mywHoa07Uymua3k1q1bMm/ePBk2bJiEhIRITExMgutFR4u8/74IiDz3nMj58yI9e/aUvn37yurVq6VixYoiInL27FnJmzevnDx5Uho3biwFChSQwMBA+eyzz+TatWtpuWseIyV5ndRk/l8Cy3Yl90vd+dBkdt6334r4+YkULy6yd69dduHCBcmbN68cPXo0dr2YmBhp2LChjBw50qFIPVMaF2rNbZUkEREiQUG2anTrJnLrll1es2ZNCQkJkcjISClatKjMmTNHREQqVqwoCxYskKJFi8q6descjNwzpCSv79pGbYz5pzFmN1DKGLMrzuMwsCv1z++VN4uJgX79oF07qFIFQkPhscfse1u2bKFChQoULVo0dn1jDO3atWP16tUORZxxaW4rsL2+f/nlF65evZrg+yJCWFgYX3wxh/Llr7FggR0SdMwYyJzZrlOiRAnCwsLw9fVl7ty5vPHGG9SoUYOdO3fSvn17Xn/9dZ555pkEt6+SJrE26u+AZcBnQJ84y6+IyL3fDa/SratXoX17mDsXunSBkSMh7sBGAQEBnDhxAhGJNzLR8ePHyZcvnwMRZ3ia2xlYdHQ0ffr0Yfz48RQoUIA//viDN998kw8++CA2P69cuULz5s3ZuzcnFy5M5ObNaJ5+uh9duvTBmJyx2+rVqxf16tWjZMmSNGnShHXr1tG2bVuqVq3KrFmzuO+++5zazXTjrmfUInJJRI6ISBsROQpcx97GkdMYU/Run1UZx/Hjdm7Z+fPtuN3jxsUv0gCVKlXC19eXr7766vblVfbv388XX3xBZ509Ps1pbmdMJ0+eZMeOHfTv359t27axb98+Dhw4wM6dO1m0aBEjR46MXfff//43N24049y5H8iU6RZ1635Itmxr6du3b7xtli1blpkzZ/LJJ5+QI0cOKlSoQPXq1Vm5cqUW6dSSlOvj2Dlqf8UOinAYiAF+Tu71dnc+tB0rbYWGihQuLJIrl8jixXdf9+DBg1K2bFkpUaKEVK9eXfz9/WXcuHFpE6gXIW3bqDW3M4Dw8HBp3ry5BAQESJkyZcQYIx988EG8dTZt2iSlS5cWEdtpzM9vkIBIjhzb5csvv5PJkyfLk08+KX5+fnfsbBYRESGRkZFu3x9vlJK8TuoIyx8DVYBgEXnaGPMc0Cb1DheU03777TfWrVtHvnz5qF+/foID8vv6+hIdHQ3AhQsXWLXKnw4d4L77YOVKKFPm7t/xyCOPsGPHDnbs2EF4eDiVKlUiZ86cd/+QcjfN7QygW7du5MuXj+PHj+Pr60u2bNn49ttvqVy5Mi+++CLw50Ak165Bx47CzZvvkSPH9/z+eyP8/csD8Pzzz1OkSBG2bt2a4L3QequleyR1wJNIETkP+BhjfERkDVAupV9ujKlvjNlvjDlojOmTwPt+xphZrve3GmOKpfQ7VXwiwttvv02VKlUICQlh6NChlChRgt27d8eu88knn2CMiS3SAAEBX9CqFQQG2k5jiRXp24wxlC9fntq1a2uR9gya2+nc2bNnCQ4OZvjw4WTLlo3MmTNTtmxZgoKC4l3qXrBgAeXKvUjNmjB7tqFgwc8pX/4b/P2zx64zdepUSpYsyapVeqt9WkrqGXW4sb0H1gHTjTFngahEPnNXxphMwDdAPeA4EGaMWSgie+Os1gW4KCIljDGtgcHA3YfLUfdk7ty5rFq1il9//ZW8rsmgp06dSuvWrdmzZw/GGPr16wdAcHAw1arVoW3bW8yblwWYTHBwR/z8HNwBlVKa2+ncuXPnKFCgQLyz3U8//ZQ2bdqQK1cuduzYwdq1axk4cDG+vsu4cQMWLIAdO64zaFAYHTt2pFq1amzYsIGQkBACAwMJCAhwcI8ynqSeUTfFdjb5F7Ac+A3btpUSlYCDInJIRG4BM13f89fvneJ6PhuoY3Qy01Q1ffp0evfuHVukAV555RUiIyPZuXNn7LIsWbLw2GN1qFUL5s/PwuOPTwE6MXnyGAeiVqlIczudK1myJJcuXWLXrj/vunvhhRdo0KABWbJkoWPHjsyd68ONG6vInj0LGzfaEQV79uxJtmzZyJIlC6GhoZQtW5ZvvvmGDRs2JDq8qEpdSR1C9CqAMSY3dkq81FAE+D3O6+NA5TutIyJRxphLQD7gj7grGWO6A92BePfpqsTduHGDXLniT5ZkjCFnzpzcuHEjdln27DWoVAnCw2HePDh9+gY9esDvv//+100qL6K5nf5lyZKFzz77jEaNGvHhhx9SqlQp5s6dy7p169i0aTNTpjzABx/YOeLnzYOCBe3nAgICmD17Nu3ataNw4cL89NNPnDhxglmzZpE/f35ndyqDSVKhNsa8CgzEHnnHAAZ7K0fxFHx3QkfPkox1EJGxwFiAwMDAv72v7qxRo0aMHj2aJk2akClTJgA2b97MqVOnqFChgmut5oSHTyNXLti4EcqWBWN6APDxxx87FLlKDZrbGUOnTp146KGHGDlyJBMnTqRatWqsWbOFd9+9n5kzoW3bGMaP9yFr1vife+655zhy5AibN29GRKhWrRqZb490otJMUtuoewNPiMgfia6ZdMeBB+O8fgA4eYd1jhtjfIE8gA7GkIo6d+7M3LlzqVGjBq1ateLYsWNMmzaNiRMnkjlzFj75BGAOsJnffw+iXLkzsZ/108bp9EBzO4OoXbs2tWvXBuD0aahe/SyHDhmyZRvAsmVfMWhQz9gBT7Zu3UpERARVq1YlZ86c1KxZ0+HoM7akFurfgGup/N1hQEljzMPACaA18I+/rLMQ6ABsBl4GQlz3o6lUkjVrVpYvX878+fNZu3Yt+fLlY+vWrdx/f3FeeQWmT4c6dU6zcWNDbtwIj/2cn58fo0aNcjBylUo0tzOYnTuhXr1rnD+fkyFDDvLOOx9x6NArdOjQgXPnzhEcHEyWLFkICAhgz549DBs2jA4dOjgddsaWlJutgaeBncAY4Mvbj+TevB1nuw2BA9j/WfR1LRsINHE9zwr8ABwEQoHiiW1TB0VIudOnRapUEQGRjz8Wee652vLdd9/FW2fr1q1SrFixOw58oJKPtB3wRHM7A5k3TyR7dpEsWU7L8OFr47134MAB8fHxkfHjx8fm9d69e6VQoUKyc+dOJ8JNV1KS10aScBBrjAkFNgC7se1Yt4v8lDt+yCGBgYGybds2p8PwWj/9ZHt8nj8PU6fCSy9BgQIF2LNnD4UKFYq3bs6cOTl58iS5c+d2KNr0yRizXUQC0+i7NLczABEYMgTefx8qVoQjR8qxfftiHnjggdh1goODqV+/PmfOnIk3/v6AAQMIDw9n+PDhToSebqQkr5N66TtKRN5Kzhco77FwIfzjH5A3L6xfD+XtYESUKFGC0NBQGjf+866dvXv3kitXLh20xPtpbqczly5dYtGiRdy4cYP69etToMADdO9uD7xbt4aJE6FFiwdYsWIFXbp0if1cWFgYWbJkwd/fP9727r//fg4fPpzWu6HiSOp91GuMMd2NMYWNMQG3H26NTKWZ20fbzZrB449DWNifRRrgnXfeoVevXmzcuBER4eeff+aVV17h7bffxscnqX9CykNpbqcjy5Yto3jx4syZM4cff/yRMmVqU7r0CaZOhYED4bvvIFs26NevH++//z7jx4/n999/Z8HB/RSiAAAaZElEQVSCBYwePRqA06dPx24vJiaGb7/9lnr16jm1SwqS3EZ9OIHHoeReb3fnQ9ux7s2NGyIdO9r26FatRK5dS3i9qVOnSpEiRcTX11f8/f1lyJAh2j7tJqRtG7Xmdjpx+fJlCQgIkE2bNomIyK5dIg88EClwTQYPPvS39Tdt2iSNGjWSIkWKSI0aNWTu3Lny2WefycMPPyxfffWVfPfdd1KnTh2pXbu23Lx5M613J91JSV4ndcCTh1P7AEE579w5aN4cNmyA/v3hww8hobGhIiIimDBhAnnz5qVJkybs3LmTSZMm0bZtW+6///40j1ulHs3t9GPp0qVUrVqVqlWrsngxtGkDuXP70rXrBP744zdgSLz1q1atyqJFfx/jplKlSkyZMoWIiAjatGlDu3btEpykR6WduxZqY0xtEQkxxjRP6H0RmeuesJS77dljO42dPg0zZ8LdRgQcOHAgRYoUISQkJPZS9wcffMAbb7zBnDlz0ihilZo0t9OfmzdvkjVrNoYOhXfegaeftv1Opk27yOnTt5K8nbj3WyvPkNgZdS0ghITH/hVAk9kLLVlij7Zz5oR162wv0LuZNWsWy5Yti9ce/e6771KwYEFu3LhB1r8OZ6S8gea2Fzh8+DBff/01e/bsoVSpUvTs2ZNHH3003joxMTHs3r2bggUfYMGCRsyZA2XL/srs2bnJkycHkyZN0jEPvNxdC7WIfOR6OlBE4nX7cw1moLyICAwfDr17Q7ly9mg7zt0ZdxQVFfW3S1++vr6ubeoYFd5Ic9vz7dq1i7p169K5c2fefPNNtmzZQvXq1Xn11VdZt24dly5dokyZMmzdupWYmACOHRtGdHQHMmX6lKiomTzxxEH8/f1p3LgxRYoU4euvvyZXrlw0a9aMPHnyOL176h4ktctuQtc3Z6dmIMq9bt2C7t3h7bchKMjefpWUIg0QFBTE8OHD4xXlUaNG8eyzz5ItWzY3RazSiOa2h+rXrx8ffvghgwYNokGDBgwYMICKFSsyYsQI3nvvPb744gvmzZvHpUtFuHIlmJiYivTqtZk8eYbSsGF9mjZtyvXr18mVKxc1a9Zk165dLFiwgEceeYQ1a9Y4vXvqHiTWRl0aeALI85e2rNzYkYWUm0VHR3Pz5k2yZctGcmcBPH/eDlzy44/Qt6+9TeNe7qrq378/tWvXpm7dutSpU4ft27cTFhbG6tWrkxWPcp7mtudbs2YNU6b8Oe7MyZMn2bRpE9euXaN+/fpMmjSJwMC+bNnyL0Su0a/fegYOfJHTp+tSokQJhgwZQqlSpZg5cyb79u2LnUM6JCSENm3acPToUR2v30sk9r/rUkAjIC+2Lev2ozzQzb2hZWxRUVH069ePQoUKERAQwNNPP83SpUvveTv79kHlyrBlC3z7LXz88b0VaYD8+fMTFhZG165duXz5MvXr1+fnn3+mZMmS9xyP8hia2x4uX7588aaR3b59O+XKlSN37twY48P8+Q+yYcP73H//TfLnb0ixYnbCnOLFi3PmjH1+/vx5goKCYos02M5iJUuW1LNqL5JYG/UCYIExpqqIbE6jmBTQu3dv9u7dS1hYGA899BArVqygU6dOzJs3j6pVqyZpGytWQMuWkDUrrFkDSfxYgvz8/GjTpg1t2rRJ/kaUx9Dc9nzdunXjrbfeYu7cueTOnZt8+fIRFhbGq6/25LXXDEuWvECOHKto0WIT+/YV5MsvvyQoKIh58+YxduxYVq5cybVr1yhVqtTftu3n50dUVJQDe6WSJSk3W2NvwMsNZAZWYyd3b5fcm7fd+UgPgyKEh4dLnjx55OzZs/GWjxw5Ulq0aJHo52NiRL78UsTHR6RsWZGjR90VqXIH0nbAE81tDxUZGSk9evQQf39/qVmzpuTLl0/y5SspDz64X0CkT58YqVGjpmTOnFlGjBghzz//vGTJkkUeeeQRefnllyV//vzy8ccfS7ly5eRanJGMduzYIQEBARIREeHg3mU8KcnrpI71/byIvGuMCcLOI9sCWAN8m2pHDCrWiRMnuO+++yhQoEC85ZUqVWLs2LF3/WxkJPzf/8GoUdCkiZ2mUofjVnehue2hfH19GTVqFB9++CH79+8HStG5cwGOHhWyZ+/BtGmL8fPzo0aNGgwYMACA5557jho1alC4cGHGjRtH7ty5OXDgAOXKlaNNmzb88ccfzJgxg9GjR5MjRw5nd1AlWVJbKzO7/m0IzBARneDdjR566CHOnj3LsWPH4i1fvXo1ZcuWvePnLl6EBg1skX7vPZg3T4u0SpTmtocrXLgwkZHP0qxZYSIifFm3LjNHjvyHH3/8kXLlyiEiTJo0iSlTphAVFUVYWBidO3cmb968+Pj4MHnyZMaOHcvNmzcpUqQI27dvp0WLFk7vlroHST2jXmSM+QW4DrxmjCkA3HBfWBlbjhw5+Ne//kVQUBAjRoygVKlSzJs3j88//5yQkJAEP3PggB1p7PBhmDwZdJ53lUSa2x5u5Ejo1Qseewxmz77J5s0zmD49lOjoaHbu3Mm+fftixzl44YUXKFu2LGvXruW5554DwBhDrVq1qFWrlpO7oVIgSWfUItIHqAoEikgkcA1o6s7AMrp+/frRvXt3/vnPf1KqVCkWLlzIsmXLePLJJ/+27urVtmf3hQsQEqJFWiWd5rbnioqCnj3h9dftlbIlSy7SunUVZsyYwWOPPcbOnTs5deoU27dvj/1M5syZadasGevXr3cwcpXa7lqojTHvxnlZV0SiAUTkKtDLnYFldMYYXn31VXbv3s2FCxdYvHgxgYF/n3N8zBh44QUoUgRCQ6FGDQeCVV5Hc9sziQgbN25k+vQl1Klzk2++sSMJzp8Po0d/Trly5Vi+fDlvvPEGr776Kk899RQ9evS43TEQgEOHDlGoUCEH90Klurv1NAP+l9DzhF57yiOj9AyNjBR54w0REGnYUOTSJacjUqmFNOj1rbnteQ4dOiRPPvmklChRX3LkOCZwU5o0mR87neyTTz4poaGhsetfunRJChUqJLlz55ajR49KTEyMzJ49WwoWLCgXL150ajfUHaQkrxO79G3u8Dyh1yqNhIfDiy/CV1/BW2/ZMbtz53Y6KuVlNLc9TJs2bahWrS/nzy8la9YHmT//KgcO9GHhwoUAZM2alYiIiNj1c+fOzdy5c4mIiKBGjRqUKFGCvn37snDhQvLmzevUbig3SKwzmdzheUKvVRo4eNB2Gjt4EMaNg65dnY5IeSnNbQ/yyy+/sH9/TbZvb8mjjxoWLYLixf25dOl9Jk+eTNOmTWnbti3/+c9/qFq1auyMdWvWrKF27dp88cUXiAhPPPFEsocaVp4rsUJd1hhzGXuEnc31HNdrHQ84jf34IzR3jcq8ahU8+6yj4SjvprntIaKjYcCAvISHD6F+fTs//O3JrQoUKMClS5cAeP311wkLC+ORRx6hXr16/Pzzz1y7do1ly5ZRtGhRB/dAuVtiQ4hmSqtA1N1NmAA9ekCJErBokf1XqeTS3PYMly9D69awbNl9ZM8+hv79A8mTpwJg+w9NnDiR+vXrA3YAlG+//ZY9e/YQGhpK27ZtqV27Npky6X/K9C6p91Erh0RHw7vvwrBhtnf3zJmgzU9KeYeoqCimTp3KvHnzyJQpEy1btqR169b4+Phw6JBtxjpwAEaPBn9/f5o0achrr73Gww8/zKxZszhz5gwTJ06Mt80yZcpQpkwZh/ZIOUELtQe7fBnatIGlS+2AB0OHgq/+F1PKK4gILVu25Ny5c7zxxhtERUUxbNgwQkJCaN9+PM2bQ0yMnTyndm2AlpQuXZoJEybwyy+/0LhxY9q3b0/27Nmd3hXlMP3fvoc6fNgebf/yix0StEcPpyNSSt2LkJAQ9u/fz44dO2JHDmvatCkPPvghU6bEULy4D4sXQ9zZYp966ilGjBjhUMTKU2mh9kAbNkBQkB2ZaMUKqFPH6YiUUvdqzZo1tGjRIrZIR0dD//45uHhxKCVLHmXLlofw93c4SOUVkjoph0ojU6bYwuzvD1u3apFWylvly5cvdmKdK1fswfd//wvFii2lT591WqRVkmmh9hAxMdCnD3TsCM88Y4v0o486HZVSKrn+8Y9/sHDhQr77biPVq8OSJTEUKDCAo0cbMXXqBNasWeN0iMpLOHLp2xgTAMwCigFHgJYicjGB9aKB3a6Xx0SkSVrFmBYiIyM5f/48fn756NQpMwsW2LboL7+EzJkT/7xSnkZz+0+FChXiww+X88orDyFyCR+fNmTKtIMlS5YQHh5O69atmTVrFs/qgAgqEU61UfcBVovIIGNMH9fr9xJY77qIlEvb0NxPRBg+fDhDhgwhMrIwly9PIzr6cb780tCzp0EHFlJeLEPndlzTpsE77wRSrJhw5UothgzpTPv2i/HxsRcyY2Ji+Pjjj7VQq0Q5dem7KTDF9XwK0MyhOBwxZswYpkyZwvDhW8iceQdZsz5G8eJvAF9rkVbeLkPnNthmrH//G9q3h2rVIDj4Ctev76Bjx46xRRqgbt26/PTTTw5GqryFU4W6kIicAnD9W/AO62U1xmwzxmwxxtwx4Y0x3V3rbTt37pw74k1Vw4YNo3nzOXTqVIwcOSA0NBOzZnVm2LBhToemVEpl6Ny+ehVefhk++wy6dYOVK6Fo0RzkyJGDffv2xVs3LCyMEjrEoEqK5E67ldgDCAb2JPBoCoT/Zd2Ld9jG/a5/i2Pbux5J7Hs9fSq86GgRH59PBURq1RL54w+7PDIyUowxsVPaqYyLNJjmMiUPze2EHTsmUq6ciI+PyBdfiMRN5cGDB0uFChVk7969EhMTIxs3bpRixYrJ3LlznQtYpamU5LXb2qhFpO6d3jPGnDHGFBaRU8aYwsDZO2zjpOvfQ8aYtcDTwG/uiDctXL1qL4fFxLxP3bpHWLKkGK5bLFm+fDnly5fXmW+Ux9Pc/rvQUGja1Ob44sXQoEH899955x2MMdSuXZvLly9z3333MXDgQIKCgpwJWHkVpy59LwQ6uJ53ABb8dQVjjL8xxs/1PD9QHdibZhGmsuPHoWZNmD8funf/hV27qjBr1jSOHDnC9OnT6datGwMGDHA6TKVSKsPl9syZUKsWZMsGmzf/vUgDGGN45513OHHiBGfOnOHgwYO88soraR+s8kpOFepBQD1jzK9APddrjDGBxpjxrnUeA7YZY34C1gCDRMQrkzksDCpVgl9/hYULYcyY0syaNZPvvvuOmjVrMnXqVKZPn86LL77odKhKpVSGye2YGPjoIzsef8WK9qz6iSfu/hkfHx9y5sypV87UPTH20nn6ERgYKNu2bXM6jFjffw8dOsB999npKXXSG5UYY8x2EQl0Og5P40m5fe2aHZzohx+gUyc7Hr+fn9NRKU+WkrzWkcncRAQGDIBWrSAw0B5ta5FWyvudOGGbsWbPhs8/t3PFa5FW7qSTcrjB9ev2KHvWLHs2PWaMJrJS6cH27dCkiZ2CdsECO8OdUu6mZ9Sp7NQp27Hk++9h8GCYNEmLtFLpwezZdhx+X1/YuFGLtEo7WqhT0Y4dtlPJ3r0wbx68+y460phSXk4EPv4YWrSAcuVsM9ZTTzkdlcpItFCnkrlzoUYN8PGxR9tNmzodkVIqpa5fh7Zt4YMPoF07CAmBQoWcjkplNFqoU0gEPvkEXnrJHmWHhkLZsk5HpZRKqdOn4bnnYMYM+PRTmDoVsmZ1OiqVEWlnshS4cQO6doXp0+1R9/jxmshKpQc7d9pOY+fP26tlOoCYcpKeUSfTmTP2aHv6dNt+NW2aFmml0oP586F6dXu1bMMGLdLKeVqok+Gnn2ynsV27bE/Qvn2105hS3k4EBg2C5s3tmAehofD0005HpZQW6nu2cKE92o6JgfXrbdu0Usq73bxpxzx4/307SNHatVC4sNNRKWVpoU4iERgyBJo1g8cft+N3ly/vdFRKqZQ6exZq17bNVwMHwnff2Qk2lPIU2pksCW7ehB49YPJke7Q9aZImslLpwe7dduCSs2ftIEUtWjgdkVJ/p2fUiTh3DurWtUW6f397q4YWaaW83+LFUK0aREbCunVapJXn0kJ9F3v22Okpt22zc85+9JF2GlPK24nA0KH29qtHH7WdxgJ1rjLlwbRQ38HSpfZo++ZNe7TdqpXTESmlUurWLTv2Qe/etiPo+vVQpIjTUSl1d1qo/0IEhg+37VYlStij7YoVnY5KKZVSf/wB9erBxIl2SNBZsyB7dqejUipx2pksjlu3oGdPGDfO3ks5dSrkyOF0VEqplNq71x58nzhhByn6xz+cjkippNMzapfz5+H5522R7tsXfvhBi7RS6cHy5VC1Kly9Cj/+qEVaeR8t1MC+fVC5MmzZAt9+a4cE9dFfRimvJgJffgkvvggPP2ybsSpXdjoqpe5dhi9HK1ZAlSpw5QqsWWMn11BKebfISPjnP+H//s/27t6wAYoWdToqpZInwxZqEfjqK2jY0B5th4XZy2NKKe924QLUrw9jxtghQefMgZw5nY5KqeTLkJ3JIiPtkfaoUfZoe/p0TWSl0oP9+6FRIzh2zHYGfeUVpyNSKuUyXKG+eNGOQLR6Nbz3np0QXtujlfJ+q1bZ3M6SBUJC7OQ5SqUHGapEHThg26PXrbNDgg4apEVaqfRg5Eho0AAefNB2GtMirdKTDFOmVq+2PT4vXLBH2x06OB2RUiqloqLs2Aevv24L9aZNUKyY01EplboyRKEeMwZeeMEOFRgaCjVqOB2RUiqlwsNtZ9BvvrFDgs6fD7lyOR2VUqkvXbdRR0XB22/beykbNrQzX+XO7XRUSqmUOnjQdho7dAgmTIDOnZ2OSCn3SbeF+tIlO5HGihXw1lswZAhkyuR0VEqplFqzxk6o4eNjO5DVquV0REq5lyOXvo0xLYwxPxtjYowxd5xgzhhT3xiz3xhz0BjTJ6nb/+03e0/06tV2SNChQ7VIK5UW3J3b48bZoX7vuw+2btUirTIGp9qo9wDNgXV3WsEYkwn4BmgAPA60McY8ntiGr1yxc0ifOWOPtrt2Ta2QlVJJ4Lbc/v136N4d6tSBzZvhkUdSK2SlPJsjhVpE9onI/kRWqwQcFJFDInILmAk0TWzbBw5AwYL2aPvZZ1MhWKVUkrkzt8+etQMVLV4MefKkRrRKeQdPbqMuAvwe5/VxIMEh9Y0x3YHurpc3f/nF7ClZ0s3Rpa78wB9OB3EPNF73KuV0AG6W7NweMcLsGTHCzdGlHm/7u9N43SvZee22Qm2MCQbuS+CtviKyICmbSGCZJLSiiIwFxrq+d5uI3LFtzBN5W8war3sZY7Y5HcPdaG4njcbrXt4Yb3I/67ZCLSJ1U7iJ48CDcV4/AJxM4TaVUimkua1U2vLkAU/CgJLGmIeNMVmA1sBCh2NSSqWc5rZS98Cp27OCjDHHgarAEmPMCtfy+40xSwFEJAroCawA9gHfi8jPSdj8WDeF7U7eFrPG617eFm8sze14NF73yjDxGpEEm4aUUkop5QE8+dK3UkopleFpoVZKKaU8mNcXancPWZjajDEBxphVxphfXf/632G9aGPMTtcjzTvaJPZ7GWP8jDGzXO9vNcYUS+sY/xJPYvF2NMaci/ObOjpmnTFmojHmrDFmzx3eN8aYL137s8sYUz6tY3Sa5rbb4tTcdhO35bWIePUDeAx7I/laIPAO62QCfgOKA1mAn4DHHYp3CNDH9bwPMPgO60U4+Jsm+nsBrwGjXc9bA7M8PN6OwNdOxZhAzDWB8sCeO7zfEFiGvee4CrDV6Zgd+I00t1M/Rs1t98brlrz2+jNqceOQhW7SFJjiej4FaOZQHHeTlN8r7n7MBuoYYxIayCIteNJ/3yQRkXXAhbus0hSYKtYWIK8xpnDaROcZNLfdQnPbjdyV115fqJMooSELizgUSyEROQXg+rfgHdbLaozZZozZYoxJ64RPyu8Vu47Y220uAfnSJLq/S+p/35dcl5tmG2MeTOB9T+JJf7OezJN+J83t1JfecjtZf6+ePNZ3LJOGQxamhrvFew+bKSoiJ40xxYEQY8xuEfktdSJMVFJ+rzT9TRORlFgWATNE5KYxpgf2jKG22yNLPk/6fd1Gc1tzOxHpLbeT9dt6RaEWLxuy8G7xGmPOGGMKi8gp1yWPs3fYxknXv4eMMWuBp7FtNWkhKb/X7XWOG2N8gTzc/ZKPOyUar4icj/NyHDA4DeJKiQwxzKbmtuZ2ItJbbifr7zWjXPr2pCELFwIdXM87AH87azDG+Btj/FzP8wPVgb1pFmHSfq+4+/EyECKu3hIOSDTev7QDNcGOiOXJFgLtXb1EqwCXbl9WVfFobt8bzW1nJS+vne4llwq97IKwRyk3gTPACtfy+4Glf+ltdwB75NrXwXjzAauBX13/BriWBwLjXc+rAbuxPRx3A10ciPNvvxcwEGjiep4V+AE4CIQCxR3+O0gs3s+An12/6RqgtMPxzgBOAZGuv98uQA+gh+t9A3zj2p/d3KHXc3p+aG67LU7NbffF6pa81iFElVJKKQ+WUS59K6WUUl5JC7VSSinlwbRQK6WUUh5MC7VSSinlwbRQK6WUUh7MKwY8Ue5ljLl9WwnYUZeigXOu15XEjrGrlPIymtvpg96epeIxxvTHzu7z378sN9i/lxhHAlNKpYjmtvfSS9/qjowxJYwxe4wxo4H/AQ8aY8LjvN/aGDPe9byQMWaua7KBUNeoO0opD6S57V20UKvEPA5MEJGngRN3We9LYIiIBAItgfFpEZxSKtk0t72EtlGrxPwmImFJWK8uUCrOtLX+xphsInLdfaEppVJAc9tLaKFWibka53kM8adpyxrnuUE7pyjlTTS3vYRe+lZJ5upsctEYU9IY44OdNOG2YOD12y+MMeXSOj6lVPJobns2LdTqXr0HLMfe8nE8zvLXgerGmF3GmL1ANyeCU0olm+a2h9Lbs5RSSikPpmfUSimllAfTQq2UUkp5MC3USimllAfTQq2UUkp5MC3USimllAfTQq2UUkp5MC3USimllAf7f9nn3kcO8lV/AAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAADQCAYAAADI+yJFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd0VNUWwOHfpkZ6kx4EhKjYAOkCUSmKlACigC1IEx88EBUEUVEsiD5EEQtEukhvoUjvIgRUpJdQpEMIPUAgyXl/nBEDJCSQzNyZZH9rZWVm7s2dPSPbfcu5+4gxBqWUUkp5pwxOB6CUUkqpxGmhVkoppbyYFmqllFLKi2mhVkoppbyYFmqllFLKi2mhVkoppbyYFmp1lYjUF5EZbthuVhHZLiIFU3vbSqmkaW77Ni3U6YiI7BORiyJyXkROicgcEfGPt8qnwGeudQuKyHgROSwiZ0TkVxGpepNt9xCRzSJyTkT2ikiPf5YZY6KBEcDb7vpsSqVnbs7tx0VkqWvdffGXaW57hhbq9KexMSYHUAQ4BnwDICKVgdzGmDWu9XIA64BHgHzAaGCOiORIZLsCvAzkBZ4CuohIq3jLfwaCRSRrKn8epZTlrtyOwhbjHoks19x2My3U6ZQx5hIwBSjneqkBsDze8j3GmC+NMUeMMbHGmGFAFuCeRLb3uTHmD2NMjDFmBzATeDTe8oPAKaCaez6RUgrcktthxpixwJ5Elmtuu5kW6nRKRLIBLYF/9rIfBHbcZP3y2GQOT8a2BagFbLlu0Tbg4duJVymVPO7M7ZvQ3HajTE4HoDxuhojEYE9/HQeedL2eBziX0B+ISC5gLPChMeZMMt7jA+xO4MjrXj/neh+lVOrzRG4nRnPbjfSIOv1paozJA2QFugDLRaQw9tRVzutXFpE7gFnAGmNM/6Q2LiJdsNeqG7oGmsSXEzidwviVUglza24nQXPbjbRQp1Oua1PTgFigJrARCIi/jmtwyAzgEPBqUtsUkbZAL6CO67rV9e4D/kph6Eqpm3BHbieD5rYbaaFOp8QKwo7S3gbMBQLjLc+MHZByEXjZGBN33d+XFBEjIiVdz1/A3gJSzxhzw6ATESmGHWG65vplSqnU44bcziAifkBm1+b9RCRLvPU1t91Mr1GnP7NEJBYwwN9AsDFmC4DrPsmqxpi1QA2gETaZT9vxYQA0MMasBPxdf3/I9frHQH5gXbx1fzLGdHI9fh4YncDpcKVU6nBXbtcGlsZ7n4vYUeSPuZ5rbruZGGOcjkF5CRGpD/zHGNM0Geu+C0QYY4YmY92s2NNitY0xx1MeqVLqVmhu+zYt1EoppZQXc/QatYiMEJHjIrI5keUiIoNFJFxENopIRU/HqJS6NZrXSqUupweTjcK2m0xMA6Cs66cj8L0HYlJKpcwoNK+VSjWOFmpjzArg5E1WCQLGGGsNkEdEingmOqXU7dC8Vip1efuo72LAgXjPD7peOxJ/JRHpiN0zJ3v27I/ce++9HgtQqdRgDBw4ABERAL+fMMbc6XRMbpSsvAbNbeX7Ll2C8HCIjr79vPb2Qi0JvHbD6DdXU/lhAJUqVTLr1693d1xKpZqTJ+HZZ+GPP6BnT/j8c/nb6ZjcLFl5DZrbyrctXGhzO1cuiIi4/bx2+hp1Ug5i7+n7R3HgsEOxKJXqdu6EatVg5UoYNQoGDHA6Io/QvFZp3nffQYMGUKIErFuXsm15e6EOBV52jRKtBpwxxtxwekwpX7R4MVStCqdOwZIlEBzsdEQeo3mt0qyYGOjSBTp3toX611/hrrtStk1HT32LyHhsd5sCInIQ6IttU4cx5gds67unsdOvXQBecSZSpVLXDz/YZL73Xpg1C0qVcjqi1KN5rdKr06fhuefsKe8ePaB/f8iYMeXbdbRQG2NaJ7HcAJ09FI5SbhcTA2+8Ad98A08/DePH2+tXaYnmtUqPwsOhUSPYsweGD4e2bVNv294+mEypNOP0aWjZEhYssMX6889TZ29bKeWspUvhmWcgQwZYtAhq107d7Xv7NWql0oTwcKhe3V6LDgmBgQO1SCuVFoSEQP36UKQIhIWlfpEGPaJWyu2WL4fmze3jhQvhscccDUcplQpiY+Gtt+Crr+Cpp2DCBMid2z3vpUfUSrnR8OFQty4ULAhr12qRViotOHsWGje2Rfr11+2AUHcVadBCrZRbxMbCm29C+/ZQpw789huUKeN0VEqplNqzx17GWrjQ3r0xaBBkcvO5aT31rVQqO3sWWreGuXOha1d7PdrdiayUcr+VK+1lrNhYmD8fnnjCM++rR9RKpaK9e6FGDZvE338PX3+tRVqptGDkSHt2LH9+exnLU0UatFArlWpWrYIqVeDQIVuoO3VyOiKlVErFxtrmJW3b2jEma9ZA2bKejUELtVKpYNQou4edN6/d265Tx+mIlFIpde4cNGsG//ufbQk6dy7kyeP5OLRQK5UCsbHw9tvwyiv2/sm1ayEgwOmolFIp9fff8OijtjgPGWJ/nLqMpVfPlLpN58/DCy9AaKg9zT14MGTO7HRUSqmUWr3aHklHR8Mvv0C9es7Go0fUSt2G/fvt3vbs2bZv93ffaZFWKi0YOxYef9z24F+zxvkiDVqolbplv/0GlSvDvn32tFiXLiDidFRKqZSIi4N33oGXX7Y74WvX2tntvIEWaqVuwbhxdm87Rw67t/3kk05HpJRKqagoaNHCTkvZoYO9ayNfPqej+pcWaqWSIS4O3n0XXnwRqlWzzffvu8/pqJRSKXXgANSsCTNn2pagQ4d632UsHUymVBKiouzpsGnToF07ez06Sxano1JKpVRYGAQF2RyfPRsaNHA6ooTpEbVSN3HwoL3tasYM+PJLO6WdFmmlfN+ECRAYCHfcYcedeGuRBi3USiVq3TrbaWzXLnsLVvfuOmhMKV8XFwd9+9p+/JUr26Pq++93Oqqb00KtVAImTrRH0lmz2nsqGzZ0OiKlVEpduACtWkG/frZJ0aJFUKCA01ElTQu1UvEYAx98YJO5UiW7t/3AA05HpZRKqUOH7M73lCm2Jejw4b5zGUsHkynlcvEitGkDkyZBcLAd/Zk1q9NRKaVS6vffoUkTOwVtaCg0auR0RLdGj6iVAg4ftgNLJk+GAQPslHZapJXyfVOmQK1a9par1at9r0iDFmql+OMPO2hs61aYPh169tRBY0r5OmPg44/h2WehQgV7GevBB52O6vZooVbp2tSpttlBhgzw66/2nkqllG+7eNFOmPPee/DSS7B4MRQs6HRUt08LtUqXjIFPPrFtAx9+2O5tP/yw01EppVLq6FHb5nf8eNsSdPRo8PNzOqqU0cFkKt25dAnat7d9u194AX780fcTWSkFGzbYQWORkbaTYLNmTkeUOvSIWqUrx47Zve1x4+z1q7FjtUgrlRbMmGFnvTIGVq1KO0Ua9IhapSN//QWNG9u97SlT4JlnnI5IKfWPK1eusGjRIiIjI6lVqxZ33XVXsv7OGPj8c+jd23YamzEDihRxc7Ae5ugRtYg8JSI7RCRcRHolsLyNiESIyAbXT3sn4lS+LzTU7m3HxcHKlVqk3U1zW92KrVu3EhAQwEcffcTs2bOpWLEivXr1whhz07+Ljra9D3r1gpYtYdmytFekwcEjahHJCHwL1AMOAutEJNQYs/W6VScaY7p4PECVJhgDX3xhE7lSJTuVXVpMZG+iua0Sc/LkSWbOnEl0dDQNGzbE398fYwwtW7bkvffeo23btlfXq1WrFtWqVaNp06YJbuv4cXt6e/Vq2xL03XfT7m2VTh5RVwHCjTF7jDGXgQmA3hyjUs2OHft44om/eftteO45w/LlWqQ9RHNb3SA0NJQyZcrwyy+/sGbNGsqXL8+gQYPYsGEDly5d4pVXXrm6br58+ejRowc//fRTgtvatMn2PvjzT9tJ8L330m6RBmevURcDDsR7fhComsB6z4hIbWAn0N0Yc+D6FUSkI9ARoESJEm4IVfmS2NhY2rXrxbhxzxATU428eb/mwIHJnD8/nTvuuNPp8NIDzW11jdOnT9OmTRsWLFhApUqVADh48CCVK1cmb9685MqVC7mu0ubOnZuoqKgbtjV7tp35KlcuWLHCnilL65w8ok5o/+f6CxKzgJLGmIeARcDohDZkjBlmjKlkjKl05536P+L07v33JzFhQncyZarKhAlw4sR/qVatKq+99prToaUXmtvqGrNnzyYwMPBqkQYoXrw4bdu2ZcuWLRw6dIiwsLCry+Li4hg2bBiN4vX7NAYGDrS3X91zj+19kB6KNDh7RH0Q8I/3vDhwOP4KxpjIeE9DgAEeiEv5sDlz4LPPmpA3byZ++UWoXBkgAx9++CFFixbl9OnT5MmTx+kw0zrNbXWNK1eu4JfAfZB+fn5ER0fz/fff07BhQ1555RXuuusuJkyYQKZMmWjXrh0Aly/Da6/BiBG2SdHo0ZAtm6c/hXOcPKJeB5QVkVIikgVoBYTGX0FE4l9RbAJs82B8yocYA19+aW+/ypx5H+PH73EVaSt79uxkzZqVCxcuOBdk+qG5ra7RoEED5s2bx969e6++dubMGUaNGkVQUBDNmjXj119/JXPmzGzcuJHOnTszf/58/Pz8OHEC6tWzRfr99+1c8empSIODR9TGmBgR6QLMBzICI4wxW0SkH7DeGBMKdBWRJkAMcBJo41S8yntdvgydO9sOY82bQ8GCI5k7N5Z69QZdXWfWrFkUKVKEIjqazO00t9X1ChcuTP/+/alatSrBwcFky5aNMWPGEBQURM2aNQEICAjgk08+uebvtm61O9+HDsHPP9tr0+mRJHWfmq+pVKmSWb9+vdNhKA+JjLT3RC9fDn362Ns0jh8/Ss2aNalcuTINGzZk06ZNjBgxgsmTJ/PYY485HXKSROR3Y0w6ufqWfJrbvm/79u1MmDCB6OhomjRpQvXq1RNdd948e2/0HXfY2yqrJjQc0YekJK+1M5nyWdu22b3tgwfhp59s326we+/r169nxIgRzJ07F39/f9asWcPdd9/tbMBKpXP33nsvH3zwwU3XMQa++Qa6d7fTUs6aBf7+N/2TNE8LtfJJ8+fDc8/ZPt1Ll8L1O+Z58uThjTfecCY4pdRtuXIF/vtfGDoUmja1vfhz5HA6KufppBzKp/yzt/3001CypL1F4yZnz5RSPuLkSXjqKVuke/e2c8Vrkbb0iFr5jCtXoFs3+P57ey/luHGayEqlBTt2QKNGsH8/jBkDL73kdETeRY+olU84dQoaNLBF+u23Yfp0LdJK+arIyEi2bdtGdHQ0CxfagWJnzsCSJVqkE6JH1Mrr7dxpB43t3QujRkFwsNMRKaVuR1RUFK+99hqhoaEULFiQw4eDuHhxAOXKZWDWLHs5S91Ij6iVV1u82O5tnzxp97a1SCvlu/7zn/8QFxfHwoVLuXjxC6KiviAubg5QgytXdjkdntfSQq281g8/wJNPQrFidtCYqy+CUsoHRUZGMnXqVK5cyU716ic5eDCI7t1jGDnyNDExp6hXr552DkyEFmrldWJioGtX29v3ySftfLOlSjkdlVIqMfv372f37t3crIHWm2++SXS0P/Pm9SU2thYFC/YGevLQQ/eTJUsWHnjgAaZOneq5oH3ITa9Ri8hNb0Q1xnyZuuGo9O7o0UsEBh5h585SFCo0jpo1D5ElS1fgxob+6vZpbqvUsHPnTtq0aUN4eDiZM2cmf/78hISEUPW6NmJ//vknc+ZcIDZ2NZcvZ+aFF0YxZMjb3H///Zw/f54aNWqQL18+9u3b58wH8XJJDSbL6fp9D1CZfxvrNwZWuCsolT7t2mWoUOE4Fy7489ZbO2jYsBgDB05g6dJFzJs3jwwZ9ARQKtLcVikSGRlJ5cqVEREyZ87Mk08+SbVq1WjcuDFbt26lQIECV9f96KMjnDw5niJFzhId/QSrVh1nx46HKV68OBMmTOCPP/6gVatW9OvXz8FP5MWMMUn+AAuAnPGe5wTmJedvPf3zyCOPGOV7li0zJmfOyyZjxlNm0aKYq6/HxMSYBx980CxcuNDB6DwLO3GFR/JFc1vdjri4OPPAAw+YQoUKmR07dpj9+/ebHj16mICAAPP888+bQYMGGWOMiYkx5vXXjQFjihX7y5w+bczcuXNNvnz5TN68ec1dd91lXn31VfP888+bRx991MTExCTxzr4rJXmd3EOUEsDleM8vAyVTYT9BKYYPh7p1wc/vDO3ahVCnTsaryzJmzEijRo347bffHIwwTdPcVrds4cKFHD16lKCgIAICAvD39+fzzz+nbNmyXLp0iUOHDnH2rL2t8quvoEOH80RF1SE8/HcaNGjA33//TdOmTTlw4ACLFi3irrvu4pdffiFjxoxJv3k6lNz7qMcCYSIyHTBAM2CM26JSadr+/fs5c+YMAQH38c47mfjySztorEmTBcyfvwrocc3627Zto3Hjxs4Em/ZpbqskrVy5klGjRnH69GlOnz7N6tWrARg+fDiPPPIIHTt2BOCxxx5j4MCB1KnTgerVbQ+EH36AV1/NwdNPh1C/fn0qVaqEMYbff/+dadOmERQU5ORH8w3JPfQGKgLdXD8VbvcQ3t0/enrMex0+fNjUrVvXFChQwNx9dwWTNesiA8Z07WrMlSvGnDt3zhQtWtT8+OOPJiYmxsTExJgRI0aYIkWKmLNnzzodvsfgwVPfRnNbJWHQoEGmRIkSZtCgQaZ+/fqmQIECply5cqZ8+fKmTp06xs/Pz/Tr188sX77cFC9e3BQt2tIUKBBn8uY1ZvHia7d17tw5M336dDN9+nRz/vx5Zz6QQ1KS18mej1pEagJljTEjReROIIcxZm+q7zmkkM5Z652MMdSoUYN69erxwgvv8swzWdi+3ZAtW09WrnyRhx9+GIBNmzbRrl07/v77bwBKlCjB8OHDeeihh5wM36M8PR+15rZKzLFjxwgICGDjxo3kz5+f4sWLs2PHDlq1asXevXsJCgoiIiKCOXPmkCtXLiIjmxATM4RSpYTZs6FsWac/gfdISV4n6xq1iPQF3gZ6u17KDPx0O2+o0qc///yTiIgI6tb9gJo1s3DoEMyfL/TsmZuQkJCr6z344IOEhYWxbt26q7/TU5H2NM1tlZC4uDg++ugjAgICiIqK4vHHHyckJIScOXNSqFAhgoODeeSRRzh58iTTpk3j7NkosmT5iosXvyUwUFizRot0akruYLJmQBMgCsAYc5h/b+9QKklHjx7Fz+9V6tXLQN68sHYt1KkDZcuW5ciRIzesX6JECe666y4HIk13NLfVDT7++GPmz5/Pd999x4MPPsi4ceP4+uuviYqKYuPGjRw9epTChQszduxYevf+lGLF1rJnzzP85z8wdy7kzev0J0hbkluoL7vOsRsAEcnuvpBUWhMbC7/8EsiWLT2oWvUya9dCQIBdNm3aNGpqb1AnaW6ra1y5coXBgwczZswYWrVqxfnz59m+fTtffvklBQoUoFGjRgwYMIDHHnuMtm0/5MMP63H4cHlatlzBF19cIHNmpz9B2pPcUd+TRGQokEdEOgBtgR/dF5ZKK86fhxdegNDQ7DzySBgnT77KwoW9KVy4MKNHj2bTpk0MGzbM6TDTM81tdY3Tp09jjKF06dIAzJw5k2bNmiEi7N+/n0yZMuHv78+rr47m1KnhZMmSk+7dl7Fly1fUrn2EZcuWkUPnoE1VyTqiNsb8D5gCTMV2MnrfGDPYnYEp37d/Pzz6KMyeDd98A2FhlXn//d6MGDGCnj17UqxYMVatWkXu3LmdDjXd0txW18uXLx/ZsmXjzz//BKBcuXJs27aNZs2a8fDDD3P48GHeeusvTp2air9/bgYPXkfZsvv47LPPKFmyJEOHDnX4E6RByRkaDgxIzmve8KO3cHiH1auNKVjQmFy5jJk3z+lofAue7Uymua1uEBISYu6++24zZ84cc+jQITNixAhz5513mlWrVpvevY0BY7JnX2NKlXrEVK9e3QQHB5vChQubunXrmrp16zodvldKSV4n99R3PezI0PgaJPCaUowbB+3a2ekply2D++5zOiJ1E5rb6gbt27cnV65c9OvXj3379lGhQgV+/nkmAwdWZ/p0aNHiJKGhdejSpR9vvGHnd7lw4QIPPfQQefLkcTj6tCep2bNeA/4DlBaRjfEW5QR+dWdgyvfExcH778Mnn0BgIEydCvnzOx2VSojmtkrKc889x3PPPQfAgQPQpAls3GhbgjZpcpYZM6KJjY3FGIOIcObMGS5cuMAdd9zhcORpT1JH1D8DvwD9gV7xXj9njDnptqiUz4mKgpdfhmnT7NH0d99BlixOR6VuQnNbJUtYGAQF2RyfPRsaNIAdO6IpWLAgP/30EyNHjqREiRKsXbuWFi1asGnTJqdDTnNuWqiNMWeAM0BrABEpiJ0YOIeI5DDG7Hd/iMrbHTxoE3nDBvjyS3j9dRBxOip1M5rbKjkmTIDg4FiyZDlJYOD/OHOmArGxz1K2bFmyZ8/Op59+Sv78+Tlx4gRjxoyhS5cu2pffDZLbmayxiOwC9gLLgX3YvXGVzq1bB1WqwK5dEBoK3btrkfYlmtsqIXFx0LcvtG4NGTL8Tq9e02jatCyDBg26ejo8JCSENm3aMGbMGHbv3k3Lli3Zv38/Xbt2dTj6tCe5DU8+BqoBO40xpYA66HWsdEdErv6cOHGCSZOgdm3ImhVWr4aGDZ2OUN0Gze00xBjDhg0bWL58ORcuXLitbVy4AK1aQb9+kCXLOHbuLEmfPq/Svn17Vq5cye7du5k3bx6BgYFs2LCBYsWKsWvXLjp06MCKFSvImVMb26W25BbqK8aYSCCDiGQwxiwFyqf0zUXkKRHZISLhItIrgeVZRWSia/laESmZ0vdUt65nz57IdYfJd975DS1bQqVK9hrWAw84FJxKKc3tNGL37t088sgjtGjRgl69euHv78+oUaNuaRuHDtmd7ylT4Nlnw2jR4hf8/QteXZ4lSxZefPFF5s2bB0CxYsXo06cPQ4YM4fnnnyeLDkxxi+QW6tMikgNYAYwTka+BmJS8sYhkBL7F3gpSDmgtIuWuW60dcMoYUwYYBAxIyXuq2/PFF18AMGvWLC5cMDRtegn4EBjFokVw552OhqdSRnM7DTDG0LRpU4KDg9m1axe//fYbK1eu5J133mHdunWAnWhj7dq1LF26lIsXL96wjd9/t5extm+Po3btgYSG1mbixAm0b9+eyMjIq+sdP35cmxR5WHILdRBwEegOzAN2AykdMVAFCDfG7DHGXAYmuN7n+vcd7Xo8Bagj1x/aKY/IkCEDFSs2IjAQZs7MSunSQ4FXGDz4C6dDUymjuZ0GrFmzhri4OLp27Xr17Fe5cuXo2rUrI0aMYNOmTZQrV462bdvSp08fSpQowYQJE67+/ZQpUKsWZMxoKFCgKTVqRLJr1y5y5MjBggULCAgIYOzYsfz555+MHDmSl156yamPmi4lt4VolDEmFsgGzMJOg5e8iawTVww4EO/5QddrCa5jjInBjlK94c5cEekoIutFZH1EREQKw1IJ8fOrQZUqsHUrTJ8OnTqdBWD/fh0c7Ms0t9OGTZs2cebMGcqXL0+jRo1YsGABAP7+/kRERNC4cWP69OnD5s2bWb16NYsXL6Zbt25s2bKVjz+GZ5+FChXg7bencO+9l/n000+ZPHkyfn5+nD59mvPnz9OpUyeqVKnCwIEDCfhnVh3lEcnqTCYirwL9sHvecYBgk7l0Ct47ob3n6/8HkZx1MMYMA4aBnVw+BTGpBDXnwoWx5M8Pv/4KDz8MIj0B+OabbxyOTaWE5rbv27lzJ++++y7nzp0jJCSEEydO0KFDB/r168fEiRO5++67KVq06DVHwQ899BDBwZ149tlotm2Dl16CYcOgX78/qVmzJsePH+ejjz5i48aNFC5cmNatW5MzZ062bt1K1qxZHfy06VNyT32/BdxvjClpjCltjClljElJIoPdy/aP97w4cDixdUQkE5Ab0GYMHmKM7TJm52v4iwMHClO+vFw9tZYxY0Ynw1OpQ3Pbx/Xv359u3brx0Ucf0blzZ6Kiovjvf/9Lp06dOHHiBBUrVqRYsWtPaBw9CpMnv8a2bRXo3x9GjwY/P7jvvvv49ddfWb58ObVr18bf359MmTKxf/9+WrRowSuvvML8+fMd+qTpV3IL9W7g9sb6J24dUFZESolIFqAVEHrdOqFAsOtxC2CJq7m5crNLl+xe9rvv2mkqW7UaBhy7uvyBBx4gJiZFY46Ud9Dc9nFr1qwhKCiIt956i2HDhrFixQqWLl2Kn58fo0aNol69eixevJjjx48DtjFRlSqG/ftz88YbK+nV69/eBy1atGD37t1MmzaNw4cPc+zYMbp160ZMTAxPPfUUERER5MqVy8FPm04lZ+YOoAKwARgKDP7n53ZnAom33aeBndj/WfRxvdYPaOJ67AdMBsKBMKB0UtvUGXZS7uhRY6pVszPkfPyxMXFxTkeUvuDZ2bM0t31cnTp1zOTJk6957cSJEyZ37tzmzJkzxhhj+vbta8qUKWNat55osmS5bLJmPWaqVHnVXL58+YbtHTx40LRq1cqIiPHz8zNt27Y1ERERZteuXaZIkSLm999/98jnSmtSktdikrETKyJhwCpgE/Y61j9FfnSif+SQSpUqmfXr1zsdhs/66y9o3BgiI2HMGHjmGacjSn9E5HdjTCUPvZfmto+bPn06PXr0YNasWdx3332cPn2aDh06kC9fvqtzQ48aNZrOnfdx4cJ7wHoqV/6E+fNHkTdv3kS3+8cff9C8eXPuvPNO8uXLR1hYGAMGDKBjx44e+mRpS0ryOrnTXMYYY964nTdQviM0FJ5/HvLkgZUroWJFpyNSHqC57eOaNWvGkSNHCAwMJHfu3ERERNCiRQu++uorAObMWUTnztm5cKEvrVrBwIFlefPNbHTv3v2mDVEqVqxIeHg4K1eu5Pz589SqVUunsHRIco+oPwH+xt6+Ef3P68YLZ9nRve5bZwx88QX06mU7jc2cCUWKOB1V+uXhI2rN7TQiOjqaffv2UbBgwatHysePw733buPUqfvo18+OORGBM2fOULJkScLDw8mvc9F6hCeOqJ93/e4S1CieAAAXaklEQVQd77WU3sKhvEB0NHTqBKNGQcuWMHIk6HSy6YrmdhqRNWtW7rnnnqvPN22yl7FOny7Jp5+G07t3mavLcufOTaFChTh27JgWah+QrEJtbLN+lcZEREDz5rBqFXzwAbz/vs58ld5obvuey5cvM2fOHI4cOUKNGjUoX962Zj9z5gyLFy8mc+bMXL5cnzZtspIrFzRr9pWrZWi/q9vYvn07J0+epHRp3R/zBTct1CLyhDFmiYg0T2i5MWaae8JS7rZ5s93bPnrUzjnbsqXTESlP0tz2Tdu3b6dBgwaULFmSgIAAPvvsM2rXrk39+vXp1q0b1apVZ9euxuzenZmyZc+ydGkuoqNbUqNGDTJmzEizZs3YsWMHvXv3pm/fvvj5+Tn9kVQyJHVEHQgsIeHevwbQZPZBc+bYeWZz5IAVK6ByZacjUg7Q3PZBwcHB9OzZk9deew2AS5cuUatWLWbMmMGqVWF880055s2DwMAINm2qSM6cWyhWrDSrVq3is88+o3Xr1hQtWpT//e9/NG3a1OFPo5ItOfdwAaWS85o3/Oi9lomLizNm4EBjRIypUMGYAwecjkglBM/eR6257SP27NljChcubGJiYq55/YUXXjCFCz9gatc2Box57z1jYmONCQoKMmPHjnUoWnW9lOR1cjuTTU3gtSkp301QnnL5MnTsCG++Cc2a2duvihd3OirlBTS3fcTly5fJkiULGTJc+7/tEycKEhExm7Vr4eefoV8/yJAB8uXLx/nz5x2KVqWmpK5R3wvcD+S+7lpWLmxnIeUDIiNt45Lly6FPn38TWaVfmtu+JyAggDvuuIPQ0FCCguysoXPnxrFoUT/i4i4wf/4lAgPtf7qIiAhmzZrFe++952TIKpUkdY36HqARkIdrr2WdAzq4KyiVerZts4PGDh6En36yfbuVQnPb54gIw4cPp2nTpsycGUpkZGtCQx8ne/Yj1K8/hP/+dykdOnTg0qVLfPfdd3Tp0oVSpXRQf1pw00JtjJkJzBSR6saY3zwUk0ol8+fDc8/ZWXGWLoXq1Z2OSHkLzW3f9OijjzJp0nS6dIlj69baVKlyiPnzS5Er1yDmzp3LzJkzyZw5M2PGjKFWrVpOh6tSSXIbnjQTkS3YOWvnAQ8DrxtjfnJbZOq2GQNDhsDrr8ODD9rWoCVKOB2V8lKa2z7iyJEjNGvWjg0b3iE6ujZ+fl/SunUm8uTpCkCjRo1o1KiRw1Eqd0julcr6xpiz2FNlB4EAoIfbolK37coV6NwZunaFRo1sMxMt0uomNLd9RFBQT3buHIsxjzJmDGzZ0pSvvx7EwoULnQ5NuVlyC3Vm1++ngfHGC/sAKzh1Cho0gO+/h7ffhunT7b3SSt2E5rYPGDXqEOvXDyFTpnwsWSK89BKULl2ad955h5CQEKfDU26W3FPfs0RkO/b02H9E5E7gkvvCUrdq5047aGzvXtu3OzjY6YiUj9Dc9nLffQdduxYla9bdhIXlpmTJf5f5+/sTGRnpWGzKM5J1RG2M6QVUByoZY64AF4AgdwamrM2bN7NgwQJOnDiR6DqLF0PVqnDyJCxZokVaJZ/mtveKiYEuXeylrPr148iWrS4XL267Zp3x48fz+OOPOxSh8pSbFmoR6RnvaV1jTCyAMSYK6OrOwNK7iIgIAgMDadiwIQMGDKBMmTL07dv3n85RVw0dCk8+CcWKQVgY1KzpUMDKp2hue7fTp+Hpp+Hbb+Gtt2DWrIwMGPAu9evXZ8iQIcyePZsXX3yRsLAwOnfu7HS4ys2SOqJuFe9x7+uWPZXKsah4XnnlFapUqcKePXtYvHgxO3bsYNq0aUycOBGwe9tdu9opKp98ElavBr1lUt0CzW0vFR4O1arBsmUwfLidKz5jRmjfvj0///wzYWFhDBkyhPvvv5/ffvvt6tzTKu1K6hq1JPI4oecqlRw5coTVq1czZcoUMmbMCEChQoX44IMP+OGHH3jqqVa0bAkLFsAbb8Dnn9tEVuoWaG57oaVLbRfBDBlg0SKoXdtOXwl2DulatWrp/dHpUFJH1CaRxwk9V6nkzJkz5M2b94Yp6IoWLcrhw9moXt1eiw4JgYEDtUir26K57WVCQqB+fShSxF7GKlZsN08++STFihWjWLFi1KtXj/DwcKfDVA5IqlA/LCJnReQc8JDr8T/PH/RAfOlSmTJluHLlCqtXr77m9QED1rBnz88cPw4LF0L79g4FqNICzW0vERsL3bvbSXPq1rWXsYoUuUjdunWpW7cuJ06cIDIykgYNGlC3bl0uXLjgdMjKw5JqIarHag7IlCkTX3/9Nc2bN6dbt27cc889fPHFSdas6ULZssLcuVCmjNNRKl+mue0dzp6FVq3gl19sJ8EvvoBMmeDnn6dzzz330KPHv71n3njjDRYvXszUqVN56aWXHIxaeZrOoeSlmjVrxrx589i37wA9emRgzZr2PPEEhIVl0iKtVBqwZ4/tv79wIfzwAwwaZIs0wN69e6lQocINf1OhQgX27t3r4UiV07RQe7HSpctz8OB37NnTlK5dYf78zOTJ43RUSqmUWrECqlSBI0fs5Dmvvnrt8ooVK7JgwQLi4uKuvmaMYcGCBVSsWNHD0SqnaaH2Unv3Qo0aNom//x6+/vrfvW2llO8aOdJei86fH9auhSeeuHGd+vXrc8cdd/Diiy/y119/sXHjRl5++WUyZcpEgwYNPB+0cpQWai+0apXd2z50yBbqTp2cjkgplVKxsdCmTQRt20KpUn8zdeohypb9d/mCBQto0KAB9913H8HBwXz55ZeULFmSFi1a0Lx5c4oXL868efOu3rKp0g8t1F5m9GioUwfy5rV723XqOB2RUiqlzp2De+/dxujRd/Lww6uoWbM/tWs/yPjx4wH4+eefad++PS+++CJTpkyhQoUKNGrUiJYtW7Jr1y7Cw8Pp378/uXLlcviTKCfoyVQvERcH77wDAwbY4jx5si3WSinf9vffUKdOFLt3l2XAgPP07FkTqEn37l2oVasW9erV45133mHSpElUq1YNgPvvv59MmTLxySefMGnSJGc/gHKcI0fUIpJPRBaKyC7X7wRLkojEisgG10+op+P0lPPnoXlzW6Q7dbK3amiRVr5Ic/taq1fby1gHD2YgOHgSPXv+O+/sAw88wGOPPcbEiROJioq6WqT/0bhxY8LCwjwdsvJCTp367gUsNsaUBRa7nifkojGmvOuniefC85z9++1EGrNmwTff2CntMmdO+u+U8lKa2y5jx8Ljj0POnNC69WDKlNmT4HrZsmUjJiaGw4cPX/P6li1bKF68uCdCVV7OqUIdBIx2PR4NNHUoDkf99pvd2967F+bOtVPaiXZZVr4t3ef2P5exXn7Z3rmxdi107FibkJCQa6ar3bRpE8uWLaNp06a0bduW9u3bc/z4cQC2bdvGW2+9RdeuOpGZcq5QFzLGHAFw/S6YyHp+IrJeRNaISKIJLyIdXeutj4iIcEe8qW7cOLu3nT07rFljZ8BSKg1I17kdFQUtWkD//tChg504J39+qF69Oi+//DIPPvgg3bp1o23btgQGBjJ06FDy5s1L//79CQgI4J577qFkyZI8/vjjdO3aleeee87pj6S8gFw/v3GqbVhkEVA4gUV9gNHGmDzx1j1ljLnhWpaIFDXGHBaR0sASoI4xZvfN3rdSpUpm/fr1KYzefeLi4P334ZNPIDAQpk61iazUP0Tkd2NMJafjSIzmdsIOHIAmTWDjRvjySzsN7fVnyLZu3crs2bPx8/OjRYsWFC1a9Jrl586dIyIiguLFi5MlSxYPRq/cLSV57bZR38aYuoktE5FjIlLEGHNERIoAxxPZxmHX7z0isgyoANw0mb1ZVJQ9HTZtmp1Q49tvQXNR+RrN7RuFhUFQkM3x2bMhsZ4k5cqVo1y5coluJ2fOnOTMmdNNUSpf5dSp71Ag2PU4GJh5/QoikldEsroeFwAeBbZ6LMJUdvCgnVt2xgy7tz1smBZplSalu9yeMMGeHbvjDjvuRBuHqdTmVKH+DKgnIruAeq7niEglEfnRtc59wHoR+QtYCnxmjPHJZF63zg4a27ULQkPtlHY6aEylUekmt+PioG9faN0aKle2R9X33+90VCotcqThiTEmErih55YxZj3Q3vV4NWlgXtxJkyA4GAoXtgNLHnjA6YiUcp/0ktsXLkCbNrYx0Suv2H78WbM6HZVKq7SFqJsYAx9+CC1bQqVKdm9bi7RSvu/QIXsZa8oUO3/08OFapJV7aQtRN7h40e5lT5xoj6aHDtVEViot+P13O7L77FmYORMaN3Y6IpUe6BH1bbp06RJXrly54fUjR+zAkkmTbEvQkSO1SCuVFkyZArVq2elmf/1Vi7TyHC3Ut2jz5s3UrVuX3LlzkydPHoKDgzl58iQAf/5pB5Vs3QrTp0PPnjpoTClfZwx8/DE8+yyUL28vYz30kNNRqfREC/UtiIiIoG7durRo0YKzZ89y4MABcubMSePGjZk61VCzJmTIYPe2g4KcjlYplVIXL8ILL8B778GLL8KSJVCokNNRqfRGr1HfgpEjR9KwYUM6deoEQNasWRk8+BsKFfqaFi2EatXskXThhHo2KaV8ytGj0LSp7dX96afQq5eeIVPO0EJ9C8LDw6lcufLV55cuQfv2wokTr1O9+m6WLLkbPz8HA1RKpYoNG+ygschI20mwWTOnI1LpmZ76vgUPPfQQS5YsAeDYMTupxrhxkDv3FwwZclaLtFJpwIwZ8Oij9tr0qlVapJXztFDfgpdffpn169fTrt1gKlaM4a+/4qhS5XMCA1dRsWIFp8NTSqWAMfDZZ9C8ue15EBYGFTStlRfQQn0LcuXKxbvvhjFmTEeOHj1O7tyNqVfvLBMnTnQ6NKVUCkRH254HvXvbJkXLlkGRIk5HpZSl16iTyRjbhahXrwJUqgQzZxalSJE5ToellEqh48ft6e3Vq6FfP3j3XR00pryLFupkiI6GTp1g1Ci7tz1ypJ0pRynl2zZtso1Ljh+3TYqefdbpiJS6kZ76TkJEBNSta4v0Bx/A+PFapJVKC2bPhho14MoVWLFCi7TyXlqob2LzZjs95fr1ds7Zvn31lJhSvs4YGDjQ3n4VEGAHjVWq5HRUSiVOC3Ui5s61e9vR0XZvu2VLpyNSSqXU5cvQvj289RY88wysXAnFijkdlVI3p4X6OsbAoEH2ulWZMnZvO16PE6WUjzpxAurVgxEjbEvQiRMhWzano1IqaTqYLJ7Ll6FLFwgJsfdSjhkD2bM7HZVSKqW2brU734cO2SZFzz/vdERKJZ8eUbtERkL9+rZI9+kDkydrkVYqLZg3D6pXh6goWL5ci7TyPVqogW3boGpVWLMGfvrJTmmXQb8ZpXyaMTB4MDRsCKVK2ctYVas6HZVSty7dl6P586FaNTh3DpYutVPaKaV825Ur8Npr0K2bHd29ahWUKOF0VErdnnRbqI2Bb76Bp5+2e9vr1tnTY0op33byJDz1FAwdaluCTp0KOXI4HZVSty9dDia7csXuaX//vd3bHjdOE1mptGDHDmjUCPbvt4NBX3rJ6YiUSrl0V6hPnbIdiBYvhrffthPC6/VopXzfwoU2t7NkgSVL7FSVSqUF6apE7dxpr0evWGFbgn72mRZppdKC776DBg3A398OGtMirdKSdFOmFi+2Iz5PnrR728HBTkeklEqpmBjb+6BzZ1uoV6+GkiWdjkqp1JUuCvXQofDkk7ZVYFgY1KzpdERKqZQ6fdoOBv32W9sSdMYMyJnT6aiUSn1p+hp1TAy8+aa9l/Lpp+3MV7lyOR2VUiqlwsPtoLE9e2D4cGjb1umIlHKfNFuoz5yxE2nMnw9vvAGffw4ZMzodlVIqpZYutRNqZMhgB5AFBjodkVLu5cipbxF5VkS2iEiciCQ6wZyIPCUiO0QkXER6JXf7u3fbe6IXL7YtQQcO1CKtlCe4O7dDQmyr38KFYe1aLdIqfXDqGvVmoDmwIrEVRCQj8C3QACgHtBaRcklt+Nw5O4f0sWN2b7t9+9QKWSmVDG7L7QMHoGNHqFMHfvsN7r47tUJWyrs5UqiNMduMMTuSWK0KEG6M2WOMuQxMAIKS2vbOnVCwoN3bfuyxVAhWKZVs7szt48dto6LZsyF37tSIVinf4M3XqIsBB+I9Pwgk2FJfRDoCHV1Po7dvl81ly7o5utRVADjhdBC3QON1r3ucDsDNbju3v/5aNn/9tZujSz2+9u9O43Wv285rtxVqEVkEFE5gUR9jzMzkbCKB10xCKxpjhgHDXO+73hiT6LUxb+RrMWu87iUi652O4WY0t5NH43UvX4z3dv/WbYXaGFM3hZs4CPjHe14cOJzCbSqlUkhzWynP8uaGJ+uAsiJSSkSyAK2AUIdjUkqlnOa2UrfAqduzmonIQaA6MEdE5rteLyoicwGMMTFAF2A+sA2YZIzZkozND3NT2O7kazFrvO7la/Fepbl9DY3XvdJNvGJMgpeGlFJKKeUFvPnUt1JKKZXuaaFWSimlvJjPF2p3tyxMbSKST0QWisgu1++8iawXKyIbXD8eH2iT1PclIllFZKJr+VoRKenpGK+LJ6l424hIRLzv1NGedSIyQkSOi8jmRJaLiAx2fZ6NIlLR0zE6TXPbbXFqbruJ2/LaGOPTP8B92BvJlwGVElknI7AbKA1kAf4CyjkU7+dAL9fjXsCARNY77+B3muT3BfwH+MH1uBUw0cvjbQMMcSrGBGKuDVQENiey/GngF+w9x9WAtU7H7MB3pLmd+jFqbrs3Xrfktc8fURs3tix0kyBgtOvxaKCpQ3HcTHK+r/ifYwpQR0QSamThCd703zdZjDErgJM3WSUIGGOsNUAeESnimei8g+a2W2huu5G78trnC3UyJdSysJhDsRQyxhwBcP0umMh6fiKyXkTWiIinEz4539fVdYy93eYMkN8j0d0ouf99n3GdbpoiIv4JLPcm3vRv1pt50/ekuZ360lpu39a/V2/u9X2VeLBlYWq4Wby3sJkSxpjDIlIaWCIim4wxu1MnwiQl5/vy6HeahOTEMgsYb4yJFpFO2COGJ9we2e3zpu/XbTS3NbeTkNZy+7a+W58o1MbHWhbeLF4ROSYiRYwxR1ynPI4nso3Drt97RGQZUAF7rcYTkvN9/bPOQRHJBOTm5qd83CnJeI0xkfGehgADPBBXSqSLNpua25rbSUhruX1b/17Ty6lvb2pZGAoEux4HAzccNYhIXhHJ6npcAHgU2OqxCJP3fcX/HC2AJcY1WsIBScZ73XWgJtiOWN4sFHjZNUq0GnDmn9Oq6hqa27dGc9tZt5fXTo+SS4VRds2weynRwDFgvuv1osDc60bb7cTuufZxMN78wGJgl+t3PtfrlYAfXY9rAJuwIxw3Ae0ciPOG7wvoBzRxPfYDJgPhQBhQ2uF/B0nF2x/Y4vpOlwL3OhzveOAIcMX177cd0Ano5FouwLeuz7OJREY9p+UfzW23xam57b5Y3ZLX2kJUKaWU8mLp5dS3Ukop5ZO0UCullFJeTAu1Ukop5cW0UCullFJeTAu1Ukop5cV8ouGJci8R+ee2ErBdl2KBCNfzKsb22FVK+RjN7bRBb89S1xCRD7Cz+/zvutcF++8lzpHAlFIporntu/TUt0qUiJQRkc0i8gPwB+AvIqfjLW8lIj+6HhcSkWmuyQbCXF13lFJeSHPbt2ihVkkpBww3xlQADt1kvcHA58aYSsBzwI+eCE4pdds0t32EXqNWSdltjFmXjPXqAvfEm7Y2r4jcYYy56L7QlFIpoLntI7RQq6RExXscx7XTtPnFeyzo4BSlfInmto/QU98q2VyDTU6JSFkRyYCdNOEfi4DO/zwRkfKejk8pdXs0t72bFmp1q94G5mFv+TgY7/XOwKMislFEtgIdnAhOKXXbNLe9lN6epZRSSnkxPaJWSimlvJgWaqWUUsqLaaFWSimlvJgWaqWUUsqLaaFWSimlvJgWaqWUUsqLaaFWSimlvNj/AYuuAdQVtGYsAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "for t in range(1, n_timepoints):\n", + " B_t, B_tau = model.adjacency_matrices_[t]\n", + " plt.figure(figsize=(7, 3))\n", + "\n", + " plt.subplot(1,2,1)\n", + " plt.plot([-1, 1],[-1, 1], marker=\"\", color=\"blue\", label=\"support\")\n", + " plt.scatter(B_t_true[t], B_t, facecolors='none', edgecolors='black')\n", + " plt.xlim(-1, 1)\n", + " plt.ylim(-1, 1)\n", + " plt.xlabel('True')\n", + " plt.ylabel('Estimated')\n", + " plt.title(f'B({t},{t})')\n", + "\n", + " plt.subplot(1,2,2)\n", + " plt.plot([-1, 1],[-1, 1], marker=\"\", color=\"blue\", label=\"support\")\n", + " plt.scatter(B_tau_true[t], B_tau, facecolors='none', edgecolors='black')\n", + " plt.xlim(-1, 1)\n", + " plt.ylim(-1, 1)\n", + " plt.xlabel('True')\n", + " plt.ylabel('Estimated')\n", + " plt.title(f'B({t},{t-1})')\n", + "\n", + " plt.tight_layout()\n", + " plt.show()" + ] + }, + { + "cell_type": "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.3" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/lingam/__init__.py b/lingam/__init__.py index 473e1a2..32b743e 100644 --- a/lingam/__init__.py +++ b/lingam/__init__.py @@ -10,7 +10,10 @@ from .causal_effect import CausalEffect from .var_lingam import VARLiNGAM from .varma_lingam import VARMALiNGAM +from .longitudinal_lingam import LongitudinalLiNGAM +from .bootstrap import LongitudinalBootstrapResult -__all__ = ['ICALiNGAM', 'DirectLiNGAM', 'BootstrapResult', 'MultiGroupDirectLiNGAM', 'CausalEffect', 'VARLiNGAM', 'VARMALiNGAM'] +__all__ = ['ICALiNGAM', 'DirectLiNGAM', 'BootstrapResult', 'MultiGroupDirectLiNGAM', + 'CausalEffect', 'VARLiNGAM', 'VARMALiNGAM', 'LongitudinalLiNGAM', 'LongitudinalBootstrapResult'] -__version__ = '1.2.1' +__version__ = '1.3' diff --git a/lingam/bootstrap.py b/lingam/bootstrap.py index c343933..3b842ec 100644 --- a/lingam/bootstrap.py +++ b/lingam/bootstrap.py @@ -107,14 +107,16 @@ def get_causal_direction_counts(self, n_directions=None, min_causal_effect=None, min_causal_effect = 0.0 else: if not 0.0 < min_causal_effect: - raise ValueError('min_causal_effect must be an value greater than 0.') + raise ValueError( + 'min_causal_effect must be an value greater than 0.') # Count causal directions directions = [] for am in self._adjacency_matrices: direction = np.array(np.where(np.abs(am) > min_causal_effect)) if split_by_causal_effect_sign: - signs = np.array([np.sign(am[i][j]) for i, j in direction.T]).astype('int64').T + signs = np.array([np.sign(am[i][j]) + for i, j in direction.T]).astype('int64').T direction = np.vstack([direction, signs]) directions.append(direction.T) directions = np.concatenate(directions) @@ -177,7 +179,8 @@ def get_directed_acyclic_graph_counts(self, n_dags=None, min_causal_effect=None, min_causal_effect = 0.0 else: if not 0.0 < min_causal_effect: - raise ValueError('min_causal_effect must be an value greater than 0.') + raise ValueError( + 'min_causal_effect must be an value greater than 0.') # Count directed acyclic graphs dags = [] @@ -231,7 +234,8 @@ def get_probabilities(self, min_causal_effect=None): min_causal_effect = 0.0 else: if not 0.0 < min_causal_effect: - raise ValueError('min_causal_effect must be an value greater than 0.') + raise ValueError( + 'min_causal_effect must be an value greater than 0.') shape = self._adjacency_matrices[0].shape bp = np.zeros(shape) @@ -244,3 +248,218 @@ def get_probabilities(self, min_causal_effect=None): else: return np.hsplit(bp, int(shape[1]/shape[0])) + +class LongitudinalBootstrapResult(object): + """The result of bootstrapping for LongitudinalLiNGAM.""" + + def __init__(self, adjacency_matrices, n_timepoints): + """Construct a BootstrapResult. + + Parameters + ---------- + adjacency_matrices : array-like, shape (n_sampling) + The adjacency matrix list by bootstrapping. + """ + self._adjacency_matrices = adjacency_matrices + self._n_timepoints = n_timepoints + + @property + def adjacency_matrices_(self): + """The adjacency matrix list by bootstrapping. + + Returns + ------- + adjacency_matrices_ : array-like, shape (n_sampling) + The adjacency matrix list, where ``n_sampling`` is + the number of bootstrap sampling. + """ + return self._adjacency_matrices + + def get_causal_direction_counts(self, n_directions=None, min_causal_effect=None, split_by_causal_effect_sign=False): + """Get causal direction count as a result of bootstrapping. + + Parameters + ---------- + n_directions : int, optional (default=None) + If int, then The top ``n_directions`` items are included in the result + min_causal_effect : float, optional (default=None) + Threshold for detecting causal direction. + If float, then causal directions with absolute values of causal effects less than ``min_causal_effect`` are excluded. + split_by_causal_effect_sign : boolean, optional (default=False) + If True, then causal directions are split depending on the sign of the causal effect. + + Returns + ------- + causal_direction_counts : dict + List of causal directions sorted by count in descending order. + The dictionary has the following format:: + + {'from': [n_directions], 'to': [n_directions], 'count': [n_directions]} + + where ``n_directions`` is the number of causal directions. + """ + # Check parameters + if isinstance(n_directions, (numbers.Integral, np.integer)): + if not 0 < n_directions: + raise ValueError( + 'n_directions must be an integer greater than 0') + elif n_directions is None: + pass + else: + raise ValueError('n_directions must be an integer greater than 0') + + if min_causal_effect is None: + min_causal_effect = 0.0 + else: + if not 0.0 < min_causal_effect: + raise ValueError( + 'min_causal_effect must be an value greater than 0.') + + # Count causal directions + cdc_list = [] + for t in range(self._n_timepoints): + + directions = [] + for m in self._adjacency_matrices: + am = np.concatenate([*m[t]], axis=1) + direction = np.array(np.where(np.abs(am) > min_causal_effect)) + if split_by_causal_effect_sign: + signs = np.array([np.sign(am[i][j]) + for i, j in direction.T]).astype('int64').T + direction = np.vstack([direction, signs]) + directions.append(direction.T) + directions = np.concatenate(directions) + + if len(directions) == 0: + cdc = {'from': [], 'to': [], 'count': []} + if split_by_causal_effect_sign: + cdc['sign'] = [] + cdc_list.append(cdc) + continue + + directions, counts = np.unique( + directions, axis=0, return_counts=True) + sort_order = np.argsort(-counts) + sort_order = sort_order[:n_directions] if n_directions is not None else sort_order + counts = counts[sort_order] + directions = directions[sort_order] + + cdc = { + 'from': directions[:, 1].tolist(), + 'to': directions[:, 0].tolist(), + 'count': counts.tolist() + } + if split_by_causal_effect_sign: + cdc['sign'] = directions[:, 2].tolist() + + cdc_list.append(cdc) + + return cdc_list + + def get_directed_acyclic_graph_counts(self, n_dags=None, min_causal_effect=None, split_by_causal_effect_sign=False): + """Get DAGs count as a result of bootstrapping. + + Parameters + ---------- + n_dags : int, optional (default=None) + If int, then The top ``n_dags`` items are included in the result + min_causal_effect : float, optional (default=None) + Threshold for detecting causal direction. + If float, then causal directions with absolute values of causal effects less than ``min_causal_effect`` are excluded. + split_by_causal_effect_sign : boolean, optional (default=False) + If True, then causal directions are split depending on the sign of the causal effect. + + Returns + ------- + directed_acyclic_graph_counts : dict + List of directed acyclic graphs sorted by count in descending order. + The dictionary has the following format:: + + {'dag': [n_dags], 'count': [n_dags]}. + + where ``n_dags`` is the number of directed acyclic graphs. + """ + # Check parameters + if isinstance(n_dags, (numbers.Integral, np.integer)): + if not 0 < n_dags: + raise ValueError('n_dags must be an integer greater than 0') + elif n_dags is None: + pass + else: + raise ValueError('n_dags must be an integer greater than 0') + + if min_causal_effect is None: + min_causal_effect = 0.0 + else: + if not 0.0 < min_causal_effect: + raise ValueError( + 'min_causal_effect must be an value greater than 0.') + + # Count directed acyclic graphs + dagc_list = [] + for t in range(self._n_timepoints): + + dags = [] + for m in self._adjacency_matrices: + am = np.concatenate([*m[t]], axis=1) + + dag = np.abs(am) > min_causal_effect + if split_by_causal_effect_sign: + direction = np.array(np.where(dag)) + signs = np.zeros_like(dag).astype('int64') + for i, j in direction.T: + signs[i][j] = np.sign(am[i][j]).astype('int64') + dag = signs + dags.append(dag) + + dags, counts = np.unique(dags, axis=0, return_counts=True) + sort_order = np.argsort(-counts) + sort_order = sort_order[:n_dags] if n_dags is not None else sort_order + counts = counts[sort_order] + dags = dags[sort_order] + + if split_by_causal_effect_sign: + dags = [{ + 'from': np.where(dag)[1].tolist(), + 'to': np.where(dag)[0].tolist(), + 'sign': [dag[i][j] for i, j in np.array(np.where(dag)).T]} for dag in dags] + else: + dags = [{ + 'from': np.where(dag)[1].tolist(), + 'to': np.where(dag)[0].tolist()} for dag in dags] + + dagc_list.append({ + 'dag': dags, + 'count': counts.tolist() + }) + + return dagc_list + + def get_probabilities(self, min_causal_effect=None): + """Get bootstrap probability. + + Parameters + ---------- + min_causal_effect : float, optional (default=None) + Threshold for detecting causal direction. + If float, then causal directions with absolute values of causal effects less than ``min_causal_effect`` are excluded. + + Returns + ------- + probabilities : array-like + List of bootstrap probability matrix. + """ + # check parameters + if min_causal_effect is None: + min_causal_effect = 0.0 + else: + if not 0.0 < min_causal_effect: + raise ValueError( + 'min_causal_effect must be an value greater than 0.') + + prob = np.zeros(self._adjacency_matrices.shape) + for adj_mat in self._adjacency_matrices: + prob += np.where(np.abs(adj_mat) > min_causal_effect, 1, 0) + prob = prob/len(self._adjacency_matrices) + + return prob diff --git a/lingam/direct_lingam.py b/lingam/direct_lingam.py index d0304d1..6c85767 100644 --- a/lingam/direct_lingam.py +++ b/lingam/direct_lingam.py @@ -5,6 +5,7 @@ import numpy as np from sklearn.utils import check_array +from sklearn.preprocessing import scale from .base import _BaseLiNGAM @@ -20,7 +21,7 @@ class DirectLiNGAM(_BaseLiNGAM): Journal of Machine Learning Research 14:111-152, 2013. """ - def __init__(self, random_state=None, prior_knowledge=None): + def __init__(self, random_state=None, prior_knowledge=None, measure='pwling'): """Construct a DirectLiNGAM model. Parameters @@ -29,9 +30,18 @@ def __init__(self, random_state=None, prior_knowledge=None): ``random_state`` is the seed used by the random number generator. prior_knowledge : array-like, shape (n_features, n_features), optional (default=None) Prior knowledge used for causal discovery, where ``n_features`` is the number of features. + + The elements of prior knowledge matrix are defined as follows [1]_: + + * ``0`` : :math:`x_i` does not have a directed path to :math:`x_j` + * ``1`` : :math:`x_i` has a directed path to :math:`x_j` + * ``-1`` : No prior knowledge is available to know if either of the two cases above (0 or 1) is true. + measure : {'pwling', 'kernel'}, default='pwling' + Measure to evaluate independence: 'pwling' [2]_ or 'kernel' [1]_. """ super().__init__(random_state) self._prior_knowledge = prior_knowledge + self._measure = measure def fit(self, X): """Fit the model to X. @@ -55,7 +65,8 @@ def fit(self, X): self._Aknw = check_array(self._prior_knowledge) self._Aknw = np.where(self._Aknw < 0, np.nan, self._Aknw) if (n_features, n_features) != self._Aknw.shape: - raise ValueError('The shape of prior knowledge must be (n_features, n_features)') + raise ValueError( + 'The shape of prior knowledge must be (n_features, n_features)') else: self._Aknw = None @@ -63,8 +74,14 @@ def fit(self, X): U = np.arange(n_features) K = [] X_ = np.copy(X) + if self._measure == 'kernel': + X_ = scale(X_) + for _ in range(n_features): - m = self._search_causal_order(X_, U) + if self._measure == 'kernel': + m = self._search_causal_order_kernel(X_, U) + else: + m = self._search_causal_order(X_, U) for i in U: if i != m: X_[:, i] = self._residual(X_[:, i], X_[:, m]) @@ -147,3 +164,47 @@ def _search_causal_order(self, X, U): M += np.min([0, self._diff_mutual_info(xi_std, xj_std, ri_j, rj_i)])**2 M_list.append(-1.0 * M) return Uc[np.argmax(M_list)] + + def _mutual_information(self, x1, x2, param): + """Calculate the mutual informations.""" + kappa, sigma = param + n = len(x1) + X1 = np.tile(x1, (n, 1)) + K1 = np.exp(-1/(2*sigma**2) * (X1**2 + X1.T**2 - 2*X1*X1.T)) + X2 = np.tile(x2, (n, 1)) + K2 = np.exp(-1/(2*sigma**2) * (X2**2 + X2.T**2 - 2*X2*X2.T)) + + tmp1 = K1 + n*kappa*np.identity(n)/2 + tmp2 = K2 + n*kappa*np.identity(n)/2 + K_kappa = np.r_[np.c_[tmp1 @ tmp1, K1 @ K2], + np.c_[K2 @ K1, tmp2 @ tmp2]] + D_kappa = np.r_[np.c_[tmp1 @ tmp1, np.zeros([n, n])], + np.c_[np.zeros([n, n]), tmp2 @ tmp2]] + + sigma_K = np.linalg.svd(K_kappa, compute_uv=False) + sigma_D = np.linalg.svd(D_kappa, compute_uv=False) + + return (-1/2)*(np.sum(np.log(sigma_K)) - np.sum(np.log(sigma_D))) + + def _search_causal_order_kernel(self, X, U): + """Search the causal ordering by kernel method.""" + Uc, Vj = self._search_candidate(U) + if len(Uc) == 1: + return Uc[0] + + if X.shape[0] > 1000: + param = [2e-3, 0.5] + else: + param = [2e-2, 1.0] + + Tkernels = [] + for j in Uc: + Tkernel = 0 + for i in U: + if i != j: + ri_j = X[:, i] if j in Vj and i in Uc else self._residual( + X[:, i], X[:, j]) + Tkernel += self._mutual_information(X[:, j], ri_j, param) + Tkernels.append(Tkernel) + + return Uc[np.argmin(Tkernels)] diff --git a/lingam/longitudinal_lingam.py b/lingam/longitudinal_lingam.py new file mode 100644 index 0000000..d459970 --- /dev/null +++ b/lingam/longitudinal_lingam.py @@ -0,0 +1,212 @@ +""" +Python implementation of the LiNGAM algorithms. +The LiNGAM Project: https://sites.google.com/site/sshimizu06/lingam +""" +import numbers +import numpy as np +from sklearn.utils import check_array, resample +from sklearn.linear_model import LinearRegression + +from .direct_lingam import DirectLiNGAM +from .bootstrap import LongitudinalBootstrapResult + + +class LongitudinalLiNGAM(): + """Implementation of Longitudinal LiNGAM algorithm [1]_ + + References + ---------- + .. [1] K. Kadowaki, S. Shimizu, and T. Washio. Estimation of causal structures in longitudinal data using non-Gaussianity. In Proc. 23rd IEEE International Workshop on Machine Learning for Signal Processing (MLSP2013), pp. 1--6, Southampton, United Kingdom, 2013. + """ + + def __init__(self, n_lags=1, measure='pwling', random_state=None): + """Construct a model. + + Parameters + ---------- + n_lags : int, optional (default=1) + Number of lags. + measure : {'pwling', 'kernel'}, default='pwling' + Measure to evaluate independence : 'pwling' or 'kernel'. + random_state : int, optional (default=None) + ``random_state`` is the seed used by the random number generator. + """ + self._n_lags = n_lags + self._measure = measure + self._random_state = random_state + self._causal_orders = None + self._adjacency_matrices = None + + def fit(self, X_list): + """Fit the model to datasets. + + Parameters + ---------- + X_list : list, shape [X, ...] + Longitudinal multiple datasets for training, where ``X`` is an dataset. + The shape of ``X`` is (n_samples, n_features), + where ``n_samples`` is the number of samples and ``n_features`` is the number of features. + + Returns + ------- + self : object + Returns the instance itself. + """ + # Check parameters + if not isinstance(X_list, (list, np.ndarray)): + raise ValueError('X_list must be a array-like.') + + if len(X_list) < 2: + raise ValueError('X_list must be a list containing at least two items') + + n_samples = check_array(X_list[0]).shape[0] + n_features = check_array(X_list[0]).shape[1] + X_t = [] + for X in X_list: + X = check_array(X) + if X.shape != (n_samples, n_features): + raise ValueError('X_list must be a list with the same shape') + X_t.append(X.T) + + self._T = len(X_t) # Number of time points + self._n = n_samples # Number of samples + self._p = n_features # Number of features + + M_tau, N_t = self._compute_residuals(X_t) + B_t, causal_orders = self._estimate_instantaneous_effects(N_t) + B_tau = self._estimate_lagged_effects(B_t, M_tau) + + # output B(t,t), B(t,t-τ) + self._adjacency_matrices = np.empty((self._T, 1+self._n_lags, self._p, self._p)) + self._adjacency_matrices[:, :] = np.nan + for t in range(1, self._T): + self._adjacency_matrices[t, 0] = B_t[t] + for l in range(self._n_lags): + if t-l == 0: + continue + self._adjacency_matrices[t, l+1] = B_tau[t, l] + self._causal_orders = causal_orders + return self + + def bootstrap(self, X_list, n_sampling): + """Evaluate the statistical reliability of DAG based on the bootstrapping. + + Parameters + ---------- + X_list : array-like, shape (X, ...) + Longitudinal multiple datasets for training, where ``X`` is an dataset. + The shape of ''X'' is (n_samples, n_features), + where ``n_samples`` is the number of samples and ``n_features`` is the number of features. + n_sampling : int + Number of bootstrapping samples. + + Returns + ------- + results : array-like, shape (BootstrapResult, ...) + Returns the results of bootstrapping for multiple datasets. + """ + # Check parameters + if not isinstance(X_list, (list, np.ndarray)): + raise ValueError('X_list must be a array-like.') + + if len(X_list) < 2: + raise ValueError('X_list must be a list containing at least two items') + + n_samples = check_array(X_list[0]).shape[0] + n_features = check_array(X_list[0]).shape[1] + X_t = [] + for X in X_list: + X = check_array(X) + if X.shape != (n_samples, n_features): + raise ValueError('X_list must be a list with the same shape') + X_t.append(X) + + self._T = len(X_t) # Number of time points + self._n = n_samples # Number of samples + self._p = n_features # Number of features + + adjacency_matrices = np.zeros( + (n_sampling, self._T, 1+self._n_lags, self._p, self._p)) + for i in range(n_sampling): + print('sampling:', i) + resampled_X_t = np.empty((self._T, self._n, self._p)) + indices = np.random.randint(0, self._n, size=(self._n,)) + for t in range(self._T): + resampled_X_t[t] = X_t[t][indices, :] + + model = self.fit(resampled_X_t) + adjacency_matrices[i] = model.adjacency_matrices_ + return LongitudinalBootstrapResult(adjacency_matrices, self._T) + + def _compute_residuals(self, X_t): + """Compute residuals N(t)""" + M_tau = np.zeros((self._T, self._n_lags, self._p, self._p)) + N_t = np.zeros((self._T, self._p, self._n)) + + for t in range(1, self._T): + # predictors + X_predictors = np.zeros((self._n, self._p*(1+self._n_lags))) + for tau in range(self._n_lags): + pos = self._p * tau + X_predictors[:, pos:pos+self._p] = X_t[t-(tau+1)].T + + # estimate M(t,t-τ) by regression + X_target = X_t[t].T + for i in range(self._p): + reg = LinearRegression() + reg.fit(X_predictors, X_target[:, i]) + for tau in range(self._n_lags): + pos = self._p * tau + M_tau[t, tau, i] = reg.coef_[pos:pos+self._p] + + # Compute N(t) + N_t[t] = X_t[t] + for tau in range(self._n_lags): + N_t[t] = N_t[t] - np.dot(M_tau[t, tau], X_t[t-(tau+1)]) + + return M_tau, N_t + + def _estimate_instantaneous_effects(self, N_t): + """Estimate instantaneous effects B(t,t) by applying LiNGAM""" + causal_orders = [] + B_t = np.zeros((self._T, self._p, self._p)) + for t in range(1, self._T): + model = DirectLiNGAM(measure=self._measure) + model.fit(N_t[t].T) + causal_orders.append(model.causal_order_) + B_t[t] = model.adjacency_matrix_ + return B_t, causal_orders + + def _estimate_lagged_effects(self, B_t, M_tau): + """Estimate lagged effects B(t,t-τ)""" + B_tau = np.zeros((self._T, self._n_lags, self._p, self._p)) + for t in range(self._T): + for tau in range(self._n_lags): + B_tau[t, tau] = np.dot(np.eye(self._p) - B_t[t], M_tau[t, tau]) + return B_tau + + @property + def causal_orders_(self): + """Estimated causal ordering. + + Returns + ------- + causal_order_ : array-like, shape (causal_order, ...) + The causal order of fitted models for B(t,t). + The shape of causal_order is (n_features), + where ``n_features`` is the number of features. + """ + return self._causal_orders + + @property + def adjacency_matrices_(self): + """Estimated adjacency matrices. + + Returns + ------- + adjacency_matrices_ : array-like, shape ((B(t,t), B(t,t-1), ..., B(t,t-τ)), ...) + The list of adjacency matrix B(t,t) and B(t,t-τ) for longitudinal datasets. + The shape of B(t,t) and B(t,t-τ) is (n_features, n_features), where + ``n_features`` is the number of features. + """ + return self._adjacency_matrices diff --git a/lingam/var_lingam.py b/lingam/var_lingam.py index 9acfa21..0e1f6b5 100644 --- a/lingam/var_lingam.py +++ b/lingam/var_lingam.py @@ -8,10 +8,9 @@ from sklearn.linear_model import LassoLarsIC, LinearRegression from statsmodels.tsa.vector_ar.var_model import VAR -from lingam import DirectLiNGAM - from .base import _BaseLiNGAM from .bootstrap import BootstrapResult +from .direct_lingam import DirectLiNGAM class VARLiNGAM: @@ -38,8 +37,8 @@ def __init__(self, lags=1, criterion='bic', prune=False, ar_coefs=None, lingam_m ar_coefs : array-like, optional (default=None) Coefficients of AR model. Estimating AR model is skipped if specified ``ar_coefs``. Shape must be (``lags``, n_features, n_features). - lingam_model : constructor - Constructor of a LiNGAM algorithm which inherits _BaseLiNGAM. + lingam_model : lingam object inherits 'lingam._BaseLiNGAM', optional (default=None) + LiNGAM model for causal discovery. If None, DirectLiNGAM algorithm is selected. random_state : int, optional (default=None) ``random_state`` is the seed used by the random number generator. """ @@ -71,7 +70,7 @@ def fit(self, X): lingam_model = self._lingam_model if lingam_model is None: - lingam_model = DirectLiNGAM + lingam_model = DirectLiNGAM() elif not issubclass(lingam_model, _BaseLiNGAM): raise ValueError('lingam_model must be a subclass of _BaseLiNGAM') @@ -83,7 +82,7 @@ def fit(self, X): lags = M_taus.shape[0] residuals = self._calc_residuals(X, M_taus, lags) - model = DirectLiNGAM() + model = lingam_model model.fit(residuals) B_taus = self._calc_b(X, model.adjacency_matrix_, M_taus) @@ -160,8 +159,22 @@ def bootstrap(self, X, n_sampling): def _estimate_var_coefs(self, X): """Estimate coefficients of VAR""" - var = VAR(X) - result = var.fit(maxlags=self._lags, ic=self._criterion, trend='nc') + # XXX: VAR.fit() is not searching lags correctly + if self._criterion not in ['aic', 'fpe', 'hqic', 'bic']: + var = VAR(X) + result = var.fit(maxlags=self._lags, trend='nc') + else: + min_value = float('Inf') + result = None + + for lag in range(1, self._lags + 1): + var = VAR(X) + fitted = var.fit(maxlags=lag, ic=None, trend='nc') + + value = getattr(fitted, self._criterion) + if value < min_value: + min_value = value + result = fitted return result.coefs, result.k_ar, result.resid diff --git a/lingam/varma_lingam.py b/lingam/varma_lingam.py index 3de829a..1336a11 100644 --- a/lingam/varma_lingam.py +++ b/lingam/varma_lingam.py @@ -7,10 +7,9 @@ from sklearn.utils import check_array, resample from statsmodels.tsa.statespace.varmax import VARMAX -from lingam import DirectLiNGAM - from .base import _BaseLiNGAM from .bootstrap import BootstrapResult +from .direct_lingam import DirectLiNGAM class VARMALiNGAM: @@ -42,8 +41,8 @@ def __init__(self, order=(1, 1), criterion='bic', prune=False, max_iter=100, ar_ ma_coefs : array-like, optional (default=None) Coefficients of MA of ARMA. Estimating ARMA model is skipped if specified ``ar_coefs`` and `ma_coefs`. Shape must be (``order[1]``, n_features, n_features). - lingam_model : constructor - Constructor of a LiNGAM algorithm which inherits _BaseLiNGAM. + lingam_model : lingam object inherits 'lingam._BaseLiNGAM', optional (default=None) + LiNGAM model for causal discovery. If None, DirectLiNGAM algorithm is selected. random_state : int, optional (default=None) ``random_state`` is the seed used by the random number generator. """ @@ -79,7 +78,7 @@ def fit(self, X): lingam_model = self._lingam_model if lingam_model is None: - lingam_model = DirectLiNGAM + lingam_model = DirectLiNGAM() elif not issubclass(lingam_model, _BaseLiNGAM): raise ValueError('lingam_model must be a subclass of _BaseLiNGAM') @@ -95,7 +94,7 @@ def fit(self, X): residuals = self._calc_residuals( X, phis, thetas, p, q) - model = DirectLiNGAM() + model = lingam_model model.fit(residuals) psis, omegas = self._calc_psi_and_omega( diff --git a/setup.py b/setup.py index 8dbebdb..a225c70 100644 --- a/setup.py +++ b/setup.py @@ -9,7 +9,7 @@ setuptools.setup( name='lingam', version=VERSION, - author='T.Ikeuchi, G.Haraoka, S.Shimizu', + author='T.Ikeuchi, G.Haraoka, W.Kurebayashi, S.Shimizu', description='LiNGAM Python Package', long_description=README, long_description_content_type='text/markdown', diff --git a/tests/test_longitudinal_lingam.py b/tests/test_longitudinal_lingam.py new file mode 100644 index 0000000..75686ed --- /dev/null +++ b/tests/test_longitudinal_lingam.py @@ -0,0 +1,421 @@ +import os + +import numpy as np +import pandas as pd + +from lingam.longitudinal_lingam import LongitudinalLiNGAM + + +def test_fit_success(): + # causal direction: x0 --> x1 --> x3 + x0 = np.random.uniform(size=1000) + x1 = 0.7*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + x3 = 0.3*x1 + np.random.uniform(size=1000) + X1 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + x0 = np.random.uniform(size=1000) + 0.5*x2 + x1 = 0.3*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + 0.5*x3 + x3 = 0.7*x1 + np.random.uniform(size=1000) + X2 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + x0 = np.random.uniform(size=1000) + 0.5*x2 + x1 = 0.5*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + 0.5*x3 + x3 = 0.5*x1 + np.random.uniform(size=1000) + X3 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + X_list = np.empty((3, 1000, 4)) + X_list[0] = X1 + X_list[1] = X2 + X_list[2] = X3 + + model = LongitudinalLiNGAM() + model.fit(X_list) + + # check causal ordering + cos = model.causal_orders_ + for co in cos: + assert co.index(0) < co.index(1) < co.index(3) + + # check B(t,t) + B_t = model.adjacency_matrices_[1, 0] # B(1,1) + assert B_t[1, 0] > 0.2 and B_t[3, 1] > 0.6 + B_t[1, 0] = 0.0 + B_t[3, 1] = 0.0 + assert np.sum(B_t) < 0.1 + + B_t = model.adjacency_matrices_[2, 0] # B(2,2) + assert B_t[1, 0] > 0.4 and B_t[3, 1] > 0.4 + B_t[1, 0] = 0.0 + B_t[3, 1] = 0.0 + assert np.sum(B_t) < 0.1 + + # check B(t,t-τ) + B_tau = model.adjacency_matrices_[1, 1] # B(1,0) + assert B_tau[0, 2] > 0.4 and B_tau[2, 3] > 0.4 + + B_tau = model.adjacency_matrices_[1, 1] # B(2,1) + assert B_tau[0, 2] > 0.4 and B_tau[2, 3] > 0.4 + + # fit by list + X_list = [X1, X2, X3] + model = LongitudinalLiNGAM() + model.fit(X_list) + + +def test_fit_invalid_data(): + # Different features + x0 = np.random.uniform(size=1000) + x1 = 0.7*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + x3 = 0.3*x1 + np.random.uniform(size=1000) + X1 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + x0 = np.random.uniform(size=1000) + 0.5*x2 + x1 = 0.3*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + 0.5*x3 + X2 = pd.DataFrame(np.array([x0, x1, x2]).T, + columns=['x0', 'x1', 'x2']) + + x0 = np.random.uniform(size=1000) + 0.5*x2 + x1 = 0.5*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + 0.5*x3 + x3 = 0.5*x1 + np.random.uniform(size=1000) + X3 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + X_list = [X1, X2, X3] + + try: + model = LongitudinalLiNGAM() + model.fit(X_list) + except ValueError: + pass + else: + raise AssertionError + + # Not list data + X = 1 + try: + model = LongitudinalLiNGAM() + model.fit(X) + except ValueError: + pass + else: + raise AssertionError + + # Include not-array data + x0 = np.random.uniform(size=1000) + x1 = 2.0*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + x3 = 4.0*x1 + np.random.uniform(size=1000) + X1 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + X_list = [X1, 1] + + try: + model = LongitudinalLiNGAM() + model.fit(X) + except ValueError: + pass + else: + raise AssertionError + + # Include non-numeric data + x0 = np.random.uniform(size=1000) + x1 = 0.7*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + x3 = 0.3*x1 + np.random.uniform(size=1000) + X1 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + x0 = np.random.uniform(size=1000) + 0.5*x2 + x1 = np.array(['X']*1000) # <== non-numeric + x2 = np.random.uniform(size=1000) + 0.5*x3 + x3 = np.array(['X']*1000) # <== non-numeric + X2 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + x0 = np.random.uniform(size=1000) + 0.5*x2 + x1 = 0.5*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + x3 = 0.5*x1 + np.random.uniform(size=1000) + X3 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + X_list = [X1, X2, X3] + + try: + model = LongitudinalLiNGAM() + model.fit(X_list) + except ValueError: + pass + else: + raise AssertionError + + # Include NaN values + x0 = np.random.uniform(size=1000) + x1 = 0.7*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + x3 = 0.3*x1 + np.random.uniform(size=1000) + X1 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + x0 = np.random.uniform(size=1000) + 0.5*x2 + x1 = 0.3*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + 0.5*x3 + x3 = 0.7*x1 + np.random.uniform(size=1000) + X2 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + X2.iloc[100, 0] = np.nan # set nan + + x0 = np.random.uniform(size=1000) + 0.5*x2 + x1 = 0.5*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + 0.5*x3 + x3 = 0.5*x1 + np.random.uniform(size=1000) + X3 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + X_list = [X1, X2, X3] + + try: + model = LongitudinalLiNGAM() + model.fit(X_list) + except ValueError: + pass + else: + raise AssertionError + + # Include infinite values + x0 = np.random.uniform(size=1000) + x1 = 0.7*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + x3 = 0.3*x1 + np.random.uniform(size=1000) + X1 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + x0 = np.random.uniform(size=1000) + 0.5*x2 + x1 = 0.3*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + 0.5*x3 + x3 = 0.7*x1 + np.random.uniform(size=1000) + X2 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + X2.iloc[100, 0] = np.inf # set inf + + x0 = np.random.uniform(size=1000) + 0.5*x2 + x1 = 0.5*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + 0.5*x3 + x3 = 0.5*x1 + np.random.uniform(size=1000) + X3 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + X_list = [X1, X2, X3] + + try: + model = LongitudinalLiNGAM() + model.fit(X_list) + except ValueError: + pass + else: + raise AssertionError + +def test_bootstrap_success(): + # causal direction: x0 --> x1 --> x3 + x0 = np.random.uniform(size=1000) + x1 = 0.7*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + x3 = 0.3*x1 + np.random.uniform(size=1000) + X1 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + x0 = np.random.uniform(size=1000) + 0.5*x2 + x1 = 0.3*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + 0.5*x3 + x3 = 0.7*x1 + np.random.uniform(size=1000) + X2 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + x0 = np.random.uniform(size=1000) + 0.5*x2 + x1 = 0.5*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + 0.5*x3 + x3 = 0.5*x1 + np.random.uniform(size=1000) + X3 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + X_list = np.empty((3, 1000, 4)) + X_list[0] = X1 + X_list[1] = X2 + X_list[2] = X3 + + model = LongitudinalLiNGAM() + model.bootstrap(X_list, n_sampling=3) + + # fit by list + X_list = [X1, X2, X3] + model = LongitudinalLiNGAM() + model.bootstrap(X_list, n_sampling=3) + +def test_bootstrap_invalid_data(): + # Different features + x0 = np.random.uniform(size=1000) + x1 = 0.7*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + x3 = 0.3*x1 + np.random.uniform(size=1000) + X1 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + x0 = np.random.uniform(size=1000) + 0.5*x2 + x1 = 0.3*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + 0.5*x3 + X2 = pd.DataFrame(np.array([x0, x1, x2]).T, + columns=['x0', 'x1', 'x2']) + + x0 = np.random.uniform(size=1000) + 0.5*x2 + x1 = 0.5*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + 0.5*x3 + x3 = 0.5*x1 + np.random.uniform(size=1000) + X3 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + X_list = [X1, X2, X3] + + try: + model = LongitudinalLiNGAM() + model.bootstrap(X_list, n_sampling=3) + except ValueError: + pass + else: + raise AssertionError + + # Not list data + X = 1 + try: + model = LongitudinalLiNGAM() + model.bootstrap(X, n_sampling=3) + except ValueError: + pass + else: + raise AssertionError + + # Include not-array data + x0 = np.random.uniform(size=1000) + x1 = 2.0*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + x3 = 4.0*x1 + np.random.uniform(size=1000) + X1 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + X_list = [X1, 1] + + try: + model = LongitudinalLiNGAM() + model.bootstrap(X, n_sampling=3) + except ValueError: + pass + else: + raise AssertionError + + # Include non-numeric data + x0 = np.random.uniform(size=1000) + x1 = 0.7*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + x3 = 0.3*x1 + np.random.uniform(size=1000) + X1 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + x0 = np.random.uniform(size=1000) + 0.5*x2 + x1 = np.array(['X']*1000) # <== non-numeric + x2 = np.random.uniform(size=1000) + 0.5*x3 + x3 = np.array(['X']*1000) # <== non-numeric + X2 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + x0 = np.random.uniform(size=1000) + 0.5*x2 + x1 = 0.5*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + x3 = 0.5*x1 + np.random.uniform(size=1000) + X3 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + X_list = [X1, X2, X3] + + try: + model = LongitudinalLiNGAM() + model.bootstrap(X_list, n_sampling=3) + except ValueError: + pass + else: + raise AssertionError + + # Include NaN values + x0 = np.random.uniform(size=1000) + x1 = 0.7*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + x3 = 0.3*x1 + np.random.uniform(size=1000) + X1 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + x0 = np.random.uniform(size=1000) + 0.5*x2 + x1 = 0.3*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + 0.5*x3 + x3 = 0.7*x1 + np.random.uniform(size=1000) + X2 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + X2.iloc[100, 0] = np.nan # set nan + + x0 = np.random.uniform(size=1000) + 0.5*x2 + x1 = 0.5*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + 0.5*x3 + x3 = 0.5*x1 + np.random.uniform(size=1000) + X3 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + X_list = [X1, X2, X3] + + try: + model = LongitudinalLiNGAM() + model.bootstrap(X_list, n_sampling=3) + except ValueError: + pass + else: + raise AssertionError + + # Include infinite values + x0 = np.random.uniform(size=1000) + x1 = 0.7*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + x3 = 0.3*x1 + np.random.uniform(size=1000) + X1 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + x0 = np.random.uniform(size=1000) + 0.5*x2 + x1 = 0.3*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + 0.5*x3 + x3 = 0.7*x1 + np.random.uniform(size=1000) + X2 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + X2.iloc[100, 0] = np.inf # set inf + + x0 = np.random.uniform(size=1000) + 0.5*x2 + x1 = 0.5*x0 + np.random.uniform(size=1000) + x2 = np.random.uniform(size=1000) + 0.5*x3 + x3 = 0.5*x1 + np.random.uniform(size=1000) + X3 = pd.DataFrame(np.array([x0, x1, x2, x3]).T, + columns=['x0', 'x1', 'x2', 'x3']) + + X_list = [X1, X2, X3] + + try: + model = LongitudinalLiNGAM() + model.bootstrap(X_list, n_sampling=3) + except ValueError: + pass + else: + raise AssertionError diff --git a/tests/test_var_lingam.py b/tests/test_var_lingam.py index bcf3148..9345146 100644 --- a/tests/test_var_lingam.py +++ b/tests/test_var_lingam.py @@ -101,7 +101,7 @@ def test_fit_success(): assert np.sum(b0) < 0.1 b1 = model.adjacency_matrices_[1] - assert b1[0, 0] < -0.3 and b1[0, 1] < -0.3 and b1[0, 2] < -0.3 \ + assert b1[0, 0] < -0.3 and b1[0, 1] < -0.3 and b1[0, 2] < -0.25 \ and b1[1, 2] > 0.05 and b1[2, 1] < -0.1 and b1[2, 2] < -0.3 b1[0, 0] = 0.0 diff --git a/tests/test_varma_lingam.py b/tests/test_varma_lingam.py index 8e9edb3..78b02f8 100644 --- a/tests/test_varma_lingam.py +++ b/tests/test_varma_lingam.py @@ -129,10 +129,10 @@ def generate_data(n=5, T=800, initial_data=None): expon = 0.1 ext = np.empty((n, T)) for i in range(n): - ext[i, :] = np.random.normal(size=(1, T)); - ext[i, :] = np.multiply(np.sign(ext[i, :]), abs(ext[i, :]) ** expon); - ext[i, :] = ext[i, :] - np.mean(ext[i, :]); - ext[i, :] = ext[i, :] / np.std(ext[i, :]); + ext[i, :] = np.random.normal(size=(1, T)) + ext[i, :] = np.multiply(np.sign(ext[i, :]), abs(ext[i, :]) ** expon) + ext[i, :] = ext[i, :] - np.mean(ext[i, :]) + ext[i, :] = ext[i, :] / np.std(ext[i, :]) # observed signals y y = np.zeros((n, T))