{"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":"# What is about ?\n\n**Briefly:** Exercises with different blending strategies on different types of the models.\nBased on solutions for the Kaggle competition: \"Open Problems - Multimodal Single-Cell Integration\"\nhttps://www.kaggle.com/competitions/open-problems-multimodal \n\n**Goals:** What are the reasonable blending stategies ? You need to balance: 1) better score 2) not to overfit 3) better simple than complicated  \n\n**Blending strategies:** The standard blending stategies include: \n\n    1) just average everything (most simple, and that is reasonable if models are of the same \"type\", and there are no params here - so no overfit (if solutions by themselves were not overfitted)). \n    2) order by scores and blend top-K scored \n    3) Sequential blend - with scalar weight optimization on each step: \n    Order the models by score from top to worst. \n    On each step - blend just result of previously blended models and one new model:\n    (  w * blend_of_previous + (1-w) * current_model )\n    \n    4) Find weights for linear blend by some optimization software (can put various natural constraints - for example to consider coefficients to be just [0,1,3..5], or positve weight or whatever) \n    4.1) Use scipy.optimizer\n    4.2) Use Optuna \n    \n    5) Find weights by some linear regression model. (We consider - individual models for each target - for multitarget problem !) \n    6) Simple Stacking - use non-linear model - e.g. LightGBM - to combine predictions. (We consider - individual models for each target - for multitarget problem !) \n    \n    7) GROUPwise blending - first group the models by the groups - for example: NN-models, GBDT-models, Linear models \n    (or use other grouping e.g. by clustering correlations). \n    In each group use simple blend - or average all, or average some topK models.\n    Then blend these results - that can be done again by aveagae all, or averaging topK , or even optimizing weights.\n    \n    That approach was considered previously in:\n    https://www.kaggle.com/code/alexandervc/mmscel-blend/notebook\n    For grouping by clustering correlations see discussion:\n    https://www.kaggle.com/competitions/tabular-playground-series-nov-2022/discussion/367413\n    \n    \n### Conclusions Briefly:\n    \n    The simple average approach seems quite often may fail to produce good results.\n    Average top1-3 seems often better.\n    \n    Sequential blend - is unexpectedly good - is simple, but it is only a bit worse than total weights optimization, but typically more interpreatable since includes only few models (at least for weight grid with step 0.1) - also we see striking phenomena that is produces sharp scores uplifts only for few models (which can be quite far from top1 model). \n    \n    GROUPwise blend - seems also same good - it is simple, but produce quite good results and interpretable. \n    \n    Scipy.optimize - works well for moderate number of models like 15 -  3 minutes.\n    But for 65 models - it will take about 1h 40 minutes - and seems no way to put timelimit. \n    (Pay attention that in current notebook - the loss function calculation is relatively slow). \n    OPTUNA - is in general worse than scipy.optimize - since it is not intended to such problems,\n    but the adavantage is that one can put a timelimit - that is quite desirable in the case of large number of models.\n    \n    TARGETWISE weights adjusting (for MULTITARGET problem) gives boost in the score. But it might be overfit - that is not explored here - we just only look on local CV scores.\n    Unexpectdely sometimes Ridge is better than LightGBM for such task. \n    \n**Overall:** one may consider using simple and interprerable blending strategies - sequential blend (N4), and groupwise blend (N7) - the results expected to be not far from optimal. And may be fine-tuning by scipy.optimize in some cases. \n\n\n### Conclusions In details:\n\n    1) Just pure average blend of ALL models - typically is worse than top1 model - somewhat unexpected - may be because we considered mostly homegenous groups of models without much diversity and with some not so small difference between top1 and the worst models.\n    \n    2) Simple average blending very few top models -  top1-2 or top1-3 models - typically gives the best results  - among  simple averages blends - i.e. no need to take many models. \n    \n    3) Step by step sequential blending - with scalar weight optimization on each step - going from top1 to the worst : \n       (  w * blend_of_previous + (1-w) * current_model ) .\n       Gives quite good(!!!) results in several respects: \n       \n       3.1) It scores almost as well as full optimization of all weights:  \\sum w_i * model_i - by scipy.optimize \n       \n       3.2) It is more interpretable - making many weights to be zero - (at least when we used grid with step 0.1) \n       \n       3.3) Even more interpretable - we see cute jumps of score for some few models which are top8 and top16 - \n       and thus taking only these models (joinlty with top1/2 ) we can achieve almost the same good results) \n       Thus we can reduce number of models to 3-5 out of dozens. \n       \n       3.4) Since it is more simple than optimization of all weights - it less prone to overfit. \n       \n       3.5) It can process 67 models in 10 minutes, while scipy.optimize - works orders of magnitude longer on 67 models. \n       (Well, for say 10 models - scipy.optimize - is faster) \n    \n    4) Scipy.optimize is typically better than Optuna at least for number of models less than 20, \n    it is faster and better results - not surprising since Optuna should work for non-smooth functions,\n    and for smooth functions classical optimization methods may work better. \n    But some detail: one can put timelimit to Optuna, while that seems to be impossible for scipy.optimize,\n    and so for larget number of models it may works toooooooo much time.... (So we suggest - sequential blend - described above).\n    \n    5,6) For the MULTItarget task (the main example for current notebook) - we can optimize weights separately for each target.\n    Typically we see strong improvement of results on our CV (but it quite might be prone to overfit ).\n    We compared blending by Ridge and LightGBM (stacking). (LightGBM was taken with the default params).\n    Difference was not that much strong, and Ridge was sometimes even better. \n    \n    7) GROUPwise blending - first group the models by the groups - for example: NN-models, GBDT-models, Linear models \n    (or use other grouping e.g. by clustering correlations). \n    In each group use simple blend - or average all, or average some topK models.\n    Then blend these results - that can be done again by aveagae all, or averaging topK , or even optimizing weights.\n    \n    That approach was considered previously in:\n    https://www.kaggle.com/code/alexandervc/mmscel-blend/notebook\n    For grouping by clustering correlations see discussion:\n    https://www.kaggle.com/competitions/tabular-playground-series-nov-2022/discussion/367413\n\n\n**Setup. Solutions - already calculated - just load them.** We produced and stored several dozens of solutions by boostings, NN, Ridge, etc... Both - out-of-fold and submission\nparts are stored in several public Kaggle datasets. In the present notebook we just load them, blend by various methods, compare the results. \n\nImportant detail: the competition was MULTItarget i.e. one needs to predict 140 targets simulatenously based on the same features. It was regression task. Actually there were two tasks in the competition - here we consider only the \"CITE-seq\" part. (The other part - \"Multiome\" is not considered here.) \n\n\nPrevious version of the notebook: \nhttps://www.kaggle.com/code/alexandervc/mmscel-blend/notebook#Prepare-data-for-submission\n\n\n### Versions\n\n\n#### 38,39,40, 41 - ALL(!) - 67 (or 65) models\n\n    Ridge blend targetwisely score: 0.899115\n    LGB targetwisely stacking score: 0.898957 \n    (Very high results - most probably overfit) \n\n    scipy.optimize best score: 0.894814\n    Wall time: 1h 43min 24s (!!!)\n    Optuna Best score: 0.894147\n    (3 minutes limit) \n\n\n    Best - Blend weighted sequential 15  Corr Score = 0.894546 blend_weight 1.0\n    (10min 32s it is quite sparse - only few solutions used )\n\n    Best - Blend top2\t0.893095\t0.122666\t0.217363\t\n    \n    BEST SOLO:\t\t0.892991\t0.135627\t0.214018\n    Blend Almost All \t0.887342\n    Blend ALL \t\t0.875747\n    WORST SOLO:\t\t0.875742\t0.053751\t0.248517\n\n#### 37 MLP_sklearn_with_bagging_etc  - MLP2  - 10 strong models - CV score from 0.890939 to 0.892438\n\n    CV score from 0.890939 to 0.892438\n\n    Ridge blend targetwisely score: 0.895625\n    LGB targetwisely stacking score: 0.895364\n\n    scipy.optimize best score: 0.893157\n    Optuna Best score: \t\t  0.893052\n\n    Blend weighted sequential 9\t0.893152\t0.127608\t0.218122\t\n    Blend weighted sequential 2\t0.893132\t0.127608\t0.218122\n\n    Blend top2 (best)\t0.893095\t0.127608\t0.218122\t\n\n    BEST SOLO:     0.892991\t0.135627\t0.214018\n    Blend ALL:\t0.890783\t0.127608\t0.218122\t\n    WORST SOLO: \t0.890729\t0.123751\t0.218543\n\n#### 36 MLP_Keras mlp-ver6-seed  - 7 strong models - CV score from 0.891866 to 0.892531 \n\n    CV score from 0.891866 to 0.892531\n\n    Ridge blend targetwisely score: 0.894712\n    LGB stacking score: 0.895025\n\n    scipy.optimize score: 0.893291\n    Optuna Best score: 0.893290\n\n    Blend weighted sequential 5\t0.893282\t0.132704\t0.215689\n    Blend weighted sequential 2\t0.893083\t0.132704\t0.215689\t\n\n    Blend top2\t0.892984\t0.132704\t0.215689\n    Blend top1\t0.892531\t0.132704\t0.215689\n\n    BEST SOLO mlp-ver6-seed0\t0.892531\t0.132896\t0.214939\n    Blend ALL:               0.891937\t0.132704\t0.215689\t\n    WORST SOLO mlp-ver6-outputs\t0.891866\t0.129989\t0.216268\n\n#### 35 LGBweak - 7 models LGB on 160 features (so relatively weak models) \n    \n    CV scores - from 0.887988 to 0.888052\n\n    Ridge blend  score: 0.890178\n    LGB stacking score: 0.890073\n\n    scipy.optimize score: 0.889464\n    Optuna Best score:  0.889463\n\n    Blend weighted sequential 6\t0.889459\t0.102188\t0.223915\n    Blend weighted sequential 4\t0.889236\t0.102188\t0.223915\t\n    Blend weighted sequential 1\t0.888940\t0.102188\t0.223915\n\n    Blend top1\t0.888940\t0.102188\t0.223915\t\n\n    Blend ALL\t\t\t0.888086\t0.102188\t0.223915\t\n    BEST SOLO Seed_62\t0.888052\t0.102284\t0.223896\n    WORST SOLO Seed_60\t0.887988\t0.102431\t0.224025\n\n\n#### 34  XGB models (12 or 13) \n    \n     Models are stronger in than previous LGB models, but have less diversity. \n     Scores from: from 0.891013 to 0.891780\n     Blend results are basically WORSE than for LGB (despite models are stronger in solo) \n     \n    Ridge targetwisely \t\t\t0.893154\n    LGB targetwisely \t\t\t0.892955\n\n    Optuna \t\t\t\t\t\t0.891973\n    Scipy.optimize              0.891978\n\n    Blend weighted sequential 3\t0.891932\t0.122666\t0.217363\n    Blend weighted sequential 2\t0.891906\t0.122666\t0.217363\n\n    Blend top2\t     \t\t0.891780\t0.122666\t0.217363\n    Top solo XGBseed20\t\t\t0.891411\t0.122191\t0.217178\t\n    Blend All           \t\t0.891024\n    WORST SOLO XGBseed446\t\t0.891013\t0.121317\t0.217974\t\n     \n\n\n#### 1-33 we use 15 LGB models\n\n    CV scores from 0.887956 to 0.891065\n        \n    from the 4 datasets:\n    https://www.kaggle.com/datasets/alexandervc/data-for-multimodal-singlecell-integration\n    https://www.kaggle.com/datasets/alexandervc/data3-multimodal-singlecell-integration\n    https://www.kaggle.com/datasets/alexandrgusev/lgb-nfeat550 \n    https://www.kaggle.com/datasets/alexandrgusev/train-size-08-features-0700-seed-639\n    \n    4) results: \n    WithOUT / WITH ROWwise scale to mean0 and std1  (results without - before version 32, results WITH - version 32)\n    0.894905 / 0.894776  - for each target blend (stack) - separately - by LGB - default params - takes 7-8 mins\n    0.894357 / 0.894282 - for each target blend - separately - by Ridge  \n    BIG UPLIFT - but quite may be overfitting  \n    0.892622 / 0.892618  - Scipy.optimizer on all 15 models - 3 minutes\n    0.892519 / 0.892514 - sequentially adding models one by one to the blend of the previous models \n    0.892411 / 0.892408 - sequentially top9 models (including tough uplift) \n    0.892288 - top1-2 & top9 (which was found to create tough uplift)\n    MAIN UPLIFT is HERE - there is one model (top9) - which blended - give significant uplift  \n    0.891289 / 0.891406- best of average sequential blend of topN - just top1&top2\n    0.891065 / same - top solo model\n    0.887982 / 0.887982 - simple average blend of all 15 models - it is almost as worst as the worst solo model \n    0.887956 / same - worst solo model \n\n\n#### 32 several LightGBM models. New option - rowwise - normalization of OOF predictions \n\n    Results - do not change much \n\n    Trick specific to that correlation metric used in that competition \n    flag_rescale_predicts_to_mean0_std1 = True # Rescaling predictions ROWwisely(!!!) to mean 0 and std1 - since correlation metric does not depend on such rescaling\n    # We were predicting targets which we transformed like that, so we may expect predictions are not far from satisfying that condition,\n    # but not exactly \n\n\n#### 1-31 several LightGBM models\n    \n**Details:**\n    \n    Models with slighly different training sets, features, random sets. \n    The solutions were created by the versions 99,100,101,102,  104,105 of the notebook: https://www.kaggle.com/code/alexandervc/mmscel-cv-modeling-advanced?scriptVersionId=111042323 And some copies of that notebook. Each running about 10 hours. \n    In the present notebook - we use these already precalculated solutions - results were stored in Kaggle datasets we can just load them.\n    In the present notebook we mostly work with the stored OOF predictions - so we can blend and caculate score comparing to given Y_true, without making a submission. Submission scores may be analyzed later.\n    \n**Conclusions/Observations:** \n\n    1) simple average blending just the first and the second scored gives better results comparing to including other models in simple blend by average\n    \n    2) Top9 solution unexpectedly improves blend quite bigger than the others. Reason is unclear. Top1+Top9 is better than Top1+Top2. And Sequential blend uplifts strongly for solution top9, not for the others.  \n    \n    3) scipy.optimize.minimize for 15 models works just 3 minutes - quite fast. It gives the best score (but we should be careful with possible overfit). \n    \n    4) results: \n    WithOUT / WITH ROWwise scale to mean0 and std1  (results without - before version 32, results WITH - version 32)\n    0.894905 / 0.894776  - for each target blend (stack) - separately - by LGB - default params - takes 7-8 mins\n    0.894357 / 0.894282 - for each target blend - separately - by Ridge  \n    BIG UPLIFT - but quite may be overfitting  \n    0.892622 / 0.892618  - Scipy.optimizer on all 15 models - 3 minutes\n    0.892519 / 0.892514 - sequentially adding models one by one to the blend of the previous models \n    0.892411 / 0.892408 - sequentially top9 models (including tough uplift) \n    0.892288 - top1-2 & top9 (which was found to create tough uplift)\n    MAIN UPLIFT is HERE - there is one model (top9) - which blended - give significant uplift  \n    0.891289 / 0.891406- best of average sequential blend of topN - just top1&top2\n    0.891065 / same - top solo model\n    0.887982 / 0.887982 - simple average blend of all 15 models - it is almost as worst as the worst solo model \n    0.887956 / same - worst solo model \n    \n    5) Optuna. As expected Optuna works worse for such kind of tasks comparing to standard optmizers from scipy.\n    The reason is clear - typically scores depend quite smoothly from the weights - so classical optimizers\n    based on gradient descent - work well, while optimizers like Optuna which are designed for more complicated tasks,\n    \"overthink\" and thus may achieve good results only very very slowly. \n\n\n### References \n\nThe competition on Kaggle - devoted how to create better blend strategies:\nhttps://www.kaggle.com/competitions/tabular-playground-series-nov-2022/overview\n(At the discussion forum one may find lots interesting ideas, e.g. see that one:\nhttps://www.kaggle.com/competitions/tabular-playground-series-nov-2022/discussion/367413 ) \n\n\n\nhttps://alexanderdyakonov.wordpress.com/2017/03/10/c%d1%82%d0%b5%d0%ba%d0%b8%d0%bd%d0%b3-stacking-%d0%b8-%d0%b1%d0%bb%d0%b5%d0%bd%d0%b4%d0%b8%d0%bd%d0%b3-blending/\nhttps://alexanderdyakonov.wordpress.com/2019/04/19/%d0%b0%d0%bd%d1%81%d0%b0%d0%bc%d0%b1%d0%bb%d0%b8-%d0%b2-%d0%bc%d0%b0%d1%88%d0%b8%d0%bd%d0%bd%d0%be%d0%bc-%d0%be%d0%b1%d1%83%d1%87%d0%b5%d0%bd%d0%b8%d0%b8/\n(From the former Kaggle Top1 and leading professor in Machine Learning - Alexander Dyakonov)\n\nGoogle for \"kaggle-ensembling-guide\":\n\nhttps://www.kaggle.com/code/mhill17/learning-the-kaggle-ensembling-guide\n\nhttps://www.kaggle.com/code/amrmahmoud123/1-guide-to-ensembling-methods/notebook\nhttps://www.kaggle.com/getting-started/117090\n\nhttps://www.kaggle.com/code/vipulgandhi/a-comprehensive-guide-to-ensemble-learning\n\nOther resources: \n\nhttps://www.projectpro.io/article/a-comprehensive-guide-to-ensemble-learning-methods/432\n\nhttps://towardsdatascience.com/ensemble-learning-stacking-blending-voting-b37737c4f483\n\nVery oftern cited (but site is always down):\nhttp://mlwave.com/kaggle-ensembling-guide/  -- there are often troubles accessing that webpage, may be consider github: https://github.com/MLWave/Kaggle-Ensemble-Guide \n\nhttps://opendatascience.com/a-survey-of-popular-ensemble-methods-part-1/\n\nhttps://machinelearningmastery.com/blending-ensemble-machine-learning-with-python/\n\nhttps://medium.com/@stevenyu530_73989/stacking-and-blending-intuitive-explanation-of-advanced-ensemble-methods-46b295da413c\n\nhttps://neptune.ai/blog/ensemble-learning-guide\n\nhttps://www.kaggle.com/code/amrmahmoud123/1-guide-to-ensembling-methods/notebook\n","metadata":{}},{"cell_type":"markdown","source":"# Key params","metadata":{}},{"cell_type":"code","source":"flag_run_stacking_by_lgb = True # It takes 7 minutes with default LGB params and 40 minutes with other  \n\nflag_rescale_predicts_to_mean0_std1 = True  # Specific feature of the MmSCel competition\n# Rescaling predictions ROWwisely(!!!) to mean 0 and std1 - since correlation metric does not depend on such rescaling\n# We were predicting targets which we transformed like that, so we may expect predictions are not far from satisfying that condition,\n# but not exactly \n\nlist_allowed_model_types =  ['_'] #  ['LGB' , 'XGB', 'LGB_NFeat160', 'MLP_ver6' , 'MLP2'  ] # \n# ['_'] - will include ALL predictions  ","metadata":{"execution":{"iopub.status.busy":"2022-11-27T18:51:36.668689Z","iopub.execute_input":"2022-11-27T18:51:36.669912Z","iopub.status.idle":"2022-11-27T18:51:36.675356Z","shell.execute_reply.started":"2022-11-27T18:51:36.669857Z","shell.execute_reply":"2022-11-27T18:51:36.674303Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Models Catalog","metadata":{}},{"cell_type":"markdown","source":"\n#### All 65 models\n\n    Include those listed below as well as number of others. \n    Typally not so strong - Ridge, one model by Catboost, 1DCNN, and very weak KNN.  \n\n#### MLP_sklearn_with_bagging_etc  - MLP2  - 10 strong models - CV score from 0.890939 to 0.892438\n\n\n    /kaggle/input/blend-results/MLP n_blends = 60 train_size = 0.7/MLP2_NFeat651Y_pred_oof_private_like.csv\n    /kaggle/input/blend-results/MLP n_blends = 10 train_size = 0.35/MLP2_NFeat651Y_pred_oof_private_like.csv\n    /kaggle/input/blend-results/MLP n_blends = 1 train_size = 1/MLP2_NFeat651Y_pred_oof_private_like.csv\n    /kaggle/input/blend-results/MLP n_blends = 10 train_size = 0.45/MLP2_NFeat651Y_pred_oof_private_like.csv\n    /kaggle/input/blend-results/MLP n_blends = 10 train_size = 0.75/MLP2_NFeat651Y_pred_oof_private_like.csv\n    /kaggle/input/blend-results/MLP n_blends = 10 train_size = 0.9/MLP2_NFeat651Y_pred_oof_private_like.csv\n    /kaggle/input/blend-results/MLP n_blends = 65 train_size = 0.9/MLP2_NFeat651Y_pred_oof_private_like.csv\n    /kaggle/input/blend-results/MultiMLP blends 72 train_size 0.9 part of x_train/MLP2_NFeat699Y_pred_oof_private_like.csv\n    /kaggle/input/blend-results/MLP n_blends = 60 train_size = 0.9/MLP2_NFeat651Y_pred_oof_private_like.csv\n    /kaggle/input/data-multimodal-singlecell-integration/MLP2_NFeat651/MLP2_NFeat651Y_pred_oof_private_like.csv\n\n\n#### MLP_Keras mlp-ver6-seed  - 7 strong models - CV score from 0.891866 to 0.892531\n\n    /kaggle/input/mlp-ver6-seed7/MLP_ver6_shevY_pred_oof_private_like.csv\n    /kaggle/input/mlp-ver6-seed4/MLP_ver6_shevY_pred_oof_private_like.csv\n    /kaggle/input/mlp-ver6-outputs/MLP_ver6_shevY_pred_oof_private_like.csv\n    /kaggle/input/mlp-ver6-seed0/MLP_ver6_shevY_pred_oof_private_like.csv\n    /kaggle/input/mlp-ver6-seed6/MLP_ver6_shevY_pred_oof_private_like.csv\n    /kaggle/input/mlp-ver6-seed2/MLP_ver6_shevY_pred_oof_private_like.csv\n    /kaggle/input/mlp-ver6-seed9/MLP_ver6_shevY_pred_oof_private_like.csv\n\n####  LGBweak  - 7 models on 160 features (relatively weak models). Cv scores: from 0.887988 to 0.888052\n\n    /kaggle/input/mmscel-cvmodeling-advanced-seeds-60-61-62/Seed_60/LGB_NFeat160Y_pred_oof_private_like.csv\n    /kaggle/input/mmscel-cvmodeling-advanced-seeds-60-61-62/Seed_61/LGB_NFeat160Y_pred_oof_private_like.csv\n    /kaggle/input/mmscel-cvmodeling-advanced-seeds-60-61-62/Seed_62/LGB_NFeat160Y_pred_oof_private_like.csv\n    /kaggle/input/data-multimodal-singlecell-integration/LGBweak2/LGB_NFeat160Y_pred_oof_private_like.csv\n    /kaggle/input/three-lgbms/rs71/LGB_NFeat160Y_pred_oof_private_like.csv\n    /kaggle/input/three-lgbms/rs72/LGB_NFeat160Y_pred_oof_private_like.csv\n    /kaggle/input/three-lgbms/rs70/LGB_NFeat160Y_pred_oof_private_like.csv\n\n####  XGB models (13) cv scores: from 0.891013 to 0.891411\n\n    '/kaggle/input/xgb-predictions-for-multimodal-competition/XGB_NFeat2719Y_pred_oof_private_like.csv',\n    '/kaggle/input/091seed112/XGB_NFeat2719Y_pred_oof_private_like.csv',\n    '/kaggle/input/xgb-1211/XGB_NFeat2719Y_pred_oof_private_like.csv',\n    '/kaggle/input/data-multimodal-singlecell-integration/XGBseed446/XGB_NFeat2719Y_pred_oof_private_like.csv',\n    '/kaggle/input/data-multimodal-singlecell-integration/XGBseed112/XGB_NFeat2719Y_pred_oof_private_like.csv',\n    '/kaggle/input/data-multimodal-singlecell-integration/XGBseed444/XGB_NFeat2719Y_pred_oof_private_like.csv',\n    '/kaggle/input/data-multimodal-singlecell-integration/XGBseed20/XGB_NFeat2291Y_pred_oof_private_like.csv',\n    '/kaggle/input/data-multimodal-singlecell-integration/XGBseed21/XGB_NFeat2291Y_pred_oof_private_like.csv',\n    '/kaggle/input/088-200xgb-nfeat2719y-pred-submission-kaggle-way/XGB_NFeat2719Y_pred_oof_private_like.csv',\n    '/kaggle/input/087-300-xgb-nfeat2719y-pred-submission-kaggle-way/XGB_NFeat2719Y_pred_oof_private_like (1).csv',\n    '/kaggle/input/data4-for-mmscel/XGB_NFeat2719Y_pred_oof_private_like.csv',\n    '/kaggle/input/data3-for-mmscel/XGB_NFeat2719Y_pred_oof_private_like.csv',\n    '/kaggle/input/data5-for-mmscel/XGB_NFeat2719Y_pred_oof_private_like.csv',\n\n#### LGB models (15) - on 500-700 features - thus reasonably strong comparing to LGBWeak on 150 features \n######    CV score from 0.887956 to 0.891065\n\n    /kaggle/input/train-size-08-features-0700-seed-639/LGB_NFeat700Y_pred_oof_private_like.csv\n    /kaggle/input/lgb-nfeat550/550_0.7_310/LGB_NFeat550Y_pred_oof_private_like.csv\n    /kaggle/input/lgb-nfeat550/550_0.4_105/LGB_NFeat550Y_pred_oof_private_like.csv\n    /kaggle/input/lgb-nfeat550/550_1_238/LGB_NFeat550Y_pred_oof_private_like.csv\n    /kaggle/input/lgb-nfeat550/550_0.8_442/LGB_NFeat550Y_pred_oof_private_like.csv\n    /kaggle/input/data3-multimodal-singlecell-integration/results (31)/LGB_NFeat499Y_pred_oof_private_like.csv\n    /kaggle/input/data3-multimodal-singlecell-integration/results (30)/LGB_NFeat499Y_pred_oof_private_like.csv\n    /kaggle/input/data3-multimodal-singlecell-integration/results (35)/LGB_NFeat499Y_pred_oof_private_like.csv\n    /kaggle/input/data3-multimodal-singlecell-integration/results (26)/LGB_NFeat600Y_pred_oof_private_like.csv\n    /kaggle/input/data3-multimodal-singlecell-integration/results (33)/LGB_NFeat499Y_pred_oof_private_like.csv\n    /kaggle/input/data3-multimodal-singlecell-integration/results (25)/LGB_NFeat600Y_pred_oof_private_like.csv\n    /kaggle/input/data3-multimodal-singlecell-integration/results (28)/LGB_NFeat570Y_pred_oof_private_like.csv\n    /kaggle/input/data3-multimodal-singlecell-integration/results (32)/LGB_NFeat499Y_pred_oof_private_like.csv\n    /kaggle/input/data3-multimodal-singlecell-integration/results (34)/LGB_NFeat499Y_pred_oof_private_like.csv\n    /kaggle/input/data3-multimodal-singlecell-integration/results (27)/LGB_NFeat580Y_pred_oof_private_like.csv\n    Total found: 15 OOF predictions. Shown only the first:  20\n        ","metadata":{}},{"cell_type":"markdown","source":"# Install/import","metadata":{}},{"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\ncount_loc = 0\nshow_how_many = 20\nprint('The first ', show_how_many , 'filenames with OOF predictions  will be shown:')\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        if 'pred_oof_private_like' in filename:\n            if np.sum( [t in filename for t in list_allowed_model_types] ) > 0: \n                count_loc += 1\n                if count_loc <= show_how_many :\n                    print(os.path.join(dirname, filename))\nprint('Total found:',count_loc, 'OOF predictions. Shown only the first: ', show_how_many  )\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","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-11-27T18:51:36.680665Z","iopub.execute_input":"2022-11-27T18:51:36.681727Z","iopub.status.idle":"2022-11-27T18:51:36.861827Z","shell.execute_reply.started":"2022-11-27T18:51:36.681665Z","shell.execute_reply":"2022-11-27T18:51:36.860575Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import time\nt0start = time.time()\nfrom datetime import datetime\ncurrent_date_and_time = datetime.now()\nprint(current_date_and_time)\n\nimport pandas as pd\nimport numpy as np\nimport os\nimport sys\n\nimport matplotlib.pyplot as plt\n#plt.style.use('dark_background')\nimport seaborn as sns\n\nif 0:\n    #If you see a urllib warning running this cell, go to \"Settings\" on the right hand side, \n    #and turn on internet. Note, you need to be phone verified.\n    !pip install --quiet tables\n\n\n    import h5py\n    !pip install hdf5plugin~=2.0 # https://forum.hdfgroup.org/t/cant-open-directory-usr-local-hdf5-lib-plugin/9738/4\n    import hdf5plugin\n\n    # !pip install scanpy\n    # import scanpy as sc\n    # import anndata\n\n    DATA_DIR = \"/kaggle/input/open-problems-multimodal/\"\n    FP_CELL_METADATA = os.path.join(DATA_DIR,\"metadata.csv\")\n\n    FP_CITE_TRAIN_INPUTS = os.path.join(DATA_DIR,\"train_cite_inputs.h5\")\n    FP_CITE_TRAIN_TARGETS = os.path.join(DATA_DIR,\"train_cite_targets.h5\")\n    FP_CITE_TEST_INPUTS = os.path.join(DATA_DIR,\"test_cite_inputs.h5\")\n\n    FP_MULTIOME_TRAIN_INPUTS = os.path.join(DATA_DIR,\"train_multi_inputs.h5\")\n    FP_MULTIOME_TRAIN_TARGETS = os.path.join(DATA_DIR,\"train_multi_targets.h5\")\n    FP_MULTIOME_TEST_INPUTS = os.pa","metadata":{"execution":{"iopub.status.busy":"2022-11-27T18:51:36.863980Z","iopub.execute_input":"2022-11-27T18:51:36.864312Z","iopub.status.idle":"2022-11-27T18:51:36.880886Z","shell.execute_reply.started":"2022-11-27T18:51:36.864282Z","shell.execute_reply":"2022-11-27T18:51:36.879579Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load Y_true: CITEseq_targets_rescaled","metadata":{}},{"cell_type":"code","source":"%%time\ndf_Y = pd.read_csv('/kaggle/input/data-for-multimodal-singlecell-integration/CITEseq_targets_rescaled.csv',index_col = 0)\nY_true = df_Y.values\nlist_all_targets = list(df_Y.columns)\nn_targets = len(list_all_targets)\nn_samples_train = Y_true.shape[0] # 70988\n\nprint(list_all_targets[:10] )\nprint(Y_true.shape)\ndf_Y","metadata":{"execution":{"iopub.status.busy":"2022-11-27T18:51:36.882236Z","iopub.execute_input":"2022-11-27T18:51:36.883167Z","iopub.status.idle":"2022-11-27T18:51:40.547688Z","shell.execute_reply.started":"2022-11-27T18:51:36.883132Z","shell.execute_reply":"2022-11-27T18:51:40.546402Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load Y_pred from various models ","metadata":{}},{"cell_type":"code","source":"%%time\ndict_df = {}\nlist_dirs_files = []\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        if '_pred_oof_private_like.csv' in filename:\n            if np.sum( [t in filename for t in list_allowed_model_types] ) > 0: \n                key_loc = dirname.split('/')[-1]\n                df = pd.read_csv(os.path.join(dirname, filename) , index_col = 0 )\n                \n                if df.shape[0] != n_samples_train: continue \n                if df.shape[1] != n_targets: continue \n                    \n                if flag_rescale_predicts_to_mean0_std1:\n                    t = df.values\n                    t -= t.mean(axis=1).reshape(-1, 1)\n                    t /= t.std(axis=1).reshape(-1, 1)\n                    df = pd.DataFrame(t, index = df.index, columns  = df.columns )\n                dict_df[key_loc] = df.copy()\n                print(key_loc, 'Shape:', dict_df[key_loc].shape, filename, dirname,)\n                print('      ',dict_df[key_loc].columns[:5], dict_df[key_loc].iloc[0,:5].values ) \n\nprint('Total:', len(dict_df) )\n            \n# for dn in list_dirs:\n#     for fn in os.listdir(dn):\n#         if '_pred_oof_private_like.csv' in fn:\n","metadata":{"execution":{"iopub.status.busy":"2022-11-27T18:51:40.551160Z","iopub.execute_input":"2022-11-27T18:51:40.551637Z","iopub.status.idle":"2022-11-27T18:54:47.069081Z","shell.execute_reply.started":"2022-11-27T18:51:40.551592Z","shell.execute_reply":"2022-11-27T18:54:47.067988Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('Total:', len(dict_df) )\n\n","metadata":{"execution":{"iopub.status.busy":"2022-11-27T18:54:47.070297Z","iopub.execute_input":"2022-11-27T18:54:47.070608Z","iopub.status.idle":"2022-11-27T18:54:47.076577Z","shell.execute_reply.started":"2022-11-27T18:54:47.070578Z","shell.execute_reply":"2022-11-27T18:54:47.075533Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Correlation scoring function\n\nCorrelation score - but it is tricky(!) - it is computed along the targets dimension, and average over the samples \n\n\nSee e.g. https://mathoverflow.net/questions/434254/how-to-maximize-certain-function-of-hundreds-variables-related-to-correlations-b\n","metadata":{}},{"cell_type":"code","source":"import gc\nrescale_Y_to_mean0_std1 = True\n\ndef get_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    Returns the average of each sample's Pearson correlation coefficient\n    \"\"\"\n\n    ############################\n    # Faster way:\n    ############################\n    \n    # Input should be matrices - does not make sense for vectors \n    if  len(y_pred.shape)< 2: return -10 # Some result to inform for incorrect input\n    if  y_pred.shape[1] < 2: return -10 # Some result to inform for incorrect input\n\n    y2 = y_pred.copy()\n    y2 -= y2.mean(axis=1).reshape(-1, 1);    y2 /= y2.std(axis=1).reshape(-1, 1)    \n    if rescale_Y_to_mean0_std1:\n        y1 = y_true # Already rescaled \n    else:\n        y1 = y_true.copy(); \n        y1 -= y1.mean(axis=1).reshape(-1, 1);    y1 /= y1.std(axis=1).reshape(-1, 1) \n        \n    c = (y1*y2).mean().mean()# Correlation for rescaled matrices is just matrix product and average \n    \n    c = (y1*y2).mean().mean()# Correlation for rescaled matrices is just matrix product and average \n    \n    # Memory control:\n    if not rescale_Y_to_mean0_std1:\n        del y1\n    del y2\n    gc.collect()\n    \n    return c\n\n    ############################\n    # Slower way:\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)","metadata":{"execution":{"iopub.status.busy":"2022-11-27T18:54:47.078249Z","iopub.execute_input":"2022-11-27T18:54:47.078712Z","iopub.status.idle":"2022-11-27T18:54:47.092004Z","shell.execute_reply.started":"2022-11-27T18:54:47.078664Z","shell.execute_reply":"2022-11-27T18:54:47.090539Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Calculate correlation scores for all input solutions","metadata":{}},{"cell_type":"code","source":"%%time\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.metrics import r2_score\n\ndf_stat = pd.DataFrame()\nIX = 0\nlist_names = list(dict_df.keys())\nfor i,key_loc in enumerate(list( dict_df.keys())  ):\n    key_loc = list_names[i]\n    df = dict_df[key_loc]\n    print(key_loc, df.shape)\n    #display(df.head(2))\n    s = get_score(Y_true,df.values)\n    df_stat.loc[IX,'Model'] = key_loc\n    df_stat.loc[IX,'Corr'] = s\n    df_stat.loc[IX,'r2'] = r2_score(Y_true,df.values)\n    df_stat.loc[IX,'MSE'] = mean_squared_error(Y_true,df.values)\n    \n#     m = mask_holdout; postfix = ' HO'\n#     df_stat.loc[IX,'Corr'+ postfix] = get_score(Y_true[m],df[m].values)\n#     df_stat.loc[IX,'r2'+ postfix] = r2_score(Y_true[m],df[m].values)\n#     df_stat.loc[IX,'MSE'+ postfix] = mean_squared_error(Y_true[m],df[m].values)\n#     m = mask_main;  postfix = ' Main'\n#     df_stat.loc[IX,'Corr'+postfix] = get_score(Y_true[m],df[m].values)\n#     df_stat.loc[IX,'r2'+postfix] = r2_score(Y_true[m],df[m].values)\n#     df_stat.loc[IX,'MSE'+postfix] = mean_squared_error(Y_true[m],df[m].values)\n    \n    IX += 1\n    \ndisplay(df_stat)","metadata":{"execution":{"iopub.status.busy":"2022-11-27T18:54:47.093387Z","iopub.execute_input":"2022-11-27T18:54:47.093717Z","iopub.status.idle":"2022-11-27T18:55:46.075265Z","shell.execute_reply.started":"2022-11-27T18:54:47.093686Z","shell.execute_reply":"2022-11-27T18:55:46.074020Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Statistics on solutions","metadata":{}},{"cell_type":"code","source":"import matplotlib.pyplot as plt\ndf_stat_sorted = df_stat.sort_values('Corr',ascending = False)# ['Corr']\ndisplay( df_stat_sorted )\ndisplay(df_stat_sorted.describe())\nfig = plt.figure(figsize = (20,4))\nplt.plot(df_stat_sorted.set_index('Model')['Corr'],'*-')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-11-27T18:55:46.077174Z","iopub.execute_input":"2022-11-27T18:55:46.077934Z","iopub.status.idle":"2022-11-27T18:55:47.120779Z","shell.execute_reply.started":"2022-11-27T18:55:46.077884Z","shell.execute_reply":"2022-11-27T18:55:47.119518Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Blend from top scored solution to bottom scored solution ","metadata":{}},{"cell_type":"code","source":"%%time\ndf_stat_with_blend = df_stat.copy()\nIX = max(df_stat_with_blend.index)+1\nfor i,model_id in  enumerate(df_stat_sorted['Model']):\n    if i == 0:\n        Y_blend_current = dict_df[model_id]\n    else:\n        Y_blend_current = (Y_blend_current/i + dict_df[model_id])/(i+1)\n    #print(Y_blend_current.shape, type(Y_blend_current))\n    s = get_score(Y_true,Y_blend_current.values)\n    key_loc = 'Blend top'+str(i)\n    print(key_loc, ' Corr Score = %.6f'%(s) )\n    df_stat_with_blend.loc[IX,'Model'] = key_loc\n    df_stat_with_blend.loc[IX,'Corr'] = s\n    df_stat_with_blend.loc[IX,'r2'] = r2_score(Y_true,df.values)\n    df_stat_with_blend.loc[IX,'MSE'] = mean_squared_error(Y_true,df.values)\n    IX += 1\n\ndf_stat_with_blend['Blend'] = df_stat_with_blend['Model'].apply( lambda x: 1 if 'Blend' in x else 0 )\nmask_blend = df_stat_with_blend['Blend'] == 1\n\nplt.figure(figsize = (20,4))\nplt.plot(df_stat_with_blend[mask_blend].sort_values('Corr',ascending = False ).set_index('Model')['Corr'],'*-')\nplt.title('Sorted Blend results')\nplt.show()\n\nplt.figure(figsize = (20,4))\nplt.plot(df_stat_with_blend[mask_blend].set_index('Model')['Corr'],'*-')\nplt.title('UN Sorted Blend results')\nplt.show()\n\ndisplay(df_stat_with_blend.sort_values('Corr',ascending = False ))    ","metadata":{"execution":{"iopub.status.busy":"2022-11-27T18:55:47.122357Z","iopub.execute_input":"2022-11-27T18:55:47.123457Z","iopub.status.idle":"2022-11-27T18:56:53.073702Z","shell.execute_reply.started":"2022-11-27T18:55:47.123415Z","shell.execute_reply":"2022-11-27T18:56:53.072465Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Blend top 2 with optimization of weight","metadata":{}},{"cell_type":"code","source":"model_id1 = df_stat_sorted['Model'].iat[0]\nmodel_id2 = df_stat_sorted['Model'].iat[1]\n\n\nY_blend_current = np.zeros_like( dict_df[model_id1])\nlist_scores_loc = []\nfor blend_weight in np.linspace(0,1,11):\n    Y_blend_current = blend_weight* dict_df[model_id1] + (1-blend_weight)* dict_df[model_id2]     \n    s = get_score(Y_true,Y_blend_current.values)\n    key_loc = 'Blend weight = '+str(blend_weight)\n    print(key_loc, ' Corr Score = %.6f'%(s) )\n    list_scores_loc.append(s)\n\nprint('Best score:', np.max(list_scores_loc) )    \n    \nplt.figure(figsize = (20,4))\nplt.plot(np.linspace(0,1,11),  list_scores_loc, '*-')\nplt.title('Dependence of blend score on weight',fontsize = 20 )\nplt.ylabel('Score',fontsize = 20)\nplt.xlabel('Weight:   w * (solution1) + (1-w) * solution2 ',fontsize = 15)\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2022-11-27T18:56:53.077592Z","iopub.execute_input":"2022-11-27T18:56:53.078058Z","iopub.status.idle":"2022-11-27T18:57:02.226570Z","shell.execute_reply.started":"2022-11-27T18:56:53.078025Z","shell.execute_reply":"2022-11-27T18:57:02.225328Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3D visualization with two weights changing ","metadata":{}},{"cell_type":"code","source":"%%time\n# from sklearn.model_selection import ParameterGrid\n# grid = {'w1': [0, 0.25,0.5, 0.75, 1], 'w2': [0, 0.25,0.5, 0.75, 1]}\n# list(ParameterGrid(grid))\n# import itertools\n# l = [0, 0.25,0.5, 0.75, 1]\n# grid = list(itertools.product(l, l))\n# grid\n\n\nif len(df_stat_sorted) > 2:\n    model_id0 = df_stat_sorted['Model'].iat[0]\n    model_id1 = df_stat_sorted['Model'].iat[1]\n    model_id2 = df_stat_sorted['Model'].iat[2]\n\n    list_1 = np.linspace(0, 1, 5)\n    list_2 = np.linspace(0, 1, 5)\n\n    Z = np.zeros( (len(list_1), len(list_2) )  )\n    for i,w1 in enumerate(list_1):\n        for j,w2 in enumerate(list_2):\n            w0_normalized = 1/(1+w1+w2)\n            w1_normalized = w1/(1+w1+w2)\n            w2_normalized = w2/(1+w1+w2)\n            \n            Y_blend = w0_normalized * dict_df[model_id1] + w1_normalized * dict_df[model_id1] + w2_normalized * dict_df[model_id2]\n            \n            s = get_score(Y_true, Y_blend.values)\n        \n            Z[i,j] = s\n            \n            \nX, Y = np.meshgrid(list_1, list_2)\nfig = plt.figure()\nax = plt.axes(projection='3d')\nax.plot_surface(X, Y, Z, rstride=1, cstride=1,\n                cmap='viridis', edgecolor='none')\n\n# ax.contour3D(X, Y, Z, 50, cmap='binary')\n\n\nax.set_xlabel('W1')\nax.set_ylabel('W2')\nax.view_init(60, 135)\n\n#fig            ","metadata":{"execution":{"iopub.status.busy":"2022-11-27T18:57:02.228196Z","iopub.execute_input":"2022-11-27T18:57:02.228543Z","iopub.status.idle":"2022-11-27T18:57:24.092360Z","shell.execute_reply.started":"2022-11-27T18:57:02.228512Z","shell.execute_reply":"2022-11-27T18:57:24.091251Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.model_selection import ParameterGrid\n\nif len(df_stat_sorted) > 2:\n    model_id0 = df_stat_sorted['Model'].iat[0]\n    model_id1 = df_stat_sorted['Model'].iat[1]\n    model_id2 = df_stat_sorted['Model'].iat[2]\n\n    grid = {'w1': [0, 0.25,0.5, 0.75, 1], 'w2': [0, 0.25,0.5, 0.75, 1]}\n    list(ParameterGrid(grid))\n\n    list_x = []; list_y = []; list_z = []\n    for params in list(ParameterGrid(grid)):\n        #print(params)\n        w1 = params['w1']; w2 = params['w2']\n        list_x.append(w1); list_y.append(w2);\n\n        w0_normalized = 1/(1+w1+w2)\n        w1_normalized = w1/(1+w1+w2)\n        w2_normalized = w2/(1+w1+w2)\n\n        Y_blend = w0_normalized * dict_df[model_id1] + w1_normalized * dict_df[model_id1] + w2_normalized * dict_df[model_id2]\n\n        s = get_score(Y_true, Y_blend.values)\n        list_z.append(s)      ","metadata":{"execution":{"iopub.status.busy":"2022-11-27T18:57:24.094153Z","iopub.execute_input":"2022-11-27T18:57:24.094933Z","iopub.status.idle":"2022-11-27T18:57:45.680015Z","shell.execute_reply.started":"2022-11-27T18:57:24.094885Z","shell.execute_reply":"2022-11-27T18:57:45.679016Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import plotly.graph_objects as go\n\nz = np.array(list_z).reshape( (len(list_1), len(list_2)) )\nz\n\nfig = go.Figure(data=[go.Surface(z=z, x=list_1, y=list_2)])\nfig.update_layout(title='Blend score', autosize=False,\n                  width=500, height=500,\n                  margin=dict(l=65, r=50, b=65, t=90))\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2022-11-27T18:57:45.681720Z","iopub.execute_input":"2022-11-27T18:57:45.682430Z","iopub.status.idle":"2022-11-27T18:57:45.700187Z","shell.execute_reply.started":"2022-11-27T18:57:45.682374Z","shell.execute_reply":"2022-11-27T18:57:45.699061Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nimport plotly.express as px\n\nv4color_loc = list_z\nfig = px.scatter_3d(x=list_x, y=list_y, z=list_z, color=v4color_loc, opacity=0.9)\n#fig = px.line_3d(x=list_x, y=list_y, z=list_z, color=v4color_loc)#, opacity=0.9)\nfig.update_traces(marker_size = 5)\n\nfig.show()\n","metadata":{"execution":{"iopub.status.busy":"2022-11-27T18:57:45.701596Z","iopub.execute_input":"2022-11-27T18:57:45.701945Z","iopub.status.idle":"2022-11-27T18:57:45.773452Z","shell.execute_reply.started":"2022-11-27T18:57:45.701914Z","shell.execute_reply":"2022-11-27T18:57:45.772311Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Sequential blending with search for optimal weight on each step","metadata":{}},{"cell_type":"code","source":"\ndef get_optimal_weight(Y_true, Y_pred1, Y_pred2):\n\n    s_best = -np.inf\n    \n    for blend_weight in np.linspace(0,1,11):\n        Y_blend_current = blend_weight * Y_pred1 + (1-blend_weight)* Y_pred2\n        if hasattr(Y_blend_current, 'values'):\n            s = get_score(Y_true,Y_blend_current.values)\n        else:\n            s = get_score(Y_true,Y_blend_current)\n        if s > s_best:\n            s_best = s\n            blend_weight_best = blend_weight\n    return blend_weight_best\n            \n","metadata":{"execution":{"iopub.status.busy":"2022-11-27T18:57:45.775093Z","iopub.execute_input":"2022-11-27T18:57:45.775549Z","iopub.status.idle":"2022-11-27T18:57:45.784955Z","shell.execute_reply.started":"2022-11-27T18:57:45.775504Z","shell.execute_reply":"2022-11-27T18:57:45.783590Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndf_stat_with_blend = df_stat.copy()\nIX = max(df_stat_with_blend.index)+1\nfor i,model_id in  enumerate(df_stat_sorted['Model']):\n    if i == 0:\n        Y_blend_current = dict_df[model_id]\n        blend_weight = 0\n    else:\n        blend_weight = get_optimal_weight(Y_true, Y_blend_current, dict_df[model_id])\n        Y_blend_current = blend_weight * Y_blend_current + (1 - blend_weight ) * dict_df[model_id]\n    #print(Y_blend_current.shape, type(Y_blend_current))\n    s = get_score(Y_true,Y_blend_current.values)\n    key_loc = 'Blend weighted sequential '+str(i)\n    print(key_loc, ' Corr Score = %.6f'%(s), 'blend_weight', blend_weight )\n    df_stat_with_blend.loc[IX,'Model'] = key_loc\n    df_stat_with_blend.loc[IX,'Corr'] = s\n    df_stat_with_blend.loc[IX,'r2'] = r2_score(Y_true,df.values)\n    df_stat_with_blend.loc[IX,'MSE'] = mean_squared_error(Y_true,df.values)\n    IX += 1\n\ndf_stat_with_blend['Blend'] = df_stat_with_blend['Model'].apply( lambda x: 1 if 'Blend' in x else 0 )\nmask_blend = df_stat_with_blend['Blend'] == 1\n\n\nplt.figure(figsize = (20,4))\nplt.plot(df_stat_with_blend[mask_blend].set_index('Model')['Corr'].values ,'*-', label = 'blend')\nplt.plot(df_stat_sorted.set_index('Model')['Corr'].values ,'*-', label = 'solo')\nplt.legend()\nplt.title('Blend and solo results', fontsize = 20)\nplt.show()\n\n\ndisplay(df_stat_with_blend.sort_values('Corr',ascending = False ))        ","metadata":{"execution":{"iopub.status.busy":"2022-11-27T18:57:45.786464Z","iopub.execute_input":"2022-11-27T18:57:45.786821Z","iopub.status.idle":"2022-11-27T19:08:18.331114Z","shell.execute_reply.started":"2022-11-27T18:57:45.786789Z","shell.execute_reply":"2022-11-27T19:08:18.329828Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Blend top 1 with selected ","metadata":{}},{"cell_type":"code","source":"if len(df_stat_sorted) > 8:\n    model_id1 = df_stat_sorted['Model'].iat[0]\n    model_id2 = df_stat_sorted['Model'].iat[8]\n\n\n    Y_blend_current = np.zeros_like( dict_df[model_id1])\n    list_scores_loc = []\n    for blend_weight in np.linspace(0,1,11):\n        Y_blend_current = blend_weight* dict_df[model_id1] + (1-blend_weight)* dict_df[model_id2]     \n        s = get_score(Y_true,Y_blend_current.values)\n        key_loc = 'Blend weight = '+str(blend_weight)\n        print(key_loc, ' Corr Score = %.6f'%(s) )\n        list_scores_loc.append(s)\n\n    print('Best score:', np.max(list_scores_loc) )    \n    plt.figure(figsize = (20,4))\n    plt.plot(np.linspace(0,1,11),  list_scores_loc, '*-')\n    plt.title('Dependence of blend score on weight',fontsize = 20 )\n    plt.ylabel('Score',fontsize = 20)\n    plt.xlabel('Weight:   w * (solution1) + (1-w) * solution2 ',fontsize = 15)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-11-27T19:08:18.332955Z","iopub.execute_input":"2022-11-27T19:08:18.333318Z","iopub.status.idle":"2022-11-27T19:08:27.462152Z","shell.execute_reply.started":"2022-11-27T19:08:18.333285Z","shell.execute_reply":"2022-11-27T19:08:27.461150Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Toy Example scipy optimize","metadata":{}},{"cell_type":"code","source":"# https://docs.scipy.org/doc/scipy/tutorial/optimize.html\nfrom scipy.optimize import minimize#, fsolve\nbnds = ((1, None), (0, None))\ndef f0(x):\n    return x[0]**2+x[1]**2\nres_scipy = minimize( f0,  [10,1] , bounds=bnds,)\nres_scipy","metadata":{"execution":{"iopub.status.busy":"2022-11-27T19:08:27.463551Z","iopub.execute_input":"2022-11-27T19:08:27.463924Z","iopub.status.idle":"2022-11-27T19:08:27.474837Z","shell.execute_reply.started":"2022-11-27T19:08:27.463892Z","shell.execute_reply":"2022-11-27T19:08:27.473704Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Use scipy optimize to find weights","metadata":{}},{"cell_type":"markdown","source":"## Example with three best solutions found before","metadata":{}},{"cell_type":"code","source":"%%time\n\n# https://docs.scipy.org/doc/scipy/tutorial/optimize.html\nfrom scipy.optimize import minimize#, fsolve\n\nlist_preds = []\nmodel_id = df_stat_sorted['Model'].iat[0]; list_preds.append( dict_df[model_id] )\nmodel_id = df_stat_sorted['Model'].iat[1]; list_preds.append( dict_df[model_id] )\nif len(df_stat_sorted) > 8:\n    model_id = df_stat_sorted['Model'].iat[8]; list_preds.append( dict_df[model_id] )\n\nprint('Number of models for blend:', len( list_preds ) )\n\n    \ndef func2optm(list_blend_weights):\n    Y_blend = np.zeros(Y_true.shape,dtype = float )   \n    for i,y_pred in enumerate(list_preds):\n        Y_blend += list_blend_weights[i]*list_preds[i]\n    score_loc = get_score(Y_true,Y_blend.values)\n    return -score_loc\n\nx0 = np.ones(len(list_preds)) * (1.0/len(list_preds)) \nbnds = []\nfor i,_ in enumerate(list_preds):\n    bnds.append((0,None))\n    \nres_scipy = minimize(fun = func2optm, x0=  x0, bounds = bnds  )    \nprint(res_scipy)\n\nprint()\nblend_weights = res_scipy.x / np.sum(res_scipy.x)\nprint('scipy.optimize Found blend weights:' , blend_weights )\nblend_score = - func2optm(blend_weights)\nprint('scipy.optimize Found blend score' , blend_score )\nprint()","metadata":{"execution":{"iopub.status.busy":"2022-11-27T19:08:27.476374Z","iopub.execute_input":"2022-11-27T19:08:27.476915Z","iopub.status.idle":"2022-11-27T19:08:50.668563Z","shell.execute_reply.started":"2022-11-27T19:08:27.476881Z","shell.execute_reply":"2022-11-27T19:08:50.667127Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Blend with scipy.optimize.minimize  for all models\n\n15 models - 3 mins 18 secs - quite fast","metadata":{}},{"cell_type":"code","source":"%%time\n\n# https://docs.scipy.org/doc/scipy/tutorial/optimize.html\nfrom scipy.optimize import minimize#, fsolve\n\nlist_preds = []\nfor model_id in df_stat_sorted['Model']:\n    list_preds.append( dict_df[model_id] )\n\nprint('Number of models for blend:', len( list_preds ) )\n    \ndef func2optm(list_blend_weights):\n    Y_blend = np.zeros(Y_true.shape,dtype = float )   \n    for i,y_pred in enumerate(list_preds):\n        Y_blend += list_blend_weights[i]*list_preds[i]\n    score_loc = get_score(Y_true,Y_blend.values)\n    return -score_loc\n\nx0 = np.ones(len(list_preds)) * (1.0/len(list_preds)) \nbnds = []\nfor i,_ in enumerate(list_preds):\n    bnds.append((0,None))\n    \nres_scipy = minimize(fun = func2optm, x0=  x0, bounds = bnds  )    \nprint(res_scipy)\n\nprint()\nblend_weights = res_scipy.x / np.sum(res_scipy.x)\nprint('scipy.optimize Found blend weights:' , blend_weights )\nblend_score = - func2optm(blend_weights)\nprint('scipy.optimize best score:' , '%.6f'%blend_score)\nprint('%.6f'%blend_score)\nprint(blend_score)","metadata":{"execution":{"iopub.status.busy":"2022-11-27T19:08:50.670449Z","iopub.execute_input":"2022-11-27T19:08:50.670814Z","iopub.status.idle":"2022-11-27T20:52:15.311678Z","shell.execute_reply.started":"2022-11-27T19:08:50.670782Z","shell.execute_reply":"2022-11-27T20:52:15.310747Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import gc\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-11-27T20:52:15.314856Z","iopub.execute_input":"2022-11-27T20:52:15.315486Z","iopub.status.idle":"2022-11-27T20:52:15.805044Z","shell.execute_reply.started":"2022-11-27T20:52:15.315450Z","shell.execute_reply":"2022-11-27T20:52:15.803821Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Optimization with Optuna\n\nAs expected Optuna works worse for such kind of tasks comparing to standard optmizers from scipy.\nThe reason is clear - typically scores depend quite smoothly from the weights - so classical optimizers\nbased on gradient descent - work well, while optimizers like Optuna which are designed for more complicated tasks,\n\"overthink\" and thus may achieve good results only very very slowly. ","metadata":{}},{"cell_type":"code","source":"import optuna","metadata":{"execution":{"iopub.status.busy":"2022-11-27T20:52:15.807295Z","iopub.execute_input":"2022-11-27T20:52:15.807785Z","iopub.status.idle":"2022-11-27T20:52:15.816421Z","shell.execute_reply.started":"2022-11-27T20:52:15.807717Z","shell.execute_reply":"2022-11-27T20:52:15.815358Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"list_preds = []\nfor model_id in df_stat_sorted['Model']:\n    list_preds.append( dict_df[model_id] )\n\ndef objectives(trial):\n    list_blend_weights = []\n    \n    for i in range(0,len(list_preds) ):\n        w = trial.suggest_float(\"W\"+str(i), 0, 1, step = 0.01) \n        list_blend_weights.append(w)\n        \n    Y_blend = np.zeros(Y_true.shape,dtype = float )   \n    for i,y_pred in enumerate(list_preds):\n        Y_blend += list_blend_weights[i]*list_preds[i]\n    score_loc = get_score(Y_true,Y_blend.values)\n    return -score_loc","metadata":{"execution":{"iopub.status.busy":"2022-11-27T20:52:15.818351Z","iopub.execute_input":"2022-11-27T20:52:15.818872Z","iopub.status.idle":"2022-11-27T20:52:15.831851Z","shell.execute_reply.started":"2022-11-27T20:52:15.818826Z","shell.execute_reply":"2022-11-27T20:52:15.830562Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"study = optuna.create_study(direction='minimize', sampler=optuna.samplers.TPESampler())\noptuna.logging.set_verbosity(optuna.logging.INFO)\n\n#study.enqueue_trial({'K2': 0.5524100000000001, 'K3': 0.30762, 'K4': 0.026800000000000157})\n\noptuna.logging.set_verbosity(optuna.logging.WARNING)\nimport warnings\nwarnings.filterwarnings(\"ignore\", category=FutureWarning)","metadata":{"execution":{"iopub.status.busy":"2022-11-27T20:52:15.834064Z","iopub.execute_input":"2022-11-27T20:52:15.834623Z","iopub.status.idle":"2022-11-27T20:52:15.846689Z","shell.execute_reply.started":"2022-11-27T20:52:15.834575Z","shell.execute_reply":"2022-11-27T20:52:15.845623Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Start main optimization","metadata":{}},{"cell_type":"code","source":"%%time\nn_trials = 1000\nstudy.optimize(objectives, n_trials=n_trials, timeout = 180)","metadata":{"execution":{"iopub.status.busy":"2022-11-27T20:52:15.848166Z","iopub.execute_input":"2022-11-27T20:52:15.848650Z","iopub.status.idle":"2022-11-27T20:55:20.114955Z","shell.execute_reply.started":"2022-11-27T20:52:15.848601Z","shell.execute_reply":"2022-11-27T20:55:20.113633Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"list_blend_weights = []\nfor key_loc in study.best_params:\n    list_blend_weights.append(study.best_params[key_loc])\nY_blend = np.zeros(Y_true.shape,dtype = float )   \nfor i,y_pred in enumerate(list_preds):\n    Y_blend += list_blend_weights[i]*list_preds[i]\nscore_loc = get_score(Y_true,Y_blend.values)\n\nprint(); print('Optuna Best params:')\nprint(list(np.round(list_blend_weights,6))  )\nprint(); print('Optuna Best score:', '%.6f'%score_loc )\nprint('%.6f'%score_loc)\nprint(score_loc)","metadata":{"execution":{"iopub.status.busy":"2022-11-27T20:55:20.116320Z","iopub.execute_input":"2022-11-27T20:55:20.116668Z","iopub.status.idle":"2022-11-27T20:55:25.174053Z","shell.execute_reply.started":"2022-11-27T20:55:20.116637Z","shell.execute_reply":"2022-11-27T20:55:25.173129Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n# Result for LGB - 15 models without Y_oof rescaling:\n\n# Worse than scipy optimize:\n\n# Optuna: Best score:  0.8925689750637171\n\n# Scipy:  Best score   0.8926221034415455\n\n# Best params:\n# [1.0, 0.22, 0.47, 0.8, 0.77, 0.44, 0.05, 0.1, 1.0, 0.9, 0.43, 0.08, 0.59, 0.33, 0.1]\n","metadata":{"execution":{"iopub.status.busy":"2022-11-27T20:55:25.175302Z","iopub.execute_input":"2022-11-27T20:55:25.176038Z","iopub.status.idle":"2022-11-27T20:55:25.180701Z","shell.execute_reply.started":"2022-11-27T20:55:25.176002Z","shell.execute_reply":"2022-11-27T20:55:25.179208Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Ridge blend for each target - individual weights for each target","metadata":{"execution":{"iopub.status.busy":"2022-11-25T08:54:35.426989Z","iopub.execute_input":"2022-11-25T08:54:35.427377Z","iopub.status.idle":"2022-11-25T08:54:35.433745Z","shell.execute_reply.started":"2022-11-25T08:54:35.427347Z","shell.execute_reply":"2022-11-25T08:54:35.432801Z"}}},{"cell_type":"code","source":"%%time\nfrom sklearn.linear_model import RidgeCV\nfrom sklearn.model_selection import cross_val_predict\nfrom sklearn.model_selection import KFold\nfrom sklearn.metrics import r2_score\nfrom sklearn.model_selection import cross_val_score\nfrom sklearn.linear_model import Ridge\n\n\nkf = KFold(n_splits=3, random_state = 0, shuffle = True)\n\n\nn_targets = len(list_all_targets)\nn_models = len( df_stat_sorted['Model']  )\nprint('n_samples_train', n_samples_train ) # \n\nlist_models_ids = list( df_stat_sorted['Model'] ) \n\ndf_data_from_all_models_for_one_target = pd.DataFrame(np.zeros( (n_samples_train,  n_models)  ), columns =  list_models_ids   )\n\ndf_stat_targetwise = pd.DataFrame(); IX = 0 \nY_blend = np.zeros(Y_true.shape,dtype = float )   \n\nfor target_id in range(n_targets):\n    for i,model_id in  enumerate(df_stat_sorted['Model']):\n        data_from_current_model = dict_df[model_id].iloc[:,target_id]\n        df_data_from_all_models_for_one_target.iloc[:,i] = data_from_current_model.values\n    Y_true_for_one_target  = Y_true[:, target_id] \n    model_loc = RidgeCV(alphas=[1e-1, 1, 1e1, 1e2, 1e3, 1e4]).fit(df_data_from_all_models_for_one_target, Y_true_for_one_target)\n    \n    model_loc2 = Ridge( alpha = model_loc.alpha_ )\n\n    Y_blend[:,target_id] = cross_val_predict(model_loc2  ,df_data_from_all_models_for_one_target, Y_true_for_one_target , cv = kf  )\n\n    s_oof =  r2_score( Y_true_for_one_target, Y_blend[:,target_id] )\n    list_cv_scores = cross_val_score(model_loc2, df_data_from_all_models_for_one_target, Y_true_for_one_target, cv=kf)\n    \n    s = model_loc.score(df_data_from_all_models_for_one_target, Y_true_for_one_target)\n    #print(list_all_targets[target_id], 'target_id', target_id, model_loc.alpha_, 'r2 score %0.3f'%s )\n    df_stat_targetwise.loc[IX,'Target'] = list_all_targets[target_id]\n    df_stat_targetwise.loc[IX,'OOF r2 score'] = s_oof\n    df_stat_targetwise.loc[IX,'CV (mean) r2 score'] = np.mean(list_cv_scores)\n    df_stat_targetwise.loc[IX,'Train r2 score'] = s\n    df_stat_targetwise.loc[IX,'Alpha'] = model_loc.alpha_\n    \n    IX += 1 \n    \n    \ndisplay( df_stat_targetwise.sort_values('CV (mean) r2 score', ascending = False).head(20) )\nscore_loc = get_score(Y_true,Y_blend)\nprint(); print('Ridge blend targetwisely score:', '%.6f'%score_loc )\nprint('%.6f'%score_loc)\nprint(score_loc)\n","metadata":{"execution":{"iopub.status.busy":"2022-11-27T20:55:25.188831Z","iopub.execute_input":"2022-11-27T20:55:25.189622Z","iopub.status.idle":"2022-11-27T20:59:31.453336Z","shell.execute_reply.started":"2022-11-27T20:55:25.189583Z","shell.execute_reply":"2022-11-27T20:59:31.451854Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"display( df_stat_targetwise.sort_values('CV (mean) r2 score', ascending = False).head(20) )\nscore_loc = get_score(Y_true,Y_blend)\nprint(); print('Ridge blend targetwisely score:', '%.6f'%score_loc )\nprint('%.6f'%score_loc)\nprint(score_loc)\n","metadata":{"execution":{"iopub.status.busy":"2022-11-27T20:59:31.455105Z","iopub.execute_input":"2022-11-27T20:59:31.455591Z","iopub.status.idle":"2022-11-27T20:59:32.204238Z","shell.execute_reply.started":"2022-11-27T20:59:31.455543Z","shell.execute_reply":"2022-11-27T20:59:32.202851Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"r2_score(Y_true,Y_blend) , mean_squared_error(Y_true,Y_blend)\nprint( 'r2_score, mse: ')\nprint( r2_score(Y_true,Y_blend) , mean_squared_error(Y_true,Y_blend) )\nprint( '%.3f'%r2_score(Y_true,Y_blend) , '%.3f'%mean_squared_error(Y_true,Y_blend) )    ","metadata":{"execution":{"iopub.status.busy":"2022-11-27T20:59:32.206089Z","iopub.execute_input":"2022-11-27T20:59:32.206444Z","iopub.status.idle":"2022-11-27T20:59:32.919392Z","shell.execute_reply.started":"2022-11-27T20:59:32.206412Z","shell.execute_reply":"2022-11-27T20:59:32.918424Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nscore_loc = get_score(Y_true,Y_blend)\nprint(); print('Best score:')\nprint(score_loc)\nprint('%.6f'%score_loc)\n","metadata":{"execution":{"iopub.status.busy":"2022-11-27T20:59:32.920854Z","iopub.execute_input":"2022-11-27T20:59:32.921222Z","iopub.status.idle":"2022-11-27T20:59:33.672334Z","shell.execute_reply.started":"2022-11-27T20:59:32.921189Z","shell.execute_reply":"2022-11-27T20:59:33.671052Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.hist(df_stat_targetwise['Train r2 score'],bins = 20 )\nplt.hist(df_stat_targetwise['OOF r2 score'],bins = 20 )\nplt.hist(df_stat_targetwise['CV (mean) r2 score'],bins = 20 )\n\n\nplt.title('r2 scores')\nplt.show()\ndf_stat_targetwise.describe()","metadata":{"execution":{"iopub.status.busy":"2022-11-27T20:59:33.674064Z","iopub.execute_input":"2022-11-27T20:59:33.674444Z","iopub.status.idle":"2022-11-27T20:59:34.009377Z","shell.execute_reply.started":"2022-11-27T20:59:33.674410Z","shell.execute_reply":"2022-11-27T20:59:34.008248Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Stacking - LightGBM ","metadata":{}},{"cell_type":"code","source":"%%time\n\nif flag_run_stacking_by_lgb:\n    import lightgbm as lgb\n\n    # Best params found for targets prediction - not for blend\n    # They give slightly worse results for blend, and work much longer - about 40 minutes \n    #     params = {'n_estimators': 500, 'reg_alpha': 6.734991732483901, 'reg_lambda': 3.0151837674804667, \n    #               'colsample_bytree': 0.9, 'subsample': 1.0, 'max_depth': 6, 'learning_rate': 0.03814860248977263, \n    #               'num_leaves': 859, 'min_child_samples': 13,\n    #              'random_state':0,# random.randint(0,10_000)#117\n    #              }\n    params = {}\n\n    from sklearn.linear_model import RidgeCV\n    from sklearn.model_selection import cross_val_predict\n    from sklearn.model_selection import KFold\n    from sklearn.metrics import r2_score\n    from sklearn.model_selection import cross_val_score\n    from sklearn.linear_model import Ridge\n    import time\n\n\n    kf = KFold(n_splits=3, random_state = 0, shuffle = True)\n\n\n    n_targets = len(list_all_targets)\n    n_models = len( df_stat_sorted['Model']  )\n    print('n_samples_train', n_samples_train ) # \n\n    list_models_ids = list( df_stat_sorted['Model'] ) \n\n    df_data_from_all_models_for_one_target = pd.DataFrame(np.zeros( (n_samples_train,  n_models)  ), columns =  list_models_ids   )\n\n    df_stat_targetwise = pd.DataFrame(); IX = 0 \n    Y_blend = np.zeros(Y_true.shape,dtype = float )   \n\n    t0 = time.time()\n    for target_id in range(n_targets):\n        if (target_id %5) == 0: print('target N:', target_id, 'Secs passed: %.1f'%(time.time()-t0 ))\n        for i,model_id in  enumerate(df_stat_sorted['Model']):\n            data_from_current_model = dict_df[model_id].iloc[:,target_id]\n            df_data_from_all_models_for_one_target.iloc[:,i] = data_from_current_model.values\n        Y_true_for_one_target  = Y_true[:, target_id] \n        #model_loc = RidgeCV(alphas=[1e-1, 1, 1e1, 1e2, 1e3, 1e4]).fit(df_data_from_all_models_for_one_target, Y_true_for_one_target)\n\n        model_loc = lgb.LGBMRegressor(**params ).fit(df_data_from_all_models_for_one_target, Y_true_for_one_target)\n\n\n        model_loc2 = lgb.LGBMRegressor(**params )\n\n        Y_blend[:,target_id] = cross_val_predict(model_loc2  ,df_data_from_all_models_for_one_target, Y_true_for_one_target , cv = kf  )\n\n        s_oof =  r2_score( Y_true_for_one_target, Y_blend[:,target_id] )\n        list_cv_scores = cross_val_score(model_loc2, df_data_from_all_models_for_one_target, Y_true_for_one_target, cv=kf)\n\n        s = model_loc.score(df_data_from_all_models_for_one_target, Y_true_for_one_target)\n        #print(list_all_targets[target_id], 'target_id', target_id, model_loc.alpha_, 'r2 score %0.3f'%s )\n        df_stat_targetwise.loc[IX,'Target'] = list_all_targets[target_id]\n        df_stat_targetwise.loc[IX,'OOF r2 score'] = s_oof\n        df_stat_targetwise.loc[IX,'CV (mean) r2 score'] = np.mean(list_cv_scores)\n        df_stat_targetwise.loc[IX,'Train r2 score'] = s\n        #df_stat_targetwise.loc[IX,'Alpha'] = model_loc.alpha_\n\n        IX += 1 \n\n\n    print(str(model_loc))\n    display( df_stat_targetwise.sort_values('CV (mean) r2 score', ascending = False).head(20) )\n    score_loc = get_score(Y_true,Y_blend)\n    print(); print('Best score:')\n    print(score_loc)\n    print('%.6f'%score_loc)\n","metadata":{"execution":{"iopub.status.busy":"2022-11-27T20:59:34.011027Z","iopub.execute_input":"2022-11-27T20:59:34.011344Z","iopub.status.idle":"2022-11-27T21:26:52.778900Z","shell.execute_reply.started":"2022-11-27T20:59:34.011315Z","shell.execute_reply":"2022-11-27T21:26:52.777618Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nif flag_run_stacking_by_lgb: \n    score_loc = get_score(Y_true,Y_blend)\n    print(); print('LGB targetwisely stacking score:', '%.6f'%score_loc)\n    print('%.6f'%score_loc)\n    print(score_loc)\n","metadata":{"execution":{"iopub.status.busy":"2022-11-27T21:26:52.781150Z","iopub.execute_input":"2022-11-27T21:26:52.781626Z","iopub.status.idle":"2022-11-27T21:26:53.513729Z","shell.execute_reply.started":"2022-11-27T21:26:52.781582Z","shell.execute_reply":"2022-11-27T21:26:53.512595Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if flag_run_stacking_by_lgb: \n    print( r2_score(Y_true,Y_blend) , mean_squared_error(Y_true,Y_blend) )\n    print( '%.3f'%r2_score(Y_true,Y_blend) , '%.3f'%mean_squared_error(Y_true,Y_blend) )    ","metadata":{"execution":{"iopub.status.busy":"2022-11-27T21:26:53.515224Z","iopub.execute_input":"2022-11-27T21:26:53.515573Z","iopub.status.idle":"2022-11-27T21:26:53.982506Z","shell.execute_reply.started":"2022-11-27T21:26:53.515541Z","shell.execute_reply":"2022-11-27T21:26:53.981249Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n# Results for initial group of 15 LGB models (not LGBweak )\n#\n# ====== WITH rescaling Y_OOF to mean=0, and std=1\n\n# Ridge: 0.894282\n# LGB:   0.894776\n\n# r2, mse\n# Ridge:  0.173 0.200\n# LGB:    0.171 0.199\n\n# ====== WITHOUT rescaling Y_OOF to mean=0, and std=1\n# Ridge: 0.8943565503334199\n# LGB:   0.8949053501818368 - default params\n# LGB2   0.8947510476432431 - our distinguished params  - worse \n\n# r2, mse\n# Ridge: (0.17464035399184563, 0.19966963915552022)\n# LGB2:  (0.17077422156026673, 0.19919473803475118)\n","metadata":{"execution":{"iopub.status.busy":"2022-11-27T21:26:53.983833Z","iopub.execute_input":"2022-11-27T21:26:53.984264Z","iopub.status.idle":"2022-11-27T21:26:53.990119Z","shell.execute_reply.started":"2022-11-27T21:26:53.984229Z","shell.execute_reply":"2022-11-27T21:26:53.989139Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if flag_run_stacking_by_lgb: \n    plt.hist(df_stat_targetwise['Train r2 score'],bins = 20 )\n    plt.hist(df_stat_targetwise['OOF r2 score'],bins = 20 )\n    plt.hist(df_stat_targetwise['CV (mean) r2 score'],bins = 20 )\n\n\n    plt.title('r2 scores')\n    plt.show()\n    df_stat_targetwise.describe()","metadata":{"execution":{"iopub.status.busy":"2022-11-27T21:26:53.991618Z","iopub.execute_input":"2022-11-27T21:26:53.992014Z","iopub.status.idle":"2022-11-27T21:26:54.293204Z","shell.execute_reply.started":"2022-11-27T21:26:53.991981Z","shell.execute_reply":"2022-11-27T21:26:54.291972Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('%.1f seconds passed total '%(time.time()-t0start) )","metadata":{"execution":{"iopub.status.busy":"2022-11-27T21:26:54.294968Z","iopub.execute_input":"2022-11-27T21:26:54.295806Z","iopub.status.idle":"2022-11-27T21:26:54.301950Z","shell.execute_reply.started":"2022-11-27T21:26:54.295734Z","shell.execute_reply":"2022-11-27T21:26:54.300527Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}