<reponame>robbisg/mvpa_itab_wu
from scipy.io import loadmat
import h5py
import pandas as pd
import seaborn as sns
from pyitab.results.dataframe import apply_function


diag_kws = {'bins': 50, 'color': 'teal', 'alpha': 0.4, 'edgecolor':None}
plot_kws = {'color': 'teal', 'edgecolor': None, 'alpha': 0.1}


path = "/media/robbis/DATA/meg/reftep/derivatives/phastimate/"

columns = ['phases32', 'hjort', 'predictedyule32', 
           'amplitudes200', 'amplitudes32', 'mep1']


def read_channels(mat, key='chlabels'):
    test = mat['chlabels']
    channels = list()
    for st in test[0]:
        obj = mat[st]
        str1 = ''.join(chr(i) for i in obj[:])
        channels.append(str1)
    return np.array(channels)


def scaler(x, b):
    factor = b / (2 * np.pi)
    return factor * (x + np.pi ) - .5


cmap = pl.cm.magma

columns = ['prediction', 'phases']

for i in range(9):
    sub = "sub-%03d" % (i+1)
    filename = os.path.join(path, sub, sub+"_space-sensor_window-500_atlas-subject_band-mu_prediction.mat")
    mat = h5py.File(filename, 'r')
    vector = []

    channels = read_channels(mat)
    idx = np.nonzero(channels == 'C3')[0]
    idx = int(idx)

    mep = np.log10(mat['AmpsMclean'][0])
    median = np.median(mep)

    vector = []
    fig, axes = pl.subplots(2, 2)

    ax = axes[0, :]
    for j, k in enumerate(columns):

        prediction = mat[k]
        data = prediction[()][:, idx, :]

        high_mep = data[:, mep > median]
        low_mep = data[:, mep <= median]

        ds = pd.DataFrame(data.T)
        mep_type = np.ones_like(mep) * -1
        mep_type[mep > median] = 1

        ds['type'] = mep_type
        dsm = ds.melt(value_vars=range(ds.shape[1] - 1), 
                      value_name=k, 
                      var_name='time', 
                      id_vars='type')

        sns.lineplot(x='time', 
                     y=k, 
                     hue='type', 
                     data=dsm, 
                     ax=ax[j], 
                     ci='sd',
                     palette=['salmon', 'teal'])

        mm = np.max([np.max(np.abs(high_mep.mean(1))), np.max(np.abs(low_mep.mean(1)))]) + 1
        ax[j].vlines(251-64, mm, -mm, linestyle='dashed', color='gray')
        ax[j].vlines(251-32, mm, -mm, linestyle='dashed', color='red')

        stim_value = data[-32, :]
        vector.append(stim_value)

    # fig, ax = pl.subplots(1, 1)
    ax3 = axes[1, 0]
    vector = np.vstack(vector)
    df = pd.DataFrame(vector.T, columns=columns)

    df['mep'] = mep
    
    b = 8
    bins = pd.cut(df['phases'], bins=b)
    df['bin'] = [v.mid for v in bins]
    sns.boxenplot(data=df, x='bin', y='mep', color='g', showfliers=False, ax=ax3)
    df['phasez'] = scaler(df['phases'], b)
    sns.scatterplot(data=df, x='phasez', y='mep', edgecolor=None, ax=ax3, color='teal', alpha=0.5)

    ax4 = axes[1, 1]
    for b in range(4, 11):
        bins = pd.cut(df['phases'], bins=b)
        df['bin'] = [v.mid for v in bins]
    
        dfm = apply_function(df, keys=['bin'], attr='mep', fx= lambda x: np.median(x))
        c = (b - 4) / 11
        ax4.plot(dfm['bin'], dfm['mep'], '-o', color=cmap(c))
        if b == 8:
            ax4.plot(dfm['bin'], dfm['mep'], '-o', lw=5, color='red')
            ax3.plot(scaler(dfm['bin'], b), dfm['mep'], '-o', lw=5, color='red')
        
        if b == 10:
            ax4.plot(dfm['bin'], dfm['mep'], '-o', lw=5, color='black')
        
    ax4.set_ylabel('mep')
    ax4.set_xlabel('phase')

    mat.close()
from scipy.io import loadmat
import h5py
import pandas as pd
import seaborn as sns
from pyitab.results.dataframe import apply_function


diag_kws = {'bins': 50, 'color': 'teal', 'alpha': 0.4, 'edgecolor':None}
plot_kws = {'color': 'teal', 'edgecolor': None, 'alpha': 0.1}


path = "/media/robbis/DATA/meg/reftep/derivatives/phastimate/"

columns = ['phases32', 'hjort', 'predictedyule32', 
           'amplitudes200', 'amplitudes32', 'mep1']


def read_channels(mat, key='chlabels'):
    test = mat['chlabels']
    channels = list()
    for st in test[0]:
        obj = mat[st]
        str1 = ''.join(chr(i) for i in obj[:])
        channels.append(str1)
    return np.array(channels)


def scaler(x, b):
    factor = b / (2 * np.pi)
    return factor * (x + np.pi ) - .5


cmap = pl.cm.magma

columns = ['prediction', 'phases']

for i in range(9):
    sub = "sub-%03d" % (i+1)
    filename = os.path.join(path, sub, sub+"_space-sensor_window-500_atlas-subject_band-mu_prediction.mat")
    mat = h5py.File(filename, 'r')
    vector = []

    channels = read_channels(mat)
    idx = np.nonzero(channels == 'C3')[0]
    idx = int(idx)

    mep = np.log10(mat['AmpsMclean'][0])
    median = np.median(mep)

    vector = []
    fig, axes = pl.subplots(2, 2)

    ax = axes[0, :]
    for j, k in enumerate(columns):

        prediction = mat[k]
        data = prediction[()][:, idx, :]

        high_mep = data[:, mep > median]
        low_mep = data[:, mep <= median]

        ds = pd.DataFrame(data.T)
        mep_type = np.ones_like(mep) * -1
        mep_type[mep > median] = 1

        ds['type'] = mep_type
        dsm = ds.melt(value_vars=range(ds.shape[1] - 1), 
                      value_name=k, 
                      var_name='time', 
                      id_vars='type')

        sns.lineplot(x='time', 
                     y=k, 
                     hue='type', 
                     data=dsm, 
                     ax=ax[j], 
                     ci='sd',
                     palette=['salmon', 'teal'])

        mm = np.max([np.max(np.abs(high_mep.mean(1))), np.max(np.abs(low_mep.mean(1)))]) + 1
        ax[j].vlines(251-64, mm, -mm, linestyle='dashed', color='gray')
        ax[j].vlines(251-32, mm, -mm, linestyle='dashed', color='red')

        stim_value = data[-32, :]
        vector.append(stim_value)

    # fig, ax = pl.subplots(1, 1)
    ax3 = axes[1, 0]
    vector = np.vstack(vector)
    df = pd.DataFrame(vector.T, columns=columns)

    df['mep'] = mep
    
    b = 8
    bins = pd.cut(df['phases'], bins=b)
    df['bin'] = [v.mid for v in bins]
    sns.boxenplot(data=df, x='bin', y='mep', color='g', showfliers=False, ax=ax3)
    df['phasez'] = scaler(df['phases'], b)
    sns.scatterplot(data=df, x='phasez', y='mep', edgecolor=None, ax=ax3, color='teal', alpha=0.5)

    ax4 = axes[1, 1]
    for b in range(4, 11):
        bins = pd.cut(df['phases'], bins=b)
        df['bin'] = [v.mid for v in bins]
    
        dfm = apply_function(df, keys=['bin'], attr='mep', fx= lambda x: np.median(x))
        c = (b - 4) / 11
        ax4.plot(dfm['bin'], dfm['mep'], '-o', color=cmap(c))
        if b == 8:
            ax4.plot(dfm['bin'], dfm['mep'], '-o', lw=5, color='red')
            ax3.plot(scaler(dfm['bin'], b), dfm['mep'], '-o', lw=5, color='red')
        
        if b == 10:
            ax4.plot(dfm['bin'], dfm['mep'], '-o', lw=5, color='black')
        
    ax4.set_ylabel('mep')
    ax4.set_xlabel('phase')

    mat.close()
# fig, ax = pl.subplots(1, 1)