binned_spiking_activity_human_epilepsy/python_code/analysis_PlosCB/create_figure_1.py

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

colors = ['steelblue', 'slategray', 'firebrick', 'darkblue', 'darkgreen']


def cm2inch(value):
    return value/2.54


def boxplot_focal_contra_onlypaired(pkl_file_pre, pkl_file_inter,
                                 resultpath, excludeInconsistencies, **kwargs):

    mr_results_pre = pd.read_pickle(pkl_file_pre)
    mr_results_inter = pd.read_pickle(pkl_file_inter)
    mr_results = pd.concat([mr_results_pre, mr_results_inter],
                           ignore_index=True)
    mr_results = mr_results[(mr_results['fit_method'] == 'offset')]
    conditions = list()

    for key in kwargs.keys():
        conditions.append(mr_results[key] == kwargs[key])
    if len(conditions) > 0:
        mr_results = mr_results[mr_utilities.conjunction(*conditions)]

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

    mr_results['logtau'] = mr_results.apply(lambda row: np.log(row['tau']),
                                            axis=1)
    mr_results['logepsilon'] = mr_results.apply(lambda row: abs(np.log(1-row[
        'm'])), axis=1)
    mr_results['epsilon'] = mr_results.apply(lambda row: 1-row['m'], axis=1)
    quantity = 'epsilon'
    # ================================================================== #
    # compute pvalue (paired, just uses data with pair)
    # ================================================================== #
    truediff = []
    paired_recordings = []
    for rec in mr_results.recording.unique():
        recframe = mr_results[mr_results['recording']==rec]
        if not len(recframe['m'].tolist()) == 2:
            continue
        focalval = recframe[recframe['binning_focus']=='focal'][quantity].item()
        contraval = recframe[recframe['binning_focus'] == 'contra-lateral'][
            quantity].item()
        truediff.append(focalval - contraval)
        paired_recordings.append(rec)

    mr_results_paired = mr_results[mr_results['recording'].isin(paired_recordings)]
    p_hierarchical = hierarchical_model.hierarchical_model_hemispheres(
        mr_results_paired)

    # ================================================================== #
    # Plotting
    # ================================================================== #
    fig, ax = plt.subplots(1, 1, figsize=(cm2inch(4), cm2inch(3)))

    gb = sns.boxplot(ax=ax, x="binning_focus", y=quantity,
                  data=mr_results_paired,  order=['focal', 'contra-lateral'],
                  dodge=True, palette=['darkblue', 'slategray'])

    for patch in ax.artists:
        r, g, b, a = patch.get_facecolor()
        patch.set_facecolor((r, g, b, .4))
        patch.set_edgecolor('gray')

    sns.stripplot(ax=ax, x="binning_focus", y=quantity,
                  data=mr_results_paired, alpha=0.9, order=['focal',
                  'contra-lateral'], size=2, dodge=True, palette=[
                  'darkblue', 'slategray'], jitter=False)
    gb.set_yscale("log")
    plt.grid(axis='x', b=False)
    n_SOZ = len(mr_results_paired[mr_results_paired['binning_focus']=='focal']['m'].tolist())
    n_nSOZ = len(
        mr_results_paired[mr_results_paired['binning_focus'] == 'contra-lateral'][
            'm'].tolist())
    ax.set_xticklabels(['ipsi\n({})'.format(n_SOZ), 'contra\n({})'.format(
        n_nSOZ)])
    ax.set_xlabel('')
    quantity_label = r'$\hat{m}$'
    ymin = 1.4 * np.max(mr_results[quantity].tolist())
    ymax = 0.6*np.min(mr_results[quantity].tolist())
    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])
    plt.ylabel(quantity_label)
    sns.despine()
    # ================================================================== #
    # Annotate pvalues
    # ================================================================== #
    x1, x2 = 0, 1
    y, h = ymax * (1.1), -0.1 * ymax
    ax.plot([x1, x1, x2, x2], [y, y + h, y + h, y], lw=1.5, c='black')
    ax.text((x1 + x2) * .5, y + 4*h, r"$p$"+"={:.3f}".format(
        p_hierarchical), ha='center', va='bottom', color='black', fontsize=10)

    # ================================================================== #
    # Final figure adjustments + saving
    # ================================================================== #
    sns.despine(fig)
    fig.subplots_adjust(bottom=None, right=None,
                        wspace=0.8, hspace=0.5, left=0.15, top=1.2)
    kmin = mr_results["kmin"].unique()[0]
    kmax = mr_results["kmax"].unique()[0]
    resultfile = '{}Fig1_boxplot_{}_kmin{}_kmax{}_onlypaired.pdf'.format(
        resultpath,quantity, kmin, kmax)
    plt.savefig(resultfile, bbox_inches='tight')
    plt.show()
    plt.close()


def boxplot_focal_contra_Engel(pkl_file_pre, pkl_file_inter,
                               resultpath, excludeInconsistencies, **kwargs):

    mr_results_pre = pd.read_pickle(pkl_file_pre)
    mr_results_inter = pd.read_pickle(pkl_file_inter)
    mr_results = pd.concat([mr_results_pre, mr_results_inter],
                           ignore_index=True)

    patients1a = [1, 6, 8, 9, 10, 11, 16, 17, 18, 19]
    mr_results = mr_results[(mr_results['patient'].isin(patients1a))]
    mr_results = mr_results[(mr_results['fit_method'] == 'offset')]

    conditions = list()
    for key in kwargs.keys():
        conditions.append(mr_results[key] == kwargs[key])
    if len(conditions) > 0:
        mr_results = mr_results[mr_utilities.conjunction(*conditions)]

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

    mr_results['logtau'] = mr_results.apply(lambda row: np.log(row['tau']),
                                            axis=1)
    mr_results['logepsilon'] = mr_results.apply(lambda row: abs(np.log(1-row[
        'm'])), axis=1)
    mr_results['epsilon'] = mr_results.apply(lambda row: 1-row[
        'm'], axis=1)
    quantity = 'epsilon'
    # ================================================================== #
    # compute pvalue (unpaired, all data)
    # ================================================================== #
    p_hierarchical = hierarchical_model.hierarchical_model_hemispheres(
        mr_results)

    # ================================================================== #
    # Plotting
    # ================================================================== #
    fig, ax = plt.subplots(1, 1, figsize=(cm2inch(4), cm2inch(3)))

    gb = sns.boxplot(ax=ax, x="binning_focus", y=quantity,
                  data=mr_results,  order=['focal', 'contra-lateral'],
                  dodge=True, palette=['darkblue', 'slategray'])

    for patch in ax.artists:
        r, g, b, a = patch.get_facecolor()
        patch.set_facecolor((r, g, b, .4))
        patch.set_edgecolor('gray')

    sns.stripplot(ax=ax, x="binning_focus", y=quantity,
                  data=mr_results, alpha=0.9, order=['focal',
                                                     'contra-lateral'],
                  size=2, dodge=True, palette=['darkblue', 'slategray'], jitter=False)
    gb.set_yscale("log")
    plt.grid(axis='x', b=False)
    n_SOZ = len(mr_results[mr_results['binning_focus']=='focal']['m'].tolist())
    n_nSOZ = len(mr_results[mr_results['binning_focus'] == 'contra-lateral'][
            'm'].tolist())
    ax.set_xticklabels(['ipsi\n({})'.format(n_SOZ), 'contra\n({})'.format(
        n_nSOZ)])
    ax.set_xlabel('')
    quantity_label = r'$\hat{m}$'
    ymin = 1.4 * np.max(mr_results[quantity].tolist())
    ymax = 0.6*np.min(mr_results[quantity].tolist())
    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])
    plt.ylabel(quantity_label)
    sns.despine()
    # ================================================================== #
    # Annotate pvalues
    # ================================================================== #
    x1, x2 = 0, 1
    y, h = ymax * (1.1), -0.1 * ymax
    ax.plot([x1, x1, x2, x2], [y, y + h, y + h, y], lw=1.5, c='black')
    ax.text((x1 + x2) * .5, y + 4*h, r"$p$"+"={:.3f}".format(p_hierarchical),
            ha='center', va='bottom', color='black', fontsize=10)

    # ================================================================== #
    # Final figure adjustments + saving
    # ================================================================== #
    sns.despine(fig)
    fig.subplots_adjust(bottom=None, right=None,
                        wspace=0.8, hspace=0.5, left=0.15, top=1.2)
    kmin = mr_results["kmin"].unique()[0]
    kmax = mr_results["kmax"].unique()[0]
    resultfile = '{}Fig1_boxplot_{}_kmin{}_kmax{}1aEngel.pdf'.format(
        resultpath,quantity, kmin, kmax)
    plt.savefig(resultfile, bbox_inches='tight')
    plt.show()
    plt.close()


def boxplot_focal_contra(pkl_file_pre, pkl_file_inter, resultpath,
                         excludeInconsistencies, **kwargs):

    mr_results_pre = pd.read_pickle(pkl_file_pre)
    mr_results_inter = pd.read_pickle(pkl_file_inter)
    mr_results = pd.concat([mr_results_pre, mr_results_inter],
                           ignore_index=True)

    mr_results = mr_results[(mr_results['fit_method'] == 'offset')]
    conditions = list()
    for key in kwargs.keys():
        conditions.append(mr_results[key] == kwargs[key])
    if len(conditions) > 0:
        mr_results = mr_results[mr_utilities.conjunction(*conditions)]

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

    mr_results['logtau'] = mr_results.apply(lambda row: np.log(row['tau']),
                                            axis=1)
    mr_results['logepsilon'] = mr_results.apply(lambda row: abs(np.log(1-row[
                                'm'])), axis=1)
    mr_results['epsilon'] = mr_results.apply(lambda row: 1-row['m'], axis=1)
    quantity = 'epsilon'
    # ================================================================== #
    # compute p-value (unpaired, all data)
    # ================================================================== #
    p_hierarchical = hierarchical_model.hierarchical_model_hemispheres(mr_results)

    # ================================================================== #
    # Plotting
    # ================================================================== #
    fig, ax = plt.subplots(1, 1, figsize=(cm2inch(4), cm2inch(3)))

    gb = sns.boxplot(ax=ax, x="binning_focus", y=quantity,
                  data=mr_results,  order=['focal', 'contra-lateral'],
                  dodge=True, palette=['darkblue', 'slategray'])

    for patch in ax.artists:
        r, g, b, a = patch.get_facecolor()
        patch.set_facecolor((r, g, b, .4))
        patch.set_edgecolor('gray')

    sns.stripplot(ax=ax, x="binning_focus", y=quantity,
                  data=mr_results, alpha=0.9, order=['focal',
                                                     'contra-lateral'],
                  size=2, dodge=True, palette=['darkblue', 'slategray'],
                  jitter=False)

    gb.set_yscale("log")
    plt.grid(axis='x', b=False)
    n_SOZ = len(mr_results[mr_results['binning_focus'] == 'focal'][
                    'm'].tolist())
    n_nSOZ = len(mr_results[mr_results['binning_focus'] == 'contra-lateral'][
            'm'].tolist())
    ax.set_xticklabels(['ipsi\n({})'.format(n_SOZ), 'contra\n({})'.format(
                        n_nSOZ)])
    ax.set_xlabel('')

    quantity_label = r'$\hat{m}$'
    ymin = 1.4 * np.max(mr_results[quantity].tolist())
    ymax = 0.6*np.min(mr_results[quantity].tolist())
    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])
    plt.ylabel(quantity_label)
    sns.despine()

    # ================================================================== #
    # Annotate pvalues
    # ================================================================== #
    x1, x2 = 0, 1
    y, h = ymax * (1.1), -0.1 * ymax
    ax.plot([x1, x1, x2, x2], [y, y + h, y + h, y], lw=1.5, c='black')
    ax.text((x1 + x2) * .5, y + 4*h, r"$p$"+"={:.3f}".format(p_hierarchical),
            ha='center', va='bottom', color='black', fontsize=10)

    # ================================================================== #
    # Final figure adjustments + saving
    # ================================================================== #
    sns.despine(fig)
    fig.subplots_adjust(bottom=None, right=None,
                        wspace=0.8, hspace=0.5, left=0.15, top=1.2)
    kmin = mr_results["kmin"].unique()[0]
    kmax = mr_results["kmax"].unique()[0]
    resultfile = '{}Fig1_boxplot_{}_kmin{}_kmax{}.pdf'.format(
        resultpath,quantity, kmin, kmax)
    plt.savefig(resultfile, bbox_inches='tight')
    plt.show()
    plt.close()


def patientwise_pointplots(pkl_file_pre, pkl_file_inter,
                            result_path, excludeInconsistencies, **kwargs):
    # -----------------------------------------------------------------#
    # Read data
    # -----------------------------------------------------------------#
    mr_results_inter = pd.read_pickle(pkl_file_inter)
    mr_results_inter['rec_type'] = mr_results_inter.apply(lambda row: 'inter',
                                                          axis=1)
    mr_results_pre = pd.read_pickle(pkl_file_pre)
    mr_results_pre['rec_type'] = mr_results_pre.apply(lambda row: 'pre',
                                                          axis=1)
    mr_results = pd.concat([mr_results_pre, mr_results_inter],
                           ignore_index=True)
    mr_results = mr_results[(mr_results['fit_method'] == 'offset')]
    conditions = list()
    for key in kwargs.keys():
        conditions.append(mr_results[key] == kwargs[key])
    if len(conditions) > 0:
        mr_results = mr_results[mr_utilities.conjunction(*conditions)]

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

    mr_results['logepsilon'] = mr_results.apply(lambda row: abs(
        np.log(1-row['m'])), axis=1)
    mr_results['epsilon'] = mr_results.apply(lambda row: 1-row['m'], axis=1)
    quantity = 'epsilon'
    # -----------------------------------------------------------------#
    # Run through all patients that also have preictal data
    # get their interictal and preictal results plot them as a pointplot
    # -----------------------------------------------------------------#
    patients = mr_results['patient'].unique()
    patients_sorted = np.sort(patients)

    fig, ax = plt.subplots(5, 4, figsize=(12, 1.5 * 5))
    v = [0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3]
    h = [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4]
    for ip, patient in enumerate(patients_sorted):
        mr_patient = mr_results[mr_results['patient'] == patient]
        recordings = mr_patient.recording.unique()

        all_SOZ = mr_patient[mr_patient['binning_focus'] == 'focal'][
                    quantity].tolist()
        all_nSOZ = mr_patient[mr_patient['binning_focus'] == 'contra-lateral'][
                    quantity].tolist()
        stat, p_unpaired = stats.mannwhitneyu(all_SOZ, all_nSOZ, alternative='two-sided')

        for rec in recordings:
            mr_rec = mr_patient[mr_patient['recording'] == rec]
            if len(mr_rec[quantity].tolist()) == 2:
                q_SOZ = mr_rec[mr_rec['binning_focus'] == 'focal'][
                    quantity].item()
                q_nSOZ = mr_rec[mr_rec['binning_focus'] == 'contra-lateral'][
                    quantity].item()
                color = 'darkorange'
                if mr_rec['rec_type'].tolist()[0] == 'pre':
                    color = 'indianred'
                ax[h[ip], v[ip]].plot([1, 2], [q_SOZ, q_nSOZ], c=color,
                                      marker='.', markersize=6)
            elif len(mr_rec[mr_rec['binning_focus']=='focal'].index) == 1:

                q_SOZ = mr_rec[mr_rec['binning_focus']=='focal'][
                    quantity].item()
                color = 'darkorange'
                if mr_rec['rec_type'].tolist()[0] == 'pre':
                    color = 'indianred'
                ax[h[ip], v[ip]].plot([1], [q_SOZ], c=color,
                                      marker='.', markersize=6)
            elif len(mr_rec[mr_rec['binning_focus']=='contra-lateral'].index) \
                    == 1:
                q_nSOZ = mr_rec[mr_rec['binning_focus']=='contra-lateral'][
                    quantity].item()
                color = 'darkorange'
                if mr_rec['rec_type'].tolist()[0] == 'pre':
                    color = 'indianred'
                ax[h[ip], v[ip]].plot([2], [q_nSOZ], c=color,
                                      marker='.', markersize=6)

        if len(all_SOZ)+len(all_nSOZ) >= 8:
            ax[h[ip], v[ip]].text(0.5, 0.9,'p={:.3f}'.format(p_unpaired),
                 horizontalalignment='center', verticalalignment='center',
                 transform=ax[h[ip], v[ip]].transAxes, color='black',
                 fontsize=10)

        ax[h[ip], v[ip]].set_title('Patient {}'.format(patient))
        ax[h[ip], v[ip]].set_yscale("log")
        quantity_label = r'$\hat{m}$'
        ax[h[ip], v[ip]].set_ylabel(quantity_label)
        ymin = 1.1
        ymax = 1e-3
        ax[h[ip], v[ip]].set_ylim(ymin, ymax)
        ax[h[ip], v[ip]].set_xlim(0.7, 2.3)
        T = [1, 2]
        ax[h[ip], v[ip]].set_xticks(T)
        ax[h[ip], v[ip]].set_xticklabels(['ipsi', 'contra'])
        T = [1, 1e-1, 1e-2, 1e-3]
        ax[h[ip], v[ip]].set_yticks(T)
        ax[h[ip], v[ip]].set_yticklabels([0, 0.9, 0.99, 0.999])
        ax[h[ip], v[ip]].set_xlabel(' ')
        ax[h[ip], v[ip]].set_xticklabels(['ipsi', 'contra'])

    sns.despine(fig)
    fig.subplots_adjust(bottom=None, right=None,
                        wspace=0.7, hspace=0.7, left=0.15, top=1.2)

    kmin = mr_results_pre['kmin'].unique()[0]
    kmax = mr_results_pre['kmax'].unique()[0]

    plotfile_mx = '{}S2_{}_patients_kmin{}_kmax{}.pdf'.format(
        result_path, quantity, kmin, kmax)
    fig.savefig(plotfile_mx, bbox_inches='tight')
    plt.show()


def boxplot_brainregions_focus_onlypaired(pkl_file_pre, pkl_file_inter,
                                          result_path,
                                          excludeInconsistencies, **kwargs):

    regions = ['H', 'A', 'PHC', 'EC']
    mr_results_pre = pd.read_pickle(pkl_file_pre)
    mr_results_inter = pd.read_pickle(pkl_file_inter)
    mr_results = pd.concat([mr_results_pre, mr_results_inter],
                           ignore_index=True)

    mr_results = mr_results[(mr_results['fit_method'] == 'offset')]
    conditions = list()
    for key in kwargs.keys():
        conditions.append(mr_results[key] == kwargs[key])
    if len(conditions) > 0:
        mr_results = mr_results[mr_utilities.conjunction(*conditions)]

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

    # -----------------------------------------------------------------#
    # Write the brainregion, independent of hemisphere, into frame
    # -----------------------------------------------------------------#
    def give_region(row):
        return row['binning'][1:]

    # -----------------------------------------------------------------#
    # Write the brainregion and focus into frame
    # -----------------------------------------------------------------#
    def give_region_and_focus(row):
        regionname = row['binning'][1:]
        focus = row['binning_focus'][0] # f for focal and c for contralateral
        return '{}{}'.format(focus, regionname)

    mr_results['region'] = mr_results.apply(lambda row: give_region(row), axis=1)
    mr_results['regionfocus'] = mr_results.apply(lambda row: give_region_and_focus(row), axis=1)
    mr_results['logtau'] = mr_results.apply(lambda row: np.log(row['tau']),
                                            axis=1)
    mr_results['logepsilon'] = mr_results.apply(lambda row: abs(np.log(
        1-row['m'])), axis=1)
    mr_results['epsilon'] = mr_results.apply(lambda row: 1-row['m'], axis=1)
    quantity = 'epsilon'
    # -----------------------------------------------------------------#
    # Filter to recordings and regions for which paired data exists
    # -----------------------------------------------------------------#
    mr_results_paired = pd.DataFrame(columns=mr_results.columns)
    differences = dict({'A':[], 'H':[], 'PHC':[], 'EC':[]})

    for rec in mr_results.recording.unique():
        recframe = mr_results[mr_results['recording']==rec]
        for area in recframe['region'].unique():
            area_frame = recframe[recframe['region'] == area]
            if not len(area_frame['m'].tolist()) == 2:
                continue
            focalval = area_frame[area_frame['binning_focus'] == 'focal'][
                quantity].item()
            contraval = \
            area_frame[area_frame['binning_focus'] == 'contra-lateral'][
                quantity].item()
            differences[area].append(focalval - contraval)
            mr_results_paired = mr_results_paired.append(area_frame,
                                                       ignore_index=True)
    mr_results = mr_results_paired

    # -----------------------------------------------------------------#
    # Compute p-values of pairwise comparison
    # -----------------------------------------------------------------#
    pvals = dict()
    medians = dict()

    for area in mr_results['region'].unique():
        region_differences = differences[area]
        stat, p_paired = stats.wilcoxon(region_differences)  # wilcoxon is by
        # default two-sided
        pvals[area] = p_paired
        medians['f{}'.format(area)] = np.median(focalval)
        medians['c{}'.format(area)] = np.median(contraval)

    global_pvals, pvals_regions = hierarchical_model.hierarchical_model_regions(
        mr_results)
    pvals = pvals_regions

    # -----------------------------------------------------------------#
    # Boxplot of regions
    # -----------------------------------------------------------------#
    fig, ax = plt.subplots(1, 1, figsize=(cm2inch(10*1.2), cm2inch(3*1.2)))
    gb = sns.boxplot(ax=ax, x="region", y=quantity, hue='binning_focus',
                  data=mr_results,  order=regions, hue_order=['focal',
                                                     'contra-lateral'],
                  dodge=True, palette=['darkblue', 'slategray'], linewidth=1)
    for patch in ax.artists:
        r, g, b, a = patch.get_facecolor()
        patch.set_facecolor((r, g, b, .4))
        patch.set_edgecolor('gray')

    sns.stripplot(ax=ax, x="region", y=quantity, hue='binning_focus',
                  data=mr_results, alpha=0.9, order=regions,
                  hue_order=['focal', 'contra-lateral'], jitter=False,
                  size=2, palette=['darkblue', 'slategray'], dodge=True)

    ax.get_legend().set_visible(False)
    plt.grid(axis='x', b=False)
    gb.set_yscale("log")
    labels = []
    for ireg, reg in enumerate(regions):
        mr_results_r = mr_results[mr_results['region']==reg]
        n_reg_f = len(mr_results_r[mr_results_r['binning_focus']=='focal'][
                          'm'].tolist())
        n_reg_c = len(mr_results_r[mr_results_r[
                                       'binning_focus']=='contra-lateral'][
                          'm'].tolist())
        labels.append('{}\n({}, {})'.format(reg, n_reg_f, n_reg_c))
    ax.set_xticklabels(labels)
    quantity_label = r'$\hat{m}$'
    ymin = 1.4 * np.max(mr_results[quantity].tolist())
    ymax = 0.4 * np.min(mr_results[quantity].tolist())
    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])
    plt.ylabel(quantity_label)
    sns.despine()
    # -----------------------------------------------------------------#
    # Annotate p-values of pairwise comparison
    # -----------------------------------------------------------------#
    y_max = 1.3*ymax
    y_min = ymin
    for ir, region in enumerate(regions):
        ax.annotate("", xy=(-0.2+ir*1.0, 1.5*y_max), xycoords='data',
                    xytext=(0.2+ir*1.0, 1.5*y_max), textcoords='data',
                    arrowprops=dict(arrowstyle="-", ec='black',
                                    connectionstyle="bar,fraction=0.3"))
        ax.text(-0.35+ir*1.0, y_max* 0.9- abs(y_max - y_min) * 0.0,
                'p'+'={:.3f}'.format(pvals[region]), fontsize=8.5)
    # -----------------------------------------------------------------#
    # Adjust some plot settings and save figures
    # -----------------------------------------------------------------#
    ax.set_xlabel('brain area')
    ax.set_ylim((ymin, ymax))
    ax.set_title('overall model: ' + str(global_pvals), fontsize=4)
    fig.subplots_adjust(bottom=None, right=None,
                        wspace=0.2, hspace=0.5, left=0.15, top=1.2)
    kmin = mr_results['kmin'].unique()[0]
    kmax = mr_results['kmax'].unique()[0]
    plotfile_mx = '{}Fig1_regions_{}_kmin{}_kmax{}_exclude_{}_onlypaired.pdf'.format(
        result_path, quantity, kmin, kmax, excludeInconsistencies)
    fig.savefig(plotfile_mx, bbox_inches='tight')
    plt.show()
    plt.close('all')


def boxplot_brainregions_focus(pkl_file_pre, pkl_file_inter, result_path,
                               excludeInconsistencies, **kwargs):

    regions = ['H', 'A', 'PHC', 'EC']
    mr_results_pre = pd.read_pickle(pkl_file_pre)
    mr_results_inter = pd.read_pickle(pkl_file_inter)
    mr_results = pd.concat([mr_results_pre, mr_results_inter],
                           ignore_index=True)

    mr_results = mr_results[(mr_results['fit_method'] == 'offset')]
    conditions = list()
    for key in kwargs.keys():
        conditions.append(mr_results[key] == kwargs[key])
    if len(conditions) > 0:
        mr_results = mr_results[mr_utilities.conjunction(*conditions)]

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

    # -----------------------------------------------------------------#
    # Write the brainregion, independent of hemisphere, into frame
    # -----------------------------------------------------------------#
    def give_region(row):
        return row['binning'][1:]

    # -----------------------------------------------------------------#
    # Write the brainregion and focus into frame
    # -----------------------------------------------------------------#
    def give_region_and_focus(row):
        regionname = row['binning'][1:]
        focus = row['binning_focus'][0] # f for focal and c for contralateral
        return '{}{}'.format(focus, regionname)

    mr_results['region'] = mr_results.apply(lambda row: give_region(row), axis=1)
    mr_results['regionfocus'] = mr_results.apply(lambda row: give_region_and_focus(row), axis=1)
    mr_results['logtau'] = mr_results.apply(lambda row: np.log(row['tau']),axis=1)
    mr_results['logepsilon'] = mr_results.apply(lambda row: -np.log(1-row['m']), axis=1)
    mr_results['epsilon'] = mr_results.apply(lambda row: 1-row['m'], axis=1)
    quantity = 'epsilon'
    # -----------------------------------------------------------------#
    # Compute p-values of pairwise comparison
    # -----------------------------------------------------------------#
    pvals = dict()
    medians = dict()

    for area in mr_results['region'].unique():
        area_frame = mr_results[mr_results['region'] == area]
        f_area_frame = area_frame[area_frame['binning_focus'] == 'focal']
        c_area_frame = area_frame[area_frame['binning_focus'] ==
                                 'contra-lateral']
        focalval = f_area_frame[quantity].tolist()
        contraval = c_area_frame[quantity].tolist()

        stat, p = stats.mannwhitneyu(focalval, contraval,
                                     alternative='two-sided')
        pvals[area] = p
        medians['f{}'.format(area)] = np.median(focalval)
        medians['c{}'.format(area)] = np.median(contraval)

    global_pvals, pvals_regions = hierarchical_model.hierarchical_model_regions(
        mr_results)
    pvals = pvals_regions
    # -----------------------------------------------------------------#
    # Boxplot of regions
    # -----------------------------------------------------------------#
    fig, ax = plt.subplots(1, 1, figsize=(cm2inch(10*1.2), cm2inch(3*1.2)))

    gb = sns.boxplot(ax=ax, x="region", y=quantity, hue='binning_focus',
                  data=mr_results,  order=regions, hue_order=['focal',
                                                     'contra-lateral'],
                  dodge=True, palette=['darkblue', 'slategray'], linewidth=1)
    gb.set_yscale("log")
    for patch in ax.artists:
        r, g, b, a = patch.get_facecolor()
        patch.set_facecolor((r, g, b, .4))
        patch.set_edgecolor('gray')

    sns.stripplot(ax=ax, x="region", y=quantity, hue='binning_focus',
                  data=mr_results, alpha=0.9, order=regions,
                  hue_order=['focal', 'contra-lateral'], jitter=False,
                  size=2, palette=['darkblue', 'slategray'], dodge=True)

    ax.get_legend().set_visible(False)
    plt.grid(axis='x', b=False)

    labels = []
    for ireg, reg in enumerate(regions):
        mr_results_r = mr_results[mr_results['region']==reg]
        n_reg_f = len(mr_results_r[mr_results_r['binning_focus']=='focal'][
                          'm'].tolist())
        n_reg_c = len(mr_results_r[mr_results_r[
                                       'binning_focus']=='contra-lateral'][
                          'm'].tolist())
        labels.append('{}\n({}, {})'.format(reg, n_reg_f, n_reg_c))
    ax.set_xticklabels(labels)

    quantity_label = r'$\hat{m}$'
    ymin = 1.4 * np.max(mr_results[quantity].tolist())
    ymax = 0.4 * np.min(mr_results[quantity].tolist())
    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])
    plt.ylabel(quantity_label)
    sns.despine()

    # -----------------------------------------------------------------#
    # Annotate p-values of pairwise comparison
    # -----------------------------------------------------------------#
    y_max = 1.3*ymax
    y_min = ymin
    for ir, region in enumerate(regions):

        ax.annotate("", xy=(-0.2+ir*1.0, 1.5*y_max), xycoords='data',
                    xytext=(0.2+ir*1.0, 1.5*y_max), textcoords='data',
                    arrowprops=dict(arrowstyle="-", ec='black',
                                    connectionstyle="bar,fraction=0.3"))
        ax.text(-0.35+ir*1.0, y_max* 0.9- abs(y_max - y_min) * 0.0,
                'p'+'={:.3f}'.format(pvals[region]), fontsize=8.5)

    # -----------------------------------------------------------------#
    # Adjust some plot settings and save figures
    # -----------------------------------------------------------------#
    ax.set_xlabel('brain area')
    ax.set_ylim((ymin, ymax))
    ax.set_title('overall model: ' + str(global_pvals), fontsize=4)
    fig.subplots_adjust(bottom=None, right=None,
                        wspace=0.2, hspace=0.5, left=0.15, top=1.2)
    kmin = mr_results['kmin'].unique()[0]
    kmax = mr_results['kmax'].unique()[0]
    plotfile_mx = '{}Fig1_regions_{}_kmin{}_kmax{}_exclude_{}.pdf'.format(
        result_path, quantity, kmin, kmax, excludeInconsistencies)
    fig.savefig(plotfile_mx, bbox_inches='tight')
    plt.show()
    plt.close('all')


def boxplot_brainregions_focus_Engel(pkl_file_pre,pkl_file_inter,result_path,
                                     excludeInconsistencies, **kwargs):

    regions = ['H', 'A', 'PHC', 'EC']
    mr_results_pre = pd.read_pickle(pkl_file_pre)
    mr_results_inter = pd.read_pickle(pkl_file_inter)
    mr_results = pd.concat([mr_results_pre, mr_results_inter],
                           ignore_index=True)
    patients1a = [1, 6, 8, 9, 10, 11, 16, 17, 18, 19]
    mr_results = mr_results[(mr_results['patient'].isin(patients1a))]
    mr_results = mr_results[(mr_results['fit_method'] == 'offset')]
    conditions = list()
    for key in kwargs.keys():
        conditions.append(mr_results[key] == kwargs[key])
    if len(conditions) > 0:
        mr_results = mr_results[mr_utilities.conjunction(*conditions)]

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

    # -----------------------------------------------------------------#
    # Write the brainregion, independent of hemisphere, into frame
    # -----------------------------------------------------------------#
    def give_region(row):
        return row['binning'][1:]

    # -----------------------------------------------------------------#
    # Write the brainregion and focus into frame
    # -----------------------------------------------------------------#
    def give_region_and_focus(row):
        regionname = row['binning'][1:]
        focus = row['binning_focus'][0] # f for focal and c for contralateral
        return '{}{}'.format(focus, regionname)

    mr_results['region'] = mr_results.apply(lambda row: give_region(row), axis=1)
    mr_results['regionfocus'] = mr_results.apply(lambda row: give_region_and_focus(row), axis=1)
    mr_results['logtau'] = mr_results.apply(lambda row: np.log(row['tau']),axis=1)
    mr_results['logepsilon'] = mr_results.apply(lambda row: abs(np.log(1-row['m'])), axis=1)
    mr_results['epsilon'] = mr_results.apply(lambda row: 1-row['m'], axis=1)
    quantity = 'epsilon'
    # -----------------------------------------------------------------#
    # Compute p-values of pairwise comparison
    # -----------------------------------------------------------------#
    pvals = dict()
    medians = dict()

    for area in mr_results['region'].unique():
        area_frame = mr_results[mr_results['region'] == area]
        f_area_frame = area_frame[area_frame['binning_focus'] == 'focal']
        c_area_frame = area_frame[area_frame['binning_focus'] ==
                                 'contra-lateral']
        focalval = f_area_frame[quantity].tolist()
        contraval = c_area_frame[quantity].tolist()

        stat, p = stats.mannwhitneyu(focalval, contraval,
                                     alternative='two-sided')
        pvals[area] = p
        medians['f{}'.format(area)] = np.median(focalval)
        medians['c{}'.format(area)] = np.median(contraval)

    global_pvals, pvals_regions = hierarchical_model.hierarchical_model_regions(
        mr_results)
    pvals = pvals_regions

    # -----------------------------------------------------------------#
    # Boxplot of regions
    # -----------------------------------------------------------------#
    fig, ax = plt.subplots(1, 1, figsize=(cm2inch(10*1.2), cm2inch(3*1.2)))
    gb = sns.boxplot(ax=ax, x="region", y=quantity, hue='binning_focus',
                  data=mr_results,  order=regions, hue_order=['focal',
                                                     'contra-lateral'],
                  dodge=True, palette=['darkblue', 'slategray'], linewidth=1)
    for patch in ax.artists:
        r, g, b, a = patch.get_facecolor()
        patch.set_facecolor((r, g, b, .4))
        patch.set_edgecolor('gray')

    sns.stripplot(ax=ax, x="region", y=quantity, hue='binning_focus',
                  data=mr_results, alpha=0.9, order=regions,
                  hue_order=['focal', 'contra-lateral'], jitter=False,
                  size=2, palette=['darkblue', 'slategray'], dodge=True)

    ax.get_legend().set_visible(False)
    plt.grid(axis='x', b=False)
    gb.set_yscale("log")
    labels = []
    for ireg, reg in enumerate(regions):
        mr_results_r = mr_results[mr_results['region']==reg]
        n_reg_f = len(mr_results_r[mr_results_r['binning_focus']=='focal'][
                          'm'].tolist())
        n_reg_c = len(mr_results_r[mr_results_r['binning_focus']==
                                   'contra-lateral']['m'].tolist())
        labels.append('{}\n({}, {})'.format(reg, n_reg_f, n_reg_c))
    ax.set_xticklabels(labels)

    quantity_label = r'$\hat{m}$'
    ymin = 1.4 * np.max(mr_results[quantity].tolist())
    ymax = 0.4 * np.min(mr_results[quantity].tolist())
    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])
    plt.ylabel(quantity_label)
    sns.despine()

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

    # -----------------------------------------------------------------#
    # Adjust some plot settings and save figures
    # -----------------------------------------------------------------#
    ax.set_xlabel('brain area')
    ax.set_ylim((ymin, ymax))
    ax.set_title('overall model: ' + str(global_pvals), fontsize=4)
    fig.subplots_adjust(bottom=None, right=None,
                        wspace=0.2, hspace=0.5, left=0.15, top=1.2)
    kmin = mr_results['kmin'].unique()[0]
    kmax = mr_results['kmax'].unique()[0]

    plotfile_mx = '{}Fig1_regions_{}_kmin{}_kmax{}_exclude_{}_Engel1a.pdf'.format(
        result_path, quantity, kmin, kmax, excludeInconsistencies)
    fig.savefig(plotfile_mx, bbox_inches='tight')
    plt.show()
    plt.close('all')


def boxplot_hemispheres():
    result_path = '../Results/preictal/singleunits/'
    pkl_file = result_path + 'mr_results.pkl'
    ksteps = (1, 400)
    pkl_file_pre = pkl_file
    result_path_inter = '../Results/interictal/singleunits/'
    pkl_file_inter = result_path_inter + 'mr_results.pkl'

    boxplot_focal_contra(pkl_file_pre, pkl_file_inter,
                         result_path, excludeInconsistencies=True,
                         kmin=ksteps[0], kmax=ksteps[1])
    boxplot_focal_contra_Engel(pkl_file_pre, pkl_file_inter,
                         result_path, excludeInconsistencies=True,
                         kmin=ksteps[0], kmax=ksteps[1])

    boxplot_focal_contra_onlypaired(pkl_file_pre, pkl_file_inter,
                         result_path, excludeInconsistencies=True,
                         kmin=ksteps[0], kmax=ksteps[1])


def pointplot_patients_hemispheres():
    pkl_file_pre = '../Results/preictal/singleunits/mr_results.pkl'
    pkl_file_inter = '../Results/interictal/singleunits/mr_results.pkl'
    result_path_pointplot = '../Results/preictal/singleunits/'

    patientwise_pointplots(pkl_file_pre, pkl_file_inter,
                           result_path_pointplot,
                           excludeInconsistencies=True,
                           kmin=1, kmax=400,
                           fit_method='offset')


def boxplot_regions():
    result_path = '../Results/preictal/singleunits_regions/'
    pkl_file = result_path + 'mr_results.pkl'

    result_path_inter = '../Results/interictal/singleunits_regions/'
    pkl_file_inter = result_path_inter + 'mr_results.pkl'
    ksteps = (1, 400)

    boxplot_brainregions_focus(pkl_file, pkl_file_inter,
                                            result_path,
                                            kmin=ksteps[0],
                                            kmax=ksteps[1],
                                            excludeInconsistencies=True,
                                            fit_method='offset')

    boxplot_brainregions_focus_Engel(pkl_file, pkl_file_inter,
                                            result_path,
                                            kmin=ksteps[0],
                                            kmax=ksteps[1],
                                            excludeInconsistencies=True,
                                            fit_method='offset')

    boxplot_brainregions_focus_onlypaired(pkl_file,
                                            pkl_file_inter,
                                            result_path,
                                            kmin=ksteps[0],
                                            kmax=ksteps[1],
                                            excludeInconsistencies=True,
                                            fit_method='offset')


def main():
    boxplot_hemispheres()
    pointplot_patients_hemispheres()
    boxplot_regions()


if __name__ == '__main__':
    main()