Coherence graph
===============

Demonstration code for computing a coherence connectivity graph using CoherenceGraph().

Link to script: `feature_extraction/coherence_graph.py <https://gitlab.com/timhermans/nnsa_public/-/blob/main/examples/feature_extraction/coherence_graph.py>`_



.. code-block:: python

    import os
    
    import numpy as np
    import matplotlib.pyplot as plt
    
    from nnsa import CoherenceGraph, read_result_from_file, assert_equal, print_object_summary
    from nnsa.utils.plotting import scale_figsize, save_fig_as
    
    plt.close('all')
    fig_width = 18


Create some test signals. 



.. code-block:: python

    
    # Create a random multi-channel signal.
    signal = np.random.rand(8, 10000)
    fs = 250


Parameters. Descriptions of the parameters are documented in the default_parameters() code. 



.. code-block:: python

    
    # Print the default parameters of CoherenceGraph():
    print('\nDefault parameters:\n', CoherenceGraph().default_parameters())
    
    # Create an instance of the CoherenceGraph class with custom parameters, overruling some defaults:
    cg = CoherenceGraph(segmentation={'segment_length': 15})
    
    # See if the custom parameters were accepted:
    print('\nCustom parameters:\n', cg.parameters)


Now that we have initialized a CoherenceGraph object with certain parameters we can run the process method, which computes the coherence between channels which can be interpreted as a graph. 



.. code-block:: python

    # Code was ported from MATLAB and constsist of 6 steps:
    # 1) Segmentation
    # 2) Optional filtering
    # 3) Artefact exclusion
    # 4) PSD estimation
    # 5) Computation of coherence in specified frequency bands per segment
    # 6) Averaging the coherence in the specified frequency bands
    
    result = cg.process(data_matrix=signal, fs=fs)


The returned object is a CoherenceGraphResult object, which is a high-level interface for manipulating/visualizing/saving the result: 



.. code-block:: python

    print('\nSummary of result object:')
    print_object_summary(result)


CoherenceGraphResult attributes. The mean attributes are mean_coh and im_coh. Both are a dictionary, with frequency bands as keys and arrays as values. 



.. code-block:: python

    
    # The arrays in mean_coh have dimensions (n_channels, n_channels, n_segments) and the values are the mean coherence
    # in the frequency band corresponding to the key.
    print('result.mean_coh.keys(): {}'.format(result.mean_coh.keys()))
    print('result.mean_coh.values()[0].shape: {}'.format(list(result.mean_coh.values())[0].shape))
    
    # The arrays in im_coh have dimensions (n_channels, n_channels, n_segments) and the values are the largest (by
    # magnitude) imaginary part in the frequency band corresponding to the key.
    print('result.im_coh.keys(): {}'.format(result.im_coh.keys()))
    print('result.im_coh.values()[0].shape: {}'.format(list(result.im_coh.values())[0].shape))


Plot the results. 



.. code-block:: python

    # Plot graphs.
    fig = plt.figure(
        tight_layout=True,
        figsize=scale_figsize((4, 12), width=fig_width, unit='cm'),
    )
    result.plot(which='delta_1')
    

.. figure:: ../examples/feature_extraction/figs/coherence_graph_1.png


See https://nnsa.readthedocs.io/en/latest/nnsa.feature_extraction.html#nnsa.feature_extraction.connectivity.CoherenceGraphResult for a list of methods for computation of graph metrics. 



.. code-block:: python

    # Compute average path lengths.
    print(result.compute_average_path_lengths('mean_coh', graph_type='functional'))


We can save the result to an hdf5 file. 



.. code-block:: python

    
    # Include the extension of the filename.
    filename = 'temp_result.hdf5'
    result.save_to_file(filename)
    
    # We can reload the result back to a CoherenceGraphResult object:
    result_loaded = read_result_from_file(filename)
    
    # Remove temporary file.
    os.remove(filename)
    
    # Verify that the loaded object is equal to the original object:
    assert_equal(actual=result_loaded, desired=result)  # Will raise an AssertionError if not equal.
    print('Loaded result object is equal to original object.')