step1_compute-coavitcation.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute co-activation\n",
    "\n",
    "## NOTE: Data is not currently shared. This code is only provided for reference.\n",
    "\n",
    "Ccomputing co-activation requires a Guassian resampling to provide fewer datapoints per seocnd (a method for downsampling), so the analysis can be perfomred more handily.\n",
    "Futher, we need to define the *target-pairs* and then quantify coactivations for those muscles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2022-11-11 19:30:34,386\tINFO worker.py:1509 -- Started a local Ray instance. View the dashboard at \u001b[1m\u001b[32m127.0.0.1:8265 \u001b[39m\u001b[22m\n"
     ]
    }
   ],
   "source": [
    "# Import libraries and setup values\n",
    "import os\n",
    "import pickle as pkl\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import ray\n",
    "from scipy.io import loadmat\n",
    "\n",
    "import EssentialEMGFuncs as ess\n",
    "\n",
    "ray.init(dashboard_host='127.0.0.1', ignore_reinit_error=True)\n",
    "save_data = False\n",
    "load_savedData = True\n",
    "resample_data, use_resampledData = False, True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup specific variables\n",
    "SUBJS = {'young': [\"PS04\", \"PS05\", \"PS06\", \"PS07\", \"PS15\", \"PS16\", \"PS17\", \"PS18\", \"PS19\", \"PS20\", \"PS21\", \"PS22\",\n",
    "                   \"PS23\", \"PS24\", \"PS25\", \"PS26\", \"PS27\"],\n",
    "         'old': [\"PS08\", \"PS09\", \"PS10\", \"PS28\", \"PS29\", \"PS31\", \"PS32\", \"PS33\", \"PS34\", \"PS35\", \"PS36\"]\n",
    "         }\n",
    "CONDS = ['LME', 'LEI', 'RME', 'REI']\n",
    "BLOCKS = {'pre': [0], 'pert': [1, 3], 'catch': [2], 'post': [4]}\n",
    "DATA_PATH = './data/emg/'\n",
    "SUBJ_PATH = {s: DATA_PATH + s + '/' for a in SUBJS for s in SUBJS[a]}\n",
    "RESMAPLED_DATA = 'gaussian-resmped_young-old_raw-data.pickle'\n",
    "COACTIVATION_DATA = 'coactivation_young-and-old_resampled_perturb-recovery-steps.pickle'\n",
    "STEP_DATA = 'stepping_time_young-and-old.pickle'\n",
    "\n",
    "if not os.path.isdir(DATA_PATH + 'group_response'):\n",
    "    os.mkdir(DATA_PATH + 'group_response')\n",
    "for a in SUBJS:\n",
    "    for s in SUBJS[a]:\n",
    "        if not os.path.isdir(SUBJ_PATH[s] + 'figs'):\n",
    "            os.mkdir(SUBJ_PATH[s] + 'figs')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load individual EMG data\n",
    "if not load_savedData:\n",
    "    par_read_feather = ray.remote(ess.read_normal_profile)\n",
    "    normal_envelope, read_idx = [{c: {s: [] for a in SUBJS for s in SUBJS[a]} for c in CONDS} for _ in range(2)]\n",
    "    for a in SUBJS:\n",
    "        for s in SUBJS[a]:\n",
    "            for c in CONDS:\n",
    "                read_idx[c][s] = par_read_feather.remote(f'{SUBJ_PATH[s]}{s}_{c}_normalEnvelope.feather')\n",
    "    for a in SUBJS:\n",
    "        for s in SUBJS[a]:\n",
    "            for c in CONDS:\n",
    "                try:\n",
    "                    normal_envelope[c][s] = ray.get(read_idx[c][s]).copy()\n",
    "                    del read_idx[c][s]\n",
    "                except Exception:\n",
    "                    print(f'cannot load {s}_{c}, disregard if you are aware of this unavalibity')\n",
    "                    normal_envelope[c][s] = pd.DataFrame()\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gaussian resampling\n",
    "It turns out that 2000 time points is a lot for this co-contraction analysis.\n",
    "So, let's have it with 200 time points. We will use Gaussian smoothing from SMART for this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if resample_data and not load_savedData:\n",
    "    kernel_size, num_dp = 0.01, 200\n",
    "    resamp_envelope = {c: {s: pd.DataFrame(data=np.zeros((num_dp, len(normal_envelope[c][s].columns))),\n",
    "                           columns=normal_envelope[c][s].columns) for a in SUBJS for s in SUBJS[a]} for c in CONDS}\n",
    "    for a in SUBJS:\n",
    "        for s in SUBJS[a]:\n",
    "            for c in CONDS:\n",
    "                if not normal_envelope[c][s].empty:\n",
    "                    for m in resamp_envelope[c][s].columns.get_level_values(0).unique().to_list():\n",
    "                        resamp_envelope[c][s][m] = ess.gaussian_smoothing_df(\n",
    "                            normal_envelope[c][s][m], normal_envelope[c][s]['time'], kernel_size, num_dp)\n",
    "                    print(f'completed resampling for {s}_{c}')\n",
    "\n",
    "    if save_data:\n",
    "        pkl.dump(resamp_envelope, open(\n",
    "            DATA_PATH + 'group_response/' + RESMAPLED_DATA, 'wb+'))\n",
    "elif load_savedData:\n",
    "    resamp_envelope = pkl.load(open(DATA_PATH + 'group_response/' + RESMAPLED_DATA, 'rb'))\n",
    "\n",
    "if use_resampledData:\n",
    "    normal_envelope = resamp_envelope\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load step times\n",
    "step_time = {c: {s: [] for a in SUBJS for s in SUBJS[a]} for c in CONDS}\n",
    "for a in SUBJS:\n",
    "    for s in SUBJS[a]:\n",
    "        for c in CONDS:\n",
    "            data = loadmat(SUBJ_PATH[s] + 'times/' + s + '_' + c + '_time.mat')\n",
    "            step_time[c][s] = pd.DataFrame(\n",
    "                data=data['strideTime'], columns=['step_index', 'type', 'start', 'pert', 'other', 'end'])\n",
    "            step_time[c][s] = ess.fill_perturb_time(step_time[c][s])\n",
    "\n",
    "if save_data:\n",
    "    pkl.dump(step_time, open(\n",
    "        DATA_PATH + 'group_response/' + STEP_DATA, 'wb+'))\n",
    "\n",
    "# Separate the steps into pre, post, perturbations and catches\n",
    "idx = pd.IndexSlice\n",
    "seperated_envelope = {c: {s: {b: [] for b in BLOCKS} for a in SUBJS for s in SUBJS[a]} for c in CONDS}\n",
    "for a in SUBJS:\n",
    "    for s in SUBJS[a]:\n",
    "        for c in CONDS:\n",
    "            last_pert = step_time[c][s][[c in BLOCKS['pert'] for c in step_time[c][s]['type']]].index[-1]\n",
    "            step_time[c][s].loc[last_pert + 1:, 'type'] = 4\n",
    "            if not bool(normal_envelope[c][s].empty):\n",
    "                col0_length = len(normal_envelope[c][s].columns.get_level_values(0).unique())\n",
    "                for b in BLOCKS:\n",
    "                    seperated_envelope[c][s][b] = normal_envelope[c][s].loc[idx[:], idx[:, np.tile(\n",
    "                        [t in BLOCKS[b] for t in step_time[c][s]['type']], (1, col0_length)).squeeze()]]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rename the muscles we need\n",
    "MUSCLES = {'LTA': ['l tibia'], 'LSO': ['l sol'], 'LRF': ['l Vast', 'l rect'], 'LST': ['l semi'],\n",
    "           'LAD': ['l deltoid ant', 'l ant deltoid'], 'LPD': ['l deltoid post', 'l post deltoid'], 'RTA': ['r tibia'],\n",
    "           'RSO': ['r sol'], 'RRF': ['r Vast', 'r rect'], 'RST': ['r semi'], 'RAD': ['r deltoid ant', 'r ant deltoid'],\n",
    "           'RPD': ['r deltoid post', 'r post deltoid']}  # the full name for the muscles paired with their keys\n",
    "\n",
    "for a in SUBJS:\n",
    "    for s in SUBJS[a]:\n",
    "        for c in CONDS:\n",
    "            for b in BLOCKS:\n",
    "                try:\n",
    "                    seperated_envelope[c][s][b].sort_index(axis=1, level=0, inplace=True)\n",
    "                    cols = seperated_envelope[c][s][b].columns.get_level_values(0).unique().to_list()\n",
    "                    for i, m in enumerate(cols):\n",
    "                        if any(t.upper() in m for mk in MUSCLES for t in MUSCLES[mk]):\n",
    "                            cols[i] = [mk for mk in MUSCLES for t in MUSCLES[mk] if t.upper() in m][0]\n",
    "                    seperated_envelope[c][s][b].columns.set_levels(cols, level=0, inplace=True)\n",
    "                    seperated_envelope[c][s][b].sort_index(axis=1, level=0, inplace=True)\n",
    "                except Exception:\n",
    "                    print(f'there was a problem in {a}:{s}_{c}_{b}')\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute co-activation\n",
    "Simialr to the Guassian resmapling, the coactivation analysis needs using the `Ray` toolbox, so we can benefit form all the computational capacity fo the processor.\n",
    "Co-activation requres *target-pairs* to compute if they are working agasints each other. This is the general case of co-contraction, which is the co-activatoin of the agonist and antagonist muscles acting on the same joint. As our study is related to the co-aontraction analysis, we defeined the *target-pairs* as the muscle pairs that act on the same joint. For each step, the first muscle is agonist and the second muscle in the pair is antagonist (See table 1 in the paper)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "target_pair = {\n",
    "    'left_step': {\n",
    "        'LTA-LSO': ['LSO', 'LTA'],\n",
    "        'LRF-LST': ['LRF', 'LST'],\n",
    "        'LAD-LPD': ['LPD', 'LAD'],\n",
    "        'RTA-RSO': ['RTA', 'RSO'],\n",
    "        'RRF-RST': ['RST', 'RRF'],\n",
    "        'RAD-RPD': ['RAD', 'RPD']\n",
    "    },\n",
    "    'right_step': {\n",
    "        'LTA-LSO': ['LTA', 'LSO'],\n",
    "        'LRF-LST': ['LST', 'LRF'],\n",
    "        'LAD-LPD': ['LAD', 'LPD'],\n",
    "        'RTA-RSO': ['RSO', 'RTA'],\n",
    "        'RRF-RST': ['RRF', 'RST'],\n",
    "        'RAD-RPD': ['RPD', 'RAD']\n",
    "    }\n",
    "}\n",
    "tp_names = list(target_pair['left_step'])\n",
    "TESTS = {'perturbed': ['start', 'other'], 'recovery': ['other', 'end']}\n",
    "par_compute_coActivation = ray.remote(ess.compute_coActivation)\n",
    "coact_idx = {t: {c: {s: {b: {p: [] for p in tp_names} for b in BLOCKS} for a in SUBJS for s in SUBJS[a]}\n",
    "                 for c in CONDS} for t in TESTS}\n",
    "\n",
    "for t in TESTS:\n",
    "    for a in SUBJS:\n",
    "        for s in SUBJS[a]:\n",
    "            for c in CONDS:\n",
    "                for b in BLOCKS:\n",
    "                    try:\n",
    "                        if c in ['LEI', 'LME']:\n",
    "                            if t == 'perturbed':\n",
    "                                t_p = target_pair['left_step']\n",
    "                            else:\n",
    "                                t_p = target_pair['right_step']\n",
    "                        else:\n",
    "                            if t == 'perturbed':\n",
    "                                t_p = target_pair['right_step']\n",
    "                            else:\n",
    "                                t_p = target_pair['left_step']\n",
    "                        strides = seperated_envelope[c][s][b].columns.get_level_values(1).unique().to_list()\n",
    "                        s_t = step_time[c][s].iloc[strides, :]\n",
    "                        for p in t_p:\n",
    "                            coact_idx[t][c][s][b][p] = par_compute_coActivation.remote(\n",
    "                                seperated_envelope[c][s][b][t_p[p]], seperated_envelope[c][s][b]['time'], s_t,\n",
    "                                fillnopert=False, event=TESTS[t][0], duration=TESTS[t][1], mode='fixed')\n",
    "                    except Exception:\n",
    "                        print(f'there was a problem in {a}:{s}_{c}_{b}')\n",
    "# Sanity check\n",
    "# for t in TESTS:\n",
    "#     for a in SUBJS:\n",
    "#         for s in SUBJS[a]:\n",
    "#             for c in CONDS:\n",
    "#                 for b in BLOCKS:\n",
    "#                     if c in ['LEI', 'LME']:\n",
    "#                         if t != 'perturbed':\n",
    "#                             t_p = target_pair['left_step']\n",
    "#                         else:\n",
    "#                             t_p = target_pair['right_step']\n",
    "#                     else:\n",
    "#                         if t != 'perturbed':\n",
    "#                             t_p = target_pair['right_step']\n",
    "#                         else:\n",
    "#                             t_p = target_pair['left_step']\n",
    "#                     strides = seperated_envelope[c][s][b].columns.get_level_values(1).unique().to_list()\n",
    "#                     s_t = step_time[c][s].iloc[strides, :]\n",
    "#                     for p in t_p:\n",
    "#                         coactivation[t][c][s][b][p] = ess.compute_coActivation(\n",
    "#                             seperated_envelope[c][s][b][t_p[p]], seperated_envelope[c][s][b]['time'], s_t,\n",
    "#                             fillnopert=False, event=TESTS[t][0], duration=TESTS[t][1], mode='wasted_contraction')\n",
    "\n",
    "coactivation = {t: {c: {s: {b: pd.DataFrame(columns=tp_names) for b in BLOCKS} for a in SUBJS for s in SUBJS[a]}\n",
    "                    for c in CONDS} for t in TESTS}\n",
    "for t in TESTS:\n",
    "    for a in SUBJS:\n",
    "        for s in SUBJS[a]:\n",
    "            for c in CONDS:\n",
    "                for b in BLOCKS:\n",
    "                    try:\n",
    "                        for p in tp_names:\n",
    "                            coactivation[t][c][s][b][p] = ray.get(coact_idx[t][c][s][b][p]).copy()\n",
    "                    except Exception:\n",
    "                        print(f'there was a problem in {a}:{s}_{c}_{b}')\n",
    "if save_data:\n",
    "    pkl.dump(coactivation, open(DATA_PATH + 'group_response/' + COACTIVATION_DATA, 'wb+'))\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "torch-beta",
   "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.13"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "7cc401fccc55e2825ef583b7b43bb1f971ba34a72aecd6f049901fe5d37b7bbc"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}