causality_for_spaceweather/notebooks/_build/jupyter_execute/chapter1_it/IDTxl.py

#!/usr/bin/env python
# coding: utf-8

# # IDTxl基本使用教程
# 
# 安装和依赖请查看：
# 
# [Installation and Requirements](https://github.com/pwollstadt/IDTxl/wiki/Installation-and-Requirements)

# In[1]:


get_ipython().run_line_magic('reload_ext', 'autoreload')
get_ipython().run_line_magic('autoreload', '2')
get_ipython().run_line_magic('matplotlib', 'inline')


# ## 使用IDTxl做第一个网络推断

# In[2]:


# %load test_first_example.py
# Import classes
from idtxl.multivariate_te import MultivariateTE
from idtxl.data import Data
from idtxl.visualise_graph import plot_network
import matplotlib.pyplot as plt
# a) Generate test data
data = Data()
data.generate_mute_data(n_samples=1000, n_replications=1)

# b) Initialise analysis object and define settings
network_analysis = MultivariateTE()
settings = {'cmi_estimator': 'JidtGaussianCMI',
            'max_lag_sources': 5,
            'min_lag_sources': 1}

# c) Run analysis
results = network_analysis.analyse_network(settings=settings, data=data)

# d) Plot inferred network to console and via matplotlib
results.print_edge_list(weights='max_te_lag', fdr=False)
plot_network(results=results, weights='max_te_lag', fdr=False)
plt.show()


# ## 理论介绍
# 
# 为从数据中得到有用的结果，我们需要对算法的理论有一个初步的了解

# 为什么使用多变量转移熵呢？

# 转移熵最早由Schreiber在2000年提出，它是一种度量两个随机过程之间信息流动的有向估计方法。然而系统中往往是包含多个变量的，此时二元转移熵可能会产生一些错误的结果：
# 
# 1. Spurious or redundant interactions
# 2. Synergistic interactions
# 
# 于是，我们需要一种可以考虑多个变量之间相互作用的方法。
# 
# IDTxl实现了Lizier在2012年提出的一个贪婪迭代算法，有如下好处：
# 
# 1. The algorithm infers all relevant sources of a target by iteratively including variables from a source's past that maximise a conditionsl mutual information criterion.
# 2. IDTxl builds the set of parent sources for each target node in the network. This iterative conditioning is designed to both remove redundancies and capture synergistic interactions in building each parent set, thus addressing the two aforementioned shortcomings of bivariate analysis.
# 3. THe inclusion of source variables requires repeated statistical testing of each investigated past variable in each iteration. IDTxl handles this multiple testing, in particular the family-wise error rate(FWER), by implementing a so-called maximum statistic test.
# 4. Further statistical tests are aimed at pruning the selected parent set and providing control of the FMER slao at the network level.
# 5. At the same time, IDTcl automates the optimization of parameters necessary for the estimation of mTE.
# 
# 前四步就是2012文章里的算法所说，最后一个在哪里，这个很重要！！

# 关于更详细的介绍，可以去读2012年那篇文章，或者查看 [Theoretical Introduction](https://github.com/pwollstadt/IDTxl/wiki/Theoretical-Introduction)

# ## 两个主要应用

# 1. Network inference
# 2. Node dynamics

# ## 网络推断的几种方法

# ###  Multivariate Transfer Entropy / Granger causality

# In[3]:


# %load ./IDTxl/demos/demo_multivariate_te.py
# Import classes
from idtxl.multivariate_te import MultivariateTE
from idtxl.data import Data
from idtxl.visualise_graph import plot_network
import matplotlib.pyplot as plt

# a) Generate test data
data = Data()
data.generate_mute_data(n_samples=1000, n_replications=5)

# b) Initialise analysis object and define settings
network_analysis = MultivariateTE()
settings = {'cmi_estimator': 'JidtGaussianCMI',
            'max_lag_sources': 5,
            'min_lag_sources': 1}

# c) Run analysis
results = network_analysis.analyse_network(settings=settings, data=data)

# d) Plot inferred network to console and via matplotlib
results.print_edge_list(weights='max_te_lag', fdr=False)
plot_network(results=results, weights='max_te_lag', fdr=False)
plt.show()


# In[4]:


# a) Generate test data
data = Data()
data.generate_mute_data(n_samples=1000, n_replications=5)


# In[5]:


data.data.shape


# In[6]:


# b) Initialise analysis object and define settings
network_analysis = MultivariateTE()
settings = {'cmi_estimator': 'JidtGaussianCMI',
            'max_lag_sources': 16,
            'min_lag_sources': 1}

# c) Run analysis
results = network_analysis.analyse_single_target(settings=settings,
                                                 data=data,
                                                 target=0,
                                                 sources=[1, 3])


# In[7]:


results.print_edge_list(weights='max_te_lag', fdr=False)
plot_network(results=results, weights='max_te_lag', fdr=False)
plt.show()


# ###  Bivariate TE / Granger causality

# In[8]:


# %load ./IDTxl/demos/demo_bivariate_te.py
# Import classes
from idtxl.bivariate_te import BivariateTE
from idtxl.data import Data
from idtxl.visualise_graph import plot_network
import matplotlib.pyplot as plt

# a) Generate test data
data = Data()
data.generate_mute_data(n_samples=1000, n_replications=5)

# b) Initialise analysis object and define settings
network_analysis = BivariateTE()
settings = {'cmi_estimator': 'JidtGaussianCMI',
            'max_lag_sources': 5,
            'min_lag_sources': 1}

# c) Run analysis
results = network_analysis.analyse_network(settings=settings, data=data)

# d) Plot inferred network to console and via matplotlib
results.print_edge_list(weights='max_te_lag', fdr=False)
plot_network(results=results, weights='max_te_lag', fdr=False)
plt.show()


# ### Multivariate mutual information

# In[9]:


# %load ./IDTxl/demos/demo_multivariate_mi.py
# Import classes
from idtxl.multivariate_mi import MultivariateMI
from idtxl.data import Data
from idtxl.visualise_graph import plot_network
import matplotlib.pyplot as plt

# a) Generate test data
data = Data()
data.generate_mute_data(n_samples=1000, n_replications=5)

# b) Initialise analysis object and define settings
network_analysis = MultivariateMI()
settings = {'cmi_estimator': 'JidtGaussianCMI',
            'max_lag_sources': 5,
            'min_lag_sources': 1}

# c) Run analysis
results = network_analysis.analyse_network(settings=settings, data=data)

# d) Plot inferred network to console and via matplotlib
results.print_edge_list(weights='max_te_lag', fdr=False)
plot_network(results=results, weights='max_te_lag', fdr=False)
plt.show()


# ### Bivariate mutual information

# In[10]:


# %load ./IDTxl/demos/demo_bivariate_mi.py
# Import classes
from idtxl.bivariate_mi import BivariateMI
from idtxl.data import Data
from idtxl.visualise_graph import plot_network
import matplotlib.pyplot as plt

# a) Generate test data
data = Data()
data.generate_mute_data(n_samples=1000, n_replications=5)

# b) Initialise analysis object and define settings
network_analysis = BivariateMI()
settings = {'cmi_estimator': 'JidtGaussianCMI',
            'max_lag_sources': 5,
            'min_lag_sources': 1}

# c) Run analysis
results = network_analysis.analyse_network(settings=settings, data=data)

# d) Plot inferred network to console and via matplotlib
results.print_edge_list(weights='max_te_lag', fdr=False)
plot_network(results=results, weights='max_te_lag', fdr=False)
plt.show()


# ## 网络推断实现细节

# ### 数据处理

# IDTxl uses its own class to handle data

# In[11]:


from idtxl.data import Data
import numpy as np


# The Data class holds up to 3D data, where 
# 
# 1. one dimension represents processes, 
# 2. one dimension represents samples, 
# 3. one dimension represents replications. 
# 
# For example, in a neuroscience setting, processes would correspond to EEG channels, samples to time steps, and replications to trials (repetitions of the same experiment).

# Initialise Data Object

# To import your own data set, create a Data object passing a 1D, 2D, or 3D numpy array. To specify which axis of the numpy array represents which dimension of the data, pass a one to three-letter string as `dim_order`, where 'p' stands for processes, 's' for samples, and 'r' for replications.

# ** Example:  3D numpy array **

# In[12]:


# Initialise a data object holding data with 50000 samples,
# 10 processes, and 5 replications.
mydata = np.arange(50000).reshape((1000, 10, 5))
data = Data(mydata, dim_order='spr')


# ** Example: 2D numpy array **

# In[13]:


# Initialise a data object holding data with 5 processes and
# 5000 samples
mydata = np.arange(5000).reshape((5, 1000))
dat = Data(mydata, dim_order='ps')


# In[14]:


mydata.shape


# In[15]:


dat.data.shape


# ** Example: 1D numpy array **

# In[16]:


# Initialise a data object holding data with 1 process and
# 1000 samples
mydata = np.arange(1000)
dat = Data(mydata, dim_order='s')


# **Example: pandas data frame**

# In[17]:


import pandas as pd
df = pd.DataFrame(
    {'col1': np.random.rand(100), 'col2': np.random.rand(100)})
data = Data(df, dim_order='sp')
# Adding data with properties: 2 processes, 100 samples, 1 replications
# overwriting existing data


# In[18]:


data.data.shape


# IDTxl also provides the functionality to [generate synthetic test data](https://github.com/pwollstadt/IDTxl/wiki/Synthetic-Test-Data)

# ### 设置模型参数

# #### 选择合适的CMI 参数

# CMI估计方法的选取很重要，对不同的数据类型应该选择不同的参数，主要是在setting里面设置`cmi_estimator`这个参数。
# 
# 主要数据类型：
# 
# 1. discrete data
# 2. jointly Gaussian continuous data
# 3. non-linear continuous data
# 
# 对于非线性的连续数据有两种估计选择：
# 
# 1. `JidtKraskovCMI`
# 2. `OpenCLKraskovCMI`

# In[19]:


import idtxl.estimators_jidt
help(idtxl.estimators_jidt)


# Different CMI estimators are available for discrete data, jointly Gaussian continuous data, and non-linear continuous data. For non-linear continuous data, two alternatives are available:
# 1. A multithreaded Java-implementation (JidtKraskovCMI) 
# 2. An OpenCL-implementation (OpenCLKraskovCMI, to be run on GPUs). 

# In[20]:


settings = {'cmi_estimator': 'JidtKraskovCMI'}


# In[21]:


import idtxl.estimators_opencl
help(idtxl.estimators_opencl)


# #### 设置Source and Target lags
# 

# 由于计算资源有限，我们不可能考虑所有的时间延迟，所以可以设置一个延迟计算范围，甚至可以设置延迟计算间隔，仅仅计算有效的延迟，合理利用先验知识，具体就是在`setting`里面设置如下参数：
# 
# `max_lag_sources`, `min_lag_sources`, `max_lag_target`, `min_lag_target`, `tau_sources`, `tau_target`
# 
# 看情况待定吧，一般要记得设置 `max_lag_sources`，否则就会计算很多无用的东西。

# In[22]:


settings = {'max_lag_sources': 16,
            'min_lag_sources': 6}


# For transfer entropy estimation, additionally, a maximum lag for the target can be set. If the lag is not specified, it will be set to be equal to the maximum source lag.

# The minimum target lag is always set to 1 sample, to ensure so-called self-prediction optimality (Wibral, 2013, PLOS 8(2)). This is further discussed in the Wiki's Theoretical Introduction.

# In[23]:


settings = {'max_lag_sources': 16,
            'min_lag_sources': 6,
            'max_lag_target': 9}


# ### 重要性统计检验

# The `stats` module provides the implementation of several statistical significance tests.

# 在网格推断中主要有五个重要性统计检验：
# 1. `Maximum test`
# 2. `Minimum test`
# 3. `Omnibus test`
# 4. `Sequential maximum test`
# 5. `FDR correction`
# 
# 关于它们的详细介绍和完整参数设置，请查看源码或说明文档：[Documentation](https://pwollstadt.github.io/IDTxl/html/index.html)， [stats section](https://github.com/pwollstadt/IDTxl/wiki/Theoretical-Introduction#statistical-tests-used-in-the-mte-algorithm)

# 参数设置还是在`setting`里面，主要就是设置打乱次数和一个水平指标，这个指标代表什么含义呢？需要去查一下

# In[24]:


settings = {'n_perm_max_stat': 200,
            'alpha_max_stat': 0.05,
            'permute_in_time': False}


# ### Conditioning on specific variables

# 这个是为了对某些特定的变量做一个条件测试，设置方法如下所示，一般我们不会特别指明

# In order to enforce the inclusion of specific variables in the condition set, add the `add_conditionals` key in the analysis settings dictionary and provide the list of variables to include. Each variable is defined as as tuple `(process index, lag)`, e.g. `(0,3)` corresponds to the variable in process 0 with a lag of 3 samples with respect to the target.

# In[25]:


settings = {'add_conditionals': [(0,3), (1,2)]}


# By setting `add_conditionals` to `faes`, the algorithm adds the current value of all sources to the conditional set. This is meant to correct for instantaneous mixing in the source.

# In[26]:


settings = {'add_conditionals': 'faes'}


# ### 分析过程

# 参数设置好了，我们就可以开始进行分析了

# #### 单个目标
# 
# 主要就是使用`analyse_single_target`这个函数，可以设置`target`和`sources`。

# In order to restrict the analysis to a specific target process in the system, run the `analyse_single_target` method and provide `target` argument in addition to the required `data` and `settings`. 
# 
# By default, when a target process is specified(e.g. process 0), all the other processes will be considered as potential sources.

# In[27]:


data = Data()
data.generate_mute_data(n_samples=1000, n_replications=5)


# In[28]:


settings = {'cmi_estimator': 'JidtKraskovCMI',
            'max_lag_sources': 5,
            'min_lag_sources':1,
            'permute_in_time': False,
            'add_conditionals': [(0,3), (1,2)]}


# `add_conditionals` 这个参数不能随便进行设置，刚才设置成`fae`的时候网络分析会报错

# In[29]:


results = network_analysis.analyse_single_target(settings=settings,
                                                 data=data,
                                                 target=0,
                                                 sources=[1, 3])


# 这个例子计算起来好慢啊

# 这是什么错误呢？

# 问题好像是出现在`add_conditionals`上面

# #### 全网络分析

# To analyse the whole network, use the `analyse_network` method and provide the required `data` and `setting` arguments.

# In[ ]:


results = network_analysis.analyse_network(settings=settings,
                                           data=data)


# THe algorithm will loop over all processes in the network and first analyse each of them separately as a target. The results will finally be combined in a [Results object](https://github.com/pwollstadt/IDTxl/wiki/The-Results-Class). By default, an FDR correction will be performed when combining the results form different target[1995](http://www.math.tau.ac.il/~ybenja/MyPapers/benjamini_hochberg1995.pdf).  
# 
# In order to skip the FDR correction, refer to the tutorial on [Statistical Significance Tests](https://github.com/pwollstadt/IDTxl/wiki/Statistical-Significance-Tests).

# 除此之外，我们还可以限制target variables 和source variables

# #### The result class

# IDTxl所有的计算结果保存在result当中，比如`ActiveInformationStorage()`，`MUltivaraiteTE()`,还有单目标及全网络分析的结果以及参数设置属性列表也都在这里，直接查看属性即可

# 一些初步的应用：
# 1. 使用`result.settings`查看参数设置
# 2. 使用`results.print_edge_list(weights='max_te_lag', fdr=False)`打印出重要的信息流动和延迟
# 3. 使用`get_adjacency_matrix()`可以查看关联关系矩阵
# 4. 可以查看对某个单独目标的贡献
# 5. 可以查看数据的性质
# 6. 可以用来可视化

# 还有一些高级功能，如`combine_results`和`PartialInformationDecomposition`等等，可以查看Documentation以及相关的论文，目前还用不到

# ### 可视化网络分析

# **下面我们来考虑下一个话题：Node dynamics !**

# ## 节点动力学分析

# ### Active Information Storage

# 这个例子其实就是对于单个的随机过程来进行分析的，就是看看$X_n^{(k)}$对$X_{n+1}$的贡献，需要提前设置所要考察的最大的k值。

# In[ ]:


# %load ./IDTxl/demos/demo_active_information_storage.py
# Import classes
from idtxl.active_information_storage import ActiveInformationStorage
from idtxl.data import Data

# a) Generate test data
data = Data()
data.generate_mute_data(n_samples=1000, n_replications=5)

# b) Initialise analysis object and define settings
network_analysis = ActiveInformationStorage()
settings = {'cmi_estimator':  'JidtGaussianCMI',
            'max_lag': 5}

# c) Run analysis
results = network_analysis.analyse_network(settings=settings, data=data)

# d) Plot list of processes with significant AIS to console
print(results.get_significant_processes(fdr=False))


# ### Partial Information Decomposition

# 下面这个是一个异或问题的例子，输入是两个变量，输出是一个变量

# In[ ]:


# %load ./IDTxl/demos/demo_partial_information_decomposition.py
# Import classes
import numpy as np
from idtxl.partial_information_decomposition import (
                                        PartialInformationDecomposition)
from idtxl.data import Data

# a) Generate test data
n = 100
alph = 2
x = np.random.randint(0, alph, n)
y = np.random.randint(0, alph, n)
z = np.logical_xor(x, y).astype(int)
data = Data(np.vstack((x, y, z)), 'ps', normalise=False)

# b) Initialise analysis object and define settings for both PID estimators
pid = PartialInformationDecomposition()
settings_tartu = {'pid_estimator': 'TartuPID', 'lags_pid': [0, 0]}
settings_sydney = {
    'alph_s1': alph,
    'alph_s2': alph,
    'alph_t': alph,
    'max_unsuc_swaps_row_parm': 60,
    'num_reps': 63,
    'max_iters': 1000,
    'pid_estimator': 'SydneyPID',
    'lags_pid': [0, 0]}

# c) Run Tartu estimator
results_tartu = pid.analyse_single_target(
    settings=settings_tartu, data=data, target=2, sources=[0, 1])

# d) Run Sydney estimator
pid = PartialInformationDecomposition()
results_sydney = pid.analyse_single_target(
    settings=settings_sydney, data=data, target=2, sources=[0, 1])

# e) Print results to console
print('\nLogical XOR')
print('Estimator            Sydney\t\tTartu\t\tExpected\n')
print('Uni s1               {0:.4f}\t\t{1:.4f}\t\t{2:.2f}'.format(
    results_sydney.get_single_target(2)['unq_s1'],
    results_tartu.get_single_target(2)['unq_s1'],
    0))
print('Uni s2               {0:.4f}\t\t{1:.4f}\t\t{2:.2f}'.format(
    results_sydney.get_single_target(2)['unq_s2'],
    results_tartu.get_single_target(2)['unq_s2'],
    0))
print('Shared s1_s2         {0:.4f}\t\t{1:.4f}\t\t{2:.2f}'.format(
    results_sydney.get_single_target(2)['shd_s1_s2'],
    results_tartu.get_single_target(2)['shd_s1_s2'],
    0))
print('Synergy s1_s2        {0:.4f}\t\t{1:.4f}\t\t{2:.2f}'.format(
    results_sydney.get_single_target(2)['syn_s1_s2'],
    results_tartu.get_single_target(2)['syn_s1_s2'],
    1))


# The result is actually interesting.
# 
# 可以看到结果还是很不错的，0.9705已经很接近于期望值1喽。

# 这个要去读一下相关的文献了：
# 1. Non negative information decomposition
# 2. Local information measure

# ### 实现细节

# 这部分内容的理论基础是：Local information Measure：具体可以查看2008，2010，2014年的那几篇文献，我们先来看看程序是什么样子~

# 具体同网路推断差不多：数据处理，设置模型参数，选择合适的估计方法，进行重要性统计检验，结果分析等等，具体参考[Documentation](https://github.com/pwollstadt/IDTxl/wiki)

# ## 注意事项

# 1. 参数太多，设置起来要小心

# In[ ]: