{"cells":[{"metadata":{},"cell_type":"markdown","source":"\nThis is a simple kernel which demonstrates that it would have been **possible to get a top 25 result (private LB 2.40) using a simple CNN without any feature engineering and without utilizing the p4677 information about the test data**. The approach is different from all other ideas I've seen here: Everything is based on the *signal power density*, which is simply the acoustic signal amplitude squared (after subtracting the mean). I am making this kernel public for two reasons:\n* I haven't seen any other kernel using the power density\n* This seems to be the only way to use neural nets directly on the data: the power density is always > 0 and can therefore be averaged over many data points (unlike the original signal, where positive and negative amplitudes would cancel, leading to a loss of information). A 150,000 points window can therefore be efficiently downsampled to several 100 data points.\n* I achieved a (unfortunately late :) ) top 25 result without much effort. The classic story: I didn't like my public LB for this method and didn't choose it... Maybe someone else is interested in improving this result? Using the p4677 info and optimizing the CNN via CV might lead to a (late) top result...\n\nAll that was done to achieve this result is the following:\n1. split the train data into non-overlapping 150k data point windows, analogous to the test data\n2. remove the global mean from train/test, square the result, reshape the squared data into (300,500) windows and take the average over each 300 point window, resulting in downsampled 300 point power density windows.\n3. train a CNN (3 convolution layers, 2 dense layers)\n4. multiply the predicted ttf by a factor which accounts for the fact that the mean power density in test is ~5% higher than in train (more on that below).\n\nI probably didn't even use the best possible CNN topology / parameters.\nHope that this solution is of interest! Thanks to the organizers for this great and frustrating :) competition, and congratulation to the winners!\n"},{"metadata":{"trusted":true,"_uuid":"7a6d5ef92752b663ef29b835642045c69f5892df","_kg_hide-input":false},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nfrom keras.models import Sequential\nfrom keras.layers import Dense, Dropout, Conv1D, MaxPooling1D, Flatten\nfrom keras.regularizers import l1","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Read the train data and prepare train_X and train_y:"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":false},"cell_type":"code","source":"df_train = pd.read_csv('../input/train.csv', dtype = {'acoustic_data': np.int16, 'time_to_failure': np.float32} ) # float32 is enough :)\n\ntrain_X = np.resize(df_train['acoustic_data'].values, (len(df_train['acoustic_data']) // 150000, 150000)).astype(np.float32) # rearrange into 150k windows\ntrain_X = (train_X - np.mean(train_X)) ** 2 # calculate power density\ntrain_X = np.mean(train_X.reshape(-1,300,500), axis=2) # downsample 500x\n\ntrain_y = np.resize(df_train['time_to_failure'].values, (len(df_train['time_to_failure']) // 150000, 150000))[:,-1] # train_y is ttf on right window edge\n\ndel df_train # free some memory\n\nprint(train_X.shape, train_y.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Read the submission file and the test data and prepare test_X:"},{"metadata":{"trusted":true},"cell_type":"code","source":"df_subm = pd.read_csv('../input/sample_submission.csv')\n\ntest_X = []\n\nfor fname in df_subm['seg_id'].values:\n    test_X.append(pd.read_csv('../input/test/' + fname + '.csv').acoustic_data.values.astype(np.int16))\ntest_X = np.array(test_X).astype(np.float32)\n\ntest_X = (test_X - np.mean(test_X)) ** 2  # calculate power density\ntest_X = np.mean(test_X.reshape(-1,300,500), axis=2) # downsample 500x\n\nprint(test_X.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This is what the individual windows look like: some examples from train(left) and test(right):"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, axes = plt.subplots(1,2)\nfig.set_size_inches(16,5)\naxes[0].grid(True)\naxes[1].grid(True)\n\naxes[0].plot(train_X[28]);\naxes[0].plot(train_X[38]);\naxes[1].plot(test_X[28]);\naxes[1].plot(test_X[38]);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now something important: as we will see below, the model (like all other models in this competition, due to the underlying physics) is actually unable to truly predict ttf for individual quake cycles. For ~ the first half of the cycle, the acoustic data is always the same in every cycle. There is no information in the data which tells whether there will be a minor quake in 3 sec or a major one in 5 sec. So the model *has* to predict a kind of average ttf over all cycles (see plot of ttf below).\n\nNow the test data has a higher mean power density than the train data. This means that the cycles in test will on average be longer. It is therefore reasonable to multiply the predicted ttfs by this ratio, which is ~1.05:"},{"metadata":{"trusted":true},"cell_type":"code","source":"ratio_mean_train_test = np.mean(test_X) / np.mean(train_X)\nprint('Ratio of mean power in train/test : ', np.mean(test_X) / np.mean(train_X))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now prepare the X data for the CNN:"},{"metadata":{"trusted":true},"cell_type":"code","source":"def prepare_X_for_cnn(data):\n    data = np.log10(data) # take log10 to handle the huge peaks\n    data -= np.mean(data) # remove mean\n    data /= np.std(data) # set std to 1.0\n    data = np.expand_dims(data, axis=-1) # reshaping for CNN/RNN\n    return data\n\ntrain_X = prepare_X_for_cnn(train_X)\ntest_X  = prepare_X_for_cnn(test_X)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Another look at the same windows as above after the CNN preparation. I've added another window to train which contains a major quake. The log10 scales everything such that a CNN can handle it:"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, axes = plt.subplots(1,2)\nfig.set_size_inches(16,5)\naxes[0].grid(True)\naxes[1].grid(True)\n\naxes[0].plot(train_X[28]);\naxes[0].plot(train_X[38]);\naxes[0].plot(train_X[29]);\naxes[1].plot(test_X[28]);\naxes[1].plot(test_X[38]);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now the CNN: The parameters have been set according to the results of a quick RandomGridCV. I'm sure that there are better topologies:"},{"metadata":{"trusted":true},"cell_type":"code","source":"def model_cnn():\n    model = Sequential([\n        Conv1D(filters=16, kernel_size=3, activation='relu', input_shape = train_X.shape[1:]),\n        MaxPooling1D(2),\n        Conv1D(filters=128, kernel_size=3, activation='relu'),\n        MaxPooling1D(2),\n        Conv1D(filters=16, kernel_size=3, activation='relu'),\n        MaxPooling1D(2),\n        Flatten(),\n        Dropout(0.1),\n        Dense(16, activation='relu', kernel_regularizer=l1(0.01)),\n        Dense(16, activation='relu', kernel_regularizer=l1(0.01)),\n        Dense(1, activation='linear') # regression\n    ])\n    model.compile(\n        loss='mse',\n        optimizer='adam',\n        metrics=['mae']\n    )\n    return model\n\nmodel = model_cnn()\nmodel.summary()\nmodel.save_weights('/tmp/model_weights_init.h5')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Same CV approach here: 16 epochs, average of 64 runs as the net isn't really stable after 16 epochs. batch_size etc. could still be tuned..."},{"metadata":{"trusted":true},"cell_type":"code","source":"test_y_pred = []\nnum_iter = 64\n\nfor i in range(num_iter):\n    model.load_weights('/tmp/model_weights_init.h5')\n    model.fit(train_X, train_y, epochs=16,  verbose=0)\n    test_y_pred.append(model.predict(test_X))\n    \ntest_y_pred = np.array(test_y_pred).reshape(num_iter,-1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Average the ttfs: "},{"metadata":{"trusted":true},"cell_type":"code","source":"test_y_pred_avg = np.mean(test_y_pred, axis=0)\ntest_y_pred_avg","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This are the predicted train ttfs. Note that the predictions are an average over all cycles. Training the CNN for more epochs would lead to overfitting, as it would simply memorize all train data points:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, axes = plt.subplots(1,1)\nfig.set_size_inches(16,5)\naxes.grid(True)\n\naxes.plot(model.predict(train_X));\naxes.plot(train_y);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Multiply by the corrective factor (see above) and submit:"},{"metadata":{"trusted":true},"cell_type":"code","source":"df_subm['time_to_failure'] = test_y_pred_avg * ratio_mean_train_test\ndf_subm.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"That's all. Hope this was helpful, maybe even for the organizers :). Remember: the CNN didn't get any info on mean, std, etc. directly but learned all that *from the shape of the peaks*."}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}