{

 "cells": [

  {

   "cell_type": "markdown",

   "metadata": {},

   "source": [

    "## Part 1: Pulse Rate Algorithm\n",

    "\n",

    "### Contents\n",

    "Fill out this notebook as part of your final project submission.\n",

    "\n",

    "**You will have to complete both the Code and Project Write-up sections.**\n",

    "- The [Code](#Code) is where you will write a **pulse rate algorithm** and already includes the starter code.\n",

    "   - Imports - These are the imports needed for Part 1 of the final project. \n",

    "     - [glob](https://docs.python.org/3/library/glob.html)\n",

    "     - [numpy](https://numpy.org/)\n",

    "     - [scipy](https://www.scipy.org/)\n",

    "- The [Project Write-up](#Project-Write-up) to describe why you wrote the algorithm for the specific case.\n",

    "\n",

    "\n",

    "### Dataset\n",

    "You will be using the **Troika**[1] dataset to build your algorithm. Find the dataset under `datasets/troika/training_data`. The `README` in that folder will tell you how to interpret the data. The starter code contains a function to help load these files.\n",

    "\n",

    "1. Zhilin Zhang, Zhouyue Pi, Benyuan Liu, ‘‘TROIKA: A General Framework for Heart Rate Monitoring Using Wrist-Type Photoplethysmographic Signals During Intensive Physical Exercise,’’IEEE Trans. on Biomedical Engineering, vol. 62, no. 2, pp. 522-531, February 2015. Link\n",

    "\n",

    "-----"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {},

   "source": [

    "### Code"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 41,

   "metadata": {},

   "outputs": [],

   "source": [

    "import glob\n",

    "from tqdm import tqdm\n",

    "import numpy as np\n",

    "import scipy as sp\n",

    "import scipy.io\n",

    "import osPre\n",

    "import scipy.signal\n",

    "import os.path\n",

    "from tqdm import tqdm\n",

    "from sklearn.model_selection import train_test_split\n",

    "from sklearn.linear_model import LinearRegression\n",

    "from sklearn.model_selection import KFold, LeaveOneGroupOut\n",

    "from sklearn.metrics import mean_squared_error\n",

    "def LoadTroikaDataset():\n",

    "    \"\"\"\n",

    "    Retrieve the .mat filenames for the troika dataset.\n",

    "\n",

    "    Review the README in ./datasets/troika/ to understand the organization of the .mat files.\n",

    "\n",

    "    Returns:\n",

    "        data_fls: Names of the .mat files that contain signal data\n",

    "        ref_fls: Names of the .mat files that contain reference data\n",

    "        <data_fls> and <ref_fls> are ordered correspondingly, so that ref_fls[5] is the \n",

    "            reference data for data_fls[5], etc...\n",

    "    \"\"\"\n",

    "    data_dir = \"./datasets/troika/training_data\"\n",

    "    data_fls = sorted(glob.glob(data_dir + \"/DATA_*.mat\"))\n",

    "    ref_fls = sorted(glob.glob(data_dir + \"/REF_*.mat\"))\n",

    "    return data_fls, ref_fls\n",

    "\n",

    "def LoadTroikaDataFile(data_fl):\n",

    "    \"\"\"\n",

    "    Loads and extracts signals from a troika data file.\n",

    "\n",

    "    Usage:\n",

    "        data_fls, ref_fls = LoadTroikaDataset()\n",

    "        ppg, accx, accy, accz = LoadTroikaDataFile(data_fls[0])\n",

    "\n",

    "    Args:\n",

    "        data_fl: (str) filepath to a troika .mat file.\n",

    "\n",

    "    Returns:\n",

    "        numpy arrays for ppg, accx, accy, accz signals.\n",

    "    \"\"\"\n",

    "    data = sp.io.loadmat(data_fl)['sig']\n",

    "    return data[2:]\n",

    "\n",

    "\n",

    "def AggregateErrorMetric(pr_errors, confidence_est):\n",

    "    \"\"\"\n",

    "    Computes an aggregate error metric based on confidence estimates.\n",

    "\n",

    "    Computes the MAE at 90% availability. \n",

    "\n",

    "    Args:\n",

    "        pr_errors: a numpy array of errors between pulse rate estimates and corresponding \n",

    "            reference heart rates.\n",

    "        confidence_est: a numpy array of confidence estimates for each pulse rate\n",

    "            error.\n",

    "\n",

    "    Returns:\n",

    "        the MAE at 90% availability\n",

    "    \"\"\"\n",

    "    # Higher confidence means a better estimate. The best 90% of the estimates\n",

    "    #    are above the 10th percentile confidence.\n",

    "    percentile90_confidence = np.percentile(confidence_est, 10)\n",

    "\n",

    "    # Find the errors of the best pulse rate estimates\n",

    "    best_estimates = pr_errors[confidence_est >= percentile90_confidence]\n",

    "\n",

    "    # Return the mean absolute error\n",

    "    return np.mean(np.abs(best_estimates))\n",

    "\n",

    "def Evaluate():\n",

    "    \"\"\"\n",

    "    Top-level function evaluation function.\n",

    "\n",

    "    Runs the pulse rate algorithm on the Troika dataset and returns an aggregate error metric.\n",

    "\n",

    "    Returns:\n",

    "        Pulse rate error on the Troika dataset. See AggregateErrorMetric.\n",

    "    \"\"\"\n",

    "    # Retrieve dataset files\n",

    "    data_fls, ref_fls = LoadTroikaDataset()\n",

    "    errs, conFs = [], []\n",

    "    for data_fl, ref_fl in zip(data_fls, ref_fls):\n",

    "        # Run the pulse rate algorithm on each trial in the dataset\n",

    "        errors, confidence = RunPulseRateAlgorithm(data_fl, ref_fl)\n",

    "        errs.append(errors)\n",

    "        conFs.append(confidence)\n",

    "        # Compute aggregate error metric\n",

    "    errs = np.hstack(errs)\n",

    "    conFs = np.hstack(conFs)\n",

    "    return AggregateErrorMetric(errs, conFs)\n",

    "\n",

    "def RunPulseRateAlgorithm(data_fl, ref_fl):\n",

    "    \n",

    "    \"\"\"\n",

    "    Run the algorithm \n",

    "    \n",

    "    Args: \n",

    "        data_fls: Names of the .mat files that contain signal data\n",

    "        ref_fls: Names of the .mat files that contain reference data\n",

    "    \n",

    "\n",

    "    Returns:\n",

    "        np.array(error): Array with error for predictions\n",

    "        np.array(confidence): Array with confidence for predictions\n",

    "    \"\"\"\n",

    "    \n",

    "    Fs = 125 # Sample Frequency\n",

    "    window_len = 8 # Window to calculate PR\n",

    "    window_shift = 2  # Difference between windows  \n",

    "       \n",

    "    reg, scores = Regressor()\n",

    "    targets, features, sigs, subs = Data_window8(data_fl, ref_fl)\n",

    "    error, confidence = [], []\n",

    "    for i,feature in enumerate(features):\n",

    "        est = reg.predict(np.reshape(feature, (1, -1)))[0]\n",

    "        \n",

    "        ppg, accx, accy, accz = sigs[i]\n",

    "        \n",

    "        ppg = Filter(ppg)        \n",

    "        accx = Filter(accx)\n",

    "        accy = Filter(accy)\n",

    "        accz = Filter(accz)        \n",

    "        \n",

    "        n = len(ppg) * 3\n",

    "        freq = np.fft.rfftfreq(n, 1/Fs)\n",

    "        fft = np.abs(np.fft.rfft(ppg,n))\n",

    "        fft[freq <= 40/60.0] = 0.0\n",

    "        fft[freq >= 240/60.0] = 0.0\n",

    "    \n",

    "        est_Fs = est / 55.0\n",

    "        Fs_win = 30  / 60.0\n",

    "        Fs_win_e = (freq >= est_Fs - Fs_win) & (freq <= est_Fs +Fs_win)\n",

    "        conf = np.sum(fft[Fs_win_e])/np.sum(fft)\n",

    "        \n",

    "        error.append(np.abs((est-targets[i])))\n",

    "        confidence.append(conf)\n",

    "    return np.array(error), np.array(confidence)\n",

    "\n",

    "def Data_window8(data_fl, ref_fl):\n",

    "    \n",

    "    \"\"\"\n",

    "    Load and prepare the data, based on windows length and shift and filters \n",

    "    \n",

    "    Args: \n",

    "        data_fls: Names of the .mat files that contain signal data\n",

    "        ref_fls: Names of the .mat files that contain reference data\n",

    "    \n",

    "\n",

    "    Returns:\n",

    "        np.array(targets): Array with targets\n",

    "        np.array(features): Array with features\n",

    "        sigs: treated signals\n",

    "        subs: name of file (or subject)\n",

    "    \"\"\"\n",

    "  \n",

    "    Fs=125  # Sampling frequency\n",

    "    window_len = 6 # Window to calculate PR\n",

    "    window_shift = 2 # Difference between windows\n",

    "    \n",

    "    sig = LoadTroikaDataFile(data_fl)\n",

    "    ref = scipy.io.loadmat(ref_fl)[\"BPM0\"]\n",

    "    ref = np.array([x[0] for x in ref])\n",

    "    subject_name = os.path.basename(data_fl).split('.')[0]        \n",

    "    start_indxs, end_indxs = Indexator(sig.shape[1], len(ref), Fs, window_len,window_shift)\n",

    "    targets, features, sigs, subs = [], [], [], []\n",

    "    for i, s in enumerate(start_indxs):\n",

    "        start_i =  start_indxs[i]\n",

    "        end_i = end_indxs[i]\n",

    "\n",

    "        ppg = sig[0, start_i:end_i]            \n",

    "        accx = sig[1, start_i:end_i]\n",

    "        accy = sig[2, start_i:end_i]\n",

    "        accz = sig[3, start_i:end_i]\n",

    "\n",

    "        ppg = Filter(ppg)\n",

    "        accx = Filter(accx)\n",

    "        accy = Filter(accy)\n",

    "        accz = Filter(accz)\n",

    "\n",

    "        feature, ppg, accx, accy, accz = CreateFeature(ppg, accx, accy, accz)\n",

    "\n",

    "        sigs.append([ppg, accx, accy, accz])\n",

    "        targets.append(ref[i])\n",

    "        features.append(feature)\n",

    "        subs.append(subject_name)\n",

    "        \n",

    "    return (np.array(targets), np.array(features), sigs, subs)\n",

    "\n",

    "def Data_window6():\n",

    "    \n",

    "    \"\"\"\n",

    "    Load and prepare the data, based on windows length and shift and filters \n",

    "    \n",

    "    Args: \n",

    "        data_fls: Names of the .mat files that contain signal data\n",

    "        ref_fls: Names of the .mat files that contain reference data\n",

    "    \n",

    "\n",

    "    Returns:\n",

    "        np.array(targets): Array with targets\n",

    "        np.array(features): Array with features\n",

    "        sigs: treated signals\n",

    "        subs: name of file (or subject)\n",

    "    \"\"\"\n",

    "    \n",

    "    Fs=125 # Sampling rate    \n",

    "    window_len = 6 # Window to calculate PR\n",

    "    window_shift = 2 # Difference between windows\n",

    "    \n",

    "    data_fls, ref_fls = LoadTroikaDataset()\n",

    "    pbar = tqdm(list(zip(data_fls, ref_fls)), desc=\"Prepare Data\")\n",

    "    targets, features, sigs, subs = [], [], [], []\n",

    "    for data_fl, ref_fl in pbar:\n",

    "        sig = LoadTroikaDataFile(data_fl)\n",

    "        ref = scipy.io.loadmat(ref_fl)[\"BPM0\"]\n",

    "        ref = np.array([x[0] for x in ref])\n",

    "        subject_name = os.path.basename(data_fl).split('.')[0]        \n",

    "        start_indxs, end_indxs = Indexator(sig.shape[1], len(ref), Fs, window_len,window_shift)\n",

    "        for i, s in enumerate(start_indxs):\n",

    "            start_i =  start_indxs[i]\n",

    "            end_i = end_indxs[i]\n",

    "\n",

    "            ppg = sig[0, start_i:end_i]            \n",

    "            accx = sig[1, start_i:end_i]\n",

    "            accy = sig[2, start_i:end_i]\n",

    "            accz = sig[3, start_i:end_i]\n",

    "\n",

    "\n",

    "            ppg = Filter(ppg)\n",

    "            accx = Filter(accx)\n",

    "            accy = Filter(accy)\n",

    "            accz = Filter(accz)\n",

    "\n",

    "            feature, ppg, accx, accy, accz = CreateFeature(ppg, accx, accy, accz)\n",

    "\n",

    "            sigs.append([ppg, accx, accy, accz])\n",

    "            targets.append(ref[i])\n",

    "            features.append(feature)\n",

    "            subs.append(subject_name)\n",

    "        \n",

    "    return (np.array(targets), np.array(features), sigs, subs)\n",

    "\n",

    "def CreateFeature(ppg, accx, accy, accz):\n",

    "    \"\"\" Create features \n",

    "    \n",

    "        Args: \n",

    "        ppg, accx, accy, accz: signals\n",

    "    \n",

    "\n",

    "    Returns:\n",

    "        np.array([ppg_feature, acc_feature]): features from PPG and ACC signals\n",

    "        ppg, accx, accy, accz: signals\n",

    "        \n",

    "    \"\"\"\n",

    "    ppg = Filter(ppg)\n",

    "    accx = Filter(accx)\n",

    "    accy = Filter(accy)\n",

    "    accz = Filter(accz)\n",

    "    \n",

    "    \n",

    "    Fs = 125\n",

    "    n = len(ppg) * 4\n",

    "    freq = np.fft.rfftfreq(n, 1/Fs)\n",

    "    fft = np.abs(np.fft.rfft(ppg,n))\n",

    "    fft[freq <= 40/60.0] = 0.0\n",

    "    fft[freq >= 240/60.0] = 0.0\n",

    "    \n",

    "    acct = np.sqrt(accx**2 + accy**2 + accz**2) # Total signal of acc\n",

    "    \n",

    "    acc_fft = np.abs(np.fft.rfft(acct, n))\n",

    "    acc_fft[freq <= 40/60.0] = 0.0\n",

    "    acc_fft[freq >= 240/60.0] = 0.0\n",

    "    \n",

    "    ppg_feature = freq[np.argmax(fft)]\n",

    "    acc_feature = freq[np.argmax(acc_fft)]\n",

    "    \n",

    "    return (np.array([ppg_feature, acc_feature]), ppg, accx, accy, accz)\n",

    "\n",

    "def RegressionAlg(features, targets, subs):\n",

    "    \"\"\" \n",

    "    \n",

    "    Create the regression model\n",

    "    \n",

    "        Args: \n",

    "        features: Features obtained from Data_window()\n",

    "        targets: Targets obtained from Data_window()\n",

    "        subs: individuals \n",

    "    \n",

    "\n",

    "    Returns:\n",

    "        regression: result from regression\n",

    "        score: scores from the regression\n",

    "        \n",

    "   \n",

    "    \"\"\"\n",

    "    AdaBoostRegressor\n",

    "    regression = RandomForestRegressor(n_estimators=400,max_depth=16)\n",

    "    scores = []\n",

    "    lf = KFold(n_splits=5)\n",

    "    splits = lf.split(features,targets,subs)\n",

    "    for i, (train_idx, test_idx) in enumerate(splits):\n",

    "        X_train, y_train = features[train_idx], targets[train_idx]\n",

    "        X_test, y_test = features[test_idx], targets[test_idx]\n",

    "        regression.fit(X_train, y_train)\n",

    "        y_pred = regression.Predict(X_test)\n",

    "        score = Error(y_test, y_pred)\n",

    "        scores.append(score)\n",

    "        \n",

    "    return (regression, scores)\n",

    "\n",

    "def Filter(signal):\n",

    "    \n",

    "     \"\"\" \n",

    "    \n",

    "    Bandpass filter between 40 and 240 BPM\n",

    "    \n",

    "        Args: \n",

    "        signal: signal to be filtered\n",

    "    \n",

    "\n",

    "    Returns:\n",

    "        signal:filtered signal\n",

    "        b,a: components from filter\n",

    "        \n",

    "   \n",

    "    \"\"\"\n",

    "    \n",

    "    pass_band=(40/60.0, 240/60.0)\n",

    "    Fs = 125\n",

    "    b, a = scipy.signal.butter(3, pass_band, btype='bandpass', fs=Fs)\n",

    "    return scipy.signal.filtfilt(b, a, signal)\n",

    "\n",

    "def Indexator(sig_len, ref_len, Fs=125, window_len_s=10, window_shift_s=2):\n",

    "    \"\"\"\n",

    "    Find start and end index to iterate over a set of signals\n",

    "    \n",

    "    Args: \n",

    "        sig_len, ref_len: signal and reference lenght\n",

    "        Fs, window_len_s, window_shift: sample frequency, window lenght and window shift\n",

    "    \n",

    "\n",

    "    Returns:\n",

    "        start_indxs: start index\n",

    "        end_indxs: end index\n",

    "    \"\"\"\n",

    "    # Set the length of the biggest signal with regards to the reference signal\n",

    "    if ref_len < sig_len:\n",

    "        n = ref_len\n",

    "    else:\n",

    "        n = sig_len\n",

    "    \n",

    "    # Start Indexes \n",

    "    start_indxs = (np.cumsum(np.ones(n) * Fs * window_shift_s) - Fs * window_shift_s).astype(int)\n",

    "    \n",

    "    # End Indexes (same size as the start indexes array)\n",

    "    end_indxs = start_indxs + window_len_s * Fs\n",

    "    return (start_indxs, end_indxs)\n",

    "\n",

    "def Predict(reg,feature, ppg, accx, accy, accz):\n",

    "    \"\"\"\n",

    "    Create the prediction based on the regression algorithm\n",

    "    \"\"\"\n",

    "    est = reg.predict(np.reshape(feature, (1, -1)))[0]\n",

    "    \n",

    "def Error(y_test, y_pred):\n",

    "    \"\"\"\n",

    "    Calculate error score of the prediction\n",

    "    \"\"\"\n",

    "    return mean_squared_error(y_test, y_pred)\n",

    "\n",

    "\n",

    "def Regressor():\n",

    "    \"\"\"\n",

    "    Apply regression\n",

    "    \"\"\"\n",

    "    fname = \"outfile.npy\"\n",

    "    reg, scores = [], []\n",

    "    if os.path.isfile(fname):\n",

    "        [reg,scores] = np.load(fname,allow_pickle=True)\n",

    "        #\n",

    "    else:\n",

    "        targets, features, sigs, subs = Data_window6()\n",

    "        reg, scores = RegressionAlg(features, targets, subs)\n",

    "        np.save(\"outfile\", [reg,scores])\n",

    "    return reg, scores"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 42,

   "metadata": {},

   "outputs": [

    {

     "data": {

      "text/plain": [

       "9.6607033250413767"

      ]

     },

     "execution_count": 42,

     "metadata": {},

     "output_type": "execute_result"

    }

   ],

   "source": [

    "Evaluate()"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {},

   "source": [

    "-----\n",

    "### Project Write-up\n",

    "\n",

    "Answer the following prompts to demonstrate understanding of the algorithm you wrote for this specific context.\n",

    "\n",

    "> - **Code Description** - Include details so someone unfamiliar with your project will know how to run your code and use your algorithm. \n",

    "> - **Data Description** - Describe the dataset that was used to train and test the algorithm. Include its short-comings and what data would be required to build a more complete dataset.\n",

    "> - **Algorithhm Description** will include the following:\n",

    ">   - how the algorithm works\n",

    ">   - the specific aspects of the physiology that it takes advantage of\n",

    ">   - a describtion of the algorithm outputs\n",

    ">   - caveats on algorithm outputs \n",

    ">   - common failure modes\n",

    "> - **Algorithm Performance** - Detail how performance was computed (eg. using cross-validation or train-test split) and what metrics were optimized for. Include error metrics that would be relevant to users of your algorithm. Caveat your performance numbers by acknowledging how generalizable they may or may not be on different datasets.\n",

    "\n",

    "Your write-up goes here..."

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {},

   "source": [

    "**Code description**\n",

    "\n",

    "The algorithm that predicts heart rate in BPM is based on the estimation of PPG signals and ACC. As requested by Udacity, I used only the second PPG signal (signals were measured on both wrists). The available code returns the mean absolute error and confidence. \n"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {},

   "source": [

    "**Data description**\n",

    "\n",

    "Data Description The Troika data set is used to build the algorithm (https://arxiv.org/pdf/1409.5181.pdf). We had the following signals to work on:\n",

    "- ECG signal \n",

    "- PPG two signals from each wrist \n",

    "- Three channels fro the accelerometer each one corresponding to (x,y and z).\n",

    "- Data is sampled at 125Hz."

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {},

   "source": [

    "**How the algorithm works**\n",

    "\n",

    "The algorithm uses RandomForestRegressor to fit the heart rate training data.\n",

    "\n",

    "- The specific aspects of the physiology that it takes advantage of:\n",

    "\n",

    "Two PPG signals were used (from each wrist). As explained in the course, the capillaries in the wrist fill with blood when the heart's ventricles contract. When the blood returns to the heart there are fewer blood cells in the wrist. The PPG sensor emits a green light which can be absorbed by the red blood cells and the photodetector will see the various levels of blood flow in the reflected light. When the blood returns to the heart fewer blood cells can absorb the green light and the photo detector can see an increase in the reflected light.\n",

    "\n",

    "- A description of the algorithm outputs:\n",

    "\n",

    "The algorithm returns BPM and the confidence rate of this prediction. The higher the confidence rate, the higher we can trust the prediction. \n",

    "\n",

    "- Caveats on algorithm outputs\n",

    "\n",

    "The confidence rate is calculated based on the magnitude of a small area that contains the estimated spectral frequency relative to the sum magnitude of the entire spectrum.\n",

    "\n",

    "- Common failure modes \n",

    "PPG might pick higher frequencies related to motion. To deal with this problem, we take into consideration the accelerometer values.\n",

    "\n"

   ]

  },

  {

   "cell_type": "markdown",

   "metadata": {},

   "source": [

    "**Algorithm Performance**\n",

    "\n",

    "The mean absolute error was calculated and the ground truth reference signal was obtained from the ECG sensors. For cross validation I used the KFold. The error rate was around ~8 BPM on the test set. Since the data is using only limited subjects the algorithm may not be able to generalize well."

   ]

  },

  {

   "cell_type": "code",

   "execution_count": null,

   "metadata": {},

   "outputs": [],

   "source": []

  }

 ],

 "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.6.3"

  }

 },

 "nbformat": 4,

 "nbformat_minor": 4

}