5-fixed_delay_hyperbrain_networks.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## License\n",
    "This file is part of the project hyperMusic. All of hyperMusic code is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. hyperMusic is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with hyperMusic. If not, see <https://www.gnu.org/licenses/>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import glob\n",
    "import sys \n",
    "\n",
    "import bct as bct\n",
    "import matplotlib.pyplot as plt\n",
    "from mne.stats import fdr_correction\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy as sp\n",
    "import seaborn as sns\n",
    "sns.set_style(\"white\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First, we get some helper and anaysis functions going on..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "This function is used to create a distribution based on permutated paired t-tests\n",
    "condition = 'Insula Left,Parietal Right,leader,follower,beta,beta'\n",
    "for condition in conditions:\n",
    "    perm_cond = grand_grand_df[grand_grand_df['conditions'] == condition]\n",
    "    perm_ste = perm_cond[['ste', 'baseline_ste']].values\n",
    "    experimental_value, _ = sp.stats.ttest_rel(perm_ste[:, 0], perm_ste[:, 1]) # experimental value of paired t-test\n",
    "    perm_distribution = np.array([monte_carlo_dist(perm_ste) for i in range(1000)])\n",
    "    p_value = (experimental_value < perm_distribution).sum()/64.\n",
    "    if p_value == 0.0:\n",
    "        out_df.loc[out_df[out_df['conditions'] == condition].index[0], 'p_value'] = 1/64.\n",
    "    else:\n",
    "        out_df.loc[out_df[out_df['conditions'] == condition].index[0], 'p_value'] = p_value\n",
    "'''\n",
    "\n",
    "def info_source(row):\n",
    "    if row['sub_source'] == row['sub_leader']: \n",
    "        return 'leader'\n",
    "    else:\n",
    "        return 'follower'\n",
    "\n",
    "def info_target(row):\n",
    "    if row['sub_target'] == row['sub_leader']:\n",
    "        return 'leader'\n",
    "    else:\n",
    "        return 'follower'\n",
    "    \n",
    "def monte_carlo_dist_unpaired(matrix_inputs):\n",
    "    # concatenate both columns into one long vector\n",
    "    n_vector = matrix_inputs.shape[0]\n",
    "    long_vector = matrix_inputs.reshape((n_vector*2))\n",
    "\n",
    "    # permute values\n",
    "    perm_vector = np.random.permutation(long_vector)\n",
    "\n",
    "    # unpack arrays\n",
    "    a_out = perm_vector[n_vector:]\n",
    "    b_out = perm_vector[:n_vector]\n",
    "\n",
    "    # get ind t-test\n",
    "    t_out, _ = sp.stats.ttest_ind(a_out, b_out)\n",
    "    return t_out\n",
    "\n",
    "# this bit of code runs paired t-tests instead of independent tests for the MonteCarlo simulation\n",
    "def monte_carlo_dist(matrix_inputs):\n",
    "    \"\"\"Randomly swap entries in two arrays.\"\"\"\n",
    "    # indices to swap\n",
    "    swap_inds = np.random.random(size = matrix_inputs.shape[0]) < 0.5\n",
    "    \n",
    "    # unpack the matrix\n",
    "    a = matrix_inputs[:, 0]\n",
    "    b = matrix_inputs[:, 1]\n",
    "    \n",
    "    # make copies of arrays a and b for output\n",
    "    a_out = np.copy(a)\n",
    "    b_out = np.copy(b)\n",
    "    \n",
    "    # swap values\n",
    "    a_out[swap_inds] = b[swap_inds]\n",
    "    b_out[swap_inds] = a[swap_inds]\n",
    "    \n",
    "    # calculate paired t-test\n",
    "    t_out, _ = sp.stats.ttest_rel(a_out, b_out)\n",
    "\n",
    "    return t_out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Table 1: Participants' information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>pair</th>\n",
       "      <th>age_1</th>\n",
       "      <th>age_2</th>\n",
       "      <th>formal_1</th>\n",
       "      <th>formal_2</th>\n",
       "      <th>ensemble_1</th>\n",
       "      <th>ensemble_2</th>\n",
       "      <th>conservatory_1</th>\n",
       "      <th>conservatory_2</th>\n",
       "      <th>profesional_1</th>\n",
       "      <th>profesional_2</th>\n",
       "      <th>current_ensemble_1</th>\n",
       "      <th>current_ensemble_2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>P03</td>\n",
       "      <td>23</td>\n",
       "      <td>20</td>\n",
       "      <td>15</td>\n",
       "      <td>12</td>\n",
       "      <td>10</td>\n",
       "      <td>8</td>\n",
       "      <td>Y</td>\n",
       "      <td>Y</td>\n",
       "      <td>N</td>\n",
       "      <td>Y</td>\n",
       "      <td>Y</td>\n",
       "      <td>N</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>P04</td>\n",
       "      <td>29</td>\n",
       "      <td>32</td>\n",
       "      <td>10</td>\n",
       "      <td>10</td>\n",
       "      <td>7</td>\n",
       "      <td>5</td>\n",
       "      <td>Y</td>\n",
       "      <td>Y</td>\n",
       "      <td>N</td>\n",
       "      <td>N</td>\n",
       "      <td>N</td>\n",
       "      <td>N</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>P05</td>\n",
       "      <td>20</td>\n",
       "      <td>18</td>\n",
       "      <td>10</td>\n",
       "      <td>15</td>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>Y</td>\n",
       "      <td>Y</td>\n",
       "      <td>N</td>\n",
       "      <td>N</td>\n",
       "      <td>Y</td>\n",
       "      <td>Y</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>P08</td>\n",
       "      <td>21</td>\n",
       "      <td>24</td>\n",
       "      <td>15</td>\n",
       "      <td>15</td>\n",
       "      <td>4</td>\n",
       "      <td>8</td>\n",
       "      <td>Y</td>\n",
       "      <td>Y</td>\n",
       "      <td>Y</td>\n",
       "      <td>N</td>\n",
       "      <td>N</td>\n",
       "      <td>N</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>P11</td>\n",
       "      <td>38</td>\n",
       "      <td>28</td>\n",
       "      <td>15</td>\n",
       "      <td>15</td>\n",
       "      <td>20</td>\n",
       "      <td>13</td>\n",
       "      <td>Y</td>\n",
       "      <td>Y</td>\n",
       "      <td>N</td>\n",
       "      <td>Y</td>\n",
       "      <td>N</td>\n",
       "      <td>N</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>P09</td>\n",
       "      <td>21</td>\n",
       "      <td>23</td>\n",
       "      <td>14</td>\n",
       "      <td>18</td>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>Y</td>\n",
       "      <td>Y</td>\n",
       "      <td>Y</td>\n",
       "      <td>Y</td>\n",
       "      <td>Y</td>\n",
       "      <td>N</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  pair  age_1  age_2  formal_1  formal_2  ensemble_1  ensemble_2  \\\n",
       "0  P03     23     20        15        12          10           8   \n",
       "1  P04     29     32        10        10           7           5   \n",
       "2  P05     20     18        10        15           6           4   \n",
       "3  P08     21     24        15        15           4           8   \n",
       "4  P11     38     28        15        15          20          13   \n",
       "5  P09     21     23        14        18           4           2   \n",
       "\n",
       "  conservatory_1 conservatory_2 profesional_1 profesional_2  \\\n",
       "0              Y              Y             N             Y   \n",
       "1              Y              Y             N             N   \n",
       "2              Y              Y             N             N   \n",
       "3              Y              Y             Y             N   \n",
       "4              Y              Y             N             Y   \n",
       "5              Y              Y             Y             Y   \n",
       "\n",
       "  current_ensemble_1 current_ensemble_2  \n",
       "0                  Y                  N  \n",
       "1                  N                  N  \n",
       "2                  Y                  Y  \n",
       "3                  N                  N  \n",
       "4                  N                  N  \n",
       "5                  Y                  N  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# get participants data\n",
    "participants = pd.read_csv('data/participant_music.csv')\n",
    "participants = participants[participants['Pair'] != 'P01']\n",
    "participants = participants[participants['Pair'] != 'P02']\n",
    "participants = participants[participants['Pair'] != 'P06']\n",
    "participants = participants.reset_index()\n",
    "del participants['index']\n",
    "\n",
    "pairs = participants['Pair'].tolist()\n",
    "n = len(pairs)/2\n",
    "sub_df = pd.DataFrame(data=np.zeros((n, 13)), \n",
    "                      columns=['pair', 'age_1', 'age_2', 'formal_1', \n",
    "                               'formal_2', 'ensemble_1', 'ensemble_2', \n",
    "                               'conservatory_1', 'conservatory_2', \n",
    "                               'profesional_1', 'profesional_2',\n",
    "                               'current_ensemble_1', 'current_ensemble_2'])              \n",
    "sub_df['pair'] = pairs[::2]\n",
    "sub_df['age_1'] = participants['age'][::2].tolist()\n",
    "sub_df['age_2'] = participants['age'][1::2].tolist()\n",
    "sub_df['formal_1'] = participants['formal_training'][::2].tolist()\n",
    "sub_df['formal_2'] = participants['formal_training'][1::2].tolist()\n",
    "sub_df['ensemble_1'] = participants['ensemble_training'][::2].tolist()\n",
    "sub_df['ensemble_2'] = participants['ensemble_training'][1::2].tolist()\n",
    "sub_df['conservatory_1'] = participants['conservatory'][::2].tolist()\n",
    "sub_df['conservatory_2'] = participants['conservatory'][1::2].tolist()\n",
    "sub_df['current_ensemble_1'] = participants['ensemble'][::2].tolist()\n",
    "sub_df['current_ensemble_2'] = participants['ensemble'][1::2].tolist()\n",
    "sub_df['profesional_1'] = participants['profesional'][::2].tolist()\n",
    "sub_df['profesional_2'] = participants['profesional'][1::2].tolist()\n",
    "sub_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis 1: Participants descriptive statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "We tested 6 pianist dyads. \n",
      "Mean age: 24.750000 ± 5.643950\n",
      "Absolute mean age difference: 3.833333 ± 2.793842\n",
      "\n",
      "Mean formal music training: 13.666667 ± 2.460804\n",
      "Absolute mean formal music training difference: 2.000000 ± 2.081666\n",
      "\n",
      "Mean ensemble music training: 7.583333 ± 4.733891\n",
      "Absolute mean ensemble music training difference: 3.166667 ± 1.863390\n",
      "\n",
      "Conservatory education: 6/6 concordant\n",
      "Currently working as profesional musicians: 3/6 concordant\n",
      "Currently playing in an ensemble: 4/6 concordant\n"
     ]
    }
   ],
   "source": [
    "# pairs descriptive statistics: means and std's\n",
    "mean_age = sub_df[['age_1', 'age_2']].values.mean()\n",
    "std_age = sub_df[['age_1', 'age_2']].values.std()\n",
    "mean_play = sub_df[['formal_1', 'formal_2']].values.mean()\n",
    "std_play = sub_df[['formal_1', 'formal_2']].values.std()\n",
    "mean_ensemble = sub_df[['ensemble_1', 'ensemble_2']].values.mean()\n",
    "std_ensemble = sub_df[['ensemble_1', 'ensemble_2']].values.std()\n",
    "\n",
    "# pairs descriptive statistics: absolute mean differences and std's\n",
    "mean_diff = abs(sub_df[['age_2']].values - sub_df[['age_1']].values).mean()\n",
    "std_diff= abs(sub_df[['age_2']].values - sub_df[['age_1']].values).std()\n",
    "mean_diff_play = abs(sub_df[['formal_2']].values - sub_df[['formal_1']].values).mean()\n",
    "std_diff_play = abs(sub_df[['formal_2']].values - sub_df[['formal_1']].values).std()\n",
    "mean_diff_ens = abs(sub_df[['ensemble_2']].values - sub_df[['ensemble_1']].values).mean()\n",
    "std_diff_ens = abs(sub_df[['ensemble_2']].values - sub_df[['ensemble_1']].values).std()\n",
    "\n",
    "# pair descriptive statistics: characteristics\n",
    "conserv_conc = (sub_df['conservatory_1'] == sub_df['conservatory_2']).sum()\n",
    "prof_conc = (sub_df['profesional_1'] == sub_df['profesional_2']).sum()\n",
    "ensemble_conc = (sub_df['current_ensemble_1'] == sub_df['current_ensemble_2']).sum()\n",
    "\n",
    "descriptive_stats = (mean_age, std_age, mean_diff, std_diff,\n",
    "                     mean_play, std_play, mean_diff_play, std_diff_play,\n",
    "                     mean_ensemble, std_ensemble, mean_diff_ens, std_diff_ens, \n",
    "                     conserv_conc, prof_conc, ensemble_conc)\n",
    "\n",
    "print('''\n",
    "We tested 6 pianist dyads. \n",
    "Mean age: %f ± %f\n",
    "Absolute mean age difference: %f ± %f\n",
    "\n",
    "Mean formal music training: %f ± %f\n",
    "Absolute mean formal music training difference: %f ± %f\n",
    "\n",
    "Mean ensemble music training: %f ± %f\n",
    "Absolute mean ensemble music training difference: %f ± %f\n",
    "\n",
    "Conservatory education: %d/6 concordant\n",
    "Currently working as profesional musicians: %d/6 concordant\n",
    "Currently playing in an ensemble: %d/6 concordant'''%descriptive_stats)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis 2: Per pair connectivity hyperbrain networks (experimental data, baseline, and scrambled pairs)\n",
    "\n",
    "First, we take every person's csv's and dump it into grand csv files per delay"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1699200, 15)\n",
      "(86400, 9)\n",
      "(1368000, 17)\n",
      "(1699200, 15)\n",
      "(86400, 9)\n",
      "(1368000, 17)\n",
      "(1699200, 15)\n",
      "(86400, 9)\n",
      "(1368000, 17)\n"
     ]
    }
   ],
   "source": [
    "# We start by setting up a couple parameters\n",
    "delays = ['3', '30', '150'] # Only looking at 20.81ms, 200ms, 1000ms\n",
    "delays_seconds = ['20ms', '200ms', '1000ms']\n",
    "\n",
    "# loop through delays\n",
    "for i_delay in range(len(delays)):\n",
    "    \n",
    "    # Create grand matrix for the original data corresponding to that delay \n",
    "    file_list = glob.glob('/Users/hectorOrozco/Desktop/hM_analysis/6-STE_fixed_delayed/csv_delay_corrected/*'+delays[i_delay]+'.csv')\n",
    "    grand_df = pd.read_csv(file_list[0], index_col = 0)\n",
    "    for file in file_list[1:]:\n",
    "        temp = pd.read_csv(file, index_col = 0)\n",
    "        grand_df = grand_df.append(temp)\n",
    "\n",
    "    # Get info direction from the csv\n",
    "    grand_df['info_source'] = grand_df.apply(info_source, axis = 1)\n",
    "    grand_df['info_target'] = grand_df.apply(info_target, axis = 1)\n",
    "\n",
    "    # Output grand csv\n",
    "    grand_df.to_csv('results/5.fixed_delay_trial/grand_csv_fixed_delay'+delays_seconds[i_delay]+'.csv')\n",
    "    print(grand_df.shape)\n",
    "\n",
    "    # Create grand baseline matrix for corresponding delay \n",
    "    file_list = glob.glob('/Users/hectorOrozco/Desktop/hM_analysis/7-STE_fixed_delay_baseline/baseline_grand_matrices_fixed/*'+delays[i_delay]+'.csv')\n",
    "    baseline_df = pd.read_csv(file_list[0], index_col = 0)\n",
    "    for file in file_list[1:]:\n",
    "        temp = pd.read_csv(file, index_col = 0)\n",
    "        baseline_df = baseline_df.append(temp)\n",
    "\n",
    "    # Output grand csv\n",
    "    baseline_df.to_csv('results/5.fixed_delay_trial/baseline_csv_fixed_delay'+delays_seconds[i_delay]+'.csv')\n",
    "    print(baseline_df.shape)\n",
    "\n",
    "    # Create grand matrix for scrambled data corresponding to delay \n",
    "    file_list = glob.glob('/Users/hectorOrozco/Desktop/hM_analysis/8-STE_fixed_delayed_scrambled/fixed_csv_scrambled_matrices/*'+delays[i_delay]+'.csv')\n",
    "    scrambled_df = pd.read_csv(file_list[0], index_col = 0)\n",
    "\n",
    "    for file in file_list[1:]:\n",
    "        temp = pd.read_csv(file, index_col = 0)\n",
    "        scrambled_df = scrambled_df.append(temp)\n",
    "\n",
    "    # Get info direction from the csv\n",
    "    scrambled_df['info_source'] = scrambled_df.apply(info_source, axis = 1)\n",
    "    scrambled_df['info_target'] = scrambled_df.apply(info_target, axis = 1)\n",
    "\n",
    "    # Output grand csv\n",
    "    scrambled_df.to_csv('results/5.fixed_delay_trial/scrambled_csv_fixed_delay'+delays_seconds[i_delay]+'.csv')\n",
    "    print(scrambled_df.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We output at this point homophonic vs polyphonic pieces for baseline, scrambled pairs, and original data just to doublec check everything is Ok "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Experimental data\n",
    "delays = ['3', '30', '150'] # Only looking at 20.81ms, 200ms, 1000ms\n",
    "delays_seconds = ['20ms', '200ms', '1000ms']\n",
    "pairs = ['P03', 'P04', 'P05', 'P08', 'P09', 'P11']\n",
    "sources_aranged = ['Prefrontal Left', 'Insula Left', 'Motor Left', 'Temporal Left', 'Parietal Left', 'Occipital Left',\n",
    "                   'Prefrontal Right', 'Insula Right', 'Motor Right', 'Temporal Right', 'Parietal Right', 'Occipital Right']\n",
    "freq_aranged = ['delta', 'theta', 'alpha', 'beta', 'gamma']\n",
    "\n",
    "# Loop through all the delays and pairs \n",
    "for i_delay in range(len(delays)):\n",
    "    grand_df = pd.read_csv('results/5.fixed_delay_trial/grand_csv_fixed_delay'+delays_seconds[i_delay]+'.csv', index_col = 0)\n",
    "    for i_pair in range(len(pairs)):\n",
    "        pair_grand_df = grand_df[grand_df['pair'] == pairs[i_pair]].copy()\n",
    "\n",
    "        # Average homophonic and polyphonic pieces separately \n",
    "        pair_matrix = pair_grand_df.groupby(['neuro_source', 'neuro_target', 'info_source', 'info_target', 'freq_source', 'freq_target', \n",
    "                                             'duo_type'], as_index=False)['ste'].mean()\n",
    "\n",
    "        # Consider these variables as categorical so they are ordered correctly when creating the figure\n",
    "        pair_matrix.loc[:, 'neuro_source'] = pd.Categorical(pair_matrix['neuro_source'], \n",
    "                                                         categories=sources_aranged, ordered = True)\n",
    "        pair_matrix.loc[:, 'neuro_target'] = pd.Categorical(pair_matrix['neuro_target'], \n",
    "                                                         categories=sources_aranged, ordered = True)\n",
    "        pair_matrix.loc[:, 'freq_source'] = pd.Categorical(pair_matrix['freq_source'], \n",
    "                                                        categories=freq_aranged, ordered = True)\n",
    "        pair_matrix.loc[:, 'freq_target'] = pd.Categorical(pair_matrix['freq_target'], \n",
    "                                                        categories=freq_aranged, ordered = True)\n",
    "\n",
    "        homo_matrix = pair_matrix[pair_matrix['duo_type'] == 'h']\n",
    "        homo_matrix = homo_matrix.pivot_table(index = ['freq_source', 'info_source', 'neuro_source'], \n",
    "                                            columns = ['freq_target', 'info_target', 'neuro_target'], \n",
    "                                            values = 'ste')\n",
    "\n",
    "        # Create normal figure\n",
    "        fig, axs = plt.subplots(1,1, figsize=(30, 30))\n",
    "        sns.heatmap(homo_matrix, cmap='Blues', square=True, linewidths=.5, ax = axs, cbar_kws = {'shrink': 0.5})\n",
    "        plt.hlines(range(12, 121, 12), 120, 0, linewidth = 0.5, linestyle = '--')\n",
    "        plt.hlines(range(24, 121, 24), 120, 0, linewidth = 1)\n",
    "        plt.vlines(range(12, 121, 12), 120, 0, linewidth = 0.5, linestyle = '--')\n",
    "        plt.vlines(range(24, 121, 24), 120, 0, linewidth = 1)\n",
    "        axs.set_title('Homophonic, original data %s, %s delay' % ((pairs[i_pair]), delays_seconds[i_delay]), fontsize=15)\n",
    "        plt.tight_layout()\n",
    "        fig.savefig('results/5.fixed_delay_trial/'+pairs[i_pair]+'_homo_original_fixed_delay'+delays_seconds[i_delay]+'.png', \n",
    "                    bbox_inches = 'tight',pad_inches = 0.1)\n",
    "\n",
    "        # Set diagonal to zero\n",
    "        np.fill_diagonal(homo_matrix.values, 0) \n",
    "        fig, axs = plt.subplots(1,1, figsize=(30, 30))\n",
    "        sns.heatmap(homo_matrix, cmap='Blues', square=True, linewidths=.5, ax = axs, cbar_kws = {'shrink': 0.5})\n",
    "        plt.hlines(range(12, 121, 12), 120, 0, linewidth = 0.5, linestyle = '--')\n",
    "        plt.hlines(range(24, 121, 24), 120, 0, linewidth = 1)\n",
    "        plt.vlines(range(12, 121, 12), 120, 0, linewidth = 0.5, linestyle = '--')\n",
    "        plt.vlines(range(24, 121, 24), 120, 0, linewidth = 1)\n",
    "        axs.set_title('Homophonic, original data, diag = 0 %s, %s delay' % ((pairs[i_pair]), delays_seconds[i_delay]), fontsize=15)\n",
    "        fig.savefig('results/5.fixed_delay_trial/'+pairs[i_pair]+'_homo_diag0_fixed_delay'+delays_seconds[i_delay]+'.png', \n",
    "                    bbox_inches = 'tight',pad_inches = 0.1)\n",
    "\n",
    "        # Set within connections to zero\n",
    "        between_matrix = pair_matrix[pair_matrix['duo_type'] == 'h'].copy()\n",
    "        between_matrix = between_matrix[between_matrix['info_source'] != between_matrix['info_target']]\n",
    "        between_matrix = between_matrix.pivot_table(index = ['freq_source', 'info_source', 'neuro_source'], \n",
    "                                            columns = ['freq_target', 'info_target', 'neuro_target'], \n",
    "                                            values = 'ste')\n",
    "        fig, axs = plt.subplots(1,1, figsize=(30, 30))\n",
    "        sns.heatmap(between_matrix, cmap='Blues', square=True, linewidths=.5, ax = axs, cbar_kws = {'shrink': 0.5})\n",
    "        plt.hlines(range(12, 121, 12), 120, 0, linewidth = 0.5, linestyle = '--')\n",
    "        plt.hlines(range(24, 121, 24), 120, 0, linewidth = 1)\n",
    "        plt.vlines(range(12, 121, 12), 120, 0, linewidth = 0.5, linestyle = '--')\n",
    "        plt.vlines(range(24, 121, 24), 120, 0, linewidth = 1)\n",
    "        axs.set_title('Homophonic, original data, between connections, %s, %s delay' % ((pairs[i_pair]), delays_seconds[i_delay]), fontsize=15)\n",
    "        fig.savefig('results/5.fixed_delay_trial/'+pairs[i_pair]+'_homo_between_fixed_delay'+delays_seconds[i_delay]+'.png', \n",
    "                    bbox_inches = 'tight',pad_inches = 0.1)\n",
    "        plt.close('all')\n",
    "        \n",
    "        poly_matrix = pair_matrix[pair_matrix['duo_type'] == 'p']\n",
    "        poly_matrix = poly_matrix.pivot_table(index = ['freq_source', 'info_source', 'neuro_source'], \n",
    "                                            columns = ['freq_target', 'info_target', 'neuro_target'], \n",
    "                                            values = 'ste')\n",
    "\n",
    "        # Create normal figure\n",
    "        fig, axs = plt.subplots(1,1, figsize=(30, 30))\n",
    "        sns.heatmap(poly_matrix, cmap='Blues', square=True, linewidths=.5, ax = axs, cbar_kws = {'shrink': 0.5})\n",
    "        plt.hlines(range(12, 121, 12), 120, 0, linewidth = 0.5, linestyle = '--')\n",
    "        plt.hlines(range(24, 121, 24), 120, 0, linewidth = 1)\n",
    "        plt.vlines(range(12, 121, 12), 120, 0, linewidth = 0.5, linestyle = '--')\n",
    "        plt.vlines(range(24, 121, 24), 120, 0, linewidth = 1)\n",
    "        axs.set_title('Polyphonic, original data %s, %s delay' % ((pairs[i_pair]), delays_seconds[i_delay]), fontsize=15)\n",
    "        plt.tight_layout()\n",
    "        fig.savefig('results/5.fixed_delay_trial/'+pairs[i_pair]+'_poly_original_fixed_delay'+delays_seconds[i_delay]+'.png', \n",
    "                    bbox_inches = 'tight',pad_inches = 0.1)\n",
    "\n",
    "        # Set diagonal to zero\n",
    "        np.fill_diagonal(poly_matrix.values, 0) \n",
    "        fig, axs = plt.subplots(1,1, figsize=(30, 30))\n",
    "        sns.heatmap(poly_matrix, cmap='Blues', square=True, linewidths=.5, ax = axs, cbar_kws = {'shrink': 0.5})\n",
    "        plt.hlines(range(12, 121, 12), 120, 0, linewidth = 0.5, linestyle = '--')\n",
    "        plt.hlines(range(24, 121, 24), 120, 0, linewidth = 1)\n",
    "        plt.vlines(range(12, 121, 12), 120, 0, linewidth = 0.5, linestyle = '--')\n",
    "        plt.vlines(range(24, 121, 24), 120, 0, linewidth = 1)\n",
    "        axs.set_title('Polyphonic, original data, diag = 0 %s, %s delay' % ((pairs[i_pair]), delays_seconds[i_delay]), fontsize=15)\n",
    "        fig.savefig('results/5.fixed_delay_trial/'+pairs[i_pair]+'_poly_diag0_fixed_delay'+delays_seconds[i_delay]+'.png', \n",
    "                    bbox_inches = 'tight',pad_inches = 0.1)\n",
    "\n",
    "        # Set within connections to zero\n",
    "        poly_matrix = pair_matrix[pair_matrix['duo_type'] == 'h'].copy()\n",
    "        poly_matrix = poly_matrix[poly_matrix['info_source'] != poly_matrix['info_target']]\n",
    "        poly_matrix = poly_matrix.pivot_table(index = ['freq_source', 'info_source', 'neuro_source'], \n",
    "                                            columns = ['freq_target', 'info_target', 'neuro_target'], \n",
    "                                            values = 'ste')\n",
    "        fig, axs = plt.subplots(1,1, figsize=(30, 30))\n",
    "        sns.heatmap(poly_matrix, cmap='Blues', square=True, linewidths=.5, ax = axs, cbar_kws = {'shrink': 0.5})\n",
    "        plt.hlines(range(12, 121, 12), 120, 0, linewidth = 0.5, linestyle = '--')\n",
    "        plt.hlines(range(24, 121, 24), 120, 0, linewidth = 1)\n",
    "        plt.vlines(range(12, 121, 12), 120, 0, linewidth = 0.5, linestyle = '--')\n",
    "        plt.vlines(range(24, 121, 24), 120, 0, linewidth = 1)\n",
    "        axs.set_title('Polyphonic, original data, between connections, %s, %s delay' % ((pairs[i_pair]), delays_seconds[i_delay]), fontsize=15)\n",
    "        fig.savefig('results/5.fixed_delay_trial/'+pairs[i_pair]+'_poly_between_fixed_delay'+delays_seconds[i_delay]+'.png', \n",
    "                    bbox_inches = 'tight',pad_inches = 0.1)\n",
    "        plt.close('all')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Baseline\n",
    "delays = ['3', '30', '150'] # Only looking at 20.81ms, 200ms, 1000ms\n",
    "delays_seconds = ['20ms', '200ms', '1000ms']\n",
    "pairs = ['P03', 'P04', 'P05', 'P08', 'P09', 'P11']\n",
    "sources_aranged = ['Prefrontal Left', 'Insula Left', 'Motor Left', 'Temporal Left', 'Parietal Left', 'Occipital Left',\n",
    "                   'Prefrontal Right', 'Insula Right', 'Motor Right', 'Temporal Right', 'Parietal Right', 'Occipital Right']\n",
    "freq_aranged = ['delta', 'theta', 'alpha', 'beta', 'gamma']\n",
    "\n",
    "# Loop through all the delays and pairs \n",
    "for i_delay in range(len(delays)):\n",
    "    grand_df = pd.read_csv('results/5.fixed_delay_trial/baseline_csv_fixed_delay'+delays_seconds[i_delay]+'.csv', index_col = 0)\n",
    "    for i_pair in range(len(pairs)):\n",
    "        pair_grand_df = grand_df[grand_df['pair'] == pairs[i_pair]].copy()\n",
    "        # Consider these variables as categorical so they are ordered correctly when creating the figure\n",
    "        pair_grand_df.loc[:, 'neuro_source'] = pd.Categorical(pair_grand_df['neuro_source'], \n",
    "                                                         categories=sources_aranged, ordered = True)\n",
    "        pair_grand_df.loc[:, 'neuro_target'] = pd.Categorical(pair_grand_df['neuro_target'], \n",
    "                                                         categories=sources_aranged, ordered = True)\n",
    "        pair_grand_df.loc[:, 'freq_source'] = pd.Categorical(pair_grand_df['freq_source'], \n",
    "                                                        categories=freq_aranged, ordered = True)\n",
    "        pair_grand_df.loc[:, 'freq_target'] = pd.Categorical(pair_grand_df['freq_target'], \n",
    "                                                        categories=freq_aranged, ordered = True)\n",
    "\n",
    "        baseline_matrix = pair_grand_df.pivot_table(index = ['freq_source', 'sub_source', 'neuro_source'], \n",
    "                                            columns = ['freq_target', 'sub_target', 'neuro_target'], \n",
    "                                            values = 'ste')\n",
    "\n",
    "        # Create normal figure\n",
    "        fig, axs = plt.subplots(1,1, figsize=(30, 30))\n",
    "        sns.heatmap(baseline_matrix, cmap='Blues', square=True, linewidths=.5, ax = axs, cbar_kws = {'shrink': 0.5})\n",
    "        plt.hlines(range(12, 121, 12), 120, 0, linewidth = 0.5, linestyle = '--')\n",
    "        plt.hlines(range(24, 121, 24), 120, 0, linewidth = 1)\n",
    "        plt.vlines(range(12, 121, 12), 120, 0, linewidth = 0.5, linestyle = '--')\n",
    "        plt.vlines(range(24, 121, 24), 120, 0, linewidth = 1)\n",
    "        axs.set_title('Baseline, original data %s, %s delay' % ((pairs[i_pair]), delays_seconds[i_delay]), fontsize=15)\n",
    "        plt.tight_layout()\n",
    "        fig.savefig('results/5.fixed_delay_trial/'+pairs[i_pair]+'_baseline_original_fixed_delay'+delays_seconds[i_delay]+'.png', \n",
    "                    bbox_inches = 'tight',pad_inches = 0.1)\n",
    "\n",
    "        # Set diagonal to zero\n",
    "        np.fill_diagonal(baseline_matrix.values, 0) \n",
    "        fig, axs = plt.subplots(1,1, figsize=(30, 30))\n",
    "        sns.heatmap(baseline_matrix, cmap='Blues', square=True, linewidths=.5, ax = axs, cbar_kws = {'shrink': 0.5})\n",
    "        plt.hlines(range(12, 121, 12), 120, 0, linewidth = 0.5, linestyle = '--')\n",
    "        plt.hlines(range(24, 121, 24), 120, 0, linewidth = 1)\n",
    "        plt.vlines(range(12, 121, 12), 120, 0, linewidth = 0.5, linestyle = '--')\n",
    "        plt.vlines(range(24, 121, 24), 120, 0, linewidth = 1)\n",
    "        axs.set_title('Baseline, original data, diag = 0 %s, %s delay' % ((pairs[i_pair]), delays_seconds[i_delay]), fontsize=15)\n",
    "        fig.savefig('results/5.fixed_delay_trial/'+pairs[i_pair]+'_baseline_diag0_fixed_delay'+delays_seconds[i_delay]+'.png', \n",
    "                    bbox_inches = 'tight',pad_inches = 0.1)\n",
    "\n",
    "        # Set within connections to zero\n",
    "        between_matrix = pair_grand_df[pair_grand_df['sub_source'] != pair_grand_df['sub_target']]\n",
    "        between_matrix = between_matrix.pivot_table(index = ['freq_source', 'sub_source', 'neuro_source'], \n",
    "                                            columns = ['freq_target', 'sub_target', 'neuro_target'], \n",
    "                                            values = 'ste')\n",
    "        fig, axs = plt.subplots(1,1, figsize=(30, 30))\n",
    "        sns.heatmap(between_matrix, cmap='Blues', square=True, linewidths=.5, ax = axs, cbar_kws = {'shrink': 0.5})\n",
    "        plt.hlines(range(12, 121, 12), 120, 0, linewidth = 0.5, linestyle = '--')\n",
    "        plt.hlines(range(24, 121, 24), 120, 0, linewidth = 1)\n",
    "        plt.vlines(range(12, 121, 12), 120, 0, linewidth = 0.5, linestyle = '--')\n",
    "        plt.vlines(range(24, 121, 24), 120, 0, linewidth = 1)\n",
    "        axs.set_title('Baseline, original data, between connections, %s, %s delay' % ((pairs[i_pair]), delays_seconds[i_delay]), fontsize=15)\n",
    "        fig.savefig('results/5.fixed_delay_trial/'+pairs[i_pair]+'_baseline_between_fixed_delay'+delays_seconds[i_delay]+'.png', \n",
    "                    bbox_inches = 'tight',pad_inches = 0.1)\n",
    "        plt.close('all')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/hectorOrozco/anaconda/lib/python2.7/site-packages/matplotlib/pyplot.py:516: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
      "  max_open_warning, RuntimeWarning)\n"
     ]
    }
   ],
   "source": [
    "# Scrambled pairs\n",
    "delays = ['3', '30', '150'] # Only looking at 20.81ms, 200ms, 1000ms\n",
    "delays_seconds = ['20ms', '200ms', '1000ms']\n",
    "pairs = ['P03', 'P04', 'P05', 'P08', 'P09', 'P11']\n",
    "sources_aranged = ['Prefrontal Left', 'Insula Left', 'Motor Left', 'Temporal Left', 'Parietal Left', 'Occipital Left',\n",
    "                   'Prefrontal Right', 'Insula Right', 'Motor Right', 'Temporal Right', 'Parietal Right', 'Occipital Right']\n",
    "freq_aranged = ['delta', 'theta', 'alpha', 'beta', 'gamma']\n",
    "\n",
    "# Loop through all the delays and pairs \n",
    "for i_delay in range(len(delays)):\n",
    "    grand_df = pd.read_csv('results/5.fixed_delay_trial/scrambled_csv_fixed_delay'+delays_seconds[i_delay]+'.csv', index_col = 0)\n",
    "    for i_pair in range(len(pairs)):\n",
    "        pair_grand_df = grand_df[grand_df['pair_a'] == pairs[i_pair]].copy()\n",
    "\n",
    "        # Average all pieces together\n",
    "        scrambled_matrix = pair_grand_df.groupby(['neuro_source', 'neuro_target', 'info_source', 'info_target', 'freq_source', 'freq_target'], as_index=False)['ste'].mean()\n",
    "\n",
    "        # Consider these variables as categorical so they are ordered correctly when creating the figure\n",
    "        scrambled_matrix.loc[:, 'neuro_source'] = pd.Categorical(scrambled_matrix['neuro_source'], \n",
    "                                                         categories=sources_aranged, ordered = True)\n",
    "        scrambled_matrix.loc[:, 'neuro_target'] = pd.Categorical(scrambled_matrix['neuro_target'], \n",
    "                                                         categories=sources_aranged, ordered = True)\n",
    "        scrambled_matrix.loc[:, 'freq_source'] = pd.Categorical(scrambled_matrix['freq_source'], \n",
    "                                                        categories=freq_aranged, ordered = True)\n",
    "        scrambled_matrix.loc[:, 'freq_target'] = pd.Categorical(scrambled_matrix['freq_target'], \n",
    "                                                        categories=freq_aranged, ordered = True)\n",
    "\n",
    "        scrambled_matrix = scrambled_matrix.pivot_table(index = ['freq_source', 'info_source', 'neuro_source'], \n",
    "                                            columns = ['freq_target', 'info_target', 'neuro_target'], \n",
    "                                            values = 'ste')\n",
    "\n",
    "        # Create normal figure\n",
    "        fig, axs = plt.subplots(1,1, figsize=(30, 30))\n",
    "        sns.heatmap(scrambled_matrix, cmap='Blues', square=True, linewidths=.5, ax = axs, cbar_kws = {'shrink': 0.5})\n",
    "        plt.hlines(range(12, 121, 12), 120, 0, linewidth = 0.5, linestyle = '--')\n",
    "        plt.hlines(range(24, 121, 24), 120, 0, linewidth = 1)\n",
    "        plt.vlines(range(12, 121, 12), 120, 0, linewidth = 0.5, linestyle = '--')\n",
    "        plt.vlines(range(24, 121, 24), 120, 0, linewidth = 1)\n",
    "        axs.set_title('Averaged pieces, scrambled data %s, %s delay' % ((pairs[i_pair]), delays_seconds[i_delay]), fontsize=15)\n",
    "        plt.tight_layout()\n",
    "        fig.savefig('results/5.fixed_delay_trial/'+pairs[i_pair]+'_avg_scrambled_fixed_delay'+delays_seconds[i_delay]+'.png', \n",
    "                    bbox_inches = 'tight',pad_inches = 0.1)\n",
    "\n",
    "        # Set within connections to zero\n",
    "        scrambled_between_matrix = pair_grand_df[pair_grand_df['sub_source'] != pair_grand_df['sub_target']].copy()\n",
    "        scrambled_between_matrix = scrambled_between_matrix.groupby(['neuro_source', 'neuro_target', 'info_source', 'info_target', 'freq_source', 'freq_target'], as_index=False)['ste'].mean()\n",
    "        scrambled_between_matrix = scrambled_between_matrix.pivot_table(index = ['freq_source', 'info_source', 'neuro_source'], \n",
    "                                            columns = ['freq_target', 'info_target', 'neuro_target'], \n",
    "                                            values = 'ste')\n",
    "        fig, axs = plt.subplots(1,1, figsize=(30, 30))\n",
    "        sns.heatmap(scrambled_between_matrix, cmap='Blues', square=True, linewidths=.5, ax = axs, cbar_kws = {'shrink': 0.5})\n",
    "        plt.hlines(range(12, 121, 12), 120, 0, linewidth = 0.5, linestyle = '--')\n",
    "        plt.hlines(range(24, 121, 24), 120, 0, linewidth = 1)\n",
    "        plt.vlines(range(12, 121, 12), 120, 0, linewidth = 0.5, linestyle = '--')\n",
    "        plt.vlines(range(24, 121, 24), 120, 0, linewidth = 1)\n",
    "        axs.set_title('Averaged pieces, scrambled data, between connections, %s, %s delay' % ((pairs[i_pair]), delays_seconds[i_delay]), fontsize=15)\n",
    "        fig.savefig('results/5.fixed_delay_trial/'+pairs[i_pair]+'_avg_scrambled_between_fixed_delay'+delays_seconds[i_delay]+'.png', \n",
    "                    bbox_inches = 'tight',pad_inches = 0.1)\n",
    "        plt.close('all')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis 3: Are there differences between playing and baseline and actual pairs and shuffled pairs? (Average across all pieces and only look at between connections)\n",
    "\n",
    "We use permutations to evaluate the statistical significance of the difference between baseline and playing together, and playing together vs scrambled participants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 936,
   "metadata": {},
   "outputs": [],
   "source": [
    "delays = ['3', '30', '150'] # Only looking at 20.81ms, 200ms, 1000ms\n",
    "delays_seconds = ['20ms', '200ms', '1000ms']\n",
    "\n",
    "for i_delay in range(len(delays)):\n",
    "    # get grand dataframes \n",
    "    playing_df = pd.read_csv('results/5.fixed_delay_trial/grand_csv_fixed_delay'+delays_seconds[i_delay]+'.csv', index_col = 0)\n",
    "    baseline_df = pd.read_csv('results/5.fixed_delay_trial/baseline_csv_fixed_delay'+delays_seconds[i_delay]+'.csv', index_col = 0)\n",
    "    scrambled_df = pd.read_csv('results/5.fixed_delay_trial/scrambled_csv_fixed_delay'+delays_seconds[i_delay]+'.csv', index_col = 0)\n",
    "\n",
    "    # get top .33% of the between connections based on the grand average \n",
    "    playing_mean_df = playing_df.groupby(['neuro_source', 'neuro_target', 'info_source', 'info_target', 'freq_source', 'freq_target', 'pair'], \n",
    "                                        as_index=False)['ste'].mean()\n",
    "    playing_mean_df = playing_mean_df.groupby(['neuro_source', 'neuro_target', 'info_source', 'info_target', 'freq_source', 'freq_target'], \n",
    "                                              as_index=False)['ste'].mean()\n",
    "    playing_mean_df = playing_mean_df[playing_mean_df['info_source'] != playing_mean_df['info_target']] #only take between connections\n",
    "    playing_mean_df = playing_mean_df.sort_values(by=['ste'], ascending = False).reset_index()\n",
    "    playing_top_df = playing_mean_df.iloc[:48, :] # Only taking the .33% top between connections per delay\n",
    "    conditions = playing_top_df.iloc[:,1:7].apply(lambda x: \",\".join(x.astype(str)), axis=1).tolist() #'Motor Right,Motor Right,leader,follower,beta,beta'\n",
    "    playing_df = playing_df.groupby(['neuro_source', 'neuro_target', 'info_source', 'info_target', 'freq_source', 'freq_target', 'pair'], \n",
    "                                        as_index=False)['ste'].mean()\n",
    "    playing_df['temp'] = playing_df.iloc[:,:6].apply(lambda x: \",\".join(x.astype(str)), axis=1).tolist()\n",
    "    playing_top_grand_df = playing_df[playing_df.temp.isin(conditions)].reset_index(drop = True)\n",
    "\n",
    "    # get these connections from baseline dataframe tricky part: baseline doesn't have a clear direction, so we take the average between source-target and target-source\n",
    "    baseline_df['temp'] = baseline_df.iloc[:,[4, 5, 2, 3]].apply(lambda x: \",\".join(x.astype(str)), axis=1)\n",
    "    baseline_df = baseline_df[baseline_df['sub_source'] != baseline_df['sub_target']]\n",
    "    baseline_top_grand_df = baseline_df[baseline_df.temp.isin(playing_top_df.iloc[:,[1, 2, 5, 6]].apply(lambda x: \",\".join(x.astype(str)), axis=1).tolist())].reset_index(drop=True)\n",
    "\n",
    "    # We get the average a-b and b-a from the baseline values \n",
    "    baseline_top_grand_df = baseline_top_grand_df.groupby(['neuro_source', 'neuro_target', 'freq_source', 'freq_target',\n",
    "                                                           'pair'], as_index=False)['ste'].mean()\n",
    "    baseline_top_grand_df.loc[:, 'temp'] = baseline_top_grand_df.iloc[:,[0, 1, 2, 3, 4]].apply(lambda x: \",\".join(x.astype(str)), axis=1).tolist()\n",
    "\n",
    "    # We create the \"grand grand\" data frame with a column for the playing value, one for baseline,and one for scrambled\n",
    "    playing_top_grand_df.loc[:,'temp'] = playing_top_grand_df.iloc[:,[0, 1, 4, 5, 6]].apply(lambda x: \",\".join(x.astype(str)), axis=1).tolist()\n",
    "    grand_grand_df = playing_top_grand_df.copy()\n",
    "    for i_baseline in range(baseline_top_grand_df.shape[0]):\n",
    "        for i_playing in range(grand_grand_df.shape[0]):\n",
    "            if baseline_top_grand_df.loc[i_baseline, 'temp'] == grand_grand_df.loc[i_playing, 'temp']:\n",
    "                grand_grand_df.loc[i_playing, 'baseline_ste'] = baseline_top_grand_df.loc[i_baseline, 'ste']\n",
    "\n",
    "    # Now we get these connections from the scrambled file\n",
    "    conditions = grand_grand_df.iloc[:,[0, 1, 2, 3, 4, 5]].apply(lambda x: \",\".join(x.astype(str)), axis=1).tolist()\n",
    "    grand_grand_df['conditions'] = conditions\n",
    "    scrambled_df = scrambled_df.groupby(['neuro_source', 'neuro_target', 'info_source', 'info_target', 'freq_source', 'freq_target', 'pair_a'], \n",
    "                                    as_index=False)['ste'].mean()\n",
    "    scrambled_df['conditions'] = scrambled_df.iloc[:,0:6].apply(lambda x: \",\".join(x.astype(str)), axis=1).tolist() # \"Insula Left,Parietal Right,leader,follower,beta,beta\"\n",
    "    scrambled_top_grand_df = scrambled_df[scrambled_df.conditions.isin(conditions)].reset_index(drop = True)\n",
    "    scrambled_top_grand_df.loc[:,'temp'] = scrambled_top_grand_df.iloc[:,:7].apply(lambda x: \",\".join(x.astype(str)), axis=1).tolist()\n",
    "    grand_grand_df.loc[:,'temp'] = grand_grand_df.iloc[:,:7].apply(lambda x: \",\".join(x.astype(str)), axis=1).tolist()\n",
    "    grand_grand_df['scrambled_ste'] = 0.0\n",
    "    for i_scrambled in range(scrambled_top_grand_df.shape[0]):\n",
    "        for i_playing in range(grand_grand_df.shape[0]):\n",
    "            if scrambled_top_grand_df.loc[i_scrambled, 'temp'] == grand_grand_df.loc[i_playing, 'temp']:\n",
    "                grand_grand_df.loc[i_playing, 'scrambled_ste'] = scrambled_top_grand_df.loc[i_scrambled, 'ste']\n",
    "\n",
    "    # We loop across the top .33% edges; first get a list of these edges\n",
    "    conditions = np.unique(conditions)\n",
    "    grand_grand_df = grand_grand_df.drop(['temp'], axis = 1)\n",
    "    out_df = grand_grand_df.groupby(['neuro_source', 'neuro_target', 'freq_source', 'freq_target', 'info_source', 'info_target', 'conditions'], as_index=False)[('ste', 'baseline_ste', 'scrambled_ste')].mean()\n",
    "    out_df['p_value_baseline'] = np.zeros(len(out_df))\n",
    "    out_df['p_value_baseline_fdr'] = np.zeros(len(out_df))\n",
    "    out_df['t_test_baseline'] = np.zeros(len(out_df))\n",
    "    out_df['p_value_scrambled'] = np.zeros(len(out_df))\n",
    "    out_df['p_value_scrambled_fdr'] = np.zeros(len(out_df))\n",
    "    out_df['t_test_scrambled'] = np.zeros(len(out_df))\n",
    "\n",
    "    # we loop across the edgest to get what we want; first baseline\n",
    "    for condition in conditions:\n",
    "        perm_cond = grand_grand_df[grand_grand_df['conditions'] == condition]\n",
    "        perm_ste = perm_cond[['ste', 'baseline_ste']].values\n",
    "        experimental_value, _ = sp.stats.ttest_rel(perm_ste[:, 0], perm_ste[:, 1]) # experimental value of ind t-test\n",
    "        perm_distribution = np.array([monte_carlo_dist(perm_ste) for i in range(500)])\n",
    "        perm_distribution = np.unique(perm_distribution)\n",
    "        p_value = (experimental_value < perm_distribution).sum()/64.\n",
    "        if p_value == 0.0:\n",
    "            out_df.loc[out_df[out_df['conditions'] == condition].index[0], 'p_value_baseline'] = 1/5000.\n",
    "        else:\n",
    "            out_df.loc[out_df[out_df['conditions'] == condition].index[0], 'p_value_baseline'] = p_value\n",
    "        out_df.loc[out_df[out_df['conditions'] == condition].index[0], 't_test_baseline'] = experimental_value  \n",
    "\n",
    "        # scrambled\n",
    "        perm_ste = perm_cond[['ste', 'scrambled_ste']].values\n",
    "        experimental_value, _ = sp.stats.ttest_rel(perm_ste[:, 0], perm_ste[:, 1]) # experimental value of ind t-test\n",
    "        perm_distribution = np.array([monte_carlo_dist(perm_ste) for i in range(500)])\n",
    "        perm_distribution = np.unique(perm_distribution)\n",
    "        p_value = (experimental_value < perm_distribution).sum()/64.\n",
    "        out_df.loc[out_df[out_df['conditions'] == condition].index[0], 'p_value_scrambled'] = p_value\n",
    "        out_df.loc[out_df[out_df['conditions'] == condition].index[0], 't_test_scrambled'] = experimental_value  \n",
    "\n",
    "    # now we do a false discovery rate procedure to control for fwer\n",
    "    reject_fdr, pval_fdr = fdr_correction(out_df['p_value_baseline'].values, alpha = 0.05, method = 'indep')\n",
    "    out_df['p_value_baseline_fdr'] = pval_fdr  \n",
    "\n",
    "    reject_fdr, pval_fdr = fdr_correction(out_df['p_value_scrambled'].values, alpha = 0.05, method = 'indep')\n",
    "    out_df['p_value_scrambled_fdr'] = pval_fdr \n",
    "\n",
    "    out_df.to_csv('results/6.fixed_betweenVSbaseline_scrambled/playing_baseline_scrambled'+delays_seconds[i_delay]+'.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because none of the between connections are significant, we will not explore the correlations between them and the PMPQ. We will only plot the averaged matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 339,
   "metadata": {},
   "outputs": [],
   "source": [
    "delays = ['3', '30', '150'] # Only looking at 20.81ms, 200ms, 1000ms\n",
    "delays_seconds = ['20ms', '200ms', '1000ms']\n",
    "sources_aranged = ['Prefrontal Left', 'Insula Left', 'Motor Left', 'Temporal Left', 'Parietal Left', 'Occipital Left',\n",
    "                   'Prefrontal Right', 'Insula Right', 'Motor Right', 'Temporal Right', 'Parietal Right', 'Occipital Right']\n",
    "freq_aranged = ['delta', 'theta', 'alpha', 'beta', 'gamma']\n",
    "pairs = ['P03', 'P04', 'P05', 'P08', 'P09', 'P11']\n",
    "plt.close('all')\n",
    "\n",
    "# Loop through all the delays and pairs \n",
    "for i_delay in range(len(delays)):\n",
    "\n",
    "    # get grand dataframes \n",
    "    playing_df = pd.read_csv('results/5.fixed_delay_trial/grand_csv_fixed_delay'+delays_seconds[i_delay]+'.csv', index_col = 0)\n",
    "    baseline_df = pd.read_csv('results/5.fixed_delay_trial/baseline_csv_fixed_delay'+delays_seconds[i_delay]+'.csv', index_col = 0)\n",
    "    scrambled_df = pd.read_csv('results/5.fixed_delay_trial/scrambled_csv_fixed_delay'+delays_seconds[i_delay]+'.csv', index_col = 0)\n",
    "\n",
    "    # Playing\n",
    "    playing_matrix = playing_df.groupby(['neuro_source', 'neuro_target', 'info_source', 'info_target', 'freq_source', 'freq_target', 'pair'], \n",
    "                                as_index=False)['ste'].mean()\n",
    "    playing_matrix = playing_matrix.groupby(['neuro_source', 'neuro_target', 'info_source', 'info_target', 'freq_source', 'freq_target'], \n",
    "                                as_index=False)['ste'].mean()\n",
    "\n",
    "    # Consider these variables as categorical so they are ordered correctly when creating the figure\n",
    "    playing_matrix.loc[:, 'neuro_source'] = pd.Categorical(playing_matrix['neuro_source'], \n",
    "                                                     categories=sources_aranged, ordered = True)\n",
    "    playing_matrix.loc[:, 'neuro_target'] = pd.Categorical(playing_matrix['neuro_target'], \n",
    "                                                     categories=sources_aranged, ordered = True)\n",
    "    playing_matrix.loc[:, 'freq_source'] = pd.Categorical(playing_matrix['freq_source'], \n",
    "                                                    categories=freq_aranged, ordered = True)\n",
    "    playing_matrix.loc[:, 'freq_target'] = pd.Categorical(playing_matrix['freq_target'], \n",
    "                                                    categories=freq_aranged, ordered = True)\n",
    "\n",
    "    playing_matrix = playing_matrix.pivot_table(index = ['freq_source', 'info_source', 'neuro_source'], \n",
    "                                        columns = ['freq_target', 'info_target', 'neuro_target'], \n",
    "                                        values = 'ste')\n",
    "    vmax_full = playing_matrix.values.max()\n",
    "    if i_delay == 2:\n",
    "        vmax_full *= 1.3\n",
    "\n",
    "\n",
    "    # Full matrix\n",
    "    fig, axs = plt.subplots(1,1)\n",
    "    fig.set_size_inches(7.25, 7.25)\n",
    "    sns.set(font = 'sans-serif', font_scale=0.1) # you can create the labels here!\n",
    "    sns.heatmap(playing_matrix, cmap='Blues', square=True, linewidths=.1, ax = axs, cbar_kws = {'shrink': 0.5}, vmax = vmax_full)\n",
    "    plt.hlines(range(12, 120, 12), 120, 0, linewidth = 0.25, linestyle = '--')\n",
    "    plt.hlines(range(24, 120, 24), 120, 0, linewidth = 0.5)\n",
    "    plt.vlines(range(12, 120, 12), 120, 0, linewidth = 0.25, linestyle = '--')\n",
    "    plt.vlines(range(24, 120, 24), 120, 0, linewidth = 0.5)\n",
    "    sns.despine(fig = fig, ax = axs)\n",
    "    # use matplotlib.colorbar.Colorbar object\n",
    "    cbar = axs.collections[0].colorbar\n",
    "    # here set the labelsize by 20\n",
    "    cbar.ax.tick_params(labelsize=4, axis='both', which = 'both', direction = 'in', pad = 0.85)\n",
    "    axs.spines['top'].set_visible(False)\n",
    "    axs.spines['right'].set_visible(False)\n",
    "    axs.get_xaxis().tick_bottom()\n",
    "    axs.get_yaxis().tick_left()\n",
    "    axs.xaxis.set_ticks_position('none')\n",
    "    axs.yaxis.set_ticks_position('none')\n",
    "    axs.tick_params(axis='both', which = 'both', direction = 'in', pad = 0.5)\n",
    "    fig.savefig('results/6.fixed_betweenVSbaseline_scrambled/playing_avg_fixed_delay'+delays_seconds[i_delay]+'.eps', \n",
    "                        bbox_inches = 'tight',pad_inches = 0.1, dpi=300)\n",
    "\n",
    "    # Between\n",
    "    playing_between = playing_df.copy()\n",
    "    playing_between.loc[playing_between['info_source'] == playing_between['info_target'], 'ste'] = 0.0\n",
    "    playing_between = playing_between.groupby(['neuro_source', 'neuro_target', 'info_source', 'info_target', 'freq_source', 'freq_target', 'pair'], \n",
    "                                as_index=False)['ste'].mean()\n",
    "    playing_between = playing_between.groupby(['neuro_source', 'neuro_target', 'info_source', 'info_target', 'freq_source', 'freq_target'], \n",
    "                                as_index=False)['ste'].mean()\n",
    "\n",
    "    playing_between.loc[:, 'neuro_source'] = pd.Categorical(playing_between['neuro_source'], \n",
    "                                                     categories=sources_aranged, ordered = True)\n",
    "    playing_between.loc[:, 'neuro_target'] = pd.Categorical(playing_between['neuro_target'], \n",
    "                                                     categories=sources_aranged, ordered = True)\n",
    "    playing_between.loc[:, 'freq_source'] = pd.Categorical(playing_between['freq_source'], \n",
    "                                                    categories=freq_aranged, ordered = True)\n",
    "    playing_between.loc[:, 'freq_target'] = pd.Categorical(playing_between['freq_target'], \n",
    "                                                    categories=freq_aranged, ordered = True)\n",
    "\n",
    "    playing_between = playing_between.pivot_table(index = ['freq_source', 'info_source', 'neuro_source'], \n",
    "                                        columns = ['freq_target', 'info_target', 'neuro_target'], \n",
    "                                        values = 'ste')\n",
    "    vmax_between = playing_between.values.max() * 1.2\n",
    "\n",
    "    fig, axs = plt.subplots(1,1)\n",
    "    fig.set_size_inches(7.25, 7.25)\n",
    "    sns.set(font = 'sans-serif', font_scale=0.1) # you can create the labels here!\n",
    "    sns.heatmap(playing_between, cmap='Blues', square=True, linewidths=.1, ax = axs, cbar_kws = {'shrink': 0.5}, vmax = vmax_between)\n",
    "    plt.hlines(range(12, 120, 12), 120, 0, linewidth = 0.25, linestyle = '--')\n",
    "    plt.hlines(range(24, 120, 24), 120, 0, linewidth = 0.5)\n",
    "    plt.vlines(range(12, 120, 12), 120, 0, linewidth = 0.25, linestyle = '--')\n",
    "    plt.vlines(range(24, 120, 24), 120, 0, linewidth = 0.5)\n",
    "    sns.despine(fig = fig, ax = axs)\n",
    "    # use matplotlib.colorbar.Colorbar object\n",
    "    cbar = axs.collections[0].colorbar\n",
    "    # here set the labelsize by 20\n",
    "    cbar.ax.tick_params(labelsize=4, axis='both', which = 'both', direction = 'in', pad = 0.85)\n",
    "    axs.spines['top'].set_visible(False)\n",
    "    axs.spines['right'].set_visible(False)\n",
    "    axs.get_xaxis().tick_bottom()\n",
    "    axs.get_yaxis().tick_left()\n",
    "    axs.xaxis.set_ticks_position('none')\n",
    "    axs.yaxis.set_ticks_position('none')\n",
    "    axs.tick_params(axis='both', which = 'both', direction = 'in', pad = 0.5)\n",
    "    fig.savefig('results/6.fixed_betweenVSbaseline_scrambled/playing_avg_between_fixed_delay'+delays_seconds[i_delay]+'.eps', \n",
    "                        bbox_inches = 'tight',pad_inches = 0.1, dpi=300)\n",
    "    # Baseline\n",
    "    baseline_df['sub_source'] = baseline_df.iloc[:,[0, 7]].apply(lambda x: \",\".join(x.astype(str)), axis=1).tolist()    \n",
    "    baseline_df['sub_target'] = baseline_df.iloc[:,[1, 7]].apply(lambda x: \",\".join(x.astype(str)), axis=1).tolist()\n",
    "\n",
    "    sub_a_pairs = ['11A,P11', '3A,P03', '5A,P05', '5B,P04', '8A,P08', '9A,P09']\n",
    "    sub_b_pairs = ['11B,P11', '4B,P03', '5B,P05', '6B,P04', '8B,P08', '9B,P09']\n",
    "    baseline_df.loc[baseline_df['sub_source'].isin(sub_a_pairs), 'sub_source'] = 'a'\n",
    "    baseline_df.loc[baseline_df['sub_target'].isin(sub_a_pairs), 'sub_target'] = 'a'\n",
    "    baseline_df.loc[baseline_df['sub_source'].isin(sub_b_pairs), 'sub_source'] = 'b'\n",
    "    baseline_df.loc[baseline_df['sub_target'].isin(sub_b_pairs), 'sub_target'] = 'b'\n",
    "\n",
    "    # Average all edges\n",
    "    baseline_matrix = baseline_df.groupby(['neuro_source', 'neuro_target', 'sub_source', 'sub_target', 'freq_source', 'freq_target', 'pair'], \n",
    "                                as_index=False)['ste'].mean()\n",
    "    baseline_matrix = baseline_matrix.groupby(['neuro_source', 'neuro_target', 'sub_source', 'sub_target', 'freq_source', 'freq_target'], \n",
    "                                as_index=False)['ste'].mean()\n",
    "\n",
    "    # Consider these variables as categorical so they are ordered correctly when creating the figure\n",
    "    baseline_matrix.loc[:, 'neuro_source'] = pd.Categorical(baseline_matrix['neuro_source'], \n",
    "                                                     categories=sources_aranged, ordered = True)\n",
    "    baseline_matrix.loc[:, 'neuro_target'] = pd.Categorical(baseline_matrix['neuro_target'], \n",
    "                                                     categories=sources_aranged, ordered = True)\n",
    "    baseline_matrix.loc[:, 'freq_source'] = pd.Categorical(baseline_matrix['freq_source'], \n",
    "                                                    categories=freq_aranged, ordered = True)\n",
    "    baseline_matrix.loc[:, 'freq_target'] = pd.Categorical(baseline_matrix['freq_target'], \n",
    "                                                    categories=freq_aranged, ordered = True)\n",
    "    baseline_matrix = baseline_matrix.pivot_table(index = ['freq_source', 'sub_source', 'neuro_source'], \n",
    "                                        columns = ['freq_target', 'sub_target', 'neuro_target'], \n",
    "                                        values = 'ste')\n",
    "\n",
    "    # Full matrix\n",
    "    fig, axs = plt.subplots(1,1)\n",
    "    fig.set_size_inches(7.25, 7.25)\n",
    "    sns.set(font = 'sans-serif', font_scale=0.1) # you can create the labels here!\n",
    "    sns.heatmap(baseline_matrix, cmap='Blues', square=True, linewidths=.1, ax = axs, cbar_kws = {'shrink': 0.5}, vmax = vmax_full)\n",
    "    plt.hlines(range(12, 120, 12), 120, 0, linewidth = 0.25, linestyle = '--')\n",
    "    plt.hlines(range(24, 120, 24), 120, 0, linewidth = 0.5)\n",
    "    plt.vlines(range(12, 120, 12), 120, 0, linewidth = 0.25, linestyle = '--')\n",
    "    plt.vlines(range(24, 120, 24), 120, 0, linewidth = 0.5)\n",
    "    sns.despine(fig = fig, ax = axs)\n",
    "    # use matplotlib.colorbar.Colorbar object\n",
    "    cbar = axs.collections[0].colorbar\n",
    "    # here set the labelsize by 20\n",
    "    cbar.ax.tick_params(labelsize=4, axis='both', which = 'both', direction = 'in', pad = 0.85)\n",
    "    axs.spines['top'].set_visible(False)\n",
    "    axs.spines['right'].set_visible(False)\n",
    "    axs.get_xaxis().tick_bottom()\n",
    "    axs.get_yaxis().tick_left()\n",
    "    axs.xaxis.set_ticks_position('none')\n",
    "    axs.yaxis.set_ticks_position('none')\n",
    "    axs.tick_params(axis='both', which = 'both', direction = 'in', pad = 0.5)\n",
    "    fig.savefig('results/6.fixed_betweenVSbaseline_scrambled/baseline_avg_fixed_delay'+delays_seconds[i_delay]+'.eps', \n",
    "                        bbox_inches = 'tight',pad_inches = 0.1, dpi=300)\n",
    "\n",
    "    # Between\n",
    "    baseline_between = baseline_df.copy()\n",
    "    baseline_between.loc[baseline_between['sub_source'] == baseline_between['sub_target'], 'ste'] = 0.0\n",
    "    baseline_between = baseline_between.groupby(['neuro_source', 'neuro_target', 'sub_source', 'sub_target', 'freq_source', 'freq_target', 'pair'], \n",
    "                                as_index=False)['ste'].mean()\n",
    "    baseline_between = baseline_between.groupby(['neuro_source', 'neuro_target', 'sub_source', 'sub_target', 'freq_source', 'freq_target'], \n",
    "                                as_index=False)['ste'].mean()\n",
    "\n",
    "\n",
    "\n",
    "    # Consider these variables as categorical so they are ordered correctly when creating the figure\n",
    "    baseline_between.loc[:, 'neuro_source'] = pd.Categorical(baseline_between['neuro_source'], \n",
    "                                                     categories=sources_aranged, ordered = True)\n",
    "    baseline_between.loc[:, 'neuro_target'] = pd.Categorical(baseline_between['neuro_target'], \n",
    "                                                     categories=sources_aranged, ordered = True)\n",
    "    baseline_between.loc[:, 'freq_source'] = pd.Categorical(baseline_between['freq_source'], \n",
    "                                                    categories=freq_aranged, ordered = True)\n",
    "    baseline_between.loc[:, 'freq_target'] = pd.Categorical(baseline_between['freq_target'], \n",
    "                                                    categories=freq_aranged, ordered = True)\n",
    "    baseline_between = baseline_between.pivot_table(index = ['freq_source', 'sub_source', 'neuro_source'], \n",
    "                                        columns = ['freq_target', 'sub_target', 'neuro_target'], \n",
    "                                        values = 'ste')\n",
    "\n",
    "    fig, axs = plt.subplots(1,1)\n",
    "    fig.set_size_inches(7.25, 7.25)\n",
    "    sns.set(font = 'sans-serif', font_scale=0.1) # you can create the labels here!\n",
    "    sns.heatmap(baseline_between, cmap='Blues', square=True, linewidths=.1, ax = axs, cbar_kws = {'shrink': 0.5}, vmax = vmax_between)\n",
    "    plt.hlines(range(12, 120, 12), 120, 0, linewidth = 0.25, linestyle = '--')\n",
    "    plt.hlines(range(24, 120, 24), 120, 0, linewidth = 0.5)\n",
    "    plt.vlines(range(12, 120, 12), 120, 0, linewidth = 0.25, linestyle = '--')\n",
    "    plt.vlines(range(24, 120, 24), 120, 0, linewidth = 0.5)\n",
    "    sns.despine(fig = fig, ax = axs)\n",
    "    # use matplotlib.colorbar.Colorbar object\n",
    "    cbar = axs.collections[0].colorbar\n",
    "    # here set the labelsize by 20\n",
    "    cbar.ax.tick_params(labelsize=4, axis='both', which = 'both', direction = 'in', pad = 0.85)\n",
    "    axs.spines['top'].set_visible(False)\n",
    "    axs.spines['right'].set_visible(False)\n",
    "    axs.get_xaxis().tick_bottom()\n",
    "    axs.get_yaxis().tick_left()\n",
    "    axs.xaxis.set_ticks_position('none')\n",
    "    axs.yaxis.set_ticks_position('none')\n",
    "    axs.tick_params(axis='both', which = 'both', direction = 'in', pad = 0.5)\n",
    "    fig.savefig('results/6.fixed_betweenVSbaseline_scrambled/baseline_avg_between_fixed_delay'+delays_seconds[i_delay]+'.eps', \n",
    "                        bbox_inches = 'tight',pad_inches = 0.1, dpi=300)\n",
    "\n",
    "    # scrambled\n",
    "    # full\n",
    "    scrambled_matrix = scrambled_df.groupby(['neuro_source', 'neuro_target', 'info_source', 'info_target', 'freq_source', 'freq_target', 'pair_a'], \n",
    "                                as_index=False)['ste'].mean()\n",
    "    scrambled_matrix = scrambled_matrix.groupby(['neuro_source', 'neuro_target', 'info_source', 'info_target', 'freq_source', 'freq_target'], \n",
    "                                as_index=False)['ste'].mean()\n",
    "\n",
    "    # Consider these variables as categorical so they are ordered correctly when creating the figure\n",
    "    scrambled_matrix.loc[:, 'neuro_source'] = pd.Categorical(scrambled_matrix['neuro_source'], \n",
    "                                                     categories=sources_aranged, ordered = True)\n",
    "    scrambled_matrix.loc[:, 'neuro_target'] = pd.Categorical(scrambled_matrix['neuro_target'], \n",
    "                                                     categories=sources_aranged, ordered = True)\n",
    "    scrambled_matrix.loc[:, 'freq_source'] = pd.Categorical(scrambled_matrix['freq_source'], \n",
    "                                                    categories=freq_aranged, ordered = True)\n",
    "    scrambled_matrix.loc[:, 'freq_target'] = pd.Categorical(scrambled_matrix['freq_target'], \n",
    "                                                    categories=freq_aranged, ordered = True)\n",
    "\n",
    "    scrambled_matrix = scrambled_matrix.pivot_table(index = ['freq_source', 'info_source', 'neuro_source'], \n",
    "                                        columns = ['freq_target', 'info_target', 'neuro_target'], \n",
    "                                        values = 'ste')\n",
    "\n",
    "\n",
    "    # Full matrix\n",
    "    fig, axs = plt.subplots(1,1)\n",
    "    fig.set_size_inches(7.25, 7.25)\n",
    "    sns.set(font = 'sans-serif', font_scale=0.1) # you can create the labels here!\n",
    "    sns.heatmap(scrambled_matrix, cmap='Blues', square=True, linewidths=.1, ax = axs, cbar_kws = {'shrink': 0.5}, vmax = vmax_full)\n",
    "    plt.hlines(range(12, 120, 12), 120, 0, linewidth = 0.25, linestyle = '--')\n",
    "    plt.hlines(range(24, 120, 24), 120, 0, linewidth = 0.5)\n",
    "    plt.vlines(range(12, 120, 12), 120, 0, linewidth = 0.25, linestyle = '--')\n",
    "    plt.vlines(range(24, 120, 24), 120, 0, linewidth = 0.5)\n",
    "    sns.despine(fig = fig, ax = axs)\n",
    "    # use matplotlib.colorbar.Colorbar object\n",
    "    cbar = axs.collections[0].colorbar\n",
    "    # here set the labelsize by 20\n",
    "    cbar.ax.tick_params(labelsize=4, axis='both', which = 'both', direction = 'in', pad = 0.85)\n",
    "    axs.spines['top'].set_visible(False)\n",
    "    axs.spines['right'].set_visible(False)\n",
    "    axs.get_xaxis().tick_bottom()\n",
    "    axs.get_yaxis().tick_left()\n",
    "    axs.xaxis.set_ticks_position('none')\n",
    "    axs.yaxis.set_ticks_position('none')\n",
    "    axs.tick_params(axis='both', which = 'both', direction = 'in', pad = 0.5)\n",
    "    fig.savefig('results/6.fixed_betweenVSbaseline_scrambled/scrambled_avg_fixed_delay'+delays_seconds[i_delay]+'.eps', \n",
    "                        bbox_inches = 'tight',pad_inches = 0.1, dpi=300)\n",
    "    \n",
    "    # Between\n",
    "    scrambled_between = scrambled_df.copy()\n",
    "    scrambled_between.loc[scrambled_between['info_source'] == scrambled_between['info_target'], 'ste'] = 0.0\n",
    "    scrambled_between = scrambled_between.groupby(['neuro_source', 'neuro_target', 'info_source', 'info_target', 'freq_source', 'freq_target', 'pair_a'], \n",
    "                                as_index=False)['ste'].mean()\n",
    "    scrambled_between = scrambled_between.groupby(['neuro_source', 'neuro_target', 'info_source', 'info_target', 'freq_source', 'freq_target'], \n",
    "                                as_index=False)['ste'].mean()\n",
    "\n",
    "    scrambled_between.loc[:, 'neuro_source'] = pd.Categorical(scrambled_between['neuro_source'], \n",
    "                                                     categories=sources_aranged, ordered = True)\n",
    "    scrambled_between.loc[:, 'neuro_target'] = pd.Categorical(scrambled_between['neuro_target'], \n",
    "                                                     categories=sources_aranged, ordered = True)\n",
    "    scrambled_between.loc[:, 'freq_source'] = pd.Categorical(scrambled_between['freq_source'], \n",
    "                                                    categories=freq_aranged, ordered = True)\n",
    "    scrambled_between.loc[:, 'freq_target'] = pd.Categorical(scrambled_between['freq_target'], \n",
    "                                                    categories=freq_aranged, ordered = True)\n",
    "\n",
    "    scrambled_between = scrambled_between.pivot_table(index = ['freq_source', 'info_source', 'neuro_source'], \n",
    "                                        columns = ['freq_target', 'info_target', 'neuro_target'], \n",
    "                                        values = 'ste')\n",
    "\n",
    "    fig, axs = plt.subplots(1,1)\n",
    "    fig.set_size_inches(7.25, 7.25)\n",
    "    sns.set(font = 'sans-serif', font_scale=0.1) # you can create the labels here!\n",
    "    \n",
    "    sns.heatmap(scrambled_between, cmap='Blues', square=True, linewidths=.1, ax = axs, cbar_kws = {'shrink': 0.5}, vmax = vmax_between)\n",
    "    plt.hlines(range(12, 120, 12), 120, 0, linewidth = 0.25, linestyle = '--')\n",
    "    plt.hlines(range(24, 120, 24), 120, 0, linewidth = 0.5)\n",
    "    plt.vlines(range(12, 120, 12), 120, 0, linewidth = 0.25, linestyle = '--')\n",
    "    plt.vlines(range(24, 120, 24), 120, 0, linewidth = 0.5)\n",
    "    sns.despine(fig = fig, ax = axs)\n",
    "    # use matplotlib.colorbar.Colorbar object\n",
    "    cbar = axs.collections[0].colorbar\n",
    "    # here set the labelsize by 20\n",
    "    cbar.ax.tick_params(labelsize=4, axis='both', which = 'both', direction = 'in', pad = 0.85)\n",
    "    axs.spines['top'].set_visible(False)\n",
    "    axs.spines['right'].set_visible(False)\n",
    "    axs.get_xaxis().tick_bottom()\n",
    "    axs.get_yaxis().tick_left()\n",
    "    axs.xaxis.set_ticks_position('none')\n",
    "    axs.yaxis.set_ticks_position('none')\n",
    "    axs.tick_params(axis='both', which = 'both', direction = 'in', pad = 0.5)\n",
    "    fig.savefig('results/6.fixed_betweenVSbaseline_scrambled/scrambled_avg_between_fixed_delay'+delays_seconds[i_delay]+'.eps', \n",
    "                        bbox_inches = 'tight',pad_inches = 0.1, dpi=300)\n",
    "\n",
    "\n",
    "    plt.close('all')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis 4: Are there differences in terms of graph theory between the different pieces?\n",
    "\n",
    "We average pieces and compare the trial based graph theory statistics between homophonic and polyphonic pieces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 362,
   "metadata": {},
   "outputs": [],
   "source": [
    "delays = ['3', '30', '150'] # Only looking at 20.81ms, 200ms, 1000ms\n",
    "delays_seconds = ['20ms', '200ms', '1000ms']\n",
    "pairs = ['P03', 'P04', 'P05', 'P08', 'P09', 'P11']\n",
    "pieces = ['mel', 'can', 'clo', 'val']\n",
    "duo_type = ['h', 'p']\n",
    "\n",
    "# Loop through all the delays\n",
    "for i_delay in range(len(delays)):\n",
    "    playing_df = pd.read_csv('results/5.fixed_delay_trial/grand_csv_fixed_delay'+delays_seconds[i_delay]+'.csv', index_col = 0)\n",
    "    grand_df = playing_df.groupby(['neuro_source', 'neuro_target', 'info_source', 'info_target', 'freq_source', 'freq_target', 'duo_type', \n",
    "                                         'piece', 'pair', 'trial'], as_index=False)['ste'].mean()  \n",
    "    # Load the pmpq file -> This will help us and save us from doing weird shit down the line\n",
    "    graph_df_original = pd.read_csv('data/pmpq/pmpq.csv')\n",
    "    graph_df = graph_df_original[graph_df_original['pair'] != 'P01'].copy()\n",
    "    graph_df = graph_df[graph_df['pair'] != 'P02']\n",
    "    graph_df = graph_df[graph_df['pair'] != 'P06']\n",
    "    graph_df = graph_df[graph_df['scale'] != 'anxiety']\n",
    "    graph_df = graph_df.groupby(['pair', 'trial', 'piece', 'scale'], as_index=False)['score'].mean()\n",
    "    graph_df['duo_type'] = 'trash_kween'\n",
    "    graph_df.loc[:, 'duo_type'] = [graph_df.loc[:, 'piece'][i][0] for i in range(len(graph_df))]\n",
    "    graph_df.loc[:, 'piece'] = [graph_df.loc[:, 'piece'][i][1:] for i in range(len(graph_df))]\n",
    "\n",
    "    for i_pair in range(len(pairs)):\n",
    "        for i_piece in range(len(pieces)):\n",
    "            for i_trial in range(5):\n",
    "                if pairs[i_pair] == 'P05' and i_trial == 3 and pieces[i_piece] == 'can':\n",
    "                    continue\n",
    "                elif pairs[i_pair] == 'P05' and i_trial == 4 and pieces[i_piece] == 'can':\n",
    "                    continue\n",
    "\n",
    "                grand_df_piece = grand_df[grand_df['piece'] == pieces[i_piece]]\n",
    "                grand_df_piece = grand_df_piece[grand_df_piece['pair'] == pairs[i_pair]]\n",
    "                current_trial = grand_df_piece[grand_df_piece['trial'] == (i_trial+1)].copy()\n",
    "                current_trial = current_trial.pivot_table(index = ['freq_source', 'info_source', 'neuro_source'], \n",
    "                                                            columns = ['freq_target', 'info_target', 'neuro_target'], \n",
    "                                                            values = 'ste')\n",
    "                current_trial = current_trial.values\n",
    "                if pieces[i_piece] == 'mel' or pieces[i_piece] == 'val':\n",
    "                    duo_type = 'h'\n",
    "                elif pieces[i_piece] == 'clo' or pieces[i_piece] == 'can':\n",
    "                    duo_type = 'p'\n",
    "                else:\n",
    "                    duo_type = 'whut'\n",
    "\n",
    "                # First, we get the connection length matrix -- we need this one for some calculations downstream\n",
    "                current_trial_length = bct.weight_conversion(current_trial, 'lengths', copy=True) # this will get us the connection length matrix needed for some of the above functions! \n",
    "\n",
    "                # Clustering\n",
    "                avg_clustering_coef = bct.clustering_coef_wd(current_trial).mean() #  The weighted clustering coefficient is the average “intensity” of triangles around a node -> We take the average! \n",
    "                graph_df = graph_df.append({'pair': pairs[i_pair], 'trial': i_trial+1, 'piece': pieces[i_piece], 'scale': 'avg_clustering_coeff', 'score': avg_clustering_coef, 'duo_type': duo_type}, ignore_index=True)\n",
    "\n",
    "                transitivity = bct.transitivity_wd(current_trial) # Transitivity is the ratio of ‘triangles to triplets’ in the network. (A classical version of the clustering coefficient)\n",
    "                graph_df = graph_df.append({'pair': pairs[i_pair], 'trial': i_trial+1, 'piece': pieces[i_piece], 'scale': 'transitivity', 'score': transitivity, 'duo_type': duo_type}, ignore_index=True)\n",
    "\n",
    "                # Core\n",
    "                assortativity = bct.assortativity_wei(current_trial) # The assortativity coefficient is a correlation coefficient between the strengths (weighted degrees) of all nodes on two opposite ends of a link. A positive assortativity coefficient indicates that nodes tend to link to other nodes with the same or similar strength.\n",
    "                graph_df = graph_df.append({'pair': pairs[i_pair], 'trial': i_trial+1, 'piece': pieces[i_piece], 'scale': 'assortativity', 'score': assortativity, 'duo_type': duo_type}, ignore_index=True)\n",
    "\n",
    "                # Degree\n",
    "                avg_strengths = bct.strengths_dir(current_trial).mean() # Node strength is the sum of weights of links connected to the node. The instrength is the sum of inward link weights and the outstrength is the sum of outward link weights.\n",
    "                graph_df = graph_df.append({'pair': pairs[i_pair], 'trial': i_trial+1, 'piece': pieces[i_piece], 'scale': 'avg_node_strengths', 'score': avg_strengths, 'duo_type': duo_type}, ignore_index=True)\n",
    "\n",
    "                # Distance\n",
    "                # We first get the distance matrix to then get the characteristic path length\n",
    "                current_trial_distance, _ = bct.distance_wei(current_trial_length) #The distance matrix contains lengths of shortest paths between all pairs of nodes. An entry (u,v) represents the length of shortest path from node u to node v. The average shortest path length is the characteristic path length of the network.\n",
    "                char_path_length, efficiency, _, _, _ = bct.charpath(current_trial_distance, include_diagonal=False) # The characteristic path length is the average shortest path length in the network. The global efficiency is the average inverse shortest path length in the network.\n",
    "                graph_df = graph_df.append({'pair': pairs[i_pair], 'trial': i_trial+1, 'piece': pieces[i_piece], 'scale': 'char_path_length', 'score': char_path_length, 'duo_type': duo_type}, ignore_index=True)\n",
    "                graph_df = graph_df.append({'pair': pairs[i_pair], 'trial': i_trial+1, 'piece': pieces[i_piece], 'scale': 'efficiency', 'score': efficiency, 'duo_type': duo_type}, ignore_index=True)\n",
    "\n",
    "                # now that we have all the statistics we want, we output the data frame just in case!\n",
    "                graph_df.to_csv('results/7.graph_theory_fixed_delay/fixed_graph_theory_pertrial'+delays_seconds[i_delay]+'.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we test for differences between polyphonic and homophonic duos at each delay using permutations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 938,
   "metadata": {},
   "outputs": [],
   "source": [
    "delays = ['3', '30', '150'] # Only looking at 20.81ms, 200ms, 1000ms\n",
    "delays_seconds = ['20ms', '200ms', '1000ms']\n",
    "graph_stats = ['avg_clustering_coeff', 'avg_node_strengths', 'char_path_length', 'efficiency']\n",
    "\n",
    "for i_delay in range(len(delays)):\n",
    "    graph_df = pd.read_csv('results/7.graph_theory_fixed_delay/fixed_graph_theory_pertrial'+delays_seconds[i_delay]+'.csv', index_col= 0)\n",
    "\n",
    "    # create output dataframe with all data\n",
    "    out_data = {'graph_stat': graph_stats, 'score_h': np.zeros(len(graph_stats)), \n",
    "               'score_p': np.zeros(len(graph_stats)), 'p_value': np.zeros(len(graph_stats)), \n",
    "               'p_value_fdr': np.zeros(len(graph_stats)), 't_test': np.zeros(len(graph_stats))}\n",
    "    out_df = pd.DataFrame(data = out_data)\n",
    "\n",
    "    for graph_stat in graph_stats:\n",
    "        perm_cond = graph_df[graph_df['scale'] == graph_stat].copy()\n",
    "        perm_cond = perm_cond.groupby(['pair', 'duo_type'], as_index=False)['score'].mean()\n",
    "        perm_h = perm_cond.loc[perm_cond['duo_type'] == 'h', 'score'].values[:, None]\n",
    "        perm_p = perm_cond.loc[perm_cond['duo_type'] == 'p', 'score'].values[:, None]            \n",
    "        perm_graph = np.concatenate((perm_p, perm_h), axis = 1)\n",
    "        experimental_value, _ = sp.stats.ttest_rel(perm_graph[:, 0], perm_graph[:, 1]) # experimental value of ind t-test\n",
    "        perm_distribution = np.array([monte_carlo_dist(perm_graph) for i in range(500)])\n",
    "        perm_distribution = np.unique(perm_distribution)\n",
    "        p_value = (np.abs(experimental_value) < np.abs(perm_distribution)).sum()/64. # this is a two-sided test    \n",
    "        out_df.loc[out_df[out_df['graph_stat'] == graph_stat].index[0], 'p_value'] = p_value\n",
    "        out_df.loc[out_df[out_df['graph_stat'] == graph_stat].index[0], 't_test'] = experimental_value\n",
    "        out_df.loc[out_df[out_df['graph_stat'] == graph_stat].index[0], 'score_h'] = perm_h.mean()\n",
    "        out_df.loc[out_df[out_df['graph_stat'] == graph_stat].index[0], 'score_p'] = perm_p.mean()\n",
    "\n",
    "    # fdr to correct for multiple comparisons\n",
    "    reject_fdr, pval_fdr = fdr_correction(out_df['p_value'].values, alpha = 0.05, method = 'indep')\n",
    "    out_df['p_value_fdr'] = pval_fdr     \n",
    "    out_df.to_csv('results/7.graph_theory_fixed_delay/graph_theory_stats'+delays_seconds[i_delay]+'.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 942,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close('all')\n",
    "# Now we plot this graph theory statistics on their own \n",
    "titles = ['a', 'b', 'c', 'd']\n",
    "for i_delay in range(len(delays_seconds)):\n",
    "    graph_df = pd.read_csv('results/7.graph_theory_fixed_delay/fixed_graph_theory_pertrial'+delays_seconds[i_delay]+'.csv', index_col= 0)\n",
    "    graph_statistics = ['avg_clustering_coeff', 'avg_node_strengths', 'char_path_length', 'efficiency']\n",
    "    graph_title = ['Average Clustering Coefficient', 'Average Node Strength', 'Characteristic Path Length', 'Efficiency']\n",
    "    graph_df = graph_df[graph_df.scale.isin(graph_statistics)].reset_index()\n",
    "\n",
    "    colors_palette = sns.cubehelix_palette(3, start=1, rot=.7)\n",
    "    fig, axs = plt.subplots(4,1)\n",
    "    fig.set_size_inches(3.54, 7.25)\n",
    "    fig.subplots_adjust(hspace = 1.03, wspace = 0.4)\n",
    "    axs = axs.ravel()\n",
    "\n",
    "    for i_graph in range(len(graph_statistics)):\n",
    "        sns.violinplot(x=\"duo_type\", y=\"score\", data=graph_df[graph_df['scale'] == graph_statistics[i_graph]], \n",
    "                       palette = colors_palette, ax = axs[i_graph])\n",
    "        axs[i_graph].spines['top'].set_visible(False)\n",
    "        axs[i_graph].spines['right'].set_visible(False)\n",
    "        axs[i_graph].get_xaxis().tick_bottom()\n",
    "        axs[i_graph].get_yaxis().tick_left()\n",
    "        #axs[i_graph].xaxis.set_ticks_position('none')\n",
    "        #axs[i_graph].yaxis.set_ticks_position('none')\n",
    "        axs[i_graph].tick_params(axis='both', which = 'both', length = 0.2, width = 0.2)\n",
    "        axs[i_graph].set_title('%s'%titles[i_graph], loc = 'left', fontname = 'sans-serif', fontsize = 15)\n",
    "        axs[i_graph].set_ylabel(graph_title[i_graph])\n",
    "        axs[i_graph].set_xticklabels([\"Homophonic\", \"Polyphonic\"])\n",
    "        axs[i_graph].xaxis.labelpad = 0.4\n",
    "        axs[i_graph].yaxis.labelpad = 0.8\n",
    "        h_where = graph_df.loc[graph_df['scale'] == graph_statistics[i_graph], 'score'].values.max()*1.3\n",
    "        axs[i_graph].text(.5, h_where*1.001, '*', fontsize=12)\n",
    "        axs[i_graph].hlines(h_where, 0, 1, linewidth =1)\n",
    "        sns.set(font = 'sans-serif', font_scale = 0.5, style = 'white') # you can create the labels here!\n",
    "    fig.savefig('results/7.graph_theory_fixed_delay/graph_stats_fixed_delay'+delays_seconds[i_delay]+'.eps', \n",
    "                                bbox_inches = 'tight',pad_inches = 0.1, dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The values are significant, so we plot them and calculate the correlations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 943,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We plot the correlations between this statistics now \n",
    "graph_statistics = ['avg_clustering_coeff', 'avg_node_strengths', 'char_path_length', 'efficiency']\n",
    "pmpq_scales = ['quality', 'synchrony', 'synergy']\n",
    "graph_title = ['Average Clustering Coefficient', 'Average Node Strength', 'Characteristic Path Length', 'Efficiency']\n",
    "pmpq_title = ['Quality', 'Synchrony', 'Synergy']\n",
    "colors_palette = sns.cubehelix_palette(2, start=1, rot=.7)\n",
    "\n",
    "for i_delay in range(len(delays_seconds)):\n",
    "    # wrangle the pmpq csv to make plotter easier\n",
    "    graph_df = pd.read_csv('results/7.graph_theory_fixed_delay/fixed_graph_theory_pertrial'+delays_seconds[i_delay]+'.csv', index_col = 0)\n",
    "    graph_df = graph_df[~graph_df.scale.isin(['anxiety'])]\n",
    "    graph_df = graph_df[~graph_df.scale.isin(['transitivity', 'assortativity'])].reset_index(drop = True)\n",
    "    graph_df.loc[:, 'conditions'] = graph_df.iloc[:,[0, 1, 2]].apply(lambda x: \",\".join(x.astype(str)), axis=1).tolist()\n",
    "    graph_df = graph_df[~graph_df['conditions'].isin(['P05,4,can', 'P05,5,can'])].reset_index(drop = True)\n",
    "    graph_df_correlations = graph_df.loc[:2, :].copy()\n",
    "    graph_df_correlations['graph_stats'] = 'place_holder'\n",
    "    graph_df_correlations['stat'] = 0.0\n",
    "\n",
    "    for i_scale in range(len(pmpq_scales)):\n",
    "        for i_stat in range(len(graph_statistics)):\n",
    "            temp_pmpq = graph_df[graph_df['scale'] == pmpq_scales[i_scale]].sort_values(by = ['pair', 'piece', 'trial']).reset_index(drop = True).copy()\n",
    "            temp_stat = graph_df[graph_df['scale'] == graph_statistics[i_stat]].sort_values(by = ['pair', 'piece', 'trial']).reset_index(drop = True).copy()\n",
    "            temp_pmpq.loc[:, 'graph_stats'] = temp_stat['scale'].tolist()\n",
    "            temp_pmpq.loc[:, 'stat'] = temp_stat['score'].tolist()\n",
    "\n",
    "            graph_df_correlations = graph_df_correlations.append(temp_pmpq)\n",
    "            graph_df_correlations = graph_df_correlations.reset_index(drop = True)\n",
    "            if i_scale == 0 and i_stat == 0:\n",
    "                graph_df_correlations = graph_df_correlations.drop(index = [0, 1, 2]).reset_index(drop = True)\n",
    "    graph_df_correlations.loc[:, 'conditions'] = graph_df_correlations.iloc[:,[3, 7]].apply(lambda x: \",\".join(x.astype(str)), axis=1).tolist()\n",
    "\n",
    "    plt.close('all')\n",
    "\n",
    "    graph_conditions = np.unique(graph_df_correlations.loc[:, 'conditions'].tolist())\n",
    "\n",
    "    out_data = {'duo_type': ['place_holder'], 'pmpq_scale': ['place_holder'], 'graph_stats': ['place_holder'], 'corr': 0.0, 'p_value': 0.0}\n",
    "    out_df = pd.DataFrame(out_data)\n",
    "\n",
    "    fig, axs = plt.subplots(3,4)\n",
    "    fig.set_size_inches(7.25, 5)\n",
    "    fig.subplots_adjust(hspace = .5, wspace = 0.5)\n",
    "    axs = axs.ravel()\n",
    "\n",
    "    i_graph = 0\n",
    "    for i_scales in range(len(pmpq_scales)):\n",
    "        for i_stats in range(len(graph_statistics)):\n",
    "            # Get the beatiful plot\n",
    "            graph_temp = graph_df_correlations[graph_df_correlations.conditions.isin([graph_conditions[i_graph]])].reset_index(drop = False).copy()\n",
    "            sns.scatterplot(x=\"score\", y=\"stat\", hue = 'duo_type', data=graph_temp, palette = colors_palette, ax = axs[i_graph], s= 12)\n",
    "            axs[i_graph].spines['top'].set_visible(False)\n",
    "            axs[i_graph].spines['right'].set_visible(False)\n",
    "            axs[i_graph].get_xaxis().tick_bottom()\n",
    "            axs[i_graph].get_yaxis().tick_left()\n",
    "            axs[i_graph].tick_params(axis='both', which = 'both', length = 0.2, width = 0.2, pad = .005)\n",
    "            axs[i_graph].set_ylabel(graph_title[i_stats])\n",
    "            axs[i_graph].set_xlabel(pmpq_title[i_scales])\n",
    "            sns.set(font = 'sans-serif', font_scale = 0.5, style = 'white') # you can create the labels here!\n",
    "            axs[i_graph].xaxis.labelpad = 0.3\n",
    "            axs[i_graph].yaxis.labelpad = 0.3\n",
    "            axs[i_graph].legend(loc='upper right', frameon=False, scatterpoints=1, fontsize=3, markerscale = 0.1, markerfirst = True)\n",
    "            if graph_temp['graph_stats'].tolist()[0] == 'avg_clustering_coeff' or graph_temp['graph_stats'].tolist()[0] == 'efficiency':\n",
    "                axs[i_graph].set_ylim(bottom = 0.0, top = 0.006) \n",
    "            i_graph += 1\n",
    "\n",
    "\n",
    "            # Get the sweet taste of defeat because we all know these correlations are not there\n",
    "            poly_pmpq = graph_temp[graph_temp['duo_type'] == 'p'].reset_index(drop = True)['score'].values\n",
    "            poly_graph = graph_temp[graph_temp['duo_type'] == 'p'].reset_index(drop = True)['stat'].values\n",
    "            poly_corr, poly_p = sp.stats.pearsonr(poly_pmpq, poly_graph)\n",
    "            out_df = out_df.append({'corr': poly_corr, 'duo_type': 'p', 'graph_stats': graph_temp.graph_stats.tolist()[0], 'p_value': poly_p, 'pmpq_scale': graph_temp.scale.tolist()[0]}, ignore_index = True)\n",
    "\n",
    "            homo_pmpq = graph_temp[graph_temp['duo_type'] == 'h'].reset_index(drop = True)['score'].values\n",
    "            homo_graph = graph_temp[graph_temp['duo_type'] == 'h'].reset_index(drop = True)['stat'].values\n",
    "            homo_corr, homo_p = sp.stats.pearsonr(homo_pmpq, homo_graph)\n",
    "            out_df = out_df.append({'corr': homo_corr, 'duo_type': 'h', 'graph_stats': graph_temp.graph_stats.tolist()[0], 'p_value': homo_p, 'pmpq_scale': graph_temp.scale.tolist()[0]}, ignore_index = True)\n",
    "\n",
    "    fig.savefig('results/7.graph_theory_fixed_delay/graph_stats_pmpq_corrs_fixed_delay'+delays_seconds[i_delay]+'.eps', \n",
    "                    bbox_inches = 'tight',pad_inches = 0.1, dpi=300)\n",
    "    out_df = out_df.drop(index = 0).reset_index(drop = True)\n",
    "    _, out_df['p_value_fdr'] = fdr_correction(out_df['p_value'].values, alpha = 0.05, method = 'indep') \n",
    "    out_df.to_csv('results/7.graph_theory_fixed_delay/correlations_pmpq_graph_fixed'+delays_seconds[i_delay]+'.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 956,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Surprisingly, some of these correlations ARE significant, so, I'm going to plot them on their own!\n",
    "graph_statistics = ['avg_clustering_coeff', 'avg_node_strengths', 'char_path_length', 'efficiency']\n",
    "pmpq_scales = ['quality']\n",
    "graph_title = ['Average Clustering Coefficient', 'Average Node Strength', 'Characteristic Path Length', 'Efficiency']\n",
    "pmpq_title = ['Quality']\n",
    "colors_palette = sns.cubehelix_palette(4, start=1, rot=.7)\n",
    "\n",
    "for i_delay in range(len(delays_seconds)):\n",
    "    # wrangle the pmpq csv to make plotter easier\n",
    "    graph_df = pd.read_csv('results/7.graph_theory_fixed_delay/fixed_graph_theory_pertrial'+delays_seconds[i_delay]+'.csv', index_col = 0)\n",
    "    graph_df = graph_df[~graph_df.scale.isin(['anxiety'])]\n",
    "    graph_df = graph_df[~graph_df.scale.isin(['transitivity', 'assortativity'])].reset_index(drop = True)\n",
    "    graph_df.loc[:, 'conditions'] = graph_df.iloc[:,[0, 1, 2]].apply(lambda x: \",\".join(x.astype(str)), axis=1).tolist()\n",
    "    graph_df = graph_df[~graph_df['conditions'].isin(['P05,4,can', 'P05,5,can'])].reset_index(drop = True)\n",
    "    graph_df_correlations = graph_df.loc[:2, :].copy()\n",
    "    graph_df_correlations['graph_stats'] = 'place_holder'\n",
    "    graph_df_correlations['stat'] = 0.0\n",
    "    \n",
    "    # only take the significant ones!\n",
    "    for i_scale in range(len(pmpq_scales)):\n",
    "        for i_stat in range(len(graph_statistics)):\n",
    "            temp_pmpq = graph_df[graph_df['scale'] == pmpq_scales[i_scale]].sort_values(by = ['pair', 'piece', 'trial']).reset_index(drop = True).copy()\n",
    "            temp_stat = graph_df[graph_df['scale'] == graph_statistics[i_stat]].sort_values(by = ['pair', 'piece', 'trial']).reset_index(drop = True).copy()\n",
    "            temp_pmpq.loc[:, 'graph_stats'] = temp_stat['scale'].tolist()\n",
    "            temp_pmpq.loc[:, 'stat'] = temp_stat['score'].tolist()\n",
    "\n",
    "            graph_df_correlations = graph_df_correlations.append(temp_pmpq)\n",
    "            graph_df_correlations = graph_df_correlations.reset_index(drop = True)\n",
    "            if i_scale == 0 and i_stat == 0:\n",
    "                graph_df_correlations = graph_df_correlations.drop(index = [0, 1, 2]).reset_index(drop = True)\n",
    "    graph_df_correlations.loc[:, 'conditions'] = graph_df_correlations.iloc[:,[3, 7]].apply(lambda x: \",\".join(x.astype(str)), axis=1).tolist()\n",
    "    graph_df_correlations = graph_df_correlations[graph_df_correlations['duo_type']== 'p']\n",
    "\n",
    "    plt.close('all')\n",
    "\n",
    "    graph_conditions = np.unique(graph_df_correlations.loc[:, 'conditions'].tolist())\n",
    "\n",
    "    \n",
    "\n",
    "    fig, axs = plt.subplots(4,1)\n",
    "    fig.set_size_inches(3.54, 7)\n",
    "    fig.subplots_adjust(hspace = .5, wspace = 0.5)\n",
    "    axs = axs.ravel()\n",
    "\n",
    "    i_graph = 0\n",
    "    for i_scales in range(len(pmpq_scales)):\n",
    "        for i_stats in range(len(graph_statistics)):\n",
    "            # Get the beatiful plot\n",
    "            graph_temp = graph_df_correlations[graph_df_correlations.conditions.isin([graph_conditions[i_graph]])].reset_index(drop = False).copy()\n",
    "            sns.regplot(x=\"score\", y=\"stat\", data=graph_temp, color = colors_palette[i_graph], ci = None, fit_reg = True, ax = axs[i_graph], scatter_kws = {'s':8}, line_kws = {'linewidth':1.5})\n",
    "            axs[i_graph].spines['top'].set_visible(False)\n",
    "            axs[i_graph].spines['right'].set_visible(False)\n",
    "            axs[i_graph].get_xaxis().tick_bottom()\n",
    "            axs[i_graph].get_yaxis().tick_left()\n",
    "            axs[i_graph].tick_params(axis='both', which = 'both', length = 0.2, width = 0.2, pad = .005)\n",
    "            axs[i_graph].set_ylabel(graph_title[i_stats])\n",
    "            axs[i_graph].set_xlabel(pmpq_title[i_scales])\n",
    "            sns.set(font = 'sans-serif', font_scale = 0.5, style = 'white') # you can create the labels here!\n",
    "            if graph_temp['graph_stats'].tolist()[0] == 'avg_clustering_coeff' or graph_temp['graph_stats'].tolist()[0] == 'efficiency':\n",
    "                axs[i_graph].set_ylim(bottom = 0.0, top = 0.005) \n",
    "            i_graph += 1\n",
    "\n",
    "\n",
    "    fig.savefig('results/7.graph_theory_fixed_delay/sign_graph_stats_pmpq_corrs_fixed_delay'+delays_seconds[i_delay]+'.eps', \n",
    "                    bbox_inches = 'tight',pad_inches = 0.1, dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis 4: Do the networks exhibit small world properties as a function of time?\n",
    "\n",
    "We first look at the MAQ results, both pre and post (how much people liked each other before and after)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 966,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First test for the difference in the maq\n",
    "# Just to see what happens, we get the slops of the change between first trial and last trial and correlate that with how much they liked each other\n",
    "\n",
    "# we first read the maq file and wrangle it\n",
    "maq_df = pd.read_csv('data/questionnaires/5.MAQ.csv')\n",
    "maq_df = maq_df[maq_df['pair'] != 'P01']\n",
    "maq_df = maq_df[maq_df['pair'] != 'P02']\n",
    "maq_df = maq_df[maq_df['pair'] != 'P06']\n",
    "maq_df = maq_df.reset_index(drop = True)\n",
    "maq_df = maq_df.drop(columns=['Timestamp', 'Email address', '1. I enjoyed playing music with my music partner',\n",
    "                             '2. I would like to play again with my music partner', \n",
    "                             '3. When I was the \"follower\", I had no trouble musically accompanying my music partner'])\n",
    "\n",
    "maq_df['score'] = maq_df.iloc[:, 3:6].values.mean(axis = 1)\n",
    "maq_df = maq_df.drop(columns = ['4. I would like to become friends with my music partner', '5. If my music partner needed help, I would help them', \n",
    "                      '6. I would trust my music partner with a secret'])\n",
    "maq_df = maq_df.groupby(by = ['time', 'pair']).mean()\n",
    "maq_df = maq_df.reset_index()\n",
    "\n",
    "# We check differences between pre and post maq\n",
    "pre_maq = maq_df.loc[maq_df.time == 'pre', 'score'].values\n",
    "post_maq = maq_df.loc[maq_df.time == 'post', 'score'].values\n",
    "pre_post_maq = np.concatenate((np.expand_dims(pre_maq, axis = 1), np.expand_dims(post_maq, axis = 1)), axis = 1)\n",
    "t_exp_maq, _ = sp.stats.ttest_rel(pre_post_maq[:, 1], pre_post_maq[:, 0]) # experimental value of ind t-test\n",
    "perm_distribution = np.array([monte_carlo_dist(pre_post_maq) for i in range(5000)])\n",
    "perm_distribution = np.unique(perm_distribution)\n",
    "p_value_maq = (t_exp_maq < perm_distribution).sum()/64.\n",
    "f=open(\"results/8.fixed_small_world/pre_post_maq_ttest.txt\",\"w+\")\n",
    "out = 'The t-test value pre and post maq is %f and the pvalue is %f'%(t_exp_maq, p_value_maq)\n",
    "f.write(out)\n",
    "f.close()\n",
    "\n",
    "# plot it gurl \n",
    "colors_palette = sns.cubehelix_palette(3, start=1, rot=.7)\n",
    "fig, axs = plt.subplots(1,1)\n",
    "fig.set_size_inches(3.54, 3.54) \n",
    "sns.set(font = 'sans-serif', font_scale = 1, style = 'white') # you can create the labels here!\n",
    "sns.pointplot(x=\"time\", y=\"score\", data=maq_df, order = ['pre', 'post'], color = colors_palette[0], ax = axs)\n",
    "axs.spines['top'].set_visible(False)\n",
    "axs.spines['right'].set_visible(False)\n",
    "axs.get_xaxis().tick_bottom()\n",
    "axs.get_yaxis().tick_left()\n",
    "axs.tick_params(axis='both', which = 'both', length = 0.2, width = 0.2, pad = .005)\n",
    "axs.set_ylabel('Music Affiliation Questionnaire')\n",
    "axs.set_xlabel('Time')\n",
    "axs.text(.5, 5.15*1.001, '*', fontsize=12)\n",
    "axs.hlines(5.15, 0, 1, linewidth =1)\n",
    "fig.savefig('results/8.fixed_small_world/maq_pre_post.eps', bbox_inches = 'tight',pad_inches = 0.1, dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1018,
   "metadata": {},
   "outputs": [],
   "source": [
    "delays_seconds = ['20ms', '200ms', '1000ms']\n",
    "titles = ['a', 'b', 'c', 'd', 'e', 'f']\n",
    "plt.close('all')\n",
    "delays_seconds = ['20ms', '200ms', '1000ms']\n",
    "d = {'delay': ['20ms', '200ms', '1000ms']*2, 'corr': np.zeros(6), 'p_value': np.zeros(6), 'p_value_fdr': np.zeros(6), 'scale': 'place_holder'}\n",
    "out_df = pd.DataFrame(data=d)\n",
    "i_need_vacation = [1, 2, 3]\n",
    "\n",
    "# Now we create the small world evolution graphs at each delay \n",
    "plt.close('all')\n",
    "plt.clf()\n",
    "colors_palette = sns.cubehelix_palette(6, start=1, rot=.7)\n",
    "fig, axs = plt.subplots(3,2)\n",
    "fig.set_size_inches(7.25, 7)\n",
    "fig.subplots_adjust(hspace = 0.4, wspace = 0.4)\n",
    "axs = axs.ravel()\n",
    "\n",
    "for i_delay in range(len(delays_seconds)):\n",
    "    # We first get the piece order by creating an agregate of pair + piece\n",
    "    play_order = pd.read_csv('data/play_order.csv')\n",
    "    play_order = play_order[play_order['pair'] != 'P01']\n",
    "    play_order = play_order[play_order['pair'] != 'P02']\n",
    "    play_order = play_order[play_order['pair'] != 'P06']\n",
    "    play_order = play_order.reset_index(drop = True)\n",
    "    play_order.loc[:, 'piece'] = [play_order.loc[:, 'piece'][i][1:] for i in range(len(play_order))]\n",
    "    piece_1 = play_order[play_order['order'] == 1]\n",
    "    piece_2 = play_order[play_order['order'] == 2]\n",
    "    piece_3 = play_order[play_order['order'] == 3]\n",
    "    piece_4 = play_order[play_order['order'] == 4]\n",
    "    piece_1 = piece_1.iloc[:,[0, 1]].apply(lambda x: \",\".join(x.astype(str)), axis=1).tolist()\n",
    "    piece_2 = piece_2.iloc[:,[0, 1]].apply(lambda x: \",\".join(x.astype(str)), axis=1).tolist()\n",
    "    piece_3 = piece_3.iloc[:,[0, 1]].apply(lambda x: \",\".join(x.astype(str)), axis=1).tolist()\n",
    "    piece_4 = piece_4.iloc[:,[0, 1]].apply(lambda x: \",\".join(x.astype(str)), axis=1).tolist()\n",
    "    pieces_order = [piece_1, piece_2, piece_3, piece_4]\n",
    "\n",
    "    # now we wrangle the graph theory stats to get both clustering coefficient and characteristic path length\n",
    "    graph_statistics = ['avg_clustering_coeff', 'char_path_length']\n",
    "    graph_df = pd.read_csv('results/7.graph_theory_fixed_delay/fixed_graph_theory_pertrial'+delays_seconds[i_delay]+'.csv', index_col = 0)\n",
    "    graph_df = graph_df[graph_df.scale.isin(graph_statistics)].reset_index(drop = True)\n",
    "    graph_df.loc[:, 'conditions'] = graph_df.iloc[:,[0, 2]].apply(lambda x: \",\".join(x.astype(str)), axis=1).tolist()\n",
    "    graph_df['order'] = 0\n",
    "    for i_time in range(len(pieces_order)):\n",
    "        graph_df.loc[graph_df['conditions'].isin(pieces_order[i_time]), 'order'] = (i_time+1)\n",
    "\n",
    "    # get the small world coefficient and test correlation against time    \n",
    "    cc = graph_df.loc[graph_df['scale'] == 'avg_clustering_coeff', 'score'].values\n",
    "    cp = graph_df.loc[graph_df['scale'] == 'char_path_length', 'score'].values\n",
    "    smoll_world = graph_df[graph_df['scale'] == 'avg_clustering_coeff'].copy().reset_index(drop = True)\n",
    "    smoll_world['scale'] = 'small_world'\n",
    "    smoll_world['score'] = cc/cp\n",
    "    \n",
    "    smoll_cor, smoll_p = sp.stats.pearsonr(smoll_world.order.values, smoll_world.score.values)\n",
    "    out_df.loc[i_delay, 'scale'] = 'small_world'\n",
    "    out_df.loc[i_delay, 'corr'] = smoll_cor\n",
    "    out_df.loc[i_delay, 'p_value'] = smoll_p\n",
    "    \n",
    "    # plot it!\n",
    "    sns.set(font = 'sans-serif', font_scale = 0.7, style = 'white') # you can create the labels here!\n",
    "    sns.regplot(x=\"order\", y=\"score\", data=smoll_world, color = colors_palette[2*i_delay], ax = axs[2*i_delay], ci = None, fit_reg = True, scatter_kws = {'s':8}, line_kws = {'linewidth':1.5})\n",
    "    axs[2*i_delay].spines['top'].set_visible(False)\n",
    "    axs[2*i_delay].spines['right'].set_visible(False)\n",
    "    axs[2*i_delay].get_xaxis().tick_bottom()\n",
    "    axs[2*i_delay].get_yaxis().tick_left()\n",
    "    axs[2*i_delay].tick_params(axis='both', which = 'both', length = 0.2, width = 0.2, pad = .005)\n",
    "    axs[2*i_delay].set_ylabel('Small World Coefficient')\n",
    "    axs[2*i_delay].set_ylim(bottom = 0.0, top = 0.000019) \n",
    "    axs[2*i_delay].set_xticks(np.arange(1, 5))\n",
    "    axs[2*i_delay].set_yticks(np.arange(0.000005, 0.000016, 0.000005))\n",
    "    axs[2*i_delay].set_xlabel('Experimental Block')    \n",
    "    axs[2*i_delay].set_title('%s'%titles[i_delay], loc = 'left', fontname = 'sans-serif', fontsize = 15)\n",
    "    \n",
    "    # now we get the changes of how  much they liked each other vs the small world correlation\n",
    "    slope_maq = pre_post_maq[:, 1] - pre_post_maq[:, 0]\n",
    "    first_trial = smoll_world[smoll_world.order == 1].groupby(['pair'])['score'].mean().values\n",
    "    last_trial = smoll_world[smoll_world.order == 4].groupby(['pair'])['score'].mean().values\n",
    "    first_last_trial = np.concatenate((np.expand_dims(first_trial, axis = 1), np.expand_dims(last_trial, axis = 1)), axis = 1)\n",
    "    slope_smoll = (first_last_trial[:, 1] - first_last_trial[:, 0]) / 3\n",
    "    d = {'slope_maq': slope_maq, 'slope_smoll': slope_smoll}\n",
    "    slopes_df = pd.DataFrame(data = d)\n",
    "    sns.regplot(x='slope_maq', y='slope_smoll', data = slopes_df, color = colors_palette[i_delay+i_need_vacation[i_delay]], ci = None, fit_reg = True, ax = axs[i_delay+i_need_vacation[i_delay]], scatter_kws = {'s':8}, line_kws = {'linewidth':1.5})\n",
    "    axs[i_delay+i_need_vacation[i_delay]].spines['top'].set_visible(False)\n",
    "    axs[i_delay+i_need_vacation[i_delay]].spines['right'].set_visible(False)\n",
    "    axs[i_delay+i_need_vacation[i_delay]].get_xaxis().tick_bottom()\n",
    "    axs[i_delay+i_need_vacation[i_delay]].get_yaxis().tick_left()\n",
    "    axs[i_delay+i_need_vacation[i_delay]].tick_params(axis='both', which = 'both', length = 0.2, width = 0.2, pad = .005)\n",
    "    axs[i_delay+i_need_vacation[i_delay]].set_ylabel('Small World Slope')\n",
    "    axs[i_delay+i_need_vacation[i_delay]].set_xlabel('MAQ Slope')\n",
    "    axs[i_delay+i_need_vacation[i_delay]].set_xticks(np.arange(0,2.1, .5))\n",
    "    axs[i_delay+i_need_vacation[i_delay]].set_title('%s'%titles[i_delay+3], loc = 'left', fontname = 'sans-serif', fontsize = 15)\n",
    "    axs[i_delay+i_need_vacation[i_delay]].set_ylim(bottom = slopes_df.slope_smoll.min()*2, top = 0.000002) \n",
    "    sns.set(font = 'sans-serif', font_scale = 0.5, style = 'white') # you can create the labels here!\n",
    "    \n",
    "    smoll_maq_cor, smoll_maq_p = sp.stats.pearsonr(slope_maq, slope_smoll)\n",
    "    out_df.loc[i_delay+3, 'scale'] = 'smoll_maq_corr'\n",
    "    out_df.loc[i_delay+3, 'corr'] = smoll_maq_cor\n",
    "    out_df.loc[i_delay+3, 'p_value'] = smoll_maq_p\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "_, out_df['p_value_fdr'] = fdr_correction(out_df['p_value'].values, alpha = 0.05, method = 'indep') \n",
    "out_df.to_csv('results/8.fixed_small_world/corr_time_smoll_world.csv')\n",
    "fig.savefig('results/8.fixed_small_world/small_world_vs_time.eps', bbox_inches = 'tight',pad_inches = 0.1, dpi=300)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}