binned_spiking_activity_human_epilepsy/python_code/analysis_PlosCB/create_figure_2.py

import os
import pandas as pd
import numpy as np
import matplotlib.pylab as plt
import seaborn as sns
import mr_utilities
from scipy import stats

mycol = ['darkgreen', 'mediumaquamarine',  'silver', 'darkblue','steelblue']


def cm2inch(value):
    return value/2.54


def mt_examples(patient, save_path, fit_method, kmin, kmax,
                windowsize, windowstep, dt, excludeInconsistencies=True):

    recordings = mr_utilities.list_preictrecordings(patient)[0]
    fig, ax = plt.subplots(2, 1, figsize=(cm2inch(13.5), cm2inch(5)))
    binnings = ['right', 'left']
    quantity = 'm'
    t_onset = 600
    for irec, rec in enumerate(recordings):
        for ibin, binning in enumerate(binnings):
            mt_file = '{}{}/mt_results_{}_{}_kmin{}_kmax{}_winsize{}_' \
                      'winstep{}_dt{}.pkl'.format(save_path, rec, binning,
                       fit_method, kmin, kmax, windowsize, windowstep, dt)

            if not os.path.isfile(mt_file):
                print('{} not existing'.format(mt_file))
                continue

            # -----------------------------------------------------------------#
            # read m(t) data
            # -----------------------------------------------------------------#
            mt_frame = pd.read_pickle(mt_file)
            if excludeInconsistencies:
                rejected_inconsistencies = ['H_nonzero', 'H_linear']
                for incon in rejected_inconsistencies:
                    if not mt_frame.empty:
                        mt_frame = mt_frame[[incon not in mt_frame[
                        'inconsistencies'].tolist()[i] for i in range(
                        len(mt_frame['inconsistencies'].tolist()))]]
            if mt_frame.empty:
                continue

            mt_frame['logtau'] = mt_frame.apply(lambda row: np.log(row['tau']),
                                                axis=1)
            mt_frame['epsilon'] = mt_frame.apply(lambda row: 1-row['m'],
                                                 axis=1)
            kmin = mt_frame.kmin.unique()[0]
            kmax = mt_frame.kmax.unique()[0]
            fit_method = mt_frame.fit_method.unique()[0]
            winsize = mt_frame.windowsize.unique()[0]
            t0_list = mt_frame.t0.tolist()
            qt = mt_frame[quantity].tolist()
            t0_list = [t0_ms * 0.001 for t0_ms in t0_list]
            half_win_sec = int((winsize / 2) / 1000 * dt)
            t_middle = [t0 + half_win_sec for t0 in t0_list]

            # -----------------------------------------------------------------#
            # Plotting
            # -----------------------------------------------------------------#
            ax[ibin].plot(t_middle, qt, '-', linewidth=1,
                          label='Rec {}'.format(irec+1), color=mycol[irec])
            ax[ibin].axhline(y=1, color='black', linestyle='--', alpha=0.3)
            ax[ibin].set_ylabel(r"$\hat{m}$")
            ax[ibin].set_xlim(
                (t0_list[0] - 10, t0_list[-1] + 2 * half_win_sec + 10))
            binning_focus = mr_utilities.binning_focus(binning, int(patient))
            if binning_focus == 'focal':
                focus_label = 'ipsi'
            else:
                focus_label = 'contra'
            ax[ibin].set_title("{}".format(focus_label), loc='left')

    ax[0].legend(loc='center left', bbox_to_anchor=(1.04, -0.2),
                 fontsize='small')
    ax[0].xaxis.grid(False)
    ax[1].xaxis.grid(False)
    ax[1].set_xlabel("time (s)")
    ax[0].annotate("", xy=(t_onset, 1), xytext=(t_onset, 1.1),
                  arrowprops=dict(arrowstyle="->", color='black'))
    ax[1].annotate("", xy=(t_onset, 1), xytext=(t_onset, 1.04),
                      arrowprops=dict(arrowstyle="->", color='black'))
    ax[0].annotate("seizure\n onset", xy=(t_onset, 1),
                      xytext=(t_onset, 1.03))
    sns.despine()
    plt.subplots_adjust(bottom=None, right=None,
                        wspace=0.3, hspace=1.3, left=None, top=None)
    plt.savefig(
        '{}/Fig2_mt_{}_{}_winsize{}_winstep{}_kmin{}_kmax{}.pdf'.format(
            save_path, patient, fit_method,
            windowsize, windowstep, kmin, kmax), bbox_inches='tight')
    plt.close()


def boxplot_mx(result_path, x, ksteps, fit_method,
                     excludeInconsistencies=True, Engel=False):
    # -----------------------------------------------------------------#
    # Read and filter data
    # -----------------------------------------------------------------#
    mx_results = '{}mt_{}_results_allrecords_kmin{}_kmax{}_{}.pkl'.format(
        result_path, x, ksteps[0], ksteps[1], fit_method)
    mx_allresults = pd.read_pickle(mx_results)

    if excludeInconsistencies:
        rejected_inconsistencies = ['H_nonzero', 'H_linear']
        for incon in rejected_inconsistencies:
            mx_allresults = mx_allresults[[incon not in mx_allresults[
                'inconsistencies'].tolist()[i] for i in range(
                    len(mx_allresults['inconsistencies'].tolist()))]]

    x = mx_allresults['partno'].max()
    mx_allresults['logtau'] = mx_allresults.apply(lambda row: np.log(row['tau']),
                                                  axis=1)
    mx_allresults['logepsilon'] = mx_allresults.apply(lambda row: abs(np.log(
        1-row['m'])), axis=1)
    mx_allresults['epsilon'] = mx_allresults.apply(lambda row:
        1-row['m'], axis=1)
    measure = 'epsilon'

    patients1a = [1, 6, 8, 9, 10, 11, 16, 17, 18, 19]
    if Engel:
        mx_allresults = mx_allresults[(mx_allresults['patient'].isin(patients1a))]
    # -----------------------------------------------------------------#
    # write estimates to lists
    # -----------------------------------------------------------------#
    focal_diff_list = []
    contralateral_diff_list = []
    all_m_first_focal = []
    all_m_last_focal = []
    all_m_first_contra = []
    all_m_last_contra = []

    for ip, patient in enumerate(mx_allresults['patient'].unique()):
        mx_patient = mx_allresults[mx_allresults['patient'] == patient]
        recordings = mx_patient.recording.unique()
        focal_diff_patient = []
        contralateral_diff_patient = []
        m_contra_first_patient = []
        m_contra_last_patient = []
        m_focal_first_patient = []
        m_focal_last_patient = []

        for recording in recordings:

            mx_recording = mx_patient[mx_patient['recording'] == recording]
            mx_focal = mx_recording[mx_recording['binning_focus'] == 'focal']
            mx_contralateral = mx_recording[
                mx_recording['binning_focus'] == 'contra-lateral']

            if not mx_focal[mx_focal['partno'] == 1].empty and not mx_focal[
                mx_focal['partno'] == x].empty:
                m_focal_first = mx_focal[mx_focal['partno'] == 1][
                    measure].item()
                m_focal_last = mx_focal[mx_focal['partno'] == x][measure].item()
                focal_diff = m_focal_last - m_focal_first
                focal_diff_log = np.log(m_focal_last) - np.log(m_focal_first)
                focal_diff_patient.append(focal_diff_log)
                m_focal_first_patient.append(m_focal_first)
                m_focal_last_patient.append(m_focal_last)

            if not mx_contralateral[mx_contralateral['partno'] == 1].empty \
                and not mx_contralateral[mx_contralateral['partno'] == x].empty:
                m_contralateral_first = mx_contralateral[mx_contralateral['partno'] == 1][measure].item()
                m_contralateral_last = mx_contralateral[mx_contralateral['partno'] == x][measure].item()
                contralateral_diff = m_contralateral_last - m_contralateral_first
                contralateral_diff_log = np.log(m_contralateral_last) - \
                                     np.log(m_contralateral_first)
                contralateral_diff_patient.append(contralateral_diff_log)
                m_contra_first_patient.append(m_contralateral_first)
                m_contra_last_patient.append(m_contralateral_last)

        if not len(focal_diff_patient) == 0:
            focal_diff_list.append(focal_diff_patient)
            all_m_first_focal.append(m_focal_first_patient)
            all_m_last_focal.append(m_focal_last_patient)
        if not len(contralateral_diff_patient) == 0:
            contralateral_diff_list.append(contralateral_diff_patient)
            all_m_first_contra.append(m_contra_first_patient)
            all_m_last_contra.append(m_contra_last_patient)

    # -----------------------------------------------------------------#
    # Permutation test
    # -----------------------------------------------------------------#
    focal_diff_list = [item for sublist in focal_diff_list for item in sublist]
    contralateral_diff_list = [item for sublist in contralateral_diff_list for item in sublist]
    all_last_contra = [item for sublist in all_m_last_contra for item in sublist]
    all_first_contra = [item for sublist in all_m_first_contra for item in
                         sublist]
    all_last_focal = [item for sublist in all_m_last_focal for item in
                         sublist]
    all_first_focal = [item for sublist in all_m_first_focal for item in
                         sublist]

    # -----------------------------------------------------------------#
    # Wilcoxon paired test
    # -----------------------------------------------------------------#
    stat, p_contra_w = stats.wilcoxon(contralateral_diff_list)
    stat, p_focal_w = stats.wilcoxon(focal_diff_list)

    # -----------------------------------------------------------------#
    # Create dataframe for plotting
    # -----------------------------------------------------------------#
    quantity = measure
    data1 = {
        quantity: all_first_focal,
        'rec_type': ['10 - 5 min' for _ in range(len(all_first_focal))],
        'binning_focus': ['SOZ' for _ in range(len(all_first_focal))]
    }

    df = pd.DataFrame(data1)

    data2 = pd.DataFrame({
        quantity: all_last_focal,
        'rec_type': ['5 - 0 min' for _ in range(len(all_last_focal))],
        'binning_focus': ['SOZ' for _ in range(len(all_last_focal))]
    })
    df = df.append(data2, ignore_index=True)

    data3 = pd.DataFrame({
        quantity: all_first_contra,
        'rec_type': ['10 - 5 min' for _ in range(len(all_first_contra))],
        'binning_focus': ['nSOZ' for _ in range(len(all_first_contra))]
    })
    df = df.append(data3, ignore_index=True)

    data4 = pd.DataFrame({
        quantity: all_last_contra,
        'rec_type': ['5 - 0 min' for _ in range(len(all_last_contra))],
        'binning_focus': ['nSOZ' for _ in range(len(all_last_contra))]
    })
    df = df.append(data4, ignore_index=True)

    # -----------------------------------------------------------------#
    # Plotting
    # -----------------------------------------------------------------#
    fig, ax = plt.subplots(1, 1, figsize=(cm2inch(5.5), cm2inch(4.5)))
    g = sns.boxplot(ax=ax, x='binning_focus', y=quantity, hue='rec_type',
                  data=df, hue_order = ['10 - 5 min', '5 - 0 min'], order=[
                  'SOZ', 'nSOZ'], dodge=True, color='slategray', palette=[
                  'lightgray', 'slategray'], linewidth=1.1, width=0.6)

    g.set_yscale("log")
    plt.grid(axis='x', b=False)
    plt.rcParams['legend.title_fontsize'] = 9
    plt.legend(loc="upper left", ncol=1,
               fontsize='x-small', bbox_to_anchor=(0, 1.7), title='time to '
               'seizure onset')
    n_f_pre = len(all_first_focal)
    n_f_inter = len(all_last_focal)
    n_c_pre = len(all_first_contra)
    n_c_inter = len(all_last_contra)
    ax.set_xticklabels(['ipsi\n({}, {})'.format(n_f_pre, n_f_inter),
                        'contra\n({}, {})'.format(n_c_pre, n_c_inter)])

    # -----------------------------------------------------------------#
    # quantity labels
    # -----------------------------------------------------------------#
    quantity_label = r'$\hat{m}$'
    ymax = 1e-3
    ymin = 1
    ax.set_ylim(ymin, ymax)
    T = [1, 1e-1, 1e-2, 1e-3]
    ax.set_yticks(T)
    ax.set_yticklabels([0, 0.9, 0.99, 0.999])

    # -----------------------------------------------------------------#
    # annotate p-values
    # -----------------------------------------------------------------#
    y_max = 2.4*ymax
    y_min = ymin
    ax.annotate("", xy=(-0.2, y_max), xycoords='data',
                xytext=(0.2, y_max), textcoords='data',
                arrowprops=dict(arrowstyle="-", ec='black',
                                connectionstyle="bar,fraction=0.2"))

    ax.text(-0.25, y_max*0.6, 'p={:.3f}'.format(p_focal_w), fontsize=8)
    ax.annotate("", xy=(0.8, y_max), xycoords='data',
                xytext=(1.2, y_max), textcoords='data',
                arrowprops=dict(arrowstyle="-", ec='black',
                                connectionstyle="bar,fraction=0.2"))

    ax.text(0.75, y_max*0.6, 'p={:.3f}'.format(p_contra_w), fontsize=8)
    ax.set_ylim(ymin, ymax)
    plt.ylabel(quantity_label)
    plt.xlabel('')
    sns.despine()

    # -----------------------------------------------------------------#
    # Saving
    # -----------------------------------------------------------------#
    kmin = mx_allresults['kmin'].unique()[0]
    kmax = mx_allresults['kmax'].unique()[0]
    plotfile_mx = '{}Fig2_mx{}_{}_kmin{}_kmax{}_Engel{}.pdf'.format(
                   result_path, x, measure, kmin, kmax, Engel)
    plt.savefig(plotfile_mx, bbox_inches='tight')


def pointplot_mx_regions_Engel(region, result_path, x, ksteps,
                               fit_method, excludeInconsistencies=True):
    # -----------------------------------------------------------------#
    # Read and filter data
    # -----------------------------------------------------------------#
    mx_results = '{}mt_{}_results_allrecords_kmin{}_kmax{}_{}.pkl'.format(
        result_path, x, ksteps[0], ksteps[1], fit_method)
    mx_allresults = pd.read_pickle(mx_results)

    if excludeInconsistencies:
        rejected_inconsistencies = ['H_nonzero', 'H_linear']
        for incon in rejected_inconsistencies:
            mx_allresults = mx_allresults[[incon not in mx_allresults[
                'inconsistencies'].tolist()[i] for i in range(
                    len(mx_allresults['inconsistencies'].tolist()))]]

    x = mx_allresults['partno'].max()
    mx_allresults['logtau'] = mx_allresults.apply(lambda row: np.log(row['tau']),
                                                  axis=1)
    mx_allresults['logepsilon'] = mx_allresults.apply(lambda row: abs(np.log(
        1-row['m'])), axis=1)
    mx_allresults['epsilon'] = mx_allresults.apply(lambda row:
        1-row['m'], axis=1)
    measure = 'epsilon'

    def give_region(row):
        return row['binning'][1:]

    mx_allresults['region'] = mx_allresults.apply(lambda row: give_region(row),
                                                  axis=1)
    if not region == 'all':
        mx_allresults = mx_allresults[mx_allresults['region'] == region]
    # -----------------------------------------------------------------#
    # write estimates to lists
    # -----------------------------------------------------------------#
    fig, ax = plt.subplots(1, 1, figsize=(cm2inch(10*1.2), cm2inch(4*1.2)))
    df = None
    patients1a = [1, 6, 8, 9, 10, 11, 16, 17, 18, 19]
    mx_allresults = mx_allresults[(mx_allresults['patient'].isin(patients1a))]
    regions = ['H', 'A', 'PHC', 'EC']
    labels = []
    pvals = dict()
    for i, reg in enumerate(regions):
        mx_region = mx_allresults[mx_allresults['region'] == reg]
        recordings = mx_region.recording.unique()
        focal_diff_list = []
        contralateral_diff_list = []
        all_first_focal = []
        all_last_focal = []
        all_first_contra = []
        all_last_contra = []
        recID = []
        for recording in recordings:
            mx_recording = mx_region[mx_region['recording'] == recording]
            mx_focal = mx_recording[mx_recording['binning_focus'] == 'focal']
            mx_contralateral = mx_recording[
                mx_recording['binning_focus'] == 'contra-lateral']

            if not mx_focal[mx_focal['partno'] == 1].empty and not mx_focal[
                mx_focal['partno'] == x].empty:
                m_focal_first = mx_focal[mx_focal['partno'] == 1][
                    measure].item()
                m_focal_last = mx_focal[mx_focal['partno'] == x][measure].item()
                focal_diff = m_focal_last - m_focal_first
                all_first_focal.append(m_focal_first)
                all_last_focal.append(m_focal_last)
                focal_diff_list.append(focal_diff)
                recID.append(recording)

            if not mx_contralateral[mx_contralateral['partno'] == 1].empty \
                    and not \
                    mx_contralateral[
                mx_contralateral['partno'] == x].empty:
                m_contralateral_first = mx_contralateral[mx_contralateral['partno'] == 1][measure].item()
                m_contralateral_last = mx_contralateral[mx_contralateral['partno'] == x][measure].item()
                contralateral_diff = m_contralateral_last - m_contralateral_first

                all_first_contra.append(m_contralateral_first)
                all_last_contra.append(m_contralateral_last)
                contralateral_diff_list.append(contralateral_diff)

        # -----------------------------------------------------------------#
        # Wilcoxon paired test
        # -----------------------------------------------------------------#
        stat, p_contra_w = stats.wilcoxon(contralateral_diff_list)
        stat, p_focal_w = stats.wilcoxon(focal_diff_list)
        pvals[reg] = p_focal_w

        # -----------------------------------------------------------------#
        # Create dataframe for plotting
        # -----------------------------------------------------------------#
        quantity = measure
        if df is None:
            data1 = {
                'region': reg,
                quantity: all_first_focal,
                'recording': recID,
                'rec_type': ['10 - 5 min' for _ in range(len(all_first_focal))],
                'binning_focus': ['SOZ' for _ in range(len(all_first_focal))]
            }

            df = pd.DataFrame(data1)
        else:
            data1 = pd.DataFrame({
                'region': reg,
                quantity: all_first_focal,
                'recording': recID,
                'rec_type': ['10 - 5 min' for _ in range(len(all_last_focal))],
                'binning_focus': ['SOZ' for _ in range(len(all_last_focal))]
            })
            df = df.append(data1, ignore_index=True)

        data2 = pd.DataFrame({
            'region': reg,
            quantity: all_last_focal,
            'recording': recID,
            'rec_type': ['5 - 0 min' for _ in range(len(all_last_focal))],
            'binning_focus': ['SOZ' for _ in range(len(all_last_focal))]
        })
        df = df.append(data2, ignore_index=True)

        n_f_first = len(all_first_focal)
        n_f_last = len(all_last_focal)
        labels.append('{}\n({}, {})'.format(reg, n_f_first, n_f_last))

        x0 = (-0.2+i) * np.ones(len(all_first_focal))
        x1 = (0.2+i) * np.ones(len(all_last_focal))
        plt.plot([x0, x1], [all_first_focal, all_last_focal],
                 color='grey', linewidth=0.5, linestyle='--')

    # -----------------------------------------------------------------#
    # Plotting
    # -----------------------------------------------------------------#
    df = df[df['binning_focus']=='SOZ']
    g = sns.stripplot(ax=ax, x='region', y=quantity, hue='rec_type',
                      data=df, hue_order=['10 - 5 min', '5 - 0 min'],
                      dodge=True, order=regions, jitter=False,
                      color='slategray', palette=['slategray', 'darkblue'],
                      size=2.5, alpha=0.7)

    ax.set_yscale("log")
    ax.grid(axis='x', b=False)
    plt.rcParams['legend.title_fontsize'] = 9
    plt.legend(loc='center left', bbox_to_anchor=(1, 0.5),
               fontsize='x-small', title='time to seizure onset')
    ax.set_xticklabels(labels)
    # -----------------------------------------------------------------#
    # quantity labels
    # -----------------------------------------------------------------#
    quantity_label = r'$\hat{m}$'
    ymax = 1e-3
    ymin = 1
    ax.set_ylim(ymin, ymax)
    T = [1, 1e-1, 1e-2, 1e-3]
    ax.set_yticks(T)
    ax.set_yticklabels([0, 0.9, 0.99, 0.999])

    # -----------------------------------------------------------------#
    # Annotate p-values of pairwise comparison
    # -----------------------------------------------------------------#
    ax.set_ylim(ymin, ymax)
    plt.ylabel(quantity_label)
    plt.xlabel('brain area')
    sns.despine()

    # -----------------------------------------------------------------#
    # Saving
    # -----------------------------------------------------------------#
    kmin = mx_allresults['kmin'].unique()[0]
    kmax = mx_allresults['kmax'].unique()[0]
    plotfile_mx = '{}Fig2_pointmxallregions{}_{}_kmin{}_kmax{}_Engel1a.pdf'.format(
                   result_path, x, measure, kmin, kmax, region)
    plt.savefig(plotfile_mx, bbox_inches='tight')


def boxplot_mx_regions(region, result_path, x, ksteps, fit_method,
                       excludeInconsistencies=True):
    # -----------------------------------------------------------------#
    # Read and filter data
    # -----------------------------------------------------------------#
    mx_results = '{}mt_{}_results_allrecords_kmin{}_kmax{}_{}.pkl'.format(
        result_path, x, ksteps[0], ksteps[1], fit_method)
    mx_allresults = pd.read_pickle(mx_results)

    if excludeInconsistencies:
        rejected_inconsistencies = ['H_nonzero', 'H_linear']
        for incon in rejected_inconsistencies:
            mx_allresults = mx_allresults[[incon not in mx_allresults[
                'inconsistencies'].tolist()[i] for i in range(
                len(mx_allresults['inconsistencies'].tolist()))]]

    x = mx_allresults['partno'].max()
    mx_allresults['logtau'] = mx_allresults.apply(lambda row: np.log(row['tau']),
                                                  axis=1)
    mx_allresults['logepsilon'] = mx_allresults.apply(lambda row: abs(np.log(
        1-row['m'])), axis=1)
    mx_allresults['epsilon'] = mx_allresults.apply(lambda row:
        1-row['m'], axis=1)
    measure = 'epsilon'

    def give_region(row):
        return row['binning'][1:]

    mx_allresults['region'] = mx_allresults.apply(lambda row: give_region(row),
                                                  axis=1)

    if not region == 'all':
        mx_allresults = mx_allresults[mx_allresults['region'] == region]

    # -----------------------------------------------------------------#
    # write estimates to lists
    # -----------------------------------------------------------------#
    df = None
    regions = ['H', 'A', 'PHC', 'EC']
    labels = []
    pvals = dict()
    for reg in regions:
        mx_region = mx_allresults[mx_allresults['region'] == reg]
        recordings = mx_region.recording.unique()
        focal_diff_list = []
        contralateral_diff_list = []
        all_first_focal = []
        all_last_focal = []
        all_first_contra = []
        all_last_contra = []
        recID = []
        for recording in recordings:
            mx_recording = mx_region[mx_region['recording'] == recording]
            mx_focal = mx_recording[mx_recording['binning_focus'] == 'focal']
            mx_contralateral = mx_recording[
                mx_recording['binning_focus'] == 'contra-lateral']

            if not mx_focal[mx_focal['partno'] == 1].empty and not mx_focal[
                mx_focal['partno'] == x].empty:
                m_focal_first = mx_focal[mx_focal['partno'] == 1][
                    measure].item()
                m_focal_last = mx_focal[mx_focal['partno'] == x][measure].item()
                focal_diff = m_focal_last - m_focal_first
                focal_diff_log = np.log(m_focal_last) - np.log(m_focal_first)
                all_first_focal.append(m_focal_first)
                all_last_focal.append(m_focal_last)
                focal_diff_list.append(focal_diff_log)
                recID.append(recording)

            if not mx_contralateral[mx_contralateral['partno'] == 1].empty \
                    and not mx_contralateral[mx_contralateral['partno'] == x].empty:
                m_contralateral_first = mx_contralateral[mx_contralateral['partno'] == 1][measure].item()
                m_contralateral_last = mx_contralateral[mx_contralateral['partno'] == x][measure].item()
                contralateral_diff = m_contralateral_last - m_contralateral_first
                contralateral_diff_log = np.log(m_contralateral_last) - \
                                     np.log(m_contralateral_first)
                all_first_contra.append(m_contralateral_first)
                all_last_contra.append(m_contralateral_last)
                contralateral_diff_list.append(contralateral_diff_log)

        # -----------------------------------------------------------------#
        # Wilcoxon paired test
        # -----------------------------------------------------------------#
        stat, p_contra_w = stats.wilcoxon(contralateral_diff_list)
        stat, p_focal_w = stats.wilcoxon(focal_diff_list)
        pvals[reg] = p_focal_w

        # -----------------------------------------------------------------#
        # Create dataframe for plotting
        # -----------------------------------------------------------------#
        quantity = measure
        if df is None:
            data1 = {
                'region': reg,
                quantity: all_first_focal,
                'recording': recID,
                'rec_type': ['10 - 5 min' for _ in range(len(all_first_focal))],
                'binning_focus': ['SOZ' for _ in range(len(all_first_focal))]
            }
            df = pd.DataFrame(data1)
        else:
            data1 = pd.DataFrame({
                'region': reg,
                quantity: all_first_focal,
                'recording': recID,
                'rec_type': ['10 - 5 min' for _ in range(len(all_last_focal))],
                'binning_focus': ['SOZ' for _ in range(len(all_last_focal))]
            })
            df = df.append(data1, ignore_index=True)

        data2 = pd.DataFrame({
            'region': reg,
            quantity: all_last_focal,
            'recording': recID,
            'rec_type': ['5 - 0 min' for _ in range(len(all_last_focal))],
            'binning_focus': ['SOZ' for _ in range(len(all_last_focal))]
        })
        df = df.append(data2, ignore_index=True)

        data3 = pd.DataFrame({
            'region': reg,
            quantity: all_first_contra,
            'recording': ['-' for _ in range(len(all_first_contra))],
            'rec_type': ['10 - 5 min' for _ in range(len(all_first_contra))],
            'binning_focus': ['nSOZ' for _ in range(len(all_first_contra))]
        })
        df = df.append(data3, ignore_index=True)

        data4 = pd.DataFrame({
            'region': reg,
            quantity: all_last_contra,
            'recording': ['-' for _ in range(len(all_first_contra))],
            'rec_type': ['5 - 0 min' for _ in range(len(all_last_contra))],
            'binning_focus': ['nSOZ' for _ in range(len(all_last_contra))]
        })
        df = df.append(data4, ignore_index=True)
        n_f_first = len(all_first_focal)
        n_f_last = len(all_last_focal)
        labels.append('{}\n({}, {})'.format(reg, n_f_first, n_f_last))

    # -----------------------------------------------------------------#
    # Plotting
    # -----------------------------------------------------------------#
    fig, ax = plt.subplots(1, 1, figsize=(cm2inch(10*1.2), cm2inch(4*1.2)))
    df = df[df['binning_focus'] == 'SOZ']
    g = sns.boxplot(ax=ax, x='region', y=quantity, hue='rec_type',
                    data=df, hue_order = ['10 - 5 min', '5 - 0 min'],
                    order=regions, dodge=True, color='slategray', palette=[
                    'lightgray','slategray'], linewidth=1.1, width=0.8)

    g.set_yscale("log")
    ax.grid(axis='x', b=False)
    plt.rcParams['legend.title_fontsize'] = 9
    plt.legend(loc='center left', bbox_to_anchor=(1, 0.5),
               fontsize='x-small', title='time to seizure onset')
    ax.set_xticklabels(labels)

    # -----------------------------------------------------------------#
    # quantity labels
    # -----------------------------------------------------------------#
    quantity_label = r'$\hat{m}$'
    ymax = 1e-3
    ymin = 1
    ax.set_ylim(ymin, ymax)
    T = [1, 1e-1, 1e-2, 1e-3]
    ax.set_yticks(T)
    ax.set_yticklabels([0, 0.9, 0.99, 0.999])

    # -----------------------------------------------------------------#
    # Annotate p-values of pairwise comparison
    # -----------------------------------------------------------------#
    y_max = 1.3 * ymax
    y_min = ymin
    for ir, region in enumerate(regions):
        if not np.isnan(pvals[region]):
            ax.annotate("", xy=(-0.2 + ir * 1.0, 2.0 * y_max), xycoords='data',
                        xytext=(0.2 + ir * 1.0, 2.0 * y_max),
                        textcoords='data',
                        arrowprops=dict(arrowstyle="-", ec='black',
                                        connectionstyle="bar,fraction=0.3"))
            ax.text(-0.35 + ir * 1.0, 0.9*y_max * 1. - abs(y_max - y_min) *
                    0.0,
                    'p' + '={:.3f}'.format(pvals[region]), fontsize=8.5)

    ax.set_ylim(ymin, ymax)
    plt.ylabel(quantity_label)
    plt.xlabel('brain area')
    sns.despine()
    # -----------------------------------------------------------------#
    # Saving
    # -----------------------------------------------------------------#
    kmin = mx_allresults['kmin'].unique()[0]
    kmax = mx_allresults['kmax'].unique()[0]
    plotfile_mx = '{}Fig2_mxallregions{}_{}_kmin{}_kmax{}.pdf'.format(
        result_path, x, measure, kmin, kmax, region)
    plt.savefig(plotfile_mx, bbox_inches='tight')


def binning_label(patient, soz):
    focifile = '../Data/patients.txt'
    f = open(focifile, 'r')
    foci = [line.rstrip('\n').split('\t') for line in f.readlines()]

    focus = 'NA'
    for i, idf in enumerate(foci[1:]):
        if int(idf[0]) == int(patient):
            focus = idf[1]

    if soz == "SOZ":
        if focus == "L":
            binning = "left"
        elif focus == "R":
            binning = "right"
        else:
            raise Exception("No clear focus found.")

    elif soz == "nSOZ":
        if focus == "L":
            binning = "right"
        elif focus == "R":
            binning = "left"
        else:
            raise Exception("No clear focus found.")
    return binning


def get_mt(patient, recording, soz, save_path, fit_method, kmin, kmax,
            windowsize, windowstep, dt, excludeInconsistencies=True):

    binning = binning_label(patient, soz)
    mt_file = '{}{}/mt_results_{}_{}_kmin{}_kmax{}_winsize{}_winstep{}_dt{}.pkl'.format(
        save_path, recording, binning, fit_method, kmin, kmax,
        windowsize, windowstep, dt)

    if not os.path.isfile(mt_file):
        raise Exception('{} not existing'.format(mt_file))

    # -----------------------------------------------------------------#
    # read m(t) data
    # -----------------------------------------------------------------#
    mt_frame = pd.read_pickle(mt_file)
    if excludeInconsistencies:
        original_length = len(mt_frame.m.tolist())
        rejected_inconsistencies = ['H_nonzero', 'H_linear', 'not_converged']
        for incon in rejected_inconsistencies:
            if not mt_frame.empty:
                mt_frame = mt_frame[[incon not in mt_frame[
                    'inconsistencies'].tolist()[i] for i in range(
                        len(mt_frame['inconsistencies'].tolist()))]]
        if len(mt_frame.m.tolist()) < 0.95 * original_length or len(
                mt_frame.m.tolist()) < 200:  # in general, it should be 260
            # windows long. if it is much shorter, it is unclear where to
            # divide it.
            raise Exception('Too many inconsistencies')
    mt_list = mt_frame.m.tolist()
    return mt_list


def mad(data):
    return np.median(np.abs(data-np.median(data)))


def variance_mt(save_path, data_path, fit_method, kmin, kmax,
                windowsize, windowstep, dt, logepsilon=True, Engel=False):

    patients = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,
                           15, 16, 17, 18, 19, 20]
    if Engel:
        patients = [1, 6, 8, 9, 10, 11, 16, 17, 18, 19]
    laterality = ["SOZ", "nSOZ"]
    # -----------------------------------------------------------------#
    # Extract m(t) values an compute variances
    # -----------------------------------------------------------------#
    df = None
    for j, soz in enumerate(laterality):
        variance_first = []
        variance_last = []
        diff_list = []
        for ir, patient in enumerate(patients):
            recordings = mr_utilities.list_preictrecordings(patient)[0]
            for irec, rec in enumerate(recordings):
                try:
                    m_list = get_mt(patient, rec, soz, data_path, fit_method,
                                    kmin, kmax, windowsize, windowstep, dt,
                                    excludeInconsistencies=True)
                    m_first = m_list[:int(len(m_list)/2)]
                    m_last = m_list[-int(len(m_list) / 2):]

                    if logepsilon:
                        m_first = [np.log(1-m) for m in m_first]
                        m_last = [np.log(1 - m) for m in m_last]

                    variance_first.append(mad(m_first))
                    variance_last.append(mad(m_last))
                    diff_list.append(mad(m_first)-mad(m_last))
                except Exception:
                    continue

        # -----------------------------------------------------------------#
        # create boxplot of variance
        # -----------------------------------------------------------------#
        stat, p = stats.wilcoxon(diff_list)
        print('p-value:'+str(p))
        quantity = 'variance'

        data1 = pd.DataFrame({
            quantity: variance_first,
            'rec_type': ['10 - 5 min' for _ in range(len(variance_first))],
            'binning_focus': [soz for _ in range(len(variance_first))],
            'pval': [p for _ in range(len(variance_last))]
        })
        if df is None:
            df = data1
        else:
            df = df.append(data1, ignore_index=True)

        data2 = pd.DataFrame({
            quantity: variance_last,
            'rec_type': ['5 - 0 min' for _ in range(len(variance_last))],
            'binning_focus': [soz for _ in range(len(variance_last))],
            'pval': [p for _ in range(len(variance_last))]
        })
        df = df.append(data2, ignore_index=True)

    # -----------------------------------------------------------------#
    # Plotting
    # -----------------------------------------------------------------#
    fig, ax = plt.subplots(1, 1, figsize=(cm2inch(5.5), cm2inch(4.5)))
    g = sns.boxplot(ax=ax, x='binning_focus', y='variance', hue='rec_type',
                    data=df, hue_order=['10 - 5 min', '5 - 0 min'], order=[
                    'SOZ', 'nSOZ'], dodge=True, color='slategray', palette=[
                    'lightgray', 'slategray'], linewidth=1.1, width=0.6)

    g.set_yscale("log")
    plt.grid(axis='x', b=False)
    plt.rcParams['legend.title_fontsize'] = 9
    plt.legend(loc="upper left", ncol=1, fontsize='x-small',
               bbox_to_anchor=(0, 1.7), title='time to seizure onset')

    n_f_pre = len(df[(df['binning_focus'] == 'SOZ') &
                     (df['rec_type'] == '5 - 0 min')].index)
    n_f_inter = len(df[(df['binning_focus'] == 'SOZ') &
                       (df['rec_type'] == '10 - 5 min')].index)
    n_c_pre = len(df[(df['binning_focus'] == 'nSOZ') &
                     (df['rec_type'] == '5 - 0 min')].index)
    n_c_inter = len(df[(df['binning_focus'] == 'nSOZ') &
                       (df['rec_type'] == '10 - 5 min')].index)
    ax.set_xticklabels(['ipsi\n({}, {})'.format(n_f_pre, n_f_inter),
                        'contra\n({}, {})'.format(n_c_pre, n_c_inter)])

    # -----------------------------------------------------------------#
    # annotate p-values
    # -----------------------------------------------------------------#
    y_max = 2.4 * np.max(variance_first)
    y_min = 0.4 * np.min(variance_first)
    ax.annotate("", xy=(-0.2, y_max), xycoords='data',
                xytext=(0.2, y_max), textcoords='data',
                arrowprops=dict(arrowstyle="-", ec='black',
                                connectionstyle="bar,fraction=0.2"))

    if quantity == 'epsilon':
        ax.text(-0.25, y_max * 0.6,
                'p={:.3f}'.format(df[df['binning_focus']=='SOZ'].pval.tolist()[0]), fontsize=8)
    else:
        ax.text(-0.2, y_max + abs(y_max - y_min) * 0.1,
                'p={:.3f}'.format(df[df['binning_focus']=='SOZ'].pval.tolist()[0]), fontsize=8)

    ax.annotate("", xy=(0.8, y_max), xycoords='data', xytext=(1.2, y_max),
                textcoords='data', arrowprops=dict(arrowstyle="-", ec='black',
                connectionstyle="bar,fraction=0.2"))

    if quantity == 'epsilon':
        ax.text(0.75, y_max * 0.6,
                'p={:.3f}'.format(df[df['binning_focus']=='nSOZ'].pval.tolist()[0]), fontsize=8)
    else:
        ax.text(0.8, y_max + abs(y_max - y_min) * 0.1,
                'p={:.3f}'.format(df[df['binning_focus']=='nSOZ'].pval.tolist()[0]), fontsize=8)

    plt.xlabel('')
    sns.despine()
    plt.ylabel(r'Variability in $\hat{m}$')
    plt.yscale('log')
    plt.savefig(
     '{}Fig2_variance_mt_Engel{}_{}_winsize{}_winstep{}_kmin{}_kmax{}_epsilon{}.pdf'.format(
            save_path, Engel, fit_method, windowsize, windowstep, kmin,
            kmax, logepsilon), bbox_inches='tight')
    plt.show()


def timeresolved_example():
    result_path = '../Results/preictal/singleunits/'
    ksteps = (1, 400)
    kmin = ksteps[0]
    kmax = ksteps[1]
    windowsize = 20000
    windowstep = 500
    dt = 4
    patients = [13]
    for patient in patients:
        mt_examples(patient, result_path, 'offset', kmin, kmax,
                    windowsize, windowstep, dt, excludeInconsistencies=True)


def boxplot_2parts_hemispheres():
    x = 2
    ksteps = (1, 400)
    fit_method = 'offset'
    result_path = '../Results/preictal/singleunits/'
    boxplot_mx(result_path, x, ksteps, fit_method,
                excludeInconsistencies=True, Engel=False)
    boxplot_mx(result_path, x, ksteps, fit_method,
                excludeInconsistencies=True, Engel=True)


def boxplot_2parts_regions():
    x = 2
    ksteps = (1, 400)
    fit_method = 'offset'
    result_path = '../Results/preictal/singleunits_regions/'
    region = 'all'
    boxplot_mx_regions(region, result_path, x, ksteps, fit_method,
                              excludeInconsistencies=True)

    pointplot_mx_regions_Engel(region, result_path, x, ksteps,
                                      fit_method, excludeInconsistencies=True)


def visualize_mt_variance():
    save_path = '../Results/preictal/singleunits/'
    data_path = '../Results/preictal/singleunits/'
    variance_mt(save_path, data_path, fit_method='offset', kmin=1, kmax=400,
                windowsize=20000, windowstep=500, dt=4,
                Engel=False, logepsilon=True)


def main():
    timeresolved_example()
    #boxplot_2parts_hemispheres()
    #boxplot_2parts_regions()
    visualize_mt_variance()


if __name__ == '__main__':
    main()