{"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":"# Multiome Time Aspect 2\n\nIn this notebook we look at several Multiome model variants for predicting day 7:\n- Baseline: training on days 2,3,4 (without using the day as a feature)\n- Training on days 2,3,4 with day 4 oversampled three times (without using the day as a feature)\n- Training on days 2,3,4 with day 4 oversampled five times (without using the day as a feature)\n- Training only on day 4 (without using the day as a feature)\n- Training on days 2,3,4, using the day as a feature\n- Training on days 2,3,4, using the day (clipped to (2, 4.5)) as a feature\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"import os, gc, pickle, datetime, scipy.sparse, math, matplotlib\nimport seaborn as sns\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom colorama import Fore, Back, Style\n\nfrom sklearn.model_selection import GroupKFold, ShuffleSplit, train_test_split\nfrom sklearn.preprocessing import StandardScaler, scale, MinMaxScaler\nfrom sklearn.decomposition import TruncatedSVD\n\nimport tensorflow as tf\nimport tensorflow.keras.backend as K\nfrom tensorflow.keras.models import Model, load_model\nfrom tensorflow.keras.callbacks import ReduceLROnPlateau, LearningRateScheduler, EarlyStopping\nfrom tensorflow.keras.layers import Dense, Input, Concatenate, Dropout, BatchNormalization\nfrom tensorflow.keras.utils import plot_model\nimport keras_tuner\n\nDATA_DIR = \"/kaggle/input/open-problems-multimodal/\"\nFP_CELL_METADATA = os.path.join(DATA_DIR,\"metadata.csv\")\n\nFP_CITE_TRAIN_INPUTS = os.path.join(DATA_DIR,\"train_cite_inputs.h5\")\nFP_CITE_TRAIN_TARGETS = os.path.join(DATA_DIR,\"train_cite_targets.h5\")\nFP_CITE_TEST_INPUTS = os.path.join(DATA_DIR,\"test_cite_inputs.h5\")\n\nFP_MULTIOME_TRAIN_INPUTS = os.path.join(DATA_DIR,\"train_multi_inputs.h5\")\nFP_MULTIOME_TRAIN_TARGETS = os.path.join(DATA_DIR,\"train_multi_targets.h5\")\nFP_MULTIOME_TEST_INPUTS = os.path.join(DATA_DIR,\"test_multi_inputs.h5\")\n\nFP_SUBMISSION = os.path.join(DATA_DIR,\"sample_submission.csv\")\nFP_EVALUATION_IDS = os.path.join(DATA_DIR,\"evaluation_ids.csv\")\n\nCV_RUNS = 20\n\nsummary_list = []\nnp.set_printoptions(linewidth=150, edgeitems=5)\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-02T15:40:58.043302Z","iopub.execute_input":"2022-12-02T15:40:58.043657Z","iopub.status.idle":"2022-12-02T15:40:58.055212Z","shell.execute_reply.started":"2022-12-02T15:40:58.043627Z","shell.execute_reply":"2022-12-02T15:40:58.054264Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# The scoring function\n\nThis competition has a special metric: For every row, it computes the Pearson correlation between y_true and y_pred, and then all these correlation coefficients are averaged.","metadata":{}},{"cell_type":"code","source":"def correlation_score(y_true, y_pred):\n    \"\"\"Scores the predictions according to the competition rules. \n    \n    It is assumed that the predictions are not constant.\n    \n    Parameters:\n    y_true: (sparse) matrix of shape (n_samples, n_outputs)\n    y_pred: array of shape (n_samples, n_outputs)\n    \n    Returns the average of each sample's Pearson correlation coefficient\"\"\"\n    if type(y_true) == pd.DataFrame: y_true = y_true.values\n    if type(y_pred) == pd.DataFrame: y_pred = y_pred.values\n    corrsum = 0\n    if type(y_true) == np.ndarray:\n        for i in range(y_true.shape[0]):\n            corrsum += np.corrcoef(y_true[i].ravel(), y_pred[i])[1, 0]\n    else:\n        for i in range(y_true.shape[0]):\n            corrsum += np.corrcoef(y_true[i].toarray().ravel(), y_pred[i])[1, 0]\n    return corrsum / y_true.shape[0]\n","metadata":{"execution":{"iopub.status.busy":"2022-12-02T15:33:42.268181Z","iopub.execute_input":"2022-12-02T15:33:42.268815Z","iopub.status.idle":"2022-12-02T15:33:42.275840Z","shell.execute_reply.started":"2022-12-02T15:33:42.268777Z","shell.execute_reply":"2022-12-02T15:33:42.274867Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Training data loading and preprocessing\n\nThe metadata is used only for the `GroupKFold`: ","metadata":{}},{"cell_type":"code","source":"metadata_df = pd.read_csv(FP_CELL_METADATA, index_col='cell_id')\nmetadata_df = metadata_df[metadata_df.technology == \"multiome\"]\nmetadata_df.shape\n\n# Load the training cell_ids and reindex the metadata to prepare the GroupKFold\nwith np.load(\"../input/multimodal-single-cell-as-sparse-matrix/train_multi_inputs_idxcol.npz\", allow_pickle=True) as idx:\n    train_index = idx['index']\n    train_columns = idx['columns']\nmeta = metadata_df.reindex(train_index)","metadata":{"execution":{"iopub.status.busy":"2022-12-02T15:33:42.277741Z","iopub.execute_input":"2022-12-02T15:33:42.278442Z","iopub.status.idle":"2022-12-02T15:33:43.139810Z","shell.execute_reply.started":"2022-12-02T15:33:42.278407Z","shell.execute_reply":"2022-12-02T15:33:43.138885Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nn_input_components = 96\nwith open(f'../input/msci-m-preprocessing/X_complete_256.pickle', 'rb') as f: \n    X = pickle.load(f)\nX = X[:,:n_input_components].copy()\nprint(f\"X shape after SVD: {X.shape}\")\n","metadata":{"execution":{"iopub.status.busy":"2022-12-02T15:33:43.142201Z","iopub.execute_input":"2022-12-02T15:33:43.142545Z","iopub.status.idle":"2022-12-02T15:33:44.337876Z","shell.execute_reply.started":"2022-12-02T15:33:43.142510Z","shell.execute_reply":"2022-12-02T15:33:44.336603Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nY = scipy.sparse.load_npz(\"../input/multimodal-single-cell-as-sparse-matrix/train_multi_targets_values.sparse.npz\")\nprint(f\"Y shape before SVD: {Y.shape}\")\n","metadata":{"execution":{"iopub.status.busy":"2022-12-02T15:33:44.339479Z","iopub.execute_input":"2022-12-02T15:33:44.339944Z","iopub.status.idle":"2022-12-02T15:34:08.032743Z","shell.execute_reply.started":"2022-12-02T15:33:44.339907Z","shell.execute_reply":"2022-12-02T15:34:08.031662Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nn_output_components = 80\n\ndef load_output_svd(title):\n    with open(f'../input/msci-m-preprocessing/osvd_256_{title}.pickle', 'rb') as f: \n        output_svd = pickle.load(f)\n    output_svd.components_ = output_svd.components_[:n_output_components]\n    output_svd.explained_variance_ = output_svd.explained_variance_[:n_output_components]\n    output_svd.explained_variance_ratio_ = output_svd.explained_variance_ratio_[:n_output_components]\n    output_svd.singular_values_ = output_svd.singular_values_[:n_output_components]\n    with open(f'../input/msci-m-preprocessing/Y_256_{title}.pickle', 'rb') as f: \n        Yp = pickle.load(f)[:,:n_output_components]\n    return output_svd, Yp\n\noutput_svd, Yp = load_output_svd(f\"full\")\n","metadata":{"execution":{"iopub.status.busy":"2022-12-02T15:34:08.034145Z","iopub.execute_input":"2022-12-02T15:34:08.035212Z","iopub.status.idle":"2022-12-02T15:34:09.151487Z","shell.execute_reply.started":"2022-12-02T15:34:08.035164Z","shell.execute_reply":"2022-12-02T15:34:09.150365Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# The model\n\nOur model is a sequential network consisting of a few dense layers. The hyperparameters will be tuned with KerasTuner.\n","metadata":{}},{"cell_type":"code","source":"best_hp = keras_tuner.HyperParameters()\nbest_hp.values = {'reg1': 0.035538207011211874, 'reg3': 1.5375567281431701e-09, 'units1': 64, 'units2': 64, 'units3': 256, 'units4': 256}\n\nLR_START = 0.01\nBATCH_SIZE = 512 # 256\n\ndef my_model(hp, n_inputs=X.shape[1]):\n    \"\"\"Sequential neural network\n    \n    Returns a compiled instance of tensorflow.keras.models.Model.\n    \"\"\"\n    activation = 'relu'\n    reg1 = hp.Float(\"reg1\", min_value=1e-8, max_value=1e1, sampling=\"log\")\n    reg3 = hp.Float(\"reg3\", min_value=1e-14, max_value=1e-5, sampling=\"log\")\n    \n    inputs = Input(shape=(n_inputs, ))\n    x0 = Dense(hp.Choice('units1', [16, 32, 64]), kernel_regularizer=tf.keras.regularizers.l2(reg1),\n              activation=activation,\n             )(inputs)\n    x1 = Dense(hp.Choice('units2', [16, 32, 64]),\n              activation=activation,\n             )(x0)\n    x1 = BatchNormalization()(x1)\n    x2 = Dense(hp.Choice('units3', [32, 64, 128, 256]),\n              activation=activation,\n             )(x1)\n    x3 = Dense(hp.Choice('units4', [32, 64, 128, 256]),\n              activation=activation,\n             )(x2)\n    x = Concatenate()([x1, x2, x3])\n    x = Dense(n_output_components, kernel_regularizer=tf.keras.regularizers.l2(reg3),\n              #activation=activation,\n             )(x)\n    regressor = Model(inputs, x)\n    regressor.compile(optimizer=tf.keras.optimizers.Adam(learning_rate=LR_START),\n                      metrics=['MSE'],\n                      loss='MSE',\n                      steps_per_execution=300,\n                     )\n    \n    return regressor\n\ndisplay(plot_model(my_model(best_hp),\n                   show_layer_names=False, show_shapes=True, dpi=72))\n","metadata":{"execution":{"iopub.status.busy":"2022-12-02T15:34:09.153114Z","iopub.execute_input":"2022-12-02T15:34:09.153491Z","iopub.status.idle":"2022-12-02T15:34:12.713013Z","shell.execute_reply.started":"2022-12-02T15:34:09.153455Z","shell.execute_reply":"2022-12-02T15:34:12.711737Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Validation\n","metadata":{}},{"cell_type":"code","source":"%%time\n# Cross-validation\nVERBOSE = 0 # set to 2 for more output, set to 0 for less output\nEPOCHS = 9000\n\nnp.random.seed(1)\ntf.random.set_seed(1)\n\nscore_list = []\nfor i, variant in enumerate(['all_days',\n                             'day_4_oversampled_3times', 'day_4_oversampled_5times', 'day_4_only',\n                             'day_feature', 'day_feature_clipped']):\n    if variant == 'all_days':\n        idx_tr = np.arange(len(X))[meta.day != 7]\n        sample_weight = None\n        XX = X\n    elif variant == 'day_4_oversampled_3times':\n        idx_tr = np.arange(len(X))[meta.day != 7]\n        sample_weight = np.ones(len(idx_tr))\n        sample_weight[meta.day[idx_tr] == 4] = 3\n        XX = X\n    elif variant == 'day_4_oversampled_5times':\n        idx_tr = np.arange(len(X))[meta.day != 7]\n        sample_weight = np.ones(len(idx_tr))\n        sample_weight[meta.day[idx_tr] == 4] = 5\n        XX = X\n    elif variant == 'day_4_only':\n        idx_tr = np.arange(len(X))[meta.day == 4]\n        sample_weight = None\n        XX = X\n    elif variant == 'day_feature':\n        idx_tr = np.arange(len(X))[meta.day != 7]\n        sample_weight = np.ones(len(idx_tr))\n        sample_weight[meta.day[idx_tr] == 4] = 5\n        XX = np.hstack([X, meta[['day']].values])\n    elif variant == 'day_feature_clipped':\n        idx_tr = np.arange(len(X))[meta.day != 7]\n        sample_weight = np.ones(len(idx_tr))\n        sample_weight[meta.day[idx_tr] == 4] = 5\n        XX = np.hstack([X, meta[['day']].values.clip(2, 4.5)])\n        \n    idx_va = np.arange(len(X))[meta.day == 7]\n    start_time = datetime.datetime.now()\n    model = None\n    gc.collect()\n    X_tr = XX[idx_tr]\n    X_va = XX[idx_va]\n    y_va_pred = np.zeros((len(idx_va), Y.shape[1]), dtype=np.float32)\n        \n    for run in range(CV_RUNS):\n        # Construct and compile the model\n        model = my_model(best_hp, X_tr.shape[1])\n\n        lr = ReduceLROnPlateau(monitor=\"val_loss\", factor=0.5, \n                               patience=4, verbose=VERBOSE)\n        es = EarlyStopping(monitor=\"val_loss\",\n                           patience=12, \n                           verbose=0,\n                           mode=\"min\", \n                           restore_best_weights=True)\n        callbacks = [lr, es, tf.keras.callbacks.TerminateOnNaN()]\n\n        # Train the model\n        history = model.fit(X_tr, Yp[idx_tr],\n                            sample_weight=sample_weight,\n                            validation_data=(X_va, Yp[idx_va]), \n                            validation_batch_size=len(idx_va),\n                            epochs=EPOCHS,\n                            verbose=VERBOSE,\n                            batch_size=BATCH_SIZE,\n                            shuffle=True,\n                            callbacks=callbacks)\n\n        history = history.history\n        callbacks, lr = None, None\n\n        # Validate the model\n        y_va_pred_1 = model.predict(X_va, batch_size=len(X_va))\n        y_va_pred_1 = output_svd.inverse_transform(y_va_pred_1)\n        corrscore = correlation_score(Y[idx_va], y_va_pred_1)\n        print(f\"Variant {variant}.{run}: {es.stopped_epoch:3} epochs, corr =  {corrscore:.5f}\")    \n        y_va_pred += y_va_pred_1\n        y_va_pred_1, model, es = None, None, None\n        \n    X_tr, X_va = None, None\n    corrscore = correlation_score(Y[idx_va], y_va_pred)\n    print(f\"Variant {variant+':':30}                                corr =  {corrscore:.5f}\")\n    score_list.append((variant, corrscore))\n","metadata":{"execution":{"iopub.status.busy":"2022-12-02T15:42:53.120349Z","iopub.execute_input":"2022-12-02T15:42:53.120719Z","iopub.status.idle":"2022-12-02T15:52:05.440848Z","shell.execute_reply.started":"2022-12-02T15:42:53.120688Z","shell.execute_reply":"2022-12-02T15:52:05.439666Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"result_df = pd.DataFrame(score_list, columns=['variant', 'corrscore'])\nmatplotlib.rcParams['font.size'] = 16\nplt.figure(figsize=(10, 4))\nplt.barh(np.arange(len(result_df)), result_df.corrscore, color='lightgreen')\nplt.yticks(np.arange(len(result_df)), result_df.variant)\nplt.xlim(0.604, 0.612)\nplt.gca().invert_yaxis()\nplt.xlabel('Correlation score (higher is better)')\nplt.title('Day 7 predictions when trained on ...')\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-02T15:52:08.754110Z","iopub.execute_input":"2022-12-02T15:52:08.755199Z","iopub.status.idle":"2022-12-02T15:52:08.968748Z","shell.execute_reply.started":"2022-12-02T15:52:08.755147Z","shell.execute_reply":"2022-12-02T15:52:08.967686Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"score_list","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-12-02T15:52:21.373729Z","iopub.execute_input":"2022-12-02T15:52:21.374754Z","iopub.status.idle":"2022-12-02T15:52:21.383296Z","shell.execute_reply.started":"2022-12-02T15:52:21.374716Z","shell.execute_reply":"2022-12-02T15:52:21.382268Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}