Power analysis
==============

Demonstration code for performing a power analysis using PowerAnalysis().

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



.. code-block:: python

    import os
    import matplotlib.pyplot as plt
    import numpy as np
    
    from nnsa import PowerAnalysis, EdfReader, print_object_summary, SUPPORTED_RESULT_FILE_TYPES, read_result_from_file, \
        assert_equal
    
    plt.close('all')


Parameters. 



.. code-block:: python

    # Print the default parameters of PowerAnalysis():
    print(PowerAnalysis().default_parameters())
    
    # Descriptions of the parameters are documented in the default_parameters() code.
    
    # Create an instance of the PowerAnalysis class with custom parameters, overruling some defaults:
    power_analysis = PowerAnalysis(artefact_criteria={'max_diff': 200},
                                   psd={'window_length': 2.5},
                                   segmentation={'segment_length': 15, 'overlap': 1})
    
    # See if the custom parameters were accepted:
    print('\nCustom parameters:')
    print(power_analysis.parameters)


Main method: power_analysis. 



.. code-block:: python

    # Now that we have initialized a PowerAnalysis object with certain parameters, we can run the power_analysis method,
    # which performs the actual power analysis:
    # Create 10 random signals to simulate 10 channel EEG data (amplitude 50).
    fs = 250
    data_matrix = np.random.rand(10, 100000)*100 - 50
    result = power_analysis.power_analysis(data_matrix, fs=fs, unit='uV', verbose=1)  # Not that some args are optional.
    
    # The returned object is a PowerAnalysisResult object, which is a high-level interface for
    # manipulating/visualizing/saving the result:
    print('\nSummary of result object:')
    print_object_summary(result)


PowerAnalysisResult attributes. 



.. code-block:: python

    # The main attribute is 'power', which contains the frequency power in certain bands:
    print('result.power: {}'.format(result.power))
    
    # power is a dictionary with the frequency bands as keys and its values are 2D arrays corresponding to
    # (channels, segments).
    freq_bands = list(result.power.keys())
    print('Frequency bands: {}'.format(freq_bands))
    print('power[\'{}\'].shape: {}'.format(freq_bands[0], result.power[freq_bands[0]].shape))
    
    # The time (in seconds) corresponding to the segments in the data in power:
    print('segment_times: {}'.format(result.segment_times))


PowerAnalysisResult methods. 



.. code-block:: python

    # Print a summary of the power analysis:
    result.summary()
    
    # Plot the power for all channels as function of segments.
    result.plot()


Save the result. 



.. code-block:: python

    # We can save the result to any supported file. To list the supported file types/extensions:
    print('Supported result file types: {}'.format(SUPPORTED_RESULT_FILE_TYPES))
    
    # E.g. save as csv (the file type to save is automatically detected from the extension of the filename).
    filename = 'temp_result.csv'
    result.save_to_file(filename)
    
    # We can reload the result back to a PowerAnalysisResult 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.')


PowerAnalysis of EegDataset. 



.. code-block:: python

    # Specify a file and open a reader to read the data (e.g. EdfReader to read an EDF(+) file).
    filepath = 'C:/data_temp/test.edf'
    
    # Read EEG channels in an EegDataset object.
    with EdfReader(filepath) as r:
        ds = r.read_eeg_dataset(dtype=np.float32)  # We can specify the dtype, lower precision, means less memory usage.
    
    # We can do the power analysis for multichannel EEG via the wrapper method 'power_analysis' of the EegDataset class.
    # Specify the parameters to the PowerAnalysis class as kwargs. Preprocessing and power analysis can be done in one line:
    result_ds = ds.reference('Cz').resample(fs_new=128).filter_saved_filter('bandpassfir_a').power_analysis(segmentation={'overlap': 0})
    
    # A data info string is appended to power_result object, which summarizes the source and preprocessing of the data.
    print('result_ds.data_info: {}'.format(result_ds.data_info))
    
    # Note that the dtype of the result is the same as the dtype of the data (e.g. np.float32).
    print('dtype result: {}'.format(result_ds.power[list(result_ds.power.keys())[0]].dtype))
    
    # Print a summary.
    result_ds.summary()
    
    # Plot the power for all channels as function of segments.
    plt.figure()
    result_ds.plot()