3-feature-engineering.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Feature Engineering\n",
    "Notebook to create all of the required features for all windows, and some preliminary feature exploration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-17T16:17:51.187759Z",
     "start_time": "2021-04-17T16:17:50.408242Z"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import os\n",
    "import gc\n",
    "from matplotlib import pyplot as plt\n",
    "import graphviz\n",
    "import random\n",
    "import pickle\n",
    "import os\n",
    "import gc\n",
    "import ast\n",
    "plt.style.use(\"seaborn\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-14T18:35:32.183079Z",
     "start_time": "2021-04-14T18:34:19.243042Z"
    }
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv(os.path.join('Data', 'data_window_mapping.csv'), \n",
    "                 usecols = ['file_num','window_1','interval','pause_threshold','af_flag'], # 'window_2','window_3'\n",
    "                 dtype = {                 # Specify integer types for memory efficiency\n",
    "                    'file_num':'uint16',\n",
    "                    'interval':'Int16',    # Nullable integer type.\n",
    "                    'af_flag':'uint8',\n",
    "                    'window_1': 'uint32',\n",
    "#                     'window_2': 'Int32',\n",
    "#                     'window_3': 'Int32',\n",
    "                    'pause_threshold': 'bool'               \n",
    "                    }\n",
    "                )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Get count of pauses within each window."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-14T17:35:06.289991Z",
     "start_time": "2021-04-14T17:35:00.637796Z"
    }
   },
   "outputs": [],
   "source": [
    "# 1. Filter to records exceeding the pause threshold, 2. Count the number of these per 30-beat segment in each file.\n",
    "pause_counts = df[df['pause_threshold']].groupby(by=['file_num', 'window_1']).size() \n",
    "pause_counts.hist(bins=31,range=(0,30),color='tab:grey',alpha=0.8, histtype='bar', ec='grey')   # Plot a histogram. Min is 1 and max is 30.\n",
    "plt.xlabel(\"Number of pause counts per 30-beat window\", fontsize=16)\n",
    "plt.ylabel(\"# of windows\", fontsize=16)\n",
    "plt.tick_params(axis='both', which='major', labelsize=12)\n",
    "plt.tick_params(axis='both', which='minor', labelsize=10)\n",
    "plt.tight_layout()\n",
    "\n",
    "# Count the total number of 30-beat segments across files to compare the % with Pauses.\n",
    "index_counts = df.groupby(by=['file_num', 'window_1']).size()\n",
    "print(\"{0:.1f}% of 30-beat segments contain at least one pause\".format(100*len(pause_counts)/len(index_counts)))\n",
    "print(\"{0:.1f}% of 30-beat segments contain more than three pauses\".format(100*np.sum(pause_counts>3)/len(index_counts)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-14T17:38:47.203026Z",
     "start_time": "2021-04-14T17:38:46.869809Z"
    }
   },
   "outputs": [],
   "source": [
    "pause_counts.to_csv(os.path.join('Data', 'pause_counts_window1.csv'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define the variables and helper functions required to calculate the summary statistics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-14T18:08:59.144005Z",
     "start_time": "2021-04-14T18:08:59.128283Z"
    }
   },
   "outputs": [],
   "source": [
    "# Define the bins outside the function rather than repeatedly calling np.linspace.\n",
    "hist_bins = np.linspace(300.0, 3000.0, 55)    # Use a small bin width to start, which can be aggregated later if necessary.\n",
    "\n",
    "# Bins for the successive differences\n",
    "# Wider bins between +-1000 and +-100 than between (-100, 100).  (Widths of 100 vs. 5)\n",
    "diff_hist_bins = np.concatenate((np.linspace(-1000,-200, 9), \n",
    "                       np.linspace(-100,100,41), \n",
    "                       np.linspace(200,1000,9)\n",
    "                      ))\n",
    "def hist_func(x):\n",
    "    # Input is a series through the pd.groupby().agg() call.\n",
    "    return np.histogram(x, \n",
    "                        bins=hist_bins,                     # Use the predefined bins above.\n",
    "                        weights=(np.ones(len(x)) / len(x))  # Return as proportion falling into the bin.\n",
    "                       )[0].tolist()\n",
    "\n",
    "def diff_hist_func(x):\n",
    "    # Input is a series through the pd.groupby().agg() call.\n",
    "    xd = x.diff()   # Get difference between successive rows. (Rows for the same ECG recording will be supplied separately to others)\n",
    "    xd = xd[~pd.isnull(xd)]\n",
    "    return np.histogram(xd, \n",
    "                        bins=diff_hist_bins,               # Use the predefined bins above.\n",
    "                        weights=(np.ones(len(xd))/len(xd)) # Return as proportion falling into the bin.\n",
    "                       )[0].tolist()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate the summary statistics for each set of windows.\n",
    "Takes approx. 40 mins to run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-08T19:33:28.134025Z",
     "start_time": "2021-04-08T18:50:45.804127Z"
    }
   },
   "outputs": [],
   "source": [
    "# Calculate summary statistics for each window in 'window_1' - done for each window group (1 to 3).\n",
    "df_window1 = (\n",
    "    # Calculate stats only where interval is non-null. Group by the file number & window number.\n",
    "    df[~pd.isnull(df['interval'])].groupby(['file_num','window_1'])[['interval', 'af_flag']]\n",
    "    .agg(\n",
    "        mean_interval=(\"interval\", np.nanmean),\n",
    "        median_interval=(\"interval\", np.nanmedian),\n",
    "        max_interval=(\"interval\", np.nanmax),\n",
    "        min_interval=(\"interval\", np.nanmin),\n",
    "        var_interval=(\"interval\", np.nanvar),\n",
    "        # Compute all percentiles simultaneously for efficiency (sort once)\n",
    "        percentiles_interval=(\"interval\", lambda x: np.percentile(x, q=[5,25,75,95]).tolist()),    \n",
    "        hists=(\"interval\", hist_func),\n",
    "        diff_hists=(\"interval\", diff_hist_func),\n",
    "        af_flag=(\"af_flag\", lambda x: np.sum(x)/len(x))\n",
    "    )\n",
    ")\n",
    "\n",
    "df_window1.to_csv(os.path.join('Data', 'df_window1.csv'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-14T18:48:33.777113Z",
     "start_time": "2021-04-14T18:48:33.731657Z"
    }
   },
   "outputs": [],
   "source": [
    "gc.collect() # Free up memory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### RMSSD and Turning Points\n",
    "In the following cell, the features for Root Mean Squared Successive Difference (RMSSD) and Turning Points are calculated. The turning points among groups of 3 successive beats are identified by first getting the difference between successive beats, and then multiplying it by the same value on the next row. If the product is positive (which means successive increases or successive decreases) then there is no turning point, if the product is negative then there is a turning point. These turning points are placed at the point they occur in the data by  using periods=-1 for the shift call."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-14T18:47:48.198057Z",
     "start_time": "2021-04-14T18:35:32.217119Z"
    }
   },
   "outputs": [],
   "source": [
    "df['interval_diff'] = df.groupby('file_num')['interval'].diff(periods=1)\n",
    "df['interval_diff_shift'] = df.groupby('file_num')['interval_diff'].shift(periods=-1)\n",
    "# If the product is negative then it is a turning point (but avoid replacing the NANs with True or False).\n",
    "df['turning_point'] = np.where(\n",
    "                            np.isnan(df['interval_diff']*df['interval_diff_shift']), \n",
    "                            np.nan, \n",
    "                            df['interval_diff']*df['interval_diff_shift'] < 0\n",
    "                        )\n",
    "df_window1_supp = (\n",
    "    df.astype({'interval_diff': 'float32', 'turning_point': 'float32'})         # Convert from int to float for the numeric operations below.\n",
    "    .groupby(['file_num','window_1'])[['interval_diff', 'turning_point']]   # Use only the relevant columns after grouping.\n",
    "    .agg(\n",
    "        # Root -> Mean -> Square -> of interval differences (converted to float)\n",
    "        RMSSD=(\"interval_diff\", lambda x: np.sqrt(np.nanmean(np.square(pd.to_numeric(x, errors='coerce'))))),\n",
    "        # Get the count of turning points divided by the total possible number for the window.\n",
    "        turning_point_ratio=(\"turning_point\", lambda x: np.nansum(x)/np.count_nonzero(~np.isnan(x)))\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-14T18:47:54.623728Z",
     "start_time": "2021-04-14T18:47:48.312195Z"
    }
   },
   "outputs": [],
   "source": [
    "df_window1_supp.to_csv(os.path.join('Data', 'df_window1_supp.csv'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-14T18:47:54.686860Z",
     "start_time": "2021-04-14T18:47:54.672865Z"
    }
   },
   "outputs": [],
   "source": [
    "df_window1_supp.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del df, df_window1, df_window1_supp\n",
    "gc.collect()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Final Cleaning Up and Writing to Numpy Array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Main feature data. \n",
    "df1 = pd.read_csv(os.path.join('Data', 'df_window1.csv')) # Some columns stored as string of lists.\n",
    "df1_supp = pd.read_csv(os.path.join('Data', 'df_window1_supp.csv'))   # Supplementary feature data - RMSSD & TPR.\n",
    "df1_pause_counts = pd.read_csv(os.path.join('Data', 'pause_counts_window1.csv'))   # The number of pause periods per window.\n",
    "\n",
    "print(\"Shapes of dataframes are df1: {},   df1_supp: {},   pause_counts: {}\".format(df1.shape, df1_supp.shape, df1_pause_counts.shape))\n",
    "\n",
    "#### Remove windows with more than 3 pause periods\n",
    "\n",
    "pause_gt_3 = df1_pause_counts[df1_pause_counts['0'] > 3]   # Get the windows with more than 3 pause periods.\n",
    "print(\"Pauses > 3: {} | Total with pauses: {} | % of total > 3: {}\".format( \\\n",
    "    len(pause_gt_3), len(df1_pause_counts), round(float(len(pause_gt_3))/len(df1_pause_counts), 4)))\n",
    "print(\"\\n\", \"-\"*40, \"\\n\")\n",
    "\n",
    "# Get all the rows of data except windows with more than 3 pause periods.\n",
    "df1 = df1[~ df1[['file_num','window_1']].apply(tuple,1).isin(pause_gt_3[['file_num','window_1']].apply(tuple,1))]\n",
    "df1_supp = df1_supp[~ df1_supp[['file_num','window_1']].apply(tuple,1).isin(pause_gt_3[['file_num','window_1']].apply(tuple,1))]\n",
    "\n",
    "print(\"DF shapes with pauses removed\\ndf1: {}  |  df1_supp: {}\".format(df1.shape, df1_supp.shape))\n",
    "\n",
    "## Percentiles and histogram have been stored as strings of lists - convert back to lists.\n",
    "# df_pct = pd.DataFrame([ast.literal_eval(l) for l in df1['percentiles_interval']], index=df1.index)\n",
    "# df_hists = pd.DataFrame([ast.literal_eval(l) for l in df1['hists']], index=df1.index)\n",
    "# df_diff_hists = pd.DataFrame([ast.literal_eval(l) for l in df1['diff_hists']], index=df1.index)\n",
    "\n",
    "## Create column names for the hist & percentile data.\n",
    "hist_bins = np.linspace(300.0, 3000.0, 55) # Pre-defined histogram bins.\n",
    "diff_hist_bins = np.concatenate((np.linspace(-1000,-200, 9), \n",
    "                       np.linspace(-100,100,41), \n",
    "                       np.linspace(200,1000,9)\n",
    "                      ))\n",
    "\n",
    "hist_cols = ['hist_' + str(int(n)) for n in hist_bins[:-1]]\n",
    "diff_hist_cols = ['diff_hist_' + str(int(n)) for n in diff_hist_bins[:-1]]\n",
    "pct_cols = ['percentile_' + n for n in ['05','25','75','95']]\n",
    "\n",
    "df1[pct_cols] = df_pct\n",
    "df1[hist_cols] = df_hists\n",
    "df1[diff_hist_cols] = df_diff_hists\n",
    "\n",
    "df1.reset_index(drop=True, inplace=True)\n",
    "df1_supp.reset_index(drop=True, inplace=True)\n",
    "\n",
    "f1[['RMSSD','turning_point_ratio']] = df1_supp[['RMSSD','turning_point_ratio']]\n",
    "\n",
    "# df1.to_csv(os.path.join('Data','df1_total.csv'))   # Save to file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Final features and save to numpy arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1 = pd.read_csv(os.path.join('Data','df1_total.csv')) \n",
    "\n",
    "df1['hist_lt_300'] = np.round(1-np.sum(df1[hist_cols], axis=1),6)\n",
    "\n",
    "hist_cols.append('hist_lt_300')\n",
    "\n",
    "#### A very small number of nans in the variance and turning point ratio columns. These occurred where the RR interval was the same number for the entire window. Hence variance and turning point ratio calculations failed. Since these are clearly low variance windows, NAs should be replaced with zero.\n",
    "\n",
    "df1[['var_interval','turning_point_ratio']] = df1[['var_interval','turning_point_ratio']].fillna(0)  # Fill nans in the selected columns.\n",
    "\n",
    "# Drop redundant columns and clean up memory.\n",
    "df1.drop(['Unnamed: 0','percentiles_interval','hists','diff_hists'], axis=1, inplace=True)\n",
    "# gc.collect()\n",
    "\n",
    "#### Calculate (Normalised) Shannon Entropy. Treat each histogram bin as one class.\n",
    "\n",
    "df1['entropy'] = entropy(df1[hist_cols], base=2.0, axis=1)/np.log2(len(hist_cols))   # Calculate Entropy.\n",
    "\n",
    "df1.to_csv(os.path.join('Data','df1_final.csv'))   # Save to file.\n",
    "\n",
    "# Save as numpy file format to save time on writing and loading.\n",
    "y = np.array(df1[['af_flag']] > 0.6)   # If more than 60% of beats are AF then window is labelled AF.\n",
    "ids=np.array(df1[['file_num','window_1']])\n",
    "X = np.array(df1.drop(['af_flag','file_num','window_1'], axis=1))\n",
    "\n",
    "np.save(os.path.join('Data','Final', 'y.npy'), y)\n",
    "np.save(os.path.join('Data','Final', 'ids.npy'), ids)\n",
    "np.save(os.path.join('Data','Final', 'X.npy'), X)\n",
    "\n",
    "col_names = df1.drop(['af_flag','file_num','window_1'], axis=1).columns.tolist()\n",
    "with open (os.path.join('Data','Final', 'col_names.txt'), 'w') as f:\n",
    "    for item in col_names:\n",
    "        f.write(\"%s\\n\" % item)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plots for feature analysis and data exploration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### Plot the distributions of entropy, RMSSD and turning point ratio for the AF and Non-AF classes.\n",
    "\n",
    "def compare_hist(col_name, ylabel):\n",
    "    if (col_name == 'entropy') | (col_name == 'turning_point_ratio'):\n",
    "        max_col = np.max(df1[col_name])\n",
    "        bins = np.linspace(0,0.1+np.floor(max_col/0.1)*0.1, int(np.floor(max_col/0.1))+2)\n",
    "    elif col_name == 'RMSSD':\n",
    "        bins=np.linspace(0,800,17)\n",
    "    else:\n",
    "        bins=None\n",
    "\n",
    "    fig, ax = plt.subplots(2,1, figsize=(8,8))\n",
    "\n",
    "    ax[0].hist(df1[df1['af_flag'] >= 0.6][col_name],bins=bins, label='AF',\\\n",
    "               color='tab:orange',alpha=0.8, histtype='bar', ec='grey')\n",
    "    ax[0].set_ylabel('# of windows',fontsize=16)\n",
    "    ax[0].tick_params(axis='both', which='major', labelsize=12)\n",
    "    ax[0].tick_params(axis='both', which='minor', labelsize=10)\n",
    "\n",
    "    ax[1].hist(df1[df1['af_flag'] < 0.6][col_name],bins=bins, label='Non-AF',\\\n",
    "               color = 'tab:blue',alpha=0.8, histtype='bar', ec='grey')\n",
    "    ax[1].set_ylabel('# of windows',fontsize=16)\n",
    "    ax[1].set_xlabel(ylabel,fontsize=16)\n",
    "    ax[1].tick_params(axis='both', which='major', labelsize=12)\n",
    "    ax[1].tick_params(axis='both', which='minor', labelsize=10)\n",
    "    \n",
    "    if col_name == 'RMSSD':\n",
    "        ax[1].yaxis.offsetText.set_fontsize(14)\n",
    "\n",
    "    fig.legend(fontsize=16)\n",
    "    fig.tight_layout()\n",
    "    \n",
    "# compare_hist('turning_point_ratio', 'Turning Point Ratio')\n",
    "# compare_hist('entropy', 'Normalised Entropy')\n",
    "# compare_hist('RMSSD', 'Root Mean Square Successive Differences (RMSSD) [ms]')\n",
    "\n",
    "### Examine the distribution of AF Flag proportion of beats per window.\n",
    "\n",
    "print(\"Number of windows with AF_Flag between 0-100%:\\t{}\".format(np.sum((df1['af_flag']>0) & (df1['af_flag']< 1))))\n",
    "print(\"Number of windows with AF_Flag between 60-100%:\\t{}\".format(sum((df1['af_flag']>=0.6) & (df1['af_flag']< 1))))\n",
    "print(\"Number of windows with AF_Flag = 100%:\\t\\t{}\".format(sum(df1['af_flag'] == 1)))\n",
    "print(\"Number of windows with AF_Flag = 0%:\\t\\t{}\".format(sum(df1['af_flag'] == 0)))\n",
    "\n",
    "df1[(df1['af_flag']>0) & (df1['af_flag']< 1)]['af_flag'].hist(\n",
    "    bins=np.linspace(0,1,11),color='tab:grey',alpha=0.8, histtype='bar', ec='grey')\n",
    "\n",
    "plt.xlabel(\"Proportion of beats labelled as AF\", fontsize=16)\n",
    "plt.ylabel(\"# of windows\", fontsize=16)\n",
    "plt.tick_params(axis='both', which='major', labelsize=12)\n",
    "plt.tick_params(axis='both', which='minor', labelsize=10)\n",
    "fig.tight_layout()\n",
    "\n",
    "print(\"{} files have at least one episode of AF. This represents {:.1f}% of all files\".format( \\\n",
    "    df1[df1['af_flag']>=0.6]['file_num'].nunique(), 100*df1[df1['af_flag']>=0.6]['file_num'].nunique()/df1['file_num'].nunique()))\n",
    "\n",
    "af_counts = df1[df1['af_flag']>=0.6].groupby('file_num').size()   # Get count of AF episodes per file number.\n",
    "\n",
    "af_counts.hist(range=(0,5500), bins=11, color = 'tab:grey',alpha=0.9, histtype='bar', ec='grey')\n",
    "\n",
    "plt.xlabel('# of AF Windows per Recording File', fontsize=16)\n",
    "plt.ylabel('# Recording Files', fontsize=16)\n",
    "plt.tick_params(axis='both', which='major', labelsize=12)\n",
    "plt.tick_params(axis='both', which='minor', labelsize=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "------\n",
    "## Previous Attempts \n",
    "More verbose and less efficient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-15T17:21:20.252976Z",
     "start_time": "2021-02-15T17:19:27.560872Z"
    }
   },
   "outputs": [],
   "source": [
    "def stats_func(x):\n",
    "    if(sum(pd.notnull(x['interval'])) > 0):\n",
    "        return pd.Series(simple_summary_stats(x[~pd.isnull(x['interval'])]))\n",
    "\n",
    "tmp = df[df['file_num'] < 10]\n",
    "tmp_stats = (\n",
    "    tmp.groupby(['file_num','window_1'])[['interval', 'af_flag', 'pause_threshold']]\n",
    "    .apply(stats_func)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-15T16:49:33.896371Z",
     "start_time": "2021-02-15T16:49:33.874566Z"
    }
   },
   "outputs": [],
   "source": [
    "g1 = df.groupby(['file_num','window_1'])[['interval', 'af_flag', 'pause_threshold']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-15T17:10:57.255943Z",
     "start_time": "2021-02-15T16:49:43.477351Z"
    }
   },
   "outputs": [],
   "source": [
    "stats = g1.agg(\n",
    "    mean_interval=(\"interval\", np.nanmean),\n",
    "    median_interval=(\"interval\", np.nanmedian),\n",
    "    max_interval=(\"interval\", np.nanmax),\n",
    "    min_interval=(\"interval\", np.nanmin),\n",
    "    var_interval=(\"interval\", np.nanvar),\n",
    "    q05_interval=(\"interval\", lambda x: x.quantile(q=0.05)),\n",
    "    q25_interval=(\"interval\", lambda x: x.quantile(q=0.25)),\n",
    "    q75_interval=(\"interval\", lambda x: x.quantile(q=0.75)),\n",
    "    q95_interval=(\"interval\", lambda x: x.quantile(q=0.95)),\n",
    "    af_flag=(\"af_flag\", lambda x: np.sum(x)/len(x)),\n",
    "    pause_threshold_num=(\"pause_threshold\", np.sum)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-15T17:16:05.648609Z",
     "start_time": "2021-02-15T17:16:05.617075Z"
    }
   },
   "outputs": [],
   "source": [
    "def simple_summary_stats(df):\n",
    "    \"\"\" \n",
    "    Calculates:\n",
    "        1. The simple summary statistics for the window \n",
    "        2. The proportion of beats annotated as AF\n",
    "        3. The number of beats with an RRI exceeding the pause threshold.\n",
    "    Input: \n",
    "        df        The dataframe for one group (window)\n",
    "        window    Which window to group by. Needs to be called three separate times to capture the three windows.\n",
    "    Returns:\n",
    "        Grouped df with the aggregations (features) as columns.\n",
    "    \"\"\"\n",
    "    return {\n",
    "        'mean_interval': np.mean(df['interval']),\n",
    "        'median_interval': np.median(df['interval']),\n",
    "        'max_interval': np.max(df['interval']),\n",
    "        'min_interval': np.min(df['interval']),\n",
    "        'var_interval': np.var(df['interval']),\n",
    "        'q05_interval': np.percentile(df['interval'], q=5),\n",
    "        'q25_interval': np.percentile(df['interval'], q=25),\n",
    "        'q75_interval': np.percentile(df['interval'], q=75),\n",
    "        'q95_interval': np.percentile(df['interval'], q=95),\n",
    "        'af_flag': np.sum(df['af_flag']/len(df)),\n",
    "        'pause_threshold_num': np.sum(df['pause_threshold'])\n",
    "    }\n",
    "\n",
    "def hists(x, bins):\n",
    "    \"\"\"\n",
    "    Return the proportion of values in certain bins, given a series x.\n",
    "    \"\"\"\n",
    "    # Null values should not be captured by the histogram.\n",
    "    x_clean = x[~pd.isnull(x)]\n",
    "    return pd.Series(\n",
    "        np.histogram(\n",
    "            x_clean,\n",
    "            bins=bins,\n",
    "            weights=(np.ones(len(x_clean)) / len(x_clean))  # Normalise by dividing by the number of data points.\n",
    "        )[0]                                           # Function returns the proportions & the bin limits, take only the proportion per bin.\n",
    "    )\n",
    "\n",
    "def diff_hists(x, bins):\n",
    "    \"\"\"\n",
    "    Find the difference between successive beats (rows) and then calculate the histogram.\n",
    "    \"\"\"\n",
    "    return(hists(x.diff(), bins))\n",
    "\n",
    "def calculate_features(x):\n",
    "    \"\"\"\n",
    "    Function to be used in call to .apply. Defines the bins and calculates the histograms for the RRI and delta RRI.\n",
    "    \"\"\"\n",
    "    if x.index[0] % 10000 == 0: print(x.index[0]) \n",
    "    if(sum(pd.notnull(x['interval'])) > 0):\n",
    "        d = {}\n",
    "        # Find the histogram breakdown of the RR intervals.\n",
    "        hist_vals = hists(pd.Series(x['interval']), bins=hist_bins)\n",
    "        # Add each histogram as a feature to the dictionary to be returned.\n",
    "        for i, val in enumerate(hist_vals):\n",
    "            d['hist_' + str(bins[i])] = val\n",
    "\n",
    "        # Find the histogram breakdown of the differenced RR intervals.\n",
    "        diff_hist_vals = diff_hists(pd.Series(x['interval']), bins=diff_hist_bins)\n",
    "        for i, val in enumerate(diff_hist_vals):\n",
    "            d['diff_hist_' + str(bins[i])] = val     \n",
    "\n",
    "#         d.update(simple_summary_stats(x[~pd.isnull(x['interval'])]))\n",
    "        return pd.Series(d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-15T16:16:00.805358Z",
     "start_time": "2021-02-15T15:26:26.309191Z"
    }
   },
   "outputs": [],
   "source": [
    "features_df_1 = (\n",
    "df\n",
    "    .groupby(['file_num','window_1'])[['interval', 'af_flag', 'pause_threshold']]\n",
    "    .apply(calculate_features)\n",
    ")\n",
    "## Interrupted after 49min so likely to take a lot longer than that."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "oldHeight": 264.667,
   "position": {
    "height": "380px",
    "left": "456px",
    "right": "20px",
    "top": "166px",
    "width": "705px"
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "varInspector_section_display": "none",
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}