{"cells":[{"metadata":{"_uuid":"e154a47bf09b8770980486e87786317a1b3038e1"},"cell_type":"markdown","source":"### This kernel presents some changes on the amazing Bruno Aquino's kernel. \n### The main difference is related to the Matthews Correlation Coefficent metric on Keras but there are also some minor improvements on parallelization and reproducibility.\n### Please upvote if it was helpful."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import pandas as pd\nimport pyarrow.parquet as pq # Used to read the data\nimport os \nimport numpy as np\nfrom joblib import Parallel, delayed\nimport tensorflow as tf\nfrom tensorflow import set_random_seed\nfrom keras.layers import * # Keras is the most friendly Neural Network library, this Kernel use a lot of layers classes\nfrom keras.models import Model\nfrom tqdm import tqdm # Processing time measurement\nfrom sklearn.model_selection import train_test_split \nfrom sklearn.metrics import matthews_corrcoef\nfrom keras import backend as K # The backend give us access to tensorflow operations and allow us to create the Attention class\nfrom keras import optimizers # Allow us to access the Adam class to modify some parameters\nfrom sklearn.model_selection import GridSearchCV, StratifiedKFold # Used to use Kfold to train our model\nfrom keras.callbacks import * # This object helps the model to train in a smarter way, avoiding overfitting","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6e6379386e44afc69bee8895a52da22199e888fb"},"cell_type":"code","source":"# select how many folds will be created\nN_SPLITS = 5\n# it is just a constant with the measurements data size\nsample_size = 800000\n\nMAX_THREADS = 2\nRANDOM_SEED = 2019","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c60f7af0118ba87d89ef4cc7e8299fe0b746e96e"},"cell_type":"code","source":"np.random.seed(RANDOM_SEED)\nset_random_seed(RANDOM_SEED)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"35ec9ae6f4489e1ef589f006269452c1f53dc911"},"cell_type":"markdown","source":"### The standart way keras calculates an epoch's loss/metric is by taking the average value of the loss/metric on batches. While it works well on losses calculated by average (like cross entropy), it creates problems in metrics that must be calculated over the entire dataset (like F1 score or Matthews Correlation).  \n### To overcome that problems, we are going to use a statefull metric that is calculated over a whole epoch. A stateful metric on keras is a special layer that allow the running of cumulative operations on each batch. The following implementation requires tensorflow as backend."},{"metadata":{"trusted":true,"_uuid":"c3340ee96becb5ca8f075d9c44b7df383ddba5ee"},"cell_type":"code","source":"class StatefullMCC(Layer):\n    def __init__(self, thresholds, **kwargs):\n        super(StatefullMCC, self).__init__(**kwargs)\n        self.thresholds = thresholds\n        self.stateful = True\n        self.name='matthews_correlation'\n\n    def reset_states(self):\n        K.get_session().run(tf.variables_initializer(self.local_variable))\n\n    def metric_variable(self, shape, dtype, validate_shape=True, name=None):\n        return tf.Variable(\n                np.zeros(shape),\n                dtype=dtype,\n                trainable=False,\n                collections=[tf.GraphKeys.LOCAL_VARIABLES],\n                validate_shape=validate_shape,\n                name=name,\n                )\n\n    def initialize_vars(self, y_true, y_pred):\n        self.tp = [None]*len(self.thresholds)\n        self.tn = [None]*len(self.thresholds)\n        self.fp = [None]*len(self.thresholds)\n        self.fn = [None]*len(self.thresholds)\n        for i,_ in enumerate(self.thresholds):\n            self.tp[i] = self.metric_variable(shape=[1], dtype=tf.int64, validate_shape=False, name='tp%d'%i)\n            self.tn[i] = self.metric_variable(shape=[1], dtype=tf.int64, validate_shape=False, name='tn%d'%i)\n            self.fp[i] = self.metric_variable(shape=[1], dtype=tf.int64, validate_shape=False, name='fp%d'%i)\n            self.fn[i] = self.metric_variable(shape=[1], dtype=tf.int64, validate_shape=False, name='fn%d'%i)\n\n            tp_op = tf.assign_add(self.tp[i], tf.count_nonzero(y_true * y_pred[i], axis=0))\n            self.add_update(tp_op)\n            tn_op = tf.assign_add(self.tn[i], tf.count_nonzero((1-y_true)*(1-y_pred[i]), axis=0))\n            self.add_update(tn_op)\n            fp_op = tf.assign_add(self.fp[i], tf.count_nonzero((1-y_true)*y_pred[i], axis=0))\n            self.add_update(fp_op)\n            fn_op = tf.assign_add(self.fn[i], tf.count_nonzero(y_true*(1-y_pred[i]), axis=0))\n            self.add_update(fn_op)\n\n        self.local_variable = tf.get_collection(tf.GraphKeys.LOCAL_VARIABLES)\n\n    def __call__(self, y_true, y_pred):\n        rounded_preds = [K.cast(K.greater_equal(y_pred, t), 'float32') for t in self.thresholds]\n        self.initialize_vars(y_true, rounded_preds)\n        mcc_vec = []\n        for i,_ in enumerate(self.thresholds):\n            num = tf.cast((self.tp[i] * self.tn[i] - self.fp[i] * self.fn[i]), 'float32')\n            den = K.sqrt(tf.cast((self.tp[i] + self.fp[i]) * (self.tp[i] + self.fn[i]) * (self.tn[i] + self.fp[i]) * (self.tn[i] + self.fn[i]), 'float32'))\n            mcc_vec.append(num/(den + tf.constant(K.epsilon())))\n        return K.max(mcc_vec)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"eda7ea366117d1ce8e5fce69e5bba333821d8b48"},"cell_type":"code","source":"# https://www.kaggle.com/suicaokhoailang/lstm-attention-baseline-0-652-lb\n\nclass Attention(Layer):\n    def __init__(self, step_dim,\n                 W_regularizer=None, b_regularizer=None,\n                 W_constraint=None, b_constraint=None,\n                 bias=True, **kwargs):\n        self.supports_masking = True\n        self.init = initializers.glorot_uniform(RANDOM_SEED)\n\n        self.W_regularizer = regularizers.get(W_regularizer)\n        self.b_regularizer = regularizers.get(b_regularizer)\n\n        self.W_constraint = constraints.get(W_constraint)\n        self.b_constraint = constraints.get(b_constraint)\n\n        self.bias = bias\n        self.step_dim = step_dim\n        self.features_dim = 0\n        super(Attention, self).__init__(**kwargs)\n\n    def build(self, input_shape):\n        assert len(input_shape) == 3\n\n        self.W = self.add_weight((input_shape[-1],),\n                                 initializer=self.init,\n                                 name='{}_W'.format(self.name),\n                                 regularizer=self.W_regularizer,\n                                 constraint=self.W_constraint)\n        self.features_dim = input_shape[-1]\n\n        if self.bias:\n            self.b = self.add_weight((input_shape[1],),\n                                     initializer='zero',\n                                     name='{}_b'.format(self.name),\n                                     regularizer=self.b_regularizer,\n                                     constraint=self.b_constraint)\n        else:\n            self.b = None\n\n        self.built = True\n\n    def compute_mask(self, input, input_mask=None):\n        return None\n\n    def call(self, x, mask=None):\n        features_dim = self.features_dim\n        step_dim = self.step_dim\n\n        eij = K.reshape(K.dot(K.reshape(x, (-1, features_dim)),\n                        K.reshape(self.W, (features_dim, 1))), (-1, step_dim))\n\n        if self.bias:\n            eij += self.b\n\n        eij = K.tanh(eij)\n\n        a = K.exp(eij)\n\n        if mask is not None:\n            a *= K.cast(mask, K.floatx())\n\n        a /= K.cast(K.sum(a, axis=1, keepdims=True) + K.epsilon(), K.floatx())\n\n        a = K.expand_dims(a)\n        weighted_input = x * a\n        return K.sum(weighted_input, axis=1)\n\n    def compute_output_shape(self, input_shape):\n        return input_shape[0],  self.features_dim","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"# just load train data\ndf_train = pd.read_csv('../input/metadata_train.csv')\n# set index, it makes the data access much faster\ndf_train = df_train.set_index(['id_measurement', 'phase'])\ndf_train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"bcbcc947c46a59a07fcd1ca50df6f92c0d173941"},"cell_type":"code","source":"def get_features(dataset='train', split_parts=10):\n    if dataset == 'train':\n        cache_file = 'X.npy'\n        meta_file = '../input/metadata_train.csv'\n    elif dataset == 'test':\n        cache_file = 'X_test.npy'\n        meta_file = '../input/metadata_test.csv'\n    if os.path.isfile(cache_file):\n        X = np.load(cache_file)\n        y = None\n        if dataset == 'train':\n            y = np.load('y.npy')\n    else:\n        meta_df = pd.read_csv(meta_file)\n\n        data_measurements = meta_df.pivot(index='id_measurement', columns='phase', values='signal_id')\n        data_measurements = data_measurements.values\n        data_measurements = np.array_split(data_measurements, split_parts, axis=0)\n        X = Parallel(n_jobs=min(split_parts, MAX_THREADS), verbose=1)(delayed(prep_data)(p, dataset) for p in data_measurements)\n        try:\n            y = meta_df.loc[meta_df['phase']==0, 'target'].values\n        except:\n            y = None\n        X = np.concatenate(X, axis=0)\n\n        if dataset == 'train':\n            np.save(\"X.npy\",X)\n            np.save(\"y.npy\",y)\n        elif dataset == 'test':\n            np.save(\"X_test.npy\",X)\n    return X, y\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"26df6c7fbfecd537404866faec13d1238ae3ebc6"},"cell_type":"code","source":"# in other notebook I have extracted the min and max values from the train data, the measurements\nmax_num = 127\nmin_num = -128","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7b0717b14bcfcba1f48d33c8161ae51c778687af"},"cell_type":"code","source":"# This function standardize the data from (-128 to 127) to (-1 to 1)\n# Theoretically it helps in the NN Model training, but I didn't tested without it\ndef min_max_transf(ts, min_data, max_data, range_needed=(-1,1)):\n    if min_data < 0:\n        ts_std = (ts + abs(min_data)) / (max_data + abs(min_data))\n    else:\n        ts_std = (ts - min_data) / (max_data - min_data)\n    if range_needed[0] < 0:    \n        return ts_std * (range_needed[1] + abs(range_needed[0])) + range_needed[0]\n    else:\n        return ts_std * (range_needed[1] - range_needed[0]) + range_needed[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c6137bbbe75c3a1509a5f98e08805dbbd492aa37"},"cell_type":"code","source":"# This is one of the most important peace of code of this Kernel\n# Any power line contain 3 phases of 800000 measurements, or 2.4 millions data \n# It would be praticaly impossible to build a NN with an input of that size\n# The ideia here is to reduce it each phase to a matrix of <n_dim> bins by n features\n# Each bean is a set of 5000 measurements (800000 / 160), so the features are extracted from this 5000 chunk data.\ndef transform_ts(ts, n_dim=160, min_max=(-1,1)):\n    # convert data into -1 to 1\n    ts_std = min_max_transf(ts, min_data=min_num, max_data=max_num)\n    # bucket or chunk size, 5000 in this case (800000 / 160)\n    bucket_size = int(sample_size / n_dim)\n    # new_ts will be the container of the new data\n    new_ts = []\n    # this for iteract any chunk/bucket until reach the whole sample_size (800000)\n    for i in range(0, sample_size, bucket_size):\n        # cut each bucket to ts_range\n        ts_range = ts_std[i:i + bucket_size]\n        # calculate each feature\n        mean = ts_range.mean()\n        std = ts_range.std() # standard deviation\n        std_top = mean + std # I have to test it more, but is is like a band\n        std_bot = mean - std\n        # I think that the percentiles are very important, it is like a distribuiton analysis from eath chunk\n        percentil_calc = np.percentile(ts_range, [0, 1, 25, 50, 75, 99, 100]) \n        max_range = percentil_calc[-1] - percentil_calc[0] # this is the amplitude of the chunk\n        relative_percentile = percentil_calc - mean # maybe it could heap to understand the asymmetry\n        # now, we just add all the features to new_ts and convert it to np.array\n        new_ts.append(np.concatenate([np.asarray([mean, std, std_top, std_bot, max_range]),percentil_calc, relative_percentile]))\n    return np.asarray(new_ts)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7460e718a605803f1d9e4fbec61750a0deb02a47"},"cell_type":"code","source":"def prep_data(signal_ids, dataset=\"train\"):\n    signal_ids_all = np.concatenate(signal_ids)\n    if dataset == \"train\":\n        praq_data = pq.read_pandas('../input/train.parquet', columns=[str(i) for i in signal_ids_all]).to_pandas()\n    elif dataset == \"test\":\n        praq_data = pq.read_pandas('../input/test.parquet', columns=[str(i) for i in signal_ids_all]).to_pandas()\n    else:\n        raise ValueError(\"Unknown dataset\")\n    X = []\n    for sids in signal_ids:\n        data = praq_data[[str(s) for s in sids]].values.T\n        X_signal = [transform_ts(signal) for signal in data]\n        X_signal = np.concatenate(X_signal, axis=1)\n        X.append(X_signal)\n    X = np.asarray(X)\n    return X","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"52dc826ab9ee1dd56c9fb29bd5c1b2d26b5928bf"},"cell_type":"code","source":"X, y = get_features(\"train\", split_parts=6)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"51ad0e25b00536de6170168499923d82ae1d735f","scrolled":true},"cell_type":"code","source":"# The X shape here is very important. It is also important undertand a little how a LSTM works\n# X.shape[0] is the number of id_measuremts contained in train data\n# X.shape[1] is the number of chunks resultant of the transformation, each of this date enters in the LSTM serialized\n# This way the LSTM can understand the position of a data relative with other and activate a signal that needs\n# a serie of inputs in a specifc order.\n# X.shape[3] is the number of features multiplied by the number of phases (3)\nprint(X.shape, y.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"289bc7d1ab8048a60025801b457f8df1d848acbc"},"cell_type":"code","source":"# This is NN LSTM Model creation\ndef model_lstm(input_shape):\n    # The shape was explained above, must have this order\n    inp = Input(shape=(input_shape[1], input_shape[2],))\n    \n    init_glorot_uniform = initializers.glorot_uniform(seed=RANDOM_SEED)\n    init_orthogonal = initializers.orthogonal(seed=RANDOM_SEED)\n    \n    # This is the LSTM layer\n    # Bidirecional implies that the 160 chunks are calculated in both ways, 0 to 159 and 159 to zero\n    # although it appear that just 0 to 159 way matter, I have tested with and without, and tha later worked best\n    # 128 and 64 are the number of cells used, too many can overfit and too few can underfit\n    x = Bidirectional(CuDNNLSTM(128, return_sequences=True, kernel_initializer=init_glorot_uniform, recurrent_initializer=init_orthogonal))(inp)\n    # The second LSTM can give more fire power to the model, but can overfit it too\n    x = Bidirectional(CuDNNLSTM(64, return_sequences=True, kernel_initializer=init_glorot_uniform, recurrent_initializer=init_orthogonal))(x)\n    # Attention is a new tecnology that can be applyed to a Recurrent NN to give more meanings to a signal found in the middle\n    # of the data, it helps more in longs chains of data. A normal RNN give all the responsibility of detect the signal\n    # to the last cell. Google RNN Attention for more information :)\n    x = Attention(input_shape[1])(x)\n    # A intermediate full connected (Dense) can help to deal with nonlinears outputs\n    x = Dense(64, activation=\"relu\", kernel_initializer=init_glorot_uniform)(x)\n    # A binnary classification as this must finish with shape (1,)\n    x = Dense(1, activation=\"sigmoid\", kernel_initializer=init_glorot_uniform)(x)\n    model = Model(inputs=inp, outputs=x)\n    # Pay attention in the addition of matthews_correlation metric in the compilation, it is a success factor key\n    model.compile(loss='binary_crossentropy', optimizer='adam', metrics=[StatefullMCC(np.linspace(0.45,0.55,11))])\n    \n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8d6f4ca319c383b1b4f671a37c5a324136e7a466"},"cell_type":"code","source":"# Here is where the training happens\n\n# First, create a set of indexes of the 5 folds\nsplits = list(StratifiedKFold(n_splits=N_SPLITS, shuffle=True, random_state=RANDOM_SEED).split(X, y))\npreds_val = []\ny_val = []\n# Then, iteract with each fold\n# If you dont know, enumerate(['a', 'b', 'c']) returns [(0, 'a'), (1, 'b'), (2, 'c')]\nfor idx, (train_idx, val_idx) in enumerate(splits):\n    K.clear_session() # I dont know what it do, but I imagine that it \"clear session\" :)\n    print(\"Beginning fold {}\".format(idx+1))\n    # use the indexes to extract the folds in the train and validation data\n    train_X, train_y, val_X, val_y = X[train_idx], y[train_idx], X[val_idx], y[val_idx]\n    # instantiate the model for this fold\n    model = model_lstm(train_X.shape)\n    # This checkpoint helps to avoid overfitting. It just save the weights of the model if it delivered an\n    # validation matthews_correlation greater than the last one.\n    ckpt = ModelCheckpoint('weights_{}.h5'.format(idx), save_best_only=True, save_weights_only=True, verbose=1, monitor='val_matthews_correlation', mode='max')\n    # Train, train, train\n    model.fit(train_X, train_y, batch_size=128, epochs=50, validation_data=[val_X, val_y], callbacks=[ckpt])\n    # loads the best weights saved by the checkpoint\n    model.load_weights('weights_{}.h5'.format(idx))\n    # Add the predictions of the validation to the list preds_val\n    preds_val.append(model.predict(val_X, batch_size=512))\n    # and the val true y\n    y_val.append(val_y)\n\n# concatenates all and prints the shape    \npreds_val = np.concatenate(preds_val)[...,0]\ny_val = np.concatenate(y_val)\npreds_val.shape, y_val.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d28151fd0be9fd9762f3f55e307d82f89bfbd291"},"cell_type":"code","source":"# The output of this kernel must be binary (0 or 1), but the output of the NN Model is float (0 to 1).\n# So, find the best threshold to convert float to binary is crucial to the result\n# this piece of code is a function that evaluates all the possible thresholds from 0 to 1 by 0.01\ndef threshold_search(y_true, y_proba):\n    thresholds = np.linspace(0.0,1.0,101)\n    scores = [matthews_corrcoef(y_true, (y_proba > t).astype(np.uint8)) for t in thresholds]\n    best_idx = np.argmax(scores)\n    return thresholds[best_idx], scores[best_idx]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6fee7f722ed08bc1453a822a4371ed2d48e08abc"},"cell_type":"code","source":"best_threshold, best_score = threshold_search(y_val, preds_val)\nprint(best_threshold, best_score)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ae9bd3fa9d8c0781c0708846bb7f2a9f9e6cbd3c"},"cell_type":"code","source":"%%time\n# Now load the test data\n# This first part is the meta data, not the main data, the measurements\nmeta_test = pd.read_csv('../input/metadata_test.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3eb186d032f79c99ffba05dd1a7fabb77e13cec5"},"cell_type":"code","source":"meta_test = meta_test.set_index(['signal_id'])\nmeta_test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6f8e94387f625bff0a9a6289e1ee038908bc5856"},"cell_type":"code","source":"%%time\nX_test_input, _ = get_features(\"test\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"af9aa6b2b8f8a2beda1a02ff998e3072fcad8d06"},"cell_type":"code","source":"X_test_input.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"cfd265d3e07c4cc1679d2c4d55fe7de631c813e7"},"cell_type":"code","source":"submission = pd.read_csv('../input/sample_submission.csv')\nprint(len(submission))\nsubmission.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2f7342296138f6bfd3e9cedd029e1035de3b98fc"},"cell_type":"code","source":"preds_test = []\nfor i in range(N_SPLITS):\n    model.load_weights('weights_{}.h5'.format(i))\n    pred = model.predict(X_test_input, batch_size=300, verbose=1)\n    pred_3 = []\n    for pred_scalar in pred:\n        for i in range(3):\n            pred_3.append(pred_scalar)\n    preds_test.append(pred_3)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9f76c471eaf983707d446c5081ab3d50c4e40ea5"},"cell_type":"code","source":"preds_test = (np.squeeze(np.mean(preds_test, axis=0)) > best_threshold).astype(np.int)\npreds_test.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b35723f85d494b4b6ec630dd7c79135a110a4062"},"cell_type":"code","source":"submission['target'] = preds_test\nsubmission.to_csv('submission.csv', index=False)\nsubmission.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d7600d0093a9880003240ef9ce0a1f1303e4d982"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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}