{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# LOAD LIBRARIES\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nfrom sklearn.metrics import average_precision_score\nfrom scipy.optimize import minimize_scalar\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.cluster import KMeans\nfrom scipy import signal\nimport os\nimport matplotlib.pyplot as plt\nfrom matplotlib import mlab\nimport torch\nimport torch.nn as nn\nimport torch.optim as optim\nimport torch.utils.data as thd\nfrom torch.utils.data import Dataset, DataLoader\nfrom torch.nn.utils.rnn import pad_sequence, pack_padded_sequence\nimport torch.nn.functional as F\nfrom tqdm import tqdm\nimport pdb\nimport sys\nimport time\nimport line_profiler\nfrom skimage.restoration import denoise_tv_bregman\nimport scipy.ndimage as ndimage\nimport warnings\n\n\nDEVICE = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-07-06T17:44:36.780738Z","iopub.execute_input":"2023-07-06T17:44:36.781048Z","iopub.status.idle":"2023-07-06T17:44:40.593902Z","shell.execute_reply.started":"2023-07-06T17:44:36.781017Z","shell.execute_reply":"2023-07-06T17:44:40.592761Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# SET UP SCRIPT PARAMETERS AND DIRECTORIES\n\n# script parameters\nn_folds = 6\nsep_models = True  # if true then separate models for defog and tdcs data\n\n# directory paths\nhome_dir = '/kaggle/input/tlvmc-parkinsons-freezing-gait-prediction/'\ntrain_dir = home_dir + 'train/'\ntest_dir = home_dir + 'test/'\n\n# feature engineering\nwindow_size = 200\nwindow_future = 75\nwindow_past = window_size - window_future\nacc_cols = ['AccV','AccML','AccAP']\nfs_dict = {'tdcsfog':128, 'defog':100}\nNFFT = 8","metadata":{"execution":{"iopub.status.busy":"2023-07-06T17:44:40.596993Z","iopub.execute_input":"2023-07-06T17:44:40.597927Z","iopub.status.idle":"2023-07-06T17:44:40.605108Z","shell.execute_reply.started":"2023-07-06T17:44:40.597882Z","shell.execute_reply":"2023-07-06T17:44:40.603866Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# FOCAL LOSS FROM CHATGPT\nclass FocalLoss(nn.Module):\n    def __init__(self, alpha=None, gamma=2, reduction='mean'):\n        super(FocalLoss, self).__init__()\n        self.alpha = alpha\n        self.gamma = gamma\n        self.reduction = reduction\n\n    def forward(self, inputs, targets):\n        ce_loss = F.cross_entropy(inputs, targets, reduction='none')\n        pt = torch.exp(-ce_loss)\n        if self.alpha is not None:\n            alpha_t = self.alpha[targets]\n            focal_loss = alpha_t * (1 - pt) ** self.gamma * ce_loss\n        else:\n            focal_loss = (1 - pt) ** self.gamma * ce_loss\n\n        if self.reduction == 'mean':\n            return focal_loss.mean()\n        elif self.reduction == 'sum':\n            return focal_loss.sum()\n        else:\n            return focal_loss","metadata":{"execution":{"iopub.status.busy":"2023-07-06T17:44:40.606497Z","iopub.execute_input":"2023-07-06T17:44:40.607134Z","iopub.status.idle":"2023-07-06T17:44:40.619352Z","shell.execute_reply.started":"2023-07-06T17:44:40.607086Z","shell.execute_reply":"2023-07-06T17:44:40.618331Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# interpolated average precision loss\ndef iap_loss(y_pred, y_true):\n    # y_pred: tensor of shape (batch_size, num_classes)\n    # y_true: tensor of shape (batch_size, num_classes), one-hot encoded\n    \n#     y_true = torch.nn.functional.one_hot(y_true.squeeze()))\n    \n    batch_size, num_classes = y_pred.shape\n    eps = 1e-7\n    \n    # Compute precision, recall, and average precision\n    tp = (y_pred * y_true).sum(dim=0) # true positives\n    fp = ((1 - y_true) * y_pred).sum(dim=0) # false positives\n    fn = (y_true * (1 - y_pred)).sum(dim=0) # false negatives\n    precision = tp / (tp + fp + eps)\n    recall = tp / (tp + fn + eps)\n    ap = torch.mean(precision[torch.argsort(recall, descending=True)])\n    \n    # Compute interpolated precision and loss\n    interp_precision = torch.zeros((num_classes, 11), device=y_pred.device)\n    interp_recall = torch.linspace(0, 1, 11, device=y_pred.device)\n    for i in range(num_classes):\n        for j in range(10, -1, -1):\n            interp_precision[i, j] = torch.max(precision[i:][recall[i:] >= interp_recall[j]])\n    interp_ap = torch.mean(interp_precision[:, 0])\n    loss = 1 - interp_ap\n    \n    return loss","metadata":{"execution":{"iopub.status.busy":"2023-07-06T17:44:40.623122Z","iopub.execute_input":"2023-07-06T17:44:40.623557Z","iopub.status.idle":"2023-07-06T17:44:40.634481Z","shell.execute_reply.started":"2023-07-06T17:44:40.623528Z","shell.execute_reply":"2023-07-06T17:44:40.633299Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# FEATURE ENGINEERING AND ANALYSIS FUNCTIONS\n\n# MEAN ABSOLUTE PRECISION\ndef mAP(preds, actuals, num_grps):\n    return [average_precision_score(y_score=preds[grp,], y_true=actuals[grp,]) for grp in range(num_grps)]\n\ndef normalize_columns(arr):\n    means = np.mean(arr, axis=0)\n    stds = np.std(arr, axis = 0)\n    stds[np.where(stds == 0)] == 1\n    \n    return (arr-means)/stds\n\n# TV FILTER\ndef total_variance_filter(y, lmbda):\n    n = len(y)\n    def objective(x):\n        return np.sum((x - y)**2) + lmbda * np.sum(np.abs(np.diff(x)))\n    x0 = np.copy(y)\n    result = minimize_scalar(objective)\n    return result.x\n\n# MA FILTER\ndef hp_lp_filter(y, size=75):\n    y_ma = ndimage.uniform_filter1d(y, size=75)\n    y_hp = y - y_ma\n    return np.column_stack([y_hp,y_ma])\n\n# SPECTROGRAM\ndef spec_interp(y, fs, nfft=64):\n    interp_inds = np.array(range(len(y)))/(len(y)-1)\n    # resample 100hz data to 120hz to align spectrograms\n    if fs == 100:\n        t_new = np.arange(0, len(y), 100 / 120)\n        y = np.interp(t_new, np.arange(len(y)), y)\n        fs = 120\n    \n    spec, freqs, t = mlab.specgram(y, NFFT=nfft, Fs=fs, noverlap=nfft//2, detrend='linear')\n    spec = 10 * np.log10(spec + 1e-25)\n    t_inds = np.array(range(len(t)))/(len(t)-1)\n    spec_db_interp = np.array([np.interp(interp_inds, t_inds, spec[i,:]) for i in range(len(spec))])\n    \n    return spec_db_interp\n\n# FILE CLUSTERING\n","metadata":{"execution":{"iopub.status.busy":"2023-07-06T17:44:40.636216Z","iopub.execute_input":"2023-07-06T17:44:40.637089Z","iopub.status.idle":"2023-07-06T17:44:40.654350Z","shell.execute_reply.started":"2023-07-06T17:44:40.637051Z","shell.execute_reply":"2023-07-06T17:44:40.651631Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# CUSTOM DATASET CLASS FOR ACCELEROMETER DATA\nclass AccDataset(Dataset):\n    def __init__(self, files_df, patients_df, valid_protocols=['tdcsfog'], tr=True):\n        self.patients_df = patients_df\n        self.valid_protocols = valid_protocols\n        self.files_df = files_df[files_df['Protocol'].isin(self.valid_protocols)] \n        self.files_df.reset_index(inplace=True, drop=True)\n        self.cv_col = cv_col\n        self.tr = tr\n\n        self.acc_feats = {}\n        self.labels = {}\n        self.tr_items = []\n        print('Loading and Processing Data:')\n        for idx in range(len(self.files_df)):\n            if idx % 100 == 0:\n                print(\"{}/{}\".format(idx,len(self.files_df)))\n            filepath = self.files_df.at[idx,'Filepath']\n            protocol = self.files_df.at[idx,'Protocol']\n            acc_feats, labels, tr_items = self.load_and_process(filepath, protocol)\n            self.acc_feats[filepath] = acc_feats\n            self.labels[filepath] = labels\n            self.tr_items.append(tr_items)\n                \n        self.tr_items = pd.concat(self.tr_items)\n        self.feature_dim = np.shape(acc_feats)[0]\n        \n        if not tr:\n            # do initialization steps for inference that would otherwise happen in set_cv_idx\n            self.cv_subjects = self.patients_df.Subject\n            self.cv_files_df = self.files_df\n            self.cv_tr_items = self.tr_items.to_numpy()\n\n    def __len__(self):\n        return len(self.cv_tr_items)\n    \n    def __getitem__(self, idx):\n        return self.process_sample(idx)\n    \n    def set_cv_idx(self, cv_idx=[None], val=False, protocols=['defog','tdcsfog'], cv_col=['combo_cv'], return_idx=False):\n        self.cv_idx = cv_idx\n        self.protocols = protocols\n        self.cv_col = cv_col\n        self.return_idx = return_idx\n        \n        if val:\n            self.cv_subjects = self.patients_df[self.patients_df[self.cv_col].isin(cv_idx).values].Subject\n        else:\n            self.cv_subjects = self.patients_df[[not s for s in self.patients_df[self.cv_col].isin(cv_idx).values]].Subject\n            \n        # training files for this cv    \n        self.cv_files_df = self.files_df[self.files_df['Subject'].isin(self.cv_subjects)]\n        self.cv_files_df = self.cv_files_df[self.cv_files_df['Protocol'].isin(protocols)] \n        \n        # training samples for this cv\n        self.cv_tr_items = self.tr_items[self.tr_items['Filepath'].isin(self.cv_files_df.Filepath)].to_numpy()\n        \n        # label frequency for this cv\n        self.lab_freq = np.zeros(4)\n        for k,v in self.labels.items():\n            if k in self.cv_files_df['Filepath'].values:\n                unique, counts = np.unique(v, return_counts=True)\n                for i in range(len(unique)):\n                    self.lab_freq[unique[i]] += counts[i]\n        self.lab_freq = np.sum(self.lab_freq)/self.lab_freq\n        self.lab_freq = self.lab_freq/np.sum(self.lab_freq)\n    \n    def load_and_process(self, filepath, protocol):\n        sample_df = pd.read_csv(filepath)\n        \n        if self.tr:\n            if protocol == 'defog':\n                # create 4th class for no FOG, and remove other class labels unless Valid and Task are both True\n                sample_df['StartHesitation'] = sample_df.apply(lambda x: x['StartHesitation'] if x['Valid'] == True and x['Task'] == True else 0, axis=1)\n                sample_df['Turn'] = sample_df.apply(lambda x: x['Turn'] if x['Valid'] == True and x['Task'] == True else 0, axis=1)\n                sample_df['Walking'] = sample_df.apply(lambda x: x['Walking'] if x['Valid'] == True and x['Task'] == True else 0, axis=1)\n                sample_df['NoFOG'] = sample_df.apply(lambda x: 0 if x['StartHesitation'] == 1 or x['Turn'] == 1 or x['Walking'] == 1 else 1, axis=1)\n                \n            elif protocol == 'tdcsfog':\n                sample_df['NoFOG'] = sample_df.apply(lambda x: 0 if x['StartHesitation'] == 1 or x['Turn'] == 1 or x['Walking'] == 1 else 1, axis=1)\n\n            labels = torch.t(torch.tensor(sample_df[['StartHesitation','Turn','Walking','NoFOG']].values, dtype=torch.int64))  \n            labels = torch.argmax(labels,dim=0)#.to(torch.int64)\n            labels = torch.unsqueeze(labels,0)\n        else: \n            labels = torch.empty(1,(len(sample_df)), dtype=torch.int64)\n        \n        acc_feats = self.acc_gen_features(sample_df, protocol)\n        item_df = pd.DataFrame({'Filepath':filepath, 'Ind':list(range(len(sample_df)))})\n            \n        return acc_feats, labels, item_df\n    \n    def acc_gen_features(self, sample_df, protocol):\n        acc_feats = sample_df[['AccV','AccML','AccAP']].values\n        fs = fs_dict[protocol]\n        \n        # change units\n        if protocol is not 'tdcsfog':\n            # change units from g to m/s^2\n            sample_df[acc_cols] = sample_df[acc_cols].applymap(lambda x: x*9.807)\n        \n        # filter acceleration\n        acc_filt = np.column_stack([hp_lp_filter(acc_feats[:,i]) for i in range(3)])\n        acc_feats = np.column_stack([acc_feats, acc_filt])\n\n        # spectrogram\n        acc_spec = np.column_stack([np.transpose(spec_interp(acc_feats[:,i], fs, nfft=NFFT)) for i in range(3)])\n        acc_feats = np.column_stack([acc_feats, acc_spec])\n\n        # time features\n        t = np.array(range(len(acc_feats)))/fs\n        t_prop = np.array(range(len(acc_feats)))/(len(acc_feats)-1)\n        acc_feats = np.column_stack([acc_feats, t, t_prop])\n        \n        acc_feats = normalize_columns(acc_feats)\n        return torch.t(torch.tensor(acc_feats,  dtype=torch.float32))\n    \n    def process_sample(self, tr_idx):\n        # RAM is limiting factor, so only create feature/label tensor windows when a sample is explicitly requested\n        \n        fp = self.cv_tr_items[tr_idx]\n        acc_feats = self.acc_feats[fp[0]]\n        labels = self.labels[fp[0]]\n        idx = fp[1]\n        x_t = self.window_tensor(acc_feats, idx)\n        y_t = labels[:,idx:(idx+1)]\n\n        if self.tr and not self.return_idx:\n            return x_t, y_t\n        else:\n            k = self.cv_files_df[self.cv_files_df['Filepath'] == fp[0]].Id.values + '_' + str(idx)\n            return x_t, y_t, k.tolist() \n    \n    def window_tensor(self, x_t, pred_idx, w_width=window_size, w_future=window_future, w_past=window_past):\n        num_t = x_t.shape[1]\n        start_idx = max(0, pred_idx - w_past)\n        end_idx = min(num_t-1, pred_idx + w_future)\n\n        # if window extends outside data, pad with closest data \n        i_w = x_t[:,start_idx:end_idx]\n        if pred_idx-w_past < 0:\n            past_pad = x_t[:,0]\n            pad_tensor = past_pad.repeat(abs(pred_idx - w_past), 1)\n            i_w = torch.cat([pad_tensor.transpose(0,1),i_w], dim=1)\n#             past_pad = x_t[:,[0]]\n#             pad_tensor = np.repeat(past_pad, abs(pred_idx - w_past),1)\n#             i_w = np.concatenate((pad_tensor, i_w),axis=1)\n        elif pred_idx + w_future > (num_t - 1):\n            fut_pad = x_t[:,-1]\n            pad_tensor = fut_pad.repeat(pred_idx + w_future - (num_t - 1), 1)\n            i_w = torch.cat([i_w, pad_tensor.transpose(0,1)], dim=1)\n#             fut_pad = x_t[:,[-1]]\n#             pad_tensor = np.repeat(fut_pad, pred_idx + w_future - (np.shape(x_t)[1] - 1), 1)\n#             i_w = np.concatenate((i_w, pad_tensor),axis=1)\n                \n        return i_w\n","metadata":{"execution":{"iopub.status.busy":"2023-07-06T17:44:40.660110Z","iopub.execute_input":"2023-07-06T17:44:40.662609Z","iopub.status.idle":"2023-07-06T17:44:40.733500Z","shell.execute_reply.started":"2023-07-06T17:44:40.662573Z","shell.execute_reply":"2023-07-06T17:44:40.732324Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# PYTORCH MODEL 2\n# Similar to model 1, but no linear layers (I think those would be good for classifying the entire sequence, but bad for classifying each timestep)\nclass cnn(nn.Module):\n    def __init__(self, input_dim, conv_layers, num_classes, kernel_size, stride, padding, incl_pool=False, dropout_p=0.25, pool_k=2, pool_s=2):\n        super(cnn, self).__init__()\n        self.input_dim = input_dim\n        self.hidden_dim = conv_layers\n        self.num_classes = num_classes\n        self.kernel_size = kernel_size\n        self.stride = stride\n        self.padding = padding\n        self.pool_k = pool_k\n        self.pool_s = pool_s\n        self.pool_reduction = 1\n        \n            # Define the convolutional layers\n        self.conv_layers = nn.ModuleList()\n        self.out_layers = nn.ModuleList()\n        self.pool_layers = nn.ModuleList()\n        for i in range(len(conv_layers)):\n            in_channels = input_dim if i == 0 else conv_layers[i-1]\n            out_channels = conv_layers[i]\n            conv_layer = nn.Conv1d(in_channels, out_channels, kernel_size[i], stride[i], padding[i])\n            self.conv_layers.append(conv_layer)\n            self.conv_layers.append(nn.ReLU())\n            self.conv_layers.append(nn.BatchNorm1d(num_features=out_channels))\n            self.conv_layers.append(nn.Dropout(p=dropout_p[i]))\n            if incl_pool:\n                self.conv_layers.append(nn.MaxPool1d(kernel_size=pool_k, stride=pool_s, padding=0))\n#                 self.conv_layers.append(nn.AvgPool1d(kernel_size=pool_k, stride=pool_s, padding=int(np.floor(pool_k/2))))\n                self.pool_reduction = self.pool_reduction/pool_s\n            \n            self.final_dropout = nn.Dropout(p=dropout_p[-1])\n            self.global_avg_pooling = nn.AdaptiveAvgPool1d(output_size=1)\n            self.fc_layer = nn.Linear(conv_layers[-1], self.num_classes)\n    \n    def forward(self, x):\n        # get input sequence shapes\n        batch_size, input_dim, seq_len = x.size()\n        for i_layer in self.conv_layers:\n            x = i_layer(x)\n        \n        x = self.global_avg_pooling(x)\n        x = self.final_dropout(x)\n        x = x.view(batch_size,-1)\n        x = self.fc_layer(x)\n        x = nn.functional.softmax(x, dim=1)\n\n        return x\n        ","metadata":{"execution":{"iopub.status.busy":"2023-07-06T17:44:40.735292Z","iopub.execute_input":"2023-07-06T17:44:40.736198Z","iopub.status.idle":"2023-07-06T17:44:40.760348Z","shell.execute_reply.started":"2023-07-06T17:44:40.736149Z","shell.execute_reply":"2023-07-06T17:44:40.759202Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class Model(nn.Module):\n    def __init__(self, input_dim, conv_layers, num_classes, kernel_size, padding, dropout_p=0.25):\n        super(Model, self).__init__()\n        self.conv_blocks = nn.Sequential()\n        for _ in range(len(conv_layers)):\n            in_channels = input_dim if _ == 0 else conv_layers[_-1]\n            self.conv_blocks.add_module(\"conv1d_\" + str(_), nn.Conv1d(in_channels=in_channels, out_channels=conv_layers[_], kernel_size=kernel_size, padding=padding))\n            self.conv_blocks.add_module(\"batchnorm_\" + str(_), nn.BatchNorm1d(conv_layers[_]))\n            self.conv_blocks.add_module(\"relu_\" + str(_), nn.ReLU())\n            self.conv_blocks.add_module(\"dropout_\" + str(_), nn.Dropout(dropout_p))\n\n        self.global_avg_pooling = nn.AdaptiveAvgPool1d(output_size=1)\n        self.fc = nn.Linear(conv_layers[_], num_classes)\n\n    def forward(self, x):\n        #x = x.permute(0, 2, 1)  # swap dimensions for Conv1D to work\n        x = self.conv_blocks(x)\n        x = self.global_avg_pooling(x)\n        x = x.view(x.size(0), -1)\n        x = self.fc(x)\n\n        return x","metadata":{"execution":{"iopub.status.busy":"2023-07-06T17:44:40.764331Z","iopub.execute_input":"2023-07-06T17:44:40.765081Z","iopub.status.idle":"2023-07-06T17:44:40.781102Z","shell.execute_reply.started":"2023-07-06T17:44:40.765008Z","shell.execute_reply":"2023-07-06T17:44:40.780054Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# CREATE DATAFRAME OF METADATA FOR TRAINING FILES\n\n# build dataframe with training files\ntr_types = ['defog','notype','tdcsfog']\ntr_type_dir = [train_dir + i for i in tr_types]\ntr_files_df = []\nfor dirname, _, filenames in os.walk(train_dir):\n    for filename in filenames:\n        #print(os.path.join(dirname, filename))\n        tr_files_df.append(\n            {\n            'Filename': filename,\n            'Filepath': os.path.join(dirname, filename),\n            'Protocol': tr_types[tr_type_dir.index(dirname)]\n            }\n        )\ntr_files_df = pd.DataFrame(tr_files_df)\n\n# add metadata for different protocols\ntdcs_meta = pd.read_csv(home_dir + 'tdcsfog_metadata.csv')\ndefog_meta = pd.read_csv(home_dir + 'defog_metadata.csv')\ndefog_meta['Test'] = 0 # add dummy 'Test' column so it can be combined with tdcs_meta, whose unique test values are (1,2,3)\nall_meta = pd.concat([tdcs_meta,defog_meta])\nall_meta['Filename'] = [i + '.csv' for i in all_meta.Id]\n\n# files df\ntr_files_df = pd.merge(tr_files_df,all_meta,how='left',on='Filename')","metadata":{"execution":{"iopub.status.busy":"2023-07-06T17:44:40.783001Z","iopub.execute_input":"2023-07-06T17:44:40.783422Z","iopub.status.idle":"2023-07-06T17:44:41.216316Z","shell.execute_reply.started":"2023-07-06T17:44:40.783345Z","shell.execute_reply":"2023-07-06T17:44:41.215309Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# CREATE DATAFRAME OF METADATA FOR TEST FILES\n\n# build dataframe with test files\ntest_types = ['defog','tdcsfog']\ntest_type_dir = [test_dir + i for i in tr_types]\ntest_files_df = []\nfor dirname, _, filenames in os.walk(test_dir):\n    for filename in filenames:\n        test_files_df.append(\n            {\n            'Filename': filename,\n            'Filepath': os.path.join(dirname, filename),\n            'Protocol': tr_types[test_type_dir.index(dirname)]\n            }\n        )\ntest_files_df = pd.DataFrame(test_files_df)\n\n# files df\ntest_files_df = pd.merge(test_files_df,all_meta,how='left',on='Filename')","metadata":{"execution":{"iopub.status.busy":"2023-07-06T17:44:41.220875Z","iopub.execute_input":"2023-07-06T17:44:41.221180Z","iopub.status.idle":"2023-07-06T17:44:41.240012Z","shell.execute_reply.started":"2023-07-06T17:44:41.221152Z","shell.execute_reply":"2023-07-06T17:44:41.239047Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# CREATE DATAFRAME OF METADATA FOR PATIENTS\npatients_df = pd.read_csv(home_dir + 'subjects.csv')\n\n# get count of number/type of files by patient\nfile_counts = tr_files_df.groupby(['Subject','Protocol']).size().reset_index(name='Count')\nfile_counts = file_counts.pivot_table(values = 'Count', index='Subject', columns='Protocol',fill_value =0)\nfile_counts['total_files'] = file_counts[['defog','notype','tdcsfog']].sum(axis=1)\nfile_counts['total_typed'] = file_counts[['defog','tdcsfog']].sum(axis=1)\n\n# patient df\npatients_df = pd.merge(patients_df,file_counts,how='left',on='Subject')\npatients_df[['defog','notype','tdcsfog','total_files','total_typed']] = patients_df[['defog','notype','tdcsfog','total_files','total_typed']].fillna(0)","metadata":{"execution":{"iopub.status.busy":"2023-07-06T17:44:41.242341Z","iopub.execute_input":"2023-07-06T17:44:41.243170Z","iopub.status.idle":"2023-07-06T17:44:41.288685Z","shell.execute_reply.started":"2023-07-06T17:44:41.243129Z","shell.execute_reply":"2023-07-06T17:44:41.287613Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"visit_feats = ['Visit','Test','Medication','Protocol']\npatient_feats = ['Age','Sex','YearsSinceDx','UPDRSIII_On','UPDRSIII_Off','NFOGQ']","metadata":{"execution":{"iopub.status.busy":"2023-07-06T17:44:41.291976Z","iopub.execute_input":"2023-07-06T17:44:41.292276Z","iopub.status.idle":"2023-07-06T17:44:41.299500Z","shell.execute_reply.started":"2023-07-06T17:44:41.292231Z","shell.execute_reply":"2023-07-06T17:44:41.298313Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# # DEFINE CV FOLDS AT THE PATIENT LEVEL\n\n# STRATIFIED K-FOLD ON PATIENTS + TOTAL FOG TIME\nevents_df = pd.read_csv(home_dir + 'events.csv')\nevents_df['Total_time'] = events_df['Completion'] - events_df['Init']\nevents_df = pd.merge(events_df, tr_files_df[['Id','Subject']], on='Id', how='inner')\nevents_agg = events_df.groupby(['Subject','Type']).agg({'Total_time': 'sum'}).reset_index()\nevents_pivot = events_agg.pivot_table(index='Subject', columns='Type', values='Total_time', fill_value=0).reset_index()\npatients_df = pd.merge(patients_df, events_pivot, on='Subject', how='left').fillna(0)\n\n# create separate fold table with unique rows per Subject\nfold_df = patients_df.groupby(['Subject']).agg({'StartHesitation': 'sum',\n                                                'Walking': 'sum',\n                                                'Turn': 'sum',\n                                                'defog': 'sum',\n                                                'tdcsfog': 'sum',\n                                                'notype': 'sum'}).reset_index()\n\n# create folds for defog model\nfold_df = fold_df.sort_values(by=['StartHesitation','Walking','Turn','defog'],ascending=True)\nfold_df['defog_cv'] = fold_df.reset_index().index % n_folds\n\n# create folds for tdcs model\nfold_df = fold_df.sort_values(by=['StartHesitation','Walking','Turn','tdcsfog'],ascending=True)\nfold_df['tdcsf_cv'] = fold_df.reset_index().index % n_folds\n\n# create folds for combined model\nfold_df = fold_df.sort_values(by=['StartHesitation','Walking','Turn','defog','tdcsfog'],ascending=True)\n# patients_df = patients_df.sort_values(by=['defog','tdcsfog'],ascending=True)\nfold_df['combo_cv'] = fold_df.reset_index().index % n_folds\n\n# create folds for combined model\nfold_df = fold_df.sort_values(by=['StartHesitation','Walking','Turn','defog','tdcsfog','notype'],ascending=True)\n# patients_df = patients_df.sort_values(by=['defog','tdcsfog'],ascending=True)\nfold_df['all_labeled_cv'] = fold_df.reset_index().index % n_folds\n\n# Add fold labels to patients_df\npatients_df = patients_df.merge(fold_df[['Subject','defog_cv','tdcsf_cv','combo_cv','all_labeled_cv']], on='Subject', how='left')\n\n# validate each subject is only in one fold\nvalue_counts = fold_df['Subject'].value_counts()\nvalue_counts = value_counts.sort_values(ascending=False)\nunique_folds = patients_df.groupby('Subject')['all_labeled_cv'].nunique()\nunique_folds.sort_values()\n\n# check number of events in each fold\nprint(\"Check \")\nfold_check = events_df.merge(fold_df[['Subject','all_labeled_cv']],on = 'Subject')\nfold_check.groupby(['all_labeled_cv','Type']).agg({'Total_time': 'sum'})","metadata":{"execution":{"iopub.status.busy":"2023-07-06T17:44:41.301052Z","iopub.execute_input":"2023-07-06T17:44:41.301425Z","iopub.status.idle":"2023-07-06T17:44:41.395333Z","shell.execute_reply.started":"2023-07-06T17:44:41.301385Z","shell.execute_reply":"2023-07-06T17:44:41.394124Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Space to protoype feature engineering\nif False:\n    protocols = ['tdcsfog']\n    cv_col = ['tdcsf_cv']\n    cv_idx = [0]\n    cv_dataset = AccDataset(files_df=tr_files_df, patients_df=patients_df, valid_protocols=protocols)\n\n    cv_dataset.set_cv_idx(cv_idx=cv_idx, val=True, protocols=protocols, cv_col=cv_col)\n    i_acc = cv_dataset.acc_feats['/kaggle/input/tlvmc-parkinsons-freezing-gait-prediction/train/tdcsfog/a171e61840.csv']\n    y = i_acc[2,]\n    x = list(range(len(y)))\n\n    # Apply total variation denoising\n    y_denoised = denoise_tv_bregman(y, weight=1000)\n    y_ma = ndimage.uniform_filter1d(y, size=75)\n\n    # # Plot the results\n    # import matplotlib.pyplot as plt\n    fig, ax = plt.subplots()\n    ax.plot(x, y, label='Noisy')\n    ax.plot(x, y_ma, label='Denoised')\n    ax.legend()\n    plt.show()\n\n    # plot difference\n    y_hp = y - y_ma\n    fig, ax = plt.subplots()\n    ax.plot(x, y_hp, label='Noisy')\n    plt.show()\n\n    # Generate a sample signal\n    time = list(range(len(y)))\n    #signal = #np.sin(2 * np.pi * 5 * time)\n\n    # Compute the spectrogram\n    fs = 120 # Sample rate\n    NFFT = 64  # Number of FFT points\n    spec, freqs, t = mlab.specgram(y_hp, NFFT=32, Fs=fs, noverlap=NFFT//2, detrend='linear')\n    spec_db = 10 * np.log10(spec)\n    \n    # Plot the spectrogram\n    plt.pcolormesh(t, freqs, spec_db)\n    plt.ylabel('Frequency [Hz]')\n    plt.xlabel('Time [s]')\n    plt.show()\n    \n    # interpolate\n    t_inds = np.array(range(len(t)))/(len(t)-1)\n    interp_inds = np.array(range(len(y)))/(len(y)-1)\n    spec_db_interp = np.array([np.interp(interp_inds, t_inds, spec[i,:]) for i in range(len(spec))])\n    \n    plt.pcolormesh(interp_inds, freqs, spec_db_interp)\n    plt.pcolormesh(spec_interp(y,120,32))\n    plt.ylabel('Frequency [Hz]')\n    plt.xlabel('Time [s]')\n    plt.show()\n    \n    lab_freq = np.zeros(4)\n    for k,v in cv_dataset.labels.items():\n        if k in cv_dataset.cv_files_df['Filepath'].values:\n            unique, counts = np.unique(v, return_counts=True)\n            for i in range(len(unique)):\n                lab_freq[unique[i]] += counts[i]\n    lab_freq = np.sum(lab_freq)/lab_freq\n    lab_freq = lab_freq/np.sum(lab_freq)","metadata":{"execution":{"iopub.status.busy":"2023-07-06T17:44:41.397152Z","iopub.execute_input":"2023-07-06T17:44:41.397974Z","iopub.status.idle":"2023-07-06T17:44:41.416289Z","shell.execute_reply.started":"2023-07-06T17:44:41.397933Z","shell.execute_reply":"2023-07-06T17:44:41.415075Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# FUNCTION TO PERFORM TRAINING LOOP (SO TDCSFOG AND DEFOG CAN HAVE SEPARATE MODELS)\ndef CV_TRAIN():\n    # Model Parameters\n    input_dim = cv_dataset.feature_dim\n    num_classes = 4\n    conv_layers = [32,64,64,128]\n    kernel_size = [9,9,2,1]\n    stride = [1,2,2,1]\n    padding = [0,0,0,0,0]#int(np.floor(kernel_size/2))\n    dropout_p=[0,0,0,0,0,0.6]\n    incl_pool=False\n    pool_k=3\n    pool_s=3\n\n    # Training Parameters\n    batch_size = 1024\n    learning_rate = 1e-4#1.2e-4 \n    num_epochs = 1\n    max_batches_per_epoch = 12000\n    max_val_batches = 500\n    val_freq = 500\n\n    # Create DataLoader\n    cv_dataloader = DataLoader(cv_dataset, batch_size=batch_size, shuffle=True, num_workers=0)\n    \n    # List to hold attributes of best models\n    best_model_list = []\n    \n    # Train loop\n    for i_cv in range(n_folds):\n#     for i_cv in [0]:\n        print('CV: ',i_cv)\n\n        i_train_loss = []\n        i_val_map = []\n        i_val_class_map = []    \n        epoch_loss_l = []\n        val_loss = 0\n        best_val_map = 0\n\n        print('Initializing Model')\n        cv_dataset.set_cv_idx(cv_idx=[i_cv], val=True, protocols=protocols, cv_col=cv_col)\n        print('Num Val Batches: ' + str(len(cv_dataloader)))\n        cv_dataset.set_cv_idx(cv_idx=[i_cv], val=False, protocols=protocols, cv_col=cv_col)\n        num_train_batches = str(len(cv_dataloader))\n        print( 'Num Train Batches: ' + num_train_batches)\n        \n        fog_model = cnn(input_dim, conv_layers, num_classes, kernel_size, stride, padding, incl_pool, dropout_p, pool_k, pool_s)\n        fog_model.to(DEVICE)\n        if use_class_weights:\n            class_weights = torch.tensor(cv_dataset.lab_freq).float()\n        else:\n            class_weights = torch.tensor([1,1,1,1]).float()\n\n        optimizer = optim.Adam(fog_model.parameters(), lr=learning_rate, weight_decay=1e-5)\n#         loss_fn = nn.CrossEntropyLoss(weight=class_weights.to(DEVICE))\n#         loss_fn = FocalLoss(alpha=class_weights.to(DEVICE), gamma=4)\n#         loss_fn = nn.MultiMarginLoss(weight=class_weights.to(DEVICE))\n#         loss_fn = FocalLoss(gamma=5)\n        if use_focal_loss:\n            loss_fn = FocalLoss(gamma=5)\n        else:\n            loss_fn = nn.CrossEntropyLoss(weight=class_weights.to(DEVICE))\n        \n        print('Training Model')\n        start_time = time.time()\n        for epoch in range(num_epochs):\n            epoch_loss = 0\n            epoch_map = 0 \n            n_batches_done = 1\n            cv_dataset.set_cv_idx(cv_idx=[i_cv], val=False, protocols=protocols, cv_col=cv_col)\n            \n            for batch in cv_dataloader:\n                x, y = batch\n                x = x.to(DEVICE)\n                y = y.to(DEVICE)\n                y_pred = fog_model(x)\n                with warnings.catch_warnings():\n                    warnings.simplefilter(\"ignore\")\n                    loss = loss_fn(y_pred.to(DEVICE), y.view(-1))\n                    optimizer.zero_grad(set_to_none=True)\n                    loss.backward()\n                    optimizer.step()\n\n                # calculate loss\n                epoch_loss += loss.item()\n                epoch_loss_l.append(loss.item())\n                avg_epoch_loss = epoch_loss*batch_size/(n_batches_done)\n\n                # calculate mAP\n                y_one_hot = torch.squeeze(torch.nn.functional.one_hot(torch.squeeze(y,1)),1)\n                batch_mAP = mAP(torch.t(y_pred).cpu().detach().numpy(),torch.t(y_one_hot).cpu().detach().numpy(),3)\n\n                if n_batches_done % val_freq == 0 or n_batches_done == len(cv_dataloader):#n_batches_done == len(cv_dataloader):#\n                    print(\"{} seconds; Epoch: {}; Batch: {}/{}; Loss: {}; Batch mAP: {} {}\".format(round((time.time() - start_time)), epoch, n_batches_done, num_train_batches, round(avg_epoch_loss,4), np.round(np.nanmean(batch_mAP),5), np.round(batch_mAP,5)))\n                    #print('Batch ' + str(n_batches_done) + f' {(time.time() - start_time):.1f} sec' + '; Loss: ' + str(round(avg_epoch_loss,4)) + '; Batch mAP: ' + str(np.round(np.mean(batch_mAP),5)) + str(np.round(batch_mAP,5)))\n                \n                    val_loss = 0\n                    actuals =[]\n                    preds=[]\n                    fog_model.eval()\n                    with torch.no_grad():\n                        cv_dataset.set_cv_idx(cv_idx=[i_cv], val=True, protocols=protocols, cv_col=cv_col)\n#                         print('Num Val Batches: ' + str(len(cv_dataloader)))\n                        n_val_batches = 0\n                        for batch in cv_dataloader:\n                            x_val, y_val = batch\n                            x_val = x_val.to(DEVICE)\n                            y_val = y_val.to(DEVICE)\n                            y_pred_val = fog_model(x_val)\n                            actuals.append(torch.t(torch.nn.functional.one_hot(y_val.squeeze())).cpu().detach().numpy())\n                            preds.append(y_pred_val.squeeze().cpu().detach().numpy())\n                            val_loss += loss_fn(y_pred_val.to(DEVICE), y_val.view(-1))\n                            \n                            n_val_batches += 1\n                            if n_val_batches > max_val_batches:\n                                break\n                    cv_dataset.set_cv_idx(cv_idx=[i_cv], val=False, protocols=protocols, cv_col=cv_col)\n                    val_loss = val_loss/(n_val_batches*batch_size)\n                    \n                    fog_model.train()\n\n                    actuals = np.concatenate(actuals,axis=1)\n                    preds = np.transpose(np.concatenate(preds,axis=0))\n                    val_map = mAP(preds,actuals,3)\n                    i_train_loss.append(np.nanmean(epoch_loss_l))\n                    i_val_map.append(np.round(np.nanmean(val_map),decimals=5))\n                    i_val_class_map.append(val_map)\n\n#                     print(\"Epoch: {}/{},  train_loss: {},  train_mAP: {},  val_mAP: {}\".format(epoch+1,num_epochs,avg_epoch_loss,None,np.round(np.nanmean(val_loss),decimals=5)))\n                    print(\"Validation mAP: {} {};\".format(np.round(np.nanmean(val_map), decimals=5), np.round(val_map, decimals=5)))\n\n                    if np.nanmean(val_map) > best_val_map:\n                        torch.save(fog_model.state_dict(), 'best_{}_cv{}.pt'.format(protocols[0],i_cv))\n                        best_val_map = np.nanmean(val_map)\n                        best_mod_pd = pd.DataFrame({'CV': i_cv, 'Epoch': epoch, 'Batch': n_batches_done, 'Val mAP': best_val_map},index=[i_cv])\n                \n                n_batches_done += 1\n                if n_batches_done > max_batches_per_epoch:\n                    break\n        \n        best_model_list.append(best_mod_pd)\n    best_mods_pd =  pd.concat(best_model_list)\n    print(best_mods_pd)\n\n    return fog_model","metadata":{"execution":{"iopub.status.busy":"2023-07-06T17:44:41.418641Z","iopub.execute_input":"2023-07-06T17:44:41.419649Z","iopub.status.idle":"2023-07-06T17:44:41.457525Z","shell.execute_reply.started":"2023-07-06T17:44:41.419608Z","shell.execute_reply":"2023-07-06T17:44:41.456408Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# TDCSF TRAINING LOOP \n\n# initialize dataloader\n# protocols = ['tdcsfog','defog']\n# cv_col = ['combo_cv']\nprotocols = ['tdcsfog']\ncv_col = ['tdcsf_cv']\ncv_dataset = AccDataset(files_df=tr_files_df, patients_df=patients_df, valid_protocols=protocols)\ncv_dataset.set_cv_idx(cv_idx=[0], val=False, protocols=protocols, cv_col=cv_col)\nuse_class_weights = True\nuse_focal_loss = False\n\n# Train model\nwith np.errstate(divide='ignore',invalid='ignore'):\n    fog_model = CV_TRAIN()","metadata":{"execution":{"iopub.status.busy":"2023-07-06T17:44:41.458898Z","iopub.execute_input":"2023-07-06T17:44:41.460295Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# DeFOG TRAINING LOOP\n\n# initialize dataloader\ndel cv_dataset\nprotocols = ['defog']\ncv_col = ['defog_cv']\ncv_dataset = AccDataset(files_df=tr_files_df, patients_df=patients_df, valid_protocols=protocols)\ncv_dataset.set_cv_idx(cv_idx=[0], val=False, protocols=protocols, cv_col=cv_col)\nuse_class_weights = False\nuse_focal_loss = True\n\n# # Train model\nwith np.errstate(divide='ignore',invalid='ignore'):\n    fog_model = CV_TRAIN()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# fig, ax = plt.subplots(n_folds)\n# for c in range(n_folds):\n#     ax[c].plot(np.array(cv_val_map[c]), label=str('cv fold: {}'.format(c)), alpha=0.5)\n#     ax[c].set_title(str('cv fold: {}'.format(c)))\n# plt.xlabel('Epochs')\n# plt.ylabel('Val mAP')\n# #plt.title('Dataset: {}, lr: {}, batch size: {}, class weights: {}, dropout: {}, pool_k: {}, conv_layers: {}'.format(protocols,learning_rate,batch_size,cv_class_w,dropout_p,pool_k,conv_layers))\n# plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# MAKE FINAL PREDICTIONS\nfinal_preds = pd.DataFrame(columns=['Id','StartHesitation','Turn','Walking'])\n\nfor i_prot in ['defog','tdcsfog']:\n    test_dataset = AccDataset(files_df=test_files_df, patients_df=patients_df, valid_protocols=[i_prot], tr=False)\n    test_dataloader = DataLoader(test_dataset, batch_size=1024, shuffle=True)\n    for batch in test_dataloader:\n        x, _, k = batch\n        y_pred_ens = np.zeros([np.shape(x)[0],4])\n        for i_model in range(n_folds):\n            fog_model.load_state_dict(torch.load('best_{}_cv{}.pt'.format(i_prot,i_model)))\n#             fog_model.load_state_dict(torch.load('best_{}_cv{}.pt'.format('tdcsfog',i_model)))\n            with torch.no_grad():\n                fog_model.eval()\n                # Make batch predictions\n                x = x.to(DEVICE)\n                y_pred = fog_model(x)\n                y_pred = torch.squeeze(y_pred,0).cpu().detach().numpy()\n                y_pred_ens += y_pred/n_folds\n        i_pd = pd.DataFrame({'Id':np.squeeze(k,0),\n                         'StartHesitation':y_pred_ens[:,0],\n                        'Turn':y_pred_ens[:,1],\n                        'Walking':y_pred_ens[:,2]\n                        })\n        final_preds = pd.concat([final_preds, i_pd], axis=0)\n    \n# write predictions to file\nfinal_preds.to_csv('submission.csv', index=False, float_format='%.5f')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# validate preds and submissions\n# sample_sub = pd.read_csv(\"/kaggle/input/tlvmc-parkinsons-freezing-gait-prediction/sample_submission.csv\")\n# shared_ids = sum(sample_sub['Id'].isin(final_preds['Id']))","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}