{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# MSCI - CITEseq - TF / Keras Baseline\n\nSimple keras nn baseline that I intend to improve over time to match competitive models.\n\nNow with the multiome part: https://www.kaggle.com/code/lucasmorin/msci-multiome-tf-keras-nn-baseline","metadata":{"execution":{"iopub.status.busy":"2022-08-25T13:46:45.876734Z","iopub.execute_input":"2022-08-25T13:46:45.877061Z","iopub.status.idle":"2022-08-25T13:46:45.897016Z","shell.execute_reply.started":"2022-08-25T13:46:45.876992Z","shell.execute_reply":"2022-08-25T13:46:45.896432Z"}}},{"cell_type":"markdown","source":"# Imports\n\nImport base libraries, graphic libraries and modelling librairies (sklearn for Cross-validation, TF/Keras for modelling).","metadata":{}},{"cell_type":"code","source":"import numpy as np, pandas as pd\nimport glob, os, gc\n\nfrom IPython.core.display import display, HTML\nimport matplotlib.pyplot as plt, seaborn as sns\n\nfrom sklearn import preprocessing, model_selection\nimport tensorflow as tf\nfrom tensorflow import keras\nfrom keras import backend as K\nimport tensorflow_probability as tfp\n\n#set backend as float16 \nK.set_floatx('float16')\ntf.keras.mixed_precision.set_global_policy('mixed_float16')\n\nDEBUG = False\nTEST = True","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-10-05T22:04:14.253492Z","iopub.execute_input":"2022-10-05T22:04:14.253934Z","iopub.status.idle":"2022-10-05T22:04:22.651546Z","shell.execute_reply.started":"2022-10-05T22:04:14.253868Z","shell.execute_reply":"2022-10-05T22:04:22.649921Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"SUB_ALL = True\n\nif SUB_ALL:\n    df_all = pd.read_csv('../input/tf-keras-nn-base-to-work/submission_all.csv')\n    df_all.to_csv('submission_all.csv', index=False)\n    del df_all","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"(needed tor ead hdf files)","metadata":{}},{"cell_type":"code","source":"!pip install tables","metadata":{"execution":{"iopub.status.busy":"2022-09-04T20:47:46.517271Z","iopub.execute_input":"2022-09-04T20:47:46.518461Z","iopub.status.idle":"2022-09-04T20:47:58.690300Z","shell.execute_reply.started":"2022-09-04T20:47:46.518422Z","shell.execute_reply":"2022-09-04T20:47:58.689165Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# read data\n\nReading data, adding celltype as categorical integer; perform 10% sampling if DEBUG mode is enabled.","metadata":{}},{"cell_type":"code","source":"%%time \n\ntrain = pd.read_hdf(\"/kaggle/input/open-problems-multimodal/train_cite_inputs.h5\").astype('float16')\n\nmeta_data = pd.read_csv('../input/open-problems-multimodal/metadata.csv')\ntrain_meta_data = meta_data.set_index('cell_id').loc[train.index]\n\ntrain = train.values\ntrain_cat =  train_meta_data.cell_type.values\n\nlabels = pd.read_hdf(\"/kaggle/input/open-problems-multimodal/train_cite_targets.h5\").astype('float16').values","metadata":{"execution":{"iopub.status.busy":"2022-09-04T20:47:58.692108Z","iopub.execute_input":"2022-09-04T20:47:58.692518Z","iopub.status.idle":"2022-09-04T20:49:03.034034Z","shell.execute_reply.started":"2022-09-04T20:47:58.692479Z","shell.execute_reply":"2022-09-04T20:49:03.032645Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"map_cat = { 'BP':0, 'EryP':1, 'HSC':2, 'MasP':3, 'MkP':4, 'MoP':5, 'NeuP':6 }\ntrain_cat = np.array([map_cat[t] for t in train_cat])","metadata":{"execution":{"iopub.status.busy":"2022-09-04T20:49:03.036987Z","iopub.execute_input":"2022-09-04T20:49:03.037588Z","iopub.status.idle":"2022-09-04T20:49:03.054035Z","shell.execute_reply.started":"2022-09-04T20:49:03.037549Z","shell.execute_reply":"2022-09-04T20:49:03.052950Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if DEBUG:\n    idx = np.random.randint(0, train.shape[0], int(train.shape[0]/10))\n    train = train[idx]\n    train_cat = train_cat[idx]\n    labels = labels[idx]","metadata":{"execution":{"iopub.status.busy":"2022-09-04T20:49:03.057179Z","iopub.execute_input":"2022-09-04T20:49:03.057449Z","iopub.status.idle":"2022-09-04T20:49:04.713813Z","shell.execute_reply.started":"2022-09-04T20:49:03.057424Z","shell.execute_reply":"2022-09-04T20:49:04.712779Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Custom Loss\n\nI implemented the needed correlation as a custom metric. To have a decreasing loss, the standard approach is to consider 1- corr instead of corr. \nUsing only 1-corr as a metric is problematic as their might be a problem with scale. As the metric is independant of scale, the scale of the output can drift uncontrollably an cause overflow errors (exacerbated by the usage of float16). One solution is to add a bit of MSE loss. The final loss is 1 - corr + lambda * MSE where lambda is a small hand-tuned hyper-parameter.","metadata":{}},{"cell_type":"code","source":"lam = 0.03\n\ndef correlation_metric(y_true, y_pred):\n    x = tf.convert_to_tensor(y_true)\n    y = tf.convert_to_tensor(y_pred)\n    mx = K.mean(x,axis=1)\n    my = K.mean(y,axis=1)\n    mx = tf.tile(tf.expand_dims(mx,axis=1),(1,x.shape[1]))\n    my = tf.tile(tf.expand_dims(my,axis=1),(1,x.shape[1]))\n    xm, ym = (x-mx)/100, (y-my)/100\n    r_num = K.sum(tf.multiply(xm,ym),axis=1)\n    r_den = tf.sqrt(tf.multiply(K.sum(K.square(xm),axis=1), K.sum(K.square(ym),axis=1)))\n    r = tf.reduce_mean(r_num / r_den)\n    r = K.maximum(K.minimum(r, 1.0), -1.0)\n    return r\n\ndef correlation_loss(y_true, y_pred):\n    return 1 - correlation_metric(y_true, y_pred) + lam * tf.keras.losses.MeanSquaredError()(tf.convert_to_tensor(y_true),tf.convert_to_tensor(y_pred))","metadata":{"execution":{"iopub.status.busy":"2022-09-04T20:49:04.715444Z","iopub.execute_input":"2022-09-04T20:49:04.715916Z","iopub.status.idle":"2022-09-04T20:49:04.728194Z","shell.execute_reply.started":"2022-09-04T20:49:04.715873Z","shell.execute_reply":"2022-09-04T20:49:04.725517Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model","metadata":{}},{"cell_type":"markdown","source":"I start with a very vanilla MLP; I try to add a cell-type embedding layer. \nTo avoid too much drift, I scale each layer with batchnorm.\nI also add some noise to make the learning more robust.\nI initially chose 'relu' as the activation function that seems well suited to handle sparse data; 'selu' is usually better than 'relu'.","metadata":{}},{"cell_type":"code","source":"hidden_units = (128, 64, 32)\ncell_embedding_size = 2\nnoise = 0.1\n\ndef base_model():\n    \n    num_input = keras.Input(shape=(train.shape[1],), name='num_data')\n    \n    cat_input = keras.Input(shape=(1,), name='cell_id')\n\n    cell_embedded = keras.layers.Embedding(8, cell_embedding_size, input_length=1)(cat_input)\n    cell_flattened = keras.layers.Flatten()(cell_embedded)\n    \n    out = keras.layers.Concatenate()([cell_flattened, num_input])\n\n    out = keras.layers.BatchNormalization()(out)\n    out = keras.layers.GaussianNoise(noise)(out)\n    \n    for n_hidden in hidden_units:\n        out = keras.layers.Dense(n_hidden, activation='selu', kernel_regularizer = tf.keras.regularizers.L2(l2=0.01))(out)\n        out = keras.layers.BatchNormalization()(out)\n        out = keras.layers.GaussianNoise(noise)(out)\n        \n    out = keras.layers.Dense(labels.shape[1], activation='selu', name='prediction')(out)\n\n    model = keras.Model(\n        inputs = [num_input, cat_input],\n        outputs = out,\n    )\n    \n    return model","metadata":{"execution":{"iopub.status.busy":"2022-09-04T20:49:04.729607Z","iopub.execute_input":"2022-09-04T20:49:04.730525Z","iopub.status.idle":"2022-09-04T20:49:04.741657Z","shell.execute_reply.started":"2022-09-04T20:49:04.730486Z","shell.execute_reply":"2022-09-04T20:49:04.740706Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Training\n\nGeneral training loop; Data is split accordingly to CV. Then I train the model with some basic callbacks. \nThen the model is evaluated out of sample (we can check that the tf corr metric match the numpy implementation).","metadata":{}},{"cell_type":"code","source":"gc.collect()\n\nepochs = 3 if DEBUG else 1000\nn_folds = 2 if DEBUG else (3 if TEST else 3)\nn_seeds = 1 if DEBUG else (1 if TEST else 3)\n\nes = tf.keras.callbacks.EarlyStopping(\n    monitor='val_correlation_metric', min_delta=1e-05, patience=5, verbose=1,\n    mode='max', restore_best_weights = True)\n\nplateau = tf.keras.callbacks.ReduceLROnPlateau(\n    monitor='val_correlation_metric', factor=0.2, patience=3, verbose=1,\n    mode='max')\n\nkf = model_selection.ShuffleSplit(n_splits=n_folds, random_state=2020, test_size = 0.4)\n\ndf_scores = []\n\nfor fold, (cal_index, val_index) in enumerate(kf.split(range(len(train)))):\n    print(f'CV {fold}/{n_folds}')\n    \n    X_train = train[cal_index, :]\n    X_train_cat = train_cat[cal_index]\n    y_train = labels[cal_index, :]\n    \n    X_test = train[val_index, :]\n    X_test_cat = train_cat[val_index]\n    y_test = labels[val_index, :]\n    \n    \n    for seed in range(n_seeds):\n        print(f'Fold: {str(fold)} - seed: {str(seed)}')\n        key = str(fold)+'-'+str(seed)\n    \n        model = base_model()\n\n        model.compile(\n            keras.optimizers.Adam(learning_rate=1e-4),\n            loss = correlation_loss,\n            metrics = correlation_metric,\n        )\n\n        model.fit([X_train,X_train_cat], \n                  y_train, \n                  batch_size=128,\n                  epochs=epochs,\n                  validation_data=([X_test,X_test_cat], y_test),\n                  callbacks=[es, plateau],\n                  shuffle=True,\n                  verbose = 1)\n\n        output_test = model.predict([X_test, X_test_cat])\n        score = np.mean([np.corrcoef(y_test[i],output_test[i])[0,1] for i in range(len(y_test))])\n        print(f'Fold: {str(fold)} - seed: {str(seed)}: {score:.2%}')\n\n        df_scores.append((fold, seed, score))\n        model.save(f'model_cite_nn_{key}')\n    \n    tf.keras.backend.clear_session()\n    del  X_train, X_train_cat, y_train, X_test, X_test_cat, y_test\n    gc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-04T20:49:04.743398Z","iopub.execute_input":"2022-09-04T20:49:04.743829Z","iopub.status.idle":"2022-09-04T20:49:44.473177Z","shell.execute_reply.started":"2022-09-04T20:49:04.743711Z","shell.execute_reply":"2022-09-04T20:49:44.472197Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del train, labels\ngc.collect","metadata":{"execution":{"iopub.status.busy":"2022-09-04T20:49:44.480074Z","iopub.execute_input":"2022-09-04T20:49:44.480390Z","iopub.status.idle":"2022-09-04T20:49:44.494519Z","shell.execute_reply.started":"2022-09-04T20:49:44.480362Z","shell.execute_reply":"2022-09-04T20:49:44.493410Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Results","metadata":{}},{"cell_type":"code","source":"df_results = pd.DataFrame(df_scores,columns=['fold','seed','score']).pivot(index='fold',columns='seed',values='score')\n\ndf_results.loc['seed_mean']= df_results.mean(numeric_only=True, axis=0)\ndf_results.loc[:,'fold_mean'] = df_results.mean(numeric_only=True, axis=1)\ndf_results","metadata":{"execution":{"iopub.status.busy":"2022-09-04T20:49:44.498768Z","iopub.execute_input":"2022-09-04T20:49:44.499048Z","iopub.status.idle":"2022-09-04T20:49:44.528245Z","shell.execute_reply.started":"2022-09-04T20:49:44.499022Z","shell.execute_reply":"2022-09-04T20:49:44.527166Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Submission\n\nLoading and preparing test data. Inference on test data. Constitution of the first part of the submission.","metadata":{}},{"cell_type":"code","source":"%%time\n\nevaluation_ids = pd.read_csv('../input/open-problems-multimodal/evaluation_ids.csv').set_index('row_id')\nunique_ids = np.unique(evaluation_ids.cell_id)\nsubmission = pd.Series(name='target', index=pd.MultiIndex.from_frame(evaluation_ids), dtype=np.float16)\n\ndel evaluation_ids\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-04T20:49:44.529542Z","iopub.execute_input":"2022-09-04T20:49:44.529887Z","iopub.status.idle":"2022-09-04T20:51:48.086213Z","shell.execute_reply.started":"2022-09-04T20:49:44.529833Z","shell.execute_reply":"2022-09-04T20:51:48.085190Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\ntest = pd.read_hdf(\"/kaggle/input/open-problems-multimodal/test_cite_inputs.h5\").astype('float16')\nmeta_data = pd.read_csv('../input/open-problems-multimodal/metadata.csv')\ntest_meta_data = meta_data.set_index('cell_id').loc[test.index]\n\ntest = test.values\ntest_cat =  test_meta_data.cell_type.values\n\nmap_cat = { 'BP':0, 'EryP':1, 'HSC':2, 'MasP':3, 'MkP':4, 'MoP':5, 'NeuP':6}\ntest_cat = np.array([map_cat[t] for t in test_cat])","metadata":{"execution":{"iopub.status.busy":"2022-09-04T20:51:48.087566Z","iopub.execute_input":"2022-09-04T20:51:48.088209Z","iopub.status.idle":"2022-09-04T20:52:37.006000Z","shell.execute_reply.started":"2022-09-04T20:51:48.088161Z","shell.execute_reply":"2022-09-04T20:52:37.004850Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gc.collect()\n\nall_preds = []\n\nfor fold in range(n_folds):\n    for seed in range(n_seeds):\n        print(f'Preds - Fold: {str(fold)} - seed: {str(seed)}')\n        key = str(fold)+'-'+str(seed)\n        \n        model_cite = tf.keras.models.load_model(f'./model_cite_nn_{key}/', compile=False)\n\n        cite_pred = model_cite.predict([test, test_cat])\n        cite_pred = cite_pred.ravel()\n        len_cite_raveled = len(cite_pred)\n        all_preds.append(cite_pred)","metadata":{"execution":{"iopub.status.busy":"2022-09-04T20:52:37.010078Z","iopub.execute_input":"2022-09-04T20:52:37.010742Z","iopub.status.idle":"2022-09-04T20:53:27.575212Z","shell.execute_reply.started":"2022-09-04T20:52:37.010687Z","shell.execute_reply":"2022-09-04T20:53:27.574163Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del test, test_cat, cite_pred\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-04T20:53:27.578564Z","iopub.execute_input":"2022-09-04T20:53:27.579043Z","iopub.status.idle":"2022-09-04T20:53:28.411949Z","shell.execute_reply.started":"2022-09-04T20:53:27.579000Z","shell.execute_reply":"2022-09-04T20:53:28.410689Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission.iloc[:len_cite_raveled] = np.nanmean(np.array(all_preds),axis=0)","metadata":{"execution":{"iopub.status.busy":"2022-09-04T20:53:28.415279Z","iopub.execute_input":"2022-09-04T20:53:28.415568Z","iopub.status.idle":"2022-09-04T20:53:29.010454Z","shell.execute_reply.started":"2022-09-04T20:53:28.415542Z","shell.execute_reply":"2022-09-04T20:53:29.009389Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission.to_csv('submission_cite.csv')\nsubmission.head()","metadata":{"execution":{"iopub.status.busy":"2022-09-04T20:53:29.012223Z","iopub.execute_input":"2022-09-04T20:53:29.012630Z","iopub.status.idle":"2022-09-04T20:56:22.784788Z","shell.execute_reply.started":"2022-09-04T20:53:29.012591Z","shell.execute_reply":"2022-09-04T20:56:22.783778Z"},"trusted":true},"execution_count":null,"outputs":[]}]}