{"cells":[{"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\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"import os\nimport math\nimport copy\nimport random\n\nfrom keras.layers import *\nfrom keras.models import Model, load_model\nfrom keras.optimizers import Adam, Nadam, SGD\nfrom keras.callbacks import EarlyStopping, ModelCheckpoint, TensorBoard\nfrom keras.utils import to_categorical\nfrom keras.preprocessing.sequence import pad_sequences\n\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import accuracy_score\nfrom sklearn.preprocessing import normalize\nfrom sklearn.metrics import confusion_matrix\nfrom sklearn.metrics import precision_recall_curve, roc_curve\nimport matplotlib.pyplot as plt\n%matplotlib inline \n\nimport seaborn as sns\nfrom sklearn.manifold import TSNE\n\nimport tensorflow as tf","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Arrays with constants"},{"metadata":{"trusted":true},"cell_type":"code","source":"classes = np.array([6, 15, 16, 42, 52, 53, 62, 64, 65, 67, 88, 90, 92, 95, 99], dtype='int32')\nclass_names = ['class_6','class_15','class_16','class_42','class_52','class_53','class_62','class_64','class_65','class_67','class_88','class_90','class_92','class_95','class_99']\nreal_class_names = ['µLens-Single','TDE','EB','SNII','SNIax','Mira','SNIbc',' KN','M-dwarf','SNIa-91bg','AGN ','SNIa','RRL','SLSN-I','class_99']\n# LSST passbands (nm)  u    g    r    i    z    y      \npassbands = np.array([357, 477, 621, 754, 871, 1004], dtype='float32')\n\nnum_models = 1\n\nlimit = 1000000 #limit of samples","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Load Training Data"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_meta = pd.read_csv('../input/PLAsTiCC-2018/training_set_metadata.csv')\ntrain_data = pd.read_csv('../input/PLAsTiCC-2018/training_set.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true},"cell_type":"code","source":"train_meta","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Meta data for graphic of SNIa objects from ZTF and LSST data**"},{"metadata":{"trusted":true,"collapsed":true},"cell_type":"code","source":"spec_meta = train_meta[train_meta['target'] == 90].copy()\nspec_meta = spec_meta[spec_meta['hostgal_specz'] <= 0.14]\nspec_meta\n# file_name = './spec_meta_lsst.csv'\n# spec_meta.to_csv(file_name)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"spec_objects = spec_meta[spec_meta['hostgal_specz'].isin([0.068 , 0.125 , 0.0498, 0.12  , 0.0763, 0.099 , 0.0669, 0.0767,\n       0.0732, 0.074 , 0.11  , 0.0845, 0.1181, 0.062 ])]['object_id'].to_numpy()\nspec_objects","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true},"cell_type":"code","source":"spec_data = train_data[train_data['object_id'].isin(spec_objects)].copy()\nspec_data\n# file_name = './spec_data_lsst.csv'\n# spec_data.to_csv(file_name)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Graphic - relation between real redshift and photometric redshift**"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(15,8))\nsns.scatterplot(data=train_meta, x=\"hostgal_specz\", y=\"hostgal_photoz\", hue=\"hostgal_photoz_err\")\nplt.savefig('./Photoz_err(specsz).png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Use only detected where flux differs from template**"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data = train_data[train_data['detected']==1]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Select two filters**"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data_ztf2 = train_data.loc[(train_data['passband'] > 0)&(train_data['passband'] < 3)]\ntrain_data_ztf2","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Select objects with more then 5 points**"},{"metadata":{"trusted":true},"cell_type":"code","source":"count = (train_data_ztf2.groupby(['object_id']).count()>=6)\ntest_objs = count[count['mjd']].index\ntest_objs\ntrain_data_ztf2 = train_data_ztf2[train_data_ztf2['object_id'].isin(test_objs)]\ntrain_data_ztf2","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Mean number of points for objects**"},{"metadata":{"trusted":true},"cell_type":"code","source":"mean_points = train_data_ztf2.groupby('object_id').count()\nmean_points['mjd'].to_numpy().mean()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Make descriptions for 7 classes**"},{"metadata":{"trusted":true},"cell_type":"code","source":"df = train_data_ztf2.copy()\nmeta = train_meta.drop(columns = ['ra', 'decl', 'gal_l', 'gal_b', 'distmod'])\ndf = df.merge(meta, how = 'left', on = 'object_id')\ndf['target'].unique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def describe(df):\n    return pd.concat([df.describe().T,\n                      df.median().rename('median')], axis=1).T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_stat = describe(df[df['target']==88])\ndf_stat.drop(columns = ['object_id', 'target'], inplace = True)\ndf_stat.to_excel(\"class_AGN.xlsx\")\ndf_stat","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Flux Density"},{"metadata":{"trusted":true},"cell_type":"code","source":"# three objects from class 42\nobjects = train_meta.loc[train_meta['target'].isin([42])]\nidxs = np.random.randint(len(objects), size=3)\nrand_objects = {int(objects.iloc[idx]['object_id']):int(objects.iloc[idx]['target']) for idx in idxs}\nprint(\"Random objects:\",rand_objects)\ndf = train_data.loc[train_data['object_id'].isin(rand_objects.keys())]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(15,10))\ngroups = df.groupby('object_id')\nfor g in groups:\n    sns.kdeplot(g[1]['flux'], shade = True, legend = False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Graphics of mean flux and lightcurves"},{"metadata":{"trusted":true},"cell_type":"code","source":"groups = train_data_ztf2.groupby('object_id')\ndf = pd.DataFrame()\nmean_g = groups.mean()\ndf['object_id'] = mean_g.index\ndf['m_flux'] = mean_g['flux'].to_numpy()\ndf['target'] = df.merge(train_meta, how = 'left', on = 'object_id')['target']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# exclude abnormally large mean flux\ndf = df[df['m_flux']<120000]\nsns.catplot(x=\"target\", y=\"m_flux\", data=df, height=7, aspect=2)\nplt.savefig('./mean_per_classes')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"spec_df = train_data_ztf2.merge(train_meta[['object_id','target']], how = 'left', on = 'object_id')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"object_ids = spec_df.groupby('target').head(1)['object_id'].to_numpy()\nobject_ids\nspec_df = spec_df.loc[spec_df['object_id'].isin(object_ids)]\nspec_df = spec_df[spec_df['target'].isin([90,42,15,67,52,62,95,88])]\nspec_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"g = sns.relplot(x=\"mjd\", y=\"flux\",hue = \"passband\", col=\"target\",palette = 'viridis', data=spec_df)\ng.fig.set_figwidth(91)\ng.fig.set_figheight(9)\nplt.savefig('./light_curves_lsst')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Preprocessing"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_freq(df):\n    \n    all_y = np.array(df['target'], dtype = 'int32')\n\n    y_count = np.unique(all_y, return_counts=True)\n\n    freq = np.ones(len(classes)-1)\n\n    freq = y_count[1]/ all_y.shape[0]\n\n    return {y_count[0][i]:freq[i] for i in range(len(classes)-1)}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true},"cell_type":"code","source":"freq_dir = get_freq(train_meta)\nfreq_dir[99] = 1.\nfreq_dir","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_keras_data(itemslist, test = False):\n    \n    if not(test):\n        print(\"TRAINING\")\n    else:\n        print(\"TESTING\")\n    # sequence_len = 256\n    sequence_len = 30\n    keys = itemslist[0].keys()\n    X = {\n            'id': np.array([i['id'] for i in itemslist], dtype='int32'),\n            'meta': np.array([i['meta'] for i in itemslist]),\n            'band': pad_sequences([i['band'] for i in itemslist], maxlen=sequence_len, dtype='int32'),\n            'hist': pad_sequences([i['hist'] for i in itemslist], maxlen=sequence_len, dtype='float32'),\n        }\n    for key in X.keys():\n        print('key = {0}, value_shape = {1}' .format(key, X[key].shape))\n\n    if not(test):\n        Y = to_categorical([i['target'] for i in itemslist], num_classes=len(classes))\n        print('target, value_shape = {0}' .format(Y.shape))\n        return X, Y\n    else: \n        return X\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_data(data_df, meta_df, test = False):\n\n    samples = []\n    groups = data_df.groupby('object_id')\n    flag = 1 #to watch on first sample dir\n\n    for g in groups:\n\n        id = g[0]\n\n        sample = {}\n        sample['id'] = int(id)\n\n        #object_id,ra,decl,gal_l,gal_b,ddf,hostgal_specz,hostgal_photoz,hostgal_photoz_err,distmod,mwebv,target\n        #615, 349.046051,-61.943836,320.796530,-51.753706,1,0.0000,0.0000,0.0000,nan,0.017,92\n        meta = meta_df.loc[meta_df['object_id'] == id]\n\n        if not(test):\n            if 'target' in meta:\n                sample['target'] = np.where(classes == int(meta['target']))[0][0]\n            else:\n                sample['target'] = len(classes) - 1\n\n        sample['meta'] = np.zeros(5, dtype = 'float32')\n        \n    \n        sample['meta'][0] = meta['ddf']\n        sample['meta'][1] = meta['hostgal_photoz']\n        sample['meta'][2] = meta['hostgal_photoz_err']\n        sample['meta'][3] = meta['mwebv']\n        sample['meta'][4] = float(meta['hostgal_photoz']) > 0\n\n        sample['specz'] = float(meta['hostgal_specz'])\n\n        z = float(sample['meta'][1])\n\n        #object_id,mjd,passband,flux,flux_err,detected\n        #615,59750.4229,2,-544.810303,3.622952,1\n\n\n        mjd      = np.array(g[1]['mjd'],      dtype='float32')\n        band     = np.array(g[1]['passband'], dtype='int32')\n        flux     = np.array(g[1]['flux'],     dtype='float32')\n        flux_err = np.array(g[1]['flux_err'], dtype='float32')\n        detected = np.array(g[1]['detected'], dtype='float32')\n\n        mjd -= mjd[0]\n        mjd /= 100 # Earth time shift in day*100\n        mjd /= (z + 1) # Object time shift in day*100\n\n\n        received_wavelength = passbands[band] # Earth wavelength in nm\n        source_wavelength = received_wavelength / (z + 1) # Object wavelength in nm\n\n\n        sample['band'] = band + 1\n        \n        flux_max = np.max(flux)\n        flux_min = np.min(flux)\n        flux_pow = flux_max - flux_min\n\n        sample['hist'] = np.zeros((flux.shape[0], 7), dtype='float32')\n        \n        sample['hist'][:,0] = mjd\n        sample['hist'][:,1] = flux/(flux_pow+1)\n        sample['hist'][:,2] = flux_err/(flux_pow+1)\n        sample['hist'][:,3] = detected\n        sample['hist'][:,4] = np.ediff1d(mjd, to_begin = [0])\n        sample['hist'][:,5] = (source_wavelength/1000)\n        sample['hist'][:,6] = (received_wavelength/1000)\n\n\n#         sample['meta'][5] = flux_pow / 10 #exclude this feature for ZTF data!!\n        \n        if flag:\n            print(\"First sample:\")\n            print(sample.keys())\n            print(\"id=\", sample['id'])\n            # print(\"target=\", sample['target'])\n            print(\"meta=\", sample['meta'])\n            print(\"specz=\", sample['specz'])\n            print(\"band shape=\", sample['band'].shape)\n            print(\"hist shape=\", sample['hist'].shape)\n            flag = 0\n            \n        samples.append(sample)\n        \n        if len(samples) % 1000 == 0:\n            print('Converting data {0}'.format(len(samples)), end='\\r')\n\n        if len(samples) >= limit:\n            break\n    \n    print('Full data number {0}'.format(len(samples)), end='\\r')\n    \n    return samples","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Augmentation"},{"metadata":{"trusted":true},"cell_type":"code","source":"def copy_sample(s, augmentate=True):\n    \n    c = copy.deepcopy(s)\n\n    if not augmentate:\n        return c\n\n    band = []\n    hist = []\n\n    drop_rate = 0.001\n\n#   drop some records\n    for k in range(s['band'].shape[0]):\n        if random.uniform(0, 1) >= drop_rate:\n            band.append(s['band'][k])\n            hist.append(s['hist'][k])\n\n    c['hist'] = np.array(hist, dtype='float32')\n    c['hist'][...,4] = np.ediff1d(c['hist'][...,0], to_begin = [0])\n    c['band'] = np.array(band, dtype='int32')\n            \n    new_z = random.normalvariate(c['meta'][1], c['meta'][2] / 1.5)\n    new_z = np.clip(new_z,0,5)\n\n    dt = (1 + c['meta'][1]) / (1 + new_z)\n    \n    c['meta'][1] = new_z\n\n    # augmentation for flux\n    c['hist'][:,1] = np.random.normal(c['hist'][:,1], c['hist'][:,2] / 1.5)\n\n    # multiply time intervals and wavelength to apply augmentation for red shift\n    c['hist'][:,0] *= dt\n    c['hist'][:,4] *= dt\n    c['hist'][:,5] *= dt\n    \n    return c","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def augmentate(samples, count):\n\n    res = []\n    index = 0\n    for s in samples:\n\n        index += 1\n        \n        if index % 1000 == 0:\n            print('Augmenting {0}/{1}   '.format(index, len(samples)), end='\\r')\n\n        for i in range(0, count):\n            res.append(copy_sample(s))\n            \n    return res","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# RNN Model"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_model(X, Y, size=80):\n\n    hist_input = Input(shape=X['hist'][0].shape, name='hist')\n    meta_input = Input(shape=X['meta'][0].shape, name='meta')\n    band_input = Input(shape=X['band'][0].shape, name='band')\n    \n    band_emb = Embedding(8, 8)(band_input)\n    \n    hist = concatenate([hist_input, band_emb])\n    hist = TimeDistributed(Dense(40, activation='relu'))(hist)\n    \n    rnn = Bidirectional(GRU(size, return_sequences=True))(hist)\n    rnn1 = SpatialDropout1D(0.5)(rnn)\n\n    gmp = GlobalMaxPool1D()(rnn1)\n    gmp1 = Dropout(0.5)(gmp)\n\n    x = concatenate([meta_input, gmp1])\n    x = Dense(128, activation='relu')(x)\n    x = Dense(128, activation='relu')(x)\n    x = Dropout(0.5)(x)\n\n    output = Dense(15, activation='softmax')(x)\n\n    model = Model(inputs=[hist_input, meta_input, band_input], outputs=output)\n    hidden_model = Model(inputs = [hist_input, meta_input, band_input], outputs = gmp)\n\n    return model, hidden_model","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Loss function**"},{"metadata":{"trusted":true},"cell_type":"code","source":"def mywloss(y_true,y_pred):\n    freq = np.array(list(freq_dir.values()))\n    yc=tf.clip_by_value(y_pred,1e-15,1-1e-15)\n    loss=-(tf.reduce_mean(tf.reduce_mean(y_true*tf.math.log(yc),axis=0)/freq))\n    return loss","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def train_model(i, samples_train, samples_valid):\n    \n    augment_count = 25\n    samples_train += augmentate(samples_train, augment_count)\n    \n    patience = len(samples_train)//2000+10\n\n    print(\"training data\")\n    train_x, train_y = get_keras_data(samples_train)\n    del samples_train\n    print(\"testing data\")\n    valid_x, valid_y = get_keras_data(samples_valid)\n    \n    #del samples_valid\n\n    model,rnn_part = get_model(train_x, train_y)\n\n    if i == 1: model.summary()\n    model.compile(optimizer='nadam', loss=mywloss, metrics=['accuracy'])\n    rnn_part.compile(optimizer='nadam', loss=mywloss, metrics=['accuracy'])\n\n    print('Training model {0} of {1}, Patience: {2}'.format(i, num_models, patience))\n    filename = './model_{0}.hdf5'.format(i)\n    max_epochs = 300\n    \n    callbacks = [EarlyStopping(patience=patience, verbose=1), ModelCheckpoint(filename, save_best_only=True)]\n    \n    model.fit(train_x, train_y, validation_data=(valid_x, valid_y), epochs=max_epochs, batch_size=5000, callbacks=callbacks, verbose=2)\n    \n    # return history\n    model.save('./model_last.hdf5')\n    rnn_part.save('./rnn_part.hdf5')\n    \n    #model = load_model(filename, custom_objects={'mywloss': mywloss})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"samples = get_data(train_data_ztf2, train_meta)\n# samples = get_data(train_data, train_meta)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Flux Density after normalization"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(15,10))\nfor obj in samples:\n    if obj['id'] in rand_objects.keys():\n        sns.kdeplot(obj['hist'][:,1], shade = True, legend = False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Training"},{"metadata":{"trusted":true},"cell_type":"code","source":"valid_size = 0.1\n\nsamples_train, samples_valid = train_test_split(samples, test_size=0.1, random_state=42)\n\ntrain_model(num_models, samples_train, samples_valid)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Loss Curves"},{"metadata":{"trusted":true},"cell_type":"code","source":"#plotting curves\nplt.plot(history.history['loss'], label='myloss(training data)')\nplt.plot(history.history['val_loss'], label='myloss(validation data)')\nplt.title('Training on data with 2 filters')\nplt.ylabel('loss value')\nplt.xlabel('No. epoch')\nplt.legend(loc=\"upper left\")\nplt.savefig('./Loss_Curves_2filters.png')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# t-SNE"},{"metadata":{"trusted":true},"cell_type":"code","source":"valid_x, valid_y = get_keras_data(samples_valid)\nfilename = './rnn_part.hdf5'\nrnn_part = load_model(filename, custom_objects={'mywloss': mywloss})\nrnn_part.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# mean number of objects in class\nneighbors = np.unique(np.argmax(valid_y, axis = 1), return_counts=True)\nneighbors[1].mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# objects per class\nclass_density = {real_class_names[num]:count for num,count in zip(neighbors[0],neighbors[1])}\nclass_density                    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# get hidden layer features\nrnn_res = rnn_part.predict(valid_x, batch_size=1000)\nprint(\"features shape = \", rnn_res.shape)\n\n# get 2d embeddings of features\ntsne = TSNE(n_components=2,perplexity = 56, n_iter = 4000, init='pca')\n\nprint(\"start fitting t-SNE\")\nres_tsne = tsne.fit_transform(rnn_res)\nprint(\"Final feature shape = \", res_tsne.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# # color palette visualization\n# num2class_name = dict(zip(classes[:-1], real_class_names))\n# print(num2class_name)\n# sns.palplot(sns.color_palette(\"Spectral\", 14))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"color_values = sns.color_palette(\"Spectral\", 14)\n\ncolor_map = dict(zip(np.arange(14), color_values))\n\n# color_list = [color_map[class_n] for class_n in np.argmax(valid_y, axis = 1)]\nclass_list = [num2class_name[classes[num]] for num in np.argmax(valid_y, axis = 1)]\n\nfig = plt.figure(figsize=(15,10))\n# plt.scatter(res_tsne[:,0], res_tsne[:,1], c = color_list, marker='.')\n# plt.legend(prop={'size':9})\nplt.title('t-SNE embedding: layer - Bidirectional GRU')\nsns.scatterplot(x=res_tsne[:,0], y=res_tsne[:,1], hue=class_list, palette=color_values, legend=\"full\")\nplt.savefig('./t-SNE_deep_layer(2).png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Graphics of embedded features for six classes**"},{"metadata":{"trusted":true},"cell_type":"code","source":"class_list = np.array(class_list)\nnames = ['TDE', 'SNIbc', 'M-dwarf', 'EB', 'SNII', 'SNIa']\nfig, axes = plt.subplots(2, 3, figsize=(20, 10))\nfor name, ax in zip(names, axes.flatten()):\n    hue = np.where(class_list == name, class_list, 'others')\n    sns.scatterplot(x=res_tsne[:,0], y=res_tsne[:,1], hue=hue, palette=color_values[5::8], ax = ax)\nplt.savefig('./t-SNE_deep_layer_6classes(2).png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Results on validation data"},{"metadata":{"trusted":true},"cell_type":"code","source":"def multi_weighted_logloss(y_valid, y_pred, freq):\n    \"\"\"\n    @author olivier https://www.kaggle.com/ogrellier\n    multi logloss for PLAsTiCC challenge\n    \"\"\"\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, 99: 1}\n    # Normalize rows and limit y_preds to 1e-15, 1-1e-15\n    y_pred = np.clip(a=y_pred, a_min=1e-15, a_max=1-1e-15)\n    # Transform to log\n    y_p_log = np.log(y_pred)\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_valid * y_p_log, axis=0)\n    # Get the number of positives for each class\n    nb_pos = freq\n\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 / y_valid.shape[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-output":false},"cell_type":"code","source":"print(\"testing data\")\nvalid_x, valid_y = get_keras_data(samples_valid)\n\n#filename = './model_1.hdf5'\nfilename = './model_last.hdf5'\n\nmodel = load_model(filename, custom_objects={'mywloss': mywloss})\n\nfreq = np.array(list(freq_dir.values())) # frequencies of classes\npreds = model.predict(valid_x, batch_size=1000)\nloss = multi_weighted_logloss(valid_y, preds, freq)\nacc = accuracy_score(np.argmax(valid_y, axis=1), np.argmax(preds,axis=1))\nprint('MW Loss: {0:.4f}, Accuracy: {1:.4f}'.format(loss, acc))\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# class_names: 'µLens-Single','TDE','EB',' SNII',' SNIax','Mira','SNIbc',' KN','M-dwarf','SNIa-91bg','AGN ','SNIa',' RRL',' SLSN-I','class_99'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Confusion matrix**"},{"metadata":{"trusted":true},"cell_type":"code","source":"cm = confusion_matrix(np.argmax(valid_y, axis = 1), np.argmax(preds, axis = 1), labels = np.arange(14))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# P vs R barplot"},{"metadata":{"trusted":true},"cell_type":"code","source":"#precision recall\nrecall = np.diag(cm) / np.sum(cm, axis = 1)\nprecision = np.diag(cm) / np.sum(cm, axis = 0)\n\nfig, ax = plt.subplots(figsize=(15,5))\nlabels = real_class_names[:-1]\nx = np.arange(len(classes)-1)  # the label locations\nwidth = 0.35  # the width of the bars\nrects1 = ax.bar(x - width/2, precision, width, label='Precision')\nrects2 = ax.bar(x + width/2, recall, width, label='Recall')\n\n# Add some text for labels, title and custom x-axis tick labels, etc.\nax.set_ylabel('Scores')\nax.set_title('Precision/Recall per class')\nax.set_xticks(x)\nax.set_xticklabels(labels)\nax.legend()\nplt.savefig('./PvsR_barplot.png')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#normalize conf matrix\ncm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]\nannot = np.around(cm, 2)\n\nfig, ax = plt.subplots(figsize=(10,9))\nsns.heatmap(cm, xticklabels=real_class_names[:-1], yticklabels=real_class_names[:-1], cmap='Blues', annot=annot)\nax.set_xlabel('Predicted Label')\nax.set_ylabel('True Label')\nax.set_aspect('equal')\nfig.tight_layout()\nplt.savefig('./conf_matrix.png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# P vs R curve"},{"metadata":{"trusted":true},"cell_type":"code","source":"# precision recall curve\nprecision = dict()\nrecall = dict()\nfig, ax = plt.subplots(figsize=(15,10))\nfor i in range(len(classes)-1):\n    precision[i], recall[i], _ = precision_recall_curve(valid_y[:, i], preds[:, i])\n    plt.plot(recall[i], precision[i], lw=2, label='class {}'.format(real_class_names[i]))\n    \nplt.xlabel(\"recall\")\nplt.ylabel(\"precision\")\nplt.legend(loc=\"best\")\nplt.title(\"precision vs recall curve\")\nplt.savefig('./PvsR_curves.png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# ROC curve"},{"metadata":{"trusted":true},"cell_type":"code","source":"# roc curve\nfpr = dict()\ntpr = dict()\n\nfig, ax = plt.subplots(figsize=(15,10))\nfor i in range(len(classes)-1):\n    fpr[i], tpr[i], _ = roc_curve(valid_y[:, i], preds[:, i])\n    plt.plot(fpr[i], tpr[i], lw=2, label='class {}'.format(real_class_names[i]))\n\nplt.xlabel(\"false positive rate\")\nplt.ylabel(\"true positive rate\")\nplt.legend(loc=\"best\")\nplt.title(\"ROC curve\")\nplt.savefig('./ROC_curves(2).png')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tf.keras.utils.plot_model(model, to_file=\"./model(2).png\", show_shapes=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Test Data"},{"metadata":{"trusted":true},"cell_type":"code","source":"test_data = pd.read_csv('../input/PLAsTiCC-2018/test_set_batch1.csv')\ntest_target = pd.read_csv('../input/PLAsTiCC-2018/sample_submission.csv')\ntest_meta = pd.read_csv('../input/PLAsTiCC-2018/test_set_metadata.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_meta.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"max_id = test_data['object_id'].max()+1\nprint(max_id)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_target = test_target.loc[test_target['object_id'] < max_id]\ntest_target.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"samples = get_data(test_data, test_meta, test = True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"testing data\")\nvalid_x = get_keras_data(samples, test = True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"valid_y = (test_target.drop(['object_id'], axis = 1)).to_numpy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"filename = './model_1.hdf5'\n#filename = './model_last.hdf5'\n\nmodel = load_model(filename, custom_objects={'mywloss': mywloss})\n\nfreq = np.array(list(freq_dir.values())) # frequencies of classes\npreds = model.predict(valid_x, batch_size=5000)\nloss = multi_weighted_logloss(valid_y, preds, freq)\nacc = accuracy_score(np.argmax(valid_y, axis=1), np.argmax(preds,axis=1))\nprint('MW Loss: {0:.4f}, Accuracy: {1:.4f}'.format(loss, acc))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cm = confusion_matrix(np.argmax(valid_y, axis = 1), np.argmax(preds, axis = 1), labels = np.arange(15))\ncm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]\nannot = np.around(cm, 2)\n\nfig, ax = plt.subplots(figsize=(10,9))\nsns.heatmap(cm, xticklabels=real_class_names, yticklabels=real_class_names, cmap='Blues', annot=annot)\nax.set_xlabel('Predicted Label')\nax.set_ylabel('True Label')\nax.set_aspect('equal')\nfig.tight_layout()","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}