Training.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: xgboost in c:\\users\\szzcc\\.conda\\envs\\mhealth23\\lib\\site-packages (1.6.2)\n",
      "Requirement already satisfied: numpy in c:\\users\\szzcc\\.conda\\envs\\mhealth23\\lib\\site-packages (from xgboost) (1.21.6)\n",
      "Requirement already satisfied: scipy in c:\\users\\szzcc\\.conda\\envs\\mhealth23\\lib\\site-packages (from xgboost) (1.7.3)\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import math\n",
    "from Lilygo.Recording import Recording\n",
    "from Lilygo.Dataset import Dataset\n",
    "from os import listdir\n",
    "from os.path import isfile, join\n",
    "from math import sqrt\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "# feature extraction\n",
    "from scipy import signal\n",
    "from scipy.stats import entropy, kurtosis, skew\n",
    "from scipy.signal import welch\n",
    "from scipy.interpolate import interp1d\n",
    "# to create a test and train split\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.utils.class_weight import compute_class_weight\n",
    "# evaluation metrics\n",
    "from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score, confusion_matrix, ConfusionMatrixDisplay\n",
    "\n",
    "# Gradient boosting\n",
    "import xgboost as xgb\n",
    "from xgboost import XGBClassifier\n",
    "# save and load models\n",
    "import joblib\n",
    "\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.optim as optim\n",
    "from torch.utils.data import TensorDataset, DataLoader"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate features for training and testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This function aims to find the component caused by gravity from data, which means the signal around 0 Hz\n",
    "def get_gravity(data):\n",
    "    filtered_data = np.zeros_like(data)\n",
    "    # Parameters in IIR filter\n",
    "    alpla = [1, -1.979133761292768, 0.979521463540373]\n",
    "    beta = [0.000086384997973502, 0.00012769995947004, 0.000086384997973502]\n",
    "    # Formula of IIR filter\n",
    "    for i in range(2, len(data)):\n",
    "        filtered_data[i] = alpla[0] * (data[i] * beta[0] + data[i-1] * beta[1] + data[i-2] * beta[2] - filtered_data[i-1] * alpla[1] - filtered_data[i-2] * alpla[2])\n",
    "    return filtered_data\n",
    "\n",
    "# This function aims to realize a low-pass filter with cutoff frequency = 1 Hz. Because according to massive amounts of data, the general \n",
    "# minimum frequency of human walking is about 1 Hz\n",
    "def get_highpass(data):\n",
    "    filtered_data = np.zeros_like(data)  # filtered_data\n",
    "    alpla = [1, -1.905384612118461, 0.910092542787947]\n",
    "    beta = [0.953986986993339, -1.907503180919730, 0.953986986993339]\n",
    "    for i in range(2, len(data)):\n",
    "        filtered_data[i] = alpla[0] * (data[i] * beta[0] + data[i-1] * beta[1] + data[i-2] * beta[2] - filtered_data[i-1] * alpla[1] - filtered_data[i-2] * alpla[2])\n",
    "    return filtered_data\n",
    "\n",
    "# This funciton aims to realize a high-pass filter with cutoff frequency = 5 Hz. Because according to massive amounts of data, the general \n",
    "# maximum frequency of human walking is about 5 Hz\n",
    "def get_lowpass(data):\n",
    "    filtered_data = np.zeros_like(data)  # filtered_data\n",
    "    alpla = [1, -1.80898117793047, 0.827224480562408]\n",
    "    beta = [0.096665967120306, -0.172688631608676, 0.095465967120306]\n",
    "    for i in range(2, len(data)):\n",
    "        filtered_data[i] = alpla[0] * (data[i] * beta[0] + data[i-1] * beta[1] + data[i-2] * beta[2] - filtered_data[i-1] * alpla[1] - filtered_data[i-2] * alpla[2])\n",
    "    return filtered_data\n",
    "\n",
    "def pre_process(data):\n",
    "    # Find the component caused by gravity from data and remove it from the singanl\n",
    "    data_gravity = get_gravity(data)\n",
    "    data_user = data - data_gravity\n",
    "    # Get user's acceleration along the gravity direction by dot product\n",
    "    data_acc = data_user * data_gravity\n",
    "    # Add low pass and high pass filter to reduce noise in signal (possible human walking rate:1 - 5Hz)\n",
    "    data_filtered = get_highpass(data_acc)\n",
    "    data_filtered = get_lowpass(data_filtered)\n",
    "    return data_filtered\n",
    "\n",
    "def get_segment_overlap(data,sampling_rate, std_win):\n",
    "    # Calculate raw magnitude of accelerometer signal\n",
    "    # amagn_acc = [np.sqrt(a**2+trace.data['ay'].values[i]**2+trace.data['az'].values[i]**2)for i, a in enumerate(trace.data['ay'].values)]\n",
    "    # Pre-process data\n",
    "    data_seg = pre_process(data)\n",
    "    # Calculate window size\n",
    "    window_size = round(std_win*sampling_rate)\n",
    "    segment_trace = [data_seg[s:s+window_size] for s in range(0, len(data_seg)-window_size, round(window_size/2))]\n",
    "    return segment_trace\n",
    "\n",
    "def get_segment_nonoverlap(data,sampling_rate, std_win):\n",
    "    # Calculate window size\n",
    "    window_size = round(std_win*sampling_rate)\n",
    "    segment_trace = [data[s:s+window_size] for s in range(0, len(data)-window_size, window_size)]\n",
    "    return segment_trace"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### features for acticity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_temporal_features(data):\n",
    "    \"\"\"\n",
    "    returns time domain features (20) from windowed raw signal:\n",
    "      - mean, standard deviation, kurtosis, skewness;\n",
    "      - RMS, zero-crossing\n",
    "      - The following percentiles: [0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 100] (Q)\n",
    "      - Range: max(x) - min(x)\n",
    "    data: windowed signal\n",
    "    \"\"\"\n",
    "    m = np.mean(data)\n",
    "    sd = np.std(data)\n",
    "    kurt = kurtosis(data)\n",
    "    sk = skew(data)\n",
    "    rms = np.sqrt(np.mean(data**2))\n",
    "    zc = np.sum(np.diff(data>=m))\n",
    "    q = np.array([0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 1])\n",
    "    Q = np.quantile(data, q)\n",
    "    range = Q[-1]-Q[0]             \n",
    "    return m, sd, kurt, sk, rms, zc, range, Q[0], Q[1], Q[2], Q[3], Q[4], Q[5], Q[6], Q[7], Q[8], Q[9], Q[10], Q[11], Q[12]\n",
    "\n",
    "\n",
    "def extract_frequential_features(data, sr):\n",
    "    \"\"\"\n",
    "    exctracting frequncy domain features (5) from windowed Fourier transform:\n",
    "    - energy, entropy, centroid, bandwidth, max_freq\n",
    "    data: windowed signal\n",
    "    sf: sampling rate\n",
    "    \"\"\"\n",
    "    window_size = len(data)\n",
    "    data -= np.mean(data)\n",
    "    ft = np.fft.fft(data)/window_size\n",
    "    sr = int(sr)\n",
    "    # get window length\n",
    "    \n",
    "    #discarding mirror part\n",
    "    ft = ft[:window_size//2]\n",
    "    #frequencies of the transofm\n",
    "    freqs = np.fft.fftfreq(window_size, 1/sr)[1:window_size//2]\n",
    "    #the spectral density is the squared of the absolute\n",
    "    Spec = np.abs(ft)**2\n",
    "    #Energy\n",
    "    E = np.sum(Spec)/(window_size//2)\n",
    "    #density\n",
    "    P = Spec[1:]/np.sum(Spec[1:])\n",
    "    #entropy\n",
    "    H = -np.sum(P*np.log2(P))/np.log2((window_size//2))\n",
    "    #centriod \n",
    "    C = np.sum(P*freqs)\n",
    "    #Absolute distance  of frequencies from from Centroid\n",
    "    distC = np.abs((C-freqs))\n",
    "    #bandwidth is the weighted mean of the distance\n",
    "    BW = np.sum(distC*P)\n",
    "    #maximum frequency \n",
    "    max_fr = freqs[np.argmax(Spec[1:])]\n",
    "    return E, H, C, BW, max_fr\n",
    "\n",
    "# 25 features\n",
    "def get_features_activity_win(data):\n",
    "    # AC, DC = get_AC_DC(data, 200)\n",
    "    m, sd, kurt, sk, rms, zc, range, q_0, q_1, q_2, q_3, q_4, q_5, q_6, q_7, q_8, q_9, q_10, q_11, q_12 = extract_temporal_features(data)\n",
    "    E, H, C, BW, max_fr = extract_frequential_features(data, 200)\n",
    "    all_features = [m, sd, kurt, sk, rms, zc, range, q_0, q_1, q_2, q_3, q_4, q_5, q_6, q_7, q_8, q_9, q_10, q_11, q_12, E, H, C, BW, max_fr]\n",
    "    return all_features\n",
    "\n",
    "def get_features_activity_trace(trace):\n",
    "    # get segment data\n",
    "    amagn = [np.sqrt(a**2+trace.data['ay'].values[i]**2+trace.data['az'].values[i]**2)for i, a in enumerate(trace.data['ax'].values)]\n",
    "    std_win = 3 # length of window in seconds\n",
    "    sampling_rate = 200\n",
    "    segment_trace = get_segment_overlap(amagn, sampling_rate, std_win)\n",
    "    # add feature extraction\n",
    "    num_features = 25\n",
    "    featured_trace_activity = np.zeros((np.shape(segment_trace)[0], num_features))\n",
    "    for i in range(np.shape(segment_trace)[0]):\n",
    "        featured_trace_activity[i,] = get_features_activity_win(segment_trace[i])\n",
    "    return featured_trace_activity"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### features for location"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_feature_location(ax, ay, az, gx, gy, gz, sampling_rate = 200):\n",
    "    acc_data = np.array([ax, ay, az]).T\n",
    "    gyro_data = np.array([gx, gy, gz]).T\n",
    "    # Compute magnitudes\n",
    "    acc_magnitude = np.sqrt(np.sum(acc_data**2, axis=1))\n",
    "    acc_mag_mean = abs(np.mean(acc_magnitude))\n",
    "    acc_mag_std = np.std(acc_magnitude)\n",
    "\n",
    "\n",
    "    # ----ACCELERATOR TIME DOMAIN----\n",
    "    ax_mean = abs(np.mean(acc_data[:,0]))\n",
    "    ay_mean = abs(np.mean(acc_data[:,1]))\n",
    "    az_mean = abs(np.mean(acc_data[:,2]))\n",
    "    a_mean_list = [ax_mean, ay_mean, az_mean]\n",
    "    a_mean_list.sort() # Sorting list of numbers in ascending\n",
    "    Am = a_mean_list[2] # Feature Am: the maximum mean among all dimensions (represents motion range for location)\n",
    "    Bm = a_mean_list[2]/a_mean_list[1] # Feature Bm and Cm: ratio of the maximum mean in different axes (represents DoF in movement for location)\n",
    "    Cm = a_mean_list[2]/a_mean_list[0] \n",
    "\n",
    "    ax_range = acc_data[np.argmax(acc_data[:,0]), 0] - acc_data[np.argmin(acc_data[:,0]), 0]\n",
    "    ay_range = acc_data[np.argmax(acc_data[:,1]), 1] - acc_data[np.argmin(acc_data[:,1]), 1]\n",
    "    az_range = acc_data[np.argmax(acc_data[:,2]), 2] - acc_data[np.argmin(acc_data[:,2]), 2]\n",
    "    a_range_list = [ax_range, ay_range, az_range]\n",
    "    a_range_list.sort() # Sorting list of numbers in ascending\n",
    "    A = a_range_list[2] # Feature A: the maximum range among all dimensions (represents motion range for location)\n",
    "    B = a_range_list[2]/a_range_list[1] # Feature B and C: ratio of the maximum ranges in different axes (represents DoF in movement for location)\n",
    "    C = a_range_list[2]/a_range_list[0] \n",
    "    \n",
    "    \n",
    "    # ----GYROSCOPE TIME DOMAIN----\n",
    "    gx_mean = abs(np.mean(gyro_data[:,0]))\n",
    "    gy_mean = abs(np.mean(gyro_data[:,1]))\n",
    "    gz_mean = abs(np.mean(gyro_data[:,2]))\n",
    "    g_mean_list = [gx_mean, gy_mean, gz_mean]\n",
    "    g_mean_list.sort() # Sorting list of numbers in ascending\n",
    "\n",
    "    gx_range = gyro_data[np.argmax(gyro_data[:,0]), 0] - gyro_data[np.argmin(gyro_data[:,0]), 0]\n",
    "    gy_range = gyro_data[np.argmax(gyro_data[:,1]), 1] - gyro_data[np.argmin(gyro_data[:,1]), 1]\n",
    "    gz_range = gyro_data[np.argmax(gyro_data[:,2]), 2] - gyro_data[np.argmin(gyro_data[:,2]), 2]\n",
    "    g_range_list = [gx_range, gy_range, gz_range]\n",
    "    g_range_list.sort() # Sorting list of numbers in ascending\n",
    "    G = g_range_list[2] \n",
    "    H = g_range_list[2]/g_range_list[1] \n",
    "    I = g_range_list[2]/g_range_list[0]\n",
    "\n",
    "\n",
    "    # ----ACCELERATOR FREQUENCY DOMAIN----\n",
    "    freq, Pxx = welch(acc_magnitude, fs=sampling_rate) # use of the fast Fourier transform for the estimation of power spectra\n",
    "    freq_band = np.logical_and(freq >= 0.3, freq <= 15) \n",
    "    power_in_band = Pxx[freq_band] \n",
    "    freq_in_band = freq[freq_band] \n",
    "    #plt.plot(freq,Pxx)\n",
    "\n",
    "    # D and F reflects impact of strides on acceleration\n",
    "    D = np.max(power_in_band) # Feature D: the maximum energy captured by the accelerator, \n",
    "    total_power = np.sum(power_in_band) # Feature F: total power in the frequencies between 0.3 and 15 Hz:\n",
    "\n",
    "    norm_Pxx = Pxx / total_power # normalize the power spectrum\n",
    "    E = entropy(norm_Pxx) # Feature E: normalized information entropy of the discrete FFT component magnitudes\n",
    "\n",
    "    \n",
    "\n",
    "    sorted_idx = np.argsort(power_in_band)[::-1] \n",
    "    first_freq = freq_in_band[sorted_idx[0]] \n",
    "    second_freq = freq_in_band[sorted_idx[1]] \n",
    "    first_power = power_in_band[sorted_idx[0]] \n",
    "    second_power = power_in_band[sorted_idx[1]]\n",
    "\n",
    "    R1 = np.sum(power_in_band[freq_in_band  < 3]) / total_power\n",
    "    R3 = np.sum(Pxx[(freq >= 1.5) & (freq <= 2.5)]) / total_power \n",
    "\n",
    "    # ----MOVING VS: STANDING----\n",
    "    moving = False\n",
    "    # Append the features to the list\n",
    "    if acc_mag_mean > 0.04:\n",
    "        moving = True\n",
    "    all_feature = [moving, acc_mag_mean ,acc_mag_std, A, B, C, Am, Bm, Cm, D, E, G, H, I, total_power, first_freq, first_power]\n",
    "    return all_feature\n",
    "\n",
    "def get_features_location_trace(trace):\n",
    "    std_win = 10 # length of window in seconds\n",
    "    sampling_rate = 200\n",
    "    ax = get_segment_overlap(trace.data['ax'].values, sampling_rate, std_win)\n",
    "    ay = get_segment_overlap(trace.data['ay'].values, sampling_rate, std_win)\n",
    "    az = get_segment_overlap(trace.data['az'].values, sampling_rate, std_win)\n",
    "    gx = get_segment_overlap(trace.data['gx'].values, sampling_rate, std_win)\n",
    "    gy = get_segment_overlap(trace.data['gy'].values, sampling_rate, std_win)\n",
    "    gz = get_segment_overlap(trace.data['gz'].values, sampling_rate, std_win)\n",
    "    # add feature extraction\n",
    "    num_features = 17\n",
    "    num_windows = np.shape(ax)[0]\n",
    "    featured_trace_location = np.zeros((num_windows, num_features))\n",
    "    for i in range(num_windows):\n",
    "        featured_trace_location[i,] = get_feature_location(ax[i], ay[i], az[i], gx[i], gy[i], gz[i])\n",
    "    return featured_trace_location"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# resample data to size: fixed length\n",
    "def resample_data(data, fixed_length):\n",
    "    num_data_points = len(data)\n",
    "    x_old = np.linspace(0, 1, num_data_points)\n",
    "    x_new = np.linspace(0, 1, fixed_length)\n",
    "    \n",
    "    # 1D interpolation for the magnitude data\n",
    "    interpolator = interp1d(x_old, data, axis=0, kind='linear')\n",
    "    resampled_data = interpolator(x_new)\n",
    "    return resampled_data\n",
    "\n",
    "def get_features_path_idx_trace(trace):\n",
    "    resample_length = 1000\n",
    "    mx = resample_data(trace.data['phone_mx'].values, resample_length)\n",
    "    my = resample_data(trace.data['phone_my'].values, resample_length)\n",
    "    mz = resample_data(trace.data['phone_mz'].values, resample_length)\n",
    "    magn_mag = [np.sqrt(m**2+my[i]**2+mz[i]**2)for i, m in enumerate(mx)]\n",
    "    return magn_mag"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### load testing trace and save features of each trace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_dir = 'E:\\\\Sunzhichao\\\\ETHz\\\\2223Spring\\\\Mobile_Health\\\\data\\\\test\\\\'\n",
    "\n",
    "def save_testing_features(dir_traces):\n",
    "    filenames = [join(dir_traces, f) for f in listdir(dir_traces) if isfile(join(dir_traces, f)) and f[-5:] == '.json']\n",
    "    filenames.sort()\n",
    "    features_all = {}\n",
    "    for filename in filenames:\n",
    "        trace_id = ''.join([*filename][-8:-5])\n",
    "        trace_feature = {}\n",
    "        trace = Recording(filename, no_labels=True, mute=True)\n",
    "        trace_feature[\"activity\"] = np.array(get_features_activity_trace(trace))\n",
    "        trace_feature[\"location\"] = np.array(get_features_location_trace(trace))\n",
    "        trace_feature[\"path_idx\"] = np.array(get_features_path_idx_trace(trace))\n",
    "        features_all[trace_id] = trace_feature\n",
    "    np.save(\"group00_features_test.npy\", features_all)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get and save features in test set\n",
    "save_testing_features(test_dir)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train and save models"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load training data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def segment(raw_data, all_data, data_count, std_win):\n",
    "    # Calculate window size\n",
    "    sampling_rate = 200\n",
    "    # std_win = 3 #s\n",
    "    window_size = round(std_win*sampling_rate)\n",
    "    prev_data_count = len(all_data)\n",
    "    for s in range(0, len(raw_data)-window_size, round(window_size/2)):\n",
    "        all_data.append(raw_data[s:s+window_size])\n",
    "    data_count.append(len(all_data)-prev_data_count)\n",
    "\n",
    "# one hot encoding for activity labels\n",
    "def one_hot(activity_labels):\n",
    "    activitu_label_onehot = [0, 0, 0, 0]\n",
    "    for i in range(4):\n",
    "        if i in activity_labels:\n",
    "            activitu_label_onehot[i] = 1\n",
    "    return activitu_label_onehot\n",
    "\n",
    "train_dir = 'E:\\\\Sunzhichao\\\\ETHz\\\\2223Spring\\\\Mobile_Health\\\\data\\\\train\\\\'\n",
    "filenames = [join(train_dir, f) for f in listdir(train_dir) if isfile(join(train_dir, f)) and f!='.DS_Store']\n",
    "filenames.sort()\n",
    "# segmented accelerometer magnitude data for activity\n",
    "magn_acc_seg = []\n",
    "data_count_activity = []\n",
    "# segmented accelerometer and gyroscope axis data for location\n",
    "x_acc_seg = []\n",
    "y_acc_seg = []\n",
    "z_acc_seg = []\n",
    "data_count_x_acc = []\n",
    "data_count_y_acc = []\n",
    "data_count_z_acc = []\n",
    "x_gyro_seg = []\n",
    "y_gyro_seg = []\n",
    "z_gyro_seg = []\n",
    "data_count_x_gyro = []\n",
    "data_count_y_gyro = []\n",
    "data_count_z_gyro = []\n",
    "\n",
    "# resampled data for path index\n",
    "magn_magnet_phone = []\n",
    "\n",
    "# labels\n",
    "activity_labels = []\n",
    "path_labels = []\n",
    "loaction_labels = []\n",
    "for i, filename in enumerate(filenames):\n",
    "    resample_length = 1000\n",
    "    trace = Recording(filename, no_labels=False, mute=True)\n",
    "    # Calculate raw magnitude of accelerometer signal\n",
    "    amagn_acc = [sqrt(a**2+trace.data['ay'].values[i]**2+trace.data['az'].values[i]**2)for i, a in enumerate(trace.data['ax'].values)]\n",
    "    # Pre-process data\n",
    "    amagn_acc = pre_process(amagn_acc)\n",
    "    segment(amagn_acc, magn_acc_seg, data_count_activity, 3)\n",
    "    x_acc = pre_process(trace.data['ax'].values)\n",
    "    segment(x_acc, x_acc_seg, data_count_x_acc, 10)\n",
    "    y_acc = pre_process(trace.data['ay'].values)\n",
    "    segment(y_acc, y_acc_seg, data_count_y_acc, 10)\n",
    "    z_acc = pre_process(trace.data['az'].values)\n",
    "    segment(z_acc, z_acc_seg, data_count_z_acc, 10)\n",
    "    x_gyro = pre_process(trace.data['gx'].values)\n",
    "    segment(x_gyro, x_gyro_seg, data_count_x_gyro, 10)\n",
    "    y_gyro = pre_process(trace.data['gy'].values)\n",
    "    segment(y_gyro, y_gyro_seg, data_count_y_gyro, 10)\n",
    "    z_gyro = pre_process(trace.data['gz'].values)\n",
    "    segment(z_gyro, z_gyro_seg, data_count_z_gyro, 10)\n",
    "    # Calculate raw x, y, z of PHONE magnetometer  signal\n",
    "    x_magnet_phone_resample = resample_data(trace.data['phone_mx'].values, resample_length)\n",
    "    y_magnet_phone_resample = resample_data(trace.data['phone_my'].values, resample_length)\n",
    "    z_magnet_phone_resample = resample_data(trace.data['phone_mz'].values, resample_length)\n",
    "    magn_mag_phone = [sqrt(m**2+y_magnet_phone_resample[i]**2+z_magnet_phone_resample[i]**2)for i, m in enumerate(x_magnet_phone_resample)]\n",
    "    magn_magnet_phone.append(magn_mag_phone)\n",
    "\n",
    "    # one hot encoding\n",
    "    activity_labels.append(one_hot(trace.labels['activities']))\n",
    "    path_labels.append(trace.labels['path_idx'])\n",
    "    loaction_labels.append(trace.labels['board_loc'])\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### model for activity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training XGBoost classifier for activity\n",
      "accuracy:  0.9470046082949308\n"
     ]
    }
   ],
   "source": [
    "# Load accelerometer data and labels for windows of first 50 traces\n",
    "activity_labels_training = np.load('activity_labels_training.npy')\n",
    "magn_data_training = np.load('activity_data_training.npy')\n",
    "# add feature extraction\n",
    "num_features = 25\n",
    "featured_data = np.zeros((magn_data_training.shape[0], num_features))\n",
    "for i in range(magn_data_training.shape[0]):\n",
    "    featured_data[i,] = get_features_activity_win(magn_data_training[i])\n",
    "# training \n",
    "X_train, X_test, y_train, y_test = train_test_split(featured_data, activity_labels_training, test_size = 0.2, shuffle = True)\n",
    "# Compute the class weights based on the frequency of each class in the training set\n",
    "class_weights = compute_class_weight(class_weight='balanced', classes=np.unique(y_train), y=y_train)\n",
    "# Set the parameters for the XGBoost model\n",
    "params = {'max_depth': 5, 'eta': 0.3, 'objective': 'multi:softmax', 'num_class': 4}\n",
    "# Create the XGBoost DMatrix object\n",
    "dtrain = xgb.DMatrix(X_train, label=y_train)\n",
    "# Train the model\n",
    "xgb_model = xgb.train(params, dtrain)\n",
    "print(\"Training XGBoost classifier for activity\")\n",
    "# Create the XGBoost DMatrix object for the test data\n",
    "dtest = xgb.DMatrix(X_test)\n",
    "# Make predictions on the test set and evaluate the model\n",
    "y_pred = xgb_model.predict(dtest)\n",
    "print(\"accuracy: \", accuracy_score(y_test, y_pred))\n",
    "# Save the model to a file\n",
    "# joblib.dump(xgb_model, 'activity_xgboost_model.joblib')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### model for location"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training XGBoost classifier for activity\n",
      "accuracy:  0.9625975217989904\n"
     ]
    }
   ],
   "source": [
    "data_count_location = data_count_x_acc\n",
    "location_labels_expand = []\n",
    "\n",
    "for i in range(len(loaction_labels)):\n",
    "    for j in range(data_count_location[i]):\n",
    "        location_labels_expand.append(loaction_labels[i])\n",
    "# add feature extraction\n",
    "num_features = 17\n",
    "num_windows = np.shape(x_acc_seg)[0]\n",
    "featured_data = np.zeros((num_windows, num_features))\n",
    "for i in range(num_windows):\n",
    "    featured_data[i,] = get_feature_location(x_acc_seg[i], y_acc_seg[i], z_acc_seg[i], x_gyro_seg[i], y_gyro_seg[i], z_gyro_seg[i])\n",
    "\n",
    "# train model\n",
    "X_train, X_test, y_train, y_test = train_test_split(featured_data, location_labels_expand, test_size = 0.3, shuffle = True)\n",
    "# Compute the class weights based on the frequency of each class in the training set\n",
    "class_weights = compute_class_weight(class_weight='balanced', classes=np.unique(y_train), y=y_train)\n",
    "# Set the parameters for the XGBoost model\n",
    "params = {'max_depth': 5, 'eta': 0.7, 'objective': 'multi:softmax', 'num_class': 3}\n",
    "# Create the XGBoost DMatrix object\n",
    "dtrain = xgb.DMatrix(X_train, label=y_train)\n",
    "# Train the model\n",
    "xgb_model = xgb.train(params, dtrain)\n",
    "print(\"Training XGBoost classifier for activity\")\n",
    "# Create the XGBoost DMatrix object for the test data\n",
    "dtest = xgb.DMatrix(X_test)\n",
    "# Make predictions on the test set and evaluate the model\n",
    "y_pred = xgb_model.predict(dtest)\n",
    "print(\"accuracy: \", accuracy_score(y_test, y_pred))\n",
    "# Save the model to a file\n",
    "# joblib.dump(xgb_model, 'location_xgboost_model.joblib')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### model for path index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the CNN model for data of length 1000\n",
    "class CNN_1000(nn.Module):\n",
    "    def __init__(self, num_classes):\n",
    "        super(CNN_1000, self).__init__()\n",
    "        self.conv1 = nn.Conv1d(1, 32, kernel_size=5)\n",
    "        self.relu1 = nn.ReLU()\n",
    "        self.maxpool1 = nn.MaxPool1d(kernel_size=2)\n",
    "        self.conv2 = nn.Conv1d(32, 64, kernel_size=5)\n",
    "        self.relu2 = nn.ReLU()\n",
    "        self.maxpool2 = nn.MaxPool1d(kernel_size=2)\n",
    "        self.dropout = nn.Dropout(p=0.5)  # Add dropout layer with 50% dropout probability\n",
    "        self.fc1 = nn.Linear(64 * 247, 255)\n",
    "        self.relu3 = nn.ReLU()\n",
    "        self.fc2 = nn.Linear(255, num_classes)\n",
    "\n",
    "    def forward(self, x):\n",
    "        x = self.conv1(x)\n",
    "        x = self.relu1(x)\n",
    "        x = self.maxpool1(x)\n",
    "        x = self.conv2(x)\n",
    "        x = self.relu2(x)\n",
    "        x = self.maxpool2(x)\n",
    "        x = self.dropout(x)  # Apply dropout\n",
    "        x = x.view(x.size(0), -1)\n",
    "        x = self.fc1(x)\n",
    "        x = self.relu3(x)\n",
    "        x = self.fc2(x)\n",
    "        return x\n",
    "    \n",
    "# Define a custom dataset class\n",
    "class PathDataset(Dataset):\n",
    "    def __init__(self, X, y):\n",
    "        self.X = X\n",
    "        self.y = y\n",
    "\n",
    "    def __len__(self):\n",
    "        return len(self.X)\n",
    "\n",
    "    def __getitem__(self, idx):\n",
    "        return self.X[idx], self.y[idx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert resampled data and labels to tensors\n",
    "X = torch.tensor(np.expand_dims(magn_magnet_phone, axis=-1), dtype=torch.float32)\n",
    "y = torch.tensor(path_labels, dtype=torch.long)\n",
    "# Split the data into training and validation sets\n",
    "X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=42)\n",
    "\n",
    "# Create dataset and dataloader objects for training and validation sets\n",
    "train_dataset = PathDataset(X_train, y_train)\n",
    "val_dataset = PathDataset(X_val, y_val)\n",
    "train_loader = DataLoader(train_dataset, batch_size=32, shuffle=True)\n",
    "val_loader = DataLoader(val_dataset, batch_size=1, shuffle=False)\n",
    "# Set device\n",
    "device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n",
    "\n",
    "# Initialize the CNN model and other training components\n",
    "num_classes = 5\n",
    "model = CNN_1000(num_classes).to(device)\n",
    "\n",
    "criterion = nn.CrossEntropyLoss()\n",
    "optimizer = optim.Adam(model.parameters(), lr=4e-06)\n",
    "\n",
    "# Training parameters\n",
    "num_epochs = 1000\n",
    "best_val_accuracy = 0\n",
    "\n",
    "# Training loop\n",
    "for epoch in range(num_epochs):\n",
    "    train_loss = 0.0\n",
    "    model.train()\n",
    "    correct_train = 0\n",
    "    total_train = 0\n",
    "\n",
    "    for inputs, targets in train_loader:\n",
    "        inputs, targets = inputs.to(device), targets.to(device)\n",
    "        inputs = inputs.permute(0, 2, 1)\n",
    "        optimizer.zero_grad()\n",
    "        outputs = model(inputs)\n",
    "        # Compute the number of correct predictions and the total number of predictions\n",
    "        _, predicted = torch.max(outputs.data, 1)\n",
    "        total_train += targets.size(0)\n",
    "        correct_train += (predicted == targets).sum().item()\n",
    "        # Compute the total loss, including the L1 regularization term\n",
    "        loss = criterion(outputs, targets) #+ l1_regularization(model, lambda_l1)\n",
    "        \n",
    "        \n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "\n",
    "        train_loss += loss.item() * inputs.size(0)\n",
    "\n",
    "    train_loss = train_loss / len(train_loader.dataset)\n",
    "    train_accuracy = 100 * correct_train / total_train\n",
    "    \n",
    "    # Validation loop\n",
    "    model.eval()\n",
    "    correct = 0\n",
    "    total = 0\n",
    "\n",
    "    with torch.no_grad():\n",
    "        for inputs, targets in val_loader:\n",
    "            inputs, targets = inputs.to(device), targets.to(device)\n",
    "            \n",
    "            inputs = inputs.permute(0, 2, 1)\n",
    "\n",
    "            outputs = model(inputs)\n",
    "            _, predicted = torch.max(outputs.data, 1)\n",
    "            total += targets.size(0)\n",
    "            correct += (predicted == targets).sum().item()\n",
    "\n",
    "    val_accuracy = 100 *correct / total\n",
    "\n",
    "    # Save the model with the best validation accuracy\n",
    "    if val_accuracy > best_val_accuracy:\n",
    "        best_val_accuracy = val_accuracy\n",
    "        model = model.to(\"cpu\")\n",
    "        # torch.save(model, \"path_index_cnn_model.pt\")\n",
    "        model = model.to(device)\n",
    "\n",
    "    print(f\"Epoch {epoch + 1}/{num_epochs}, Loss: {train_loss:.4f},Train Accuracy: {train_accuracy:.2f}%, Validation Accuracy: {val_accuracy:.2f}%\")\n",
    "\n",
    "print(f\"Best validation accuracy: {best_val_accuracy:.4f}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "mhealth23",
   "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.12"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}