gin commit from UX410UAK

New files: 10

python_code/analysis_PlosCB/create_figure_1.py

@@ -0,0 +1,943 @@
+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()

python_code/analysis_PlosCB/create_figure_2.py

@@ -0,0 +1,920 @@
+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()

python_code/analysis_PlosCB/create_supp_figures.py

@@ -0,0 +1,492 @@
+import os
+import pandas as pd
+import mrestimator as mre
+import numpy as np
+import matplotlib.pylab as plt
+import seaborn as sns
+import matplotlib.gridspec as gridspec
+from scipy import stats as sps
+import mr_utilities
+
+
+def cm2inch(value):
+    return value/2.54
+
+
+mycol = ['darkgreen', 'mediumaquamarine',  'slategray', 'darkblue','steelblue',
+         'firebrick', 'purple', 'orange', "red", "violet", "darkred",
+         "darkviolet", "green", "salmon", "black", 'darkgreen',
+         'mediumaquamarine',  'silver', 'darkblue','steelblue', 'firebrick',
+         'purple', 'orange', "red", "violet", "darkred", "darkviolet"]
+
+
+def distribution_A_recordings(resultpath):
+    data_path = '../Data/preictal/'
+    recordings = [dI for dI in os.listdir(data_path) if os.path.isdir(
+        os.path.join(data_path, dI))]
+
+    A_mean = []
+    A_nonzero = []
+    A_nonzero_frac = []
+
+    for recording in recordings:
+
+        outputfile = '{}{}/activity_SU_4ms.pkl'.format(data_path, recording)
+
+        if not os.path.isfile(outputfile):
+            print("No binning data for recording ", recording)
+            continue
+
+        data = pd.read_pickle(outputfile)
+        binnings = data['binning'].tolist()
+
+        for binning in binnings:
+            activity = data[data['binning'] == binning]['activity'].tolist()[0]
+            A_mean.append(np.mean(activity))
+            A_nonzero.append(mr_utilities.count_nonzero(activity))
+            A_nonzero_frac.append(mr_utilities.count_nonzero(activity) /
+                                  len(activity) * 100)
+
+    A_mean_Hz = [A_m*(1000/4) for A_m in A_mean]
+    A_mean = A_mean_Hz
+    quantile = 0.50
+    A_mean_50 = np.quantile(A_mean, quantile)
+    A_nonzero_50 = np.quantile(A_nonzero, quantile)
+    A_nonzero_frac_50 = np.quantile(A_nonzero_frac, quantile)
+
+    quantile = 0.25
+    A_mean_25 = np.quantile(A_mean, quantile)
+    A_nonzero_25 = np.quantile(A_nonzero, quantile)
+    A_nonzero_frac_25 = np.quantile(A_nonzero_frac, quantile)
+
+    sns.set(style='white')
+    fig, ax = plt.subplots(1, 3, figsize=((cm2inch(21), cm2inch(2.8))))
+    bins = np.arange(0, np.max(A_mean), 20)
+    ax[0].hist(A_mean, alpha=0.6, color='darkblue')
+    ax[0].set_xlabel(r'population firing rate R (Hz)')
+    ax[0].set_ylabel(r'# recordings')
+    ax[0].text(0.8, 0.8, r'$q_{50}$'+'={:.1f} Hz'.format(A_mean_50),
+               horizontalalignment='center', verticalalignment='center',
+               transform=ax[0].transAxes, color='green', fontsize=11)
+    ax[0].text(0.8, 0.65, r'$q_{25}$'+'={:.1f} Hz'.format(A_mean_25),
+               horizontalalignment='center', verticalalignment='center',
+               transform=ax[0].transAxes, color='green', fontsize=11)
+
+    bins = np.arange(0, np.max(A_nonzero), 20)
+    ax[1].hist(A_nonzero, alpha=0.6, color='darkblue')
+    ax[1].set_xlabel(r'$n_{A_t\neq 0}$')
+    ax[1].text(0.8, 0.8, r'$q_{50}$'+'={:.0f}'.format(A_nonzero_50),
+               horizontalalignment='center', verticalalignment='center',
+               transform=ax[1].transAxes, color='green', fontsize=11)
+    ax[1].text(0.8, 0.65, r'$q_{25}$'+'={:.0f}'.format(A_nonzero_25),
+               horizontalalignment='center', verticalalignment='center',
+               transform=ax[1].transAxes, color='green', fontsize=11)
+
+    bins = np.arange(0, np.max(A_nonzero_frac), 20)
+    ax[2].hist(A_nonzero_frac, alpha=0.6, color='darkblue')
+    ax[2].set_xlabel(r'$n_{A_t\neq 0}$ / $n_{A_t}$ (%)')
+    ax[2].text(0.7, 0.8, r'$q_{50}$'+'={:.1f} %'.format(A_nonzero_frac_50),
+               horizontalalignment='center',
+               verticalalignment='center', transform=ax[2].transAxes,
+               color='green', fontsize=11)
+    ax[2].text(0.7, 0.65, r'$q_{25}$'+'={:.1f} %'.format(A_nonzero_frac_25),
+               horizontalalignment='center',
+               verticalalignment='center', transform=ax[2].transAxes,
+               color='green', fontsize=11)
+
+    ax[0].tick_params(top='off', bottom='on', left='on', right='off',
+                      labelleft='on',labelbottom='on')
+    ax[1].tick_params(top='off', bottom='on', left='on', right='off',
+                      labelleft='on', labelbottom='on')
+    ax[2].tick_params(top='off', bottom='on', left='on', right='off',
+                      labelleft='on', labelbottom='on')
+
+    sns.despine()
+    plt.subplots_adjust(bottom=None, right=None,
+                        wspace=0.6, hspace=0.3, left=0.15, top=1.2)
+    figfile = '{}S4_recording_statistics.pdf'.format(resultpath)
+    plt.savefig(figfile, bbox_inches='tight')
+    plt.show()
+
+
+def mt_all_supp(save_path,fit_method, kmin, kmax, windowsize, windowstep, dt):
+
+    patients = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
+                18, 19, 20]
+    fig = plt.figure(figsize=(cm2inch(2*13.5), cm2inch(1.2*10*5)))
+    outer = gridspec.GridSpec(10, 2, wspace=0.2, hspace=0.9)
+    laterality = ["SOZ", "nSOZ"]
+
+    for ip, patient in enumerate(patients):
+        inner = gridspec.GridSpecFromSubplotSpec(2, 1, subplot_spec=outer[ip],
+                                                 wspace=0.1, hspace=0.9)
+
+        for j, soz in enumerate(laterality):
+            ax = plt.Subplot(fig, inner[j])
+            ax = mt_plot(ax, patient, soz, save_path,fit_method, kmin, kmax,
+                          windowsize, windowstep, dt,
+                          excludeInconsistencies=True)
+
+            if j == 0:
+                ax.set_title('Patient {}'.format(int(patient)), fontsize=12)
+                ax.tick_params(top='off', bottom='on', left='off',
+                               right='off', labelleft='on',
+                               labelbottom='off')
+                ax.set_xlabel('')
+
+            fig.add_subplot(ax)
+    sns.despine(fig)
+    plt.savefig(
+        '{}/S6_allmt_{}_{}_winsize{}_winstep{}_kmin{}_kmax{}.pdf'.format(
+            save_path, patient, fit_method,
+            windowsize, windowstep, kmin, kmax), bbox_inches='tight')
+
+    plt.savefig(
+        '{}/S6_allmt_{}_{}_winsize{}_winstep{}_kmin{}_kmax{}.png'.format(
+            save_path, patient, fit_method,
+            windowsize, windowstep, kmin, kmax), bbox_inches='tight', dpi=350)
+    plt.savefig(
+        '{}/S6_allmt_{}_{}_winsize{}_winstep{}_kmin{}_kmax{}.svg'.format(
+            save_path, patient, fit_method,
+            windowsize, windowstep, kmin, kmax), bbox_inches='tight')
+    plt.show()
+    plt.close()
+
+
+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 mt_plot(ax, patient, soz, save_path, fit_method, kmin, kmax,
+            windowsize, windowstep, dt, excludeInconsistencies=True):
+
+    recordings = mr_utilities.list_preictrecordings(patient)[0]
+    binning = binning_label(patient, soz)
+
+    for irec, rec in enumerate(recordings):
+
+        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:
+            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:
+                continue
+
+        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()
+        mt = mt_frame.m.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.plot(t_middle, mt, '.', label='Rec {}'.format(irec + 1),
+                color=mycol[irec], markersize=0.5)
+        ax.axhline(y=1, color='black', linestyle='--', alpha=0.3)
+        ax.set_ylabel(r"$\hat{m}$")
+        ax.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.set_title("{}".format(focus_label), loc='left')
+        ax.set_xlim((-10, 610))
+
+    handles, labels = ax.get_legend_handles_labels()
+    if len(labels) == 0:
+        if soz == 'SOZ':
+            ax.set_title("ipsi", loc='left')
+        if soz == 'nSOZ':
+            ax.set_title("contra", loc='left')
+        ax.axhline(y=1, color='black', linestyle='--', alpha=0.3)
+        ax.set_ylabel(r"$\hat{m}$")
+        ax.set_ylim((0.8, 1.05))
+        ax.set_xlim((-10, 610))
+
+    ax.xaxis.grid(False)
+    ax.set_xlabel("time (s)")
+
+    ax.tick_params(top='off', bottom='off', left='off', right='off',
+                   labelleft='on', labelbottom='off')
+    ax.tick_params(top='off', bottom='on', left='off', right='off',
+                   labelleft='on', labelbottom='on')
+    sns.despine()
+    return ax
+
+
+def windowsize_analysis_recs(data_path, recording, binning, outputfilename,
+                             resultpath):
+
+    outputfile = '{}{}/{}'.format(data_path, recording, outputfilename)
+
+    if not os.path.isfile(outputfile):
+        print("No binning data for recording ", recording)
+        print(outputfile)
+        return
+    data = pd.read_pickle(outputfile)
+
+    dt = data['dt'].unique()[0]
+    activity = data[data['binning'] == binning]['activity'].tolist()[0]
+
+    winstep = 100
+    ksteps = (1, 400)
+    windowsize = np.arange(2500, 40000, 2500)
+    N_estimates = len(activity) - np.max(windowsize)
+    tau_estimate = []
+    tau_conf = []
+    m_estimate = []
+    m_conf = []
+    for Lw in windowsize:
+        tau_L = []
+        m_L = []
+        for iw in range(0, N_estimates, winstep):
+            activity_window = activity[iw:iw+Lw]
+            input = mre.input_handler(activity_window)
+            rk = mre.coefficients(input, steps=ksteps, dt=dt)
+            fitres_offset = mre.fit(rk, fitfunc=mre.f_exponential_offset)
+            tau_L.append(fitres_offset.tau)
+            m_L.append(fitres_offset.mre)
+        conf_int_off = sps.mstats.mquantiles(tau_L, prob=[0.025, 0.975],
+                                             alphap=0, betap=1, axis=None,
+                                             limit=())
+        conf_int_off = [abs(np.mean(tau_L) - ci) for ci in conf_int_off]
+        tau_conf.append(conf_int_off)
+        tau_estimate.append(np.mean(tau_L))
+
+        conf_int_off = sps.mstats.mquantiles(m_L, prob=[0.025, 0.975],
+                                             alphap=0, betap=1, axis=None,
+                                             limit=())
+        conf_int_off = [abs(np.mean(m_L) - ci) for ci in conf_int_off]
+        m_conf.append(conf_int_off)
+        m_estimate.append(np.mean(m_L))
+
+    windowsize_seconds = [w * 4 / 1000 for w in windowsize]
+    fig, ax = plt.subplots(1, 1, figsize=((cm2inch(9), cm2inch(0.9 * 4))))
+    ax.scatter(windowsize_seconds, tau_estimate,
+               label=r'$\hat{\tau}_{offset}$', c='steelblue',
+               marker='.', s=9)
+    ax.errorbar(windowsize_seconds, tau_estimate,
+                yerr=np.transpose(np.array(tau_conf)),
+                linestyle="None", c='steelblue', elinewidth=0.9,
+                capsize=2.5)
+
+    ax.set_ylim((-10, 800))
+    T = [20, 40, 60, 80, 100, 120, 140]
+    ax.set_xticks(T)
+    ax.set_xlabel('window size (s)')
+    ax.set_ylabel(r'$\hat{\tau}$ (ms)')
+    sns.despine()
+    plt.subplots_adjust(bottom=None, right=None,
+                        wspace=0.4, hspace=0.3, left=0.15, top=1.2)
+    figfile = '{}S5_estimates_winsizes_rec_{}_{}.pdf'.format(
+        resultpath, recording, binning)
+    plt.savefig(figfile, bbox_inches='tight')
+
+    fig, ax = plt.subplots(1, 1, figsize=((cm2inch(9), cm2inch(0.9 * 4))))
+    ax.scatter(windowsize_seconds, m_estimate,
+               label=r'$\hat{\tau}_{offset}$', c='royalblue',
+               marker='.', s=9)
+    ax.errorbar(windowsize_seconds, m_estimate,
+                yerr=np.transpose(np.array(m_conf)),
+                linestyle="None", c='royalblue', elinewidth=0.9,
+                capsize=2.5)
+
+    T = [20, 40, 60, 80, 100, 120, 140]
+    ax.set_xticks(T)
+    ax.set_ylim((0.5, 1.25))
+    ax.set_xlabel('window size (s)')
+    ax.set_ylabel(r'$\hat{m}$')
+    sns.despine()
+    plt.subplots_adjust(bottom=None, right=None,
+                        wspace=0.4, hspace=0.3, left=0.15, top=1.2)
+    figfile = '{}S5_estimates_m_winsizes_rec_{}_{}.pdf'.format(
+        resultpath, recording, binning)
+    plt.savefig(figfile, bbox_inches='tight')
+    plt.show()
+
+
+def windowsize_analysis():
+    binning = 'left'
+    outputfilename = 'activity_SU_4ms.pkl'
+    resultpath = '../Results/preictal/singleunits/'
+
+    data_path = '../Data/interictal/'
+    recording = '13ref'
+    windowsize_analysis_recs(data_path, recording, binning, outputfilename,
+                             resultpath)
+
+    data_path = '../Data/preictal/'
+    recording = '13_02'
+    windowsize_analysis_recs(data_path, recording, binning, outputfilename,
+                             resultpath)
+
+
+def all_mt():
+    result_path = '../Results/preictal/singleunits/'
+    ksteps = (1, 400)
+    kmin = ksteps[0]
+    kmax=ksteps[1]
+    windowsize = 20000
+    windowstep = 500
+    dt = 4
+    mt_all_supp(result_path, 'offset', kmin, kmax, windowsize, windowstep, dt)
+
+
+def example_autocorrelations(patient, pkl_file_inter, pkl_file_pre,
+                             resultpath, **kwargs):
+    mr_results_pre = pd.read_pickle(pkl_file_pre)
+    mr_results_inter = pd.read_pickle(pkl_file_inter)
+    mr_results_pre['rectype'] = 'preictal'
+    mr_results_inter['rectype'] = 'interictal'
+    mr_results = pd.concat([mr_results_pre, mr_results_inter],
+                           ignore_index=True)
+    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)]
+
+    mr_results = mr_results[mr_results['patient'] == patient]
+    binning_foci = ['focal', 'contra-lateral']
+    kmin = mr_results.kmin.unique()[0]
+    kmax = mr_results.kmax.unique()[0]
+
+    sns.set(style="white")
+    fig, ax = plt.subplots(2, 4, figsize=(cm2inch(24), cm2inch(6)))
+    fit_method = 'offset'
+    mr_rec_pre = mr_results[mr_results['rectype'] == 'preictal']
+    pre_recordings_choice = mr_rec_pre.recording.unique()[:3]
+
+    for si, soz in enumerate(binning_foci):
+        # plot interictal in first position
+        mr_rec_inter = mr_results[mr_results['rectype'] == 'interictal']
+        mr_rec_inter = mr_rec_inter[mr_rec_inter['binning_focus'] == soz]
+        rk = mr_rec_inter['rk'].item()[0]
+        k_steps = mr_rec_inter['rk_steps'].item()
+        dt = mr_rec_inter['dt'].item()
+        time_steps = np.array(k_steps[0]) * dt
+        ax[si, 0].plot(time_steps, rk, '-', linewidth=0.8,
+                      c='slategray', alpha=0.4)
+        fitfunc = mr_utilities.string_to_func(fit_method)
+        fitparams = mr_rec_inter['fit_params'].item()
+        fitparams = tuple(fitparams[0])
+        m = mr_rec_inter['m'].item()
+        ax[si, 0].plot(time_steps, fitfunc(time_steps, *fitparams),
+                       label=r'$m$ = {:.2f}'.format(m), color='slategray')
+
+        ax[si, 0].tick_params(top='off', bottom='on', left='on',
+                              right='off', labelleft='on', labelbottom='on')
+
+        ax[si, 0].set_xlabel(r"time lag $k$ (ms)")
+        ax[si, 0].legend(loc=1, fontsize='x-small')
+        ax[si, 0].ticklabel_format(axis='y', style='sci', scilimits=(-1, 1),
+                                   useMathText=True)
+
+        # plot 3 preictal recordings next to it
+        mr_rec_pre = mr_results[mr_results['rectype'] == 'preictal']
+        mr_rec_pre = mr_rec_pre[mr_rec_pre['binning_focus'] == soz]
+        for irec, rec in enumerate(pre_recordings_choice):
+            mr_rec = mr_rec_pre[mr_rec_pre['recording'] == rec]
+            rk = mr_rec['rk'].item()[0]
+            k_steps = mr_rec['rk_steps'].item()
+            dt = mr_rec['dt'].item()
+            time_steps = np.array(k_steps[0]) * dt
+            ax[si, irec+1].plot(time_steps, rk, '-', linewidth=0.8,
+                           c='darkblue', alpha=0.4)
+            fitfunc = mr_utilities.string_to_func(fit_method)
+            fitparams = mr_rec['fit_params'].item()
+            fitparams = tuple(fitparams[0])
+            m = mr_rec['m'].item()
+            ax[si, irec+1].plot(time_steps, fitfunc(time_steps, *fitparams),
+                           label=r'$m$ = {:.2f}'.format(m), color='darkblue')
+
+            ax[si, irec+1].tick_params(top='off', bottom='on', left='on',
+                                  right='off', labelleft='on',labelbottom='on')
+
+            ax[si, irec+1].set_xlabel(r"time lag $k$ (ms)")
+            ax[si, irec+1].legend(loc=1, fontsize='x-small')
+            ax[si, irec+1].ticklabel_format(axis='y', style='sci',
+                                       scilimits=(-1, 1), useMathText=True)
+    ax[0, 0].set_ylabel("ACF")
+    ax[1, 0].set_ylabel("ACF")
+    plt.subplots_adjust(bottom=None, right=None,
+                        wspace=0.5, hspace=1, left=0.15, top=1.2)
+    sns.despine(fig)
+    resultpath_rec = '{}S1_autocorrelations_kmin{}_kmax{}_patient{}.pdf'.format(resultpath, kmin, kmax, patient)
+    plt.savefig(resultpath_rec, bbox_inches='tight')
+    resultpath_rec = '{}S1_autocorrelations_kmin{}_kmax{}_patient{}.png'.format(resultpath, kmin, kmax, patient)
+    plt.savefig(resultpath_rec, bbox_inches='tight', dpi=300)
+    resultpath_rec = '{}S1_autocorrelations_kmin{}_kmax{}_patient{}.svg'.format(resultpath, kmin, kmax, patient)
+    plt.savefig(resultpath_rec, bbox_inches='tight')
+    plt.close()
+
+
+def create_example_ACF_plot():
+    pkl_file_pre = '../Results/preictal/singleunits/mr_results.pkl'
+    pkl_file_inter = '../Results/interictal/singleunits/mr_results.pkl'
+    resultpath = '../Results/preictal/singleunits/'
+    ksteps = (1, 400)
+    patient = 12
+    example_autocorrelations(patient, pkl_file_inter, pkl_file_pre,
+         resultpath, kmin=ksteps[0], kmax=ksteps[1])
+    patient = 13
+    example_autocorrelations(patient, pkl_file_inter, pkl_file_pre,
+         resultpath, kmin=ksteps[0], kmax=ksteps[1])
+
+
+def main():
+    #create_example_ACF_plot()
+    #windowsize_analysis()
+    all_mt()
+
+
+if __name__ == '__main__':
+    main()

python_code/analysis_PlosCB/hierarchical_model.py

@@ -0,0 +1,39 @@
+from pymer4.models import Lmer
+
+
+def hierarchical_model_regions(dataframe):
+    # overall p-value
+    df_reduced = dataframe[['logepsilon', 'patient', 'binning_focus', 'region']]
+    model = Lmer('logepsilon ~ binning_focus*region + (1|patient)',
+                     data=df_reduced)
+
+    model.fit()
+    anova_result = model.anova(force_orthogonal=True)
+    global_pvals = anova_result['P-val']
+    # separate p-values
+    pvals_regions = dict()
+    for reg in dataframe.region.unique():
+        print(reg)
+        df_region = dataframe[dataframe['region'] == reg]
+        df_reduced = df_region[['logepsilon', 'patient', 'binning_focus']]
+        # for formulas, check
+        # https://stats.stackexchange.com/questions/13166/rs-lmer-cheat-sheet
+        model = Lmer('logepsilon ~ binning_focus + (1|patient)',
+                     data=df_reduced)
+        model.fit()
+        anova_result = model.anova(force_orthogonal=True)
+        region_pval = anova_result['P-val'].item()
+        pvals_regions[reg] = region_pval
+
+    return global_pvals, pvals_regions
+
+
+def hierarchical_model_hemispheres(dataframe):
+    df_reduced = dataframe[['logepsilon', 'patient', 'binning_focus']]
+    model = Lmer('logepsilon ~ binning_focus + (1|patient)',
+                 data=df_reduced)
+    model.fit()
+    anova_result = model.anova(force_orthogonal=True)
+    pval = anova_result['P-val'].item()
+    return pval
+

python_code/analysis_PlosCB/mr_epilepsy_fullrec.py

@@ -0,0 +1,124 @@
+import os
+import pandas as pd
+import mrestimator as mre
+import numpy as np
+import math
+import mr_utilities
+
+
+def mr_analysis_confint(data_path, result_path, outputfilename,
+                        recordings='all', binnings=[], ksteps=(1, 400),
+                        fit_method='offset'):
+
+    if recordings == 'all':
+        recordinglist = [dI for dI in os.listdir(data_path) if
+                         os.path.isdir(os.path.join(data_path, dI))]
+    else:
+        recordinglist = recordings
+
+    focus_list = '../Data/patients.txt'
+    f = open(focus_list, 'r')
+    foci = [line.rstrip('\n').split('\t') for line in f.readlines()]
+
+    # -----------------------------------------------------------------------#
+    # Run through all recordings and binnings and perform MR estimation
+    # -----------------------------------------------------------------------#
+    for recording in recordinglist:
+        outputfile = '{}{}/{}'.format(data_path, recording, outputfilename)
+
+        if not os.path.isfile(outputfile):
+            print("No binning data for recording ", recording)
+            print(outputfile)
+            continue
+
+        data = pd.read_pickle(outputfile)
+        patient = int(data['patient'].unique()[0])
+
+        # read the location of the focus
+        focus = 'NA'
+        for i, idf in enumerate(foci[1:]):
+            if int(idf[0]) == int(patient):
+                focus = idf[1]
+        if focus == 'NA':
+            print('Error: patient {} not in foci list.'.format(patient))
+            continue
+
+        existing_binnings = data['binning'].tolist()
+        if not binnings:
+            binnings = existing_binnings
+
+        for binning in binnings:
+
+            binning_focus = mr_utilities.binning_focus(binning, patient)
+            dt = data['dt'].unique()[0]
+            activity = data[data['binning'] == binning]['activity'].tolist()[0]
+
+            if mr_utilities.count_nonzero(activity) < 1000:
+                continue
+
+            # split into trials for bootstrap confidence intervals
+            ntrials = 30
+            len_trial = int(np.floor(len(activity) / ntrials))
+            unusable_activity = len(activity) - ntrials * len_trial
+            activity = activity[unusable_activity:]
+            trials = np.reshape(activity, (ntrials, len_trial))
+
+            # apply MR estimator
+            prepared = mre.input_handler(trials)
+            rk = mre.coefficients(prepared, dt=dt, dtunit='ms',
+                                  steps=ksteps, method='sm', knownmean=None)
+            try:
+                maxfev = 200 * (len(rk.coefficients)+1)
+                fit_res = mre.fit(rk, fitfunc=mr_utilities.string_to_func(
+                                    fit_method), maxfev=maxfev)
+
+                if not math.isnan(fit_res.tau):
+                    inconsistencies = mr_utilities.consistency_checks(
+                                            activity, rk, fit_res)
+                    mr_savefile(result_path, patient, recording,
+                                focus, binning, binning_focus, ksteps,
+                                fit_method, rk, fit_res, inconsistencies)
+
+            except RuntimeError:
+                print('RUNTIME-ERROR', recording, binning)
+                pass
+
+
+def mr_savefile(result_path, patient, recording, focus, binning,
+                binning_focus, ksteps, fit_method, rk, fitres,
+                inconsistencies):
+
+    datafile_patient = result_path + 'mr_results.pkl'
+    m_quantiles = fitres.mrequantiles
+    tau_quantiles = fitres.tauquantiles
+    rk_stderr = rk.stderrs
+
+    save_res = dict()
+    save_res['patient'] = patient
+    save_res['recording'] = recording
+    save_res['focus'] = focus
+    save_res['binning_focus'] = binning_focus
+    save_res['binning'] = binning
+    save_res['kmin'] = ksteps[0]
+    save_res['kmax'] = ksteps[1]
+    save_res['rk'] = [list(rk.coefficients)]
+    save_res['rk_steps'] = [list(rk.steps)]
+    save_res['dt'] = rk.dt
+    save_res['fit_method'] = fit_method
+    save_res['fit_params'] = [list(fitres.popt)]
+    save_res['m'] = fitres.mre
+    save_res['ssres'] = fitres.ssres
+    save_res['ssres_weighted'] = fitres.ssres*len(fitres.popt)
+    save_res['tau'] = fitres.tau
+    save_res['inconsistencies'] = inconsistencies
+    save_res['rk_stderr'] = rk_stderr
+    save_res['m_quantiles'] = m_quantiles
+    save_res['tau_quantiles'] = tau_quantiles
+
+    if os.path.isfile(datafile_patient):
+        mr_results_patient = pd.read_pickle(datafile_patient)
+        mr_results_patient = mr_results_patient.append(save_res,
+                                                       ignore_index=True)
+    else:
+        mr_results_patient = pd.DataFrame.from_records([save_res])
+    mr_results_patient.to_pickle(datafile_patient)

python_code/analysis_PlosCB/mr_epilepsy_timeresolved.py

@@ -0,0 +1,225 @@
+import numpy as np
+import pandas as pd
+import mrestimator as mre
+import os
+import time
+import mr_utilities
+import math
+
+
+def analyse_timeresolved_mr(data_path, result_path, recordings,
+                            binnings='all', ksteps=(1,400),
+                            windowsize=20000, windowstep=500,
+                            fit_method="offset"):
+
+    for recording in recordings:
+        outputfile = '{}{}/activity_SU_4ms.pkl'.format(data_path, recording)
+
+        if os.path.isfile(outputfile):
+            data = pd.read_pickle(outputfile)
+        else:
+            print('Error: no datafile for recording {}'.format(recording))
+            continue
+
+        existing_binnings = data['binning'].unique().tolist()
+        patient = int(data['patient'].unique()[0])
+
+        if binnings == 'all':
+            binnings = existing_binnings
+
+        for binning in binnings:
+
+            dt = data['dt'].unique()[0]
+            activity_list = data[data['binning'] == binning]['activity'].tolist()
+
+            if len(activity_list) == 0:
+                continue
+            if len(activity_list) > 1:
+                print('Error: Multiple activity entries for binning.')
+                data = data[~data['activity'].apply(tuple).duplicated()]
+                activity_list = data[data['binning'] == binning][
+                    'activity'].tolist()
+                activity = activity_list[0]
+            else:
+                activity = activity_list[0]
+
+            t_onset = 600 # in seconds, given 10min recordings
+            # the recording is significantly less than 10min long,
+            # it is unclear whether seizure onset is indeed at the end of
+            # the recording, therefore better skip
+            if len(activity) < 0.8 * t_onset * 10**3 / dt:
+                continue
+
+            save_path_recording = '{}{}/'.format(result_path, recording)
+
+            timeresolved_mr(activity, windowsize, windowstep,
+                            ksteps[0], ksteps[1], fit_method,
+                            save_path_recording, patient, recording,
+                            binning, dt)
+            print('break to cool cpu...')
+            time.sleep(30)
+
+
+def timeresolved_mr(activity, windowsize, windowstep, kmin, kmax,
+                    fit_method, save_path_recording, patient, recording,
+                    binning, dt):
+
+    """ performs MR estimation in sliding windows of provided activity and
+    calls function to store results in pkl file
+    """
+
+    for t in range(0, len(activity) - windowsize, windowstep):
+        windowed_data = activity[t:t + windowsize]
+        prepared = mre.input_handler(windowed_data)
+        rk = mre.coefficients(prepared, dt=dt, dtunit='ms', steps=(kmin, kmax))
+        try:
+            fit_res = mre.fit(rk, fitfunc=mr_utilities.string_to_func(
+                fit_method))
+            t0 = t * dt
+
+            mt_savefile(save_path_recording, patient, recording, binning,
+                        dt, kmin, kmax, fit_method, windowsize, windowstep,
+                        t0=t0, activity=windowed_data, rk=rk, fitres=fit_res)
+        except (RuntimeError, ValueError):
+            continue
+
+
+def mt_savefile(result_path, patient, recording, binning, dt, kmin, kmax,
+                fit_method, windowsize, windowstep, t0, activity, rk,
+                fitres):
+
+    """ save results of m(t) estimation performed by function
+    timeresolved_mr to pkl file.
+    """
+
+    datafile = '{}mt_results_{}_{}_kmin{}_kmax{}_winsize{}_winstep{}_dt{}' \
+               '.pkl'.format(result_path, binning, fit_method, kmin, kmax,
+                             windowsize, windowstep, dt)
+
+    if math.isnan(fitres.tau):
+        inconsistencies = ['not_converged']
+    else:
+        inconsistencies = mr_utilities.consistency_checks(activity, rk, fitres)
+    try:
+        fano = np.var(activity)/np.mean(activity)
+        variance = np.var(activity)
+    except ZeroDivisionError:
+        fano = np.nan
+        variance = np.nan
+
+    save_res = dict()
+    save_res['patient'] = patient
+    save_res['recording'] = recording
+    save_res['binning'] = binning
+    save_res['kmin'] = kmin
+    save_res['kmax'] = kmax
+    save_res['windowsize'] = windowsize
+    save_res['windowstep'] = windowstep
+    save_res['t0'] = t0
+    save_res['rk'] = [list(rk.coefficients)]
+    save_res['rk_steps'] = [list(rk.steps)]
+    save_res['dt'] = dt
+    save_res['fit_method'] = fit_method
+    save_res['fit_params'] = [list(fitres.popt)]
+    save_res['m'] = fitres.mre
+    save_res['tau'] = fitres.tau
+    save_res['inconsistencies'] = inconsistencies
+    save_res['fano'] = fano
+    save_res['variance'] = variance
+    save_res['offset'] = fitres.popt[2]
+
+    if os.path.isfile(datafile):
+        mr_results_patient = pd.read_pickle(datafile)
+        mr_results_patient = mr_results_patient.append(save_res,
+                                                       ignore_index=True)
+    else:
+        mr_results_patient = pd.DataFrame.from_records([save_res])
+    mr_results_patient.to_pickle(datafile)
+
+
+def mt_x_parts(x, outputfilename, data_path, result_path, recordings,
+               binnings='all', ksteps=(1, 400), fit_method ="offset"):
+
+    mx_allresults = None
+    for recording in recordings:
+        outputfile = '{}{}/{}'.format(data_path, recording, outputfilename)
+
+        if os.path.isfile(outputfile):
+            data = pd.read_pickle(outputfile)
+        else:
+            print(' File not found {}'.format(outputfile))
+            continue
+
+        existing_binnings = data['binning'].unique().tolist()
+        patient = int(data['patient'].unique()[0])
+
+        kmin = ksteps[0]
+        kmax = ksteps[1]
+
+        if binnings == 'all':
+            binnings = existing_binnings
+
+        for binning in binnings:
+            dt = data['dt'].unique()[0]
+            activity_list = data[data['binning'] == binning]['activity'].tolist()
+            activity = activity_list[0]
+            partlength = int(np.floor(len(activity)/x))
+            tau = []
+            m = []
+            binning_focus = mr_utilities.binning_focus(binning, patient)
+
+            for i in range(x):
+                activity_part = activity[i*partlength:(i+1)*partlength]
+                nspikes = np.sum(np.array(activity_part))
+                n_nonzero = mr_utilities.count_nonzero(activity_part)
+                prepared = mre.input_handler(activity_part)
+                rk = mre.coefficients(prepared, dt=dt, dtunit='ms',
+                                      steps=(kmin, kmax))
+                try:
+                    fit_res = mre.fit(rk,
+                                      fitfunc=mr_utilities.string_to_func(
+                                          fit_method))
+
+                    if not math.isnan(fit_res.tau):
+                        tau.append(fit_res.tau)
+                        m.append(fit_res.mre)
+
+                        inconsistencies = mr_utilities.consistency_checks(
+                            activity_part, rk, fit_res)
+                        fano = np.var(activity_part)/np.mean(activity_part)
+                        variance = np.var(activity_part)
+                        mx_result = {'patient': patient,
+                                      'recording': recording,
+                                      'binning': binning,
+                                      'binning_focus': binning_focus,
+                                      'partno': i+1,
+                                      'm': fit_res.mre,
+                                      'tau': fit_res.tau,
+                                      'fano': fano,
+                                      'variance': variance,
+                                      'offset': fit_res.popt[2],
+                                      'rk': [list(rk.coefficients)],
+                                      'fit_params': [list(fit_res.popt)],
+                                      'nspikes': nspikes,
+                                      'n_nonzero': n_nonzero,
+                                      'inconsistencies': inconsistencies,
+                                      'kmin': kmin,
+                                      'kmax': kmax,
+                                      'fit_method': fit_method}
+
+                        if mx_allresults is None:
+                            mx_allresults = pd.DataFrame.from_records(
+                                [mx_result])
+
+                        else:
+                            mx_allresults = mx_allresults.append(
+                                mx_result, ignore_index=True)
+
+                except (RuntimeError, ValueError):
+                    continue
+
+    datafile_mx = '{}mt_{}_results_allrecords_kmin{}_kmax{}_{}.pkl'.format(
+        result_path, x, ksteps[0], ksteps[1], fit_method)
+
+    if mx_allresults is not None:
+        mx_allresults.to_pickle(datafile_mx)

python_code/analysis_PlosCB/mr_utilities.py

@@ -0,0 +1,237 @@
+import sys
+sys.path.insert(0, '../mrestimator/')
+import numpy as np
+import functools
+import mrestimator as mre
+from scipy import stats
+
+
+def consistency_checks(activity, rk, fitres):
+    """ Checks the consistency of MR estimates
+    We only use H_nonzero and H_linear as exclusion criteria.
+    """
+
+    tau_max = 2000  # unusually large tau
+    tau_min = 0
+    inconsistencies = []
+    if fitres.tau > tau_max or fitres.tau < tau_min:
+        inconsistencies.append('H_tau')
+    if count_nonzero(activity) < 1000:
+        inconsistencies.append('H_nonzero')
+    if check_linear(rk, fitres.rsquared):
+        inconsistencies.append('H_linear')
+    if not random_residuals(rk, fitres.fitfunc, *fitres.popt):
+        inconsistencies.append('H_residual')
+    if not random_residuals_alternative(rk, fitres.fitfunc, *fitres.popt):
+        inconsistencies.append('H_residual_alt')
+    if Ftest05(rk, fitres.rsquared, fitres.popt):
+        inconsistencies.append('H_F05')
+    if Ftest10(rk, fitres.rsquared, fitres.popt):
+        inconsistencies.append('H_F10')
+
+    return inconsistencies
+
+
+def Ftest05(rk, adjusted_rsquared, popt):
+    n_obs = len(rk.coefficients)
+    n_param = len(popt)
+    r_squared = 1.0- (1.0 - adjusted_rsquared) * (n_obs-1-n_param)/(n_obs-1)
+    F = (r_squared / (n_param-1)) / ((1-r_squared)/(n_obs-n_param))
+    df1 = n_param-1 #n_obs
+    df2 = n_obs-n_param
+    p_value = 1 - stats.f.cdf(F, df1, df2)
+    if p_value > 0.05:
+        # complex model is not significanty better in explaining the data
+        return True
+    else:
+        return False
+
+
+def Ftest10(rk, adjusted_rsquared, popt):
+    n_obs = len(rk.coefficients)
+    n_param = len(popt)
+    r_squared = 1.0- (1.0 - adjusted_rsquared) * (n_obs-1-n_param)/(n_obs-1)
+    F = (r_squared / (n_param-1)) / ((1-r_squared)/(n_obs-n_param))
+    df1 = n_param-1 #n_obs
+    df2 = n_obs-n_param
+    p_value = 1 - stats.f.cdf(F, df1, df2)
+    if p_value > 0.10:
+        # complex model is not significanty better in explaining the data
+        return True
+    else:
+        return False
+
+
+def check_linear(rk, rsquared):
+    """ Check if linear function better fits the rk than the fit that was
+    done (typically exponential offset). Done by comparing the rsquared of
+    both fits, which are weighted with the number of parameters i the
+    fitfunction.
+    :return: bool
+        True: linear fit better (Inconsistency!)
+        False: linear fit worse
+    """
+
+    try:
+        linear_fit = mre.fit(rk, fitfunc=mre.f_linear, description='linear')
+        if linear_fit.rsquared is not None and rsquared is not None:
+            if linear_fit.rsquared > rsquared:
+                return True
+            else:
+                return False
+        else:
+            return True
+    except RuntimeError:
+        return False
+
+
+def random_residuals(rk, fitfunc, *args):
+    """ Check if residuals (fit minus data) are consistent with noise. If
+    this is not the case, i.e. if there are systematic trends in the
+    residuals, the fit is inconsistent.
+    """
+
+    residuals = rk.coefficients - fitfunc(rk.steps*rk.dt, *args)
+    k_win = 20
+    for k_1 in range(len(rk.steps)-(k_win+1)):
+        sign = []
+        for k in range(k_1, k_1+k_win):
+            if residuals[k] > 1*np.median(abs(residuals)):
+                sign.append(1)
+            elif residuals[k] < -1*np.median(abs(residuals)):
+                sign.append(-1)
+            else:
+                sign.append(0)
+        if len(set(sign)) <= 1 and sign[0] != 0:
+            return False
+    return True
+
+
+def random_residuals_alternative(rk, fitfunc, *args):
+
+    """ Check if residuals (fit minus data) are consistent with noise. If
+    this is not the case, i.e. if there are systematic trends in the
+    residuals, the fit is inconsistent.
+    """
+    residuals = rk.coefficients - fitfunc(rk.steps*rk.dt, *args)
+    k_win = 10
+    k_range = 50
+    sign = []
+    for k in range(rk.steps[0], rk.steps[-k_win]):
+        res_mean = np.mean(residuals[k:k+k_win])
+        if res_mean > 1*np.median(abs(residuals)):
+            sign.append(1)
+        elif res_mean < -1*np.median(abs(residuals)):
+            sign.append(-1)
+        else:
+            sign.append(0)
+    for si in range(len(sign)-k_range):
+        if len(set(sign[si:si+k_range])) <= 1 and sign[si] != 0:
+            return False
+    return True
+
+
+def list_preictrecordings(patient):
+    listfile = '../Data/preIct_recordings_patients.txt'
+    list = np.loadtxt(listfile, dtype=str)
+    indices = np.where(list[:, 1] == str(patient))
+
+    recordings = []
+    for i in indices:
+        recordings.append(list[i, 0])
+    return recordings
+
+
+def count_nonzero(activity):
+    activity_nonzero = 0
+    for t in activity:
+        if t != 0:
+            activity_nonzero += 1
+    return activity_nonzero
+
+
+def conjunction(*conditions):
+    return functools.reduce(np.logical_and, conditions)
+
+
+def string_to_func(fitmethod):
+
+    if fitmethod == 'exp':
+        fitfunc = mre.f_exponential
+    elif fitmethod == 'offset':
+        fitfunc = mre.f_exponential_offset
+    elif fitmethod == 'complex':
+        fitfunc = mre.f_complex
+    else:
+        print(
+            'Error: fitmethod {} does not exist'.format(fitmethod))
+        exit()
+    return fitfunc
+
+
+def twosided_pvalue(value, diff_distribution):
+    diff_distribution = np.array(diff_distribution)
+    if value < 0:
+        p_value = len(diff_distribution[diff_distribution <=
+                                              value])/len(
+            diff_distribution) + len(diff_distribution[diff_distribution >=
+                                              (-1)*value])/len(
+            diff_distribution)
+    else:
+        p_value = len(diff_distribution[diff_distribution >=
+                                              value])/len(
+            diff_distribution) + len(diff_distribution[diff_distribution <=
+                                              (-1)*value])/len(
+            diff_distribution)
+    return p_value
+
+
+def bootstrap_distribution(sample, statistic=np.var, nboot=5000):
+    sample_stat = []
+    for bssample in range(nboot):
+        sample_vals = np.random.choice(sample, size=len(sample),
+                                       replace=True)
+        sample_stat.append(statistic(sample_vals))
+    return sample_stat
+
+
+def binning_focus(binning, patient):
+
+    binning_focus = 'X'
+    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 focus == 'NA':
+        print('Error: patient {} not in foci list.'.format(patient))
+        return 1
+
+    if binning.startswith('left') and focus == 'L' or \
+            binning.startswith('right') and \
+            focus == 'R':
+        binning_focus = 'focal'
+    elif binning.startswith('left') and focus == 'R' or \
+            binning.startswith('right') and \
+            focus == 'L':
+        binning_focus = 'contra-lateral'
+    elif focus == 'B' and binning.startswith('left'):
+        binning_focus = 'B_L'
+    elif focus == 'B' and binning.startswith('right'):
+        binning_focus = 'B_R'
+    elif '?' in focus and binning.startswith('right'):
+        binning_focus = 'U_R'
+    elif '?' in focus and binning.startswith('left'):
+        binning_focus = 'U_L'
+    elif (binning.startswith('L') and focus == 'L') or (
+            binning.startswith('R') and focus == 'R'):
+        binning_focus = 'focal'
+    elif (binning.startswith('L') and focus == 'R') or (
+            binning.startswith('R') and focus == 'L'):
+        binning_focus = 'contra-lateral'
+
+    return binning_focus

python_code/analysis_PlosCB/readme_code.md

@@ -0,0 +1,22 @@
+# Summary of python scripts
+
+
+## Running the Multistep Regression Analysis
+
+* run_analysis_fullrecs.py (Analysis of full recordings from all brain areas)
+* run_analysis_timeresolved.py (Timeresolved analysis, including the 2-part analysis and the sliding window analysis. Note that the sliding window analysis of all recordings is computationally expensive.)
+
+
+## Creating the figures
+
+* create_figure_1.py
+* create_figure_2.py
+* create_supp_figures.py
+
+
+## Other functions / utilities called by the above mentioned scipts
+
+* mr_epilepsy_fullrecs.py
+* mr_epilepsy_timeresolved.py
+* mr_utilities.py
+* hierarchical_model.py

python_code/analysis_PlosCB/run_analysis_fullrecs.py

@@ -0,0 +1,69 @@
+import mr_epilepsy_fullrec
+import os
+
+
+def run_mr_regions():
+
+    for rectype in ['preictal', 'interictal']:
+
+        data_path = '../Data/{}/'.format(rectype)
+        result_path = '../Results/{}/singleunits_regions/'.format(rectype)
+
+        # -----------------------------------------------------------------------#
+        # create result directories
+        # -----------------------------------------------------------------------#
+        recordings = [dI for dI in os.listdir(data_path) if os.path.isdir(
+            os.path.join(data_path, dI))]
+
+        if not os.path.exists(result_path):
+            os.makedirs(result_path)
+        for rec in recordings:
+            if not os.path.exists(result_path + rec):
+                os.makedirs(result_path + rec)
+
+        # -----------------------------------------------------------------------#
+        # run mr analysis
+        # -----------------------------------------------------------------------#
+        outputfilename = 'activity_brainregions_SU_4ms.pkl'
+
+        mr_epilepsy_fullrec.mr_analysis_confint(data_path, result_path,
+                    outputfilename=outputfilename, recordings='all',
+                    binnings=[], ksteps=(1, 400), fit_method='offset')
+
+
+def run_mr_hemispheres():
+
+    for rectype in ['preictal', 'interictal']:
+        data_path = '../Data/{}/'.format(rectype)
+        result_path = '../Results/{}/singleunits/'.format(rectype)
+
+        # -----------------------------------------------------------------------#
+        # create result directories
+        # -----------------------------------------------------------------------#
+        recordings = [dI for dI in os.listdir(data_path) if os.path.isdir(
+            os.path.join(data_path, dI))]
+
+        if not os.path.exists(result_path):
+            os.makedirs(result_path)
+        for rec in recordings:
+            if not os.path.exists(result_path + rec):
+                os.makedirs(result_path + rec)
+
+        # -----------------------------------------------------------------------#
+        # run mr analysis
+        # -----------------------------------------------------------------------#
+        outputfilename = 'activity_SU_4ms.pkl'
+
+        mr_epilepsy_fullrec.mr_analysis_confint(data_path, result_path,
+                            outputfilename=outputfilename, recordings='all',
+                            binnings=['left', 'right'], ksteps=(1, 400),
+                            fit_method='offset')
+
+
+def main():
+    run_mr_hemispheres()
+    #run_mr_regions()
+
+
+if __name__ == '__main__':
+    main()

python_code/analysis_PlosCB/run_analysis_timeresolved.py

@@ -0,0 +1,112 @@
+import os
+import mr_epilepsy_timeresolved
+
+
+def run_2parts_regions():
+    # -----------------------------------------------------------------------#
+    # preictal data 4ms
+    # -----------------------------------------------------------------------#
+    data_path = '../Data/preictal/'
+    result_path = '../Results/preictal/singleunits_regions/'
+
+    # -----------------------------------------------------------------------#
+    # create result directories
+    # -----------------------------------------------------------------------#
+    recordings = [dI for dI in os.listdir(data_path) if os.path.isdir(
+        os.path.join(data_path, dI))]
+
+    if not os.path.exists(result_path):
+        os.makedirs(result_path)
+    for rec in recordings:
+        if not os.path.exists(result_path + rec):
+            os.makedirs(result_path + rec)
+
+    # -----------------------------------------------------------------------#
+    # run mr analysis
+    # -----------------------------------------------------------------------#
+    outputfilename = 'activity_brainregions_SU_4ms.pkl'
+    x = 2
+    mr_epilepsy_timeresolved.mt_x_parts(x, outputfilename, data_path,
+       result_path, recordings=recordings, binnings='all', ksteps=(1, 400),
+       fit_method="offset")
+
+
+def run_2parts_hemispheres():
+    # -----------------------------------------------------------------------#
+    # preictal data 4ms
+    # -----------------------------------------------------------------------#
+    data_path = '../Data/preictal/'
+    result_path = '../Results/preictal/singleunits/'
+
+    # -----------------------------------------------------------------------#
+    # create result directories
+    # -----------------------------------------------------------------------#
+    recordings = [dI for dI in os.listdir(data_path) if os.path.isdir(
+        os.path.join(data_path, dI))]
+
+    if not os.path.exists(result_path):
+        os.makedirs(result_path)
+    for rec in recordings:
+        if not os.path.exists(result_path + rec):
+            os.makedirs(result_path + rec)
+
+    # -----------------------------------------------------------------------#
+    # run timeresolved mr analysis
+    # -----------------------------------------------------------------------#
+    x = 2
+    outputfilename = "activity_SU_4ms.pkl"
+    mr_epilepsy_timeresolved.mt_x_parts(x, outputfilename, data_path,
+       result_path, recordings=recordings, binnings=['left', 'right'],
+       ksteps=(1, 400), fit_method="offset")
+
+
+def run_timeresolved_hemispheres():
+    # -----------------------------------------------------------------------#
+    # preictal data 4ms
+    # -----------------------------------------------------------------------#
+    data_path = '../Data/preictal/'
+    result_path = '../Results/preictal/singleunits/'
+
+    # -----------------------------------------------------------------------#
+    # create result directories
+    # -----------------------------------------------------------------------#
+    recordings = [dI for dI in os.listdir(data_path) if os.path.isdir(
+        os.path.join(data_path, dI))]
+
+    if not os.path.exists(result_path):
+        os.makedirs(result_path)
+    for rec in recordings:
+        if not os.path.exists(result_path + rec):
+            os.makedirs(result_path + rec)
+
+    # -----------------------------------------------------------------------#
+    # run timeresolved mr analysis
+    # -----------------------------------------------------------------------#
+    ksteps = (1, 400)
+    for binning in ['left', 'right']:
+
+        for rec in recordings[60:]:
+            windowsize = 20000
+            windowstep = 500
+            binning = binning
+            fit_method = 'offset'
+            dt = 4
+            mt_file = '{}{}/mt_results_{}_{}_kmin{}_kmax{}_winsize{}_winstep{}_dt{}.pkl'.format(
+                result_path, rec, binning, fit_method, ksteps[0], ksteps[1],
+                windowsize, windowstep, dt)
+
+            if not os.path.isfile(mt_file):
+                mr_epilepsy_timeresolved.analyse_timeresolved_mr(data_path,
+                    result_path, recordings=[rec], binnings=[binning],
+                    ksteps=ksteps, windowsize=20000, windowstep=500,
+                    fit_method='offset')
+
+
+def main():
+    run_timeresolved_hemispheres()
+    #run_2parts_hemispheres()
+    #run_2parts_regions()
+
+
+if __name__ == '__main__':
+    main()