FOOOF_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Author: Jacob Woessner <woessner.jacob@gmail.com>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## HOW TO RUN\n",
    "- Set the path to the folder with the data in the variable `import_path`\n",
    "- Set project name in the variable `project_name`\n",
    "- Optional: Set home directory in the variable `home_dir` where files will be generated (default is the location of this notebook)\n",
    "- Hit `Run All` in the tab\n",
    "- The script automatically populates all the folders and files in the home directory\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>Note:</b> The script should skip files that have already been processed so if an error occurs you should be able to fix it and hit Run All without having to start from the beginning.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## INPUT:\n",
    "- Resting EEG files in the CNT format, file name should be in the format of a subject id at the beginning of the file name separated by an underscore (Ex: 1234_resting1.cnt)\n",
    "<div class=\"alert alert-block alert-warning\">\n",
    "<b>Warning:</b> The script concatenates Resting EEG files with the same subject id (Ex: 1234_resting1.cnt, 1234_resting2.cnt, 1234_resting3.cnt) so make sure the subject id is the same for all files and that eyes closed and eyes open files are not concatenated together.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## OUTPUT:\n",
    "- Log file with the name `<ENTER_NAME>.log` in the home directory\n",
    "- psd files in a psd folder in the home directory (json format)\n",
    "- individual fooof csv files in a fooof folder in the home directory\n",
    "- a fooof csv file with aperiodic parameters and error for all subjects in the home directory\n",
    "- a fooof func csv file with aperiodic, periodic and orginal function values\n",
    "- wide format spss file with fooof csv and fooof func csv combined\n",
    "- Analysis of the data (fooof plots, raincloud plots, and statistics)\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">\n",
    "<b>Warning:</b> As of January 1st, 2023, There may be a bug that appends empty data to cnt files in mne. I have included a conditional statement to check for this bug.\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mne\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "from fooof import FOOOF\n",
    "from fooof.sim.gen import gen_aperiodic\n",
    "import os\n",
    "import json\n",
    "import pyreadstat\n",
    "import datetime\n",
    "import sys\n",
    "import ptitprince as pt\n",
    "import seaborn as sns\n",
    "\n",
    "# configure backend here for plotting\n",
    "matplotlib.use('Agg')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## File paths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CNT import path\n",
    "import_cnt = 'PATH'\n",
    "\n",
    "#Project name string\n",
    "project_name = 'ENTER_NAME'\n",
    "\n",
    "# home directory: where this script is located\n",
    "home_dir = os.getcwd()\n",
    "\n",
    "# Set stdout and stderr used for logging print statements (comment out if you want to print to terminal)\n",
    "now = datetime.datetime.now()\n",
    "log_file = open(f'{home_dir}/log_{project_name}_{now.strftime(\"%Y-%m-%d_%H:%M\")}.txt', 'a')\n",
    "sys.stdout = log_file\n",
    "sys.stderr = log_file\n",
    "\n",
    "\n",
    "def make_dir_path(path):\n",
    "    os.makedirs(path, exist_ok=True)\n",
    "    return path\n",
    "\n",
    "# Make a directory called psd_{project_name} if it doesnt exist\n",
    "psd_path = make_dir_path(f'{home_dir}/psd_{project_name}')\n",
    "\n",
    "# Make a directory called fooof_{project_name} if it doesnt exist\n",
    "fooof_path = make_dir_path(f'{home_dir}/fooof_{project_name}')\n",
    "\n",
    "# Make a directory called analysis_{project_name} if it doesnt exist\n",
    "analysis_path = make_dir_path(f'{home_dir}/analysis_{project_name}')\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Processing Settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fooof fit frequency range in Hz\n",
    "freq_range = [1, 55]\n",
    "\n",
    "# fooof fit mode\n",
    "mode = 'fixed'\n",
    "\n",
    "# epoch duration in seconds\n",
    "epoch_dur = 2.048\n",
    "\n",
    "# Baseline correction\n",
    "baseline = (-0.2, 0)\n",
    "\n",
    "# psd method\n",
    "psd_method = 'welch'\n",
    "\n",
    "# Selected channels\n",
    "fooof_channels = 'all'\n",
    "\n",
    "# Bad channels to remove from fooof_channels\n",
    "bad_fooof_channels = ['X', 'Y', 'BLANK']\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis Settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select channels to be used for graphical analysis AS A LIST of lower case strings\n",
    "# e.g. ['fz', 'cz', 'pz']\n",
    "analysis_channels = ['fcz', 'cz', 'pz', 'f8', 'f7']\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate PSD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "#Get list of all files in directory\n",
    "file_list = os.listdir(import_cnt)\n",
    "print(file_list)\n",
    "#Get list of all files that end with .cnt\n",
    "subject_list = [i.split('_')[0] for i in file_list if i.endswith('.cnt')]\n",
    "cnt_list = [i for i in file_list if i.endswith('.cnt')]\n",
    "#Sort list of cnt files\n",
    "cnt_list.sort()\n",
    "subject_list = list(set(subject_list.copy()))\n",
    "subject_list.sort()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import_path = import_cnt\n",
    "export_path = psd_path\n",
    "\n",
    "fooof_dict = {}\n",
    "spec_dict = {}\n",
    "freq_dict = {}\n",
    "\n",
    "spec_list = []\n",
    "epoch_num_list = []\n",
    "for i in subject_list:\n",
    "    if (f'{i}_{project_name}_psd.json') in os.listdir(export_path):\n",
    "        print(f'{i}_{project_name}_psd.json already exists. Skipping...')\n",
    "        continue\n",
    "\n",
    "    #Get items in cnt_list that are in subject_list\n",
    "    one_two_epoch_list = []\n",
    "    for j in cnt_list:\n",
    "        if i in j:\n",
    "            file_path = os.path.join(import_path, j)\n",
    "            raw  = mne.io.read_raw(file_path, preload=True)\n",
    "\n",
    "            # Select channels\n",
    "            if fooof_channels != 'all':\n",
    "                raw.pick_channels(fooof_channels)\n",
    "            # Remove bad channels\n",
    "            if bad_fooof_channels != []:\n",
    "                raw.drop_channels(bad_fooof_channels)\n",
    "\n",
    "            # Condition to check for bug\n",
    "            if file_path.endswith('cnt'):\n",
    "                # Check if there are 10 zero values in a row in the data\n",
    "                # Check if there are any occurrences of 10 consecutive zero values\n",
    "                data = raw.to_data_frame(picks=[raw.info['ch_names'][0]])\n",
    "                has_consecutive_zeros = any(\n",
    "                    data['value'].rolling(window=10).sum() == 0)\n",
    "                # Print error message if there are 10 consecutive zeros\n",
    "                if has_consecutive_zeros:\n",
    "                    print(\n",
    "                        f'ERROR: {file_path} has 10 consecutive zero values. Please check data.')\n",
    "                    continue\n",
    "                \n",
    "            srate = raw.info['sfreq']\n",
    "            var_n_fft = int(srate * epoch_dur)\n",
    "            onsets = raw.annotations.onset\n",
    "            durations = raw.annotations.duration\n",
    "            epoch_rest = mne.make_fixed_length_epochs(raw, epoch_dur, preload=True)\n",
    "            epoch_rest.apply_baseline(baseline=baseline)\n",
    "            one_two_epoch_list.append(epoch_rest)\n",
    "    epochs = mne.concatenate_epochs(one_two_epoch_list)\n",
    "\n",
    "    psd = epochs.compute_psd(method=psd_method, n_fft=var_n_fft).average()\n",
    "    # Create json object to store psd data\n",
    "    psd_dict = {}\n",
    "    psd_dict['spec'] = psd.get_data().tolist()\n",
    "    psd_dict['freq'] = psd.freqs.tolist()\n",
    "    psd_dict['epoch_num'] = len(epochs)\n",
    "    psd_dict['subject'] = i\n",
    "    psd_dict['ch_names'] = psd.ch_names\n",
    "    # Export psd data as json\n",
    "    with open(os.path.join(psd_path, f'{i}_{project_name}_psd.json'), 'w') as f:\n",
    "        json.dump(psd_dict, f)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FOOF Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "#Get list of all files in directory\n",
    "file_list = os.listdir(psd_path)\n",
    "print(file_list)\n",
    "#Get list of all files that end with .cnt\n",
    "file_list.sort()\n",
    "\n",
    "df_fooof_list = []\n",
    "df_fooof_func_list = []\n",
    "spec_list = []\n",
    "bands = {'delta': [1, 4],\n",
    "        'theta': [4, 8],\n",
    "        'alpha': [8, 13],\n",
    "        'beta': [13, 30],\n",
    "        'gamma': [30, 50]}\n",
    "\n",
    "loop = 0\n",
    "for i in file_list:\n",
    "    if i.endswith('.json'):\n",
    "        with open(os.path.join(psd_path, i)) as f:\n",
    "            data = json.load(f)\n",
    "\n",
    "        freqs = np.array(data['freq'])\n",
    "        fooof_sub_list = []\n",
    "        fooof_func_sub_list = []\n",
    "    \n",
    "        file = i\n",
    "        sub_num = file.split('_')[0]\n",
    "\n",
    "\n",
    "        if f\"{sub_num}_{project_name}_func.csv\" in str(os.listdir(fooof_path)) and f\"{sub_num}_{project_name}.csv\" in str(os.listdir(fooof_path)):\n",
    "            print(f'{sub_num}_{project_name}_func.csv and {sub_num}_{project_name}.csv already exists. Skipping...')\n",
    "            df_fooof_func_list.append(pd.read_csv(os.path.join(fooof_path, f\"{sub_num}_{project_name}_func.csv\")))\n",
    "            df_fooof_list.append(pd.read_csv(os.path.join(fooof_path, f\"{sub_num}_{project_name}.csv\")))\n",
    "            continue\n",
    "        for index, spectrum in enumerate(data['spec']):\n",
    "            \n",
    "            \n",
    "            spectrum = np.array(spectrum)\n",
    "            # Create FOOOF object\n",
    "            fm = FOOOF(aperiodic_mode=mode, verbose=False)\n",
    "            \n",
    "            # If there is a Dataerror, skip the channel\n",
    "            try:\n",
    "                fm.fit(freqs, spectrum, freq_range=freq_range)\n",
    "                print(f'Fit: {file} at channel {data[\"ch_names\"][index]} \\n')\n",
    "            except:\n",
    "                #write to log file\n",
    "                print(f'Error in {file} at channel {data[\"ch_names\"][index]} \\n')\n",
    "                continue\n",
    "            \n",
    "            # Remove 0 Hz from the spectrum\n",
    "            freqs_no_0 = freqs[1:]\n",
    "            spectrum_no_0 = spectrum[1:]\n",
    "\n",
    "            # Generate aperiodic component\n",
    "            ap_spec = 10**gen_aperiodic(freqs_no_0, fm.aperiodic_params_) \n",
    "            p_spec = spectrum_no_0 - ap_spec\n",
    "            \n",
    "            # Create empty dataframe\n",
    "            df_fooof = pd.DataFrame()\n",
    "            df_fooof_func = pd.DataFrame()\n",
    "            # Add aperiodic series\n",
    "            df_fooof_func['pc'] = p_spec\n",
    "            df_fooof_func['apc'] = ap_spec\n",
    "            df_fooof_func['psd'] = spectrum_no_0\n",
    "            # print(f'psd: {spectrum}')\n",
    "            df_fooof_func['freqs'] = freqs_no_0\n",
    "            df_fooof_func['subject'] = data['subject']\n",
    "            df_fooof_func['channel'] = data['ch_names'][index]\n",
    "            df_fooof_func['epoch_num'] = data['epoch_num']\n",
    "            #Generate series 0-len(freqs) for plotting\n",
    "            df_fooof_func['num'] = np.arange(0, len(freqs_no_0))\n",
    "\n",
    "            offset =  fm.aperiodic_params_[0]\n",
    "            if mode == 'knee':\n",
    "                knee = fm.aperiodic_params_[1]\n",
    "                exponent = fm.aperiodic_params_[2]\n",
    "            else:\n",
    "                exponent = fm.aperiodic_params_[1]\n",
    "            df_fooof['ap_offset'] = [offset]\n",
    "            if (mode == 'knee'):\n",
    "                df_fooof['ap_knee'] = [knee]\n",
    "            df_fooof['ap_exponent'] = [exponent]\n",
    "            df_fooof['r_squared'] = [fm.r_squared_]\n",
    "            df_fooof['error'] = [fm.error_]\n",
    "            df_fooof['subject'] = data['subject']\n",
    "            df_fooof['channel'] = data['ch_names'][index]\n",
    "\n",
    "            #Add peak parameters in band\n",
    "            # Convert the peak_params_ to a dataframe\n",
    "            params_df = pd.DataFrame(fm.peak_params_, columns=['CF', 'PW', 'BW'])\n",
    "\n",
    "            # Create a function to get the peak_nums\n",
    "            def get_peak_nums(df, bands):\n",
    "                peak_nums = {}\n",
    "                for band in bands:\n",
    "                    peak_nums[band] = len(df[df['CF'].between(bands[band][0], bands[band][1], inclusive='left')])\n",
    "                return peak_nums\n",
    "\n",
    "            # Call the function to get the peak_nums dictionary\n",
    "            peak_nums = get_peak_nums(params_df, bands)\n",
    "\n",
    "            # Iterate over the bands and create the columns in df_fooof\n",
    "            for band, band_limits in bands.items():\n",
    "                band_df = params_df[params_df['CF'].between(bands[band][0], bands[band][1], inclusive='left')].reset_index(drop=True)\n",
    "                for i in range(peak_nums[band]):\n",
    "                    df_fooof[f'{band}_CF_{i}'] = [band_df.loc[i, 'CF']]\n",
    "                    df_fooof[f'{band}_PW_{i}'] = [band_df.loc[i, 'PW']]\n",
    "                    df_fooof[f'{band}_BW_{i}'] = [band_df.loc[i, 'BW']]\n",
    "            \n",
    "            df_fooof_list.append(df_fooof)\n",
    "            df_fooof_func_list.append(df_fooof_func)\n",
    "         \n",
    "            #Add data to list\n",
    "            fooof_sub_list.append(df_fooof)\n",
    "            fooof_func_sub_list.append(df_fooof_func)\n",
    "\n",
    "    fooof_sub = pd.concat(fooof_sub_list)\n",
    "    fooof_func_sub = pd.concat(fooof_func_sub_list)\n",
    "    \n",
    "    \n",
    "    # Reset index\n",
    "    fooof_sub = fooof_sub.reset_index(drop=True)\n",
    "    fooof_func_sub = fooof_func_sub.reset_index(drop=True)\n",
    "    fooof_sub.to_csv(os.path.join(fooof_path,f\"{sub_num}_{project_name}.csv\"))\n",
    "    fooof_func_sub.to_csv(os.path.join(fooof_path, f\"{sub_num}_{project_name}_func.csv\"))\n",
    "\n",
    "#Concatenate all dataframes\n",
    "df_fooof = pd.concat(df_fooof_list)\n",
    "df_fooof_func = pd.concat(df_fooof_func_list)\n",
    "\n",
    "# Reset index\n",
    "df_fooof = df_fooof.reset_index(drop=True)\n",
    "df_fooof_func = df_fooof_func.reset_index(drop=True)\n",
    "\n",
    "\n",
    "now = datetime.datetime.now()\n",
    "now = now.strftime(\"%Y-%m-%d_%H-%M\")\n",
    "\n",
    "# Export to csv\n",
    "df_fooof.to_csv(f\"fooof_comb_{project_name}_{now}.csv\")\n",
    "df_fooof_func.to_csv(f\"fooof_func_comb_{project_name}_{now}.csv\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convert to wide"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Get list of all files in directory\n",
    "file_list = os.listdir()\n",
    "print(file_list)\n",
    "#Get list of all files that end with .cnt\n",
    "file_list.sort()\n",
    "print(file_list)\n",
    "\n",
    "#Get all files that contain string 'fooof_'\n",
    "fooof_file_list = [i for i in file_list if ('fooof_comb' in i and i.endswith('.csv'))]\n",
    "fooof_func_file_list = [i for i in file_list if ('fooof_func_comb' in i and i.endswith('.csv'))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import last df_fooof and df_fooof_func file\n",
    "df_fooof = pd.read_csv(fooof_file_list[-1])\n",
    "df_fooof_func = pd.read_csv(fooof_func_file_list[-1])\n",
    "\n",
    "#Drop all columns that contain 'Unnamed'\n",
    "df_fooof = df_fooof.loc[:, ~df_fooof.columns.str.contains('^Unnamed')]\n",
    "df_fooof_func = df_fooof_func.loc[:, ~df_fooof_func.columns.str.contains('^Unnamed')]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_fooof"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_fooof_func"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Round frequencies to 2 decimals\n",
    "df_fooof_func['freq'] = df_fooof_func['freqs'].round(2) \n",
    "df_fooof_func['freq'] = df_fooof_func['freq'].astype(str) + 'H'\n",
    "df_fooof_func\n",
    "df_select = df_fooof_func[df_fooof_func['num'] <120]\n",
    "\n",
    "# Convert offset to microvolts\n",
    "df_fooof_microvolts = df_fooof.apply(lambda x: np.log10(10**x*10**12) if 'offset' in x.name else x, axis=0)\n",
    "\n",
    "#Get column names in df_fooof\n",
    "fooof_cols = df_fooof.columns.to_list()\n",
    "#Remove subject and channel\n",
    "fooof_cols.remove('subject')\n",
    "fooof_cols.remove('channel')\n",
    "fooof_cols"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "wide_df_fooof_microvolts = df_fooof_microvolts.pivot(\n",
    "                    index=['subject'], \n",
    "                    columns=['channel'], \n",
    "                    values=fooof_cols\n",
    "                    )\n",
    "\n",
    "\n",
    "value = wide_df_fooof_microvolts.columns.get_level_values(0)\n",
    "channel = wide_df_fooof_microvolts.columns.get_level_values(1)\n",
    "\n",
    "# #Set the index to source_fileS\n",
    "# wide_df_fcz = wide_df_fcz.set_index(('lab_num'))\n",
    "\n",
    "#Rename the columns\n",
    "wide_df_fooof_microvolts.columns = project_name + '_' + value + '_' + channel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "wide_df_fooof_func = df_select.pivot(\n",
    "                    index=['subject'], \n",
    "                    columns=['channel', 'freq'], \n",
    "                    values=['psd', 'apc', 'pc', 'epoch_num']\n",
    "                    )\n",
    "\n",
    "\n",
    "value = wide_df_fooof_func.columns.get_level_values(0)\n",
    "channel = wide_df_fooof_func.columns.get_level_values(1)\n",
    "num = wide_df_fooof_func.columns.get_level_values(2)\n",
    "\n",
    "# #Set the index to source_fileS\n",
    "# wide_df_fcz = wide_df_fcz.set_index(('lab_num'))\n",
    "\n",
    "#Rename the columns\n",
    "wide_df_fooof_func.columns = project_name + '_' + value + '_' + channel + '_' + num "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Multiply every column but the last one by 10^12\n",
    "wide_fooof_func_microvolts = pd.concat([wide_df_fooof_func.filter(regex='(psd)|(apc)|(pc)') * 1e12, wide_df_fooof_func.filter(regex='(epoch_num)')], axis=1)\n",
    "\n",
    "#Reset index\n",
    "wide_fooof_func_microvolts.reset_index(inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epoch_cols = [col for col in wide_fooof_func_microvolts.columns if 'epoch_num' in col]\n",
    "wide_fooof_func_microvolts = wide_fooof_func_microvolts.drop(epoch_cols[1:], axis=1)\n",
    "\n",
    "# Rename column that has epoch_num in it to epoch_num\n",
    "epoch_num_name = [col for col in wide_fooof_func_microvolts.columns if 'epoch_num' in col]\n",
    "wide_fooof_func_microvolts = wide_fooof_func_microvolts.rename(columns={epoch_num_name[0]: 'epoch_num'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wide_fooof_func_microvolts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Convert index to string\n",
    "wide_df_fooof_microvolts.index = wide_df_fooof_microvolts.index.astype(str)\n",
    "#Convert index to string\n",
    "wide_fooof_func_microvolts.index = wide_fooof_func_microvolts.index.astype(str)\n",
    "\n",
    "#Reset index\n",
    "wide_df_fooof_microvolts.reset_index(inplace=True)\n",
    "#Reset index\n",
    "wide_fooof_func_microvolts.reset_index(inplace=True)\n",
    "\n",
    "#Set index to subject\n",
    "wide_df_fooof_microvolts = wide_df_fooof_microvolts.set_index('subject')\n",
    "#Set index to subject\n",
    "wide_fooof_func_microvolts = wide_fooof_func_microvolts.set_index('subject')\n",
    "\n",
    "# Set index to string\n",
    "wide_df_fooof_microvolts.index = wide_df_fooof_microvolts.index.astype(str)\n",
    "# Set index to string\n",
    "wide_fooof_func_microvolts.index = wide_fooof_func_microvolts.index.astype(str)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Concatenate the two dataframes\n",
    "wide_df = pd.concat([ wide_fooof_func_microvolts, wide_df_fooof_microvolts], axis=1)\n",
    "#Replace offset with off\n",
    "wide_df.columns = wide_df.columns.str.replace('offset', 'off')\n",
    "#Replace exponent with exp\n",
    "wide_df.columns = wide_df.columns.str.replace('exponent', 'exp')\n",
    "\n",
    "wide_df.reset_index(inplace=True)\n",
    "wide_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wide_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "now = datetime.datetime.now()\n",
    "now = now.strftime(\"%Y-%m-%d_%H-%M\")\n",
    "pyreadstat.write_sav(wide_df, os.path.join(home_dir, f\"final_fooof_{project_name}_{now}.sav\"))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis\n",
    "- Freq and FOOOF plots\n",
    "- Raincloud plots\n",
    "- Stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "now = datetime.datetime.now()\n",
    "now = now.strftime(\"%Y-%m-%d_%H-%M\")\n",
    "#Get list of all files in directory\n",
    "file_list = os.listdir(home_dir)\n",
    "print(file_list)\n",
    "#Get list of all files that end with .cnt\n",
    "file_list.sort()\n",
    "print(file_list)\n",
    "\n",
    "#Get all files that contain string 'fooof_'\n",
    "fooof_file_list = [i for i in file_list if ('fooof_comb' in i and i.endswith('.csv'))]\n",
    "fooof_func_file_list = [i for i in file_list if ('fooof_func_comb' in i and i.endswith('.csv'))]\n",
    "spss_file_list = [i for i in file_list if ('final_fooof' in i and i.endswith('.sav'))]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import last df_fooof and df_fooof_func file\n",
    "df_fooof = pd.read_csv(fooof_file_list[-1])\n",
    "df_fooof_func = pd.read_csv(fooof_func_file_list[-1])\n",
    "final_fooof = pd.read_spss(spss_file_list[-1])\n",
    "\n",
    "# Drop columns subject and index\n",
    "final_fooof = final_fooof.drop(columns=['subject', 'index'])\n",
    "\n",
    "# Convert offset to microvolts\n",
    "df_fooof = df_fooof.apply(lambda x: np.log10(10**x*10**12) if 'offset' in x.name else x, axis=0)\n",
    "\n",
    "# Convert psd, apc, and pc to microvolts\n",
    "df_fooof_func = df_fooof_func.apply(lambda x: x*10**12 if 'psd' in x.name or 'apc' in x.name or 'pc' in x.name else x, axis=0)\n",
    "\n",
    "# set dtype to numeric in final_fooof\n",
    "final_fooof = final_fooof.apply(pd.to_numeric, errors='ignore')\n",
    "\n",
    "\n",
    "channels_regex  = [f'(_{i}_)' for i in analysis_channels]\n",
    "channels_regex = '|'.join(channels_regex)\n",
    "\n",
    "# Make column names lower case\n",
    "final_fooof.columns = final_fooof.columns.str.lower()\n",
    "\n",
    "# Set all text to lower case in df_fooof_func and df_fooof\n",
    "df_fooof_func = df_fooof_func.apply(lambda x: x.astype(str).str.lower() if 'channel' in x.name else x, axis=0)\n",
    "df_fooof = df_fooof.apply(lambda x: x.astype(str).str.lower() if 'channel' in x.name else x, axis=0)\n",
    "\n",
    "df_fooof_func = df_fooof_func.loc[df_fooof_func['channel'].isin(analysis_channels)]\n",
    "df_fooof = df_fooof.loc[df_fooof['channel'].isin(analysis_channels)]\n",
    "\n",
    "final_fooof = final_fooof.filter(regex=channels_regex)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_fooof_func_grouped = df_fooof_func.groupby(['subject', 'channel'])\n",
    "# Loop through each group \n",
    "for name, group in df_fooof_func_grouped:\n",
    "\n",
    "    offset = df_fooof.loc[(df_fooof['subject'] == name[0]) & (df_fooof['channel'] == name[1]), 'ap_offset'].values[0].round(2)\n",
    "    if mode == 'knee':\n",
    "        knee = df_fooof.loc[(df_fooof['subject'] == name[0]) & (df_fooof['channel'] == name[1]), 'ap_knee'].values[0].round(2)\n",
    "    elif mode == 'fixed':\n",
    "        knee = 'NA'\n",
    "    exponent = df_fooof.loc[(df_fooof['subject'] == name[0]) & (df_fooof['channel'] == name[1]), 'ap_exponent'].values[0].round(2)\n",
    "\n",
    "    # Plot psd, apc, and pc\n",
    "    fig, ax = plt.subplots()\n",
    "    ax = group.plot(x='freqs', y= ['psd', 'apc', 'pc'], ax=ax, title=f'Subject: {name[0]} Channel: {name[1]}')\n",
    "    ax.set_xlabel('Frequency (Hz)')\n",
    "    ax.set_ylabel('Power (uV^2/Hz)')\n",
    "    ax.legend(['PSD', 'APC', 'PC'])\n",
    "    ax.set_xlim([0, 55])\n",
    "    # Place the Following Text at the Bottom of the Figure\n",
    "    fig.text(0.5, -0.01, f'Offset: {offset} Knee: {knee} Exponent: {exponent}', ha='center', fontsize=12)\n",
    "    \n",
    "    # Save the Figure\n",
    "    fig.tight_layout()\n",
    "    fig.savefig(os.path.join(analysis_path, f'fooof_{name[0]}_{name[1]}_plot_{now}.png'), dpi=300, bbox_inches='tight')\n",
    "\n",
    "    plt.close()\n",
    "\n",
    "    # Plot psd, apc, and pc\n",
    "    fig, ax = plt.subplots()\n",
    "    ax = group.plot(x='freqs', y= ['psd', 'apc', 'pc'], ax=ax, title=f'Subject: {name[0]} Channel: {name[1]}', logy=True)\n",
    "    ax.set_xlabel('Frequency (Hz)')\n",
    "    ax.set_ylabel('Power (uV^2/Hz)')\n",
    "    ax.legend(['PSD', 'APC', 'PC'])\n",
    "    ax.set_xlim([0, 55])\n",
    "    # Place the Following Text at the Bottom of the Figure\n",
    "    fig.text(0.5, -0.01, f'Offset: {offset} Knee: {knee} Exponent: {exponent}', ha='center', fontsize=12)\n",
    "    \n",
    "    # Save the Figure\n",
    "    fig.tight_layout()\n",
    "    fig.savefig(os.path.join(analysis_path, f'fooof_{name[0]}_{name[1]}_logy_{now}.png'), dpi=300, bbox_inches='tight')\n",
    "\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "# Plot every column in final_fooof as a raincloud plot seperate raincloud plot\n",
    "for column in final_fooof.columns:\n",
    "    fig = plt.figure()\n",
    "\n",
    "    # Remove outliers from final_fooof[column]\n",
    "    z_score = 3\n",
    "    data_no_out = final_fooof[column][~((final_fooof[column] < (final_fooof[column].mean() - z_score * final_fooof[column].std())) |\n",
    "                                        (final_fooof[column] > (final_fooof[column].mean() + z_score * final_fooof[column].std())))]\n",
    "\n",
    "    ax = pt.RainCloud(data=final_fooof[column], orient='h')\n",
    "    ax.set_title(column)\n",
    "    \n",
    "    mean = round(float(final_fooof[column].mean()), 2)\n",
    "    median = round(float(final_fooof[column].median()), 2)\n",
    "    sd = round(float(final_fooof[column].std()), 2)\n",
    "\n",
    "    # Place descriptive statistics at the bottom of the figure\n",
    "    fig.text(.5, 0.02, f'Mean: {mean} Median: {median} SD: {sd}', ha='center', fontsize=12)\n",
    "\n",
    "    plt.savefig(os.path.join(analysis_path, f'fooof_{column}_raincloud_{now}.png'), dpi=300)\n",
    "    # wait for user to press enter to continue\n",
    "    \n",
    "    plt.close()\n",
    "\n",
    "    fig = plt.figure()\n",
    "\n",
    "    ax = pt.RainCloud(data=data_no_out, orient='h')\n",
    "    ax.set_title(f'{column}_noout')\n",
    "    \n",
    "    mean = round(float(data_no_out.mean()), 2)\n",
    "    median = round(float(data_no_out.median()), 2)\n",
    "    sd = round(float(data_no_out.std()), 2)\n",
    "\n",
    "    # Place descriptive statistics at the bottom of the figure\n",
    "    fig.text(.5, 0.02, f'Mean: {mean} Median: {median} SD: {sd}', ha='center', fontsize=12)\n",
    "\n",
    "    plt.savefig(os.path.join(analysis_path, f'fooof_{column}_raincloud_noout_{now}.png'), dpi=300)\n",
    "    # wait for user to press enter to continue\n",
    "    \n",
    "    plt.close()\n",
    "\n",
    "\n",
    "    "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "mne",
   "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.10.10"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}