{"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":"markdown","source":"# Dataset","metadata":{}},{"cell_type":"code","source":"import glob\nimport pandas as pd\nimport numpy as np\nfrom torch.utils.data import Dataset\nfrom pathlib import Path\nimport pywt\n\nBASE_FEATURE_COUNT = 5\nWINDOW_LENGTH = 10240\nSTRIDE_DENOMINATOR = 8\nSTRIDE_LENGTH = WINDOW_LENGTH // STRIDE_DENOMINATOR\n\ndef GetBandPowers(signals, sample_rate_hz, freq_bands):\n    # Clculate time vector.\n    sample_count = len(signals)\n    t = np.arange(0, sample_count) / sample_rate_hz\n\n    # Set up wavelet parameters\n    wavelet = 'mexh'\n    \n    freqs = np.logspace(np.log10(0.5), np.log10(freq_bands[-1][1]), num=100)  # Frequency range to analyze\n    n_cycles = freqs / 2.  # Number of cycles for each frequency\n\n    # Loop through each accelerometer and compute wavelet power\n    power = np.zeros((signals.shape[1], len(freqs), signals.shape[0]), dtype=np.float32)\n    for i in range(signals.shape[1]):\n        cwtmatr, _ = pywt.cwt(signals[:, i], n_cycles, wavelet, sampling_period=1/sample_rate_hz)\n        power[i] = (np.abs(cwtmatr)**2) / sample_rate_hz\n\n    # Sum power across frequency bands\n    band_power = np.zeros((signals.shape[1], len(freq_bands), signals.shape[0]), dtype=np.float32)\n    for i in range(signals.shape[1]):\n        for j, band in enumerate(freq_bands):\n            freq_idx = np.logical_and(freqs >= band[0], freqs < band[1])\n            band_power[i, j] = np.sum(power[i, freq_idx, :], axis=0)\n\n    return band_power\n\nclass ParkinsonsDataset(Dataset):\n    def __init__(self, frequency_bands, max_sequence_length = 10240, stride_denominator = 1):\n        self.MaxLength = max_sequence_length\n        self.FeaturesSequences = []\n        self.TimeSequences = []\n        self.FileIds = []\n\n        self.FileIdsToMaxTimestamps = {}\n        \n        self.FrequencyBands = frequency_bands\n        \n        # Has accleration units in g (9.806 m/s^2)\n        filepaths = glob.glob('/kaggle/input/tlvmc-parkinsons-freezing-gait-prediction/test/*/*.csv')\n        for filepath in filepaths:\n            df = pd.read_csv(filepath)\n            \n            #if 'tdcsfog' in filepath:\n            #if 'tdcsfog' not in filepath:\n            #    G = 9.806\n            #    df[['AccV', 'AccML', 'AccAP']] /= G\n                \n            df['NormalizedTime'] = df['Time'] / df['Time'].max()\n            df['SinNormalizedTime'] = np.sin(df['NormalizedTime'] * np.pi)\n            signals = df[['AccV', 'AccML', 'AccAP']].to_numpy()\n            sample_rate_hz = 100 if 'defog' in filepath else 128\n            if self.FrequencyBands is not None:\n                band_powers = GetBandPowers(signals, sample_rate_hz, self.FrequencyBands) # Has shape (accelerometer axis count, band count, sequence length)\n                all_features = np.concatenate((\n                    df[['NormalizedTime', 'SinNormalizedTime']].to_numpy(),\n                    band_powers.reshape(-1, band_powers.shape[-1]).T,\n                    signals\n                ), axis = 1)\n            else:\n                all_features = np.concatenate((\n                    df[['NormalizedTime', 'SinNormalizedTime']].to_numpy(),\n                    signals\n                ), axis = 1)\n            \n            times = df['Time'].to_numpy()\n\n            raw_length = len(all_features)\n            file_id = Path(filepath).stem\n            for start_index in range(0, raw_length, self.MaxLength // stride_denominator):\n                self.FeaturesSequences.append(all_features[start_index:start_index+self.MaxLength])\n                self.TimeSequences.append(times[start_index:start_index+self.MaxLength])\n                self.FileIds.append(file_id)\n                \n            self.FileIdsToMaxTimestamps[file_id] = max(times)\n\n    def __len__(self):\n        return len(self.FeaturesSequences)\n\n    def __getitem__(self, index):\n        sequence_length = len(self.FeaturesSequences[index])\n\n        feature_count = BASE_FEATURE_COUNT + 3*len(self.FrequencyBands) if self.FrequencyBands is not None else BASE_FEATURE_COUNT\n        padded_features = np.zeros((self.MaxLength, feature_count))\n        padded_features[:sequence_length,:] = self.FeaturesSequences[index]\n        \n        padded_times = np.ones((self.MaxLength,), dtype = int) * -1\n        padded_times[:sequence_length] = self.TimeSequences[index]\n\n        return padded_features, padded_times, self.FileIds[index]\n    \n#dataset = ParkinsonsDataset([(0.5, 3), (3, 10)], WINDOW_LENGTH, STRIDE_DENOMINATOR)\ndataset = ParkinsonsDataset(None, WINDOW_LENGTH, STRIDE_DENOMINATOR)\nprint(len(dataset))\nprint(dataset[0])\nprint(np.shape(dataset[0][0]))\n\nprint(dataset[-1])","metadata":{"execution":{"iopub.status.busy":"2023-06-08T04:22:43.051664Z","iopub.execute_input":"2023-06-08T04:22:43.051987Z","iopub.status.idle":"2023-06-08T04:22:46.793999Z","shell.execute_reply.started":"2023-06-08T04:22:43.051955Z","shell.execute_reply":"2023-06-08T04:22:46.792469Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model architecture","metadata":{}},{"cell_type":"code","source":"import torch\nimport torch.nn as nn\nimport torch.nn.functional as F\n\nimport torch\nimport torch.nn as nn\n\nclass Conv1dBlockSE(nn.Module):\n    def __init__(self, channels, kernel_size=3, stride=1, padding=1, reduction=16, dropout=0):\n        super(Conv1dBlockSE, self).__init__()\n        self.conv = nn.Conv1d(channels, channels, kernel_size, stride, padding, bias=False, groups=1)\n        self.bn = nn.BatchNorm1d(channels)\n        self.relu = nn.ReLU(inplace=True)\n        self.dropout = nn.Dropout1d(p=dropout)\n\n        self.pool = nn.AdaptiveAvgPool1d(1)\n        self.fc1 = nn.Linear(channels, channels // reduction)\n        self.fc2 = nn.Linear(channels // reduction, channels)\n        self.sigmoid = nn.Sigmoid()\n\n    def forward(self, x):\n        residual = x\n        x = self.conv(x)\n        x = self.bn(x)\n        x = self.relu(x)\n        x = self.dropout(x)\n\n        # Squeeze\n        se = self.pool(x)\n        se = se.view(se.size(0), -1)\n        se = self.fc1(se)\n        se = self.relu(se)\n        se = self.fc2(se)\n        se = self.sigmoid(se)\n\n        # Excitation\n        x = x * se.unsqueeze(2)\n\n        # Add residual and return\n        x += residual\n\n        return x\n\nclass NonResidualConvSE(nn.Module):\n    def __init__(self, in_channels, out_channels, kernel_size=3, stride=1, padding=1, reduction=16, dropout=0):\n        super(NonResidualConvSE, self).__init__()\n\n        self.Conv = nn.Sequential(\n            nn.Conv1d(in_channels, out_channels, kernel_size, stride, padding, bias=False),\n            nn.BatchNorm1d(out_channels),\n            nn.ReLU(),\n            nn.Dropout(p = dropout),\n        )\n\n        self.Pooling = nn.AdaptiveAvgPool1d(1)\n        self.SqueezeExcitationWeightGenerator = nn.Sequential(\n            nn.Linear(out_channels, out_channels // reduction),\n            nn.ReLU(),\n            nn.Linear(out_channels // reduction, out_channels),\n            nn.Sigmoid()\n        )\n\n    def forward(self, x):\n        x = self.Conv(x)\n        pooled_x = self.Pooling(x)\n        excitation_weights = self.SqueezeExcitationWeightGenerator(pooled_x.view(pooled_x.size(0), -1))\n        x = x * excitation_weights.unsqueeze(2)\n\n        return x\n\nclass Conv1dBlockResidual(nn.Module):\n    def __init__(self, channels, kernel_size=3, stride=1, padding=1, dropout=0):\n        super(Conv1dBlockResidual, self).__init__()\n\n        self.Layers = nn.Sequential(\n            nn.Conv1d(channels, channels, kernel_size, stride, padding, bias=False),\n            nn.BatchNorm1d(channels),\n            nn.ReLU(),\n            nn.Dropout1d(p = dropout),\n            nn.Conv1d(channels, channels, kernel_size, stride, padding, bias=False),\n            nn.BatchNorm1d(channels),\n            nn.ReLU(),\n            nn.Dropout1d(p = dropout)\n        )\n\n    def forward(self, x):\n        x = x + self.Layers(x)\n\n        return x\n\nclass Conv1dBlock(nn.Module):\n    def __init__(self, in_channels, out_channels, kernel_size=3, stride=1, padding=1, dropout=0):\n        super(Conv1dBlock, self).__init__()\n        self.conv = nn.Conv1d(in_channels, out_channels, kernel_size, stride, padding, bias=False)\n        self.bn = nn.BatchNorm1d(out_channels)\n        self.relu = nn.ReLU(inplace=True)\n        self.dropout = nn.Dropout1d(p=dropout)\n\n    def forward(self, x):\n        x = self.conv(x)\n        x = self.bn(x)\n        x = self.relu(x)\n        x = self.dropout(x)\n        return x\n\nclass Conv1dBlockPreprocessedSE(nn.Module):\n    def __init__(self, in_channels, out_channels, reduction=16, use_second_se=False, preprocessor_dropout=0, se_dropout=0):\n        super(Conv1dBlockPreprocessedSE, self).__init__()\n\n        KERNEL_SIZE=3\n        STRIDE=1\n        PADDING=1\n        self.Preprocessor = Conv1dBlock(in_channels, out_channels, KERNEL_SIZE, STRIDE, PADDING, preprocessor_dropout)\n        self.SqueezeAndExcitation1 = Conv1dBlockSE(out_channels, KERNEL_SIZE, STRIDE, PADDING, reduction, se_dropout)\n        self.SqueezeAndExcitation2 = Conv1dBlockSE(out_channels, KERNEL_SIZE, STRIDE, PADDING, reduction, se_dropout)\n\n        self.UseSecondSe = use_second_se\n\n    def forward(self, x):\n        x = self.Preprocessor(x)\n        x = self.SqueezeAndExcitation1(x)\n        if self.UseSecondSe:\n            x = self.SqueezeAndExcitation2(x)\n        return x\n\nclass SequenceClassifier(nn.Module):\n    def __init__(\n            self, \n            in_channels=3, \n            out_channels=4, \n            model_width_coef=32, \n            reduction=16, \n            use_second_se=False, \n            preprocessor_dropout=0, \n            se_dropout=0,\n            initial_dropout=0,\n            center_dropout=0):\n        super(SequenceClassifier, self).__init__()\n\n        features = model_width_coef\n\n        self.encoder1 = nn.Sequential(\n            NonResidualConvSE(in_channels, features, reduction = reduction//2, dropout=initial_dropout),\n            NonResidualConvSE(features, features, reduction = reduction//2, dropout=initial_dropout),\n            Conv1dBlockPreprocessedSE(features, features, reduction, use_second_se, preprocessor_dropout, se_dropout)\n        )\n\n        # Encoder part\n        # self.encoder1 = Conv1dBlockPreprocessedSE(in_channels, features, reduction, use_second_se, preprocessor_dropout, se_dropout)\n        self.pool1 = nn.MaxPool1d(kernel_size=2, stride=2)\n        self.encoder2 = Conv1dBlockPreprocessedSE(features, features*2, reduction, use_second_se, preprocessor_dropout, se_dropout)\n        self.pool2 = nn.MaxPool1d(kernel_size=2, stride=2)\n        self.encoder3 = Conv1dBlockPreprocessedSE(features*2, features*4, reduction, use_second_se, preprocessor_dropout, se_dropout)\n        self.pool3 = nn.MaxPool1d(kernel_size=2, stride=2)\n        self.encoder4 = Conv1dBlockPreprocessedSE(features*4, features*8, reduction, use_second_se, preprocessor_dropout, se_dropout)\n        self.pool4 = nn.MaxPool1d(kernel_size=2, stride=2)\n        self.encoder5 = Conv1dBlockPreprocessedSE(features*8, features*16, reduction, use_second_se, preprocessor_dropout, se_dropout)\n        self.pool5 = nn.MaxPool1d(kernel_size=2, stride=2)\n\n        # Bottleneck part\n        self.bottleneck = Conv1dBlock(features*16, features*32, dropout=center_dropout)\n\n        # Decoder part\n        self.upconv5 = nn.ConvTranspose1d(features*32, features*16, kernel_size=2, stride=2)\n        self.decoder5 = Conv1dBlockPreprocessedSE(features*32, features*16, reduction, use_second_se, preprocessor_dropout, se_dropout)\n        self.upconv4 = nn.ConvTranspose1d(features*16, features*8, kernel_size=2, stride=2)\n        self.decoder4 = Conv1dBlockPreprocessedSE(features*16, features*8, reduction, use_second_se, preprocessor_dropout, se_dropout)\n        self.upconv3 = nn.ConvTranspose1d(features*8, features*4, kernel_size=2, stride=2)\n        self.decoder3 = Conv1dBlockPreprocessedSE(features*8, features*4, reduction, use_second_se, preprocessor_dropout, se_dropout)\n        self.upconv2 = nn.ConvTranspose1d(features*4, features*2, kernel_size=2, stride=2)\n        self.decoder2 = Conv1dBlockPreprocessedSE(features*4, features*2, reduction, use_second_se, preprocessor_dropout, se_dropout)\n        self.upconv1 = nn.ConvTranspose1d(features*2, features, kernel_size=2, stride=2)\n        self.decoder1 = Conv1dBlockPreprocessedSE(features*2, features, reduction, use_second_se, preprocessor_dropout, se_dropout)\n\n        # Output layer\n        self.out_conv = nn.Conv1d(features, out_channels, kernel_size=1)\n\n    def forward(self, x):\n        x = x.permute(0, 2, 1)\n\n        # Encoder part\n        enc1 = self.encoder1(x)\n        enc2 = self.encoder2(self.pool1(enc1))\n        enc3 = self.encoder3(self.pool2(enc2))\n        enc4 = self.encoder4(self.pool3(enc3))\n        enc5 = self.encoder5(self.pool4(enc4))\n\n        # enc1 = self.Enc1Resnet(enc1)\n        # enc2 = self.Enc2Resnet(enc2)\n\n\n        # Bottleneck part\n        bottleneck = self.bottleneck(self.pool5(enc5))\n\n        # Decoder part\n        dec5 = self.upconv5(bottleneck)\n        dec5 = self.decoder5(torch.cat((dec5, enc5), dim=1))\n        dec4 = self.upconv4(dec5)\n        dec4 = self.decoder4(torch.cat((dec4, enc4), dim=1))\n        dec3 = self.upconv3(dec4)\n        dec3 = self.decoder3(torch.cat((dec3, enc3), dim=1))\n        dec2 = self.upconv2(dec3)\n        dec2 = self.decoder2(torch.cat((dec2, enc2), dim=1))\n        dec1 = self.upconv1(dec2)\n        dec1 = self.decoder1(torch.cat((dec1, enc1), dim = 1))\n\n        results = self.out_conv(dec1)\n        results = results.permute(0, 2, 1)\n\n        return results","metadata":{"execution":{"iopub.status.busy":"2023-06-08T04:22:46.797937Z","iopub.execute_input":"2023-06-08T04:22:46.798373Z","iopub.status.idle":"2023-06-08T04:22:46.864690Z","shell.execute_reply.started":"2023-06-08T04:22:46.798330Z","shell.execute_reply":"2023-06-08T04:22:46.862631Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Inference","metadata":{}},{"cell_type":"code","source":"import torch\nfrom torch.utils.data import DataLoader\nimport sys\nfrom math import ceil\nimport glob\n\nBATCH_SIZE = 64\ndataloader = DataLoader(dataset, batch_size=BATCH_SIZE, num_workers=1)\n\nmodel_paths =[\n    #'/kaggle/input/parkinsonsfrequencydomainmodels/Limited/split0_346_459_48_479.pth'\n    #'/kaggle/input/parkinsonsfrequencydomainmodels/Limited/split1_410_738_146_374.pth'\n    #'/kaggle/input/parkinsonsfrequencydomainmodels/Limited/split10_385_570_133_417.pth'\n    \n    '/kaggle/input/parkinsonsfrequencydomainmodels/None/split10_504_543_164_386.pth', # Strong! 0.411\n    \n    #'/kaggle/input/cherrypickedparkinsons/CherryPicked/split10_494_544_150_406.pth'\n    '/kaggle/input/cherrypickedparkinsons/CherryPicked/split10_551_573_175_376.pth' # Strong ish. 0.38\n    #'/kaggle/input/cherrypickedparkinsons/CherryPicked/split10_547_564_198_358.pth'\n    \n]\n#model_paths = glob.glob('/kaggle/input/parkinsonsfrequencydomainmodels/Limited/*.pth')\n\n#model_paths = glob.glob('/kaggle/input/cherrypickedparkinsons/CherryPicked/*.pth')\n\nmodels = []\nfor model_path in model_paths:\n    model = SequenceClassifier(in_channels = BASE_FEATURE_COUNT, model_width_coef=32).cuda()\n    state_dict = torch.load(model_path, map_location=torch.device('cuda:0'))\n    for key in list(state_dict.keys()):\n        state_dict[key.replace('_orig_mod.', '')] = state_dict.pop(key)\n\n    model.load_state_dict(state_dict)\n    model.eval()\n    models.append(model)\n\n# Should have entries like\n# FileId : {\n#     timestamp_0 : [\n#         [class 1 score 1, class 1 score 2, ...],\n#         [class 2 score 1, class 2 score 2, ...],\n#         [class 3 score 1, class 3 score 2, ...],\n#     ],\n#     timestamp_1 : [\n#         [class 1 score 1, class 1 score 2, ...],\n#         [class 2 score 1, class 2 score 2, ...],\n#         [class 3 score 1, class 3 score 2, ...],\n#     ],\n#     ...\n# }\n# \nfile_ids_to_timestamps_to_scores = {}\nwith open('submission.csv', 'w') as submission_file:\n    submission_file.write('Id,StartHesitation,Turn,Walking\\n')\n    \n    for padded_features, padded_times, padded_ids in dataloader:\n        padded_features = padded_features.float().cuda()\n        # Has shape (model count, series count, padded_series_length, event_class_count + 1)\n        all_predictions = []\n        with torch.no_grad():\n            for model in models:\n                predictions = model(padded_features)\n                predictions = torch.nn.functional.softmax(predictions, dim=2).cpu().numpy()\n                all_predictions.append(predictions)\n                \n                padded_features[:,:,-2] *= -1\n                predictions = model(padded_features)\n                predictions = torch.nn.functional.softmax(predictions, dim=2).cpu().numpy()\n                all_predictions.append(predictions)\n\n        averaged_predictions = np.mean(all_predictions, axis = 0)\n\n        for series_index, file_id in enumerate(padded_ids):\n            if file_id not in file_ids_to_timestamps_to_scores.keys():\n                file_ids_to_timestamps_to_scores[file_id] = [\n                    [[], [] ,[]]\n                    for timestamp in range(dataset.FileIdsToMaxTimestamps[file_id] + 1)\n                ]\n\n            series_predictions = averaged_predictions[series_index]\n            series_timestamps = padded_times[series_index].numpy()\n            file_timestamps_to_scores = file_ids_to_timestamps_to_scores[file_id]\n            max_file_timestamp = dataset.FileIdsToMaxTimestamps[file_id]\n            for timestep_index, timestep_predictions in enumerate(series_predictions):\n                timestamp = int(series_timestamps[timestep_index])\n                if timestamp < 0:\n                    break\n\n                timestamp_scores = file_timestamps_to_scores[timestamp]\n                timestamp_scores[0].append(timestep_predictions[1])\n                timestamp_scores[1].append(timestep_predictions[2])\n                timestamp_scores[2].append(timestep_predictions[3])\n\n                max_expected_endpoint_samples = min(ceil(dataset.FileIdsToMaxTimestamps[file_id] / STRIDE_LENGTH), STRIDE_DENOMINATOR)\n                if (timestamp == max_file_timestamp) and (len(timestamp_scores[0]) == max_expected_endpoint_samples):\n                    print('Dumping:', file_id, timestep_index, timestamp, max_file_timestamp)\n                    for timestamp_to_store in range(dataset.FileIdsToMaxTimestamps[file_id] + 1):\n                        sample_id = f'{file_id}_{timestamp_to_store}'\n                        scores = file_timestamps_to_scores[timestamp_to_store]\n                        \n                        weights = np.ones(len(scores[0]))\n                        #weights[0] = .5\n                        weights[-1] = .5\n                        \n                        #weight_count = len(scores[0])\n                        #weights = np.sin(np.linspace(0, np.pi, weight_count)) + 0.25\n                        \n                        avg_scores = np.average(scores, axis = 1, weights = weights)\n                        submission_file.write(f'{sample_id},{avg_scores[0]},{avg_scores[1]},{avg_scores[2]}\\n')\n                    \n                    file_ids_to_timestamps_to_scores.pop(file_id)\n\n\n        print(sys.getsizeof(file_ids_to_timestamps_to_scores))\n","metadata":{"execution":{"iopub.status.busy":"2023-06-08T04:22:46.866346Z","iopub.execute_input":"2023-06-08T04:22:46.867922Z","iopub.status.idle":"2023-06-08T04:23:22.575002Z","shell.execute_reply.started":"2023-06-08T04:22:46.867880Z","shell.execute_reply":"2023-06-08T04:23:22.572822Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Sanity check results","metadata":{}},{"cell_type":"code","source":"import pandas as pd\n\npredictions_df = pd.read_csv('submission.csv')\npredictions_df.head()","metadata":{"execution":{"iopub.status.busy":"2023-06-08T04:23:22.578580Z","iopub.execute_input":"2023-06-08T04:23:22.580571Z","iopub.status.idle":"2023-06-08T04:23:22.869692Z","shell.execute_reply.started":"2023-06-08T04:23:22.580525Z","shell.execute_reply":"2023-06-08T04:23:22.868637Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"predictions_df.tail()","metadata":{"execution":{"iopub.status.busy":"2023-06-08T04:23:22.871451Z","iopub.execute_input":"2023-06-08T04:23:22.871837Z","iopub.status.idle":"2023-06-08T04:23:22.886741Z","shell.execute_reply.started":"2023-06-08T04:23:22.871789Z","shell.execute_reply":"2023-06-08T04:23:22.884946Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"predictions_df.describe()","metadata":{"execution":{"iopub.status.busy":"2023-06-08T04:23:22.888760Z","iopub.execute_input":"2023-06-08T04:23:22.889555Z","iopub.status.idle":"2023-06-08T04:23:22.935329Z","shell.execute_reply.started":"2023-06-08T04:23:22.889516Z","shell.execute_reply":"2023-06-08T04:23:22.934348Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"file_ids_to_timestamps_to_scores","metadata":{"execution":{"iopub.status.busy":"2023-06-08T04:23:22.937226Z","iopub.execute_input":"2023-06-08T04:23:22.938373Z","iopub.status.idle":"2023-06-08T04:23:22.945732Z","shell.execute_reply.started":"2023-06-08T04:23:22.938330Z","shell.execute_reply":"2023-06-08T04:23:22.944447Z"},"trusted":true},"execution_count":null,"outputs":[]}]}