{"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":"# Overview\n\n### Background\n\n- CITEseq samples take input X (RNA sequence vector) to predict output Y (Protein sequence vector)\n- Import generated features (n_samples, 240) shape created using PCA and simple feature engineering\n- Feed samples to encoder-decoder NN (see structure below)\n- Train one fold using pytorch (50 epochs, AdamW optimizer, and Cosine scheduler)\n- Use network to make predictions\n\n### Improving Upon this Notebook\n\n**The best way to improve upon this notebook is likely feature engineering**\n- Look for feature importance\n- Improve dimensionality reduction technique\n- Use domain knowledge\n- This is just a baseline\n\n**Apply similar model for Multiome samples** \n- Currently, I am just borrowing another submission for multiome and this notebook only predicts CITEseq\n- However, you can expand upon this notebook to make Multiome predictions\n\n**Change NN Structure**\n- test rnn or cnn\n- try adding attention mechanism\n- change structure to adjust for new features\n\n### Neural Network Structure\n![](https://www.googleapis.com/download/storage/v1/b/kaggle-forum-message-attachments/o/inbox%2F6537187%2Fcc192e980bc9f9248d2bcae2c6accd74%2Fencoder_decoder.PNG?generation=1662083041543650&alt=media)\n\n![](https://www.googleapis.com/download/storage/v1/b/kaggle-forum-message-attachments/o/inbox%2F6537187%2F1a17ab66143625efff11e8a063e1dac1%2Fenc_dec2.PNG?generation=1662083054477703&alt=media)\n\n![](https://www.googleapis.com/download/storage/v1/b/kaggle-forum-message-attachments/o/inbox%2F6537187%2Fb373fe534194a31dfc3505f90f488b25%2FFCBlock.PNG?generation=1662087114666150&alt=media)","metadata":{}},{"cell_type":"markdown","source":"# Read Data","metadata":{}},{"cell_type":"code","source":"import os, gc, pickle\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport numpy as np\nfrom colorama import Fore, Back, Style\nfrom matplotlib.ticker import MaxNLocator\nfrom tqdm import tqdm\n\nfrom sklearn.model_selection import KFold\nfrom sklearn.metrics import mean_squared_error\n\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nfrom torch.optim import AdamW\nfrom torch.optim import lr_scheduler\nfrom torch.utils.data import DataLoader, Dataset\n\ndevice = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\nprint(device)\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\")","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-09-04T13:43:20.438123Z","iopub.execute_input":"2022-09-04T13:43:20.439193Z","iopub.status.idle":"2022-09-04T13:43:23.324324Z","shell.execute_reply.started":"2022-09-04T13:43:20.439090Z","shell.execute_reply":"2022-09-04T13:43:23.323053Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install --quiet tables","metadata":{"execution":{"iopub.status.busy":"2022-09-04T13:43:23.326717Z","iopub.execute_input":"2022-09-04T13:43:23.327680Z","iopub.status.idle":"2022-09-04T13:43:36.912187Z","shell.execute_reply.started":"2022-09-04T13:43:23.327634Z","shell.execute_reply":"2022-09-04T13:43:36.911060Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cite_train_x = np.load(\"../input/citeseq-pca-240-preprocessing/cite_train_x.npy\")\nprint(cite_train_x.shape)\n\ncite_test_x = np.load(\"../input/citeseq-pca-240-preprocessing/cite_test_x.npy\")\nprint(cite_train_x.shape)\n\ncite_train_y = pd.read_hdf(FP_CITE_TRAIN_TARGETS).values\nprint(cite_train_y.shape)","metadata":{"execution":{"iopub.status.busy":"2022-09-04T13:43:36.913835Z","iopub.execute_input":"2022-09-04T13:43:36.914225Z","iopub.status.idle":"2022-09-04T13:43:38.991495Z","shell.execute_reply.started":"2022-09-04T13:43:36.914187Z","shell.execute_reply":"2022-09-04T13:43:38.989992Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Config","metadata":{}},{"cell_type":"code","source":"class CFG:\n    tr_batch_size = 16 # 16\n    va_batch_size = 128 # 32\n    \n    optimizer = \"AdamW\"\n    lr = 1e-5\n    weight_decay = 0.1\n    betas = (0.9, 0.999)\n    epochs = 50","metadata":{"execution":{"iopub.status.busy":"2022-09-04T13:43:38.994607Z","iopub.execute_input":"2022-09-04T13:43:38.995404Z","iopub.status.idle":"2022-09-04T13:43:39.002304Z","shell.execute_reply.started":"2022-09-04T13:43:38.995365Z","shell.execute_reply":"2022-09-04T13:43:39.000818Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Dataset","metadata":{}},{"cell_type":"code","source":"class CtieseqDataset(Dataset):\n    \"\"\"\n    Train, Validation or Test dataset for CITEseq samples\n    Prepares data for simple vector to vector NN\n    \"\"\"\n    def __init__(self, X, y=None):\n        self.train = False \n        if y is not None:\n            self.train = True\n        self.X = X\n        self.y = y\n            \n    def __len__(self):\n        return len(self.X)\n    \n    def __getitem__(self, idx):\n        X = self.X[idx]\n        \n        if self.train:\n            y = self.y[idx]\n            return {\n                \"X\" : torch.tensor(X).to(device),\n                \"y\" : torch.tensor(y).to(device)\n            }\n        else:\n            return {\n                \"X\" : torch.tensor(X).to(device)\n            }","metadata":{"execution":{"iopub.status.busy":"2022-09-04T13:43:39.003835Z","iopub.execute_input":"2022-09-04T13:43:39.004309Z","iopub.status.idle":"2022-09-04T13:43:39.015563Z","shell.execute_reply.started":"2022-09-04T13:43:39.004251Z","shell.execute_reply":"2022-09-04T13:43:39.014412Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Utils","metadata":{}},{"cell_type":"code","source":"def criterion(outputs, labels):\n    \"\"\" MSE Loss function\"\"\"\n    return nn.MSELoss()(outputs, labels)\n\ndef correlation_score(y_true, y_pred):\n    \"\"\"\n    Scores the predictions according to the competition rules. \n    It is assumed that the predictions are not constant.\n    Returns the average of each sample's Pearson correlation coefficient\n    \"\"\"\n    \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    for i in range(len(y_true)):\n        corrsum += np.corrcoef(y_true[i], y_pred[i])[1, 0]\n    return corrsum / len(y_true)\n\ndef get_optimizer(model, lr, weight_decay, betas):\n    \"\"\" Gets AdamW optimizer \"\"\"\n    param_optimizer = list(model.named_parameters())\n    no_decay = ['bias', 'LayerNorm.bias', 'LayerNorm.weight']\n    optimizer_grouped_parameters = [\n        {'params': [p for n, p in param_optimizer if not any(\n            nd in n for nd in no_decay)], 'weight_decay': weight_decay},\n        {'params': [p for n, p in param_optimizer if any(\n            nd in n for nd in no_decay)], 'weight_decay': 0.0}\n    ]\n    optimizer = AdamW(optimizer_grouped_parameters,\n                      lr=lr,\n                      weight_decay=weight_decay,\n                      betas=betas,\n                     )\n    return optimizer\n\ndef get_scheduler(optimizer, T_max=300):\n    \"\"\" Gets Consine scheduler \"\"\"\n    scheduler = torch.optim.lr_scheduler.CosineAnnealingLR(optimizer,\n                                                           T_max=T_max)\n    return scheduler","metadata":{"execution":{"iopub.status.busy":"2022-09-04T13:43:39.017283Z","iopub.execute_input":"2022-09-04T13:43:39.017747Z","iopub.status.idle":"2022-09-04T13:43:39.031412Z","shell.execute_reply.started":"2022-09-04T13:43:39.017706Z","shell.execute_reply":"2022-09-04T13:43:39.030439Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model","metadata":{}},{"cell_type":"code","source":"class FCBlock(nn.Module):\n    \"\"\"\n    A Pytorch Block for a fully connected Layer\n    Includes Linear, Activation Function, and Dropout\n    \"\"\"\n    def __init__(self, input_dim, hidden_dim, dropout):\n        super().__init__()\n        self.fc = nn.Linear(input_dim, hidden_dim)\n        self.dropout = nn.Dropout(dropout)\n        \n    def forward(self, x):\n        x = self.fc(x)\n        x = F.selu(x)\n        x = self.dropout(x)\n        return x\n\nclass Encoder(nn.Module):\n    \"\"\"\n    Encoder module to generate embeddings of a RNA vector\n    \"\"\"\n    def __init__(self):\n        super().__init__()\n        self.l0 = FCBlock(240, 120, 0.05)\n        self.l1 = FCBlock(120, 60, 0.05)\n        self.l2 = FCBlock(60, 30, 0.05)\n        \n    def forward(self, x):\n        x = self.l0(x)\n        x = self.l1(x)\n        x = self.l2(x)\n        return x\n    \nclass Decoder(nn.Module):\n    \"\"\"\n    Decoder module to extract Protein sequences from RNA embeddings\n    \"\"\"\n    def __init__(self):\n        super().__init__()\n        self.l0 = FCBlock(30, 70, 0.05)\n        self.l1 = FCBlock(70, 100, 0.05)\n        self.l2 = FCBlock(100, 140, 0.05)\n        \n    def forward(self, x):\n        x = self.l0(x)\n        x = self.l1(x)\n        x = self.l2(x)\n        return x\n    \nclass CtieseqModel(nn.Module):\n    \"\"\"\n    Wrapper for the Encoder and Decoder modules\n    Converts RNA sequence to Protein sequence\n    \"\"\"\n    def __init__(self):\n        super().__init__()\n        self.encoder = Encoder()\n        self.decoder = Decoder()\n        \n    def forward(self, x):\n        embeddings = self.encoder(x)\n        outputs = self.decoder(embeddings)\n        return outputs","metadata":{"execution":{"iopub.status.busy":"2022-09-04T13:43:54.133582Z","iopub.execute_input":"2022-09-04T13:43:54.133967Z","iopub.status.idle":"2022-09-04T13:43:54.148410Z","shell.execute_reply.started":"2022-09-04T13:43:54.133918Z","shell.execute_reply":"2022-09-04T13:43:54.147344Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train Loop Functions","metadata":{}},{"cell_type":"code","source":"class AverageMeter(object):\n    \"\"\"Computes and stores the average and current value\"\"\"\n    def __init__(self):\n        self.reset()\n    def reset(self):\n        self.val = 0\n        self.avg = 0\n        self.sum = 0\n        self.count = 0\n    def update(self, val, n=1):\n        self.val = val\n        self.sum += val * n\n        self.count += n\n        self.avg = self.sum / self.count","metadata":{"execution":{"iopub.status.busy":"2022-09-04T13:43:54.658081Z","iopub.execute_input":"2022-09-04T13:43:54.658473Z","iopub.status.idle":"2022-09-04T13:43:54.664907Z","shell.execute_reply.started":"2022-09-04T13:43:54.658440Z","shell.execute_reply":"2022-09-04T13:43:54.664046Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def train_one_epoch(model, optimizer, scheduler, dataloader, device, epoch):\n    \"\"\" Trains one epoch and returns loss \"\"\"\n    model.train()\n    \n    losses = AverageMeter()\n    corr = AverageMeter()\n    \n    bar = tqdm(enumerate(dataloader), total=len(dataloader))\n    for step, data in bar:\n        X, y = data[\"X\"], data[\"y\"]\n        \n        batch_size = X.size(0)\n\n        outputs = model(X)\n\n        n = outputs.size(0)\n        loss = criterion(outputs, y)\n        losses.update(loss.item(), n)\n        loss.backward()\n        \n        outputs = outputs.detach().cpu().numpy()\n        y = y.detach().cpu().numpy()\n        corr_score = correlation_score(y, outputs)\n        corr.update(corr_score, n)\n        \n        optimizer.step()\n        optimizer.zero_grad()\n\n        if scheduler is not None:\n            scheduler.step()\n        \n        bar.set_postfix(Epoch=epoch, Train_Loss=losses.avg, Corr=corr.avg,\n                        LR=optimizer.param_groups[0]['lr'])\n    gc.collect()\n    \n    return losses.avg","metadata":{"execution":{"iopub.status.busy":"2022-09-04T13:43:54.961959Z","iopub.execute_input":"2022-09-04T13:43:54.963228Z","iopub.status.idle":"2022-09-04T13:43:54.975047Z","shell.execute_reply.started":"2022-09-04T13:43:54.963177Z","shell.execute_reply":"2022-09-04T13:43:54.974059Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"@torch.no_grad()\ndef valid_one_epoch(model, optimizer, dataloader, device, epoch):\n    \"\"\" Evaluates one epoch and returns loss \"\"\"\n    model.eval()\n    \n    losses = AverageMeter()\n    corr = AverageMeter()\n    \n    bar = tqdm(enumerate(dataloader), total=len(dataloader))\n    for step, data in bar:        \n        X, y = data[\"X\"], data[\"y\"]\n        \n        batch_size = X.size(0)\n\n        outputs = model(X)\n        \n        n = outputs.size(0)\n        loss = criterion(outputs, y)\n        losses.update(loss.item(), n)\n        \n        outputs = outputs.detach().cpu().numpy()\n        y = y.detach().cpu().numpy()\n        corr_score = correlation_score(y, outputs)\n        corr.update(corr_score, n)\n        \n        bar.set_postfix(Epoch=epoch, Valid_Loss=losses.avg, Corr=corr.avg,\n                        LR=optimizer.param_groups[0]['lr'])   \n    \n    gc.collect()\n    \n    return losses.avg","metadata":{"execution":{"iopub.status.busy":"2022-09-04T13:43:55.337183Z","iopub.execute_input":"2022-09-04T13:43:55.337793Z","iopub.status.idle":"2022-09-04T13:43:55.346412Z","shell.execute_reply.started":"2022-09-04T13:43:55.337760Z","shell.execute_reply":"2022-09-04T13:43:55.345327Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def train_one_fold(model, \n                   optimizer, \n                   scheduler, \n                   train_loader, \n                   valid_loader, \n                   fold):\n    \"\"\" Trains and saves a full fold of a pytorch model \"\"\"\n    best_epoch_loss = np.inf\n    model.to(device)\n\n    for epoch in range(CFG.epochs):\n        gc.collect()\n        train_epoch_loss = train_one_epoch(model, \n                                           optimizer, \n                                           scheduler, \n                                           dataloader=train_loader, \n                                           device=device, \n                                           epoch=epoch)\n\n        val_epoch_loss = valid_one_epoch(model,\n                                         optimizer, \n                                         valid_loader, \n                                         device=device, epoch=epoch)\n        \n        if val_epoch_loss <= best_epoch_loss:\n            print(f\"Validation Loss Improved ({best_epoch_loss} ---> {val_epoch_loss})\")\n            best_epoch_loss = val_epoch_loss\n            torch.save(model.state_dict(), f\"model_f{fold}.bin\")\n            \n    print(\"Best Loss: {:.4f}\".format(best_epoch_loss))","metadata":{"execution":{"iopub.status.busy":"2022-09-04T13:43:55.968736Z","iopub.execute_input":"2022-09-04T13:43:55.969839Z","iopub.status.idle":"2022-09-04T13:43:55.979036Z","shell.execute_reply.started":"2022-09-04T13:43:55.969789Z","shell.execute_reply":"2022-09-04T13:43:55.977981Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Run Training","metadata":{}},{"cell_type":"code","source":"kf = KFold(n_splits=5, shuffle=True, random_state=42)\nscore_list = []\nfor fold, (idx_tr, idx_va) in enumerate(kf.split(cite_train_x)):\n    print(f\"\\nfold = {fold}\")\n    X_tr = cite_train_x[idx_tr] \n    y_tr = cite_train_y[idx_tr]\n    \n    X_va = cite_train_x[idx_va]\n    y_va = cite_train_y[idx_va]\n    \n    ds_tr = CtieseqDataset(X_tr, y_tr)\n    ds_va = CtieseqDataset(X_tr, y_tr)\n    dl_tr = DataLoader(ds_tr, batch_size=CFG.tr_batch_size, shuffle=True)\n    dl_va = DataLoader(ds_va, batch_size=CFG.va_batch_size, shuffle=False)\n    \n    model = CtieseqModel()\n    optimizer = get_optimizer(model, CFG.lr, CFG.weight_decay, CFG.betas)\n    scheduler = get_scheduler(optimizer) \n\n    train_one_fold(model, optimizer, scheduler, dl_tr, dl_va, fold)\n    break","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-09-04T13:53:33.201452Z","iopub.execute_input":"2022-09-04T13:53:33.201851Z","iopub.status.idle":"2022-09-04T13:54:48.654003Z","shell.execute_reply.started":"2022-09-04T13:53:33.201815Z","shell.execute_reply":"2022-09-04T13:54:48.653000Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Predict & Submit","metadata":{}},{"cell_type":"code","source":"def predict(fold):\n    preds = list()\n    ds = CtieseqDataset(cite_test_x)\n    dl = DataLoader(ds, batch_size=32, shuffle=False)\n    model = CtieseqModel()\n    model.load_state_dict(torch.load(f\"model_f{fold}.bin\"))\n    model.eval()\n\n    bar = tqdm(enumerate(dl), total=len(dl))\n    for step, data in bar:        \n        X = data[\"X\"]\n\n        batch_size = X.size(0)\n\n        outputs = model(X)\n        preds.append(outputs.detach().cpu().numpy())\n    test_pred = np.concatenate(preds)\n    return test_pred\n\n#test_preds = np.array([predict(0), predict(1), predict(2), predict(3), predict(4)])\ntest_preds = np.array([predict(0)])\ntest_preds = np.mean(test_preds, axis=0)\nprint(test_preds.shape)","metadata":{"execution":{"iopub.status.busy":"2022-09-04T13:49:17.966343Z","iopub.execute_input":"2022-09-04T13:49:17.966749Z","iopub.status.idle":"2022-09-04T13:49:23.447583Z","shell.execute_reply.started":"2022-09-04T13:49:17.966714Z","shell.execute_reply":"2022-09-04T13:49:23.446520Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Submission\nsubmission = pd.read_csv('../input/msci-multiome-quickstart-w-sparse-matrices/submission.csv', index_col='row_id', squeeze=True)\nsubmission.iloc[:len(test_preds.ravel())] = test_preds.ravel()\nassert not submission.isna().any()\nsubmission = submission.round(6) # reduce the size of the csv\nsubmission.to_csv('submission.csv')\nsubmission","metadata":{"execution":{"iopub.status.busy":"2022-09-04T13:49:41.248633Z","iopub.execute_input":"2022-09-04T13:49:41.249023Z","iopub.status.idle":"2022-09-04T13:52:49.940775Z","shell.execute_reply.started":"2022-09-04T13:49:41.248991Z","shell.execute_reply":"2022-09-04T13:52:49.939483Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}