accelerometer_ventilated_analysis_2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "![alt text](./Cerny_logo_1.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The effect of ambulance acceleration on mechanical ventilation during neonatal transport\n",
    "\n",
    "#### Author: Dr Gusztav Belteki"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Import the required libraries and set options"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import IPython\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "import os\n",
    "import sys\n",
    "import pickle\n",
    "\n",
    "from scipy import stats\n",
    "from pandas import Series, DataFrame\n",
    "from datetime import datetime, timedelta\n",
    "from matplotlib.ticker import ScalarFormatter\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "%matplotlib inline\n",
    "matplotlib.style.use('classic')\n",
    "matplotlib.rcParams['figure.facecolor'] = 'w'\n",
    "#matplotlib.rcParams['date.autoformatter.second'] = '%H:%M:%S.%f'\n",
    "\n",
    "pd.set_option('display.max_rows', 250)\n",
    "pd.set_option('display.max_columns', 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Python version: {}\".format(sys.version))\n",
    "print(\"pandas version: {}\".format(pd.__version__))\n",
    "print(\"matplotlib version: {}\".format(matplotlib.__version__))\n",
    "print(\"seaborn version: {}\".format(sns.__version__))\n",
    "print(\"NumPy version: {}\".format(np.__version__))\n",
    "print(\"SciPy version: {}\".format(sp.__version__))\n",
    "print(\"IPython version: {}\".format(IPython.__version__))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. List and set the working directory and the directory to write out data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Topic of the Notebook which will also be the name of the subfolder containing results\n",
    "TOPIC = 'accelerometer_ventilated'\n",
    "\n",
    "# Name of the external hard drive\n",
    "DRIVE = 'GUSZTI'\n",
    "\n",
    "# Directory containing clinical and blood gas data\n",
    "CWD = os.path.join('/Users', 'guszti', 'ventilation_fabian')\n",
    "\n",
    "DIR_WRITE = os.path.join(CWD, 'Analyses', TOPIC)\n",
    "DATA_DUMP = os.path.join('/Volumes', DRIVE, 'data_dump', 'fabian')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DIR_WRITE, DATA_DUMP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Import ventilator, accelerometer and clinical data from pickle archives"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time \n",
    "with open('%s/%s.pickle' % (DATA_DUMP, 'clin_df_1_1100'), 'rb') as handle:\n",
    "    clin_df = pickle.load(handle)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "with open('%s/%s.pickle' % (DATA_DUMP, 'combined_ventilated_1min'), 'rb') as handle:\n",
    "    combined_ventilated_1min = pickle.load(handle)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "with open('%s/%s.pickle' % (DATA_DUMP, 'ventilator_mode'), 'rb') as handle:\n",
    "    ventilator_mode = pickle.load(handle)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "with open('%s/%s.pickle' % (DATA_DUMP, 'accelero_ventilated_2'), 'rb') as handle:\n",
    "    accelero_ventilated = pickle.load(handle)\n",
    "\n",
    "with open('%s/%s.pickle' % (DATA_DUMP, 'accelero_ventilated_low_pass'), 'rb') as handle:\n",
    "    accelero_ventilated_low_pass = pickle.load(handle)\n",
    "    \n",
    "with open('%s/%s.pickle' % (DATA_DUMP, 'accelero_ventilated_high_pass'), 'rb') as handle:\n",
    "    accelero_ventilated_high_pass = pickle.load(handle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Remove those recordings whose postnatal age is over term + 6 weeks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "patients = [key[1] for key in combined_ventilated_1min.keys()]\n",
    "len(patients)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clin_df = clin_df.loc[patients]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(clin_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "clin_df[clin_df['Corrected gestational Age (weeks)'] > 46]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clin_df = clin_df[clin_df['Corrected gestational Age (weeks)'] <= 46]\n",
    "len(clin_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "patients = sorted(clin_df.index)\n",
    "len(patients)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Statistics on clinical details of the final set of patients\n",
    "\n",
    "The _postnatal age_ shown in the table is the age at the end of the transfer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clin_df.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clin_df_stats = round(clin_df.describe(percentiles = [0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99]), 2)\n",
    "clin_df_stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Boxplots on selected clinical parameters "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize = (3, 4))\n",
    "dpi = 300\n",
    "filetype = 'jpg'\n",
    "xticklabels = ['Gestational age']\n",
    "\n",
    "# Define styling for each boxplot component\n",
    "medianprops = {'color': 'black', 'linewidth': 2}\n",
    "meanpointprops = {'marker':'D', 'markeredgecolor':'black', 'markerfacecolor':'black'}\n",
    "boxprops = {'color': 'black', 'linestyle': '-'}\n",
    "whiskerprops = {'color': 'black', 'linestyle': '-'}\n",
    "capprops = {'color': 'black', 'linestyle': '-'}\n",
    "flierprops = {'color': 'black', 'marker': '.'}\n",
    "\n",
    "plt.boxplot(clin_df['Gestational Age (weeks)'] ,\n",
    "    widths = 0.3, whis = [0, 100], showfliers = True, showmeans = True, medianprops=medianprops, \n",
    "    meanprops = meanpointprops, boxprops=boxprops, whiskerprops=whiskerprops, capprops=capprops, \n",
    "    flierprops = flierprops)\n",
    "\n",
    "ax.set_ylabel('weeks', size = 14)\n",
    "ax.tick_params(axis='both', which='major', labelsize=14)\n",
    "ax.set_xticklabels(xticklabels)\n",
    "plt.grid(True)\n",
    "\n",
    "fig.savefig('%s/%s.%s' % (DIR_WRITE, 'accelero_gest_age', filetype),\n",
    "    dpi = dpi, format = filetype, bbox_inches='tight', pad_inches=0.1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize = (3, 4))\n",
    "dpi = 300\n",
    "filetype = 'jpg'\n",
    "xticklabels = ['Birth weight']\n",
    "\n",
    "medianprops = {'color': 'black', 'linewidth': 2}\n",
    "meanpointprops = {'marker':'D', 'markeredgecolor':'black', 'markerfacecolor':'black'}\n",
    "boxprops = {'color': 'black', 'linestyle': '-'}\n",
    "whiskerprops = {'color': 'black', 'linestyle': '-'}\n",
    "capprops = {'color': 'black', 'linestyle': '-'}\n",
    "flierprops = {'color': 'black', 'marker': '.'}\n",
    "\n",
    "plt.boxplot(clin_df['Birth Weight'], widths = 0.3, whis = [0, 100], showfliers = True,\n",
    "             showmeans = True, meanprops = meanpointprops, medianprops=medianprops, boxprops=boxprops, \n",
    "        whiskerprops=whiskerprops, capprops=capprops, flierprops = flierprops)\n",
    "\n",
    "ax.set_xticklabels(xticklabels)\n",
    "ax.set_ylabel('grams', size = 14)\n",
    "ax.tick_params(axis='both', which='major', labelsize=14)\n",
    "plt.grid(True)\n",
    "\n",
    "fig.savefig('%s/%s.%s' % (DIR_WRITE, 'accelero_birth_weight', filetype),\n",
    "    dpi = dpi, format = filetype, bbox_inches='tight', pad_inches=0.1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "below_29 = sorted(clin_df[clin_df['Corrected gestational Age (weeks)'] < 29].index)\n",
    "len(below_29)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "between_29_32 = sorted(clin_df[(clin_df['Corrected gestational Age (weeks)'] >= 29) & \n",
    "                       (clin_df['Corrected gestational Age (weeks)'] < 33)].index)\n",
    "len(between_29_32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "between_33_36 = sorted(clin_df[(clin_df['Corrected gestational Age (weeks)'] >= 33) & \n",
    "                       (clin_df['Corrected gestational Age (weeks)'] < 37)].index)\n",
    "len(between_33_36)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "between_37_40 = sorted(clin_df[(clin_df['Corrected gestational Age (weeks)'] >= 37) & \n",
    "                       (clin_df['Corrected gestational Age (weeks)'] < 41)].index)\n",
    "len(between_37_40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "over_40 = sorted(clin_df[clin_df['Corrected gestational Age (weeks)'] >= 41].index)\n",
    "len(over_40)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Clinical problems"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# To analyse primary reason for transfer\n",
    "\n",
    "for index, row in clin_df.iterrows():\n",
    "    print(index, row['Gestational Age (weeks)'], row['Corrected gestational Age (weeks)'],\n",
    "        row['Pathology_English'], '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "writer = pd.ExcelWriter('%s/%s' % (DIR_WRITE, 'clinical_data_accelerometer_final.xlsx'))\n",
    "clin_df.to_excel(writer, 'cases')\n",
    "clin_df_stats.to_excel(writer, 'stats')\n",
    "writer.save()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. Remove the samples with  >46 week corrected gestational age from the ventilator and accelerometer data as well"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "combined_ventilated_1min = {key: value for key, value in combined_ventilated_1min.items() \n",
    "                               if key[1] in patients}\n",
    "\n",
    "accelero_ventilated = {key: value for key, value in accelero_ventilated.items() \n",
    "                               if key[1] in patients}\n",
    "\n",
    "accelero_ventilated_low_pass = {key: value for key, value in accelero_ventilated_low_pass.items() \n",
    "                               if key[1] in patients}\n",
    "\n",
    "accelero_ventilated_high_pass = {key: value for key, value in accelero_ventilated_high_pass.items() \n",
    "                               if key[1] in patients}\n",
    "\n",
    "ventilator_mode = {key: value for key, value in ventilator_mode.items() if key[1] in patients}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(combined_ventilated_1min),  len(ventilator_mode)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(accelero_ventilated), len(accelero_ventilated_low_pass), len(accelero_ventilated_high_pass),"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7. Calculate recording durations of the final dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "recording_durations_final = {}\n",
    "\n",
    "for recording in combined_ventilated_1min:\n",
    "    recording_durations_final[recording] = len(combined_ventilated_1min[recording])\n",
    "    \n",
    "recording_durations_final = DataFrame(Series(recording_durations_final), columns = ['duration (minutes)'])\n",
    "recording_durations_final.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "recording_durations_final.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "recording_durations_final.sum(), recording_durations_final.sum() / 60, recording_durations_final.sum() / 3600, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "writer = pd.ExcelWriter('%s/%s' % (DIR_WRITE, 'recording_durations_final.xlsx'))\n",
    "recording_durations_final.to_excel(writer, 'final')\n",
    "recording_durations_final.describe().to_excel(writer, 'stats')\n",
    "writer.save()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8. Analyse ventilator modes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Only consider those recordings where the ventilation mode did not change during the recording"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "for key, value in ventilator_mode.items():\n",
    "    print(key, value['Ventilator_mode'].unique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simv_only = {key[1] for key, value in ventilator_mode.items()\n",
    "           if 'SIMV' in value['Ventilator_mode'].unique() and len(value['Ventilator_mode'].unique()) == 1}\n",
    "len(simv_only)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sippv_only = {key[1] for key, value in ventilator_mode.items()\n",
    "           if 'SIPPV' in value['Ventilator_mode'].unique() and len(value['Ventilator_mode'].unique()) == 1}\n",
    "len(sippv_only)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simvpsv_only = {key[1] for key, value in ventilator_mode.items()\n",
    "           if 'SIMVPSV' in value['Ventilator_mode'].unique() and len(value['Ventilator_mode'].unique()) == 1}\n",
    "len(simvpsv_only)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many recording has multiple modes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "multiple_mode = {key[1] for key, value in ventilator_mode.items()\n",
    "           if len(value['Ventilator_mode'].unique()) > 1}\n",
    "len(multiple_mode)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many times was VG on for part or all of the recording"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "has_VG = set()\n",
    "\n",
    "for recording in combined_ventilated_1min:\n",
    "    if 'VG_set_kg_mean' in combined_ventilated_1min[recording].columns:\n",
    "        has_VG.add(recording)\n",
    "        # print(recording)\n",
    "\n",
    "len(has_VG)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "has_noVG = set()\n",
    "\n",
    "for recording in combined_ventilated_1min:\n",
    "    if 'VG_set_kg_mean' not in combined_ventilated_1min[recording].columns:\n",
    "        has_noVG.add(recording)\n",
    "    \n",
    "    elif sum(combined_ventilated_1min[recording]['VG_set_kg_mean'].isnull()):\n",
    "        has_noVG.add(recording)\n",
    "    \n",
    "len(has_noVG)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 9. Generate explanatory graphs for acceleration vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "filetype = 'jpg'\n",
    "dpi  = 200\n",
    "\n",
    "fig = plt.figure(figsize = [6,6])\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.quiver(0, 0, 0, 3.5, 2.8, 3.5, linewidth = 2, color = 'red', linestyle = '-', arrow_length_ratio = 0.1)\n",
    "ax.quiver(0, 0, 0, 3.5, 0, 0, linewidth = 2, color = 'black', arrow_length_ratio = 0.1)\n",
    "ax.quiver(0, 0, 0, 0, 2.8, 0, linewidth = 2, color = 'black', arrow_length_ratio = 0.1)\n",
    "ax.quiver(0, 0, 0, 0, 0, 3.5, linewidth = 2, color = 'black', arrow_length_ratio = 0.1)\n",
    "\n",
    "ax.quiver(0, 2.8, 0, 0, 0, 3.5, linewidth = 2, linestyle = 'dotted', color = 'black', arrow_length_ratio = 0.1,)\n",
    "ax.quiver(0, 2.8, 3.5, 3.5, 0, 0, linewidth = 2, linestyle = 'dotted', color = 'black', arrow_length_ratio = 0.1,)\n",
    "\n",
    "ax.set_xlim([0, 5]); ax.set_ylim([0,4]); ax.set_zlim([0, 4])\n",
    "ax.set_xlabel('X (m/sec$^2$)'); ax.set_ylabel('Y (m/sec$^2$)'); ax.set_zlabel('Z (m/sec$^2$)')\n",
    "ax.view_init(elev=20., azim = 200)\n",
    "\n",
    "fig.savefig('%s/%s.%s' % (DIR_WRITE, 'accel_explanation', filetype), \n",
    "    dpi = dpi, format = filetype, bbox_inches='tight', pad_inches=0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 10. Generate and export histograms of accelerations using the original accelerometer data (~100 Hz sampling)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A. Acceleration vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def acceleration(rec, direction, bins = 50, xlim = [-3, 3], ylim = [0, 100000], logyscale = False,\n",
    "               color = 'blue', show = True, save = False, dpi = 200, filetype = 'jpg'):\n",
    "    \n",
    "    fig, ax = plt.subplots(figsize = [4, 4])\n",
    "    \n",
    "    accelero_ventilated[rec][direction].hist(ax = ax, bins=bins, color = color)\n",
    "    ax.set_xlabel('m/sec$^2$', size = 12)\n",
    "    ax.set_ylabel('count', size = 12)\n",
    "    ax.set_title(direction, size = 14)\n",
    "    ax.set_xlim(*xlim)\n",
    "    ax.set_ylim(*ylim)\n",
    "    \n",
    "    if logyscale:\n",
    "        ax.set_yscale('log')\n",
    "    \n",
    "    name = rec[1] + '_log' if logyscale else rec[1]\n",
    "    \n",
    "    if save:       \n",
    "        fig.savefig('%s/%s_%s_%s.%s' % (DIR_WRITE, 'hist', name, direction, filetype), dpi = dpi, \n",
    "           format = filetype, bbox_inches='tight', pad_inches=0.1)\n",
    "        \n",
    "    if not show:\n",
    "        plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acceleration(('default__393.txt', 'AL000628'), 'X', bins = 20, logyscale = True, xlim = [-8,8],\n",
    "             ylim = [1, 1000000], save = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "acceleration(('default__393.txt', 'AL000628'), 'Y', bins = 50, logyscale = True, xlim = [-9,9],\n",
    "             ylim = [1, 1000000], save = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "acceleration(('default__393.txt', 'AL000628'), 'Z', bins = 30, logyscale = True, xlim = [-10,12],\n",
    "             ylim = [1, 1000000], save = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Absolute value of acceleration vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def acceleration_abs(rec, direction, bins = 50, xlim = [0, 8], ylim = [0, 100000], logyscale = False,\n",
    "               color = 'blue', show = True, save = False, dpi = 200, filetype = 'jpg'):\n",
    "    \n",
    "    fig, ax = plt.subplots(figsize = [4, 4])\n",
    "    \n",
    "    accelero_ventilated[rec]['%s_abs' % direction].hist(ax = ax, bins=bins, color = color)\n",
    "    ax.set_xlabel('m/sec$^2$', size = 12)\n",
    "    ax.set_ylabel('count', size = 12)\n",
    "    ax.set_title(direction, size = 14)\n",
    "    ax.set_xlim(*xlim)\n",
    "    ax.set_ylim(*ylim)\n",
    "    \n",
    "    if logyscale:\n",
    "        ax.set_yscale('log')\n",
    "    \n",
    "    name = rec[1] + '_log' if logyscale else rec[1]\n",
    "    \n",
    "    if save:       \n",
    "        fig.savefig('%s/%s_%s_%s.%s' % (DIR_WRITE, 'hist_abs', name, direction, filetype), dpi = dpi, \n",
    "        format = filetype, bbox_inches='tight', pad_inches=0.1)\n",
    "        \n",
    "    if not show:\n",
    "        plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "acceleration_abs(('default__393.txt', 'AL000628'), 'X', bins = 20, logyscale = True, xlim = [0,7],\n",
    "             ylim = [1, 1000000], save = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "acceleration_abs(('default__393.txt', 'AL000628'), 'Y', bins = 20, logyscale = True, xlim = [0,10],\n",
    "             ylim = [1, 1000000], save = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "acceleration_abs(('default__393.txt', 'AL000628'), 'Z', bins = 20, logyscale = True, xlim = [0,12],\n",
    "             ylim = [1, 1000000], save = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### X` (front-back), `Y` (left-right) and `Z` (up-down) acceleration vectors together"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def acceleration_XYZ(rec, bins = 50, xlim = [-3, 3], logyscale = False,\n",
    "               color = 'blue', show = True, save = False, dpi = 200, filetype = 'jpg'):\n",
    "    \n",
    "    fig, axes = plt.subplots(nrows=1, ncols=3, figsize = [12, 4], sharex=False, sharey=True)\n",
    "    \n",
    "    for par, ax in zip(['X', 'Y', 'Z'], axes.ravel()):\n",
    "        \n",
    "        data = accelero_ventilated[rec][par]\n",
    "        # Remove outlier data outside the X axis scale\n",
    "        data = data[(data > xlim[0]) & (data < xlim[1])]\n",
    "        data.hist(ax = ax, bins=bins, color = color)\n",
    "        ax.set_xlabel('m/sec$^2$', size = 14)\n",
    "        ax.set_ylabel('count', size = 14)\n",
    "        ax.set_title(par, size = 16)\n",
    "        ax.set_xlim(*xlim)\n",
    "    \n",
    "        if logyscale:\n",
    "            ax.set_yscale('log')\n",
    "    \n",
    "    name = rec[1] + '_log' if logyscale else rec[1]\n",
    "    \n",
    "    if save:\n",
    "        fig.savefig('%s/%s_%s.%s' % (DIR_WRITE, 'hist_XYZ', name, filetype), dpi = dpi, \n",
    "        format = filetype, bbox_inches='tight', pad_inches=0.1)\n",
    "        \n",
    "    if not show:\n",
    "        plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Example\n",
    "acceleration_XYZ(('default__393.txt', 'AL000628'), bins = 50, save = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Example\n",
    "acceleration_XYZ(('default__393.txt', 'AL000628'), bins = 50, xlim = (-10, 10), logyscale=True, save = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def accelerometer_XYZ_all(recs, bins = 50, xlim = [-3, 3], logyscale = False, color = 'blue', \n",
    "    name = 'hist_XYZ_all', save = True, show = False, dpi = 200, filetype = 'jpg'):\n",
    "    '''\n",
    "    Generates a multiplot graph with the acceleration of all recordings with each recording in a separate\n",
    "    row. This results in a very large vertical side length of the image\n",
    "    '''\n",
    "    fig, axes = plt.subplots(nrows=len(recs), ncols=3, figsize = [12, len(recs) * 4], sharex=False, sharey=False)\n",
    "    fig.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.5, hspace=0.5)\n",
    "    \n",
    "    ax_ravelled = axes.ravel()\n",
    "    i = 0\n",
    "    for rec in recs:\n",
    "        for par in ['X', 'Y', 'Z']:\n",
    "            ax = ax_ravelled[i]\n",
    "            data = accelero_ventilated[rec][par]\n",
    "            # Remove outlier data outside the X axis scale\n",
    "            data = data[(data > xlim[0]) & (data < xlim[1])]\n",
    "            data.hist(ax = ax , bins=bins, color = color)\n",
    "            ax.set_xlabel('m/sec$^2$', size = 10)\n",
    "            ax.set_ylabel('count', size = 10)\n",
    "            ax.set_title('%s - %s' % (rec[1], par), size = 10)\n",
    "            ax.set_xlim([*xlim])\n",
    "            \n",
    "            if logyscale:\n",
    "                ax.set_yscale('log')\n",
    "            \n",
    "            i+=1\n",
    "            if i % 10 == 0:\n",
    "                print(datetime.now(), 'Created %s graphs' %i)\n",
    "            \n",
    "    name = name + '_log' if logyscale else name\n",
    "    \n",
    "    if save:\n",
    "        fig.savefig('%s/%s.%s' % (DIR_WRITE, name, filetype), dpi = dpi, \n",
    "        format = filetype, bbox_inches='tight', pad_inches=0.1)\n",
    "        \n",
    "    if not show:\n",
    "        plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example plot\n",
    "accelerometer_XYZ_all([('default__369.txt', 'AL000611'), ('default__393.txt', 'AL000628')], \n",
    "                      show = True, save= False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example plot\n",
    "accelerometer_XYZ_all([('default__369.txt', 'AL000611'), ('default__393.txt', 'AL000628')], \n",
    "                      xlim = (-10, 10), logyscale=True, show = True, save= False)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "scrolled": false
   },
   "source": [
    "%%time\n",
    "\n",
    "accelerometer_XYZ_all(sorted(accelero_ventilated.keys(), key = lambda x: int(x[0][9:-4])), \n",
    "    dpi = 100, filetype = 'pdf', )"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "%%time\n",
    "\n",
    "accelerometer_XYZ_all(sorted(accelero_ventilated.keys(), key = lambda x: int(x[0][9:-4])), \n",
    "                xlim = (-10, 10), dpi = 100, filetype = 'pdf', logyscale=True,)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def accelerometer_XYZ_all_2(recs, bins = 50, xlim = [-3, 3], logyscale = False, color = 'blue',\n",
    "    name = 'hist_XYZ_all', save = True, show = False, dpi = 200, filetype = 'jpg'):\n",
    "    '''\n",
    "    Generates a multiplot graph with the acceleration of all recordings with each recording in a separate\n",
    "    row. This results in a very large vertical side length of the image\n",
    "    '''\n",
    "    fig, axes = plt.subplots(nrows=len(recs) // 4+1, ncols=12, figsize = [60, 60], sharex=False, sharey=False)\n",
    "    fig.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.6, hspace=0.6)\n",
    "    \n",
    "    ax_ravelled = axes.ravel()\n",
    "    i = 0\n",
    "    for rec in recs:\n",
    "        for par in ['X', 'Y', 'Z']:\n",
    "            ax = ax_ravelled[i]\n",
    "            data = accelero_ventilated[rec][par]\n",
    "            # Remove outlier data outside the X axis scale\n",
    "            data = data[(data > xlim[0]) & (data < xlim[1])]\n",
    "            data.hist(ax = ax , bins=bins, color = color)\n",
    "            ax.set_xlabel('m/sec$^2$', size = 10)\n",
    "            ax.set_ylabel('count', size = 10)\n",
    "            ax.set_title('%s - %s' % (rec[1], par), size = 10)\n",
    "            ax.set_xlim([*xlim])\n",
    "            \n",
    "            if logyscale:\n",
    "                ax.set_yscale('log')\n",
    "            \n",
    "            i+=1\n",
    "            if i % 10 == 0:\n",
    "                print(datetime.now(), 'Created %s graphs' %i)\n",
    "            \n",
    "    name = name + '_log' if logyscale else name\n",
    "    \n",
    "    if save:\n",
    "        fig.savefig('%s/%s.%s' % (DIR_WRITE, name, filetype), dpi = dpi, \n",
    "        format = filetype, bbox_inches='tight', pad_inches=0.1)\n",
    "        \n",
    "    if not show:\n",
    "        plt.close()"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "%%time\n",
    "\n",
    "accelerometer_XYZ_all_2(sorted(accelero_ventilated.keys(), key = lambda x: int(x[0][9:-4])), \n",
    "    dpi = 100, name = 'hist_XYZ_all_2', filetype = 'pdf', )"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "%%time\n",
    "\n",
    "accelerometer_XYZ_all_2(sorted(accelero_ventilated.keys(), key = lambda x: int(x[0][9:-4])), \n",
    "   xlim = (-10, 10), dpi = 100, name = 'hist_XYZ_all_2', filetype = 'pdf', logyscale=True,)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### B. Length (Euclidean (L2) norm) of the acceleration vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def accelerometer_euclidean_norm(rec, bins = 50, xlim = [0, 4], ylim = [0, 40000],\n",
    "        logyscale = False, color = 'blue', show = True, save = False, dpi = 200, filetype = 'jpg'):\n",
    "    \n",
    "    fig, ax = plt.subplots(figsize = [4, 4])\n",
    "    \n",
    "    data = accelero_ventilated[rec]['length']\n",
    "    # Remove outlier data outside the X axis scale\n",
    "    data = data[(data > xlim[0]) & (data < xlim[1])]\n",
    "    data.hist(ax = ax, bins=bins, color = color)\n",
    "    ax.set_xlabel('m/sec$^2$', size = 12)\n",
    "    ax.set_ylabel('count', size = 12)\n",
    "    ax.set_title('Euclidean length', size = 14)\n",
    "    ax.set_xlim(*xlim)\n",
    "    ax.set_ylim(*ylim)\n",
    "    \n",
    "    if logyscale:\n",
    "        ax.set_yscale('log')\n",
    "    \n",
    "    name = rec[1] + '_log' if logyscale else rec[1]\n",
    "    \n",
    "    if save:       \n",
    "        fig.savefig('%s/%s_%s.%s' % (DIR_WRITE, 'hist_euclidean', name, filetype), dpi = dpi, \n",
    "            format = filetype, bbox_inches='tight', pad_inches=0.1)\n",
    "        \n",
    "    if not show:\n",
    "        plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Example\n",
    "accelerometer_euclidean_norm(('default__393.txt', 'AL000628'), bins = 50, save = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "accelerometer_euclidean_norm(('default__393.txt', 'AL000628'), bins = 30,\n",
    "                             logyscale=True, xlim = (0,10), ylim = (1, 1000000), save = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 11. Analyse distributions of acceleration vector components in the whole recordings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Median of absolute values of acceleration vectors and the euclidean length in each recording overall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "to_keep = ['X_abs', 'Y_abs', 'Z_abs', 'length']\n",
    "\n",
    "accelero_stats_median = {}\n",
    "for recording in sorted(accelero_ventilated.keys()):\n",
    "    accelero_stats_median[recording] = accelero_ventilated[recording][to_keep].median()\n",
    "    \n",
    "accelero_stats_median = DataFrame(accelero_stats_median).T\n",
    "\n",
    "accelero_stats_median.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "accelero_stats_median_describe = \\\n",
    "    accelero_stats_median.describe(percentiles = [0.01, 0.05, 0.1, 0.25, 0.5, 0.75, 0.95, 0.99])\n",
    "accelero_stats_median_describe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "writer = pd.ExcelWriter('%s/%s' % (DIR_WRITE, 'median_acceleration_recordings.xlsx'))\n",
    "accelero_stats_median.to_excel(writer, 'median_accels')\n",
    "accelero_stats_median_describe.to_excel(writer, 'median_accel_stats')\n",
    "writer.save()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Distribution of the median absolute acceleration of each recording\n",
    "Medians are not normally distributed as it can be seen on the graphs; use non-parametric statistics: group medians and IQRs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "par = 'X_abs'\n",
    "data = accelero_stats_median[par]\n",
    "\n",
    "fig, ax = plt.subplots(figsize = [4,4])\n",
    "ax.hist(data, bins = 25, density = True, color = 'black', alpha = 0.7)    \n",
    "# Plot the PDF.\n",
    "xmin, xmax = ax.get_xlim()\n",
    "x = np.linspace(xmin, xmax, 100)\n",
    "p = stats.norm.pdf(x, data.mean(), data.std())\n",
    "ax.plot(x, p, linewidth=2, color = 'red')\n",
    "ax.set_xlabel('acceleration (m/sec$^2$)', size = 12)\n",
    "ax.set_ylabel('count', size = 12)\n",
    "ax.set_xlim(0, data.mean()*2)\n",
    "plt.title(par, fontsize = 12)\n",
    "plt.grid(True)\n",
    "\n",
    "fig.savefig('%s/%s.%s' % (DIR_WRITE, 'accel_%s_median_hist_with_pdf' % par, 'jpg' ),\n",
    "    dpi = 200, format = 'jpg', bbox_inches='tight', pad_inches=0.1,);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "par = 'Y_abs'\n",
    "data = accelero_stats_median[par]\n",
    "\n",
    "fig, ax = plt.subplots(figsize = [4,4])\n",
    "ax.hist(data, bins = 25, density = True, color = 'black', alpha = 0.7)    \n",
    "# Plot the PDF.\n",
    "xmin, xmax = ax.get_xlim()\n",
    "x = np.linspace(xmin, xmax, 100)\n",
    "p = stats.norm.pdf(x, data.mean(), data.std())\n",
    "ax.plot(x, p, linewidth=2, color = 'red')\n",
    "ax.set_xlabel('acceleration (m/sec$^2$)', size = 12)\n",
    "ax.set_ylabel('count', size = 12)\n",
    "ax.set_xlim(0, data.mean()*2)\n",
    "plt.title(par, fontsize = 12)\n",
    "plt.grid(True)\n",
    "\n",
    "fig.savefig('%s/%s.%s' % (DIR_WRITE, 'accel_%s_median_hist_with_pdf' % par, 'jpg' ),\n",
    "    dpi = 200, format = 'jpg', bbox_inches='tight', pad_inches=0.1,);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "par = 'Z_abs'\n",
    "data = accelero_stats_median['Z_abs']\n",
    "\n",
    "fig, ax = plt.subplots(figsize = [4,4])\n",
    "ax.hist(data, bins = 25, density = True, color = 'black', alpha = 0.7)    \n",
    "# Plot the PDF.\n",
    "xmin, xmax = ax.get_xlim()\n",
    "x = np.linspace(xmin, xmax, 100)\n",
    "p = stats.norm.pdf(x, data.mean(), data.std())\n",
    "ax.plot(x, p, linewidth=2, color = 'red')\n",
    "ax.set_xlabel('acceleration (m/sec$^2$)', size = 12)\n",
    "ax.set_ylabel('count', size = 12)\n",
    "ax.set_xlim(0, data.mean()*2)\n",
    "plt.title(par, fontsize = 12)\n",
    "plt.grid(True)\n",
    "\n",
    "fig.savefig('%s/%s.%s' % (DIR_WRITE, 'accel_Z_abs_median_hist_with_pdf', 'jpg'),\n",
    "    dpi = 200, format = 'jpg', bbox_inches='tight', pad_inches=0.1,);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "par = 'Euclidean_abs'\n",
    "data = accelero_stats_median['length']\n",
    "\n",
    "fig, ax = plt.subplots(figsize = [4,4])\n",
    "ax.hist(data, bins = 25, density = True, color = 'black', alpha = 0.7)    \n",
    "# Plot the PDF.\n",
    "xmin, xmax = ax.get_xlim()\n",
    "x = np.linspace(xmin, xmax, 100)\n",
    "p = stats.norm.pdf(x, data.mean(), data.std())\n",
    "ax.plot(x, p, linewidth=2, color = 'red')\n",
    "ax.set_xlabel('acceleration (m/sec$^2$)', size = 12)\n",
    "ax.set_ylabel('count', size = 12)\n",
    "ax.set_xlim(0, data.mean()*2)\n",
    "plt.title('Euclidean_abs', fontsize = 12)\n",
    "plt.grid(True)\n",
    "\n",
    "fig.savefig('%s/%s.%s' % (DIR_WRITE, 'accel_euclidean_median_hist_with_pdf', 'jpg'),\n",
    "    dpi = 200, format = 'jpg', bbox_inches='tight', pad_inches=0.1,);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "stats.mannwhitneyu(accelero_stats_median['X_abs'], accelero_stats_median['Y_abs'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "stats.mannwhitneyu(accelero_stats_median['X_abs'], accelero_stats_median['Z_abs'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize = (4, 4))\n",
    "dpi = 300\n",
    "filetype = 'jpg'\n",
    "xticklabels = ['X', 'Y', 'Z']\n",
    "\n",
    "medianprops = {'color': 'black', 'linewidth': 2}\n",
    "meanpointprops = {'marker':'D', 'markeredgecolor':'black', 'markerfacecolor':'black'}\n",
    "boxprops = {'color': 'black', 'linestyle': '-'}\n",
    "whiskerprops = {'color': 'black', 'linestyle': '-'}\n",
    "capprops = {'color': 'black', 'linestyle': '-'}\n",
    "flierprops = {'color': 'black', 'marker': '.'}\n",
    "\n",
    "plt.boxplot([accelero_stats_median['X_abs'], accelero_stats_median['Y_abs'], accelero_stats_median['Z_abs']], \n",
    "    widths = 0.3, whis = [5, 95], showfliers = False, showmeans = True, meanprops = meanpointprops, \n",
    "        medianprops=medianprops, boxprops=boxprops, whiskerprops=whiskerprops, capprops=capprops, \n",
    "            flierprops = flierprops)\n",
    "\n",
    "ax.hlines([1.1, 1.3], [1, 1], [2, 3])\n",
    "ax.vlines([1, 2, 1, 3,], [1.1, 1.1, 1.3, 1.3], [1.05, 1.05, 1.25, 1.25])\n",
    "ax.text(1.2, 1.15, 'p<0.001'); ax.text(1.55, 1.35, 'p<0.001')\n",
    "\n",
    "ax.set_xticklabels(xticklabels)\n",
    "ax.set_ylabel('Median acceleration of recordings \\n (m/sec$^2$)', size = 12)\n",
    "ax.set_ylim(0, 1.5)\n",
    "ax.tick_params(axis='both', which='major', labelsize=14)\n",
    "plt.grid(True)\n",
    "\n",
    "fig.savefig('%s/%s.%s' % (DIR_WRITE, 'X_Y_Z_median_boxplot', 'jpg'),\n",
    "    dpi = 200, format = 'jpg', bbox_inches='tight', pad_inches=0.1,);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 12. Combine individual DataFrames into one"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# combine DataFrames into one\n",
    "combined_ventilated_1min_all = pd.concat(combined_ventilated_1min, sort = True)\n",
    "len(combined_ventilated_1min_all)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove unnecessary columns\n",
    "\n",
    "columns_to_keep = ['FiO2_set_mean', 'FiO2_set_median', 'FiO2_set_sd', 'Leak_mean', 'Leak_median', 'Leak_sd', \n",
    "                   \n",
    "                   'MAP_mean', 'MAP_sd', 'MV_kg_mean', 'MV_kg_sd', 'MVresp_mean', 'MVresp_sd', \n",
    "                   'PIP_mean', 'PIP_sd', 'PIP_set_mean', 'PIP_set_sd',\n",
    "                   'RR_mean', 'RR_sd', 'RR_set_mean', 'RR_set_sd', 'VG_set_kg_mean', 'VG_set_kg_sd',\n",
    "                   'VTemand_kg_mean', 'VTemand_kg_sd',\n",
    "                   'VTemand_resp_kg_mean', 'VTemand_resp_kg_sd', 'VTespon_pat_kg_mean', 'VTespon_pat_kg_sd',\n",
    "                   'VTimand_kg_mean', 'VTimand_kg_sd', \n",
    "                   \n",
    "                   'Pdiff_VG_abs_mean', 'Pdiff_VG_abs_median', 'Pdiff_VG_abs_sd', \n",
    "                   'Pdiff_VG_mean', 'Pdiff_VG_median', 'Pdiff_VG_sd', \n",
    "                   'Pdiff_noVG_abs_mean', 'Pdiff_noVG_abs_median', 'Pdiff_noVG_abs_sd', \n",
    "                   'Pdiff_noVG_mean', 'Pdiff_noVG_median', 'Pdiff_noVG_sd',\n",
    "                   \n",
    "                   'RRdiff_abs_mean', 'RRdiff_abs_median', 'RRdiff_abs_sd', \n",
    "                   'RRdiff_mean', 'RRdiff_median', 'RRdiff_sd',\n",
    "                   \n",
    "                   'VTdiff_abs_mean', 'VTdiff_abs_median', 'VTdiff_abs_sd', \n",
    "                   'VTdiff_mean', 'VTdiff_median', 'VTdiff_sd',\n",
    "                   \n",
    "                   'X_pos_median', 'X_neg_median', 'Y_pos_median', 'Y_neg_median',\n",
    "                   'X_abs_median_low_pass', 'X_abs_median_high_pass', 'Y_abs_median_low_pass', \n",
    "                   'Y_abs_median_high_pass', 'Z_abs_median_low_pass', 'Z_abs_median_high_pass', \n",
    "                   'length_median_low_pass', 'length_median_high_pass']\n",
    "\n",
    "combined_ventilated_1min_all = combined_ventilated_1min_all[columns_to_keep]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "combined_ventilated_1min_all.index.rename(['accel_rec', 'patient', 'datetime'], inplace = True)\n",
    "combined_ventilated_1min_all.reset_index(inplace = True)\n",
    "combined_ventilated_1min_all.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 13. Analyse the correlation vibration  in X (front-back) or Y (side-to-side) direction with IQR in Z (vertical) direction\n",
    "\n",
    "- vibration: high-pass filter\n",
    "- sustained acceleration: low pass filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize = (6, 6))\n",
    "dpi = 300\n",
    "filetype = 'jpg'\n",
    "\n",
    "a = combined_ventilated_1min_all.dropna(subset = ['Z_abs_median_high_pass','X_abs_median_high_pass'])\n",
    "corr = stats.spearmanr(a['Z_abs_median_high_pass'], a['X_abs_median_high_pass'])\n",
    "\n",
    "ax.scatter(combined_ventilated_1min_all['Z_abs_median_high_pass'], \n",
    "           combined_ventilated_1min_all['X_abs_median_high_pass'], \n",
    "           marker = 'o', facecolor = 'blue', edgecolor = 'blue', alpha = 0.1)\n",
    "ax.text(0.1, 0.7, 'r: %.2f\\np<0.0001' % corr[0], fontsize = 14)\n",
    "ax.grid(True)\n",
    "ax.set_xlim(0, combined_ventilated_1min_all['Z_abs_median_high_pass'].max() + 0.05)\n",
    "ax.set_ylim(-0.02, combined_ventilated_1min_all['X_abs_median_high_pass'].max() + 0.05)\n",
    "ax.set_xlabel('vertical vibration (m/sec$^2$)', size = 14)\n",
    "ax.set_ylabel('front-back vibration (m/sec$^2$)', size = 14)\n",
    "ax.yaxis.set_major_formatter(ScalarFormatter())\n",
    "ax.tick_params(axis='both', which='major', labelsize=14)\n",
    "\n",
    "fig.savefig('%s/%s.%s' % (DIR_WRITE, 'Z_abs_high_pass_X_abs_high_pass', filetype),\n",
    "    dpi = dpi, format = filetype, bbox_inches='tight', pad_inches=0.1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize = (6, 6))\n",
    "dpi = 300\n",
    "filetype = 'jpg'\n",
    "\n",
    "a = combined_ventilated_1min_all.dropna(subset = ['Z_abs_median_high_pass','Y_abs_median_high_pass'])\n",
    "corr = stats.spearmanr(a['Z_abs_median_high_pass'], a['Y_abs_median_high_pass'])\n",
    "\n",
    "ax.scatter(combined_ventilated_1min_all['Z_abs_median_high_pass'], \n",
    "           combined_ventilated_1min_all['Y_abs_median_high_pass'], \n",
    "           marker = 'o', facecolor = 'blue', edgecolor = 'blue', alpha = 0.1)\n",
    "ax.text(0.1, 0.6, 'r: %.2f\\np<0.0001' % corr[0], fontsize = 14)\n",
    "ax.grid(True)\n",
    "ax.set_xlim(0, combined_ventilated_1min_all['Z_abs_median_high_pass'].max() + 0.05)\n",
    "ax.set_ylim(-0.02, combined_ventilated_1min_all['Y_abs_median_high_pass'].max() + 0.05)\n",
    "ax.set_xlabel('vertical vibration (m/sec$^2$)', size = 14)\n",
    "ax.set_ylabel('side-to-side vibration (m/sec$^2$)', size = 14)\n",
    "ax.yaxis.set_major_formatter(ScalarFormatter())\n",
    "ax.tick_params(axis='both', which='major', labelsize=14)\n",
    "\n",
    "fig.savefig('%s/%s.%s' % (DIR_WRITE, 'Z_abs_high_pass_Y_abs_high_pass', filetype),\n",
    "    dpi = dpi, format = filetype, bbox_inches='tight', pad_inches=0.1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 14. Visualize the distribution of the median vibration and sustained acceleration along the three axes of the acceleration vector\n",
    "\n",
    "1 minute medians are compared\n",
    "\n",
    "- vibration: hig pass filter\n",
    "- sustained acceleration: low pass filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "dpi = 300\n",
    "filetype = 'jpg'\n",
    "\n",
    "a = combined_ventilated_1min_all.dropna(subset = ['X_abs_median_low_pass','X_abs_median_high_pass'])\n",
    "#corr = stats.spearmanr(a['X_abs_median_low_pass'], a['X_abs_median_high_pass'])\n",
    "\n",
    "plot = sns.jointplot(x='X_abs_median_low_pass', y='X_abs_median_high_pass', \n",
    "    data = combined_ventilated_1min_all, kind='scatter', color='blue', alpha = 0.1, \n",
    "    height=6, ratio=5, space=0.2, dropna=True, xlim=[0, 2], ylim=[0,0.8])\n",
    "\n",
    "plot.ax_joint.hlines(0.15, 0, 2, color = 'black', linestyle = 'dashed')\n",
    "#plot.ax_joint.text(0.1, 0.7, 'r: %.2f\\np<0.0001' % corr[0], fontsize = 14)\n",
    "plot.ax_joint.grid('True')\n",
    "plot.ax_joint.set_xlabel('front-back sustained acceleration (m/sec$^2$)', size = 14)\n",
    "plot.ax_joint.set_ylabel('front-back vibration (m/sec$^2$)', size = 14)\n",
    "plot.ax_joint.tick_params(axis='both', which='major', labelsize=14)\n",
    "\n",
    "plt.savefig('%s/%s.%s' % (DIR_WRITE, 'X_abs_acceleration_X_abs_vibration_jointplot', filetype),\n",
    "    dpi = dpi,  format = filetype, bbox_inches='tight', pad_inches=0.1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dpi = 300\n",
    "filetype = 'jpg'\n",
    "\n",
    "a = combined_ventilated_1min_all.dropna(subset = ['Y_abs_median_low_pass','Y_abs_median_high_pass'])\n",
    "#corr = stats.spearmanr(a['Y_abs_median_low_pass'], a['Y_abs_median_high_pass'])\n",
    "\n",
    "plot = sns.jointplot(x='Y_abs_median_low_pass', y='Y_abs_median_high_pass', \n",
    "    data = combined_ventilated_1min_all, kind='scatter', color='blue', alpha = 0.1, \n",
    "    height=6, ratio=5, space=0.2, dropna=True, xlim=[0, 1.2], ylim=[0,0.8])\n",
    "\n",
    "plot.ax_joint.hlines(0.07, 0, 2, color = 'black', linestyle = 'dashed')\n",
    "#plot.ax_joint.text(0.1, 0.7, 'r: %.2f\\np<0.0001' % corr[0], fontsize = 14)\n",
    "plot.ax_joint.grid('True')\n",
    "plot.ax_joint.set_xlabel('side-to-side sustained acceleration (m/sec$^2$)', size = 14)\n",
    "plot.ax_joint.set_ylabel('side-to-side vibration (m/sec$^2$)', size = 14)\n",
    "plot.ax_joint.tick_params(axis='both', which='major', labelsize=14)\n",
    "\n",
    "plt.savefig('%s/%s.%s' % (DIR_WRITE, 'Y_abs_acceleration_Y_abs_vibration_jointplot', filetype),\n",
    "    dpi = dpi,  format = filetype, bbox_inches='tight', pad_inches=0.1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dpi = 300\n",
    "filetype = 'jpg'\n",
    "\n",
    "a = combined_ventilated_1min_all.dropna(subset = ['Z_abs_median_low_pass','Z_abs_median_high_pass'])\n",
    "#corr = stats.spearmanr(a['Z_abs_median_low_pass'], a['Z_abs_median_high_pass'])\n",
    "\n",
    "plot = sns.jointplot(x='Z_abs_median_low_pass', y='Z_abs_median_high_pass', \n",
    "    data = combined_ventilated_1min_all, kind='scatter', color='blue', alpha = 0.1, \n",
    "    height=6, ratio=5, space=0.2, dropna=True, xlim=[0, 1.2], ylim=[0, 1.2])\n",
    "\n",
    "plot.ax_joint.hlines(0.2, 0, 2, color = 'black', linestyle = 'dashed')\n",
    "#plot.ax_joint.text(0.02, 1, 'r: %.2f\\np<0.0001' % corr[0], fontsize = 14)\n",
    "plot.ax_joint.grid('True')\n",
    "plot.ax_joint.set_xlabel('vertical sustained acceleration (m/sec$^2$)', size = 14)\n",
    "plot.ax_joint.set_ylabel('vertical vibration (m/sec$^2$)', size = 14)\n",
    "plot.ax_joint.tick_params(axis='both', which='major', labelsize=14)\n",
    "\n",
    "plt.savefig('%s/%s.%s' % (DIR_WRITE, 'Z_abs_acceleration_Z_abs_vibration_jointplot', filetype),\n",
    "    dpi = dpi,  format = filetype, bbox_inches='tight', pad_inches=0.1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dpi = 300\n",
    "filetype = 'jpg'\n",
    "\n",
    "a = combined_ventilated_1min_all.dropna(subset = ['length_median_low_pass','length_median_high_pass'])\n",
    "corr = stats.spearmanr(a['length_median_low_pass'], a['length_median_high_pass'])\n",
    "\n",
    "plot = sns.jointplot(x='length_median_low_pass', y='length_median_high_pass', \n",
    "    data = combined_ventilated_1min_all, kind='scatter', color='blue', alpha = 0.1, \n",
    "    height=6, ratio=5, space=0.2, dropna=True, xlim=[0, 2], ylim=[0, 2])\n",
    "\n",
    "plot.ax_joint.hlines(0.3, 0, 2, color = 'black', linestyle = 'dashed')\n",
    "plot.ax_joint.text(0.1, 1.7, 'r: %.2f\\np<0.0001' % corr[0], fontsize = 14)\n",
    "plot.ax_joint.grid('True')\n",
    "plot.ax_joint.set_xlabel('Euclidean sustained acceleration (m/sec$^2$)', size = 14)\n",
    "plot.ax_joint.set_ylabel('Euclidean vibration (m/sec$^2$)', size = 14)\n",
    "plot.ax_joint.tick_params(axis='both', which='major', labelsize=14)\n",
    "\n",
    "plt.savefig('%s/%s.%s' % (DIR_WRITE, 'length_acceleration_length_vibration_jointplot', filetype),\n",
    "    dpi = dpi,  format = filetype, bbox_inches='tight', pad_inches=0.1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 15. Identify periods of vibration and sustained acceleration\n",
    "\n",
    "- Periods of vibration: 1-min-median of the Euclidean length of the acceleration vector after passing through a high-pass filter is >=0.3m/sec2 (high vibration)\n",
    "\n",
    "- Periods of 'sustained' acceleration: 1-min-median of the Euclidean length of the acceleration vector after passing through a high-pass filter is <0.3m/sec2 low vibration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "low_1_vibr = {}; high_1_vibr = {}\n",
    "\n",
    "for patient in combined_ventilated_1min_all['patient'].unique():\n",
    "    patient_frame = combined_ventilated_1min_all[combined_ventilated_1min_all['patient'] == patient]\n",
    "    \n",
    "    # Only consider minutes when there was at least 0.3 m/sec2 Euclidean vibration\n",
    "    patient_frame = patient_frame[patient_frame['length_median_high_pass'] >=0.3]\n",
    "    \n",
    "    low_1_vibr[patient] = patient_frame.sort_values('length_median_high_pass', ascending = True)[:1]\n",
    "    high_1_vibr[patient] = patient_frame.sort_values('length_median_high_pass', ascending = True)[-1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "low_1_accel = {}; high_1_accel = {}\n",
    "\n",
    "for patient in combined_ventilated_1min_all['patient'].unique():\n",
    "    patient_frame = combined_ventilated_1min_all[combined_ventilated_1min_all['patient'] == patient]\n",
    "    \n",
    "    # Only consider minutes when there was less than 0.3 m/sec2 Euclidean vibration\n",
    "    patient_frame = patient_frame[patient_frame['length_median_high_pass'] <0.3]\n",
    "    \n",
    "    # Only include those recordings where there is at least two minutes to compare\n",
    "    if len(patient_frame) >=2:\n",
    "        low_1_accel[patient] = patient_frame.sort_values('length_median_low_pass', ascending = True)[:1]\n",
    "        high_1_accel[patient] = patient_frame.sort_values('length_median_low_pass', ascending = True)[-1:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 16. Compare periods with vibration with periods of no acceleration "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### A. All minutes (irrespective of ventilation modes and patient characteristics)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def two_paired_test(arr_1, arr_2, par):\n",
    "    return {     'parameter': par,\n",
    "                 'low_mean':  round(arr_1[par].mean(), 2),\n",
    "                 'low_median': round(arr_1[par].median(), 2),\n",
    "                 'low_sd':    round(arr_1[par].std(), 2),\n",
    "                 'low_min': round(arr_1[par].min(), 2),\n",
    "                 'low_max': round(arr_1[par].max(), 2),\n",
    "                 'high_mean': round(arr_2[par].mean(), 2),\n",
    "                 'high_median': round(arr_2[par].median(), 2),\n",
    "                 'high_sd':   round(arr_2[par].std(), 2),\n",
    "                 'high_min': round(arr_2[par].min(), 2),\n",
    "                 'high_max': round(arr_2[par].max(), 2),\n",
    "                 'statistic': round(stats.ttest_rel(arr_1[par], arr_2[par])[0], 2),\n",
    "                 'p-value':   stats.ttest_rel(arr_1[par], arr_2[par])[1],}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "low_1_vibr_frme = pd.concat(low_1_vibr).reset_index(level = 1)\n",
    "low_1_vibr_frme.head();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "high_1_vibr_frme = pd.concat(high_1_vibr).reset_index(level = 1)\n",
    "high_1_vibr_frme.head();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "columns_to_analyse_1 = ['length_median_high_pass',]\n",
    "\n",
    "comparison_stats_1 = []\n",
    "for column in columns_to_analyse_1:\n",
    "    comparison_stats_1.append(two_paired_test(low_1_vibr_frme, high_1_vibr_frme, column))\n",
    "\n",
    "comparison_stats_frme_1 = DataFrame(comparison_stats_1, index= columns_to_analyse_1)\n",
    "comparison_stats_frme_1;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "columns_to_analyse_2 = [ 'FiO2_set_mean',   'Leak_median', 'MV_kg_mean', 'PIP_mean', 'VTemand_kg_mean', ]\n",
    "\n",
    "comparison_stats_2 = []\n",
    "\n",
    "for column in columns_to_analyse_2:\n",
    "    comparison_stats_2.append(two_paired_test(low_1_vibr_frme, high_1_vibr_frme, column))\n",
    "\n",
    "comparison_stats_frme_2 = DataFrame(comparison_stats_2, index= columns_to_analyse_2)\n",
    "comparison_stats_frme_2;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "comparison_stats_frme = pd.concat([comparison_stats_frme_1, comparison_stats_frme_2])\n",
    "comparison_stats_frme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "writer = pd.ExcelWriter('%s/%s' % (DIR_WRITE, 'vibr_stats_all_recs.xlsx'))\n",
    "comparison_stats_frme.to_excel(writer, 'stats')\n",
    "writer.save()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### B. VG ventilation\n",
    "\n",
    "During these minutes the set tidal volume `VG_set_kg_mean` has a specific value different from np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(len(combined_ventilated_1min_all), \n",
    " sum(combined_ventilated_1min_all['VG_set_kg_mean'].notnull()),\n",
    " sum(combined_ventilated_1min_all['VG_set_kg_mean'].isnull()), \n",
    " sum(combined_ventilated_1min_all['VG_set_kg_mean'] == 0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VG = combined_ventilated_1min_all[combined_ventilated_1min_all['VG_set_kg_mean'].notnull()].copy()\n",
    "len(VG)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(VG['accel_rec'].unique()), len(VG['patient'].unique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "low_1_VG = {}; high_1_VG = {}\n",
    "\n",
    "for patient in VG['patient'].unique():\n",
    "    patient_frame_VG = VG[VG['patient'] == patient]\n",
    "    \n",
    "    # Only consider minutes when there was at least 0.3 m/sec2 vibration\n",
    "    patient_frame_VG = patient_frame_VG[patient_frame_VG['length_median_high_pass'] >=0.3]\n",
    "    \n",
    "    low_1_VG[patient] = patient_frame_VG.sort_values('length_median_high_pass', ascending = True)[:1]\n",
    "    high_1_VG[patient] = patient_frame_VG.sort_values('length_median_high_pass', ascending = True)[-1:]\n",
    "\n",
    "low_1_frme_VG = pd.concat(low_1_VG)\n",
    "high_1_frme_VG = pd.concat(high_1_VG)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "columns_to_analyse_1_VG = ['length_median_high_pass',]\n",
    "\n",
    "comparison_stats_1_VG = []\n",
    "\n",
    "for column in columns_to_analyse_1_VG:\n",
    "    comparison_stats_1_VG.append(two_paired_test(low_1_frme_VG, high_1_frme_VG, column))\n",
    "\n",
    "comparison_stats_frme_1_VG = DataFrame(comparison_stats_1_VG, index = columns_to_analyse_1_VG)\n",
    "comparison_stats_frme_1_VG;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "columns_to_analyse_2_VG = ['PIP_mean', 'VTemand_kg_mean', 'VTdiff_abs_median',  'Pdiff_VG_abs_median', ]\n",
    "\n",
    "comparison_stats_2_VG = []\n",
    "for column in columns_to_analyse_2_VG:\n",
    "    comparison_stats_2_VG.append(two_paired_test(low_1_frme_VG, high_1_frme_VG, column))\n",
    "\n",
    "comparison_stats_frme_2_VG = DataFrame(comparison_stats_2_VG, index = columns_to_analyse_2_VG)\n",
    "comparison_stats_frme_2_VG;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "comparison_stats_frme_VG = pd.concat([comparison_stats_frme_1_VG, comparison_stats_frme_2_VG])\n",
    "comparison_stats_frme_VG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "writer = pd.ExcelWriter('%s/%s' % (DIR_WRITE, 'vibr_stats_VG.xlsx'))\n",
    "comparison_stats_frme_VG.to_excel(writer, 'stats')\n",
    "writer.save()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### C. Minutes with pressure limited (non-VG) ventilation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Plim = combined_ventilated_1min_all[combined_ventilated_1min_all['VG_set_kg_mean'].isnull()].copy()\n",
    "len(Plim)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "len(Plim['accel_rec'].unique()), len(Plim['patient'].unique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "low_1_Plim= {}; high_1_Plim = {}\n",
    "\n",
    "for patient in Plim['patient'].unique():\n",
    "    patient_frame_Plim = Plim[Plim['patient'] == patient]\n",
    "    \n",
    "    # Only consider minutes when there was at least 0.3 m/sec2 vibration\n",
    "    patient_frame_Plim = patient_frame_Plim[patient_frame_Plim['length_median_high_pass'] >=0.3]\n",
    "    \n",
    "    low_1_Plim[patient] = patient_frame_Plim.sort_values('length_median_high_pass', ascending = True)[:1]\n",
    "    high_1_Plim[patient] = patient_frame_Plim.sort_values('length_median_high_pass', ascending = True)[-1:]\n",
    "\n",
    "low_1_frme_Plim = pd.concat(low_1_Plim)\n",
    "high_1_frme_Plim = pd.concat(high_1_Plim)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "columns_to_analyse_1_Plim = ['length_median_high_pass',]\n",
    "\n",
    "comparison_stats_1_Plim = []\n",
    "\n",
    "for column in columns_to_analyse_1_Plim:\n",
    "    comparison_stats_1_Plim.append(two_paired_test(low_1_frme_Plim, high_1_frme_Plim, column))\n",
    "\n",
    "comparison_stats_frme_1_Plim = DataFrame(comparison_stats_1_Plim, index = columns_to_analyse_1_Plim)\n",
    "comparison_stats_frme_1_Plim;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "columns_to_analyse_2_Plim = ['PIP_mean', 'VTemand_kg_mean', 'Pdiff_noVG_abs_median',]\n",
    "\n",
    "comparison_stats_2_Plim = []\n",
    "\n",
    "for column in columns_to_analyse_2_Plim:\n",
    "    comparison_stats_2_Plim.append(two_paired_test(low_1_frme_Plim, high_1_frme_Plim, column))\n",
    "\n",
    "comparison_stats_frme_2_Plim = DataFrame(comparison_stats_2_Plim, index = columns_to_analyse_2_Plim)\n",
    "comparison_stats_frme_2_Plim;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "comparison_stats_frme_Plim = pd.concat([comparison_stats_frme_1_Plim, comparison_stats_frme_2_Plim])\n",
    "comparison_stats_frme_Plim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "writer = pd.ExcelWriter('%s/%s' % (DIR_WRITE, 'vibr_stats_Plim.xlsx'))\n",
    "comparison_stats_frme_Plim.to_excel(writer, 'stats')\n",
    "writer.save()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### D. Minutes with SIPPV ventilation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sippv = combined_ventilated_1min_all[combined_ventilated_1min_all['patient'].isin(sippv_only)].copy()\n",
    "len(sippv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "len(sippv['accel_rec'].unique()), len(sippv['patient'].unique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "low_1_sippv= {}; high_1_sippv = {}\n",
    "\n",
    "for patient in sippv['patient'].unique():\n",
    "    patient_frame_sippv = sippv[sippv['patient'] == patient]\n",
    "    \n",
    "    # Only consider minutes when there was at least 0.3 m/sec2 vibration\n",
    "    patient_frame_sippv = patient_frame_sippv[patient_frame_sippv['length_median_high_pass'] >=0.3]\n",
    "   \n",
    "    low_1_sippv[patient] = patient_frame_sippv.sort_values('length_median_high_pass', ascending = True)[:1]\n",
    "    high_1_sippv[patient] = patient_frame_sippv.sort_values('length_median_high_pass', ascending = True)[-1:]\n",
    "\n",
    "low_1_frme_sippv = pd.concat(low_1_sippv)\n",
    "high_1_frme_sippv = pd.concat(high_1_sippv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "columns_to_analyse_1_sippv = ['length_median_high_pass',]\n",
    "\n",
    "comparison_stats_1_sippv = []\n",
    "\n",
    "for column in columns_to_analyse_1_sippv:\n",
    "    comparison_stats_1_sippv.append(two_paired_test(low_1_frme_sippv, high_1_frme_sippv, column))\n",
    "\n",
    "comparison_stats_frme_1_sippv = DataFrame(comparison_stats_1_sippv, index = columns_to_analyse_1_sippv)\n",
    "comparison_stats_frme_1_sippv;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "columns_to_analyse_2_sippv = ['RR_mean',  'RRdiff_abs_median', ]\n",
    "\n",
    "comparison_stats_2_sippv = []\n",
    "\n",
    "for column in columns_to_analyse_2_sippv:\n",
    "    comparison_stats_2_sippv.append(two_paired_test(low_1_frme_sippv, high_1_frme_sippv, column))\n",
    "\n",
    "comparison_stats_frme_2_sippv = DataFrame(comparison_stats_2_sippv, index = columns_to_analyse_2_sippv)\n",
    "comparison_stats_frme_2_sippv;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "comparison_stats_frme_sippv = pd.concat([comparison_stats_frme_1_sippv, comparison_stats_frme_2_sippv])\n",
    "comparison_stats_frme_sippv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "writer = pd.ExcelWriter('%s/%s' % (DIR_WRITE, 'vibr_stats_sippv.xlsx'))\n",
    "comparison_stats_frme_sippv.to_excel(writer, 'stats')\n",
    "writer.save()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### E. Minutes with SIMV ventilation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "simv = combined_ventilated_1min_all[combined_ventilated_1min_all['patient'].isin(simv_only)].copy()\n",
    "len(simv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(simv['accel_rec'].unique()), len(simv['patient'].unique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "low_1_simv= {}; high_1_simv = {}\n",
    "\n",
    "for patient in simv['patient'].unique():\n",
    "    patient_frame_simv = simv[simv['patient'] == patient]\n",
    "    \n",
    "    # Only consider minutes when there was at least 0.3 m/sec2 vibration\n",
    "    patient_frame_simv = patient_frame_simv[patient_frame_simv['length_median_high_pass'] >=0.3]\n",
    "    \n",
    "    low_1_simv[patient] = patient_frame_simv.sort_values('length_median_high_pass', ascending = True)[:1]\n",
    "    high_1_simv[patient] = patient_frame_simv.sort_values('length_median_high_pass', ascending = True)[-1:]\n",
    "\n",
    "low_1_frme_simv = pd.concat(low_1_simv)\n",
    "high_1_frme_simv = pd.concat(high_1_simv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "columns_to_analyse_1_simv = ['length_median_high_pass', ]\n",
    "\n",
    "comparison_stats_1_simv = []\n",
    "\n",
    "for column in columns_to_analyse_1_simv:\n",
    "    comparison_stats_1_simv.append(two_paired_test(low_1_frme_simv, high_1_frme_simv, column))\n",
    "\n",
    "comparison_stats_frme_1_simv = DataFrame(comparison_stats_1_simv, index = columns_to_analyse_1_simv)\n",
    "comparison_stats_frme_1_simv;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "columns_to_analyse_2_simv = ['VTemand_kg_mean', 'VTemand_resp_kg_mean', 'VTespon_pat_kg_mean', ]\n",
    "\n",
    "comparison_stats_2_simv = []\n",
    "\n",
    "for column in columns_to_analyse_2_simv:\n",
    "    comparison_stats_2_simv.append(two_paired_test(low_1_frme_simv, high_1_frme_simv, column))\n",
    "\n",
    "comparison_stats_frme_2_simv = DataFrame(comparison_stats_2_simv, index = columns_to_analyse_2_simv)\n",
    "comparison_stats_frme_2_simv;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "comparison_stats_frme_simv = pd.concat([comparison_stats_frme_1_simv, comparison_stats_frme_2_simv])\n",
    "comparison_stats_frme_simv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "writer = pd.ExcelWriter('%s/%s' % (DIR_WRITE, 'vibr_stats_simv.xlsx'))\n",
    "comparison_stats_frme_simv.to_excel(writer, 'stats')\n",
    "writer.save()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### E. Effect of gestation\n",
    "\n",
    "<29 weeks corrected gestational age, all ventilated data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vent_data_below_29 = combined_ventilated_1min_all[combined_ventilated_1min_all['patient'].isin(below_29)].copy()\n",
    "len(vent_data_below_29)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### All received VG ventilation\n",
    "\n",
    "(len(vent_data_below_29), \n",
    " sum(vent_data_below_29['VG_set_kg_mean'].notnull()),\n",
    " sum(vent_data_below_29['VG_set_kg_mean'].isnull()), \n",
    " sum(vent_data_below_29['VG_set_kg_mean'] == 0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vent_data_below_29['VTdiff'] = vent_data_below_29['VTemand_kg_mean'] - vent_data_below_29['VG_set_kg_mean']\n",
    "vent_data_below_29['VTdiff_abs'] = np.abs(vent_data_below_29['VTemand_kg_mean'] - vent_data_below_29['VG_set_kg_mean'])\n",
    "vent_data_below_29['Pdiff_VG'] = vent_data_below_29['PIP_set_mean'] - vent_data_below_29['PIP_mean']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vent_data_below_29.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(vent_data_below_29['accel_rec'].unique()), len(vent_data_below_29['patient'].unique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "low_1_vent_data_below_29= {}; high_1_vent_data_below_29 = {}\n",
    "\n",
    "for patient in vent_data_below_29['patient'].unique():\n",
    "    patient_frame_vent_data_below_29 = vent_data_below_29[vent_data_below_29['patient'] == patient]\n",
    "    \n",
    "    # Only consider minutes when there was at least 0.3 m/sec2 vibration\n",
    "    patient_frame_vent_data_below_29 = \\\n",
    "        patient_frame_vent_data_below_29[patient_frame_vent_data_below_29['length_median_high_pass'] >=0.3]\n",
    "    \n",
    "    low_1_vent_data_below_29[patient] = patient_frame_vent_data_below_29.sort_values('length_median_high_pass', \n",
    "                                                                                ascending = True)[:1]\n",
    "    high_1_vent_data_below_29[patient] = patient_frame_vent_data_below_29.sort_values('length_median_high_pass', \n",
    "                                                                                ascending = True)[-1:]\n",
    "\n",
    "low_1_frme_vent_data_below_29 = pd.concat(low_1_vent_data_below_29)\n",
    "high_1_frme_vent_data_below_29 = pd.concat(high_1_vent_data_below_29)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "columns_to_analyse_1_vent_data_below_29 = ['length_median_high_pass', ]\n",
    "\n",
    "comparison_stats_1_vent_data_below_29 = []\n",
    "\n",
    "for column in columns_to_analyse_1_vent_data_below_29:\n",
    "    comparison_stats_1_vent_data_below_29.append(two_paired_test(low_1_frme_vent_data_below_29, \n",
    "                                                                 high_1_frme_vent_data_below_29, column))\n",
    "\n",
    "comparison_stats_frme_1_vent_data_below_29 = DataFrame(comparison_stats_1_vent_data_below_29, \n",
    "                                                       index = columns_to_analyse_1_vent_data_below_29)\n",
    "comparison_stats_frme_1_vent_data_below_29;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "columns_to_analyse_2_vent_data_below_29 = ['FiO2_set_mean',  'Leak_median', \n",
    "                                           'MV_kg_mean', 'PIP_mean',  'VTemand_kg_mean',]\n",
    "\n",
    "comparison_stats_2_vent_data_below_29 = []\n",
    "\n",
    "for column in columns_to_analyse_2_vent_data_below_29:\n",
    "    comparison_stats_2_vent_data_below_29.append(two_paired_test(low_1_frme_vent_data_below_29, \n",
    "                                                                 high_1_frme_vent_data_below_29, column))\n",
    "\n",
    "comparison_stats_frme_2_vent_data_below_29 = DataFrame(comparison_stats_2_vent_data_below_29, \n",
    "                                                       index = columns_to_analyse_2_vent_data_below_29)\n",
    "comparison_stats_frme_2_vent_data_below_29;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "comparison_stats_frme_vent_data_below_29 = pd.concat([comparison_stats_frme_1_vent_data_below_29, \n",
    "                                                       comparison_stats_frme_2_vent_data_below_29])\n",
    "comparison_stats_frme_vent_data_below_29"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "writer = pd.ExcelWriter('%s/%s' % (DIR_WRITE, 'vibr_stats_below_29weeks_recs.xlsx'))\n",
    "comparison_stats_frme_vent_data_below_29.to_excel(writer, 'stats')\n",
    "writer.save()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 17. Compare periods with low and high sustained acceleration \n",
    "Compare the 1-1 minute with the lowest and highest sustained acceleration in each recording considering those minutes when there was no vibration (1-minute median of Euclidean acceleration was less than 0.3 m/sec2 in the data passed through a high pass filter).\n",
    "This includes periods with the ambulance was not moving and also periods which were dominated by true acceleration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "len(low_1_accel), len(high_1_accel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "low_1_accel_frame = pd.concat(low_1_accel).reset_index(level = 1)\n",
    "low_1_accel_frame;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "high_1_accel_frame = pd.concat(high_1_accel).reset_index(level = 1)\n",
    "high_1_accel_frame;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(low_1_accel_frame), len(high_1_accel_frame)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "(high_1_accel_frame['length_median_low_pass'] - low_1_accel_frame['length_median_low_pass']) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Only include those recordings when the difference is at least 0.25 m.sec2\n",
    "difference = (high_1_accel_frame['length_median_low_pass'] - low_1_accel_frame['length_median_low_pass'])\n",
    "to_included_sustained_accel = difference[difference > 0.25].index\n",
    "to_included_sustained_accel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "low_1_accel_frame = low_1_accel_frame.reindex(to_included_sustained_accel)\n",
    "high_1_accel_frame = high_1_accel_frame.reindex(to_included_sustained_accel)\n",
    "len(low_1_accel_frame), len(high_1_accel_frame)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "columns_to_analyse_1 = ['length_median_low_pass', ]\n",
    "\n",
    "comparison_stats_1_accel = []\n",
    "for column in columns_to_analyse_1:\n",
    "    comparison_stats_1_accel.append(two_paired_test(low_1_accel_frame, high_1_accel_frame, column))\n",
    "\n",
    "comparison_stats_frme_1_accel = DataFrame(comparison_stats_1_accel, index= columns_to_analyse_1)\n",
    "comparison_stats_frme_1_accel;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "columns_to_analyse_2 = ['FiO2_set_median', 'Leak_median', 'MV_kg_mean','PIP_mean',  'VTemand_kg_mean',  ]\n",
    "\n",
    "comparison_stats_2_accel = []\n",
    "for column in columns_to_analyse_2:\n",
    "    comparison_stats_2_accel.append(two_paired_test(low_1_accel_frame, high_1_accel_frame, column))\n",
    "\n",
    "comparison_stats_frme_2_accel = DataFrame(comparison_stats_2_accel, index= columns_to_analyse_2)\n",
    "comparison_stats_frme_2_accel;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "comparison_stats_frme_accel = pd.concat([comparison_stats_frme_1_accel, comparison_stats_frme_2_accel])\n",
    "comparison_stats_frme_accel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "writer = pd.ExcelWriter('%s/%s' % (DIR_WRITE, 'sustained_acccel_stats_all_recs.xlsx'))\n",
    "comparison_stats_frme_accel.to_excel(writer, 'stats')\n",
    "writer.save()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 18. Export the filtered DataFrames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('%s/%s.pickle' % (DATA_DUMP, 'combined_ventilated_1min'), 'wb') as handle:\n",
    "    pickle.dump(combined_ventilated_1min, handle, protocol=pickle.HIGHEST_PROTOCOL) \n",
    "\n",
    "with open('%s/%s.pickle' % (DATA_DUMP, 'low_1_vibr_frame'), 'wb') as handle:\n",
    "    pickle.dump(low_1_vibr_frme, handle, protocol=pickle.HIGHEST_PROTOCOL) \n",
    "\n",
    "with open('%s/%s.pickle' % (DATA_DUMP, 'high_1_vibr_frame'), 'wb') as handle:\n",
    "    pickle.dump(high_1_vibr_frme, handle, protocol=pickle.HIGHEST_PROTOCOL) \n",
    "                            \n",
    "with open('%s/%s.pickle' % (DATA_DUMP, 'low_1_accel_frame'), 'wb') as handle:\n",
    "    pickle.dump(low_1_accel_frame, handle, protocol=pickle.HIGHEST_PROTOCOL) \n",
    "\n",
    "with open('%s/%s.pickle' % (DATA_DUMP, 'high_1_accel_frame'), 'wb') as handle:\n",
    "    pickle.dump(high_1_accel_frame, handle, protocol=pickle.HIGHEST_PROTOCOL)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figures for the paper"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Figure 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Figure 1 ABC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rec = 'default__393.txt', 'AL000628'\n",
    "dpi = 300 \n",
    "filetype = 'jpg'\n",
    "xlim = [-3, 3]\n",
    "bins = 50\n",
    "logyscale = False\n",
    "color = 'blue'\n",
    "\n",
    "fig, axes = plt.subplots(nrows=1, ncols=3, figsize = [12, 3], sharex=False, sharey=False)\n",
    "fig.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.5, hspace=0.5)\n",
    "    \n",
    "ax_ravelled = axes.ravel()\n",
    "\n",
    "for i, par in enumerate(['X', 'Y', 'Z']):\n",
    "    ax = axes[i]\n",
    "    data = accelero_ventilated[rec][par]\n",
    "    # Remove outlier data outside the X axis scale\n",
    "    data = data[(data > xlim[0]) & (data < xlim[1])]\n",
    "    data.hist(ax = ax , bins=bins, color = color)\n",
    "    ax.set_xlabel('m/sec$^2$', size = 10)\n",
    "    ax.set_ylabel('count', size = 10)\n",
    "    ax.set_title(par, size = 10)\n",
    "    ax.set_xlim([*xlim])\n",
    "    ax.set_ylim([0, 60000])\n",
    "            \n",
    "    if logyscale:\n",
    "        ax.set_yscale('log')\n",
    "    \n",
    "fig.savefig('%s/%s.%s' % (DIR_WRITE, 'Figure_1ABC', filetype), \n",
    "    dpi = dpi, format = filetype, bbox_inches='tight', pad_inches=0.1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Figure 1D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dpi = 300\n",
    "filetype = 'jpg'\n",
    "xticklabels = ['X', 'Y', 'Z']\n",
    "\n",
    "medianprops = {'color': 'black', 'linewidth': 2}\n",
    "meanpointprops = {'marker':'D', 'markeredgecolor':'black', 'markerfacecolor':'black'}\n",
    "boxprops = {'color': 'black', 'linestyle': '-'}\n",
    "whiskerprops = {'color': 'black', 'linestyle': '-'}\n",
    "capprops = {'color': 'black', 'linestyle': '-'}\n",
    "flierprops = {'color': 'black', 'marker': '.'}\n",
    "\n",
    "fig, ax = plt.subplots(figsize = (4, 4))\n",
    "\n",
    "plt.boxplot([accelero_stats_median['X_abs'], accelero_stats_median['Y_abs'], accelero_stats_median['Z_abs']], \n",
    "    widths = 0.3, whis = [5, 95], showfliers = False, showmeans = True, meanprops = meanpointprops, \n",
    "        medianprops=medianprops, boxprops=boxprops, whiskerprops=whiskerprops, capprops=capprops, \n",
    "            flierprops = flierprops)\n",
    "\n",
    "ax.hlines([1.1, 1.3], [1, 1], [2, 3], color = 'black')\n",
    "ax.vlines([1, 2, 1, 3,], [1.1, 1.1, 1.3, 1.3], [1.05, 1.05, 1.25, 1.25], color = 'black')\n",
    "ax.text(1.2, 1.15, 'p<0.001'); ax.text(1.55, 1.35, 'p<0.001')\n",
    "\n",
    "ax.set_xticklabels(xticklabels)\n",
    "ax.set_ylabel('Median acceleration of recordings \\n (m/sec$^2$)', size = 12)\n",
    "ax.set_ylim(0, 1.5)\n",
    "ax.tick_params(axis='both', which='major', labelsize=12)\n",
    "plt.grid(True)\n",
    "\n",
    "fig.savefig('%s/%s.%s' % (DIR_WRITE, 'Figure_1D', 'jpg'),\n",
    "    dpi = dpi, format = filetype, bbox_inches='tight', pad_inches=0.1,);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Figure 1E"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filetype = 'jpg'\n",
    "dpi  = 300\n",
    "\n",
    "fig = plt.figure(figsize = [6,6])\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.quiver(0, 0, 0, 3.5, 2.8, 3.5, linewidth = 2, color = 'black', linestyle = 'dashed', arrow_length_ratio = 0.1)\n",
    "ax.quiver(0, 0, 0, 3.5, 0, 0, linewidth = 2, color = 'black', arrow_length_ratio = 0.1)\n",
    "ax.quiver(0, 0, 0, 0, 2.8, 0, linewidth = 2, color = 'black', arrow_length_ratio = 0.1)\n",
    "ax.quiver(0, 0, 0, 0, 0, 3.5, linewidth = 2, color = 'black', arrow_length_ratio = 0.1)\n",
    "\n",
    "ax.quiver(0, 2.8, 0, 0, 0, 3.5, linewidth = 2, linestyle = 'dotted', color = 'black', arrow_length_ratio = 0.1,)\n",
    "ax.quiver(0, 2.8, 3.5, 3.5, 0, 0, linewidth = 2, linestyle = 'dotted', color = 'black', arrow_length_ratio = 0.1,)\n",
    "\n",
    "ax.set_xlim([0, 5]); ax.set_ylim([0,4]); ax.set_zlim([0, 4])\n",
    "ax.set_xlabel('X (m/sec$^2$)'); ax.set_ylabel('Y (m/sec$^2$)'); ax.set_zlabel('Z (m/sec$^2$)')\n",
    "ax.view_init(elev=20., azim = 200)\n",
    "\n",
    "fig.savefig('%s/%s.%s' % (DIR_WRITE, 'Figure_1E', filetype), \n",
    "    dpi = dpi, format = filetype, bbox_inches='tight', pad_inches=0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Figure 1F"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rec = 'default__393.txt', 'AL000628'\n",
    "bins = 50\n",
    "xlim = [0, 4]\n",
    "ylim = [0, 40000]\n",
    "logyscale = False \n",
    "color = 'blue' \n",
    "dpi = 300 \n",
    "filetype = 'jpg'\n",
    "\n",
    "fig, ax = plt.subplots(figsize = [4, 4])\n",
    "    \n",
    "data = accelero_ventilated[rec]['length']\n",
    "# Remove outlier data outside the X axis scale\n",
    "data = data[(data > xlim[0]) & (data < xlim[1])]\n",
    "data.hist(ax = ax, bins=bins, color = color)\n",
    "ax.set_xlabel('m/sec$^2$', size = 12)\n",
    "ax.set_ylabel('count', size = 12)\n",
    "ax.set_title('Euclidean Length', size = 12)\n",
    "ax.set_xlim(*xlim)\n",
    "ax.set_ylim(*ylim)\n",
    "    \n",
    "if logyscale:\n",
    "    ax.set_yscale('log')\n",
    "    \n",
    "fig.savefig('%s/%s.%s' % (DIR_WRITE, 'Figure_1F', filetype), \n",
    "    dpi = dpi, format = filetype, bbox_inches='tight', pad_inches=0.1)\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Figure 1 combined"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dpi = 300\n",
    "filetype = 'pdf'\n",
    "\n",
    "fig = plt.figure(figsize = [9,6], dpi = dpi)\n",
    "fig.subplots_adjust(left=0.1, bottom=0.1, right=0.9, top=0.9, hspace=0.4, wspace=0.5)\n",
    "\n",
    "\n",
    "#Fig 1A,B,C\n",
    "\n",
    "rec = 'default__393.txt', 'AL000628'\n",
    "xlim = [-3, 3]\n",
    "bins = 50\n",
    "logyscale = False\n",
    "color = 'black'\n",
    "\n",
    "for i, par in enumerate(['X', 'Y', 'Z']):\n",
    "    ax = fig.add_subplot(2,3, i+1)\n",
    "    data = accelero_ventilated[rec][par]\n",
    "    # Remove outlier data outside the X axis scale\n",
    "    data = data[(data > xlim[0]) & (data < xlim[1])]\n",
    "    data.hist(ax = ax , bins=bins, color = color)\n",
    "    ax.set_xlabel('m/sec$^2$', size = 10)\n",
    "    ax.set_ylabel('count', size = 10)\n",
    "    ax.set_title(par, size = 10)\n",
    "    ax.tick_params(axis='both', which='major', labelsize=10)\n",
    "    ax.set_xlim([*xlim])\n",
    "    ax.set_ylim([0, 60000])\n",
    "    plt.grid(False)\n",
    "    if logyscale:\n",
    "        ax.set_yscale('log')\n",
    "\n",
    "        \n",
    "# Figure 1D\n",
    "\n",
    "ax = fig.add_subplot(2,3,4,)\n",
    "\n",
    "xticklabels = ['X', 'Y', 'Z']\n",
    "medianprops = {'color': 'black', 'linewidth': 2}\n",
    "meanpointprops = {'marker':'D', 'markeredgecolor':'black', 'markerfacecolor':'black'}\n",
    "boxprops = {'color': 'black', 'linestyle': '-'}\n",
    "whiskerprops = {'color': 'black', 'linestyle': '-'}\n",
    "capprops = {'color': 'black', 'linestyle': '-'}\n",
    "flierprops = {'color': 'black', 'marker': '.'}\n",
    "\n",
    "plt.boxplot([accelero_stats_median['X_abs'], accelero_stats_median['Y_abs'], accelero_stats_median['Z_abs']], \n",
    "    widths = 0.3, whis = [5, 95], showfliers = False, showmeans = True, meanprops = meanpointprops, \n",
    "    medianprops=medianprops, boxprops=boxprops, whiskerprops=whiskerprops, capprops=capprops, \n",
    "    flierprops = flierprops)\n",
    "\n",
    "ax.hlines([1.1, 1.3], [1, 1], [2, 3], color = 'black')\n",
    "ax.vlines([1, 2, 1, 3,], [1.1, 1.1, 1.3, 1.3], [1.05, 1.05, 1.25, 1.25], color = 'black')\n",
    "ax.text(1.2, 1.15, 'p<0.001'); ax.text(1.55, 1.35, 'p<0.001')\n",
    "ax.set_xticklabels(xticklabels)\n",
    "ax.tick_params(axis='both', which='major', labelsize=10)\n",
    "ax.set_ylabel('Median acceleration\\n (m/sec$^2$)', size = 10)\n",
    "ax.set_ylim(0, 1.5)\n",
    "ax.tick_params(axis='both', which='major', labelsize=10)\n",
    "plt.grid(False)\n",
    "\n",
    "\n",
    "# Figure 1E\n",
    "\n",
    "ax = fig.add_subplot(2,3,5, projection='3d')\n",
    "\n",
    "ax.quiver(0, 0, 0, 3.5, 2.8, 3.5, linewidth = 1, color = 'black', linestyle = 'dashed', arrow_length_ratio = 0.15)\n",
    "ax.quiver(0, 0, 0, 3.5, 0, 0, linewidth = 1, color = 'black', arrow_length_ratio = 0.15)\n",
    "ax.quiver(0, 0, 0, 0, 2.8, 0, linewidth = 1, color = 'black', arrow_length_ratio = 0.15)\n",
    "ax.quiver(0, 0, 0, 0, 0, 3.5, linewidth = 1, color = 'black', arrow_length_ratio = 0.15)\n",
    "\n",
    "ax.quiver(0, 2.8, 0, 0, 0, 3.5, linewidth = 1, linestyle = 'dotted', color = 'black', arrow_length_ratio = 0.15,)\n",
    "ax.quiver(0, 2.8, 3.5, 3.5, 0, 0, linewidth = 1, linestyle = 'dotted', color = 'black', arrow_length_ratio = 0.15,)\n",
    "\n",
    "ax.set_xticks(range(0,6),)\n",
    "ax.set_yticks(range(0,5),)\n",
    "ax.set_zticks(range(0,5),)\n",
    "ax.tick_params(axis='both', which='major', labelsize=10)\n",
    "ax.set_xlim([0, 5]); ax.set_ylim([0,4]); ax.set_zlim([0, 4])\n",
    "ax.set_xlabel('X (m/sec$^2$)', fontsize = 10)\n",
    "ax.set_ylabel('Y (m/sec$^2$)', fontsize = 10)\n",
    "ax.set_zlabel('Z (m/sec$^2$)', fontsize = 10)\n",
    "ax.view_init(elev=20., azim = 200)\n",
    "\n",
    "\n",
    "# Figure 1F\n",
    "\n",
    "ax = fig.add_subplot(2,3,6)\n",
    "\n",
    "rec = 'default__393.txt', 'AL000628'\n",
    "bins = 50\n",
    "xlim = [0, 4]\n",
    "ylim = [0, 40000]\n",
    "logyscale = False  \n",
    "   \n",
    "data = accelero_ventilated[rec]['length']\n",
    "# Remove outlier data outside the X axis scale\n",
    "data = data[(data > xlim[0]) & (data < xlim[1])]\n",
    "data.hist(ax = ax, bins=bins, color = color)\n",
    "ax.set_xlabel('m/sec$^2$', size = 10)\n",
    "ax.set_ylabel('count', size = 10)\n",
    "ax.tick_params(axis='both', which='major', labelsize=10)\n",
    "ax.set_title('Euclidean Length', size = 10)\n",
    "ax.set_xlim(*xlim)\n",
    "ax.set_ylim(*ylim)\n",
    "plt.grid(False)\n",
    "if logyscale:\n",
    "    ax.set_yscale('log')\n",
    "    \n",
    "fig.text(0.03, 0.97, 'A', fontsize = 16)\n",
    "fig.text(0.33, 0.97, 'B', fontsize = 16)\n",
    "fig.text(0.62, 0.97, 'C', fontsize = 16)\n",
    "\n",
    "fig.text(0.03, 0.47, 'D', fontsize = 16) \n",
    "fig.text(0.33, 0.47, 'E', fontsize = 16)\n",
    "fig.text(0.62, 0.47, 'F', fontsize = 16)\n",
    "    \n",
    "fig.savefig('%s/%s.%s' % (DIR_WRITE, 'Figure_1', filetype), \n",
    "    dpi = dpi, format = filetype, bbox_inches='tight', pad_inches=0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Figure 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Figure 2A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style('whitegrid', {'axes.grid' : False})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dpi = 300\n",
    "filetype = 'jpg'\n",
    "\n",
    "a = combined_ventilated_1min_all.dropna(subset = ['X_abs_median_low_pass','X_abs_median_high_pass'])\n",
    "#corr = stats.spearmanr(a['X_abs_median_low_pass'], a['X_abs_median_high_pass'])\n",
    "\n",
    "plot = sns.jointplot(x='X_abs_median_low_pass', y='X_abs_median_high_pass', \n",
    "    data = combined_ventilated_1min_all, kind='scatter', color='black', alpha = 0.1, \n",
    "    height=6, ratio=5, space=0.2, dropna=True, xlim=[0, 2], ylim=[0,0.8])\n",
    "\n",
    "plot.ax_joint.hlines(0.15, 0, 2, color = 'black', linestyle = 'dashed')\n",
    "#plot.ax_joint.text(0.1, 0.7, 'r: %.2f\\np<0.0001' % corr[0], fontsize = 14)\n",
    "plot.ax_joint.grid('True')\n",
    "plot.ax_joint.set_xlabel('front-back sustained acceleration (m/sec$^2$)', size = 14)\n",
    "plot.ax_joint.set_ylabel('front-back vibration (m/sec$^2$)', size = 14)\n",
    "plot.ax_joint.tick_params(axis='both', which='major', labelsize=14)\n",
    "\n",
    "plt.savefig('%s/%s.%s' % (DIR_WRITE, 'X_abs_acceleration_X_abs_vibration_jointplot', filetype),\n",
    "    dpi = dpi,  format = filetype, bbox_inches='tight', pad_inches=0.1);\n",
    "\n",
    "plt.savefig('%s/%s.%s' % (DIR_WRITE, 'Figure_2A', filetype),\n",
    "    dpi = dpi,  format = filetype, bbox_inches='tight', pad_inches=0.1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Figure 2B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dpi = 300\n",
    "filetype = 'jpg'\n",
    "\n",
    "a = combined_ventilated_1min_all.dropna(subset = ['Y_abs_median_low_pass','Y_abs_median_high_pass'])\n",
    "#corr = stats.spearmanr(a['Y_abs_median_low_pass'], a['Y_abs_median_high_pass'])\n",
    "\n",
    "plot = sns.jointplot(x='Y_abs_median_low_pass', y='Y_abs_median_high_pass', \n",
    "    data = combined_ventilated_1min_all, kind='scatter', color='black', alpha = 0.1, \n",
    "    height=6, ratio=5, space=0.2, dropna=True, xlim=[0, 1.2], ylim=[0,0.8])\n",
    "\n",
    "plot.ax_joint.hlines(0.07, 0, 2, color = 'black', linestyle = 'dashed')\n",
    "#plot.ax_joint.text(0.1, 0.7, 'r: %.2f\\np<0.0001' % corr[0], fontsize = 14)\n",
    "plot.ax_joint.grid('True')\n",
    "plot.ax_joint.set_xlabel('side-to-side sustained acceleration (m/sec$^2$)', size = 14)\n",
    "plot.ax_joint.set_ylabel('side-to-side vibration (m/sec$^2$)', size = 14)\n",
    "plot.ax_joint.tick_params(axis='both', which='major', labelsize=14)\n",
    "\n",
    "plt.savefig('%s/%s.%s' % (DIR_WRITE, 'Figure_2B', filetype),\n",
    "    dpi = dpi, format = filetype, bbox_inches='tight', pad_inches=0.1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Figure 2C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dpi = 300\n",
    "filetype = 'jpg'\n",
    "\n",
    "a = combined_ventilated_1min_all.dropna(subset = ['Z_abs_median_low_pass','Z_abs_median_high_pass'])\n",
    "#corr = stats.spearmanr(a['Z_abs_median_low_pass'], a['Z_abs_median_high_pass'])\n",
    "\n",
    "plot = sns.jointplot(x='Z_abs_median_low_pass', y='Z_abs_median_high_pass', \n",
    "    data = combined_ventilated_1min_all, kind='scatter', color='black', alpha = 0.1, \n",
    "    height=6, ratio=5, space=0.2, dropna=True, xlim=[0, 1.2], ylim=[0, 1.2])\n",
    "\n",
    "plot.ax_joint.hlines(0.2, 0, 2, color = 'black', linestyle = 'dashed')\n",
    "#plot.ax_joint.text(0.02, 1, 'r: %.2f\\np<0.0001' % corr[0], fontsize = 14)\n",
    "plot.ax_joint.grid('True')\n",
    "plot.ax_joint.set_xlabel('vertical sustained acceleration (m/sec$^2$)', size = 14)\n",
    "plot.ax_joint.set_ylabel('vertical vibration (m/sec$^2$)', size = 14)\n",
    "plot.ax_joint.tick_params(axis='both', which='major', labelsize=14)\n",
    "\n",
    "plt.savefig('%s/%s.%s' % (DIR_WRITE, 'Figure_2C', filetype),\n",
    "    dpi = dpi, format = filetype, bbox_inches='tight', pad_inches=0.1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Figure 2D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dpi = 300\n",
    "filetype = 'jpg'\n",
    "\n",
    "a = combined_ventilated_1min_all.dropna(subset = ['length_median_low_pass','length_median_high_pass'])\n",
    "#corr = stats.spearmanr(a['length_median_low_pass'], a['length_median_high_pass'])\n",
    "\n",
    "plot = sns.jointplot(x='length_median_low_pass', y='length_median_high_pass', \n",
    "    data = combined_ventilated_1min_all, kind='scatter', color='black', alpha = 0.1, \n",
    "    height=6, ratio=5, space=0.2, dropna=True, xlim=[0, 2], ylim=[0, 2])\n",
    "\n",
    "plot.ax_joint.hlines(0.3, 0, 2, color = 'black', linestyle = 'dashed')\n",
    "#plot.ax_joint.text(0.1, 1.7, 'r: %.2f\\np<0.0001' % corr[0], fontsize = 14)\n",
    "plot.ax_joint.grid('True')\n",
    "plot.ax_joint.set_xlabel('Euclidean sustained acceleration (m/sec$^2$)', size = 14)\n",
    "plot.ax_joint.set_ylabel('Euclidean vibration (m/sec$^2$)', size = 14)\n",
    "plot.ax_joint.tick_params(axis='both', which='major', labelsize=14)\n",
    "\n",
    "plt.savefig('%s/%s.%s' % (DIR_WRITE, 'Figure_2D', filetype),\n",
    "    dpi = dpi, format = filetype, bbox_inches='tight', pad_inches=0.1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Supplementary Figure 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "recs = sorted(accelero_ventilated.keys())\n",
    "bins = 50 \n",
    "xlim = [-3, 3] \n",
    "logyscale = False\n",
    "color = 'blue'\n",
    "dpi = 200 \n",
    "filetype = 'pdf'\n",
    "\n",
    "fig, axes = plt.subplots(nrows=len(recs) // 4+1, ncols=12, figsize = [60, 60], sharex=False, sharey=False)\n",
    "fig.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.6, hspace=0.6)\n",
    "    \n",
    "ax_ravelled = axes.ravel()\n",
    "i = 0\n",
    "for rec in recs:\n",
    "    for par in ['X', 'Y', 'Z']:\n",
    "        ax = ax_ravelled[i]\n",
    "        data = accelero_ventilated[rec][par]\n",
    "        # Remove outlier data outside the X axis scale\n",
    "        data = data[(data > xlim[0]) & (data < xlim[1])]\n",
    "        data.hist(ax = ax , bins=bins, color = color)\n",
    "        ax.set_xlabel('m/sec$^2$', size = 10)\n",
    "        ax.set_ylabel('count', size = 10)\n",
    "        ax.set_title('%s - %s' % (rec[1], par), size = 10)\n",
    "        ax.set_xlim([*xlim])\n",
    "        if logyscale:\n",
    "            ax.set_yscale('log') \n",
    "        i+=1\n",
    "    \n",
    "fig.savefig('%s/%s.%s' % (DIR_WRITE, 'Supplementary_Figure_1', filetype), dpi = dpi, \n",
    "    format = filetype, bbox_inches='tight', pad_inches=0.1)\n",
    "\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "recs = sorted(accelero_ventilated.keys())\n",
    "bins = 50 \n",
    "xlim = [-3, 3] \n",
    "logyscale = False\n",
    "color = 'blue'\n",
    "dpi = 200 \n",
    "filetype = 'pdf'\n",
    "\n",
    "fig, axes = plt.subplots(nrows=len(recs), ncols=3, figsize = [8, 3 * len(recs)], sharex=False, sharey=False)\n",
    "fig.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.6, hspace=0.6)\n",
    "    \n",
    "ax_ravelled = axes.ravel()\n",
    "i = 0\n",
    "for rec in recs:\n",
    "    for par in ['X', 'Y', 'Z']:\n",
    "        ax = ax_ravelled[i]\n",
    "        data = accelero_ventilated[rec][par]\n",
    "        # Remove outlier data outside the X axis scale\n",
    "        data = data[(data > xlim[0]) & (data < xlim[1])]\n",
    "        data.hist(ax = ax , bins=bins, color = color)\n",
    "        ax.set_xlabel('m/sec$^2$', size = 10)\n",
    "        ax.set_ylabel('count', size = 10)\n",
    "        ax.set_title('%s - %s' % (rec[1], par), size = 10)\n",
    "        ax.set_xlim([*xlim])\n",
    "        if logyscale:\n",
    "            ax.set_yscale('log') \n",
    "        i+=1\n",
    "    \n",
    "fig.savefig('%s/%s.%s' % (DIR_WRITE, 'Supplementary_Figure_1_mod', filetype), dpi = dpi, \n",
    "    format = filetype, bbox_inches='tight', pad_inches=0.1)\n",
    "\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}