{"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":"<div class=\"alert alert-block alert-success\" style=\"font-size:30px\">\n[LB:0.811] Normalized Ensembles for Pearson's Correlation Score Function\n</div>\n\n\n<div class=\"alert alert-block alert-danger\" style=\"text-align:center; font-size:20px;\">\n    ❤️ Dont forget to ▲upvote▲ if you find this notebook usefull!  ❤️\n</div>\n\n\nIn this notebook I want to share a more robust ensembling method that can be used for both final multimodel ensembles and for ensembling of the same model trained on multiple folds. MSCI competition uses **correlation coefficient** as a scoring function, which affects the choice of ensembling logic. \n\n<span style='font-size:18px'>TL;DR: standardize your outputs per cell_id before adding base submissions to the ensemble!</span>\n\n\n\n\n\nLets explore the following claims:\n\n## Statement 1. Correlation loss is insensitive to linear transformations of predictions\n\nYou can have 2 solutions with the same score but with vastly different predictions for any given cell_id and gene_id. And this is because of this property. Quick proof. Let $X$ be our solution and $Y$ are ground truth labels. First recall that:\n$$\ncorr(X, Y) = \\frac {cov(X, Y)} {\\sigma_X * \\sigma_Y} = \\frac {\\mathbb{E}[(X-\\mathbb{E}(X))*(Y-\\mathbb{E}(Y))]} {\\sigma_X * \\sigma_Y}\n$$\nNow take another solution $X'=C_1+C_2*X$\n\n$$\ncorr(X', Y)=corr(C_1+C_2*X, Y) = \\frac {\\mathbb{E}[(C_1+C_2*X-\\mathbb{E}[C_1+C_2*X])*(Y-\\mathbb{E}(Y))]} {\\sigma_{[C_1+C_2*X]} * \\sigma_Y}\\\\\n=\\frac {\\mathbb{E}[(C_1+C_2*X-C_1+C_2*\\mathbb{E}(X))*(Y-\\mathbb{E}(Y))]} {C_2*\\sigma_{X} * \\sigma_Y}\\\\\n=\\frac {C_2*\\mathbb{E}[(X-\\mathbb{E}(X))*(Y-\\mathbb{E}(Y))]} {C_2*\\sigma_{X} * \\sigma_Y}\\\\\n=\\frac {\\mathbb{E}[(X-\\mathbb{E}(X))*(Y-\\mathbb{E}(Y))]} {\\sigma_{X} * \\sigma_Y}\\\\\n=corr(X, Y)\n$$\nSo we see that multiplying by $C_2$ and adding $C_1$ doesn't affect the score. \n\nFrom practical standpoint this means that we could have 2 similar solutions which we want to weight with coefficients $w_1$ and $w_2$, but the difference in magnitude of these solutions could be huge (e.g. $C_2=123$) which would make correct weighting impossible. This is an unlikely scenario if MSE metric was used to train base models, but it's totally possible if models were optimized directly with correlation score loss function!\n\n## Statement 2. Per-cell_id standardization helps to rescale base submissions\nUnder assumption that two base submissions are similar and demonstrate similar performance we could rescale them in the way that they become comparable and weighting in a regular way becomes adequate:\n\n$$\nX'=\\frac {X-\\mathbb{E}X} {\\sigma_X}\n$$\n\n## Statement 3. Weighting coefficients don't have to add up to 1!\nThis is one of the benefit of the loss function that is agnostic to linear transformations. You don't have to weight base submissions as usual with $\\sum_i w_i=1$. Any coefficients will do the job!\n\n\n## Statement 4. Only collect predictions for one of the technologies (CITEseq, Multiome) from every base solution \nThis is another hack unrelated to the correlation score function.\nMost of public notebooks build models for a single technology (CITEseq or Multiome) and paste the rest of predictions from the best availble public notebook for the other technology. \nThis results in less control you have over base predictions. E.g. you might end up with a good ensemble for CITEseq, but all Multiome predictions would actually come from a single source notebook which is suboptimal!\n\nIn this notebook we carefully pick up only relevant predictions from every base submission. ","metadata":{}},{"cell_type":"markdown","source":"<div class=\"alert alert-block alert-success\" style=\"font-size:30px\">\nA toy example demonstrating benefits of normalization\n</div>","metadata":{}},{"cell_type":"code","source":"import numpy as np \nimport pandas as pd \nimport glob\nfrom tqdm.notebook import tqdm\nimport os","metadata":{"execution":{"iopub.status.busy":"2022-09-09T06:15:23.951143Z","iopub.execute_input":"2022-09-09T06:15:23.951520Z","iopub.status.idle":"2022-09-09T06:15:23.956494Z","shell.execute_reply.started":"2022-09-09T06:15:23.951489Z","shell.execute_reply":"2022-09-09T06:15:23.955433Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# our groud truth targets\ntargets = np.random.randn(100000)","metadata":{"execution":{"iopub.status.busy":"2022-09-09T06:15:23.959550Z","iopub.execute_input":"2022-09-09T06:15:23.960328Z","iopub.status.idle":"2022-09-09T06:15:23.970177Z","shell.execute_reply.started":"2022-09-09T06:15:23.960249Z","shell.execute_reply":"2022-09-09T06:15:23.969428Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# submission1 = targets + some random noise\nsubmission1 = targets + 0.5 * np.random.randn(100000)\nsubmission1 ","metadata":{"execution":{"iopub.status.busy":"2022-09-09T06:15:23.971484Z","iopub.execute_input":"2022-09-09T06:15:23.971976Z","iopub.status.idle":"2022-09-09T06:15:23.984685Z","shell.execute_reply.started":"2022-09-09T06:15:23.971950Z","shell.execute_reply":"2022-09-09T06:15:23.983498Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# submission2 = targets + same amount of random noise + linear transformation\nsubmission2 = 4 * (targets + 0.5 * np.random.randn(100000))\nsubmission2 ","metadata":{"execution":{"iopub.status.busy":"2022-09-09T06:15:23.986482Z","iopub.execute_input":"2022-09-09T06:15:23.986781Z","iopub.status.idle":"2022-09-09T06:15:23.995208Z","shell.execute_reply.started":"2022-09-09T06:15:23.986755Z","shell.execute_reply":"2022-09-09T06:15:23.994379Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# correlation with target of submission1 and submission2 is quite similar\nnp.corrcoef(submission1, targets)","metadata":{"execution":{"iopub.status.busy":"2022-09-09T06:15:23.997427Z","iopub.execute_input":"2022-09-09T06:15:23.997750Z","iopub.status.idle":"2022-09-09T06:15:24.008031Z","shell.execute_reply.started":"2022-09-09T06:15:23.997722Z","shell.execute_reply":"2022-09-09T06:15:24.006577Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.corrcoef(submission2, targets)","metadata":{"execution":{"iopub.status.busy":"2022-09-09T06:15:24.009253Z","iopub.execute_input":"2022-09-09T06:15:24.010442Z","iopub.status.idle":"2022-09-09T06:15:24.018626Z","shell.execute_reply.started":"2022-09-09T06:15:24.010383Z","shell.execute_reply":"2022-09-09T06:15:24.017168Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Let's evaluate the standard average ensemble\nnp.corrcoef((submission1 + submission2) / 2, targets)","metadata":{"execution":{"iopub.status.busy":"2022-09-09T06:15:24.019863Z","iopub.execute_input":"2022-09-09T06:15:24.020413Z","iopub.status.idle":"2022-09-09T06:15:24.030734Z","shell.execute_reply.started":"2022-09-09T06:15:24.020373Z","shell.execute_reply":"2022-09-09T06:15:24.029624Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now let's standardise submissions first","metadata":{}},{"cell_type":"code","source":"# Let's standardise first. You can see the gain of 0.92 -> 0.94 after applying normalization trick!\n\ndef std(x):\n    return (x - np.mean(x)) / np.std(x)\n\nnp.corrcoef(std(submission1) + std(submission2), targets)","metadata":{"execution":{"iopub.status.busy":"2022-09-09T06:15:24.032250Z","iopub.execute_input":"2022-09-09T06:15:24.032816Z","iopub.status.idle":"2022-09-09T06:15:24.043771Z","shell.execute_reply.started":"2022-09-09T06:15:24.032776Z","shell.execute_reply":"2022-09-09T06:15:24.042685Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<div class=\"alert alert-block alert-success\" style=\"font-size:30px\">\nExample of ensembling with rescaling of base solutions\n</div>","metadata":{}},{"cell_type":"code","source":"SUBMISSIONS = {\n    \n    # LB: 0.81 https://www.kaggle.com/code/xiafire/lb-t15-msci-multiome-catboostregressor\n    '../input/lb-t15-msci-multiome-catboostregressor/submission.csv': 1.,\n    \n    # LB 0.81 https://www.kaggle.com/code/sskknt/msci-citeseq-keras-quickstart-dropout\n    '../input/msci-citeseq-keras-quickstart-dropout/submission.csv': 1.,         \n    \n    # LB: 0.809 https://www.kaggle.com/code/ambrosm/msci-citeseq-keras-quickstart\n    '../input/msci-citeseq-keras-quickstart/submission.csv': 0.7,\n        \n    # LB: 0.808 https://www.kaggle.com/code/fabiencrom/msci-multiome-torch-quickstart-submission\n    '../input/msci-multiome-torch-quickstart-submission/submission.csv': 0.5,\n    \n    # LB: 0.804 https://www.kaggle.com/code/xiafire/fork-of-msci-multiome-randomsampling-sp-6b182b\n    '../input/fork-of-msci-multiome-randomsampling-sp-6b182b/submission.csv': 0.3,\n        \n    # LB: 0.803 - https://www.kaggle.com/code/jsmithperera/multiome-quickstart-w-sparse-m-tsvd-32\n    '../input/multiome-quickstart-w-sparse-m-tsvd-32/submission.csv': 0.2,\n        \n    # LB: 0.803 - https://www.kaggle.com/code/fabiencrom/msci-multiome-quickstart-w-sparse-matrices\n    '../input/msci-multiome-quickstart-w-sparse-matrices/submission.csv': 0.2,            \n    \n    # LB: 0.803 - https://www.kaggle.com/code/ambrosm/msci-citeseq-quickstart/notebook\n    '../input/msci-citeseq-quickstart/submission.csv': 0.2,\n        \n    \n    # 0.797 - https://www.kaggle.com/code/ravishah1/citeseq-rna-to-protein-encoder-decoder-nn\n    #'../input/citeseq-rna-to-protein-encoder-decoder-nn/submission.csv': 0.5,\n        \n    # LB: 0.792 - https://www.kaggle.com/code/swimmy/lgbm-baseline-msci-citeseq\n    #'../input/lgbm-baseline-msci-citeseq/submission.csv': 0.2\n}","metadata":{"execution":{"iopub.status.busy":"2022-09-09T06:15:24.046756Z","iopub.execute_input":"2022-09-09T06:15:24.047400Z","iopub.status.idle":"2022-09-09T06:15:24.052416Z","shell.execute_reply.started":"2022-09-09T06:15:24.047335Z","shell.execute_reply":"2022-09-09T06:15:24.051734Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cell_ids = pd.read_parquet('../input/multimodal-single-cell-as-sparse-matrix/evaluation.parquet').cell_id","metadata":{"execution":{"iopub.status.busy":"2022-09-09T06:15:24.054922Z","iopub.execute_input":"2022-09-09T06:15:24.055181Z","iopub.status.idle":"2022-09-09T06:15:36.460710Z","shell.execute_reply.started":"2022-09-09T06:15:24.055156Z","shell.execute_reply":"2022-09-09T06:15:36.459801Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def gen_std_submission(path, cell_ids):\n    \"\"\"\n    Standardize submission per cell_id\n    \"\"\"\n    df = pd.read_csv(path)\n    df['cell_id'] = cell_ids    \n    vals = []\n    for idx, g in tqdm(df.groupby('cell_id', sort=False), desc=f'Standardizing {path}', miniters=1000):\n        vals.append(std(g.target).values)\n    vals = np.concatenate(vals)\n    return vals","metadata":{"execution":{"iopub.status.busy":"2022-09-09T06:15:36.461565Z","iopub.execute_input":"2022-09-09T06:15:36.461828Z","iopub.status.idle":"2022-09-09T06:15:36.467533Z","shell.execute_reply.started":"2022-09-09T06:15:36.461801Z","shell.execute_reply":"2022-09-09T06:15:36.466844Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def gen_ensemble(technology):\n    ensemble = None\n    for path in tqdm([path for path in SUBMISSIONS.keys() if technology in path], desc='Process submission'):\n        weight = SUBMISSIONS[path]\n        if ensemble is None:\n            ensemble = gen_std_submission(path, cell_ids) * weight\n        else:\n            ensemble += gen_std_submission(path, cell_ids) * weight\n    return ensemble","metadata":{"execution":{"iopub.status.busy":"2022-09-09T06:15:36.468617Z","iopub.execute_input":"2022-09-09T06:15:36.468900Z","iopub.status.idle":"2022-09-09T06:15:36.482543Z","shell.execute_reply.started":"2022-09-09T06:15:36.468866Z","shell.execute_reply":"2022-09-09T06:15:36.481758Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"PRED_SEGMENTS = [(0, 6812820), (6812820, 65744180)]\nensemble = []\nfor tech, (from_idx, to_idx) in tqdm(list(zip(['citeseq', 'multiome'], PRED_SEGMENTS)), desc='Technology'):    \n    ensemble.append(gen_ensemble(tech)[from_idx: to_idx])\n    \n    \nensemble = np.concatenate(ensemble)","metadata":{"execution":{"iopub.status.busy":"2022-09-09T06:15:36.483774Z","iopub.execute_input":"2022-09-09T06:15:36.484210Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_submit = pd.read_parquet('../input/multimodal-single-cell-as-sparse-matrix/sample_submission.parquet')\ndf_submit['target'] = ensemble\ndf_submit.to_csv('submission.csv', index=False)\ndf_submit","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<div class=\"alert alert-block alert-danger\" style=\"text-align:center; font-size:20px;\">\n    ❤️ Dont forget to ▲upvote▲ if you find this notebook usefull!  ❤️\n</div>","metadata":{"execution":{"iopub.status.busy":"2022-09-04T18:45:04.379501Z","iopub.execute_input":"2022-09-04T18:45:04.379886Z","iopub.status.idle":"2022-09-04T18:45:04.386762Z","shell.execute_reply.started":"2022-09-04T18:45:04.379855Z","shell.execute_reply":"2022-09-04T18:45:04.385279Z"}}}]}