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):

    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_utilities.exclude_inconsistencies(mr_results)
    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 pvalue (paired, just uses data with pair)
    # ================================================================== #
    paired_recordings = []
    for rec in mr_results.recording.unique():
        recframe = mr_results[mr_results['recording'] == rec]
        if not len(recframe['m'].tolist()) == 2:
            continue
        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)
    resultfile = '{}Fig1_boxplot_onlypaired.pdf'.format(resultpath)
    plt.savefig(resultfile, bbox_inches='tight')
    plt.show()
    plt.close()


def boxplot_focal_contra_Engel(pkl_file_pre, pkl_file_inter, resultpath):

    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_utilities.exclude_inconsistencies(mr_results)

    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 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)
    resultfile = '{}Fig1_boxplot_1aEngel.pdf'.format(resultpath)
    plt.savefig(resultfile, bbox_inches='tight')
    plt.show()
    plt.close()


def boxplot_focal_contra(pkl_file_pre, pkl_file_inter, resultpath):

    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_utilities.exclude_inconsistencies(mr_results)
    mr_results['epsilon'] = mr_results.apply(lambda row: 1-row['m'], axis=1)
    mr_results['logepsilon'] = mr_results.apply(lambda row: np.log(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)
    resultfile = '{}Fig1_boxplot.pdf'.format(resultpath)
    plt.savefig(resultfile, bbox_inches='tight')
    plt.show()
    plt.close()


def patientwise_pointplots(pkl_file_pre, pkl_file_inter, result_path):
    # -----------------------------------------------------------------#
    # 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_utilities.exclude_inconsistencies(mr_results)

    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 the number of recordings is sufficiently large, output p-value
        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)
    plotfile_mx = '{}S2_patients.pdf'.format(result_path)
    fig.savefig(plotfile_mx, bbox_inches='tight')
    plt.show()


def boxplot_brainregions_focus_onlypaired(pkl_file_pre, pkl_file_inter,
                                          result_path):

    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_utilities.exclude_inconsistencies(mr_results)

    # -----------------------------------------------------------------#
    # 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['epsilon'] = mr_results.apply(lambda row: 1-row['m'], axis=1)
    mr_results['logepsilon'] = mr_results.apply(
        lambda row: np.log(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)
    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
            mr_results_paired = mr_results_paired.append(area_frame,
                                                         ignore_index=True)
    mr_results = mr_results_paired

    # -----------------------------------------------------------------#
    # Compute p-values
    # -----------------------------------------------------------------#
    global_pvals, pvals = hierarchical_model.hierarchical_model_regions(
        mr_results)

    # -----------------------------------------------------------------#
    # 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)
    plotfile_mx = '{}Fig1_regions_onlypaired.pdf'.format(result_path)
    fig.savefig(plotfile_mx, bbox_inches='tight')
    plt.show()
    plt.close('all')


def boxplot_brainregions_focus(pkl_file_pre, pkl_file_inter, result_path):

    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_utilities.exclude_inconsistencies(mr_results)

    # -----------------------------------------------------------------#
    # 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['epsilon'] = mr_results.apply(lambda row: 1-row['m'], axis=1)
    mr_results['logepsilon'] = mr_results.apply(
        lambda row: np.log(1 - row['m']), axis=1)
    quantity = 'epsilon'
    # -----------------------------------------------------------------#
    # Compute p-values
    # -----------------------------------------------------------------#
    global_pvals, pvals = hierarchical_model.hierarchical_model_regions(
                                    mr_results)
    # -----------------------------------------------------------------#
    # 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)
    plotfile_mx = '{}Fig1_regions.pdf'.format(result_path)
    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):

    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_utilities.exclude_inconsistencies(mr_results)

    # -----------------------------------------------------------------#
    # 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['epsilon'] = mr_results.apply(lambda row: 1-row['m'], axis=1)
    mr_results['logepsilon'] = mr_results.apply(
        lambda row: np.log(1 - row['m']), axis=1)
    quantity = 'epsilon'

    # -----------------------------------------------------------------#
    # Compute p-values
    # -----------------------------------------------------------------#
    global_pvals, pvals = hierarchical_model.hierarchical_model_regions(
        mr_results)

    # -----------------------------------------------------------------#
    # 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)

    plotfile_mx = '{}Fig1_regions_Engel1a.pdf'.format(result_path)
    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'
    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)
    boxplot_focal_contra_Engel(pkl_file_pre, pkl_file_inter,
                         result_path)
    boxplot_focal_contra_onlypaired(pkl_file_pre, pkl_file_inter,
                         result_path)


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)


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'

    boxplot_brainregions_focus(pkl_file, pkl_file_inter, result_path)
    boxplot_brainregions_focus_Engel(pkl_file, pkl_file_inter, result_path)
    boxplot_brainregions_focus_onlypaired(pkl_file, pkl_file_inter,
                                          result_path)


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


if __name__ == '__main__':
    main()