{"cells":[{"metadata":{"_uuid":"5516e5311b44cec16b10fa7c866b9fa33c3736e4"},"cell_type":"markdown","source":"# PLAsTiCC Astronomical Classification 2018"},{"metadata":{"_uuid":"6aa08fddc51d403113ab19cb6ab1cb80ffecb4cb"},"cell_type":"markdown","source":"I use this kernel for analysis of signals form through passband. I convert `flux`, `flux errors` and `detected` values to images with 3 channels and 6 rows. Quantity of columns depends on dataset - I binnarize time sequencies with a fixed number of bins."},{"metadata":{"trusted":true,"_uuid":"02bcbb646e8d3775efe991cc53081189237dc5ec"},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom scipy.stats import skew\nfrom matplotlib import pylab as plt","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"94f35fcb5e1e9665426d166ce5bdff1dfe792f37"},"cell_type":"code","source":"from sklearn.model_selection import StratifiedKFold","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"cbe7be57d563e48f0dc228eccc532de9b3503022"},"cell_type":"code","source":"import gc\nimport time\nimport warnings\nwarnings.simplefilter(action = 'ignore')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8089995971fcbe4c8fd97fb07a38a4db562d6738"},"cell_type":"code","source":"import tensorflow as tf\nfrom tensorflow import keras","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2d8d6154d78fa5420e2c7a736bdfdb68f8de6eec"},"cell_type":"markdown","source":"## Loading dataset"},{"metadata":{"_uuid":"3e876ca2a8545570a153bc671b2daf87bafc2edf"},"cell_type":"markdown","source":"For illustration purpose I use only galactic train subset."},{"metadata":{"trusted":true,"_uuid":"208cacc71487587788bf04a5badae415affdd257"},"cell_type":"code","source":"train_all = pd.read_csv('../input/training_set.csv')\ntrain_all.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"590463f685e54f110746898d025998c01958d0ca"},"cell_type":"code","source":"train_meta = pd.read_csv('../input/training_set_metadata.csv')\ntrain_meta.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1a6959389e64908f71b9cb2d2d9f522f59d5192a"},"cell_type":"code","source":"galactic_objects = list(train_meta[np.isnan(train_meta['distmod'])]['object_id'])\nlen(galactic_objects)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"4be4b8928e858b148bcdf6beed503a772bd73b05"},"cell_type":"code","source":"galactic_set = train_all[train_all['object_id'].isin(galactic_objects)]\ngalactic_set.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"bb6e8d9a008b396e1e7e803758b13bfad1f9a8c2"},"cell_type":"code","source":"galactic_target = train_meta[train_meta['object_id'].isin(galactic_objects)][['object_id', 'target']]\ngalactic_target.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"30d50b9989fd913afc13c0e8e10ceb283654b579"},"cell_type":"code","source":"galactic_classes = sorted(galactic_target['target'].unique())\ngalactic_classes","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"521fcb553554a2a4fa24beafb2acc2da519183c1"},"cell_type":"code","source":"del train_all, train_meta\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"aae273bd8d7b863df2fe9e8a8e7a093d402db11f"},"cell_type":"markdown","source":"## Functions"},{"metadata":{"_uuid":"7e86276d34386fe5d6bb6eead1b973b0a7fbabd5"},"cell_type":"markdown","source":"This block contains all of functions I use for preparing data, creating and training model, visualizing results."},{"metadata":{"_uuid":"bfe958a83246b792f79c2b59d40f44e8261a9d13"},"cell_type":"markdown","source":"### for preparing data"},{"metadata":{"trusted":true,"_uuid":"a5c6e605d979099b1e0f04e3409c273a26179c38"},"cell_type":"code","source":"# This function calculate some statistics for timeseries. \n# In this model I use only the median of the secont column for defining step size of binnarizing.\n# Others statistics using in another models.\n\ndef get_stats(df):\n    groups = df.groupby('passband')\n    res = groups['mjd'].apply(np.count_nonzero).values\n    res = np.vstack((res, groups['mjd'].apply(np.asarray).apply(lambda x: np.median(x[1:] - x[:-1])))) # This\n    #res = np.vstack((res, groups['flux'].apply(np.mean)))\n    #res = np.vstack((res, groups['flux'].apply(np.std)))\n    #res = np.vstack((res, groups['flux'].apply(skew)))\n    #res = np.vstack((res, groups['flux_err'].apply(np.mean)))\n    #res = np.vstack((res, groups['flux_err'].apply(np.std)))\n    #res = np.vstack((res, groups['flux_err'].apply(skew)))\n    #res = np.vstack((res, groups['detected'].apply(np.mean)))\n    #res = np.vstack((res, groups['detected'].apply(np.std)))\n    \n    return np.transpose(res)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"058e572746f744aadfcaeaa33517c092240bd0b8"},"cell_type":"code","source":"# This function converts numpy array from source dataset into 3-channels binned array with fixed width.\n# Columns in array must contain values: mjd, passband, flux, flux_err, detected.\n# Used in Batch Generator wich adds zeros to the equal length.\n\ndef to_binned_timeseries(ndar, step):\n    warnings.simplefilter(action = 'ignore')\n    \n    # the first time for object\n    start = np.min(ndar[:, 0])\n    # sequence duration for object\n    mjd_lendth = np.max(ndar[:, 0]) - start\n    # count of bins for object timeseries\n    timeseries_lendth = int(mjd_lendth / step) + 1\n    # matrix for counts in each bin for each row\n    cnt = np.zeros((6, timeseries_lendth))\n    # matrix for result with 3 channels: flux, flux_err, detected\n    # corresponds to data_format = 'channels_last' for CPU\n    res = np.zeros((6, timeseries_lendth, 3))\n    \n    # loop for rows in sourse array for calculating summs\n    for i in range(ndar.shape[0]):\n        row = ndar[i, :]\n        col_num = int((row[0] - start) / step)\n        cnt[int(row[1]), col_num] += 1\n        res[int(row[1]), col_num, 0] += row[2]\n        res[int(row[1]), col_num, 1] += row[3]\n        res[int(row[1]), col_num, 2] += row[4]\n        \n    # get mean values exclude nans\n    res[:, :, 0] /= cnt\n    res[:, :, 1] /= cnt\n    res[:, :, 2] /= cnt\n    \n    # normalizing flux channels by rows\n    for channel in range(2):\n        means = np.reshape([np.mean(res[i, ~np.isnan(res[i, :, channel]), channel]) for i in range(6)]*timeseries_lendth, \n                           (6, timeseries_lendth), order = 'F')\n        stds = np.reshape([np.std(res[i, ~np.isnan(res[i, :, channel]), channel]) for i in range(6)]*timeseries_lendth, \n                          (6, timeseries_lendth), order = 'F')\n        res[:, :, channel] = (res[:, :, channel] - means) / stds\n        \n    # replacing nans to zeros\n    res = np.nan_to_num(res)\n        \n    return res","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"98a3df5c7ec91e7f64392f1a3d0282c73adc548c"},"cell_type":"code","source":"#Calculating of constant for this model and dataset\n\nMAX_LENDTH = -1\nfor obj in galactic_objects:\n    ndar = galactic_set[galactic_set['object_id'] == obj][['mjd', 'passband', 'flux', 'flux_err', 'detected']]\n    stats = get_stats(ndar)\n    data = to_binned_timeseries(ndar.values, np.median(stats[:, 1]))\n    if data.shape[1] > MAX_LENDTH:\n        MAX_LENDTH = data.shape[1]\n        \nprint('Count of columns in `image`:', MAX_LENDTH)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b85d0a4245c656795c0e0fed7834b02e31bccdb4"},"cell_type":"markdown","source":"### for CNN-model"},{"metadata":{"trusted":true,"_uuid":"856f252e61965185dccd52683838f7f086582b10"},"cell_type":"code","source":"class BatchGenerator(keras.utils.Sequence):\n    \n    def __init__(self, X, y, batch_size = 32, predict = False):\n        self.X = X\n        self.index = list(X['object_id'].unique())\n        self.y = y\n        self.batch_size = batch_size\n        self.predict = predict\n\n        if not predict:\n            self.on_epoch_end()\n        \n    def __getitem__(self, index_batch):\n        idx = self.index[index_batch * self.batch_size : (index_batch + 1) * self.batch_size]\n        batch = np.zeros((len(idx), 6, MAX_LENDTH, 3))\n        if not self.predict:\n            target = np.zeros((len(idx), self.y.shape[1]))\n        \n        for i, obj in enumerate(idx):\n            ndar = self.X[self.X['object_id'] == obj][['mjd', 'passband', 'flux', 'flux_err', 'detected']]\n            stats = get_stats(ndar) # for defining step size\n            data = to_binned_timeseries(ndar.values, np.median(stats[:, 1]))\n            \n            # adding zeros to MAX_LENDTH\n            if data.shape[1] < MAX_LENDTH:\n                data = np.concatenate((data, np.zeros((6, MAX_LENDTH - data.shape[1], 3))), axis = 1)\n                \n            batch[i] = data\n            if not self.predict:\n                target[i] = self.y.loc[obj].values\n\n        if self.predict:\n            return batch\n        else:\n            return batch, target\n        \n    def on_epoch_end(self):\n        if not self.predict:\n            np.random.shuffle(self.index)\n        \n    def __len__(self):\n        if self.predict:\n            return int(np.ceil(len(self.index) / self.batch_size))\n        else:\n            return int(len(self.index) / self.batch_size)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"848ee2c42da06263dd956f90a43c45632710fdb8"},"cell_type":"code","source":"def get_model(class_cnt, input_shape, dropout = .5):\n    inputs = keras.layers.Input(shape = input_shape)\n    \n    # Convolutional block\n    \n    x = inputs\n    \n    x = keras.layers.Conv2D(filters = 8, kernel_size = 1, padding = 'same', use_bias = False, \n                            kernel_initializer = keras.initializers.he_normal(seed = 0),\n                            kernel_regularizer = keras.regularizers.l2(0.01))(x)\n    x = keras.layers.BatchNormalization(momentum = 0.9)(x)\n    x = keras.layers.LeakyReLU(alpha = .3)(x)\n    \n    x = keras.layers.Conv2D(filters = 16, kernel_size = 3, padding = 'same', use_bias = False, \n                            kernel_initializer = keras.initializers.he_normal(seed = 0),\n                            kernel_regularizer = keras.regularizers.l2(0.01))(x)\n    x = keras.layers.BatchNormalization(momentum = 0.9)(x)\n    x = keras.layers.LeakyReLU(alpha = .2)(x)\n    \n    x = keras.layers.Conv2D(filters = 32, kernel_size = 3, padding = 'same', use_bias = False, \n                            kernel_initializer = keras.initializers.he_normal(seed = 0),\n                            kernel_regularizer = keras.regularizers.l2(0.01))(x)\n    x = keras.layers.BatchNormalization(momentum = 0.9)(x)\n    x = keras.layers.LeakyReLU(alpha = .1)(x)\n    \n    x = keras.layers.GlobalAveragePooling2D()(x)\n    \n    # Dense block\n    \n    x = keras.layers.Flatten()(x)\n    \n    x = keras.layers.Dense(class_cnt * 4, \n                           kernel_initializer = keras.initializers.he_normal(seed = 0),\n                           kernel_regularizer = keras.regularizers.l2(0.01))(x)\n    x = keras.layers.BatchNormalization()(x)\n    x = keras.layers.LeakyReLU(alpha = 0)(x)\n    \n    x = keras.layers.Dropout(dropout)(x)\n    \n    x = keras.layers.Dense(class_cnt * 2, \n                           kernel_initializer = keras.initializers.he_normal(seed = 0),\n                           kernel_regularizer = keras.regularizers.l2(0.01))(x)\n    x = keras.layers.BatchNormalization()(x)\n    x = keras.layers.LeakyReLU(alpha = 0)(x)\n    \n    x = keras.layers.Dropout(dropout)(x)\n    \n    outputs = keras.layers.Dense(class_cnt, \n                                 kernel_initializer = keras.initializers.he_normal(seed = 0), \n                                 activation = 'softmax')(x)\n    \n    return keras.Model(inputs, outputs)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ee08ee1e21f50d2c8e8d28c0a8c09cf74a710388"},"cell_type":"code","source":"OPTIMIZER = keras.optimizers.Adam(lr = 0.0005)\n\n# The behavior of this metric completely coincides with the custom function. Only absolute values differ.\nLOSS = 'categorical_crossentropy' \n\nMETRICS = ['categorical_accuracy']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ec2863abea40cfc300197c2116da48f2c461c794"},"cell_type":"code","source":"# Cross-validation for Keras model\ndef cv_scores(num_folds, classes, model_file_prefix, \n              X_train, y_train, \n              early_stopping = -1,\n              n_epoch = 50, batch_size = 32, rs = 0):\n    \n    def lr_schedule_cosine(x):\n        return .001 * (np.cos(np.pi * x / n_epoch) + 1.) / 2\n    \n    warnings.simplefilter('ignore')\n    \n    print(\"Starting cross-validation at {} with random_state {}\".format(time.ctime(), rs))\n\n    folds = StratifiedKFold(n_splits = num_folds, shuffle = True, random_state = rs)\n        \n    # Create arrays to store results\n    train_pred = pd.DataFrame(columns = classes, index = y_train['object_id'])\n    valid_pred = pd.DataFrame(columns = classes, index = y_train['object_id'])\n    \n    y = pd.get_dummies(y_train.set_index('object_id')['target']).reset_index()\n        \n    histories = {}\n\n    # Cross-validation cycle\n    for n_fold, (train_idx, valid_idx) in enumerate(folds.split(np.zeros(y_train.shape[0]), \n                                                                y_train.set_index('object_id'))):\n        print('--- Fold {} started at {}'.format(n_fold, time.ctime()))\n        \n        # Preparing data\n        train_y = y.iloc[train_idx]\n        train_objects = train_y['object_id'].values\n        train_x = X_train[X_train['object_id'].isin(train_objects)]\n        \n        valid_y = y.iloc[valid_idx]\n        valid_objects = valid_y['object_id'].values\n        valid_x = X_train[X_train['object_id'].isin(valid_objects)]\n        \n        # Defining new model\n        train_gen = BatchGenerator(train_x, train_y.set_index('object_id'), batch_size = batch_size)\n        valid_gen = BatchGenerator(valid_x, valid_y.set_index('object_id'), batch_size = batch_size)\n        \n        model = get_model(len(classes), (6, MAX_LENDTH, 3))\n            \n        model.compile(optimizer = OPTIMIZER, loss = LOSS, metrics = METRICS)\n        \n        model_file = model_file_prefix + '_fold_' + str(n_fold) + '.h5'\n        callbacks = [\n                keras.callbacks.LearningRateScheduler(lr_schedule_cosine),\n                keras.callbacks.ModelCheckpoint(filepath = model_file, \n                                                monitor = 'val_loss', \n                                                save_best_only = True, save_weights_only = True)\n        ]\n        \n        if early_stopping > 0:\n            callbacks.append(keras.callbacks.EarlyStopping(monitor = 'val_loss', patience = early_stopping))\n\n        # Fitting model\n        model.fit_generator(train_gen, validation_data = valid_gen, callbacks = callbacks, epochs = n_epoch)\n        \n        histories[n_fold] = model.history.history\n        \n        # Prediction for train and valid data\n\n        train_gen = BatchGenerator(train_x, None, batch_size = 1, predict = True)\n        valid_gen = BatchGenerator(valid_x, None, batch_size = 1, predict = True)\n        \n        model.load_weights(model_file)\n        \n        train_pred.loc[train_objects] = pd.DataFrame(model.predict_generator(train_gen), \n                                                  columns = classes, index = train_objects) \n        valid_pred.loc[valid_objects] = pd.DataFrame(model.predict_generator(valid_gen), \n                                                  columns = classes, index = valid_objects)\n        \n        del train_x, train_y, valid_x, valid_y\n        gc.collect()\n        \n    return train_pred, valid_pred, histories","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ba0d671efd3aaa21efc7189254480fdd253f70f8"},"cell_type":"code","source":"# Custom score function for galactic subset\n\ndef weighted_multiclass_logloss(y_true, y_pred):\n    class_weights = [1, 1, 1, 1, 1]\n    \n    y_pred_clip = np.clip(a = y_pred, a_min = 1e-15, a_max = 1 - 1e-15)\n    \n    loss = np.sum(y_true * y_pred_clip.applymap(np.log), axis = 0)\n    loss /= np.sum(y_true, axis = 0)\n    loss *= class_weights\n    return -(np.sum(loss) / np.sum(class_weights))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"502a90652eb176e53075b238aaf4fe76b1e33d71"},"cell_type":"code","source":"# Function for visualizing history\n\ndef plot_history(hist):\n    n_folds = len(hist)\n    _, axes = plt.subplots(n_folds, 2, figsize = (25, 7 * n_folds))\n    \n    for row in range(n_folds):\n        n_epoch = len(hist[row][\"loss\"])\n        axes[row, 0].plot(range(1, n_epoch + 1), hist[row][\"loss\"], label = \"Train loss\")\n        axes[row, 0].plot(range(1, n_epoch + 1), hist[row][\"val_loss\"], label = \"Valid loss\")\n        axes[row, 0].legend()\n    \n        axes[row, 1].plot(range(1, n_epoch + 1), hist[row]['categorical_accuracy'], label = 'Train accuracy')\n        axes[row, 1].plot(range(1, n_epoch + 1), hist[row]['val_categorical_accuracy'], label = 'Valid accuracy')\n        axes[row, 1].legend()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"74f4878eb340b5db53c1c041c07e440b00140aca"},"cell_type":"markdown","source":"## Train model with cross-validation"},{"metadata":{"scrolled":true,"trusted":true,"_uuid":"7749c817bdb69d3f7d4283602bbfbe2c98cda5c0"},"cell_type":"code","source":"gal_train_pred, gal_valid_pred, gal_histories = cv_scores(num_folds = 2, \n                                                          n_epoch = 100, \n                                                          model_file_prefix = 'galactic',\n                                                          classes = galactic_classes, \n                                                          X_train = galactic_set, \n                                                          y_train = galactic_target)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"22fdc0c8bb0bd1f3d36aeb9899fc0b45ea98c742"},"cell_type":"code","source":"plot_history(gal_histories)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"68e6bc68167e96948ebcca6342d77643fe15f4d1"},"cell_type":"code","source":"gal_train_pred.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"eaf48adc8ea63dbae5273027b2d7626e40693deb"},"cell_type":"code","source":"gal_valid_pred.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6e30f5225e35d010eace8baf34a4cbfdf326f6b2"},"cell_type":"code","source":"y = pd.get_dummies(galactic_target.set_index('object_id')['target'])\ny.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8315c023591fdc6a3f0fc94d313510cf41cb1b10"},"cell_type":"code","source":"print('Custom score for train: ', weighted_multiclass_logloss(y, gal_train_pred))\nprint('Custom score for valid: ', weighted_multiclass_logloss(y, gal_valid_pred))","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.5"}},"nbformat":4,"nbformat_minor":1}