{"cells":[{"metadata":{"_uuid":"1fa3f89331fab7bfe8fc7f54cca431bd96e02d66"},"cell_type":"markdown","source":"## Keras + RNN(GRU) to handle passbands as timeseries\nMost of the kernels I saw here were not treating passbands as timeseries, but calculating meta feature such as mean/min/max of flux.\nHere I'm making simple kernel which uses all the timeseries data as input of RNN(GRU).\nI hope this helps your exploration. I'd appreciate your insights and feedback."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport sys\nimport time\nimport tensorflow as tf\nfrom keras.utils import to_categorical\nfrom keras.callbacks import EarlyStopping, TensorBoard\nfrom keras.layers import GRU, Dense, Activation, Dropout, CuDNNGRU, concatenate, Input\nfrom keras.callbacks import ReduceLROnPlateau,ModelCheckpoint\nfrom keras.models import Sequential, Model\nfrom keras.optimizers import Adam\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.model_selection import StratifiedKFold\nfrom collections import Counter\nimport matplotlib.pyplot as plt\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nimport os\nprint(os.listdir(\"../input\"))\n\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"train = pd.read_csv('../input/training_set.csv')\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"080a2678472db3cce9f0e5e7f9b547e9c7e00001"},"cell_type":"markdown","source":"## Standardize\nStandardize input before we process data. NN likes properly scaled data."},{"metadata":{"trusted":true,"_uuid":"ce4c94327d5e13dd826aaff09ce23a8c6b9ec1ef"},"cell_type":"code","source":"ss = StandardScaler()\ntrain[['mjd', 'flux', 'flux_err']] = ss.fit_transform(train[['mjd', 'flux', 'flux_err']])\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b81257eb21c6486d30cc5eb15a05cc1c567a61db"},"cell_type":"markdown","source":"## Magical transformation\nAs you see above, train data can be grouped by object_id, passband, list of [flux] order by mjd.\n\nMeaning 1 train input row would be<br>\n``\nobject_id = 615, passband=1, flux0 = 0.228287, flux1 = xxx, ... flux72 = xxx\n``\n\nSimilary we can have similar row for flux_err.<br>\n``\nobject_id = 615, passband=1, flux_err0 = 0.0.005226, flux_err = xxx, ... flux_err = xxx\n``\n\nThe following methods do the transformation!"},{"metadata":{"trusted":true,"_uuid":"38d3c824aad15a17881738d5d39307e07a2dc7fa"},"cell_type":"code","source":"def get_timeseries_by_column_name(train, passband, column_name):\n    train_passband = train[train.passband == passband]\n    ## This is where magic happens\n    train_column_timeseries = train_passband.groupby(['object_id', 'passband'])[column_name].apply(lambda df: df.reset_index(drop=True)).unstack()\n    train_column_timeseries.fillna(0, inplace=True)\n    train_column_timeseries  = train_column_timeseries.reset_index()\n    train_column_timeseries['feature_id'] = column_name\n    return train_column_timeseries\n\ndef get_timeseries(train, passband):\n    df = pd.concat([get_timeseries_by_column_name(train, passband, column_name) for column_name in ['mjd', 'flux', 'flux_err', 'detected']])\n    df = df.sort_values(['object_id', 'passband'])\n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"bd134ab1642d875fdbd5767546ab35c6dafd3762"},"cell_type":"code","source":"# One example for flux where passband=0\nget_timeseries_by_column_name(train, passband=0, column_name='flux').head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d54727b145e78207edef22014bebd2c7a34b33be"},"cell_type":"code","source":"# for each passband we get input rows for mjd, flux, flux_err, detected\ntrain_list = []\nfor passband in range(0, 6):\n    train_passband =  get_timeseries(train, passband=passband)\n    train_list.append(train_passband)\n    print(\"train_list[%d]\" % passband, train_list[passband].shape)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9b238dee407fcbe17c81523805843373f12870d9"},"cell_type":"markdown","source":"As you can see the shapes above, each passband __may have__ different length of time series. eg) passband0 has len = 75, but passband1 = 61."},{"metadata":{"_uuid":"994aff83da25bd928bce8e59bdf1df936744a4d5"},"cell_type":"markdown","source":"## Transform for RNN input\n Input to RNN is (batch_size, timesteps, input_dim).\n - batch_size is batch_size :)\n - timesteps is 75 for passband0,\n - input_dim is number of feature the timeseries has, here we have 4 ['mjd', 'flux', 'flux_err', 'detected']."},{"metadata":{"trusted":true,"_uuid":"f9a832583b79c3f1df6ffb3cafa6286ffeab2cb4"},"cell_type":"code","source":"num_features = len(['mjd', 'flux', 'flux_err', 'detected'])\ndrop_features= [ \"feature_id\", \"object_id\", \"passband\"]\nX_train_list = []\nfor passband in range(0, 6):\n    num_columns = len(train_list[passband].columns) - len(drop_features)\n    print(\"num_columns[%d]:\" % passband, num_columns)\n    X_train_list.append(train_list[passband].drop(drop_features, axis=1).values.reshape(-1, num_features, num_columns).transpose(0, 2, 1))\n    print(\"X_train_list[%d].shape:\" % passband, X_train_list[passband].shape)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6d61b2f0975655e8471603f5d045ac9fb1e87b3c"},"cell_type":"markdown","source":"## Merge target value from metadata"},{"metadata":{"trusted":true,"_uuid":"1abbead77165840e740ede35386c936b358b262b"},"cell_type":"code","source":"meta_train = pd.read_csv('../input/training_set_metadata.csv')\nmeta_train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"db6c8c26ef049bd7ce692393cf153f0564410d9f"},"cell_type":"code","source":"# List all classes\nclasses = sorted(meta_train.target.unique())\nclasses","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"bbdbf039787a75805465dacd60e4e040efd78afc"},"cell_type":"code","source":"class_map = dict()\nfor i,val in enumerate(classes):\n    class_map[val] = i\nclass_map","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f781a650d97a4461e591a796d1d6a8dd38f55b47"},"cell_type":"code","source":"# We only need target for earch object_id, so using feature_id=='flux'\nmerged_train = train_list[0][train_list[0].feature_id == 'flux'].merge(meta_train, on='object_id', how='left')\nmerged_train = merged_train.drop(['ra',\t'decl',\t'gal_l',\t'gal_b',\t'ddf',\t'hostgal_specz',\t'hostgal_photoz',\t'hostgal_photoz_err',\t'distmod',\t'mwebv'], axis=1)\nmerged_train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"0b35c27a3554a88096f7a3369b103818a69a732c"},"cell_type":"code","source":"targets = merged_train.target\ntarget_map = np.zeros((targets.shape[0],))\ntarget_map = np.array([class_map[val] for val in targets])\nY = to_categorical(target_map)\nY.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5b47d28ec9f85edbfe04e092f0b4f3a85c4be3ed"},"cell_type":"code","source":"def multi_weighted_logloss(y_ohe, y_p):\n    \"\"\"\n    @author olivier https://www.kaggle.com/ogrellier\n    multi logloss for PLAsTiCC challenge\n    \"\"\"\n    classes = [6, 15, 16, 42, 52, 53, 62, 64, 65, 67, 88, 90, 92, 95]\n    class_weight = {6: 1, 15: 2, 16: 1, 42: 1, 52: 1, 53: 1, 62: 1, 64: 2, 65: 1, 67: 1, 88: 1, 90: 1, 92: 1, 95: 1}\n    # Normalize rows and limit y_preds to 1e-15, 1-1e-15\n    y_p = np.clip(a=y_p, a_min=1e-15, a_max=1-1e-15)\n    # Transform to log\n    y_p_log = np.log(y_p)\n    # Get the log for ones, .values is used to drop the index of DataFrames\n    # Exclude class 99 for now, since there is no class99 in the training set \n    # we gave a special process for that class\n    y_log_ones = np.sum(y_ohe * y_p_log, axis=0)\n    # Get the number of positives for each class\n    nb_pos = y_ohe.sum(axis=0).astype(float)\n    # Weight average and divide by the number of positives\n    class_arr = np.array([class_weight[k] for k in sorted(class_weight.keys())])\n    y_w = y_log_ones * class_arr / nb_pos    \n    loss = - np.sum(y_w) / np.sum(class_arr)\n    return loss\n\n# https://www.kaggle.com/c/PLAsTiCC-2018/discussion/69795\ndef mywloss(y_true,y_pred):  \n    yc=tf.clip_by_value(y_pred,1e-15,1-1e-15)\n    loss=-(tf.reduce_mean(tf.reduce_mean(y_true*tf.log(yc),axis=0)/wtable))\n    return loss\n\ndef plot_loss_acc(history):\n    plt.plot(history.history['loss'][1:])\n    plt.plot(history.history['val_loss'][1:])\n    plt.title('model loss')\n    plt.ylabel('val_loss')\n    plt.xlabel('epoch')\n    plt.legend(['train','Validation'], loc='upper left')\n    plt.show()\n    \n    plt.plot(history.history['acc'][1:])\n    plt.plot(history.history['val_acc'][1:])\n    plt.title('model Accuracy')\n    plt.ylabel('val_acc')\n    plt.xlabel('epoch')\n    plt.legend(['train','Validation'], loc='upper left')\n    plt.show()\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"92c130ae6cfee6849441ddfc86087584f6fd6b5c"},"cell_type":"markdown","source":"## RNN(GRU) begins here :)\nNote that\n- I'm using Keras functional API which supports multples inputs.\n- I intentionally __chose small values for GRU parameters and epochs for demo purpose__, as GRU is slow to train."},{"metadata":{"trusted":true,"_uuid":"70b8fdd12e217c131cca662844115115160fe8be"},"cell_type":"code","source":"n_classes = len(classes)\n\ndef weight_variable(shape, name=None):\n    return np.random.normal(scale=.01, size=shape)\n\ndef build_model():\n    def basic_layer(input_):\n        output = GRU(64,\n                     kernel_initializer=weight_variable,\n                     dropout=0.5,\n                     recurrent_dropout=0.5,\n                     return_sequences=True)(input_)\n        output = Dropout(0.5)(output)\n        output = GRU(32, return_sequences=True, dropout=0.5, recurrent_dropout=0.5)(output)\n        output = GRU(16, dropout=0.5, recurrent_dropout=0.5)(output)\n        output = Dense(32)(output)\n        return output\n    \n    # Keras functional API supports multiple inputs!\n    input0 = Input(shape=(None, num_features), dtype='float32', name='passband0')\n    input1 = Input(shape=(None, num_features), dtype='float32', name='passband1')\n    input2 = Input(shape=(None, num_features), dtype='float32', name='passband2')\n    input3 = Input(shape=(None, num_features), dtype='float32', name='passband3')\n    input4 = Input(shape=(None, num_features), dtype='float32', name='passband4')\n    input5 = Input(shape=(None, num_features), dtype='float32', name='passband5')\n    \n    merged_output = concatenate([basic_layer(input0),basic_layer(input1),basic_layer(input2),basic_layer(input3),basic_layer(input4),basic_layer(input5)])\n    merged_output = Dense(64)(merged_output)\n    final_output = Dense(len(classes), activation='softmax')(merged_output)\n    return Model(inputs=[input0, input1, input2, input3, input4, input5], outputs=[final_output])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1e732ae6746341a30383919037e4c0bbb057d401"},"cell_type":"code","source":"y_count = Counter(target_map)\nwtable = np.zeros((len(classes),))\nfor i in range(len(classes)):\n    wtable[i] = y_count[i] / target_map.shape[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"305d1e36ee958915d9f2faccedd0f5808091bb04"},"cell_type":"code","source":"batch_size = 512\ny_map = target_map\ny_categorical = Y\nfolds = StratifiedKFold(n_splits=5, shuffle=True, random_state=1)\nstart = time.time()\nclfs = []\noof_preds = np.zeros((len(X_train_list[0]), len(classes)))\nepochs = 10 #5000\n\nfor fold_, (trn_, val_) in enumerate(folds.split(y_map, y_map)):\n    checkPoint = ModelCheckpoint('./keras.model',monitor='val_loss',mode = 'min', save_best_only=True, verbose=0)\n    X0 = X_train_list[0]\n    X1 = X_train_list[1]\n    X2 = X_train_list[2]\n    X3 = X_train_list[3]\n    X4 = X_train_list[4]\n    X5 = X_train_list[5]\n    x_train0, x_train1, x_train2, x_train3, x_train4, x_train5, y_train = X0[trn_], X1[trn_], X2[trn_], X3[trn_], X4[trn_], X5[trn_], y_categorical[trn_]\n    x_valid0, x_valid1, x_valid2, x_valid3, x_valid4, x_valid5, y_valid = X0[val_], X1[val_], X2[val_], X3[val_], X4[val_], X5[val_], y_categorical[val_]\n    \n    model = build_model()    \n    optimizer = Adam(lr=0.001, beta_1=0.9, beta_2=0.999)#\n    early_stopping = EarlyStopping(monitor='loss', patience=30, verbose=1)\n    model.compile(loss=mywloss, optimizer=optimizer, metrics=['accuracy'])\n    history = model.fit([x_train0, x_train1, x_train2, x_train3, x_train4, x_train5], y_train,\n                    validation_data=[[x_valid0, x_valid1, x_valid2, x_valid3, x_valid4, x_valid5], y_valid], \n                    epochs=epochs,\n                        batch_size=batch_size,\n                    shuffle=True,verbose=1,callbacks=[checkPoint, early_stopping])      \n    plot_loss_acc(history)\n    \n    print('Loading Best Model')\n    model.load_weights('./keras.model')\n    # # Get predicted probabilities for each class\n    oof_preds[val_, :] = model.predict([x_valid0, x_valid1, x_valid2, x_valid3, x_valid4, x_valid5],batch_size=batch_size)\n    print(multi_weighted_logloss(y_valid, model.predict([x_valid0, x_valid1, x_valid2, x_valid3, x_valid4, x_valid5],batch_size=batch_size)))\n    clfs.append(model)\n    \nprint('MULTI WEIGHTED LOG LOSS : %.5f ' % multi_weighted_logloss(y_categorical,oof_preds))\n\nelapsed_time = time.time() - start\nprint(\"elapsed_time:\", elapsed_time)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a2aa9b4ee897889ea605e5b9a5db1de0b7db9ae2"},"cell_type":"markdown","source":"## Ideas for improvement\n- Add one more input for metadata.\n- Derive some features as timeseries and feed into the RNN?"}],"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}