{"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":"# *KAGGLE CHALLENGE: LANL Earthquake Prediction*\n\nUn projet de Matthieu Dagommer, Paul Boulgakoff, Godefroy Bichon, Germain L'Hostis\n\nVersions utilisées:\n\nPython: 3.10.4\nTorch: 1.11","metadata":{"id":"JqdSSxEJLtds"}},{"cell_type":"code","source":"import numpy as np\nimport math\nimport matplotlib\nimport matplotlib.pyplot as plt \nimport time\n\nimport torch as th\nth.cuda.empty_cache()\n\nimport torch.autograd as autograd\nfrom torch.autograd import Variable\nimport torch.nn.functional as F\nimport torch.nn as nn\n\n#from torchsummary import summary \n\nimport gzip\nimport pickle\n\nimport pandas as pd\nimport scipy\nimport scipy.stats as stats\nfrom scipy.stats import kurtosis, skew\nimport csv\nimport os\nimport pickle\nimport gc\n\nfrom sklearn.preprocessing import StandardScaler, MinMaxScaler","metadata":{"execution":{"iopub.status.busy":"2022-05-29T21:55:07.561656Z","iopub.execute_input":"2022-05-29T21:55:07.562024Z","iopub.status.idle":"2022-05-29T21:55:07.570234Z","shell.execute_reply.started":"2022-05-29T21:55:07.561976Z","shell.execute_reply":"2022-05-29T21:55:07.569166Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"### Ensuring CPU is connected\n\nprint(th.cuda.is_available())\nprint(th.cuda.get_device_name())\ndevice = th.device('cuda' if th.cuda.is_available() else 'cpu')","metadata":{"execution":{"iopub.status.busy":"2022-05-29T21:55:08.018761Z","iopub.execute_input":"2022-05-29T21:55:08.019125Z","iopub.status.idle":"2022-05-29T21:55:08.024306Z","shell.execute_reply.started":"2022-05-29T21:55:08.019095Z","shell.execute_reply":"2022-05-29T21:55:08.023274Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"os.mkdir('./Models')","metadata":{"execution":{"iopub.status.busy":"2022-05-29T21:55:16.312163Z","iopub.execute_input":"2022-05-29T21:55:16.313214Z","iopub.status.idle":"2022-05-29T21:55:16.330959Z","shell.execute_reply.started":"2022-05-29T21:55:16.313175Z","shell.execute_reply":"2022-05-29T21:55:16.32987Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# *General Hyperparameters*","metadata":{}},{"cell_type":"code","source":"### Setting General Hyperparameters \n\nbatch_size = 100 # number of patches per batch\nvalid_rate = 0.1 # fraction of data dedicated to validation\noverlap_rate = 0.2 # overlap\n\n# Parameters\n\nnrows = 50_000_000\nseed = 1\npatch_size = 150_000","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"th.manual_seed(seed)\nnp.random.seed(seed)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# *Data Preparation*\n\n","metadata":{}},{"cell_type":"code","source":"### Loading Data\n\n#rootpath = os.getcwd() + \"/\"\nrootpath = \"../input/LANL-Earthquake-Prediction/\"\n\ntrain_data = pd.read_csv(rootpath + \"train.csv\", usecols = ['acoustic_data', 'time_to_failure'], \\\n                         dtype = {'acoustic_data': np.int16, 'time_to_failure': np.float64}, nrows = nrows)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X = th.squeeze(th.tensor(train_data['acoustic_data'].values, dtype = th.int16))\nY = th.squeeze(th.tensor(train_data['time_to_failure'].values, dtype = th.float64))\n\ndel train_data\ngc.collect()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"### Function to create patched sequences with some overlap out of the training data\n\n# Overlap is expressed as a fraction of patch size\ndef patching(patch_size, X, Y, overlap_rate):\n    \n    overlap = int(patch_size*overlap_rate)\n    L = X.shape[0] # total length of the training acoustic signal\n    n_patch = int(np.floor((L-patch_size)/(patch_size-overlap)+1))\n    ids_no_seism = []\n\n    X_patch = th.zeros(n_patch,patch_size)\n    Y_patch = th.zeros(n_patch)\n    \n    for i in range (n_patch):\n        X_patch[i,:] = X[i*(patch_size-overlap):(i+1)*patch_size-i*overlap] \n        Y_patch[i] = Y[(i+1)*patch_size - i*overlap] \n    \n        # Removing patches with no seism\n        if th.min(Y[i*(patch_size-overlap):(i+1)*patch_size - i*overlap]) > 0.001:\n            ids_no_seism.append(i)\n    \n    X_patch = X_patch[ids_no_seism]\n    Y_patch = Y_patch[ids_no_seism]\n    \n    return(n_patch, X_patch, Y_patch)","metadata":{"id":"nn_5lrssRSpJ","outputId":"e0eb114e-3c57-458a-bec2-53709bbd85d8"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"### Normalizing Data\n\nss = StandardScaler()\nmm = MinMaxScaler()\n    \nX = th.squeeze(th.from_numpy(ss.fit_transform(np.array(X).reshape(-1, 1))))\nY = th.squeeze(th.from_numpy(mm.fit_transform(np.array(Y).reshape(-1, 1))))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"### Distribution of data in patches\n\n_, X_patch, Y_patch = patching(patch_size, X, Y, overlap_rate = overlap_rate)\n\ndel X; del Y\ngc.collect()\n\nprint(X_patch.shape)\nprint(\"There are \", X_patch.shape[0], \"time series available for training and validation after patching with overlap. \\n\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"### Initializating Data Sets\n\n\n# Shuffling data\nidx = np.arange(X_patch.shape[0])\nshuffled_idx = np.random.shuffle(idx)\nX_patch = th.squeeze(X_patch[shuffled_idx,:])\nY_patch = th.squeeze(Y_patch[shuffled_idx])\n\nN_samples = X_patch.shape[0]\nN_valid = int(valid_rate*N_samples) # number of validation patches\n\nX_valid = X_patch[:N_valid,:].cpu()\nY_valid = Y_patch[:N_valid].cpu()\n\n\n# Inputs of nn.Conv1d must have the following shape: (N, C_in, *),\n#where N: number of samples (batch size), C_in: number of input channels, *: can be any dimension (150_000 for our time series)\n\n# Add channel dimension\nX_valid = th.unsqueeze(X_valid, 1)\nX_train = X_patch[N_valid:,:].cpu()\nY_train = Y_patch[N_valid:].cpu()\nN_train = X_train.shape[0]\n\nn_batch = int(X_train.shape[0] / batch_size) # Number of batch per epoch\n\nX_train = X_train.reshape([N_train, 1, patch_size])","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# *Function for Results plotting*","metadata":{}},{"cell_type":"code","source":"### Plotting Graphs\n\ndef plot_and_save_results(train_losses, valid_losses, best_mvd, best_mtd, min_tl, min_vl, mm, model_name):\n\n    best_epoch = np.fromiter(valid_losses, dtype=float).argmin()\n    \n    fig, ax = plt.subplots(1, 2, figsize = (20,10))\n    plt.rcParams['font.size'] = '20'\n    ax[0].set(title = \"Losses: Training and Validation\") \n    ax[0].set_xlabel(\"epochs\", fontsize = 20)\n    ax[0].set_ylabel(\"MSE\", fontsize = 20)\n    ax[0].plot(train_losses,\"r\", label = \"Training\", linewidth = 3)\n    ax[0].plot(valid_losses, \"b\", label = \"Validation\", linewidth = 3)\n    ax[0].legend(loc = \"upper right\", fontsize = 18)\n    ax[0].axvline(x=int(best_epoch), color = 'black', linestyle =\"--\", linewidth = 3)\n    ax[0].annotate(\"Best epoch: {}\\nMSE_train: {:.3f}\\nMSE_valid: {:.3f}\".format(int(best_epoch), min_tl, min_vl), \\\n                   xy = (0.5,0.5), xycoords = 'axes fraction')\n    \n    best_mvd_plot = mm.inverse_transform(best_mvd.reshape(-1, 1))\n    best_mtd_plot = mm.inverse_transform(best_mtd.reshape(-1, 1))\n\n    N_valid = best_mvd_plot.shape[0]\n    \n    mean = float(np.mean(best_mvd_plot))\n    std_dev = float(np.std(best_mvd_plot))\n    kurt = float(kurtosis(best_mvd_plot))\n    skewn = float(skew(best_mvd_plot))\n    q1 = float(np.quantile(best_mvd_plot, 0.25))\n    median = float(np.quantile(best_mvd_plot, 0.5))\n    q3 = float(np.quantile(best_mvd_plot, 0.75))\n    mae = np.absolute(best_mvd_plot).sum() / N_valid\n    mse = np.square(best_mvd_plot).sum() / N_valid\n    \n    text = \"mean: {:.3f}\\nstd: {:.3f}\\nkurt: {:.3f}\\nskew: {:.3f}\\nq1: {:.3f}\\nmed: {:.3f}\\nq3: {:.3f}\\niqr: {:.3f}\\nmae: {:.3f}\\nmse: {:.3f}\\n\".format(mean, std_dev, kurt, skewn, q1, median, q3, q3-q1, mae, mse)\n    \n    ax[1].hist(best_mvd_plot, alpha = 0.3, label = \"validation set\", bins = 100, density = True, range = (-16, 16))\n    ax[1].set(title = \"TTF error distributions at best epoch\")\n    ax[1].set_xlabel(\"Error (seconds)\", fontsize = 20)\n    ax[1].set_ylabel(\"Density\", fontsize = 20)\n    ax[1].hist(best_mtd_plot, alpha = 0.3, label = \"training set\", bins = 100, density = True, range = (-16, 16))\n    ax[1].annotate(text, xy =(-15, 0.1))\n    ax[1].legend()\n\n\n    plt.gcf()\n    plt.savefig('Models/' + model_name + '/' + model_name + \"_plot.jpg\")\n    plt.show()\n\n    model_features = {\"N_samples\": N_samples, \"N_train\": N_train, \"N_valid\": N_valid, \\\n                      \"overlap_rate\": overlap_rate, \"learning_rate\": learning_rate, \"num_epochs\": num_epochs, \\\n                     \"seed\": seed, \"batch_size\": batch_size, \"train_losses\": train_losses, \"valid_losses\": valid_losses, \\\n                     \"valid_differences\": valid_differences, \"best_mtd\": best_mtd, \"best_mvd\": best_mvd, \"min_tl\": min_tl, \\\n                      \"min_vl\": min_vl}\n\n    pickle.dump(model_features, open('Models/' + model_name + '/' + model_name + \".p\", \"wb\" ))\n    \n    model_features_display = {\"N_samples\": N_samples, \"N_train\": N_train, \"N_valid\": N_valid, \\\n                  \"overlap_rate\": overlap_rate, \"learning_rate\": learning_rate, \"num_epochs\": num_epochs, \\\n                 \"seed\": seed, \"batch_size\": batch_size}\n    \n    pickle.dump(model_features_display, open('Models/' + model_name + '/' + model_name + \"display.p\", \"wb\" ))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# *1D Convolutional Neural Network*\n\nThis section contains two 1D CNN architecture that we designed during this project, with their respective hyperparameters and the training loop.\nTo use the LSTM instead, you can skip this section and go to the next.","metadata":{}},{"cell_type":"code","source":"### Hyperparameters specific to 1D CNN\n\nlearning_rate=0.0001\nnum_epochs=500\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"### 1D Convolutional Network #1\n\nclass NN(nn.Module):\n    \n    def __init__ (self,num_classes=1):\n        super(NN, self).__init__()\n        self.fc1=nn.Conv1d(in_channels=1, out_channels=1, \n                             kernel_size=10, padding=1, stride=5)\n        self.fc2=nn.AdaptiveMaxPool1d(200)\n        self.fct3=nn.Linear(200,50)\n        self.fct4=nn.Linear(50,num_classes)\n    \n    def forward (self,x):\n        x=self.fc1(x)\n        x=self.fc2(x)\n        x=self.fct3(x)\n        x=self.fct4(x)\n        return x\n    ","metadata":{"id":"bOtetCjsDj8t"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"### 1D Convolutional Network #2\n\nclass NN_2(nn.Module):\n    def __init__(self,num_classes=1):\n        super(NN_2, self).__init__()\n        self.seq_1=nn.Sequential(nn.Conv1d(in_channels = 1, out_channels = 16, kernel_size = 10, stride = 5),\n                                 nn.ReLU(),\n                                 nn.MaxPool1d(kernel_size = 2),\n                                 nn.Conv1d(in_channels = 16, out_channels = 32, kernel_size = 10, stride = 5, \\\n                                           padding = 'valid'),\n                                 nn.ReLU(),\n                                 nn.MaxPool1d(kernel_size = 2),\n                                 nn.Conv1d(in_channels = 32, out_channels = 64, kernel_size = 10, stride = 5, \\\n                                           padding = 'valid'),\n                                 nn.ReLU()\n                                )\n        \n        self.seq_2=nn.Sequential(nn.Dropout(p = 0.4),\n                                 nn.Linear(in_features = 64, out_features = 1)\n                                )\n    def forward(self,x):\n        x=self.seq_1(x)\n        x = th.mean(x, dim = 2, keepdim = True)\n        x = th.squeeze(x)\n        x=self.seq_2(x)\n        return (x)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"### Initialize Network, Define Loss and optimizer\n\nmodel = NN_2()\n#print(model)\nmodel.to(device)\nloss=nn.MSELoss()\noptimizer=th.optim.Adam(model.parameters(),learning_rate)\n","metadata":{"id":"hH24ktcBDvrq","outputId":"55f4b5f1-f013-4e34-da59-49e8d0d19e99"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"### Plot Model Graph \n\n#from torchviz import make_dot\n#yhat = model(X_valid.to(device))\n#make_dot(yhat, params=dict(list(model.named_parameters()))).render(\"cnn_2\", format=\"png\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"### Train network\n\nth.cuda.empty_cache()\n\ntrain_losses, valid_losses = [], []\n\nbest_mvd, best_mtd = [], [] # validation, training differences at best epoch\n\nmin_vl = 1000\nmin_tl = 1000\n\nN_train_per_epoch = n_batch*batch_size \n\nstart = time.perf_counter()\n\nfor epoch in range(num_epochs):\n\n    th.cuda.empty_cache()\n    \n    # Generating Random Batches every epoch\n    random_idx = np.random.randint(0, N_train, N_train_per_epoch)\n    X_train_batch = th.reshape(X_train[random_idx], [n_batch, batch_size, 1, patch_size])\n    Y_train_batch = th.reshape(Y_train[random_idx], [n_batch, batch_size, 1, 1])\n    \n    running_loss = 0\n    _loss = 0\n\n    for ids in range (n_batch):\n        \n        # Set gradients to 0\n        optimizer.zero_grad()\n\n        # Model Prediction\n        predicted_ttfs = th.squeeze(model(X_train_batch[ids].cuda())).cpu()\n        # Loss function\n        _loss = loss(predicted_ttfs, th.squeeze(Y_train_batch[ids]))\n        \n        del predicted_ttfs\n        gc.collect()\n        \n        # Backpropagation\n        _loss.backward()\n        running_loss += _loss.item()\n        \n        # Updating Model Weights\n        optimizer.step()\n        \n    optimizer.zero_grad()\n    th.cuda.empty_cache()\n    \n    model.eval()\n    with th.no_grad():\n    \n        # Check Accuracy with Validation Set\n        \n        Y_valid_pred = th.squeeze(model(X_valid.cuda())).cpu()\n        valid_loss = loss(Y_valid_pred, Y_valid)\n        valid_differences = Y_valid_pred[:] - Y_valid[:]\n        valid_differences = valid_differences.numpy()\n        \n        valid_losses.append(valid_loss.item())\n        train_losses.append(running_loss / n_batch)\n        \n        if valid_loss < min_vl:\n            \n            min_vl = valid_loss\n            min_tl = running_loss\n            \n            best_mvd = valid_differences\n            \n            Y_final = th.squeeze(model(X_train[:int(N_train/10)].to(device))).cpu()\n            \n            best_mtd = Y_train[:int(N_train/10)] - Y_final[:]\n            best_mtd = best_mtd.cpu().detach().numpy()\n\n    model.train()\n    \n    if epoch%10 == 0:\n        print(\"Epoch: {}\\t\".format(epoch),\n                \"train Loss: {:.5f}.. \".format(train_losses[-1]),\n                \"valid Loss: {:.5f}.. \".format(valid_losses[-1])) \n\nprint(\"---------- Best : {:.3f}\".format(min(valid_losses)), \" at epoch \" \n    , np.fromiter(valid_losses, dtype=float).argmin(), \" / \", epoch + 1)\n\n\nend = time.perf_counter()\nprint(\"\\ntime elapsed: {:.3f}\".format(end - start))","metadata":{"id":"QVS0pEmAD4P4","scrolled":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#model_name = input(\"Choose a name for the model: \")\nmodel_name = \"Cnn1d\"\nos.mkdir('Models/' + model_name)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_and_save_results(train_losses, valid_losses, best_mvd, best_mtd, min_tl, min_vl, mm, model_name)\nth.save(model.state_dict(), 'Models/' + model_name + '/' + model_name)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# *Retrieve Test Data*","metadata":{}},{"cell_type":"code","source":"### Retrieve Test Data\n\n#rootpath = os.getcwd() + \"/\"\nrootpath = \"../input/LANL-Earthquake-Prediction/\"\n\nX_test_ = []\n\nfor filename in os.listdir(rootpath + \"test\"):\n    temp_df = pd.read_csv(rootpath + \"test/\" + filename)\n    X_test_.append(temp_df)\n\npatch_size = X_test_[0].shape[0]\nsample_submission = pd.read_csv(rootpath + \"sample_submission.csv\")\n\nX_test = th.zeros((len(X_test_), patch_size))\nfor i in range(len(X_test_)):\n    X_test[i,:] = th.tensor(X_test_[i][\"acoustic_data\"], dtype = th.float32)\n\nY_test = th.tensor(sample_submission[\"time_to_failure\"])\n\nN_test = X_test.shape[0]\nX_test = X_test.reshape([N_test, 1, patch_size])","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# *Submission*","metadata":{}},{"cell_type":"code","source":"### Load Model\n\nmodel = NN_2().cuda()\nmodel.load_state_dict(th.load('Models/Cnn1d_2/Cnn1d_2'))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model.eval()\nwith th.no_grad():\n    Y_test_predicted = model(X_test.cuda())","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Y_test_predicted = np.array(Y_test_predicted.cpu()).squeeze()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_submission['time_to_failure'] = Y_test_predicted","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_submission.to_csv('submission.csv')","metadata":{},"execution_count":null,"outputs":[]}]}