From 96a179f82f181c3c772860ecf3397ed5a6d91af7 Mon Sep 17 00:00:00 2001 From: Patricia Wollstadt Date: Sat, 24 Nov 2018 17:20:33 +0100 Subject: [PATCH 1/3] Add class to handle adjacency matrices Add class to handle adjacency matrices. Zero-lag analysis may result in networks that contain edges with a delay of zero. The previous handling of adjacency matrices as numpy arrays where non-zero entries represent weighted edges didn't allow for edges with a delay of zero. The AdjacencyMatrix class contains a boolean matrix for edges and a weight matrix of varying type, allowing for weights of value zero. Adapt unit tests and functions expecting adjacency matrices as input accordingly. Update test data. Fixes #20. Fix minor bugs in network plotting. --- idtxl/idtxl_io.py | 20 +- idtxl/network_comparison.py | 4 +- idtxl/results.py | 195 ++++++++++++++------ idtxl/visualise_graph.py | 35 ++-- test/data/mute_results_0.p | Bin 2052 -> 2199 bytes test/data/mute_results_1.p | Bin 2389 -> 2525 bytes test/data/mute_results_2.p | Bin 2052 -> 2052 bytes test/data/mute_results_3.p | Bin 2547 -> 2787 bytes test/data/mute_results_4.p | Bin 2389 -> 2611 bytes test/data/mute_results_full.p | Bin 3996 -> 4092 bytes test/generate_test_data.py | 8 +- test/systemtest_multivariate_te.py | 3 +- test/systemtest_multivariate_te_discrete.py | 3 +- test/systemtest_network_comparison.py | 33 ++-- test/test_idtxl_io.py | 4 - test/test_network_comparison.py | 21 ++- test/test_results.py | 111 ++++++++--- 17 files changed, 292 insertions(+), 145 deletions(-) diff --git a/idtxl/idtxl_io.py b/idtxl/idtxl_io.py index 01ef1748..338ddc0a 100644 --- a/idtxl/idtxl_io.py +++ b/idtxl/idtxl_io.py @@ -375,7 +375,7 @@ class for directed graphs (DiGraph). Multiple options for the weight are available (see documentation of method get_adjacency_matrix for details). Args: - adjacency_matrix : 2D numpy array + adjacency_matrix : AdjacencyMatrix instances adjacency matrix to be exported, returned by get_adjacency_matrix() method of Results() class weights : str @@ -389,9 +389,9 @@ class for directed graphs (DiGraph). Multiple options for the weight are """ # use 'weights' parameter (string) as networkx edge property name and use # adjacency matrix entries as edge property values - custom_type = [(weights, type(adjacency_matrix[0, 0]))] - custom_npmatrix = np.matrix(adjacency_matrix, dtype=custom_type) - return nx.from_numpy_matrix(custom_npmatrix, create_using=nx.DiGraph()) + G = nx.DiGraph() + G.add_weighted_edges_from(adjacency_matrix.get_edge_list(), weights) + return G def export_networkx_source_graph(results, target, sign_sources=True, fdr=True): @@ -410,7 +410,7 @@ def export_networkx_source_graph(results, target, sign_sources=True, fdr=True): target : int target index sign_sources : bool [optional] - add only sources with significant information contribution + add sources with significant information contribution only (default=True) fdr : bool [optional] return FDR-corrected results (default=True) @@ -494,7 +494,7 @@ def export_brain_net_viewer(adjacency_matrix, mni_coord, file_name, **kwargs): https://doi.org/10.1371/journal.pone.0068910 Args: - adjacency_matrix : 2D numpy array + adjacency_matrix : AdjacencyMatrix instance adjacency matrix to be exported, returned by get_adjacency_matrix() method of Results() class mni_coord : numpy array @@ -516,15 +516,13 @@ def export_brain_net_viewer(adjacency_matrix, mni_coord, file_name, **kwargs): """ # Check input and get default settings for plotting. The default for # node labels is a list of '-' (no labels). - n_nodes = adjacency_matrix.shape[0] - n_edges = np.sum(adjacency_matrix > 0) + n_nodes = adjacency_matrix.n_nodes() + n_edges = adjacency_matrix.n_edges() labels = kwargs.get('labels', ['-' for i in range(n_nodes)]) node_color = kwargs.get('node_color', np.ones(n_nodes)) node_size = kwargs.get('node_size', np.ones(n_nodes)) if n_edges == 0: Warning('No edges in results file. Nothing to plot.') - assert adjacency_matrix.shape[0] == adjacency_matrix.shape[1], ( - 'Adjacency matrix must be quadratic.') assert mni_coord.shape[0] == n_nodes and mni_coord.shape[1] == 3, ( 'MNI coordinates must have shape [n_nodes, 3].') assert len(labels) == n_nodes, ( @@ -551,6 +549,6 @@ def export_brain_net_viewer(adjacency_matrix, mni_coord, file_name, **kwargs): with open('{0}.edge'.format(file_name), 'w') as text_file: for i in range(n_nodes): for j in range(n_nodes): - print('{0}\t'.format(adjacency_matrix[i, j]), + print('{0}\t'.format(adjacency_matrix._edge_matrix[i, j]), file=text_file, end='') print('', file=text_file) diff --git a/idtxl/network_comparison.py b/idtxl/network_comparison.py index 04c2b313..f3bed906 100755 --- a/idtxl/network_comparison.py +++ b/idtxl/network_comparison.py @@ -129,8 +129,8 @@ def compare_links_within(self, settings, link_a, link_b, network, data): settings=self.settings, union_network=self.union, results={ - 'cmi_diff_abs': {link_a[1]: [np.abs(self.cmi_diff)], - link_b[1]: [np.abs(self.cmi_diff)]}, + 'cmi_diff_abs': {link_a[1]: np.abs(self.cmi_diff), + link_b[1]: np.abs(self.cmi_diff)}, 'a>b': {link_a[1]: [te_a > te_b], link_b[1]: [te_a > te_b]}, 'pval': {link_a[1]: [pvalue], link_b[1]: [pvalue]}, 'cmi_surr': self.cmi_surr, diff --git a/idtxl/results.py b/idtxl/results.py index ef09ab77..a38481af 100644 --- a/idtxl/results.py +++ b/idtxl/results.py @@ -1,8 +1,13 @@ """Provide results class for IDTxl network analysis.""" +import sys +import warnings import copy as cp import numpy as np from . import idtxl_utils as utils +warnings.simplefilter(action='ignore', category=FutureWarning) +MIN_INT = -sys.maxsize - 1 # minimum integer for initializing adj. matrix + class DotDict(dict): """Dictionary with dot-notation access to values. @@ -49,6 +54,77 @@ def __setstate__(self, state): # self.__dict__ = self +class AdjacencyMatrix(): + """Adjacency matrix representing inferred networks. + + Attributes: TODO + is_edge : bool + matrix of edges in the network + weights : int | float | bool + weight for each edge + """ + def __init__(self, n_nodes, weight_type): + self._edge_matrix = np.zeros((n_nodes, n_nodes), dtype=bool) + self._weight_matrix = np.zeros((n_nodes, n_nodes), dtype=weight_type) + if np.issubdtype(weight_type, np.integer): + self._weight_type = np.integer + elif np.issubdtype(weight_type, np.float): + self._weight_type = np.float + elif weight_type is bool: + self._weight_type = weight_type + else: + raise RuntimeError('Unknown weight data type {0}.'.format( + weight_type)) + + def n_nodes(self): + """Return number of nodes.""" + return self._edge_matrix.shape[0] + + def n_edges(self): + return self._edge_matrix.sum() + + def add_edge(self, i, j, weight): + """Add weighted edge (i, j) to adjacency matrix.""" + if not np.issubdtype(type(weight), self._weight_type): + raise TypeError( + 'Can not add weight of type {0} to adjacency matrix of type ' + '{1}.'.format(type(weight), self._weight_type)) + self._edge_matrix[i, j] = True + self._weight_matrix[i, j] = weight + + def add_edge_list(self, i_list, j_list, weights): + """Add multiple weighted edges (i, j) to adjacency matrix.""" + if len(i_list) != len(j_list): + raise RuntimeError( + 'Lists with edge indices must be of same length.') + if len(i_list) != len(weights): + raise RuntimeError( + 'Edge weights must have same length as edge indices.') + for i, j, weight in zip(i_list, j_list, weights): + self.add_edge(i, j, weight) + + def print_matrix(self): + """Print weight and edge matrix.""" + print(self._edge_matrix) + print(self._weight_matrix) + + def get_edge_list(self): + """Return list of weighted edges. + + Returns + list of tuples + each entry represents one edge in the graph: (i, j, weight) + """ + edge_list = np.zeros(self.n_edges(), dtype=object) # list of tuples + ind = 0 + for i in range(self.n_nodes()): + for j in range(self.n_nodes()): + if self._edge_matrix[i, j]: + edge_list[ind] = (i, j, self._weight_matrix[i, j]) + ind += 1 + return edge_list + + class Results(): """Parent class for results of network analysis algorithms. @@ -81,19 +157,15 @@ def __init__(self, n_nodes, n_realisations, normalised): def _print_edge_list(self, adjacency_matrix, weights): """Print edge list to console.""" - link_found = False - for s in range(self.data_properties.n_nodes): - for t in range(self.data_properties.n_nodes): - if adjacency_matrix[s, t]: - link_found = True - if weights == 'binary': - print('\t{0} -> {1}'.format( - s, t, weights, adjacency_matrix[s, t])) - else: - print('\t{0} -> {1}, {2}: {3}'.format( - s, t, weights, adjacency_matrix[s, t])) - - if not link_found: + edge_list = adjacency_matrix.get_edge_list() + if edge_list.size > 0: + for e in edge_list: + if weights == 'binary': + print('\t{0} -> {1}'.format(e[0], e[1])) + else: + print('\t{0} -> {1}, {2}: {3}'.format( + e[0], e[1], weights, e[2])) + else: print('No significant links found in the network.') def _check_result(self, process, settings): @@ -183,12 +255,12 @@ class ResultsSingleProcessAnalysis(Results): e.g., estimation of active information storage. Note that for convenience all dictionaries in this class can additionally - be accessed using dot-notation: + be accessed using dot-notation: >>> res_network.settings.cmi_estimator - - or - + + or + >>> res_network.settings['cmi_estimator']. Attributes: @@ -436,11 +508,11 @@ class ResultsNetworkInference(ResultsNetworkAnalysis): MultivariateTE or Bivariate TE. Note that for convenience all dictionaries in this class can additionally - be accessed using dot-notation: - + be accessed using dot-notation: + >>> res_network.settings.cmi_estimator - - or + + or >>> res_network.settings['cmi_estimator']. @@ -572,10 +644,11 @@ def get_adjacency_matrix(self, weights, fdr=True): fdr : bool [optional] return FDR-corrected results (default=True) + + Returns: + AdjacencyMatrix instance """ - adjacency_matrix = np.zeros( - (self.data_properties.n_nodes, self.data_properties.n_nodes), - dtype=int) + adjacency_matrix = AdjacencyMatrix(self.data_properties.n_nodes, int) if weights == 'max_te_lag': for t in self.targets_analysed: @@ -583,26 +656,38 @@ def get_adjacency_matrix(self, weights, fdr=True): delays = self.get_target_delays(target=t, criterion='max_te', fdr=fdr) - if sources.size: - adjacency_matrix[sources, t] = delays + adjacency_matrix.add_edge_list( + sources, np.ones(len(sources), dtype=int) * t, delays) elif weights == 'max_p_lag': for t in self.targets_analysed: sources = self.get_target_sources(target=t, fdr=fdr) delays = self.get_target_delays(target=t, criterion='max_p', fdr=fdr) - if sources.size: - adjacency_matrix[sources, t] = delays + adjacency_matrix.add_edge_list( + sources, np.ones(len(sources), dtype=int) * t, delays) elif weights == 'vars_count': for t in self.targets_analysed: single_result = self.get_single_target(target=t, fdr=fdr) - for s in single_result.selected_vars_sources: - adjacency_matrix[s[0], t] += 1 + sources = np.zeros(len(single_result.selected_vars_sources)) + weights = np.zeros(len(single_result.selected_vars_sources)) + for i, s in enumerate(single_result.selected_vars_sources): + sources[i] = s[0] + weights[i] += 1 + adjacency_matrix.add_edge_list( + sources, np.ones(len(sources), dtype=int) * t, weights) elif weights == 'binary': for t in self.targets_analysed: single_result = self.get_single_target(target=t, fdr=fdr) - for s in single_result.selected_vars_sources: - adjacency_matrix[s[0], t] = 1 + sources = np.zeros( + len(single_result.selected_vars_sources), dtype=int) + weights = np.zeros( + len(single_result.selected_vars_sources), dtype=int) + for i, s in enumerate(single_result.selected_vars_sources): + sources[i] = s[0] + weights[i] = 1 + adjacency_matrix.add_edge_list( + sources, np.ones(len(sources), dtype=int) * t, weights) else: raise RuntimeError('Invalid weights value') return adjacency_matrix @@ -642,12 +727,12 @@ class ResultsPartialInformationDecomposition(ResultsNetworkAnalysis): algorithms. Note that for convenience all dictionaries in this class can additionally - be accessed using dot-notation: - + be accessed using dot-notation: + >>> res_pid._single_target[2].source_1 - - or - + + or + >>> res_pid._single_target[2].['source_1']. Attributes: @@ -787,43 +872,43 @@ def get_adjacency_matrix(self, weights='comparison'): connectivity for significant links - 'diff_abs': absolute difference + Returns: + AdjacencyMatrix instance """ # Note: right now, the network comparison work on the uncorrected # networks only. This may have to change in the future, in which case # the value for 'fdr' when accessing single target results or adjacency # matrices has to be taken from the analysis settings. if weights == 'comparison': - adjacency_matrix = np.zeros( - (self.data_properties.n_nodes, self.data_properties.n_nodes), - dtype=bool) + adjacency_matrix = AdjacencyMatrix( + self.data_properties.n_nodes, int) for t in self.targets_analysed: sources = self.get_target_sources(t) - for (i, s) in enumerate(sources): - adjacency_matrix[s, t] = self.ab[t][i] + for i, s in enumerate(sources): + adjacency_matrix.add_edge(s, t, int(self.ab[t][i])) elif weights == 'union': - adjacency_matrix = np.zeros( - (self.data_properties.n_nodes, self.data_properties.n_nodes), - dtype=int) + adjacency_matrix = AdjacencyMatrix( + self.data_properties.n_nodes, int) for t in self.targets_analysed: sources = self.get_target_sources(t) - if sources.size: - adjacency_matrix[sources, t] = 1 + adjacency_matrix.add_edge_list( + sources, np.ones(len(sources), dtype=int) * t, + np.ones(len(sources), dtype=int)) elif weights == 'diff_abs': - adjacency_matrix = np.zeros( - (self.data_properties.n_nodes, self.data_properties.n_nodes), - dtype=float) + adjacency_matrix = AdjacencyMatrix( + self.data_properties.n_nodes, float) for t in self.targets_analysed: sources = self.get_target_sources(t) for (i, s) in enumerate(sources): - adjacency_matrix[s, t] = self.cmi_diff_abs[t][i] + print(self.cmi_diff_abs) + adjacency_matrix.add_edge(s, t, self.cmi_diff_abs[t][i]) elif weights == 'pvalue': - adjacency_matrix = np.ones( - (self.data_properties.n_nodes, self.data_properties.n_nodes), - dtype=float) + adjacency_matrix = AdjacencyMatrix( + self.data_properties.n_nodes, float) for t in self.targets_analysed: sources = self.get_target_sources(t) for (i, s) in enumerate(sources): - adjacency_matrix[s, t] = self.pval[t][i] + adjacency_matrix.add_edge(s, t, self.pval[t][i]) else: raise RuntimeError('Invalid weights value') diff --git a/idtxl/visualise_graph.py b/idtxl/visualise_graph.py index 3b73fe31..70b7288b 100644 --- a/idtxl/visualise_graph.py +++ b/idtxl/visualise_graph.py @@ -86,7 +86,7 @@ def plot_selected_vars(results, target, sign_sources=True, target : int index of target process sign_sources : bool [optional] - add only sources significant information contribution + plot sources with significant information contribution only (default=True) display_edge_labels : bool [optional] display TE value on edge lables (default=False) @@ -167,13 +167,16 @@ def _plot_adj_matrix(adj_matrix, mat_color='gray_r', diverging=False, # https://stackoverflow.com/questions/25500541/ # matplotlib-bwr-colormap-always-centered-on-zero if diverging: - max_val = np.max(abs(adj_matrix)) + max_val = np.max(abs(adj_matrix._weight_matrix)) min_val = -max_val else: - max_val = np.max(adj_matrix) - min_val = -np.min(adj_matrix) - plt.imshow(adj_matrix, cmap=mat_color, interpolation='nearest', - vmin=min_val, vmax=max_val) + max_val = np.max(adj_matrix._weight_matrix) + min_val = -np.min(adj_matrix._weight_matrix) + + adj_matrix_masked = np.ma.masked_where( + np.invert(adj_matrix._edge_matrix), adj_matrix._weight_matrix) + plt.imshow(adj_matrix_masked, cmap=mat_color, + interpolation='nearest', vmin=min_val, vmax=max_val) # Set the colorbar and make colorbar match the image in size using the # fraction and pad parameters (see https://stackoverflow.com/a/26720422). @@ -187,7 +190,8 @@ def _plot_adj_matrix(adj_matrix, mat_color='gray_r', diverging=False, elif cbar_label == 'binary': cbar_label = 'Binary adjacency matrix' elif cbar_label == 'p-value': - cbar_ticks = np.arange(0, 1.001, 0.1) + # cbar_ticks = np.arange(0, 1.001, 0.1) + cbar_ticks = np.arange(0, max_val * 1.01, max_val / 5) else: cbar_ticks = np.arange(min_val, max_val + 0.01 * max_val, cbar_stepsize) @@ -195,8 +199,8 @@ def _plot_adj_matrix(adj_matrix, mat_color='gray_r', diverging=False, cbar.set_label(cbar_label, rotation=90) # Set x- and y-ticks. - plt.xticks(np.arange(adj_matrix.shape[1])) - plt.yticks(np.arange(adj_matrix.shape[0])) + plt.xticks(np.arange(adj_matrix.n_nodes())) + plt.yticks(np.arange(adj_matrix.n_nodes())) ax = plt.gca() ax.xaxis.tick_top() return cbar @@ -280,18 +284,21 @@ def plot_network_comparison(results): cbar_label = 'A != B' elif results.settings.tail_comp == 'one': cbar_label = 'A > B' - _plot_adj_matrix(results.get_adjacency_matrix('comparison').astype(int), + adj_matrix_comparison = results.get_adjacency_matrix('comparison') + _plot_adj_matrix(adj_matrix_comparison, mat_color='OrRd', cbar_label=cbar_label, cbar_stepsize=1) ax.set_title('Comparison {0}'.format(cbar_label), y=1.1) ax = plt.subplot(235) # plot abs. differences adjacency matrix - _plot_adj_matrix(results.get_adjacency_matrix('diff_abs'), + adj_matrix_diff = results.get_adjacency_matrix('diff_abs') + _plot_adj_matrix(adj_matrix_diff, mat_color='BuGn', cbar_label='norm. CMI diff [a.u.]', cbar_stepsize=0.1) ax.set_title('CMI diff abs (A - B)', y=1.1) - ax = plt.subplot(236) # plot p-value adjacency matrix - _plot_adj_matrix(results.get_adjacency_matrix('pvalue'), mat_color='gray', - cbar_label='p-value', cbar_stepsize=0.05) + ax = plt.subplot(236) # plot p-value adjacency matrix + adj_matrix_pval = results.get_adjacency_matrix('pvalue') + _plot_adj_matrix(adj_matrix_pval, mat_color='Greys', + cbar_label='p-value') ax.set_title('p-value [%]', y=1.1) return graph_union, fig diff --git a/test/data/mute_results_0.p b/test/data/mute_results_0.p index dfb5572d4ffeeed8bdcf4b0142f8dccdb45e3348..e3b5632adac3f8288f58a4e94a93cb8c0355eab4 100644 GIT binary patch literal 2199 zcma)6TXz&i5Z+A)gykv(?{`r*iV+czo7sd7F3BXL45FhF$DWK6vmSczp3qt9y1xR`DF?oTPWEs;j>G>Z>`|UG+UV8w_wN z;y{Xmm7XS*ry@Q*8ctD8!zwL6uKzsb=b)=-_vUgrL8bJ=nt*@aABp6MUzG(|h4DQ_ zdjken>%K#U^y@@M99G*K(Kd!<4Nk;#G*g}!g*EmjG=@%t^14$evyPBNLa(u_!|sQv zeJ$Eml6HJl4}CWlP%wU4Y1U~dz`Ay{dJjm=k%ZT%g!M*`PM(=X6!R(-upzT!H5H!O zHyU{zuOP;fI_SohC~RukYY@u}Y;I``QyzD!5$Cij{V0U{TJ~lPBA$oV(DRj% z1Om30b^+~@#2KuuSY}<@l$Gv7+j_JGQI!PF3<+Yyxot_4OM+U&eaY%h7}w!`j9iP6 z(m=|I=R0-{C*W~S!-?ZZoZ<0tctWT5s2kVd$t6t*2^x$n5BI4` z>4)p;xuL7;$;)Pg!qfI@M8-+zMRj<_-i9v($!bbpbuVJhqZ*OO!1j*3nK+(Z&a1?i zR?Q^x+`XKcP@Z3QDh|_|I7lciI~1A~vvu5NHn{$lo>A%yQNN*t7mTi$^|)F$9+AXp z@CX}~KB9k-sTJ8#)U<_870~5dGwNUQP&x4m+G27wlXh+~K z1~fyma}Hi98abVRS5Zf_mE<68l6vs!BeJ)|_S zJ;(?){RYF{DvBQL?I<4D+bAC035?-^EgFw)9@%O=-im^p_%YtJ?qg!v=A+Su-9Qz9&{m;lrIPR8v6?FB$+XwK~eEc8;qFlpo#;KOp> z1!oaHvVdrLEP%^W&{O-Mnu7K$pa`1%8}5|5RNj21?`GxB>8qvegqj5AvJyMVbPbGBJ}Rl42_e_Uz)+JCzQLV|Q-HsZ!wSZe?98!xDE^e08&F2s{C zqePIW_7)_wbC(>o{-?BFW-&njVsZr}Rm-ebBV;p+bx)3++! acNSc0kLml)FJ*Q31?4A4Ed!OgIxiwhvJ#{@O2qxn|$y~Yn zX)|>#5^D=#l?fJBSPnABZbSA44CI|7nM{T>L{wo?GD&DigQjF246eeiyp^%C%o)KD z$&~WLnuMW$LCi;L-mj<}?8cPCd1pTc)1bjAiRz*6#u7$K#ovF`&u`8Fd|9X0NkD5v zQBh+G_Ea_QSVa-HnZC`n&7F`mn5g$M?zR4K4&t5UxbBBUDym?gb<)|> z!^e!UzaRS%5I+jxKwoxUm}CG4G0QM!S!cqH@ZeC}+1D#DM4(W?6V`yK+Kfy26&es3 z3lY_*ViFE_k}=dG;VWJzVO)nJR@TX39!2A{yHPyan;gj~7RUh{#j_zi^Ozgg;K`1` z4kC@PPJNRS%ZXy5XA+f zo`mO2mzOS{C$o!-@VuF5N={nM7_zg&XgZ-6cBZP!4B4`=bLx`nz(bu z6tgC%T#MGp8ocCWF?14MHs^{nKVX6=c*UGgz^jJnu5t{HCiJx(AkAwSi8>7s3|==g zs$$y=-Y}KxwqwNzYC-vR25+`a;gz~in8XvRA^~{=asf?ZbqkwYT7YAHTGArW&~Kp% zw8pf!v|Q7PylrJ;KI}S&@g0$9iGU?^5-qrQxc0$ubdLE_$U6fVs2GgfMlg~EsBq?` z`A7pA$c}g-mb3$yIo_+aIcqg{rF!m_)-4wcwO9il99xb2G8q zd#muiw)L(SOlfzw;52>=wBQUg>P$<(il_UFGKJ3AW`fBu2#Othfe%*UY~IRR1J%#mtoG-&i8alfLr6kbP;?KdyX7bf3SkQfL!-XhcXG%r2@t|~nbdw`L6-|1VA6`tm@^n|M4ArEo zn(3-H1FirUI6sGf$6aFvk1DEvEYNJq$5r&Tp>nNjs}On9?xcDlNeGty%$f}ok}&cx zxQ@F-FTkFQJqX(|UvrJ5`n?3v4QQk$g1l~0dBk-atQJ=2sr!;@cM7EW@o9 z+(s_=RX&~$^*f#Q`!f8HWZmi?YD>Yk`*@z(|Ho5Dr~I)DKP4&e^r!q;!CmAicH`JB F`~&*crd$93 diff --git a/test/data/mute_results_1.p b/test/data/mute_results_1.p index c0020f20e50077557e6015a0578807042d044682..a3c1249fb19a1450a87d3a7b813b15f9fa46f324 100644 GIT binary patch literal 2525 zcma)7>30)V6i-qLC2SRn3%HA#3RXlwDA}V;P2HRqF zxy17==B4HET_9|;;|{I22Z>AUB1uD{b0;LXLn@=RB@zjbFlO2%56(Ql)?sVSB1^+0 zjBm{-GcdSVHdW#=v#c@4g~`fHwDn?H$<@3eXE@it9;PU>&}gd_;+EC2wng=r#-L4F z(quR7&^{IIJV_Mqlxhx;szt|jNs*b34d>OAb_~)C z1Fg1Ua-(K>FuS0n(atnKhBXJvOp8o$=?=6_N1NpoHA@}VEFW>s9TVhgR>^TqM$4-0 zm*ISjoQjb{A+_odT%aWJ`%<`2*louxYnF*PFH&03K8PR)%o9FFiD?m6WiVfySqv8o zB({U`szrv01()P|)>ZfZxp;SX=<>@dqb|`54=xonY}&X^UDMYGmx=TS(ftx!J|-y9 ztP0f{!@VM(J$#TqheUNTc~-3uxKf#n$Xv}foHAUc%*8Jnl2v6U@80uW*C}ZX8CcMi zHx|d$jl5iZVbxe7*Zh|=7Rt2^C+9FciGzep!@-#Ee3xjrz42TqmV` z+NulP7#h>mitAvb%tZ7HDPNIA89|$^+Kz$7>rs56vXNArlg5SOL0A7DNg&Gh*mw^MT^8p7H@|eGg4CIUy>0Kk&;xTjS&NG8i1u4DIuk#7HJX*mlSy=HVf4S z%cP{#8o+X;1N-%eg)$}v`dlwAa`OPJz#`J509N8F89+C_S^`*wuT%i5iI`mT*;vrG z3_uUoFQQ~C)h^r<(pARm$JAU`tEN4iv*xBD8~W1*mqi+GRZ z9FtbTjf28Dz^b_7I3niy{MOqLLAH&s6}VGMjKE#Fqz>{ZY?R?{&RJ>%Hqo%1!4+^% z*v`E&DA*1?yOI0!v0VGMj`95It^2d_8}4WDfUeBJ4(c4ofEj@Y=>{rl4^S-pkPK>6 z_+hGGw!deL_IcyKWl!jT&8yBokOjPLbzMj-S{0NE6!Y1WoeLJ~yN(HI#~#6lkun`;zBg&L*HP^_p(hobN*_!)Hv8bu+w@jzR|IB}p}089oyfMJamVgQ>tfGsjq zBLG_kfI6Go7(A-eI&0f`^kXtS9!2jE(MP-cdiqPJvaplE6FNnDPL=w5`nr$O1{;*r zXuazvd751^j6`XI<}{-WYS?vI_r)Cf?SQ{P9`RIh;d!{JV*~IMKb#`UZ0xUj(^Y;s zJuL{^jX~nW9>Lt+#&*+XbU7z&Uxu!v-EMe2VEJ41_;2lA6 zq`-4g{4P)To(%8vbSSJ31jTh>9cJ($QoJV?))7I^(TT#MLQ@~nP)_4X)K6Xf(cFZ+cIUq)q*iL!M$eZ}DG2|a(qqra8myD0j&h_1`&dj>!JA36QV)0~ju Xrzp+O&1rsN@GJIu1pjvHzTWjOlmK2j literal 2389 zcmZ{lX?qh@6vva6LJ6x-#0}gvpp=TJh-kzDu}mpCNFA4Yy_31g9XpxHxpOBF5CxSX zRm2S!aKl{`ao@gyZ+Y!2e9yxN;Dw%hGm{v^d7d`S%$ak}|9Af9c5iRVb>vQOP_Qua zWSBEspJPn!2*l)wUuJ^&C6B)pVXMBwlmH${;=lZJIU>;>yuDY z3A2q=*35Q&%m{P3YyYc=8~8A{Y0biLNh1;lkw646&&XOc@Lrt>I|v!fZ(6hQYwJZH ztN0QwGJ13w$(idhK_o1&vKTZ33lsU3K)8}ui62#A05j&2Dv9PrcuS#3`$r;J3NFU> z`|!5I?5F~Zr*vg6C{d5psTW~)mrPmnF@QP_vhp2Qq6+GTaOs$CNJEm18iUKM=@=)d z?*vu2T&u)hh~$D%Cxk1s#_iiSlJ%pbutdX5=aQ0Gii7D4Q$1BZUfBVv(D{KIGUBm1 z^WduDJsTT`{u(;CRyQ*#XgC?v$s}B@>8yrp^q{rFq^negkq|+JN>q03lqS-0jxMGK zdR?cQnyVj(D)q1;T(3cuQ9FYhbo)wsBpE?3NVm%1#-{GRRCNgp@eY-NfV?iU8JEOp zbvL|a1eSGqN$FIjm*X&$$CSC0T~jPK8QF->u&o97jD)xn9!tzhl;LjXY7bW6=IB2L z^HvWAN(L)SIDlGX5IxQYtGJym6@oT$1RWosu|%YnXs8vH^PSdHxW&lqhgHK_8*V)S zt4+8~?PPX8+@1tFYTP-vLzgM9fi*_P$R2<@O<0>ZdNufUaM#$_*yg|S(!Y$1gmw1R zaQ~9>)_DBHnF8D`;T~HLCuANjAahc=q=wFe0w3bKk$X9IvEGFHk}ih%TFkS1;AHXi zuEPb`AYr4eU7PGa0;>41so9;oPDPj&fZ_z)Z)E-lPKYj{A_pV9CD2_1n_B`66Gl-W z|F5(y_He&qt}YmhVXFxb#PA?(*#RBLeA(J3Fp*_n4BM;$jBQKT?IoM~&;)G9plLCT zQ@zEhFr4Va#lo2r_lU4kMBtd)QRMX-{zZ|9ETPd89u^Chx>aqX=GiC0Zb4vgp`}8 z4arK50HZv9L4tMs$i}dPaU^%@+|uP|7Z30bD2`NiDk@n`g}2o=R8zZ5*qt=>bhm^( z5}rZzHvX)z^~yR_QuQwlj|A!B+48KWvbVMTJ}#U(L#R4CL)f?F&(?wO6F+b;_&GfG z)N$QDYaDpK8S%Zg9%T2%!wY+0KX#?_$9X0803W9vo=w&IlnE~;t-sXG^Pq&6Q=W%X zo^9{G62q$|ycWaj$dhmD%=8VF^{@$VCRyLoS=)AaTf#fGVw|}1JjO>l&g(Efs#3md z!h1=|_q$VmAmKyg8>bFCrnP>Q=XH+bY}@B!mE{u?K25TG)-CsQ316ggzf9$}?Q=Ya zuT1zlh7(Gz?8yB_Wj$%aw@KFTbk??oznAcXmV04K{dgfu{iIURmy4Dpn1rgxe3{A$Q7*7OeQ0!lwQ~r@XzlTBDvr-WFCewexzU>#o$KZ z*;Gg`ATr``z&eJuWlU>uBG#hya&-v~TF22C+AYchJ0KgjkVL|$aelz=g^B$T+7(DI zJ_UdRR zEUArdY>C40u5}P=S*P5I1ceh_UBjx&?MB2oZAdQ);i0Z|0)vR_qBV3qWhB0UlQk=k zc1hy2tW%igaL<&c9z)v^wE0nk`1U&SW32P^o+g+0&4_!F1$Gz*@GwRm!pOCRT6+y1 zv9kC+36E-f81aDk9@hDoHG=k4EP}%0+9w(He9CPJPv~zZ;7MIc7-78a({<{@Q{R^S;8?@}jj!n1pt65_WQ*&pt6mC|p&sOwjq zT~A)NTNIwR4q#=Rgl-hT3)X3TAxKtRda8R7a~?H`LFzioH#0YW#6Il^-a%q+UaNtT|I@w8KOc%31^I4%tjoo43|h^w|IoDNDu3Km8m5- zThPRWb{M&6d=2H5XcuY2Nf0Fk0k2p8$Z9D(n2-zO-I$$qbhIsS+kBc<^2Rny6pW0< zKUdHZJ+0VCjHE7{ufv-KBWvW0A!C>(YlzIzP$v3;Nh4#7bl`$DhW!faqfn`fJ}%NE zZ`I*#OkxaoUbrQ)N+yq-K zZ=|Q)!XW(sdhC2({cMt5P!ksgHpMQoX{shJ*-VWUYphgbvo$ui$L3KwkV+@k80-Dg zVOllUSR$AWOt!?TeM4<`VAF)H4t$M6P21mA?SEs!^+L`+qp!Q!qc6@M90=)*#GJ5^f{sG5%Uk?2P{n=<=pZ literal 2052 zcmZ`)*>)R65S6^dji~}We zrPim-)YV9=%!SKLu&~T>kU4n+vYXJAcaCQ=8B!BbnMuhcp&nLnRWv;#CsHRXA>CogC&-G)}u4z@r1nk&I%29KZ=Y>%}vVxp4&^ zZyM|%(g>?Gh%w@crgIbnsORCpLeE#&3JoMYX`ANFq*T4e;3;Px#t9mFQ5BvxNLUC_ zTrlcMc*b;j@xm;bT3CQ*%|ugj(rTVVcJ>%eC-nT@RCW1b^3x$Fm^Fi0@GLt!<=R;+bW34`jwr7J?Tv$F8g4S~cmq3I@6S3YaQ?006 z=(XO$2Ucbqrl+$m%m%z$f-9~PcuwdKv z{NMPPe^y4pqT8JA-x|Ap?aqUH1+W!-?DAng#SPpi9I-wvK`Ftyz%7j47`R)!zy3=B zK2vbfmw2l-e+-;fZC{k&GS(F69%&!o3G~uH7En;*-UKTEAbsJFvl#pRy)BYxv=Xv@1_{ z#Y(W8bX753wP(N;-~#98@$Zyt%-~T)^^XObP5HQrzA;p;v}_e3Z`z$yFC+=U;-6Ww zVL}o{9tKx&m*@r9cCiQHZp>F)BdJa=L9_|A)I^ZiO^Q4p<7fFeXl~1V1GyZ+52qS-~yj KD0bu69sCD11gm!d diff --git a/test/data/mute_results_3.p b/test/data/mute_results_3.p index a17374e8e847ace3f63b22f755513b04c4f2b22a..f85481fc144c0e184cbe5d02dcf71ee25ccec20f 100644 GIT binary patch literal 2787 zcma)7S#%Ud6rGtQ5QbFQogog4 z4((WtyNC8MQdIc!7F|XZyvC6TjO|Ge$6h*!tyvC%aXq4kR^3)ImTi*^Gb|I%=t+-5 zCso(cYU&0zsyPS7r_yn>Gt~{-nt)|S2Bxrd8QMmnEoWuaoYJA@T;N{k$_6$a;t znCWo7=)|JD5eU|6!aPN{Bcf{R3YooH?r zF76eOsJV8k4v>3EeR9_}{@E_73+7p8JAq5n!>}`3HFYZwm!&7-7Zt(kWCow!ac$em zstggBG$?N`w~a`E|9Ly25*vdd5?R}Nfq4E;zX5)uO!yqf}U6|FNQZgP*{$iFmE zzb1nzQp}|#n$V4|GF54}EgX~?*!@(>S7cg3z-B6@rK9mG6rZnb1QmHwxlkM^Z~i?h zgwoHS4p;ZMbebl-6$f`QN5W2~uYzk5QbgpRkq{n%k`#oE5go2=ftd*@D#fHCshIdn zintP!`RamMQba28V0O9;=XFR9WlVIuvz@TWbuBOli%7*D%*7+>K?NQ$59Sd;+2Ej1 z$=A2Q4TG=4&e~DdP?<_kC@dIWG6G3|QvxD!5U5HCz(b2jP(9e@7+DkM_oic*6eDU< zE}fxCQGrZtf+7rcdL!(TRn#dXIp>=&eniOkq3=<;D{2TI7kzQW^n1+ zJIis?NqHV0P(!tCwKMLZQdK(+&HjB{=faLdN;G23ftxsD-LO!GMG2`$Y}@r9r&6iL z|L`MzNf8E%wVG1iI{(1U9xRdJ77uPkYMbezno_>l`~ztZZj<455AGne@)Rv8_!Oa( znCixlaCJ8PmEDBWnUg#inFT$U`XpO%yuPn&C7s_xYr=FEPuDI3@O05Guz9-Dy=%3plO`WYg2RAm zF0Iw1m>}I1q<3)AkxaE5%{{=ITPs6XFsmnK4Uq{?=w{GUfKV@wdQgUkf~a*OYJoa-Pw$As z$%wkSW9FYblkhNu^#Y+|vktD3Ck#6kLK}FxjWRqEqH6)VVJD*SKg5JwS@uY<3~<`#I|M=nOv;}Mj>?~cNs$}IxstsFc4yloSB zZx^KlyyE{dhWB_nl=NPnZa~uSGxz|B+7U{6pMd8> zpBP`#RH*qQ>I?bfeq8#6{D}vj%J7*7pAQl87an{m!+sCGD!}Dy-qHaXz6mD$_T;#H z$KYTAF5mO0A7uD3i26xH4H5Do28T~a$VYg(qcZ#)r29pr8zSUm41Ps?_TkS~&DF~P E2TH%YMgRZ+ literal 2547 zcmZ`)*?ZJf6rX7+lwlQ$xPiM4D5YT00+j&^RtN+Q0uim@W^y|>G?VlsH-)mJ5@r-p z5f^aZMJw+6zI~ssJo|+2eEb1C=(#t^41+M=*JdU;_nh-v&dJ{9qHXbQPPTRneZrQ1(^>K17ZGPCwUb{90|^{JUm#u)UyA`Jo>7{n#ca6m0+?t@WzC8K1O zCaSk$P(V4i-BJK8=M`D*@#S`r=U_A@9h285U@&n888N)FYnx#Jt-XcgN5$uuJ`o?5 z$u`3wC4&>cL^+IATJ&jnyWjGSqUZax$PIHieVVJtGp;d6{cl}nXiwn)yIr5p|SQkMQqQ7i4}buhI0}cffxEn0!&h}`WU=7M19i>C`_*CGJ>Vx zQf$8kZ(Gz1OEA4ED?471IK~ihLhSCcsy-P5h-D#Hu4QvXL7V_CS7k$zBtASy;R<~; z#`(mxyfR!VQDQF!XFerX09Q$j>({L?maSe5GbEX*E-4c;QJA`z;;HEI>N=r;{^>ZlNH#OzlVHFbG6vvU$xbI+CxwRf#^J4mzA#X+C2RrG*YvYMo3BrY*TJuiAKEi1@@K{1om|GW(&YCu_p;1F%AayG7~je&|gCEinN(xJNSI z-vuiL&I7PY10}CCOD}fAYE@NN{fi&@PswmtV^(`wXY{X)HjJDufXZR5$vW6FCd=K+ zu-1JVXi2SgOu^_KZ6|t9Z9h_g`#G#P*}MuvT{bA0i{J{-^Tf}=Mi%#GcYz-F_J9VP zus8Pqt{yadT7`qe2Of%`PlJaeFbL~J&?xn0hKdsWj)IzmE$%Ix9@!AB?VrKAdkcK~ zsooQ9I~eLy6r#uautJK;5~!#K{Ryg)L{;N}ali%6mYPbFlxz!1mL44tM;#5y$x&Ax zZJT~(_RO;<3*d1WJcm9I(Y^*j63rXw3pos>^lin-*~m8av^K|}Gq5d!VGXuNutT73 z5~z0y)VuQd17YvR3;K8u73fc%ojLo=^tSWR?~UJO!ku;W`$SWZXz*y#)MJeh_H&33 z{TBS2W6D)v5l+O*5?2kkac3QntW@Hq;-YKCZX#kQ2T)z}cl;u84cD_U_;FmP;x$m) zVHP}530cLI1J$@m;NUKJGF>5u@&fErdDg`ivPG12sBB(U80~4+$riA=DqF0wZdHi< zj6jl#jKg*)9k#zgub%~Ef1c-y?VOW(PWOWzeK-_zjzB<1nOlpk>T5H&iSx^zTh z{V30NGqgkv?g^3QV+}q@vV7VI_cIQkr*OYW;npnpWdvVo@O1>=2)H}zaK9B37sr6rY(=C<6#Wk;Q!<6icZTpZGP2cGC{+qtVWjMDb41Hg(^Fv~+T?S)?HVW1?vO-O>sw)7%1KFjlG+$z!wU_->wf zFfO)ZH69+@&y$iOoGt;r*{IW6u zgEYfHt8JLvs97FN>{QZdXPO_wnuKM>M5efOJ=(^i&2sXZrS@r-k2og}2y!*6;J7BE zMb-9;Z~;b+#>k$KTD1qJC`tT22c`or|aa-%s(vCkGx3`2Yzo;_m63y`75<$bt6-(8{?d@=>NN*6`FTiC3 zf)dT@rCN2km$zjP?&Z%uQC&=)m3j$Wp^QXiu4Wrf5w27w;};FdDlwCH@Aeps) zos{xvtuAz9XiQUkT?ZRwCZfNd@)enx5wzK=?HFjh0mT<88%f1EX8-lo1h;l2oLP5d&`Qf;kx}A*G}mX#@$E6nP~! z3)Kb9Qc|i7V6IY+{d&Yg850A2t``@%sSD;|5otsK^YN7opaoww0W82*Du9JVOm67e zSkO0jK`YkJqhxza4Y((yp@`RyskyFJN_#kG%}ql#^rsCji!>}!QV7PJ2a6GTatqud zLoOrL2(2!JC7qp}%m2Yg{F4$4mg>1Weu6)6YXEIB+!nwxDM7O|m#gDz!5>&2fFwhE z03ETbCiEnrGl1J=xFdi&&*~=Tvlo{Uh(FS-`uq zTNi4KwhC1Rg|_UeY1Xf+nzxEAh&_TI+LmSO&$Jz>-!R1InFSlnvY0Ancp!$;!{Jym z6giw^KiC{jUDKNF^G{4Ul?8`EukO<(T}lasLx!T9T`p&&59LgoQP(J5kD^8yenyG< z8b!_ZoipII0j!gu6u^3o$}xZq9KeGzY>WUrBml(h!ov(U>9oS$W*+^B439?9kBR6* zE$yuxh0|H+XAtNV={a5KXl-vfM5}C2vrH>pZ{camG7Lm%whm3Rjlp*88gyUGlV2qG zW91P~6?dbDOZ|8k?BIt}M45;qA9^XM{2qHk5VjM8#D`sixhJa^0d1ypIcd8ybUtmN z3+O_^Y1@Ml;m%an>8Svomf@KIo{e3Vb=n)iJ{g`1;Q829S*I5Qcu|Iz0@zQgY#-o_ zzbwNm96QSYRo?jUntqMJ>x0;TgP(p=hJ(@Rx5VjTwfi=MLt+|{ScjX&VV?Gg4Dax? zsKKLR8Wjz`%itJJV^^%f@YSH+zvj^=WcVhE{#Ha+l=mHj@BfeTe&A_N%J5^9=BJ@)erE6s_Ie!u JcIv*~@GodUeCPlG literal 2389 zcmZ{lX?qh@6vva6LJ6x-#0}gvpp=TJh-kzDu}mpCNFA4Yy_31g9XpxHxpOBF5CxSX zRm2S!aKl{`ao@gyZ+Y!2e9yxN;Dw%hGm{v^d7d`S%$ak}|9Af9c5iRVb>vQOP_Qua zWSBEspJPn!2*l)wUuJ^&C6B)pVXMBwlmH${;=lZJIU>;>yuDY z3A2q=*35Q&%m{P3YyYc=8~8A{Y0biLNh1;lkw646&&XOc@Lrt>I|v!fZ(6hQYwJZH ztN0QwGJ13w$(idhK_o1&vKTZ33lsU3K)8}ui62#A05j&2Dv9PrcuS#3`$r;J3NFU> z`|!5I?5F~Zr*vg6C{d5psTW~)mrPmnF@QP_vhp2Qq6+GTaOs$CNJEm18iUKM=@=)d z?*vu2T&u)hh~$D%Cxk1s#_iiSlJ%pbutdX5=aQ0Gii7D4Q$1BZUfBVv(D{KIGUBm1 z^WduDJsTT`{u(;CRyQ*#XgC?v$s}B@>8yrp^q{rFq^negkq|+JN>q03lqS-0jxMGK zdR?cQnyVj(D)q1;T(3cuQ9FYhbo)wsBpE?3NVm%1#-{GRRCNgp@eY-NfV?iU8JEOp zbvL|a1eSGqN$FIjm*X&$$CSC0T~jPK8QF->u&o97jD)xn9!tzhl;LjXY7bW6=IB2L z^HvWAN(L)SIDlGX5IxQYtGJym6@oT$1RWosu|%YnXs8vH^PSdHxW&lqhgHK_8*V)S zt4+8~?PPX8+@1tFYTP-vLzgM9fi*_P$R2<@O<0>ZdNufUaM#$_*yg|S(!Y$1gmw1R zaQ~9>)_DBHnF8D`;T~HLCuANjAahc=q=wFe0w3bKk$X9IvEGFHk}ih%TFkS1;AHXi zuEPb`AYr4eU7PGa0;>41so9;oPDPj&fZ_z)Z)E-lPKYj{A_pV9CD2_1n_B`66Gl-W z|F5(y_He&qt}YmhVXFxb#PA?(*#RBLeA(J3Fp*_n4BM;$jBQKT?IoM~&;)G9plLCT zQ@zEhFr4Va#lo2r_lU4kMBtd)QRMX-{zZ|9ETPd89u^Chx>aqX=GiC0Zb4vgp`}8 z4arK50HZv9L4tMs$i}dPaU^%@+|uP|7Z30bD2`NiDk@n`g}2o=R8zZ5*qt=>bhm^( z5}rZzHvX)z^~yR_QuQwlj|A!B+48KWvbVMTJ}#U(L#R4CL)f?F&(?wO6F+b;_&GfG z)N$QDYaDpK8S%Zg9%T2%!wY+0KX#?_$9X0803W9vo=w&IlnE~;t-sXG^Pq&6Q=W%X zo^9{G62q$|ycWaj$dhmD%=8VF^{@$VCRyLoS=)AaTf#fGVw|}1JjO>l&g(Efs#3md z!h1=|_q$VmAmKyg8>bFCrnP>Q=XH+bY}@B!mE{u?K25TG)-CsQ316ggzf9$}?Q=Ya zuT1zlh7(Gz?8yB_Wj$%aw@KFTbk??oznAcXmV04K{dgfu{iIURR#a&Scg8>ncV8{Rg0=#6@XKklTe7?!N%pI8}=e|2YptRUF zk3=cjYTc@>wrXn&)Y{tCu4*^CFCd6^0c#h!+D%()i{H8Ly<{ddzaR6*x$m9rp7Y)B zo;N$HyR&JvCD$Y!x0F+Ez3HzBM{O?3jq3|Nj&us$PCZ5D-aygqR256sMIsT`Q7T)= zxHK|;by2O(cB^`-fq!i*IT6m?`K;x*Dx0@eQPP;?B7&h<1lIT^}= z)$7Q-mAAK8uCf)?84D})*@CYh2R$bV{EPWQHdS(|-gvrPt;Iu)Md;JML<*4jujl2N=bqM2G>D9XH@%c9O%$y(@p zP=rIXwM_{pn{%Y4=rkR(fKJzy6pHY-a?U0vM`yIRuIbzT!h&7Pd{ZfJDIpy@?b4YV zhHI`~W3BAyptE%Pw38}j=5`Zc(fJp~v`4@s z3$Q6CO_z4jqL>jis*Gx*#_=yHaCHjGXBRCtB1Wx8OOn&juj}MMV!PJ_ zs`04V)@M}d6;gLu7hPUy4om3;kEtb*Y$)F$G%AOd`&T14IaE%q35{^5(Fx*}Iyajs z(2Bui6_O%iWuhb&ri`eDVpU9_?=~P4^eq+x0(7}rpqVv{2fA5{CF)T-mTHgIVX5+H zy(7@b284Fwmi6{E;c7cgc}&c#Ev4O8?}D1PtS7^_*Qc4u_0*vY$u6Z%Y#QB336rjf z8P)oFwa}G`L}J5n{IoD4impnv)p6#L8$EiHNlB01 z9Ga?qON8Fy(OXS=n@3%aXc?ti3u6hiQnJL8`r4r8LfY?zEF;oS*R@4cWbdSuN!<)h zbw8zr-%rm{a{PXJOv<1iabhPTHP|+D+q8Y{N4MS}lM5*Ba+7w|5V-W?S19$i3)2K6iQ z%23$J-UB8L2HrPn?>);pT01ib;&hXuUApT7naO@Z^FBrI zN7tKEI5r+Yx9YBM^Sefi+;uf~{Q>U!cAPF>hdZ#~9^8oqx8g2+9qvYq&&wgoe$b

`TDB@*;oPqpz6sRgb<_Mw=(NpeIdwiVFg7zphVk2#s$j`eqp#-(v4?oAh+x z{T=O{Y&>@3vE4`F^j$^Y(&tdnv?RpxglJ8OctWfg6g<@rhEx4H9sWFT?0?M=t_S+hVq&dm7wd#SynoL7Y4C^G z`U{VKY0?poj)tb80sDeSzcT4XkA5ARh6e0!Jo>Fkzw_u=8Iga_=>5T@KXMxwm49ON zMxTFwR`i#$BYcUY|7z0Tg6O~N=uzk2KNS7*Wai(?oc0xy{uLbQzx9z0&A(R_{Rc<- zTsZ$;(@6f;Kf>T3^yvGa@TcDKkyywG0P|!N(v(&Dxf*hq%CZKps6y6imobojEpjXy zfKgpYYsp_XCC5PyBO~kC4e{f#nsNeU5I>O(93Rd{ISCpiCqtgZl$29A9*IxEYRajQ zLE;8BR3x4Tjgr$Lf!1rm+?1zs2>wvW8T_EX8UFMXJUHb{ZN5X?DfF{l&VxnCiy=>9*2*T%gAL|mHRS@xV1rB8z#EL5 zyK*5cN?r}HH6rl5hD zN+DNvDlHL(nQ`I5xZ;CI%1LG*#WI~`aM6g$Q*kE;`x#uE7=y=!X{VelT%y2ZiiYHZ znJE`$D(E-e*kP>f>Viv^&nzaFF?$OHU8dERl$v%a^~;Mq6Qj11cA1f311tkqbhdR2 z?l{?a=TenSpD^7%XTa!#D-~%AVOE0buewmC39OqF!s#_70={arv~?&Lwa(SWXjE?7 z5m_^X5n;CCuLmh%aE&?{r0+8y^|McA(*|>qBvU#9<|xy(c#*VxcJ1rdz}&JH3qoIu zHKNv|o-3v0I;|?l$6AT0xM#R{VKU6`VT1bcdQP{b0dJT3m={e{qoHKb7+2_Bmchnn zseopqA1e;?nUO|x~)mfFdirA_A zF0^Z*?XbGN%7TtzxIu?b3SG4wUgsO8=mtb!jbgs11=bSIVbFBwifUoii>S;TFt|zc#D4SnYMxLG%aR9=uft z!-Kb(yct1bsZEX*MNX-pNt3sB#!m0q;H~eO$(MJ=zOmB0$9m zRgVt6J}RzI6*x#sxYf$>m6lemO4f%Z3yx;!XjTW?KkBHXLoGYE?K*WN2K^GCgg!ym zt`5?-&Xvxe!vvO6{uJG1+S@xw8=Q^H*(EbO4u@6lnWfBL=ZbkIbVGu$oV zeE}=?kf1193e9k@2lwgleh)sN)#SJ6X`(;4--8Eqc+i7~@Up2>g{@s4d{BoEdGKL| zWQK{%hl$MH*nXr2kKlxA@F-4U4?e=Eqh2!3-2B{*!xh|qG|w&8poH6x5yX${@R$$r z6AEI>)Wx&p+@rD3(T&Z6$46tKw&9brU${F4pOml%xgNF<0&)F#3DnYO`anpftsMI9}ZOgdFA8B9HpqBokt1&F&$p;;lEg3nwKOT4@&cs z*L3+$J@}apCpMPdCTM*%8e{S&$PvkrgpUHr8i>)#~&J-~ViPhOnJyNTB6Tl k1-ZjGD{ 0: fp += 1 + if res.adjacency_matrix._edge_matrix[0, 1] == True: tp += 1 + if res.adjacency_matrix._edge_matrix[1, 2] == True: tp += 1 + if res.adjacency_matrix._edge_matrix[0, 2] == True: fp += 1 fn = 2 - tp print('TP: {0}, FP: {1}, FN: {2}'.format(tp, fp, fn)) diff --git a/test/systemtest_multivariate_te.py b/test/systemtest_multivariate_te.py index 9243a3f9..68172d94 100644 --- a/test/systemtest_multivariate_te.py +++ b/test/systemtest_multivariate_te.py @@ -209,7 +209,8 @@ def test_multivariate_te_lorenz_2(): # Just analyse the direction of coupling results = lorenz_analysis.analyse_single_target(settings, data, target=1) print(results._single_target) - print(results.get_adjacency_matrix('binary', fdr=False)) + adj_matrix = results.get_adjacency_matrix(weights='binary', fdr=False) + adj_matrix.print_matrix() def test_multivariate_te_mute(): diff --git a/test/systemtest_multivariate_te_discrete.py b/test/systemtest_multivariate_te_discrete.py index a19431f6..948048b3 100644 --- a/test/systemtest_multivariate_te_discrete.py +++ b/test/systemtest_multivariate_te_discrete.py @@ -227,7 +227,8 @@ def test_multivariate_te_lorenz_2(): # Just analyse the coupled direction results = lorenz_analysis.analyse_single_target(settings, data, 1) print(results._single_target) - print(results.get_adjacency_matrix(weights='binary', fdr=False)) + adj_matrix = results.get_adjacency_matrix(weights='binary', fdr=False) + adj_matrix.print_matrix() def test_multivariate_te_mute(): diff --git a/test/systemtest_network_comparison.py b/test/systemtest_network_comparison.py index cd30893d..6c9ebc32 100644 --- a/test/systemtest_network_comparison.py +++ b/test/systemtest_network_comparison.py @@ -85,24 +85,21 @@ def test_network_comparison(): def _verify_test(c_within, c_between, res): # Test values for verification - tp = res.get_adjacency_matrix('binary') > 0 # get true positives - print(c_within.get_adjacency_matrix('union')) - assert (c_between.get_adjacency_matrix('union')[tp] == 1).all(), ( - 'Missing union link in between network comparison.') - assert (c_within.get_adjacency_matrix('union')[tp] == 1).all(), ( - 'Missing union link in wihin network comparison.') - assert (c_between.get_adjacency_matrix('pvalue')[tp] < 1).all(), ( - 'Wrong p-value in between network comparison.') - assert (c_within.get_adjacency_matrix('pvalue')[tp] < 1).all(), ( - 'Wrong p-value in wihin network comparison.') - # assert (c_between.adjacency_matrix_comparison[tp]).all(), ( - # 'Wrong comparison in between network comparison.') - # assert (c_within.adjacency_matrix_comparison[tp]).all(), ( - # 'Wrong comparison in wihin network comparison.') - assert (c_between.get_adjacency_matrix('diff_abs')[tp] > 0).all(), ( - 'Missed difference in between network comparison.') - assert (c_within.get_adjacency_matrix('diff_abs')[tp] > 0).all(), ( - 'Missed difference in wihin network comparison.') + # Get true positives. + adj_mat_binary = res.get_adjacency_matrix('binary') + tp = adj_mat_binary._edge_matrix + + for comp, comp_type in zip([c_between, c_within], ['between', 'within']): + adj_mat_union = comp.get_adjacency_matrix('union') + adj_mat_pval = comp.get_adjacency_matrix('pvlaue') + adj_mat_diff = comp.get_adjacency_matrix('diff_abs') + print(adj_mat_union._weight_matrix) + assert (adj_mat_union._edge_matrix[tp]).all(), ( + 'Missing union link in {} network comparison.'.format(comp_type)) + assert (adj_mat_pval._weight_matrix[tp] < 1).all(), ( + 'Wrong p-value in {} network comparison.'.format(comp_type)) + assert (adj_mat_diff._edge_matrix[tp] > 0).all(), ( + 'Missed difference in {} network comparison.'.format(comp_type)) def _generate_dummy_data(data): diff --git a/test/test_idtxl_io.py b/test/test_idtxl_io.py index d8bdbf3c..7978b0de 100644 --- a/test/test_idtxl_io.py +++ b/test/test_idtxl_io.py @@ -97,10 +97,6 @@ def test_export_brain_net(): node_size=node_size) # Test input checks. - with pytest.raises(AssertionError): # test for square adj. matrix - io.export_brain_net_viewer(adjacency_matrix=np.ones((3, 10)), - mni_coord=mni_coord[:3, :], - file_name=outfile) with pytest.raises(AssertionError): # no. entries in mni matrix io.export_brain_net_viewer(adjacency_matrix=adj_matrix, mni_coord=mni_coord[:3, :], diff --git a/test/test_network_comparison.py b/test/test_network_comparison.py index 8458dbbd..2649b8b7 100644 --- a/test/test_network_comparison.py +++ b/test/test_network_comparison.py @@ -444,14 +444,17 @@ def test_compare_links_within(): network=res, data=data) for r in [res_indep, res_dep]: - assert (r.get_adjacency_matrix('diff_abs')[link_a[0], link_a[1]] == - r.get_adjacency_matrix('diff_abs')[link_b[0], link_b[1]]), ( + adj_mat_diff = r.get_adjacency_matrix('diff_abs') + adj_mat_comp = r.get_adjacency_matrix('comparison') + adj_mat_pval = r.get_adjacency_matrix('pvalue') + assert (adj_mat_diff._weight_matrix[link_a[0], link_a[1]] == + adj_mat_diff._weight_matrix[link_b[0], link_b[1]]), ( 'Absolute differences for link comparison not equal.') - assert (r.get_adjacency_matrix('comparison')[link_a[0], link_a[1]] == - r.get_adjacency_matrix('comparison')[link_b[0], link_b[1]]), ( + assert (adj_mat_comp._weight_matrix[link_a[0], link_a[1]] == + adj_mat_comp._weight_matrix[link_b[0], link_b[1]]), ( 'Comparison results for link comparison not equal.') - assert (r.get_adjacency_matrix('pvalue')[link_a[0], link_a[1]] == - r.get_adjacency_matrix('pvalue')[link_b[0], link_b[1]]), ( + assert (adj_mat_pval._weight_matrix[link_a[0], link_a[1]] == + adj_mat_pval._weight_matrix[link_b[0], link_b[1]]), ( 'P-value for link comparison not equal.') assert (r.targets_analysed == [link_a[1], link_b[1]]).all(), ( 'Analysed targets are not correct.') @@ -498,8 +501,10 @@ def test_tails(): network_set_b=np.array((res_2, res_3)), data_set_a=np.array((data, data)), data_set_b=np.array((data, data))) - print(c_within.get_adjacency_matrix('pvalue')) - print(c_between.get_adjacency_matrix('pvalue')) + adj_mat_within = c_within.get_adjacency_matrix('pvalue') + adj_mat_within.print_matrix() + adj_mat_between = c_between.get_adjacency_matrix('pvalue') + adj_mat_between.print_matrix() if __name__ == '__main__': diff --git a/test/test_results.py b/test/test_results.py index 26972bb5..373ed115 100644 --- a/test/test_results.py +++ b/test/test_results.py @@ -6,6 +6,7 @@ import itertools as it import copy as cp import numpy as np +from idtxl.results import AdjacencyMatrix from idtxl.multivariate_te import MultivariateTE from idtxl.bivariate_te import BivariateTE from idtxl.multivariate_mi import MultivariateMI @@ -128,13 +129,13 @@ def test_results_network_inference(): assert res.data_properties.normalised == normalisation, ( 'Incorrect value for data normalisation.') adj_matrix = res.get_adjacency_matrix('binary', fdr=False) - assert adj_matrix.shape[0] == n_nodes, ( + assert adj_matrix._edge_matrix.shape[0] == n_nodes, ( 'Incorrect number of rows in adjacency matrix.') - assert adj_matrix.shape[1] == n_nodes, ( + assert adj_matrix._edge_matrix.shape[1] == n_nodes, ( 'Incorrect number of columns in adjacency matrix.') - assert adj_matrix.shape[0] == n_nodes, ( + assert adj_matrix._edge_matrix.shape[0] == n_nodes, ( 'Incorrect number of rows in adjacency matrix.') - assert adj_matrix.shape[1] == n_nodes, ( + assert adj_matrix._edge_matrix.shape[1] == n_nodes, ( 'Incorrect number of columns in adjacency matrix.') @@ -249,10 +250,13 @@ def test_delay_reconstruction(): res_network.combine_results(nw.analyse_single_target( settings=settings, data=data, target=3)) adj_mat = res_network.get_adjacency_matrix('max_te_lag', fdr=False) - print(adj_mat) - assert adj_mat[0, 1] == delay_1, ('Estimate for delay 1 is not correct.') - assert adj_mat[0, 2] == delay_2, ('Estimate for delay 2 is not correct.') - assert adj_mat[0, 3] == delay_3, ('Estimate for delay 3 is not correct.') + adj_mat.print_matrix() + assert adj_mat._weight_matrix[0, 1] == delay_1, ( + 'Estimate for delay 1 is not correct.') + assert adj_mat._weight_matrix[0, 2] == delay_2, ( + 'Estimate for delay 2 is not correct.') + assert adj_mat._weight_matrix[0, 3] == delay_3, ( + 'Estimate for delay 3 is not correct.') for target in range(1, 4): est_mi = res_network._single_target[target].omnibus_te @@ -331,50 +335,103 @@ def test_results_network_comparison(): t = [1, 2] test = ['Within', 'Between'] for (i, res) in enumerate([res_within, res_between]): + # Get adjacency matrices. + adj_matrix_union = res.get_adjacency_matrix('union') + adj_matrix_comp = res.get_adjacency_matrix('comparison') + adj_matrix_diff = res.get_adjacency_matrix('diff_abs') + adj_matrix_pval = res.get_adjacency_matrix('pvalue') + + # Test all adjacency matrices for non-zero entries at compared edges + # and zero/False entries otherwise. + # Union network # TODO do we need the max_lag entry? - assert (res.get_adjacency_matrix('union')[s, t] == 1).all(), ( + assert adj_matrix_union._edge_matrix[s, t].all(), ( '{0}-test did not return correct union network links.'.format( test[i])) - no_diff = np.extract(np.invert(res.get_adjacency_matrix('comparison')), - res.get_adjacency_matrix('union')) - assert (no_diff == 0).all(), ( - '{0}-test did not return 0 in union network for no links.'.format( - test[i])) # Comparison - assert res.get_adjacency_matrix('comparison')[s, t].all(), ( + assert adj_matrix_comp._edge_matrix[s, t].all(), ( '{0}-test did not return correct comparison results.'.format( test[i])) - no_diff = np.extract(np.invert(res.get_adjacency_matrix('comparison')), - res.get_adjacency_matrix('comparison')) - assert (no_diff == 0).all(), ( + no_diff = np.extract(np.invert(adj_matrix_union._edge_matrix), + adj_matrix_comp._edge_matrix) + assert not no_diff.any(), ( '{0}-test did not return 0 comparison for non-sign. links.'.format( test[i])) # Abs. difference - assert (res.get_adjacency_matrix('diff_abs')[s, t] > 0).all(), ( + assert (adj_matrix_diff._edge_matrix[s, t]).all(), ( '{0}-test did not return correct absolute differences.'.format( test[i])) - no_diff = np.extract(np.invert(res.get_adjacency_matrix('comparison')), - res.get_adjacency_matrix('diff_abs')) - assert (no_diff == 0).all(), ( + no_diff = np.extract(np.invert(adj_matrix_union._edge_matrix), + adj_matrix_diff._edge_matrix) + assert not no_diff.any(), ( '{0}-test did not return 0 difference for non-sign. links.'.format( test[i])) # p-value p_max = 1 / comp_settings['n_perm_comp'] - assert (res.get_adjacency_matrix('pvalue')[s, t] == p_max).all(), ( + assert (adj_matrix_pval._weight_matrix[s, t] == p_max).all(), ( '{0}-test did not return correct p-value for sign. links.'.format( test[i])) - no_diff = np.extract(np.invert(res.get_adjacency_matrix('comparison')), - res.get_adjacency_matrix('pvalue')) - assert (no_diff == 1).all(), ( + no_diff = np.extract(np.invert(adj_matrix_union._edge_matrix), + adj_matrix_pval._edge_matrix) + assert not no_diff.any(), ( '{0}-test did not return p-vals of 1 for non-sign. links.'.format( test[i])) +def test_adjacency_matrix(): + # Test AdjacencyMatrix class + n_nodes = 5 + adj_mat = AdjacencyMatrix(n_nodes, int) + # Test adding single edge + edge = (0, 1) + weight = 1 + adj_mat.add_edge(edge[0], edge[1], weight) + assert adj_mat._weight_matrix[edge[0], edge[1]] == weight, ( + 'Weighted edge was not added.') + assert adj_mat._edge_matrix[edge[0], edge[1]], ( + 'Edge was not added.') + with pytest.raises(TypeError): + adj_mat.add_edge(0, 2, 3.5) + # Test adding edge list + i_list = [1, 1, 2, 4, 0] + j_list = [2, 3, 0, 3, 3] + weights = [2, 3, 4, 5, 6] + adj_mat.add_edge_list(i_list, j_list, weights) + for i in range(len(i_list)): + assert adj_mat._weight_matrix[i_list[i], j_list[i]] == weights[i], ( + 'Weighted edge was not added.') + assert adj_mat._edge_matrix[i_list[i], j_list[i]], ( + 'Edge was not added.') + adj_mat.add_edge_list(np.array(i_list), j_list, weights) + + # Test getting list of edges + adj_mat.get_edge_list() + + # Test comparison of entry types + a = AdjacencyMatrix(n_nodes, weight_type=bool) + a.add_edge(0, 1, True) + for t in [int, np.int, np.int32, np.int64]: + a = AdjacencyMatrix(n_nodes, weight_type=t) + a.add_edge(0, 1, int(1)) + a.add_edge(0, 1, np.int(1)) + a.add_edge(0, 1, np.int32(1)) + a.add_edge(0, 1, np.int64(1)) + for t in [float, np.float, np.float32, np.float64]: + a = AdjacencyMatrix(n_nodes, weight_type=t) + a.add_edge(0, 1, float(1)) + a.add_edge(0, 1, np.float(1)) + a.add_edge(0, 1, np.float32(1)) + a.add_edge(0, 1, np.float64(1)) + with pytest.raises(TypeError): + AdjacencyMatrix(n_nodes, weight_type=(2, 3)) + + if __name__ == '__main__': + test_adjacency_matrix() + test_console_output() test_results_network_inference() test_results_network_comparison() - test_console_output() test_pickle_results() test_delay_reconstruction() test_combine_results() From 55f672ae832b77fb9a0239b9ea175bb3c4114c89 Mon Sep 17 00:00:00 2001 From: Patricia Wollstadt Date: Sun, 2 Dec 2018 17:46:57 +0100 Subject: [PATCH 2/3] Remove redefinition of BROJA exception Fixes #22. --- idtxl/synergy_tartu.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/idtxl/synergy_tartu.py b/idtxl/synergy_tartu.py index ebed887a..2599d56d 100644 --- a/idtxl/synergy_tartu.py +++ b/idtxl/synergy_tartu.py @@ -47,9 +47,6 @@ def p_vidx(i): def q_vidx(i): return 3*i+2 -class BROJA_2PID_Exception(Exception): - pass - class Solve_w_ECOS: # (c) Abdullah Makkeh, Dirk Oliver Theis From 65aa8e34c0e440cb131fc47411af87b32485feb8 Mon Sep 17 00:00:00 2001 From: Patricia Wollstadt Date: Sun, 2 Dec 2018 18:00:50 +0100 Subject: [PATCH 3/3] Update documentation Update documentation. Add original reference for the example network used in the MuTE paper and make clear that the Sydney PID estimator estimates the BROJA measure. Add documentation for AdjacencyMatrix class. --- idtxl/data.py | 18 ++++++++++++++---- idtxl/estimators_pid.py | 13 +++++++++---- idtxl/results.py | 9 +-------- 3 files changed, 24 insertions(+), 16 deletions(-) diff --git a/idtxl/data.py b/idtxl/data.py index 6025ed86..a1934904 100755 --- a/idtxl/data.py +++ b/idtxl/data.py @@ -34,7 +34,7 @@ class Data(): >>> data_2.set_data(data_new, 's') Note: - Realisations are stored as attribute 'data'. This can only be set via + Realisations are stored as attribute 'data'. This can only be set via the 'set_data()' method. Args: @@ -828,9 +828,9 @@ def generate_mute_data(self, n_samples=1000, n_replications=10): Generate example data and overwrite the instance's current data. The network is used as an example the paper on the MuTE toolbox (Montalto, - PLOS ONE, 2014, eq. 14). The network consists of five autoregressive - (AR) processes with model orders 2 and the following (non-linear) - couplings: + PLOS ONE, 2014, eq. 14) and was orginially proposed by Baccala & + Sameshima (2001). The network consists of five autoregressive (AR) + processes with model orders 2 and the following (non-linear) couplings: 0 -> 1, u = 2 (non-linear) 0 -> 2, u = 3 @@ -838,6 +838,16 @@ def generate_mute_data(self, n_samples=1000, n_replications=10): 3 -> 4, u = 1 4 -> 3, u = 1 + References: + + - Montalto, A., Faes, L., & Marinazzo, D. (2014) MuTE: A MATLAB toolbox + to compare established and novel estimators of the multivariate + transfer entropy. PLoS ONE 9(10): e109462. + https://doi.org/10.1371/journal.pone.0109462 + - Baccala, L.A. & Sameshima, K. (2001). Partial directed coherence: a + new concept in neural structure determination. Biol Cybern 84: + 463–474. https://doi.org/10.1007/PL00007990 + Args: n_samples : int number of samples simulated for each process and replication diff --git a/idtxl/estimators_pid.py b/idtxl/estimators_pid.py index 2d0f80b4..b79d05ac 100755 --- a/idtxl/estimators_pid.py +++ b/idtxl/estimators_pid.py @@ -5,7 +5,6 @@ Bertschinger, N., Rauh, J., Olbrich, E., Jost, J., & Ay, N. (2014). Quantifying Unique Information. Entropy, 16(4), 2161–2183. http://doi.org/10.3390/e16042161 - """ import numpy as np from . import synergy_tartu @@ -17,9 +16,9 @@ class SydneyPID(Estimator): """Estimate partial information decomposition of discrete variables. - Fast implementation of the partial information decomposition (PID) - estimator for discrete data. The estimator does not require JAVA or GPU - modules to run. + Fast implementation of the BROJA partial information decomposition (PID) + estimator for discrete data (Bertschinger, 2014). The estimator does not + require JAVA or GPU modules to run. The estimator finds shared information, unique information and synergistic information between the two inputs s1 and s2 with respect to @@ -34,6 +33,12 @@ class SydneyPID(Estimator): CMI decreases; and an outer loop which decreases the size of the probability mass increment the virtualised swapping utilises. + References + + - Bertschinger, N., Rauh, J., Olbrich, E., Jost, J., & Ay, N. (2014). + Quantifying unique information. Entropy, 16(4), 2161–2183. + http://doi.org/10.3390/e16042161 + Args: settings : dict estimation parameters diff --git a/idtxl/results.py b/idtxl/results.py index a38481af..34f6a875 100644 --- a/idtxl/results.py +++ b/idtxl/results.py @@ -55,14 +55,7 @@ def __setstate__(self, state): class AdjacencyMatrix(): - """Adjacency matrix representing inferred networks. - - Attributes: TODO - is_edge : bool - matrix of edges in the network - weights : int | float | bool - weight for each edge - """ + """Adjacency matrix representing inferred networks.""" def __init__(self, n_nodes, weight_type): self._edge_matrix = np.zeros((n_nodes, n_nodes), dtype=bool) self._weight_matrix = np.zeros((n_nodes, n_nodes), dtype=weight_type)