{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install -q wtte","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Inspired by this [kernel](https://www.kaggle.com/mayer79/rnn-starter-for-huge-time-series) and [this package](https://ragulpr.github.io/2016/12/22/WTTE-RNN-Hackless-churn-modeling/) (and [this example](https://github.com/gm-spacagna/deep-ttf/blob/master/notebooks/Keras-WTT-RNN%20Engine%20failure.ipynb))"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom tqdm.auto import tqdm\n\nimport keras\nimport wtte.wtte as wtte\n\nfrom tensorflow import set_random_seed\nset_random_seed(5944)\n\n%matplotlib inline\nimport matplotlib.pyplot as plt\n\nK = keras.backend","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","collapsed":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"cell_type":"markdown","source":"# Custom fn"},{"metadata":{"trusted":true},"cell_type":"code","source":"def _keras_unstack_hack(ab):\n    \"\"\"Implements tf.unstack(y_true_keras, num=2, axis=-1).\n       Keras-hack adopted to be compatible with theano backend.\n    \"\"\"\n    ndim = len(K.int_shape(ab))\n    if ndim == 0:\n        print('can not unstack with ndim=0')\n    else:\n        a = ab[..., 0]\n        b = ab[..., 1]\n    return a, b\n\ndef weibull_loglik_discrete(y_true, ab_pred, name=None):\n    y_ = y_true[:, 0]\n    u_ = y_true[:, 1]\n    a_ = ab_pred[:, 0]\n    b_ = ab_pred[:, 1]\n\n    hazard0 = K.pow((y_ + 1e-35) / a_, b_)\n    hazard1 = K.pow((y_ + 1) / a_, b_)\n\n    return -1 * K.mean(u_ * K.log(k.exp(hazard1 - hazard0) - 1.0) - hazard1)\n\n\"\"\"\n    Not used for this model, but included in case somebody needs it\n    For math, see https://ragulpr.github.io/assets/draft_master_thesis_martinsson_egil_wtte_rnn_2016.pdf (Page 35)\n\"\"\"\ndef weibull_loglik_continuous(y_true, ab_pred, name=None):\n    y_ = y_true[:, 0]\n    u_ = y_true[:, 1]\n    a_ = ab_pred[:, 0]\n    b_ = ab_pred[:, 1]\n\n    ya = (y_ + 1e-35) / a_\n    return -1 * K.mean(u_ * (K.log(b_) + b_ * K.log(ya)) - K.pow(ya, b_))\n\n\n\"\"\"\n    Custom Keras activation function, outputs alpha neuron using exponentiation and beta using softplus\n\"\"\"\ndef activate(ab):\n    a = K.exp(ab[:, 0])\n    b = K.softplus(ab[:, 1])\n\n    a = K.reshape(a, (K.shape(a)[0], 1))\n    b = K.reshape(b, (K.shape(b)[0], 1))\n\n    return K.concatenate((a, b), axis=1)\n\ndef output_lambda(x, init_alpha=1.0, max_beta_value=5.0, max_alpha_value=None):\n    \"\"\"Elementwise (Lambda) computation of alpha and regularized beta.\n\n        Alpha: \n        (activation) \n        Exponential units seems to give faster training than \n        the original papers softplus units. Makes sense due to logarithmic\n        effect of change in alpha. \n        (initialization) \n        To get faster training and fewer exploding gradients,\n        initialize alpha to be around its scale when beta is around 1.0,\n        approx the expected value/mean of training tte. \n        Because we're lazy we want the correct scale of output built\n        into the model so initialize implicitly; \n        multiply assumed exp(0)=1 by scale factor `init_alpha`.\n\n        Beta: \n        (activation) \n        We want slow changes when beta-> 0 so Softplus made sense in the original \n        paper but we get similar effect with sigmoid. It also has nice features.\n        (regularization) Use max_beta_value to implicitly regularize the model\n        (initialization) Fixed to begin moving slowly around 1.0\n\n        Assumes tensorflow backend.\n\n        Args:\n            x: tensor with last dimension having length 2\n                with x[...,0] = alpha, x[...,1] = beta\n\n        Usage:\n            model.add(Dense(2))\n            model.add(Lambda(output_lambda, arguments={\"init_alpha\":100., \"max_beta_value\":2.0}))\n        Returns:\n            A positive `Tensor` of same shape as input\n    \"\"\"\n    a, b = _keras_unstack_hack(x)\n\n    # Implicitly initialize alpha:\n    if max_alpha_value is None:\n        a = init_alpha * K.exp(a)\n    else:\n        a = init_alpha * K.clip(x=a, min_value=K.epsilon(),\n                                max_value=max_alpha_value)\n\n    m = max_beta_value\n    if m > 1.05:  # some value >>1.0\n        # shift to start around 1.0\n        # assuming input is around 0.0\n        _shift = np.log(m - 1.0)\n\n        b = K.sigmoid(b - _shift)\n    else:\n        b = K.sigmoid(b)\n\n    # Clipped sigmoid : has zero gradient at 0,1\n    # Reduces the small tendency of instability after long training\n    # by zeroing gradient.\n    b = m * K.clip(x=b, min_value=K.epsilon(), max_value=1. - K.epsilon())\n\n    x = K.stack([a, b], axis=-1)\n\n    return x","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def extract_features(z):\n     return np.c_[z.mean(axis=1), \n                  z.min(axis=1),\n                  z.max(axis=1),\n                  z.std(axis=1)]\n\ndef create_X(x, last_index=None, n_steps=150, step_length=1000):\n    if last_index == None:\n        last_index=len(x)\n       \n    assert last_index - n_steps * step_length >= 0\n\n    # Reshaping and approximate standardization with mean 5 and std 3.\n    temp = (x[(last_index - n_steps * step_length):last_index].reshape(n_steps, -1) - 5 ) / 3\n    \n    # Extracts features of sequences of full length 1000, of the last 100 values and finally also \n    # of the last 10 observations. \n    return np.c_[extract_features(temp),\n                 extract_features(temp[:, -step_length // 10:]),\n                 extract_features(temp[:, -step_length // 100:])]\n\n\ndef weibull_loss_discrete(y_true, y_pred, name=None):\n    \"\"\"calculates a keras loss op designed for the sequential api.\n    \n        Discrete log-likelihood for Weibull hazard function on censored survival data.\n        For math, see \n        https://ragulpr.github.io/assets/draft_master_thesis_martinsson_egil_wtte_rnn_2016.pdf (Page 35)\n        \n        Args:\n            y_true: tensor with last dimension having length 2\n                with y_true[:,...,0] = time to event, \n                     y_true[:,...,1] = indicator of not censored\n                \n            y_pred: tensor with last dimension having length 2 \n                with y_pred[:,...,0] = alpha, \n                     y_pred[:,...,1] = beta\n\n        Returns:\n            A positive `Tensor` of same shape as input\n            \n    \"\"\"    \n    y,u = _keras_unstack_hack(y_true)\n    a,b = _keras_unstack_hack(y_pred)\n\n    hazard0 = K.pow((y + 1e-35) / a, b)\n    hazard1 = K.pow((y + 1.0) / a, b)\n    \n    loglikelihoods = u * K.log(K.exp(hazard1 - hazard0) - 1.0) - hazard1\n    loss = -1 * K.mean(loglikelihoods)\n    return loss","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Import data"},{"metadata":{"trusted":true},"cell_type":"code","source":"float_data = pd.read_csv('../input/train.csv',\n                         dtype={'acoustic_data': np.float32, 'time_to_failure': np.float32}).values","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Modelling"},{"metadata":{},"cell_type":"markdown","source":"## Bookkeeping"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Query \"create_X\" to figure out the number of features\nn_features = create_X(float_data[0:150000]).shape[1]\nprint(\"Our RNN is based on %i features\"% n_features)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def generator(data, min_index=0, max_index=None, batch_size=16, n_steps=150, step_length=1000, test=False):\n    if max_index is None:\n        max_index = len(data) - 1\n     \n    while True:\n        # Pick indices of ending positions\n        rows = np.random.randint(min_index + n_steps * step_length, max_index, size=batch_size)\n         \n        # Initialize feature matrices and targets\n        samples = np.zeros((batch_size, n_steps, n_features))\n        targets = np.zeros((batch_size, 2))\n        \n        for j, row in enumerate(rows):\n            samples[j] = create_X(data[:, 0], last_index=row, n_steps=n_steps, step_length=step_length)\n            targets[j, 0] = data[row - 1, 1]\n            if test:\n                targets[j, 1] = 1\n            else:\n                targets[j, 1] = 1 if targets[j, 0] > 0.03 else 0\n        yield samples, targets","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Position of second (of 16) earthquake. Used to have a clean split\n# between train and validation\nsecond_earthquake = 50085877\nfloat_data[second_earthquake, 1]\n\nbatch_size = 32\nmask_value = -99\ninit_alpha = -1.0 / np.log(1.0 - 1.0/ (5 + 1.0) )\ninit_alpha = init_alpha / n_features\n\n# Initialize generators\ntrain_gen = generator(float_data, batch_size=batch_size) # Use this for better score\n# train_gen = generator(float_data, batch_size=batch_size, min_index=second_earthquake + 1)\nvalid_gen = generator(float_data, batch_size=batch_size, max_index=second_earthquake, test=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"aux, aux2 = next(train_gen)\nprint(aux.shape, aux2.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Define Model"},{"metadata":{"trusted":true},"cell_type":"code","source":"def weibull_mode(alpha, beta):\n    assert np.all(beta > 1)\n    return alpha * ((beta-1)/beta)**(1/beta)\n\ndef weibull_quantiles(alpha, beta, p):\n    return alpha*np.power(-np.log(1.0-p),1.0/beta)\n\ndef weibull_mode_K(a, b):\n    #assert np.all(beta > 1)\n    return a * K.pow((b - 1) / b, (1 / b))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def mode_mae(y_true, y_pred, name=None):\n    y, u = _keras_unstack_hack(y_true)\n    a, b = _keras_unstack_hack(y_pred)\n    pred = weibull_mode_K(a, b)\n    \n    return K.mean(K.abs(pred - y), axis=-1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# callbacks\nhistory = keras.callbacks.History()\nww = wtte.WeightWatcher()\nnanterminator = keras.callbacks.TerminateOnNaN()\ncb = [history, ww, nanterminator]\n\n# Start building our model\nmodel = keras.Sequential()\n\n# Mask parts of the lookback period that are all zeros (i.e., unobserved) so they don't skew the model\n#model.add(tfkl.Masking(mask_value=mask_value, input_shape=(None, n_features)))\nmodel.add(keras.layers.Conv1D(64, kernel_size=32, input_shape=(None, n_features)))\nmodel.add(keras.layers.Conv1D(64, kernel_size=16))\nmodel.add(keras.layers.MaxPool1D(2))\n\nmodel.add(keras.layers.Conv1D(128, kernel_size=16))\nmodel.add(keras.layers.Conv1D(128, kernel_size=12))\nmodel.add(keras.layers.MaxPool1D(2))\n\n# recurrent layer(s)\nmodel.add(keras.layers.GRU(50, activation='tanh', recurrent_dropout=0.25, dropout=0.3, input_shape=(None, n_features)))\n#model.add(keras.layers.CuDNNLSTM(100))\n#model.add(keras.layers.CuDNNGRU(30, input_shape=(None, n_features)))\n\n# We need 2 neurons to output Alpha and Beta parameters for our Weibull distribution\n\n# model.add(keras.layers.TimeDistributed(tfkl.Dense(2)))\nmodel.add(keras.layers.Dense(2))\n\n# Apply the custom activation function mentioned above\nmodel.add(keras.layers.Activation(activate))\n\nmodel.add(keras.layers.Lambda(\n    wtte.output_lambda, \n    arguments={'init_alpha': 3, \n               'max_beta_value': 500.0, \n               'scalefactor': 0.5,}))\n\n# Use the discrete log-likelihood for Weibull survival data as our loss function\nloss = wtte.loss(kind='discrete', reduce_loss=False).loss_function\n\nmodel.compile(loss=loss,\n              optimizer=keras.optimizers.SGD(lr=.01, clipnorm=2),\n              metrics=[mode_mae])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Fitting"},{"metadata":{"trusted":true},"cell_type":"code","source":"steps_train = 4196 // batch_size\nsteps_test = steps_train // 2\n\nhistory = model.fit_generator(train_gen,\n                              steps_per_epoch=1000,\n                              epochs=50,\n                              verbose=1,\n                              callbacks=cb,\n                              validation_data=valid_gen,\n                              validation_steps=200)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Evaluating"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(history.history['loss'], label='training')\nplt.plot(history.history['val_loss'],label='validation')\nplt.legend();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ww.plot()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"np.array(y_train_true).shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from tqdm.auto import tqdm\n\ny_train_true = []\ny_train_pred = []\ny_train_pred_25 = []\ny_train_pred_75 = []\n\nvalid_gen = generator(float_data, batch_size=batch_size, max_index=second_earthquake)\n\n\nfor _ in tqdm(range(500)):\n    feats, tgt = next(train_gen)\n    y_train_true.append(tgt[-1, 0])\n    preds = model.predict(feats)\n    alpha_pred, beta_pred = preds[-1, 0], preds[-1, 1]\n    y_train_pred.append(weibull_mode(alpha_pred, beta_pred))\n    y_train_pred_25.append(weibull_quantiles(alpha_pred, beta_pred, 0.25))\n    y_train_pred_75.append(weibull_quantiles(alpha_pred, beta_pred, 0.75))\n\nplt.figure(figsize=(9, 9))\nplt.scatter(np.array(y_train_true), np.array(y_train_pred))\nplt.xlabel('ytrue')\nplt.ylabel('ypred');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.metrics import mean_absolute_error\n\nmean_absolute_error(y_train_true, y_train_pred)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"# Prediction"},{"metadata":{"trusted":true},"cell_type":"code","source":"from tqdm.auto import tqdm\n\n# Load submission file\nsubmission = pd.read_csv('../input/sample_submission.csv', index_col='seg_id', dtype={\"time_to_failure\": np.float32})\n\n# Load each test data, create the feature matrix, get numeric prediction\nfor i, seg_id in enumerate(tqdm(submission.index)):\n  #  print(i)\n    seg = pd.read_csv('../input/test/' + seg_id + '.csv')\n    x = seg['acoustic_data'].values\n    preds = model.predict(np.expand_dims(create_X(x), 0))\n    alpha_pred, beta_pred = preds[-1, 0], preds[-1, 1]\n    submission.time_to_failure[i] = weibull_mode(alpha_pred, beta_pred)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Submission"},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.to_csv('submission.csv')","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}