<jupyter_start><jupyter_text>Part 1: Pulse Rate Algorithm ContentsFill out this notebook as part of your final project submission.**You will have to complete both the Code and Project Write-up sections.**- The [Code](Code) is where you will write a **pulse rate algorithm** and already includes the starter code.   - Imports - These are the imports needed for Part 1 of the final project.      - [glob](https://docs.python.org/3/library/glob.html)     - [numpy](https://numpy.org/)     - [scipy](https://www.scipy.org/)- The [Project Write-up](Project-Write-up) to describe why you wrote the algorithm for the specific case. DatasetYou 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.1. <NAME>, <NAME>, <NAME>, ‘‘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-----  Code<jupyter_code>import matplotlib.pyplot as plt
import scipy.signal as sg
import pandas as pd
import copy
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split
import pickle
from sklearn.metrics import mean_absolute_error, mean_squared_error
from sklearn.preprocessing import normalize
import os

import glob

import numpy as np
import scipy as sp
import scipy.io


def LoadTroikaDataset():
    """
    Retrieve the .mat filenames for the troika dataset.

    Review the README in ./datasets/troika/ to understand the organization
    of the .mat files.

    Returns:
        data_fls: Names of the .mat files that contain signal data
        ref_fls: Names of the .mat files that contain reference data
        <data_fls> and <ref_fls> are ordered correspondingly,
        so that ref_fls[5] is the reference data for data_fls[5], etc...
    """
    data_dir = "./datasets/troika/training_data"
    data_fls = sorted(glob.glob(data_dir + "/DATA_*.mat"))
    ref_fls = sorted(glob.glob(data_dir + "/REF_*.mat"))
    return data_fls, ref_fls


def LoadTroikaDataFile(data_fl):
    """
    Loads and extracts signals from a troika data file.

    Usage:
        data_fls, ref_fls = LoadTroikaDataset()
        ppg, accx, accy, accz = LoadTroikaDataFile(data_fls[0])

    Args:
        data_fl: (str) filepath to a troika .mat file.

    Returns:
        numpy arrays for ppg, accx, accy, accz signals.
    """
    data = sp.io.loadmat(data_fl)["sig"]
    return data[2:]


# in your prediction file
def load_model(model_path):
    """
    Loads a saved sklearn regression model or train new model
    Args:
    Returns:
        return the saved model
    """
    path_n = model_path
    if os.path.isfile(path_n):
        fil = open(path_n, "rb")
        prediction_model = pickle.load(fil)
        print("model found")
    else:
        prediction_model = model_train()
        print("new model is created")
    return prediction_model


def Load_labels(data_fl, ts):
    """
    Load the reference PBM for a file and perform interpolation to convert it
    to the same timestamp of the measurments
    Args:
        data_fl: (str) filepath to a troika .mat file.
        ts: time stamp of the sensors measurements for the same file
    Returns:
        numpy arrays for ref BPM with the same length of sensor measurements
    """
    # as mentioed in the dataset read_me the ref value is obtained using a
    # window of 8s and 6s overlap
    # so we can assume that the sample freq of the ref data is 2s
    fs_labels = 2
    data = sp.io.loadmat(data_fl)["BPM0"]
    ts_labels = np.arange(0, len(data) / fs_labels, 1 / fs_labels)
    labels_interp = np.interp(ts, ts_labels, data[:, 0])
    return labels_interp


def butter_bandpass(data, rang, fs):
    """
    bandpass the data withn certain frequency range for all sensors
    Args:
        data: a list of sensors measurements
        range: tuple contain the min, max frequency
        fs: sample frequence for the signal
    Returns:
        return a list arrays for all the filtered signals in the time domain
    """
    results = []
    zeros, poles = sg.butter(5, rang, btype="bandpass", fs=fs)
    # filter each sensor measurement
    for key in data:
        results.append(sg.filtfilt(zeros, poles, key))
    return results


def fft_sensors(data, fs):
    """
    fast fourier transform for the each sensor data
    Args:
        data: a list of sensors measurements
        fs: sample frequence for the signal
    Returns:
        return a list of arrays for all the signals in the frequency domaoin
        data is on this format data - > sensor1_freq, sensor1_powers
                                    - > sensor2_freq, sensor2_powers
                                    - > sensorn_freq, sensorn_powers
    """
    ff_data = []
    # FFT each sensor measurement
    for key in data:
        ff_freq = np.fft.rfftfreq(len(key), 1.0 / fs)
        ff_power = np.fft.rfft(key)
        ff_data.append(np.array([ff_freq, ff_power]))
    return ff_data


def fft_inv(fftdata):
    """
    inverse fast fourier transform for each sensor data
    Args:
        data: a list of arrays for each sensor frequencies and powers
        fs: sample frequence for the signal
    Returns: list of arrays for each sensor signal in the time domain
    """
    data = {}
    # IFFT for each sensor data
    for key in fftdata:
        # Each sensor has freqs, power arrays so we chose the power here
        data[key] = np.fft.irfft(fftdata[key][1])
    return data


# Deprecated -> was used for visualization in the begining only
def vizualize_data(data_dict_fs, wind_title, xlim, vtype="freq"):
    len_dict = len(data_dict_fs)
    plt.clf
    fig, ax = plt.subplots(len_dict, 1, figsize=(10, 6 * len_dict))
    fig.canvas.set_window_title(wind_title)
    ax = ax.flatten()
    for key, axes in zip(data_dict_fs, ax):
        if type == "freq":
            power = np.abs(key[0])
            freqs = np.abs(key[1])
            axes.plot(freqs, power, label=key)
            axes.set_xlabel(key)
            axes.set_xlim(xlim)
            axes.legend()
        else:
            axes.plot(key)
            # axes.set_xlabel(key)
            axes.set_xlim(xlim)
            axes.legend()

        # axes.set_xlabel('Time (sec)')


def load_all_dataset(data_fls, ref_fls, fs, window_size=6, window_overlab=4):

    """
    Load data, labels from all files, windowing on each 6 seconds data with
    4 sec overlab each sample data is a window of 6 seconds and shifted 2 s
    from the previoues one
    Args:
        data_fls: all the files for the dataset sensor measurements
        ref_fls: all the files for the dataset sensor labels
        fs: sample frequency
    Returns: list of array for all samples in the dataset, each sample data
     has all the sensors measurements in this sample in the time domain
    """
    X = []
    Y = []
    window_len = 6 * fs
    window_shift_len = (window_size - window_overlab) * fs

    dataset = {"ppg": [], "accx": [], "accy": [], "accz": [], "labels": []}
    # Read All Files
    for data_file, ref_file in zip(data_fls, ref_fls):
        ppg, accx, accy, accz = LoadTroikaDataFile(data_file)
        dataset["ppg"].extend(ppg)
        dataset["accx"].extend(accx)
        dataset["accy"].extend(accy)
        dataset["accz"].extend(accz)
        ts = np.arange(0, len(ppg) / fs, 1 / fs)
        labels = Load_labels(ref_file, ts)
        labels /= 60.0
        dataset["labels"].extend(labels)
    # dataframe have data and labels for all the data and labels
    # all files are appended
    df = pd.DataFrame(dataset)
    # Convert dataframe to samples each one is shifted 2s
    number_samples = int(len(df) / window_shift_len)
    for window_ind in range(number_samples):
        sample_start = int(window_ind * window_shift_len)
        sample_end = int(sample_start + window_len)
        ppg = df["ppg"][sample_start:sample_end].values
        accx = df["accx"][sample_start:sample_end].values
        accy = df["accy"][sample_start:sample_end].values
        accz = df["accz"][sample_start:sample_end].values
        labels = df["labels"][sample_start:sample_end].values
        X.append(np.array([ppg, accx, accy, accz]))
        Y.append(np.array([labels]))
    return X, Y


def clean_datset(samples, labels, fs):
    """
    Clean each seample data by performing bandpass filter on range of 0.7-4
    Hz and add FFT for each sample in the dataset
    Args:
        samples: list of array for each sample. each sample data has several
         sensor measurements
        labels: list of arrays for each sample labels
    Returns: tuble of dataset time and frequency data for all sensors, labels
     for each sample
    """
    clean_x = []
    clean_y = []
    for sample, label in zip(samples, labels):
        bandpass_filtered = butter_bandpass(sample, (0.7, 4), fs)
        fft_sample = fft_sensors(bandpass_filtered, fs)
        clean_x.append((bandpass_filtered, fft_sample))
        clean_y.append([np.mean(label)])
    return (clean_x, clean_y)


def AggregateErrorMetric(pr_errors, confidence_est):
    """
    Computes an aggregate error metric based on confidence estimates.

    Computes the MAE at 90% availability.

    Args:
        pr_errors: a numpy array of errors between pulse rate estimates and
         corresponding reference heart rates.
        confidence_est: a numpy array of confidence estimates for each pulse
         rate error.

    Returns:
        the MAE at 90% availability
    """
    # Higher confidence means a better estimate. The best 90% of the estimates
    #    are above the 10th percentile confidence.
    percentile90_confidence = np.percentile(confidence_est, 10)

    # Find the errors of the best pulse rate estimates
    best_estimates = pr_errors[confidence_est >= percentile90_confidence]

    # Return the mean absolute error
    return np.mean(np.abs(best_estimates))


def remove_powerful_acc_freqs_from_ppg(
    peaks_freq_ppg, peaks_freqs_ppg_vals, peaks_freq_acc
):
    """
    in any sample, filter the ppg peak frequencies that exist in the
     acceleration peak frequencies .
    If an acceleration frequency exists within +-0.25Hz range for any of
     the ppg peak freqs, this freq is removed.
    Args:
        peaks_freq_ppg: a list of the powerful freqs in the ppg freq signal
        peaks_freqs_ppg_vals: power of each ppg peak frequency
        peaks_freq_acc:a list of the powerful freqs in the accel freq signal
    Returns:return a tuple of the most powerful ppg freq and its power after
     removing near acc freqs
    """
    # if there is only one peak in this sample, take this sample
    if len(peaks_freq_ppg) > 1:

        # sort ppg freqs according to their power
        sorted_ppg_freqs = peaks_freq_ppg[np.argsort(peaks_freqs_ppg_vals)]
        sorted_ppg_powers = peaks_freqs_ppg_vals[np.argsort(
                                                peaks_freqs_ppg_vals)]
        ppg_freqs_acc_removed = []
        ppg_freqs_acc_removed_val = []
        # for each ppg freq search for powerful freqs near the range
        for freq, power in zip(sorted_ppg_freqs, sorted_ppg_powers):
            flag = 0
            for acc_freq in peaks_freq_acc:
                if (freq + 0.15 <= acc_freq) & (acc_freq <= freq + 0.15):
                    flag = 1
                    break
            if flag == 0:
                ppg_freqs_acc_removed.append(freq)
                ppg_freqs_acc_removed_val.append(power)

        ppg_1st_power_freq_acc_removed = ppg_freqs_acc_removed[-1]
        ppg_1st_power_freq_acc_removed_val = ppg_freqs_acc_removed_val[-1]
    else:
        ppg_1st_power_freq_acc_removed = peaks_freq_ppg[-1]
        ppg_1st_power_freq_acc_removed_val = peaks_freqs_ppg_vals[-1]
    return (ppg_1st_power_freq_acc_removed, ppg_1st_power_freq_acc_removed_val)


def featurize_samples(samples, fs):
    """
    samples->(time_data,freq_data)->(ppg,accx,accy,accz)
                                    ->fppg,faccx,faccy,faccz->freqs,power
    """
    """
    time and freq featurization for each sample sensors data
    Args:
        samples: a list of the dataset samples
        fs: sample frequency
    Returns: return a list of arrays for each sample features
    """
    features = []
    for sample in samples:
        # extract each sample to each sensor time & freq data
        ppg_t = sample[0][0]
        freqs = np.abs(sample[1][0][0])
        ppg_f = np.abs(sample[1][0][1])
        accx_f = np.abs(sample[1][1][1])
        accy_f = np.abs(sample[1][2][1])
        accz_f = np.abs(sample[1][3][1])
        #  average freq power for all accel axes
        acc_all_f = (accx_f + accy_f + accz_f) / 3.0

        # Time features
        peaks, _ = sg.find_peaks(
            ppg_t, height=1.5, distance=30
        )
        # minimum 0.25s distance(4Hz)
        # periods features
        peaks_dif = np.diff(peaks) / 60.0
        min_peak = np.min(peaks_dif)
        max_peak = np.max(peaks_dif)
        mean_peak = np.mean(peaks_dif)
        std_peak = np.std(peaks_dif)

        # peaks for ppg freqs & their power
        # minimum 1/3Hz distance
        peaks_freq_ppg, _ = sg.find_peaks(ppg_f, distance=3)
        peaks_freqs_ppg_vals = ppg_f[peaks_freq_ppg]
        peaks_freq_ppg = freqs[peaks_freq_ppg]
        # take the peaks in the allowed freq only
        peaks_freqs_ppg_vals = peaks_freqs_ppg_vals[
            (peaks_freq_ppg < 4.0) & (peaks_freq_ppg > 0.67)
        ]
        peaks_freq_ppg = peaks_freq_ppg[
            (peaks_freq_ppg < 4.0) & (peaks_freq_ppg > 0.67)
        ]
        # peaks for accel freqs & their power
        peaks_freq_acc, _ = sg.find_peaks(
            acc_all_f, distance=3
        )  # minimum 1/3Hz distance
        peaks_freqs_acc_vals = acc_all_f[peaks_freq_acc]
        peaks_freq_acc = freqs[peaks_freq_acc]
        peaks_freqs_acc_vals = peaks_freqs_acc_vals[
            (peaks_freq_acc < 4.0) & (peaks_freq_acc > 0.67)
        ]
        peaks_freq_acc = peaks_freq_acc[
            (peaks_freq_acc < 4.0) & (peaks_freq_acc > 0.67)
        ]
        # freq features for ppg
        ppg_freq_mean = np.mean(peaks_freq_ppg)
        ppg_freq_std = np.std(peaks_freq_ppg)
        ppg_1st_power_freq = peaks_freq_ppg[np.argsort(
                                    peaks_freqs_ppg_vals)][-1]
        ppg_1st_power_freq_val = peaks_freqs_ppg_vals[np.argsort(
                                    peaks_freqs_ppg_vals)][-1]
        (ppg_1st_freq_acc_fil, ppg_1st_freq_acc_fil_val,) = \
            remove_powerful_acc_freqs_from_ppg(
                        peaks_freq_ppg, peaks_freqs_ppg_vals, peaks_freq_acc
                                    )

        # freq features for acc
        acc_all_freq_mean = np.mean(peaks_freq_acc)
        acc_all_freq_std = np.std(peaks_freq_acc)
        acc_all_max_strong_freq = peaks_freq_acc[np.argsort(
                                            peaks_freqs_acc_vals)][-1]
        acc_all_max_strong_freq_power = peaks_freqs_acc_vals[
            np.argsort(peaks_freqs_acc_vals)
        ][-1]

        time_features = [min_peak, max_peak, mean_peak, std_peak]
        freq_features = [
            ppg_freq_mean,
            ppg_freq_std,
            ppg_1st_power_freq,
            ppg_1st_power_freq_val,
            ppg_1st_freq_acc_fil,
            ppg_1st_freq_acc_fil_val,
            acc_all_freq_mean,
            acc_all_freq_std,
            acc_all_max_strong_freq,
            acc_all_max_strong_freq_power,
        ]
        total_features = time_features + freq_features
        features.append(np.array(total_features))

    return np.array(features)


def model_train(estimators=650, depth=14, file_path="model_final"):
    """
    train a random forest regressor on the dataset and save it in the
     provided path
    Args:
        estimators: number of trees in the model
        depth: single value for the depth of each tree of the model
        file_path: a path string to which model will be saved
    Returns: return the trained regression model
    """
    # Reading ref and sensors data, create timestamp for both
    fs = 125
    data_fls, ref_fls = LoadTroikaDataset()
    samples, labels = load_all_dataset(data_fls, ref_fls, fs)
    clean_x, clean_y = clean_datset(samples, labels, fs)
    """
    samples->(time_data,freq_data)->(ppg,accx,accy,accz)
                                    ->fppg,faccx,faccy,faccz->freqs,power
    """
    dataset_feats = featurize_samples(clean_x, fs)
    train_x, test_x, train_y, test_y = train_test_split(
        dataset_feats, clean_y, random_state=42, test_size=0.2
    )
    clf = RandomForestRegressor(
        n_estimators=estimators, max_depth=depth, random_state=42
    )
    clf.fit(train_x, np.ravel(train_y))
    y_pred = clf.predict(test_x)
    mae_val = mean_absolute_error(
        y_pred * 60.0, np.array(test_y) * 60.0
    )  # *60 toconvert from Hz t BPM
    mse_val = mean_squared_error(
        y_pred * 60, np.array(test_y) * 60
    )  # *60 to convert from Hz t BPM
    # res={"label":np.ravel(test_y)*60,"prediction":y_pred*60}
    print("Model Results : \n")
    print("mean absolute error", mae_val)
    print("mean square error", mse_val)
    with open(file_path, "wb") as f:
        pickle.dump(clf, f)
        print("model saved in the following dir: %s" % file_path)
    return clf


def Evaluate():
    """
    Top-level function evaluation function.
    Runs the pulse rate algorithm on the Troika dataset and returns an
     aggregate error metric.

    Returns:
        Pulse rate error on the Troika dataset. See AggregateErrorMetric.
    """
    # Retrieve dataset files
    data_fls, ref_fls = LoadTroikaDataset()
    errs, confs = [], []
    for data_fl, ref_fl in zip(data_fls, ref_fls):
        # Run the pulse rate algorithm on each trial in the dataset
        errors, confidence = RunPulseRateAlgorithm(data_fl, ref_fl)
        # errors=list(errors)
        # confidence=list(confidence)
        # print(errors.shape,confidence.shape)
        errs.append(errors)
        confs.append(confidence)
        # Compute aggregate error metric
    errs = np.hstack(errs)
    confs = np.hstack(confs)
    # print(len(errs),len(confs))
    return AggregateErrorMetric(errs, confs)


def RunPulseRateAlgorithm(data_fl, ref_fl):
    """
    load regresssion model, perform prediction over a file  samples
    return mean error and confidence for each sample
    Args:
        data_fl: file name for sensor measurements
        ref_fl: file name that has true PBM value
        file_path: a path string to which model will be saved
    Returns: tuple of two arrays for mean error confidence
    """
    model_path = "model_final"
    fs = 125
    samples, labels = load_all_dataset([data_fl], [ref_fl], fs)
    clean_x, clean_y = clean_datset(samples, labels, fs)
    file_samples_x = featurize_samples(clean_x, fs)
    file_samples_y = np.array(clean_y)
    # ppg, accx, accy, accz = LoadTroikaDataFile(data_fl)
    reg_model = load_model(model_path)
    samples_pred = reg_model.predict(file_samples_x)

    mae_val = np.abs(
        (samples_pred.reshape(-1, 1)) * 60.0 - np.array(file_samples_y) * 60.0
    )  # *60 toconvert from Hz t BPM
    # Compute pulse rate estimates and estimation confidence.
    confidence = []
    mae_val = np.ravel(mae_val)
    for ind in range(len(clean_y)):
        pred_freq = samples_pred[ind]
        fft_freqs = np.abs(clean_x[ind][1][0][0])
        fft_f_power = np.abs(clean_x[ind][1][0][1])
        freq_wind = \
            (fft_freqs > pred_freq - 0.5) & (fft_freqs < pred_freq + 0.5)
        conf = np.sum(fft_f_power[freq_wind]) / np.sum(fft_f_power)
        confidence.append(conf)
    # print(mae_val.shape,np.array(confidence).shape)
    # Return per-estimate mean absolute error and confidence as
    # a 2-tuple of numpy arrays.
    return mae_val, np.array(confidence)


def get_predictions(data_fl, ref_fl, model_path="model_final"):
    """
    load regresssion model, perform prediction over a file  samples
    return prediction and confidence for each sample
    Args:
        data_fl: file name for sensor measurements
        model_path: a path string to which model will be saved
    Returns: a data for predictions and confidence of each 2s measurement
    """
    fs = 125
    samples, labels = load_all_dataset([data_fl], [ref_fl], fs)
    clean_x, clean_y = clean_datset(samples, labels, fs)
    file_samples_x = featurize_samples(clean_x, fs)
    file_samples_y = np.array(clean_y)
    # ppg, accx, accy, accz = LoadTroikaDataFile(data_fl)
    reg_model = load_model(model_path)
    samples_pred = reg_model.predict(file_samples_x)
    confidence = []
    for ind in range(len(clean_y)):
        pred_freq = samples_pred[ind]
        fft_freqs = np.abs(clean_x[ind][1][0][0])
        fft_f_power = np.abs(clean_x[ind][1][0][1])
        freq_wind = \
            (fft_freqs > pred_freq - 0.5) & (fft_freqs < pred_freq + 0.5)
        conf = np.sum(fft_f_power[freq_wind]) / np.sum(fft_f_power)
        confidence.append(conf)
    # print(len(file_samples_y),len(samples_pred),len(confidence))
    results = {
        "true PPM": np.ravel(file_samples_y) * 60.0,
        "predicted_PPM": samples_pred * 60.0,
        "confidence": confidence,
    }
    return pd.DataFrame(results)

model_train(estimators=650, depth=14, file_path="model_final")
Evaluate()<jupyter_output>model found
model found
model found
model found
model found
model found
model found
model found
model found
model found
model found
model found<jupyter_text>----- Project Write-upAnswer the following prompts to demonstrate understanding of the algorithm you wrote for this specific context.> - **Code Description** - Include details so someone unfamiliar with your project will know how to run your code and use your algorithm. > - **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.> - **Algorithhm Description** will include the following:>   - how the algorithm works>   - the specific aspects of the physiology that it takes advantage of>   - a describtion of the algorithm outputs>   - caveats on algorithm outputs >   - common failure modes> - **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.Your write-up goes here...  Code Description Running the code 1. download all requred libraries are avaiable and all dependancies with Python 3.7 2. to run test if everything is working, you can use ```RunPulseRateAlgorithm``` function. To use this function you should provide 1 file that have PPG , 3-axis Accelerometer sensors measurements and 1 more file that have the true BPM values. All files should be the same format to Troika Dataset. 3. to run the training, testing process at the same at once, you can use ```Evaluate``` function which will do everything and return the mean average error of the algorithm. Troika dataset should be provided in the same directory of the notebook in a folder called 'datasets' 4. to run the project for getting BPM values for certain file you can run the following code<jupyter_code>data_fs,ref_fs=LoadTroikaDataset()
df=get_predictions(data_fs[0],ref_fs[0])
df.sample(10)<jupyter_output>model found
{'fraction_non_alphanumeric': 0.08375005252762953, 'fraction_numerical': 0.011009791150144976, 'mean_word_length': 2.7726696258719086, 'pattern_counts': {'":': 11, '<': 20, '<?xml version=': 0, '>': 41, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '390057', 'license': 'MIT', 'n_tokens_mistral': 7759, 'n_tokens_neox': 7490, 'n_words': 2259, 'path': 'part_1/.ipynb_checkpoints/pulse_rate_starter-checkpoint.ipynb', 'repo_name': 'ahmedokasha000/Pulse-Rate-Estimation-For-Wearables', 'size': '23812'}