learned_tuning_NREM_QW_separate.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import scipy.io\n",
    "from scipy.stats import ranksums\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import h5py\n",
    "from learned_tuning.learned_tuning import calculate_learned_tuning, calculate_place_field_fidelity_of_learned_tuning\n",
    "import pandas as pd\n",
    "\n",
    "\n",
    "data_dir = r'/home/kouroshmaboudi/Documents/Learned_tuning_Python/Datasets'\n",
    "sessions = [d for d in os.listdir(data_dir) if os.path.isdir(os.path.join(data_dir, d))]\n",
    "\n",
    "included_session_IDX = [x for x in range(17) if x not in (12, 13)]\n",
    "# included_session_IDX = [10]\n",
    "\n",
    "sessions = [sessions[i] for i in included_session_IDX]\n",
    "\n",
    "for session_idx, session_name in enumerate(sessions):\n",
    "\n",
    "    print(session_name)\n",
    "\n",
    "    session_dataset_path = os.path.join(data_dir, session_name)\n",
    "    session_number = included_session_IDX[session_idx]\n",
    "\n",
    "\n",
    "\n",
    "    #--------------------------------------------------------------------------------------------\n",
    "    # Load epochs information\n",
    "\n",
    "    filename = f'{session_name}.fileInfo_for_python.mat'\n",
    "    file_path = os.path.join(session_dataset_path, filename)\n",
    "\n",
    "    mat_file = scipy.io.loadmat(file_path)\n",
    "    session_info = mat_file[\"fileInfo\"]\n",
    "\n",
    "    epochs = session_info[\"behavior\"][0][0][0][0][\"time\"]\n",
    "\n",
    "\n",
    "    #---------------------------------------------------------------------------------------------\n",
    "    # Brain state detection results\n",
    "\n",
    "    if 0<=session_number<=4 or 6<=session_number<=10: # Grosmark and Giri datasets\n",
    "        filename = f'{session_name}.brainStateDetection_HMMtheta_EMG_SWS_SchmidtTrigger.mat'\n",
    "        file_path = os.path.join(session_dataset_path, filename)\n",
    "\n",
    "        mat_file = scipy.io.loadmat(file_path)\n",
    "        brainStates_bouts_label = mat_file['brainState']['bouts'][0][0][:, :-1]\n",
    "        bouts_start_end = brainStates_bouts_label[:, :-1]\n",
    "        bout_duration = bouts_start_end[:, 1] - bouts_start_end[:, 0] \n",
    "        bout_labels = brainStates_bouts_label[:, -1].astype(int)\n",
    "\n",
    "        brainStates_names = []\n",
    "        for i in range(4):\n",
    "            brainStates_names.append(mat_file['brainState']['names'][0][0][i][0][0])\n",
    "\n",
    "    else: # Miyawaki dataset\n",
    "        filename = f'{session_name}.fileInfo_for_python.mat'\n",
    "        file_path = os.path.join(session_dataset_path, filename)\n",
    "\n",
    "        mat_file = scipy.io.loadmat(file_path)\n",
    "        brainStates_bouts_label = mat_file['fileInfo']['brainStates'].item()\n",
    "\n",
    "        bouts_start_end = brainStates_bouts_label[:, :-1]\n",
    "        bout_duration = bouts_start_end[:, 1] - bouts_start_end[:, 0] \n",
    "\n",
    "        bout_labels = brainStates_bouts_label[:, -1].astype(int)\n",
    "\n",
    "        # swapping the 3s and 4s, because the 3s in Hiro's brain_state_df are QWAKE while in the other datasets they are active wake(WAKE)\n",
    "        bout_labels_temp = bout_labels.copy()\n",
    "        bout_labels_temp[bout_labels == 3] = 4\n",
    "        bout_labels_temp[bout_labels == 4] = 3\n",
    "        bout_labels = bout_labels_temp[:]\n",
    "        del bout_labels_temp\n",
    "\n",
    "    brainStates_names = ['NREM', 'REM', 'WAKE', 'QWAKE', 'Undetermined']\n",
    "\n",
    "    too_short_state_bout_index = np.where(bout_duration < 6)\n",
    "    bout_labels[too_short_state_bout_index] = 5\n",
    "\n",
    "\n",
    "\n",
    "    unique_state_labels = np.unique(bout_labels)\n",
    "    num_bouts = bouts_start_end.shape[0]\n",
    "\n",
    "    time_bins_centers = np.arange(0, bouts_start_end[-1,1]+1, 1)+0.5    \n",
    "    time_bins_brain_state = np.full(time_bins_centers.shape, np.nan)\n",
    "\n",
    "    for bout_idx, bout_timing in enumerate(bouts_start_end):\n",
    "        inside_bout_idx = np.logical_and(bout_timing[0] < time_bins_centers,  time_bins_centers < bout_timing[1]) \n",
    "        time_bins_brain_state[inside_bout_idx] = bout_labels[bout_idx]\n",
    "\n",
    "\n",
    "    chunk_size_in_hours = 1/60 # one minutes\n",
    "    chunk_size = int(chunk_size_in_hours*3600) # in seconds\n",
    "    num_chunks = len(time_bins_brain_state) // chunk_size\n",
    "    chunk_start_times = [i*chunk_size_in_hours for i in range(num_chunks)] # centers of the time chunks\n",
    "    \n",
    "    reshaped_time_bins_brain_state = time_bins_brain_state[:num_chunks*chunk_size].reshape(num_chunks, chunk_size)\n",
    "\n",
    "    bins = np.arange(1, unique_state_labels.shape[0]+2)\n",
    "    counts = np.apply_along_axis(lambda x:np.histogram(x, bins = bins)[0] , axis = 1, arr=reshaped_time_bins_brain_state)\n",
    "\n",
    "\n",
    "    # Normalize the counts to get percentages\n",
    "    percentages = (counts / counts.sum(axis=1, keepdims=True))*100\n",
    "\n",
    "    # Create a DataFrame with percentages for each category and each chunk\n",
    "    brain_state_df = pd.DataFrame(percentages, columns=brainStates_names)\n",
    "\n",
    "    \n",
    "   \n",
    "\n",
    "    #-------------------------------------------------------------------------------------------\n",
    "    # Loading the population burst evenst with all their corresponding measured lfp features\n",
    "\n",
    "    overwrite = False # in case we need to read the .mat file again, if there was a change\n",
    "\n",
    "    filename = f'{session_name}.PBEs.pkl'\n",
    "    file_path = os.path.join(session_dataset_path, filename)\n",
    "\n",
    "    if os.path.exists(file_path) and overwrite == False:\n",
    "\n",
    "        # PBEs = np.load(file_path, allow_pickle=True)\n",
    "        PBEs = pd.read_pickle(file_path)\n",
    "    else: # if it doesn't exist then read it from the .mat file\n",
    "\n",
    "        filename = f'{session_name}.PBEInfo_replayScores_with_spindle_and_deltaPowers.mat'\n",
    "        file_path = os.path.join(session_dataset_path, filename)\n",
    "        f = h5py.File(file_path, \"r\")\n",
    "        PBEInfo = f['PBEInfo_replayScores']    \n",
    "       \n",
    "\n",
    "        # Store the population burst events in a pandas DataFrame\n",
    "        num_PBEs = PBEInfo[\"fr_1msbin\"].shape[0]\n",
    "        attributes = list(PBEInfo.keys())\n",
    "\n",
    "        PBEs = pd.DataFrame(columns=[attr for attr in attributes[1:] if attr not in ['posteriorProbMat', 'postMat_nonNorm']])\n",
    "\n",
    "        # Loop over the PBEs and add each one as a row to the DataFrame\n",
    "        num_dots = int(num_PBEs * (10/100))\n",
    "        count = 0\n",
    "\n",
    "        for pbe in range(num_PBEs): #  \n",
    "  \n",
    "            for attr in PBEs.columns:\n",
    "                ref = PBEInfo[attr][pbe][0]\n",
    "                obj = f[ref]\n",
    "\n",
    "                if attr in ['epoch', 'brainState']: # convert the ascii code to string\n",
    "                    arr = np.array(obj).flatten()\n",
    "                    epoch = \"\".join(chr(code) for code in arr)\n",
    "                    PBEs.at[pbe, attr] = epoch\n",
    "                elif attr in ['fr_1msbin', 'fr_20msbin', 'posteriorProbMat', 'postMat_nonNorm']: # no need to flatten\n",
    "                    arr = np.array(obj)\n",
    "                    PBEs.at[pbe, attr] = arr\n",
    "                else: \n",
    "                    arr = np.array(obj).flatten()\n",
    "                    PBEs.at[pbe, attr] = arr\n",
    "\n",
    "        if (pbe+1) % num_dots == 1:\n",
    "            count += 1\n",
    "            message = \"Importing PBEs\" + \".\" * count\n",
    "            print(message, end=\"\\r\")\n",
    "\n",
    "        print(\"All PBEs were imported\") \n",
    "    \n",
    "        filename = f'{session_name}.PBEs.pkl'\n",
    "        file_path = os.path.join(session_dataset_path, filename)\n",
    "        PBEs.to_pickle(file_path)\n",
    "\n",
    "    num_PBEs = PBEs.shape[0]\n",
    "\n",
    "    # add a brain state label depending on the time bin of the PBE\n",
    "    for pbe in range(num_PBEs):\n",
    "        idx = np.where(np.logical_and(time_bins_centers-0.5 <= PBEs.at[pbe, 'peakT'],  time_bins_centers+0.5 > PBEs.at[pbe, 'peakT']))[0]\n",
    "        PBEs.at[pbe, 'brain_state'] = time_bins_brain_state[idx]\n",
    "\n",
    "\n",
    "\n",
    "    #----------------------------------------------------------------------------------------------\n",
    "    # Load spike data\n",
    "    filename = f'{session_name}.spikes_for_python.mat'\n",
    "    file_path = os.path.join(session_dataset_path, filename)\n",
    "\n",
    "    mat_file = scipy.io.loadmat(file_path)\n",
    "    spikes_pyr = mat_file[\"spikes_pyr\"]\n",
    "\n",
    "\n",
    "    # Load unit stability information\n",
    "\n",
    "    filename = f'{session_name}.cluster_quality_by_block'\n",
    "    file_path = os.path.join(session_dataset_path, filename)\n",
    "\n",
    "    mat_file = scipy.io.loadmat(file_path)\n",
    "    cluster_quality_by_block = mat_file['cluster_quality_by_block'][0]\n",
    "\n",
    "\n",
    "\n",
    "    #### Extracting all place fields from the imported .mat file\n",
    "    spatial_tuning_smoothed = spikes_pyr[\"spatialTuning_smoothed\"]\n",
    "\n",
    "    num_units    = spatial_tuning_smoothed[0].shape[0]\n",
    "    # num_units    = spatial_tuning_smoothed.shape[0] # for RatN only\n",
    "\n",
    "    num_pos_bins = spatial_tuning_smoothed[0][0]['uni'][0][0].size\n",
    "\n",
    "    # print(num_units, num_pos_bins)\n",
    "\n",
    "    spikes = []; # spike data and place field info of each unit\n",
    "\n",
    "    # attributes = list(spikes_pyr.dtype.names) % if we want to work on all variable in the imported .mat data structure\n",
    "    running_directions = {'LR', 'RL', 'uni'}\n",
    "    other_attributes   = {'spike_times', 'shank_id','cluster_id'}\n",
    "\n",
    "    iter = 0\n",
    "    for unit in range(num_units):\n",
    "        \n",
    "        # Create dictionaries for each unit and store the matrices\n",
    "        \n",
    "        unit_spikes = dict()\n",
    "        \n",
    "        unit_spikes['place_fields']  = {}\n",
    "        unit_spikes['peak_pos_bins'] = {}\n",
    "\n",
    "        \n",
    "        for direction in running_directions:\n",
    "            try:\n",
    "                if session_number in [6,7]:\n",
    "                    unit_spikes['place_fields'][direction] = spatial_tuning_smoothed[unit][0][direction][0][0].reshape(num_pos_bins) \n",
    "                    unit_spikes['peak_pos_bins'][direction] = spikes_pyr['peakPosBin'][unit][0][direction][0][0][0][0]\n",
    "                else:\n",
    "                    unit_spikes['place_fields'][direction] = spatial_tuning_smoothed[0][unit][direction][0][0].reshape(num_pos_bins) \n",
    "                    unit_spikes['peak_pos_bins'][direction] = spikes_pyr['peakPosBin'][0][unit][direction][0][0][0][0]\n",
    "\n",
    "\n",
    "            except ValueError:\n",
    "                if iter == 0:\n",
    "                    print(\"This session has only one running direction\")\n",
    "                iter += 1\n",
    "\n",
    "\n",
    "        if session_number in [9, 10]: # for Rat V sessions\n",
    "            unit_spikes['spike_times'] = spikes_pyr['time'][0][unit]\n",
    "            unit_spikes['shank_id']    = spikes_pyr['id'][0][unit][0][1]\n",
    "            unit_spikes['shank_id']    += 1\n",
    "            unit_spikes['cluster_id']  = spikes_pyr['id'][0][unit][0][0]\n",
    "\n",
    "        elif session_number in [6, 7]: # for RatN and RatS\n",
    "            unit_spikes['spike_times'] = spikes_pyr['time'][unit][0] \n",
    "            unit_spikes['shank_id']    = spikes_pyr['id'][unit][0][0][0]\n",
    "            unit_spikes['cluster_id']  = spikes_pyr['id'][unit][0][0][1]\n",
    "\n",
    "        elif session_number == 8: # RatU  \n",
    "            unit_spikes['spike_times'] = spikes_pyr['time'][0][unit]\n",
    "            unit_spikes['shank_id']    = spikes_pyr['id'][0][unit][0][0] # shank indices already starts at zero\n",
    "            unit_spikes['shank_id']    += 1\n",
    "            unit_spikes['cluster_id']  = spikes_pyr['id'][0][unit][0][1]\n",
    "        else: # Grosmark, Hiro, and all other sessions\n",
    "            unit_spikes['spike_times'] = spikes_pyr['time'][0][unit]\n",
    "            unit_spikes['shank_id']    = spikes_pyr['id'][0][unit][0][0] # need to go one down for the other datasets\n",
    "            unit_spikes['cluster_id']  = spikes_pyr['id'][0][unit][0][1]\n",
    "\n",
    "\n",
    "\n",
    "        spikes.append(unit_spikes) \n",
    "\n",
    "\n",
    "\n",
    "    # place fields by pooling spikes across both running directions\n",
    "    place_fields_uni = []\n",
    "    for unit in range(num_units):\n",
    "        place_fields_uni.append(spikes[unit]['place_fields']['uni'])\n",
    "    place_fields_uni = np.array(place_fields_uni)\n",
    "\n",
    "    place_fields_uni[place_fields_uni == 0] = 1e-4\n",
    "\n",
    "\n",
    "\n",
    "    # --------------------------------------------------------------------------------------------\n",
    "    # Load cluster quality data (L-ratios)\n",
    "    \n",
    "    filename = f'{session_name}.clusterQuality.mat'\n",
    "    file_path = os.path.join(session_dataset_path, filename)\n",
    "\n",
    "    mat_file = scipy.io.loadmat(file_path)\n",
    "\n",
    "    # Access data structure\n",
    "    num_shanks = len(mat_file[\"clusterQuality\"][\"Lratio\"][0])\n",
    "\n",
    "    L_ratios = list()\n",
    "    for shank in range(num_shanks):    \n",
    "        curr_shank_L_ratios = dict()\n",
    "        curr_shank_L_ratios[\"L_ratios\"] = mat_file[\"clusterQuality\"][\"Lratio\"][0][shank]\n",
    "        curr_shank_L_ratios[\"cluster_ids\"] = mat_file[\"clusterQuality\"][\"clus\"][0][shank]\n",
    "        \n",
    "        L_ratios.append(curr_shank_L_ratios)\n",
    "    # L_ratios = []\n",
    "\n",
    "\n",
    "    #------------------------------------------------------------------------------------------------------\n",
    "    # Learned tunings during Non-REM versus Quiet Wake ripples\n",
    "\n",
    "\n",
    "    # Calculate learned tunings separately for NREM and QW PBEs and calculate thier PF fidelities\n",
    "    time_bin_duration= 0.02\n",
    "    # num_PF_shuffles = 10000\n",
    "    learned_tunings_NREM_vs_QW = {}  # create an empty dictionary to store the results\n",
    "\n",
    "\n",
    "    for epoch in ['PRE', 'POST']:\n",
    "        if epoch == 'PRE':\n",
    "            epoch_idx = 0\n",
    "            epoch_duration = epochs[0,1] - epochs[0,0]\n",
    "        elif epoch == 'POST':\n",
    "            epoch_idx = 2\n",
    "            epoch_duration = 4*60*60\n",
    "            \n",
    "        # PBEs for the current epoch\n",
    "        \n",
    "        if_inside_epoch = PBEs['peakT'].between(epochs[epoch_idx,0], epochs[epoch_idx,0]+epoch_duration)\n",
    "\n",
    "        # NREM PBEs\n",
    "\n",
    "\n",
    "        \n",
    "        select_IDX = np.where(np.logical_and(if_inside_epoch, (PBEs['brain_state'] == 1)))[0]\n",
    "        epoch_NREM_PBEs = PBEs.loc[select_IDX].reset_index(drop=True)\n",
    "        num_PBEs_NREM = len(epoch_NREM_PBEs)\n",
    "\n",
    "\n",
    "        # QW PBEs\n",
    "        select_IDX = np.where(np.logical_and(if_inside_epoch, (PBEs['brain_state'] == 4)))[0]\n",
    "        epoch_QW_PBEs = PBEs.loc[select_IDX].reset_index(drop=True)\n",
    "        num_PBEs_QW = len(epoch_QW_PBEs)\n",
    "\n",
    "\n",
    "\n",
    "        # learned tunings\n",
    "        # NREM\n",
    "        learned_tunings_NREM = calculate_learned_tuning(epoch_NREM_PBEs, spikes, L_ratios, time_bin_duration)\n",
    "        \n",
    "        # learned_tuning_place_field_pearson_corr_NREM = np.full((num_units,), np.nan)\n",
    "\n",
    "        # learned_tuning_place_field_pearson_corr_NREM, _, median_LT_PF_pearson_corr_NREM = calculate_place_field_fidelity_of_learned_tuning(learned_tunings_NREM[active_units_epochs_intersect, :], place_fields_uni[active_units_epochs_intersect, :], num_PF_shuffles)\n",
    "        # The reason for why we are restricting the calculations to only active_units_epochs_intersect is for the shuffle procedure. We want to make sure that the only units with significant PFs are swapped for the each unit's PF\n",
    "\n",
    "        # QW\n",
    "        learned_tunings_QW = calculate_learned_tuning(epoch_QW_PBEs, spikes, L_ratios, time_bin_duration)\n",
    "\n",
    "        # learned_tuning_place_field_pearson_corr_QW = np.full((num_units,), np.nan)\n",
    "        # learned_tuning_place_field_pearson_corr_QW, _, median_LT_PF_pearson_corr_QW = calculate_place_field_fidelity_of_learned_tuning(learned_tunings_QW[active_units_epochs_intersect, :], place_fields_uni[active_units_epochs_intersect, :], num_PF_shuffles)\n",
    "        \n",
    "        # store the results in the dictionary\n",
    "        learned_tunings_NREM_vs_QW[epoch] = {\n",
    "            'learned_tunings_NREM': learned_tunings_NREM,\n",
    "            # 'learned_tuning_place_field_pearson_corr_NREM': learned_tuning_place_field_pearson_corr_NREM,\n",
    "            # 'median_LT_PF_pearson_corr_NREM': median_LT_PF_pearson_corr_NREM,\n",
    "            'number_of_PBEs_NREM':num_PBEs_NREM,\n",
    "            'learned_tunings_QW': learned_tunings_QW,\n",
    "            # 'learned_tuning_place_field_pearson_corr_QW': learned_tuning_place_field_pearson_corr_QW,\n",
    "            # 'median_LT_PF_pearson_corr_QW': median_LT_PF_pearson_corr_QW,\n",
    "            'number_of_PBEs_QW':num_PBEs_QW\n",
    "        }\n",
    "\n",
    "        filename = f'{session_name}.learned_tunings_NREM_vs_QW_new_brain_state.npy'\n",
    "        file_path = os.path.join(session_dataset_path, filename)\n",
    "        np.save(file_path, learned_tunings_NREM_vs_QW)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "    # # Plot the distributions of learned tunings for NREM and QW periods within each epoch\n",
    "\n",
    "    # colors = sns.color_palette(\"husl\", 2) # Set color palette\n",
    "    # sns.set_style('whitegrid') # Set style and context\n",
    "    # sns.set_context('paper')\n",
    "\n",
    "    # custom_params = {\"axes.spines.right\": False, \"axes.spines.top\": False}\n",
    "    # sns.set_theme(style=\"ticks\", rc=custom_params)\n",
    "\n",
    "    # fig, axes = plt.subplots(1,2, figsize = (10, 4))\n",
    "\n",
    "    # for i, curr_epoch in enumerate(['PRE', 'POST']):\n",
    "\n",
    "    #     learned_tuning_place_field_pearson_corr_NREM = learned_tunings_NREM_vs_QW[curr_epoch]['learned_tuning_place_field_pearson_corr_NREM']\n",
    "    #     median_LT_PF_pearson_corr_NREM_pvalue = learned_tunings_NREM_vs_QW[curr_epoch]['median_LT_PF_pearson_corr_NREM']['p_value']\n",
    "        \n",
    "    #     learned_tuning_place_field_pearson_corr_QW = learned_tunings_NREM_vs_QW[curr_epoch]['learned_tuning_place_field_pearson_corr_QW']\n",
    "    #     median_LT_PF_pearson_corr_QW_pvalue = learned_tunings_NREM_vs_QW[curr_epoch]['median_LT_PF_pearson_corr_QW']['p_value']\n",
    "\n",
    "    #     num_PBEs_NREM = learned_tunings_NREM_vs_QW[curr_epoch]['number_of_PBEs_NREM']\n",
    "    #     num_PBEs_QW = learned_tunings_NREM_vs_QW[curr_epoch]['number_of_PBEs_QW']\n",
    "\n",
    "\n",
    "    #     # Plot the distributions using ECDFs with overlaid ticks\n",
    "    #     sns.ecdfplot(learned_tuning_place_field_pearson_corr_NREM, ax = axes[i], label=f'{curr_epoch} NREM', color=colors[0], linewidth = 2)\n",
    "    #     sns.ecdfplot(learned_tuning_place_field_pearson_corr_QW, ax = axes[i], label=f'{curr_epoch} QW', color=colors[1], linewidth = 2)\n",
    "    #     axes[i].set_xlim([-1,1])\n",
    "    #     axes[i].set_xlabel(f'{curr_epoch} LT fidelity', fontsize=10)\n",
    "    #     axes[i].set_ylabel('Proportion of units', fontsize=10)\n",
    "    #     axes[i].tick_params(labelsize=8)\n",
    "    #     axes[i].legend(fontsize=8)\n",
    "\n",
    "\n",
    "    #     # Add ticks to the x-axis\n",
    "    #     y_min, y_max = axes[i].get_ylim()\n",
    "    #     for x in learned_tuning_place_field_pearson_corr_NREM:\n",
    "    #         axes[i].plot([x, x], [y_min+0.01, y_min + 0.02], '|-', color=colors[0], linewidth=0.25)\n",
    "    #     for x in learned_tuning_place_field_pearson_corr_QW:\n",
    "    #         axes[i].plot([x, x], [y_min+0.06, y_min + 0.07], '|-', color=colors[1], linewidth=0.25)\n",
    "\n",
    "    #     # Calculate and display medians\n",
    "    #     median_nrem = np.nanmedian(learned_tuning_place_field_pearson_corr_NREM)\n",
    "    #     median_qw = np.nanmedian(learned_tuning_place_field_pearson_corr_QW)\n",
    "\n",
    "\n",
    "    #     def get_pval_statement(pvalue):\n",
    "    #         if pvalue < 0.0001:\n",
    "    #             pvalue_statement = 'P<1e-4'\n",
    "    #         else:\n",
    "    #             pvalue_statement = f'P={pvalue:.4f}'\n",
    "    #         return pvalue_statement\n",
    "\n",
    "\n",
    "    #     axes[i].axvline(median_nrem, color=colors[0], linestyle='dashed', label=f'Median {curr_epoch} NREM={median_nrem:.2f},{get_pval_statement(median_LT_PF_pearson_corr_NREM_pvalue)}')\n",
    "    #     axes[i].axvline(median_qw, color=colors[1], linestyle='dashed', label=f'Median {curr_epoch} QW={median_qw:.2f},{get_pval_statement(median_LT_PF_pearson_corr_QW_pvalue)}')\n",
    "    #     axes[i].legend(fontsize=6)\n",
    "\n",
    "    #     # Perform rank-sum test\n",
    "    #     statistic, p_value = ranksums(learned_tuning_place_field_pearson_corr_NREM,\n",
    "    #                                 learned_tuning_place_field_pearson_corr_QW)\n",
    "\n",
    "    #     # Add line with p-value above the plot\n",
    "    #     axes[i].text((median_nrem + median_qw) / 2, 0.2, f'p-value = {p_value:.4f}', ha='center', fontsize=6)\n",
    "\n",
    "        \n",
    "\n",
    "    #     # Bar plot of frequency of ripples in each brain states\n",
    "    #     ax_inset = axes[i].inset_axes([0.2, 0.4, 0.2, 0.2])\n",
    "\n",
    "\n",
    "    #     ax_inset.bar([0, 1],[num_PBEs_NREM, num_PBEs_QW], color=[colors[0], colors[1]]) #unique_strings, \n",
    "    #     ax_inset.set_xticks([0, 1], ['NREM', 'QW'])\n",
    "    #     ax_inset.set_ylabel('number of PBEs', fontsize=6)\n",
    "    #     ax_inset.tick_params(labelsize=6)\n",
    "    #     # axes[i].tight_layout()\n",
    "\n",
    "\n",
    "    # filename = f'{session_name}.learned_tunings_NREM_vs_QW_with_cluster_quality_criterion.svg'\n",
    "    # file_path = os.path.join(session_dataset_path, filename)\n",
    "    # plt.savefig(file_path)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epoch_NREM_PBEs"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}