Merge pull request #115 from pwollstadt/develop

Add bug fixes related to FDR corrections to master; add dependency on statsmodels

demos/demo_active_information_storage.py

@@ -8,8 +8,7 @@
 
 # b) Initialise analysis object and define settings
 network_analysis = ActiveInformationStorage()
-settings = {'cmi_estimator':  'JidtGaussianCMI',
-            'max_lag': 5}
+settings = {"cmi_estimator": "JidtGaussianCMI", "max_lag": 5}
 
 # c) Run analysis
 results = network_analysis.analyse_network(settings=settings, data=data)

demos/demo_bivariate_mi.py

@@ -1,23 +1,26 @@
 # Import classes
+import matplotlib.pyplot as plt
+
 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}
+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)
+results.print_edge_list(weights="max_te_lag", fdr=False)
+plot_network(results=results, weights="max_te_lag", fdr=False)
 plt.show()

demos/demo_bivariate_pid.py

@@ -1,5 +1,6 @@
 # Import classes
 import numpy as np
+
 from idtxl.bivariate_pid import BivariatePID
 from idtxl.data import Data
 
@@ -9,46 +10,61 @@
 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)
+data = Data(np.vstack((x, y, z)), "ps", normalise=False)
 
 # b) Initialise analysis object and define settings for both PID estimators
 pid = BivariatePID()
-settings_tartu = {'pid_estimator': 'TartuPID', 'lags_pid': [0, 0]}
+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]}
+    "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])
+    settings=settings_tartu, data=data, target=2, sources=[0, 1]
+)
 
 # d) Run Sydney estimator
 pid = BivariatePID()
 results_sydney = pid.analyse_single_target(
-    settings=settings_sydney, data=data, target=2, sources=[0, 1])
+    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))
+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,
+    )
+)

demos/demo_bivariate_te.py

@@ -1,23 +1,26 @@
 # Import classes
+import matplotlib.pyplot as plt
+
 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}
+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)
+results.print_edge_list(weights="max_te_lag", fdr=False)
+plot_network(results=results, weights="max_te_lag", fdr=False)
 plt.show()

demos/demo_core_estimators.py

@@ -1,22 +1,24 @@
 """Demonstrate the usage of IDTxl's core estimators."""
+
 import numpy as np
-from idtxl.estimators_python import PythonKraskovCMI
+
 from idtxl.estimators_jidt import (
     JidtDiscreteAIS,
     JidtDiscreteCMI,
     JidtDiscreteMI,
     JidtDiscreteTE,
-    JidtKraskovAIS,
-    JidtKraskovCMI,
-    JidtKraskovMI,
-    JidtKraskovTE,
     JidtGaussianAIS,
     JidtGaussianCMI,
     JidtGaussianMI,
     JidtGaussianTE,
+    JidtKraskovAIS,
+    JidtKraskovCMI,
+    JidtKraskovMI,
+    JidtKraskovTE,
 )
-from idtxl.estimators_opencl import OpenCLKraskovMI, OpenCLKraskovCMI
+from idtxl.estimators_opencl import OpenCLKraskovCMI, OpenCLKraskovMI
 from idtxl.estimators_pid import SydneyPID, TartuPID
+from idtxl.estimators_python import PythonKraskovCMI
 from idtxl.idtxl_utils import calculate_mi
 
 # Generate Gaussian test data

demos/demo_hd_estimator.py

@@ -12,15 +12,17 @@
 
 >>> test_fast_emb.py
 """
+
 # Import classes
 import numpy as np
+
 from idtxl.data_spiketime import Data_spiketime
 from idtxl.embedding_optimization_ais_Rudelt import OptimizationRudelt
 
 # a) Generate test data
 print("\nTest optimization_Rudelt using BBC estimator on Rudelt data")
-data = Data_spiketime()             # initialise empty data object
-data.load_Rudelt_data()             # load Rudelt spike time data
+data = Data_spiketime()  # initialise empty data object
+data.load_Rudelt_data()  # load Rudelt spike time data
 
 # b) Initialise analysis object and define settings
 # Run optimization with the Bayesian bias criterion (BBC) (alternatively, set
@@ -29,20 +31,22 @@
 # 'embedding_number_of_bins_set', and 'embedding_scaling_exponent_set' to
 # reduce the run time of this demo. Use defaults if unsure about the settings.
 settings = {
-    'embedding_past_range_set': [0.005, 0.31548, 5.0],
-    'embedding_number_of_bins_set': [1, 3, 5],
-    'embedding_scaling_exponent_set':
-        {'number_of_scalings': 3,
-         'min_first_bin_size': 0.005,
-         'min_step_for_scaling': 0.01},
-    'number_of_bootstraps_R_max': 10,
-    'number_of_bootstraps_R_tot': 10,
-    'auto_MI_bin_size_set': [0.01, 0.025, 0.05, 0.25],
-    'estimation_method': 'bbc',
-    'debug': True,
-    'visualization': True,
-    'output_path': '.',
-    'output_prefix': 'systemtest_optimizationRudelt_image1'}
+    "embedding_past_range_set": [0.005, 0.31548, 5.0],
+    "embedding_number_of_bins_set": [1, 3, 5],
+    "embedding_scaling_exponent_set": {
+        "number_of_scalings": 3,
+        "min_first_bin_size": 0.005,
+        "min_step_for_scaling": 0.01,
+    },
+    "number_of_bootstraps_R_max": 10,
+    "number_of_bootstraps_R_tot": 10,
+    "auto_MI_bin_size_set": [0.01, 0.025, 0.05, 0.25],
+    "estimation_method": "bbc",
+    "debug": True,
+    "visualization": True,
+    "output_path": ".",
+    "output_prefix": "systemtest_optimizationRudelt_image1",
+}
 hd_estimator = OptimizationRudelt(settings)
 
 # c) Run analysis and save a plot of the results to current directory.

demos/demo_multivariate_mi.py

@@ -1,23 +1,26 @@
 # Import classes
-from idtxl.multivariate_mi import MultivariateMI
+import matplotlib.pyplot as plt
+
 from idtxl.data import Data
+from idtxl.multivariate_mi import MultivariateMI
 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}
+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)
+results.print_edge_list(weights="max_te_lag", fdr=False)
+plot_network(results=results, weights="max_te_lag", fdr=False)
 plt.show()