{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":59094,"databundleVersionId":6541963,"sourceType":"competition"},{"sourceId":6675604,"sourceType":"datasetVersion","datasetId":3841755}],"dockerImageVersionId":30558,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# What is about ?   \n\nSome modeling, CV, params tuning \n\nTake a look on https://www.kaggle.com/code/alexandervc/op2-eda-baseline-s for EDA, baselines \n\nUse CV scheme basically by cell type, more precisely the one proposed by AmbrosM \nPlease upvote: https://www.kaggle.com/code/ambrosm/scp-quickstart?scriptVersionId=144293041&cellId=8\n\nThread for CV/LB modeling scores: https://www.kaggle.com/competitions/open-problems-single-cell-perturbations/discussion/444494\n\nThe notebook planned to be discussed at the webinar: https://www.kaggle.com/competitions/open-problems-single-cell-perturbations/discussion/444825\n\n\n###  Versions and Notes \n\n\n\n    77 blending models over different configs \n\n    76  ensemble model from v75 \n    \n    75 LB0.617 CV 0.993 - Ridge on many TE features tsvd features. Only compound. CV is better than top before (CV0.996639), but LB is worse than LB0.612 . \n\n    74 cosmetic changes \n    \n    73 LB0.605 Added ensembling with aggregates by compound and cell type separetly at the final submission step (following proposals by ZXMKCD, LIUDA CHELDIEVA)\n\n    70 LB0.729 CV 1.036 - LGB optimized params. CV is not so bad, but LB is awful - much worse than just mean or zero predicts (0.664)\n\n    63,64,65 LightGBM optuna params search trials 100 , 66 - 200, 67,68-500 69 - 1000\n        CV 1.036 - best found, it is better than just mean - so better than Catboost , but we optimized more params than for catboost \n        100 trials - 20 minutes, 1000 trials - 2h40min  \n        \n    62 LB0.664 CV 1.106 - Catboost - EXACTLY the same CV and LB as for mean prediction. That means model have not learn anything and just predicted by mean. \n    \n    61 cosmetic changes\n    \n    56,57,58 Catboost 100 trials of optuna  59 - 500 trials, 60 - 1000 trials (2h 40min)\n        Outcomes: (in short - bad results)\n        Got best CV score 1.106 - which is much worse 1.03 we get with Ridge with only mild tuning. And worse 0.996 - for optimized Ridge+TE. \n        The best params: 'n_estimators': 1, 'depth': 1 which basically means it is trivial model.\n        So: Catboost directly on cat features seems to be BAD approach, we need to check tuning some other params, but hardly it will give uplift. \n    \n    55 Catboost \n       Simple modeling\n       Optuna parameter optimization. \n    \n    54 LB0.612 CV0.996639 - same as v53, but blending 50 random folds on submission, (model - new params - alpha for Ridge and smoothign for TE) \n    53 LB0.612 CV0.996639 - use newly found params: alphas for Ridge (V47) and smoothings for Target Encoder (V51,52) \n    \n    52 Search for best smoothing for compound  - component-dependent , n_components = 25\n    \n    51  Search for best smoothing for cell-type - component-dependent , n_components = 25\n    \n    48 BUG - same param as v43 - so got same result. LB 0.615 CV0.997829  New alpha params found in 48\n    \n    47 Search again for best alpha with newly found smoothings \n    \n    45 LB 0.614 CV 0.999888 - same model as V43, but submission scheme changed a bit:  50 random folds \n    \n    44 BUG Search again for best alpha with newly found smoothings \n    \n    43 LB 0.615 CV 0.999888 with optimal smoothers for each component separate  \n    \n    42  Search for best smoothing for cell-type - component-dependent , n_components = 25\n    \n    41 Search for best smoothing for compound  - component-dependent , n_components = 25\n    \n    35-40 - bugs\n    \n    34 quick save\n    \n    33 LB 0.625 - submit same model but train on entire submit part ONCE - without 50 folds. \n\n    32 LB 0.624\n        CV 1.023748 - use best found alpha for each component found before - V30,31 \n        Added inference section for submit. Blending 50 random folds. Like in V2 here.  \n        \n    \n    30,31 run extensive search of best alpha for each component alpha from 1e0 to 1e8 . Best CV: 1.030232. About 40 minutes \n    \n    28  Try different alpha for different TSVD components - got small improvement: CV 1.030249\n        Conclusion: We need smaller alpha for higher components starting from 8 - 10 component - alpha 1e4 works better than 1e6. \n       \n        Added r2 score computation.\n        Conclusion (strange): r2 and mmse are inconsistent (!) - negatively correlated and better r2 does not mean better mmse and vice versa.  \n        \n\n    So current best  CV 1.030952 (worse than AmbrosM) with params:\n    alpha = 1e6,\n    smoothing: (1e7,1e7) - same for cell_type and compound -  Checked other combinations - smaller give worse result, greater - same\n    tsvd = 25\n    \n    V27 TE-smoothing params different for cell_type/compound - does not see improvement \n\n    V26 CV 1.030952 alpha1_000_000.0 smooth1000000000 tsvd25\n\n    V25 complete rewritten from scratch - previous code mostly deleted\n        Use CV scheme basically by cell type, more precisely the one proposed by AmbrosM \n        Please upvote: https://www.kaggle.com/code/ambrosm/scp-quickstart?scriptVersionId=144293041&cellId=8\n        Make a first tuning for alpha in smooth parameter \n        alpha1000000.0 smooth10000000 - got as the best for tsvd35 Ridge TE\n\n    Ridge on each target, with target encoding both features,   with preliminary denoising by tsvd \n    Shows good local results (better than all previous) but bad LB score. May be CV/scroring scheme should be improved.  \n    V23 1.788  RidgeAlpha20  100tsvd  30folds\n    V22 1.7104 RidgeAlpha20   20tsvd  30folds LB 0.672 \n    V21 1.7506 RidgeAlpha100  20tsvd  30folds\n    V20 1.781  RidgeAlpha100  35tsvd  30folds\n    V19 1.7428 RidgeAlpha20   35tstv  30folds\n\n    \n    V9 SVR+TE worse than Ridge,LGB (local), may be tune more(?)\n    V8 LGB - with categorical/ordinal - quite worse, with TE similar to Ridge, but a bit worse on 50 folds, (while a bit better on 3 folds), may be need to tune more - todo(?)\n\n    V7 Catboost - local - not good results (yet?) , neither train score, neither test score\n\n    V5 LB 0.616 - same but 614 folds i.e. leave one out validation\n    V4 LB 0.616 - same as V2, (50 folds), but bug corrected - now reducer.inverse_transform is done on each fold - as should be \n    V3 LB 0.617 - same with 50 folds, condition - use only r2_test >0.01 models - does not help \n    V2 LB 0.617 - same with 50 folds\n    V1 LB 0.623 - target encoding compound, cell type NOT used, Ridge alpha =   100000, 3 folds","metadata":{}},{"cell_type":"markdown","source":"# Preliminaries and data load ","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 time\nt0start = time.time() \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\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\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\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":"2023-10-12T11:02:44.668418Z","iopub.execute_input":"2023-10-12T11:02:44.668684Z","iopub.status.idle":"2023-10-12T11:02:46.464348Z","shell.execute_reply.started":"2023-10-12T11:02:44.668661Z","shell.execute_reply":"2023-10-12T11:02:46.462990Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load train data","metadata":{}},{"cell_type":"code","source":"%%time\nfn = '/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet'\ndf_de_train = pd.read_parquet(fn)# , index_col = 0)\nprint(df_de_train.shape)\ndf_de_train","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:02:46.466777Z","iopub.execute_input":"2023-10-12T11:02:46.467605Z","iopub.status.idle":"2023-10-12T11:02:49.246096Z","shell.execute_reply.started":"2023-10-12T11:02:46.467562Z","shell.execute_reply":"2023-10-12T11:02:49.245115Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(df_de_train['sm_name'].nunique() )\ndf_de_train['cell_type'].value_counts() ","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:02:49.247443Z","iopub.execute_input":"2023-10-12T11:02:49.248017Z","iopub.status.idle":"2023-10-12T11:02:49.259804Z","shell.execute_reply.started":"2023-10-12T11:02:49.247993Z","shell.execute_reply":"2023-10-12T11:02:49.258699Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load test and sample submission","metadata":{}},{"cell_type":"code","source":"%%time\nfn = '/kaggle/input/open-problems-single-cell-perturbations/id_map.csv'\ndf_id_map = pd.read_csv(fn)\nprint(df_id_map.shape)\ndisplay(df_id_map)\nfn = '/kaggle/input/open-problems-single-cell-perturbations/sample_submission.csv'\ndf_sample_submit = pd.read_csv(fn, index_col = 0)\nprint(df_sample_submit.shape)\ndisplay( df_sample_submit )\n","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:02:49.262182Z","iopub.execute_input":"2023-10-12T11:02:49.262901Z","iopub.status.idle":"2023-10-12T11:02:51.805429Z","shell.execute_reply.started":"2023-10-12T11:02:49.262878Z","shell.execute_reply":"2023-10-12T11:02:51.804528Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_id_map['cell_type'].value_counts() ","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:02:51.807091Z","iopub.execute_input":"2023-10-12T11:02:51.807431Z","iopub.status.idle":"2023-10-12T11:02:51.814404Z","shell.execute_reply.started":"2023-10-12T11:02:51.807397Z","shell.execute_reply":"2023-10-12T11:02:51.813544Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# CV scheme \n\nFollowing Ambros M.:  https://www.kaggle.com/code/ambrosm/scp-quickstart?scriptVersionId=144293041&cellId=8\n\nSee also discussion: https://www.kaggle.com/competitions/open-problems-single-cell-perturbations/discussion/444494\n\nThere is alternative proposal by MT: https://www.kaggle.com/code/masato114/scp-quickstart-another-cv-strategy/notebook\n\nThat proposal puts in validation same cell types as on LB, but different compounds.  While AmbrosM proposals puts same compounds, but different cell types. \n\n","metadata":{}},{"cell_type":"markdown","source":"## Preliminaries for CV","metadata":{}},{"cell_type":"code","source":"temp = df_de_train.groupby(['cell_type']).size() \ntemp","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:02:51.815665Z","iopub.execute_input":"2023-10-12T11:02:51.816502Z","iopub.status.idle":"2023-10-12T11:02:51.827232Z","shell.execute_reply.started":"2023-10-12T11:02:51.816467Z","shell.execute_reply":"2023-10-12T11:02:51.826396Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"temp = df_de_train.groupby(['cell_type']).size() > 20\ntemp","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:02:51.828195Z","iopub.execute_input":"2023-10-12T11:02:51.828478Z","iopub.status.idle":"2023-10-12T11:02:51.838512Z","shell.execute_reply.started":"2023-10-12T11:02:51.828458Z","shell.execute_reply":"2023-10-12T11:02:51.837973Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"validation_cell_types = temp[temp].index # 4 cell types\nvalidation_cell_types","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:02:51.839324Z","iopub.execute_input":"2023-10-12T11:02:51.839541Z","iopub.status.idle":"2023-10-12T11:02:51.850590Z","shell.execute_reply.started":"2023-10-12T11:02:51.839522Z","shell.execute_reply":"2023-10-12T11:02:51.850054Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_sm_names = df_de_train.query(\"cell_type == 'B cells'\").sm_name.values # 17 compounds including the two control compounds\ntrain_sm_names","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:02:51.851390Z","iopub.execute_input":"2023-10-12T11:02:51.851742Z","iopub.status.idle":"2023-10-12T11:02:52.790639Z","shell.execute_reply.started":"2023-10-12T11:02:51.851721Z","shell.execute_reply":"2023-10-12T11:02:52.789754Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## CV scheme itself. Mean,Median,etc... predictions trivial baselines\n\n\nLB scores for submissions by mean, median, percentiles can be found here:\nhttps://www.kaggle.com/code/alexandervc/op2-eda-baseline-s?scriptVersionId=143396607&cellId=1\n\n\n    CV 1.057 V6 LB 0.666 quantile(0.7)\n    CV 1.008 V5 LB 0.657 quantile(0.6) \n            results are again better, that probably indicates some shift between train and public data \n    CV  0.997 V4 LB 0.659 - median instead of mean - results are a bit better, \n            it might mean either a bit of presense of outliers, or  public is somewhat different from train - next experiments suggests second is true \n    CV 1.106  V1 LB 0.664 - submission of train means - the simplest baseline \n    \n    Not a perfect match of CV to LB but not terrible miscorrepondence\n    \nOverall situation is strange CV of median is better than for Ridge+TE, but LB is much worse. \n","metadata":{}},{"cell_type":"code","source":"%%time\nY = df_de_train.iloc[:,5:].values\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.metrics import r2_score\nverbose = 100\nfor predict_method in ['mean', 'median', 'percentile60', 'percentile70']:\n    print( predict_method )\n    mrrmse_list = [] ; r2_list = []\n    for fold, val_cell_type in enumerate(validation_cell_types):\n        mask_va = (df_de_train.cell_type == val_cell_type) & ~df_de_train.sm_name.isin(train_sm_names)\n        mask_tr = ~mask_va # 485 or 487 training rows\n        print('fold:', fold, val_cell_type, mask_va.sum() , mask_tr.sum(), mask_va.sum()+ mask_tr.sum() ,  )\n\n        X_train_categorical = df_de_train[['cell_type','sm_name']][mask_tr]\n        X_valid_categorical = df_de_train[['cell_type','sm_name']][mask_va]    \n        Y_train = Y[mask_tr,:]\n        Y_valid = Y[mask_va,:]    \n        print( 'X_valid_categorical.shape', X_valid_categorical.shape , 'X_train_categorical.shape',X_train_categorical.shape, 'Y_train.shape', Y_train.shape, 'Y_valid.shape', Y_valid.shape  )\n\n        if predict_method == 'mean':\n            vec1 = np.mean( Y_train, axis = 0)\n        elif predict_method == 'median':\n            vec1 = np.median( Y_train, axis = 0)\n        elif predict_method == 'percentile60':\n            vec1 = np.percentile( Y_train,60, axis = 0)\n        elif predict_method == 'percentile70':\n            vec1 = np.percentile( Y_train,70, axis = 0)\n\n        Y_pred_valid = np.zeros_like( Y_valid) + vec1[ np.newaxis,:] # broadcasting vec1 to Y_pred sizes \n        mrrmse = np.sqrt(np.square(Y_pred_valid - Y_valid).mean(axis=1)).mean()\n        r2 = r2_score( Y_valid , Y_pred_valid ) \n        print(f\"# Fold {fold}: mse: {mrrmse:5.3f} r2: {r2:5.3f} val='{val_cell_type}'\")    \n        mrrmse_list.append( mrrmse )\n        r2_list.append( r2 )\n        print()\n\n    mrrmse = np.array(mrrmse_list).mean()    \n    r2 = np.array(r2_list).mean()    \n    if verbose >= 10:\n        print(f\"# {predict_method} Overall mrrmse: {mrrmse:5.3f}, r2: {r2:5.3f}\" )\n    print()\n    print('--------------------------------------------------------')\n    ","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:02:52.793307Z","iopub.execute_input":"2023-10-12T11:02:52.793552Z","iopub.status.idle":"2023-10-12T11:02:57.187044Z","shell.execute_reply.started":"2023-10-12T11:02:52.793532Z","shell.execute_reply":"2023-10-12T11:02:57.186026Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Catboost model\n\nCatboost allows to work with categorical features as it is . Here is CatBoost model with preliminary reduction of target by TSVD.\n\nResults are quite poor up to now - the best optuna 1000 rounds can find - params which really do not learn anything, but predict by average. \n\n\nMay be we need to tune more params not only three as it is done below. ","metadata":{}},{"cell_type":"markdown","source":"## Preliminaries","metadata":{}},{"cell_type":"code","source":"%%time\n\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.metrics import r2_score\nfrom sklearn.decomposition import TruncatedSVD\n\nn_components = 25\nreducer = TruncatedSVD(n_components=n_components, n_iter=7, random_state=42)\n\n","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:02:57.188379Z","iopub.execute_input":"2023-10-12T11:02:57.188770Z","iopub.status.idle":"2023-10-12T11:02:57.338361Z","shell.execute_reply.started":"2023-10-12T11:02:57.188738Z","shell.execute_reply":"2023-10-12T11:02:57.337458Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Catboost with CV. Simple example  ","metadata":{}},{"cell_type":"code","source":"%%time\n\nverbose = 10000\nimport catboost\nfrom catboost import CatBoostRegressor, Pool\nfrom sklearn.multioutput import MultiOutputRegressor\ncategorical_features = ['cell_type','sm_name']\nmodel = MultiOutputRegressor( CatBoostRegressor(cat_features=categorical_features, verbose = 0,  # Categorical features ) ) \n                        iterations=1,  # Number of boosting iterations\n                          depth=1,        # Depth of the trees\n                          learning_rate=0.01,  # Learning rate\n                          loss_function='RMSE'))  # Specify your loss function (e.g., RMSE for regression)\n#                           verbose=0) ) # Set verbose to 0 to suppress output\n\nmrrmse_list = []\nfor fold, val_cell_type in enumerate(validation_cell_types):\n    mask_va = (df_de_train.cell_type == val_cell_type) & ~df_de_train.sm_name.isin(train_sm_names)\n    mask_tr = ~mask_va # 485 or 487 training rows\n    print(mask_va.sum() , mask_tr.sum(), mask_va.sum()+ mask_tr.sum() , 'fold:', fold, val_cell_type )\n    \n\n    Y_red_train = reducer.fit_transform(Y[mask_tr,:])\n    X_train_categorical = df_de_train[['cell_type','sm_name']][mask_tr]\n    X_valid_categorical = df_de_train[['cell_type','sm_name']][mask_va]    \n    #Yr_valid = reducer.transform(Y[mask_va,:])\n    if verbose >= 1000:        print('Y_red_train.shape',Y_red_train.shape, 'X_valid_categorical.shape', X_valid_categorical.shape , 'X_train_categorical.shape',X_train_categorical.shape )\n\n    model.fit(X_train_categorical, Y_red_train )\n    Y_red_valid_pred = model.predict(X_valid_categorical) # \n    Y_valid_pred = reducer.inverse_transform( Y_red_valid_pred )\n    \n    mrrmse = np.sqrt(np.square(Y[mask_va,:] - Y_valid_pred).mean(axis=1)).mean(); mrrmse_list.append(mrrmse) \n    r2 = r2_score( Y[mask_va,:] ,Y_valid_pred  ); r2_list.append(r2)\n    if verbose >= 100:\n        print(f\"# Fold {fold}: mse: {mrrmse:5.3f} r2: {r2:5.3f} val='{val_cell_type}'\")\n    if verbose >= 10: print()\n        \nprint(model); print()        \nmrrmse = np.array(mrrmse_list).mean()    \nif verbose >= 10:\n    print(f\"# Overall {mrrmse:5.3f}\" )\n","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:02:57.339589Z","iopub.execute_input":"2023-10-12T11:02:57.340197Z","iopub.status.idle":"2023-10-12T11:03:06.701294Z","shell.execute_reply.started":"2023-10-12T11:02:57.340162Z","shell.execute_reply":"2023-10-12T11:03:06.699988Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"##  Prepare submission. Retrain on full data","metadata":{}},{"cell_type":"code","source":"%%time\ndf_de_train = pd.read_parquet('/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet')\ndf_id_map = pd.read_csv('/kaggle/input/open-problems-single-cell-perturbations/id_map.csv')\nfrom catboost import CatBoostRegressor, Pool\nfrom sklearn.multioutput import MultiOutputRegressor\nfrom sklearn.decomposition import TruncatedSVD\n\nmodel = MultiOutputRegressor( CatBoostRegressor(cat_features=['cell_type','sm_name'], verbose = 0,  # Categorical features ) ) \n                        iterations=1,  # Number of boosting iterations\n                          depth=1,        # Depth of the trees\n                          learning_rate=0.01,  # Learning rate\n                          loss_function='RMSE'))  # Specify your loss function (e.g., RMSE for regression)\nX_train_categorical = df_de_train[['cell_type','sm_name']]\nY = df_de_train.iloc[:,5:].values\nn_components = 25; reducer = TruncatedSVD(n_components=n_components, n_iter=7, random_state=42)\nY_red_train = reducer.fit_transform(Y)\nmodel.fit(X_train_categorical, Y_red_train )\n\nX_submit_categorical = df_id_map[['cell_type','sm_name']]\nY_red_submit = model.predict(X_submit_categorical )\nY_submit = reducer.inverse_transform(  Y_red_submit )\n\ndf_submit = pd.DataFrame(Y_submit, columns = df_de_train.columns[5:])\ndf_submit.index.name = 'id'; print('df_submit.shape', df_submit.shape )\ndisplay(df_submit.head(3))","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:03:06.702998Z","iopub.execute_input":"2023-10-12T11:03:06.703808Z","iopub.status.idle":"2023-10-12T11:03:10.316331Z","shell.execute_reply.started":"2023-10-12T11:03:06.703763Z","shell.execute_reply":"2023-10-12T11:03:10.310979Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nif 0:# flag_save_submit:\n    df_submit.to_csv('submission_catboost.csv')","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:03:10.318116Z","iopub.execute_input":"2023-10-12T11:03:10.318634Z","iopub.status.idle":"2023-10-12T11:03:10.326880Z","shell.execute_reply.started":"2023-10-12T11:03:10.318593Z","shell.execute_reply":"2023-10-12T11:03:10.325780Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Optuna param optimization for Catboost\n\nFirst optimization (1000+ iterations) lead to very bad conclusion - we get CV 1.106 LB 0.664 - exactly as for prediction by means. \nSo the model have not learned anything but just predicted by mean. \n\nMay be we need to optimize other params also but no much optimism \n\n    First optimization we tried: \n    'n_estimators':  trial.suggest_categorical('n_estimators', [1,2,5,10,20, 50, 100, 200, 500]) , # \n    'depth' : trial.suggest_int('depth', 1,16) ,\n    'learning_rate': trial.suggest_float('learning_rate',0.01, 0.1) ,  \n","metadata":{}},{"cell_type":"code","source":"%%time\nimport optuna\n\nimport catboost\nfrom catboost import CatBoostRegressor, Pool\nfrom sklearn.multioutput import MultiOutputRegressor\ncategorical_features = ['cell_type','sm_name']\n\nimport time\n\ndict_trial_count = {0:0}\ndef objective(trial):\n    verbose = 1\n    dict_trial_count[0] += 1\n    t0 = time.time()\n\n    params = {\n        'random_state' :  0,\n        'verbose' :  0, \n        'cat_features': categorical_features, \n        'n_estimators':  trial.suggest_categorical('n_estimators', [1,2,5,10,20, 50, 100, 200, 500]) , # \n        'depth' : trial.suggest_int('depth', 1,16) ,\n        'learning_rate': trial.suggest_float('learning_rate',0.01, 0.1) ,  \n        'loss_function': 'RMSE'\n    }\n    \n    model = MultiOutputRegressor( CatBoostRegressor(**params)) \n    \n    mrrmse_list = []\n    for fold, val_cell_type in enumerate(validation_cell_types):\n        mask_va = (df_de_train.cell_type == val_cell_type) & ~df_de_train.sm_name.isin(train_sm_names)\n        mask_tr = ~mask_va # 485 or 487 training rows\n        if verbose >= 10:\n            print(  'fold:', fold, val_cell_type, mask_va.sum() , mask_tr.sum(), mask_va.sum()+ mask_tr.sum() )\n\n\n        Yr_train = reducer.fit_transform(Y[mask_tr,:])\n        X_train_categorical = df_de_train[['cell_type','sm_name']][mask_tr]\n        X_valid_categorical = df_de_train[['cell_type','sm_name']][mask_va]    \n        #Yr_valid = reducer.transform(Y[mask_va,:])\n        if verbose >= 1000:        print('Yr_train.shape',Yr_train.shape, 'X_valid_categorical.shape', X_valid_categorical.shape , 'X_train_categorical.shape',X_train_categorical.shape )\n\n        model.fit(X_train_categorical, Yr_train )\n        Yr_valid_pred = model.predict(X_valid_categorical) # \n        Y_valid = reducer.inverse_transform( Yr_valid_pred )\n\n        mrrmse = np.sqrt(np.square(Y[mask_va,:] - Y_valid).mean(axis=1)).mean()\n        r2 = r2_score( Y[mask_va,:] ,Y_valid  )\n        if verbose >= 10:\n            print(f\"# Fold {fold}, {val_cell_type}: mse: {mrrmse:5.3f} r2: {r2:5.3f} val='{val_cell_type}'\")\n        mrrmse_list.append(mrrmse)    \n\n        if verbose >= 10: print()\n\n    \n    mrrmse = np.array(mrrmse_list).mean()    \n    if verbose >= 1:\n        print(params)\n        print(f\"# Overall {mrrmse:5.3f}\", 'trial_count = ', dict_trial_count[0] ,'timing: %.1f'%(time.time()-t0) )\n    \n    return mrrmse\n\n","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:03:10.329145Z","iopub.execute_input":"2023-10-12T11:03:10.330340Z","iopub.status.idle":"2023-10-12T11:03:11.497679Z","shell.execute_reply.started":"2023-10-12T11:03:10.330290Z","shell.execute_reply":"2023-10-12T11:03:11.496625Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nn_trials = 0\nif n_trials > 0:\n    print('Optimization starts.  n_trials = ', n_trials)\n    dict_trial_count[0] = 0\n\n\n    optuna.logging.set_verbosity(optuna.logging.WARNING)\n    import warnings\n    warnings.filterwarnings(\"ignore\", category=FutureWarning)\n\n    study = optuna.create_study()\n    study.optimize(objective, n_trials=n_trials)\n\n    # Output for best found params: \n    print(); print('Best params:')\n    print(study.best_params)  # E.g. {'x': 2.002108042}","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:03:11.498952Z","iopub.execute_input":"2023-10-12T11:03:11.499952Z","iopub.status.idle":"2023-10-12T11:03:11.507163Z","shell.execute_reply.started":"2023-10-12T11:03:11.499913Z","shell.execute_reply":"2023-10-12T11:03:11.506075Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# LightGBM\n\nLightGBM as CatBoost can support categorical features out of box without any encoding. \n\n\nOptuna hyperparameter search leads to params that better than predictions by mean:  CV 1.036. However LB score 0.729 - which is very bad - much worse than predictions by zeros/means 0.664.\n\n\n\nSubmission in notebook version: 70 LB0.729 CV 1.036 - LGB optimized params. CV is not so bad, but LB is awful - much worse than just mean or zero predicts (0.664)\n\n\nOptuna optimization in notebooks versions:     63,64,65 LightGBM optuna params search trials 100 , 66 - 200, 67,68-500 69 - 1000\n        CV 1.036 - best found, it is better than just mean - so better than Catboost , but we optimized more params than for catboost \n        100 trials - 20 minutes, 1000 trials - 2h40min  \n\n\nMay be some other params should be tuned, by we are not optimistic.\n","metadata":{}},{"cell_type":"markdown","source":"## Simple LightGBM modeling with CV score","metadata":{}},{"cell_type":"code","source":"%%time\nverbose = 100\nimport lightgbm as lgb\nfrom sklearn.multioutput import MultiOutputRegressor\nparams_best1_cv1_036 = {'random_state': 0, 'n_estimators': 20, 'reg_alpha': 6.764079452929363, 'reg_lambda': 0.41900776876588564, \n        'colsample_bytree': 0.3, 'subsample': 0.7, 'max_depth': 1, 'learning_rate': 0.08456104070184789, 'num_leaves': 682, 'min_child_samples': 102}\n\nmodel = MultiOutputRegressor( lgb.LGBMRegressor( **params_best1_cv1_036 ) )\n\nmrrmse_list = []; r2_list = []\nfor fold, val_cell_type in enumerate(validation_cell_types):\n    mask_va = (df_de_train.cell_type == val_cell_type) & ~df_de_train.sm_name.isin(train_sm_names)\n    mask_tr = ~mask_va # 485 or 487 training rows\n    print(mask_va.sum() , mask_tr.sum(), mask_va.sum()+ mask_tr.sum() , 'fold:', fold, val_cell_type )\n    \n\n    Y_red_train = reducer.fit_transform(Y[mask_tr,:])\n    X_train_categorical = df_de_train[['cell_type','sm_name']][mask_tr].astype('category')\n    X_valid_categorical = df_de_train[['cell_type','sm_name']][mask_va].astype('category')    \n    #Yr_valid = reducer.transform(Y[mask_va,:])\n    if verbose >= 1000:        print('Yr_train.shape',Y_red_train.shape, 'X_valid_categorical.shape', X_valid_categorical.shape , 'X_train_categorical.shape',X_train_categorical.shape )\n\n    model.fit(X_train_categorical, Y_red_train )\n    Y_red_valid_pred = model.predict(X_valid_categorical) # \n    Y_valid_pred = reducer.inverse_transform( Y_red_valid_pred )        \n    \n    mrrmse = np.sqrt(np.square(Y[mask_va,:] - Y_valid_pred).mean(axis=1)).mean(); mrrmse_list.append(mrrmse) \n    r2 = r2_score( Y[mask_va,:] ,Y_valid_pred  ); r2_list.append(r2)\n    if verbose >= 100:\n        print(f\"# Fold {fold}: mse: {mrrmse:5.3f} r2: {r2:5.3f} val='{val_cell_type}'\")\n    if verbose >= 10: print()\n        \nprint(model); print()        \nmrrmse = np.array(mrrmse_list).mean(); r2 =   np.array(r2_list).mean();  \nif verbose >= 10:\n    print(f\"# Overall {mrrmse:5.3f}, r2 {r2:5.3f}\" )    ","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:03:11.508550Z","iopub.execute_input":"2023-10-12T11:03:11.508828Z","iopub.status.idle":"2023-10-12T11:03:19.974230Z","shell.execute_reply.started":"2023-10-12T11:03:11.508805Z","shell.execute_reply":"2023-10-12T11:03:19.973020Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Optuna+LightGBM\n\n\n    Norebook Versions: 63,64,65 LightGBM optuna params search trials 100 , 66 - 200, 67,68-500 69 - 1000\n        CV 1.036 - best found, it is better than just mean - so better than Catboost , but we optimized more params than for catboost \n        100 trials - 20 minutes, 1000 trials - 2h40min  \n        ","metadata":{}},{"cell_type":"code","source":"%%time\nimport time\ndict_trial_count = {0:0}\ndef objective(trial):\n    verbose = 1\n    dict_trial_count[0] += 1\n    t0 = time.time()\n\n    params = {\n        'random_state' :  0,\n        'n_estimators' :  trial.suggest_categorical('n_estimators',  [1,2,5,10,20, 50, 100, 200, 500]) ,\n        \n        'reg_alpha': trial.suggest_float('reg_alpha', 1e-3, 10.0),\n        'reg_lambda': trial.suggest_float('reg_lambda', 1e-3, 10.0),\n        'colsample_bytree': trial.suggest_categorical('colsample_bytree', [0.3,0.4,0.5,0.6,0.7,0.8,0.9, 1.0]),\n        'subsample': trial.suggest_categorical('subsample', [0.4,0.5,0.6,0.7,0.8,1.0]),\n        \n        'max_depth': trial.suggest_int('max_depth', 1,16),\n        'learning_rate': trial.suggest_float('learning_rate', 0.01, 0.1),\n#        'max_leaves' : trial.suggest_int('max_leaves', 0, 1000),\n#         'min_child_samples': trial.suggest_int('min_child_samples', 1, 300),\n        'num_leaves' : trial.suggest_int('num_leaves', 2, 1000),\n        'min_child_samples': trial.suggest_int('min_child_samples', 1, 300),\n        #'cat_smooth' : trial.suggest_int('min_data_per_groups', 1, 100)        \n    }\n    \n    model = MultiOutputRegressor(  lgb.LGBMRegressor(**params) ) \n    \n    mrrmse_list = []; r2_list = []\n    for fold, val_cell_type in enumerate(validation_cell_types):\n        mask_va = (df_de_train.cell_type == val_cell_type) & ~df_de_train.sm_name.isin(train_sm_names)\n        mask_tr = ~mask_va # 485 or 487 training rows\n        if verbose >= 10:\n            print(  'fold:', fold, val_cell_type, mask_va.sum() , mask_tr.sum(), mask_va.sum()+ mask_tr.sum() )\n\n\n        Y_red_train = reducer.fit_transform(Y[mask_tr,:])\n        X_train_categorical = df_de_train[['cell_type','sm_name']][mask_tr].astype('category')\n        X_valid_categorical = df_de_train[['cell_type','sm_name']][mask_va].astype('category')    \n        #Yr_valid = reducer.transform(Y[mask_va,:])\n        if verbose >= 1000:        print('Yr_train.shape',Yr_train.shape, 'X_valid_categorical.shape', X_valid_categorical.shape , 'X_train_categorical.shape',X_train_categorical.shape )\n\n        model.fit(X_train_categorical, Y_red_train )\n        Y_red_valid_pred = model.predict(X_valid_categorical) # \n        Y_valid_pred = reducer.inverse_transform( Y_red_valid_pred )        \n\n\n        mrrmse = np.sqrt(np.square(Y[mask_va,:] - Y_valid_pred).mean(axis=1)).mean(); mrrmse_list.append(mrrmse) \n        r2 = r2_score( Y[mask_va,:] ,Y_valid_pred  ); r2_list.append(r2)\n        if verbose >= 100:\n            print(f\"# Fold {fold}: mse: {mrrmse:5.3f} r2: {r2:5.3f} val='{val_cell_type}'\")\n        if verbose >= 10: print()\n        \n\n    mrrmse = np.array(mrrmse_list).mean(); r2 =   np.array(r2_list).mean();  \n    if verbose >= 1:\n        print(params)\n        print(f\"# Overall {mrrmse:5.3f}, r2 {r2:5.3f} \", 'trial_count = ', dict_trial_count[0] ,'timing: %.1f'%(time.time()-t0) )\n    \n    return mrrmse\n\n","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:03:19.976448Z","iopub.execute_input":"2023-10-12T11:03:19.977374Z","iopub.status.idle":"2023-10-12T11:03:19.992715Z","shell.execute_reply.started":"2023-10-12T11:03:19.977326Z","shell.execute_reply":"2023-10-12T11:03:19.992125Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nn_trials = 1# 1000 - \nstudy = optuna.create_study()\n\nif n_trials > 0:\n    print('Optimization starts.  n_trials = ', n_trials)\n    dict_trial_count[0] = 0\n\n\n    optuna.logging.set_verbosity(optuna.logging.WARNING)\n    import warnings\n    warnings.filterwarnings(\"ignore\", category=FutureWarning)\n\n    study = optuna.create_study()\n    study.optimize(objective, n_trials=n_trials)\n\n    # Output for best found params: \n    print(); print('Best params:')\n    print(study.best_params)  # E.g. {'x': 2.002108042}","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:03:19.993809Z","iopub.execute_input":"2023-10-12T11:03:19.994281Z","iopub.status.idle":"2023-10-12T11:03:39.606574Z","shell.execute_reply.started":"2023-10-12T11:03:19.994258Z","shell.execute_reply":"2023-10-12T11:03:39.603203Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Prepare submission, retrain on full data. LightGBM","metadata":{}},{"cell_type":"code","source":"%%time\n\nparams_best1_cv1_036 = {'random_state': 0, 'n_estimators': 20, 'reg_alpha': 6.764079452929363, 'reg_lambda': 0.41900776876588564, 'colsample_bytree': 0.3, 'subsample': 0.7, 'max_depth': 1, 'learning_rate': 0.08456104070184789, 'num_leaves': 682, 'min_child_samples': 102}\nparams = params_best1_cv1_036\n\nflag_loc = True\n# try:\n#     params = study.best_params\n# except:\n#     flag_loc = False\n\nif flag_loc:\n    model = MultiOutputRegressor(  lgb.LGBMRegressor(**params) ) \n\n    X_train_categorical = df_de_train[['cell_type','sm_name']].astype('category')\n    Y = df_de_train.iloc[:,5:].values\n    n_components = 25; reducer = TruncatedSVD(n_components=n_components, n_iter=7, random_state=42)\n    Y_red_train = reducer.fit_transform(Y)\n    model.fit(X_train_categorical, Y_red_train )\n\n    X_submit_categorical = df_id_map[['cell_type','sm_name']].astype('category')\n    Y_red_submit = model.predict(X_submit_categorical )\n    Y_submit = reducer.inverse_transform(  Y_red_submit )\n\n    df_submit = pd.DataFrame(Y_submit, columns = df_de_train.columns[5:])\n    df_submit.index.name = 'id'; print('df_submit.shape', df_submit.shape )\n    display(df_submit.head(3))\n\n    if 1:# flag_save_submit:\n        df_submit.to_csv('submission_lgb.csv')","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:03:39.608579Z","iopub.execute_input":"2023-10-12T11:03:39.609208Z","iopub.status.idle":"2023-10-12T11:03:50.297880Z","shell.execute_reply.started":"2023-10-12T11:03:39.609164Z","shell.execute_reply":"2023-10-12T11:03:50.296943Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Ridge + TargetEncoder + tsvd\n\nOptimized values of alpha for ridge and smoothing for Target Encoder has been found and used.\n\nSee history of optimization in the first section of the present notebook or the discussion post:\nhttps://www.kaggle.com/competitions/open-problems-single-cell-perturbations/discussion/444494\n","metadata":{}},{"cell_type":"code","source":"%%time\n\nverbose = 1\nsmoothing = (1e7,1e7) # smoothing = 100\n# alpha = 1e6\nn_components = 25\nalpha_default = 1e6\nsmoothing_default = (1e7,1e7)\n\n\ndict_best_alpha_for_component = {0: 1000000.0, 1: 1000000.0, 2: 1000000.0, 3: 1000000.0, 4: 100000.0, 5: 1000000.0, 6: 1000000.0, 7: 1000000.0, 8: 100000.0, 9: 100000.0, 10: 10000.0, 11: 10000.0, 12: 100000.0, 13: 10000.0, 14: 100000.0, 15: 100000.0, 16: 100000.0, 17: 10000.0, 18: 10000.0, 19: 10000.0, 20: 10000.0, 21: 10000.0, 22: 10000.0, 23: 1000.0, 24: 10000.0}\n# From versions 31,32 of the notebook:\n# https://www.kaggle.com/code/alexandervc/op2-models-cv-tuning?scriptVersionId=144944894&cellId=19\ndict_best_alpha_for_component = {0: 1000000.0, 1: 100000.0, 2: 100000.0, 3: 1000000.0, 4: 100000.0, 5: 1000000.0, 6: 1000000.0, 7: 1000000.0, 8: 100000.0, 9: 100000.0, 10: 10000.0, 11: 10000.0, 12: 100000.0, 13: 10000.0, 14: 100000.0, 15: 100000.0, 16: 100000.0, 17: 10000.0, 18: 10000.0, 19: 10000.0, 20: 10000.0, 21: 10000.0, 22: 10000.0, 23: 1000.0, 24: 10000.0}\n# From versions 47 of the notebook:\n# https://www.kaggle.com/code/alexandervc/op2-models-cv-tuning?scriptVersionId=144999017\n    \n\ndict_best_smooth_compound_for_component = {0: 1000000000000000.0, 1: 1000000000000000.0, 2: 1000000000000000.0, 3: 1000000000000000.0, 4: 1.0, 5: 10000000000000.0, 6: 10000000000000.0, 7: 1000000000000000.0, 8: 1000000000000000.0, 9: 100.0, 10: 100.0, 11: 100.0, 12: 100.0, 13: 100.0, 14: 1000000000000000.0, 15: 1000.0, 16: 1000000000000000.0, 17: 100.0, 18: 100.0, 19: 1000000000000000.0, 20: 10000000000000.0, 21: 1000.0, 22: 10000000000000.0, 23: 100.0, 24: 1000000000000000.0}\n# Version 41: https://www.kaggle.com/code/alexandervc/op2-models-cv-tuning?scriptVersionId=144969071&cellId=20\ndict_best_smooth_compound_for_component = {0: 1000000000000000.0, 1: 1000000000000000.0, 2: 1000000000000000.0, 3: 1000000000000000.0, 4: 10.0, 5: 1000000000000000.0, 6: 10000000000000.0, 7: 1000000000000000.0, 8: 1000000000000000.0, 9: 100.0, 10: 100.0, 11: 100.0, 12: 100.0, 13: 100.0, 14: 1000000000000000.0, 15: 1000.0, 16: 1000000000000000.0, 17: 100.0, 18: 100.0, 19: 10000000000000.0, 20: 1000000000000000.0, 21: 10000000000000.0, 22: 1000000000000000.0, 23: 100.0, 24: 1000000000000000.0}\n# Version 52: https://www.kaggle.com/code/alexandervc/op2-models-cv-tuning?scriptVersionId=145015691\n\ndict_best_smooth_celltype_for_component = {0: 1.0, 1: 1.0, 2: 1.0, 3: 1000000000000000.0, 4: 10.0, 5: 100.0, 6: 1000000000000000.0, 7: 100.0, 8: 1000000000000000.0, 9: 100.0, 10: 10.0, 11: 100.0, 12: 1.0, 13: 1000.0, 14: 10000000000000.0, 15: 100.0, 16: 100.0, 17: 10.0, 18: 1.0, 19: 10000000000000.0, 20: 100.0, 21: 10.0, 22: 1.0, 23: 1.0, 24: 1.0}\n# Version 42: https://www.kaggle.com/code/alexandervc/op2-models-cv-tuning?scriptVersionId=144969161&cellId=20\ndict_best_smooth_celltype_for_component = {0: 0.0, 1: 1.0, 2: 1.0, 3: 1000000000000000.0, 4: 10.0, 5: 1000.0, 6: 1000000000000000.0, 7: 100.0, 8: 1000000000000000.0, 9: 100.0, 10: 10.0, 11: 100.0, 12: 0.0, 13: 1000.0, 14: 10000000000000.0, 15: 100.0, 16: 100.0, 17: 10.0, 18: 1.0, 19: 10000000000000.0, 20: 100.0, 21: 10.0, 22: 0.0, 23: 1.0, 24: 0.0}\n# Version 51: https://www.kaggle.com/code/alexandervc/op2-models-cv-tuning?scriptVersionId=145015609\n\ndf_stat = pd.DataFrame(); IX_stat = -1\n\nfrom sklearn.linear_model import Ridge\nimport time\nfrom sklearn.metrics import r2_score\nfrom sklearn.decomposition import TruncatedSVD\n\n# from sklearn.preprocessing import TargetEncoder\nimport category_encoders as ce\n\n\nY = df_de_train.iloc[:,5:].values\n\ncc = 0\nt0 = time.time()\n# for n_components in [25]:\nif 1:\n    if 1:\n        if 1:\n# for i_selected_component in range( n_components):#  [1,2,3]:\n#     for smoothing in [ (0,1e7), (1e0,1e7), (1e1,1e7),(1e2,1e7),(1e3,1e7),(1e4,1e7),  (1e5,1e7), (1e6,1e7), (1e7,1e7),(1e10,1e7), (1e13,1e7),(1e15,1e7) ]: # Checked other combinations - smaller give worse result, greater - same\n#         if 1:\n#         for alpha_for_selected_component in [1,1e1 ,1e2, 1e3,1e4,1e5,1e6,1e7,1e8]: # [1e3, 1e4,1e5,1e6,1e7]:\n            smoothing = (smoothing[1], smoothing[0])  \n            reducer = TruncatedSVD(n_components=n_components, n_iter=7, random_state=42)\n\n            mrrmse_list = []\n\n            for fold, val_cell_type in enumerate(validation_cell_types):\n                t0_fold = time.time()\n\n                mask_va = (df_de_train.cell_type == val_cell_type) & ~df_de_train.sm_name.isin(train_sm_names)\n                mask_tr = ~mask_va # 485 or 487 training rows\n                if verbose >= 1000:\n                    print(mask_va.sum() , mask_tr.sum(), mask_va.sum()+ mask_tr.sum() , 'fold:', fold, val_cell_type )\n\n                Yr_train = reducer.fit_transform(Y[mask_tr,:])\n                if verbose >= 1000:\n                    print('Yr_train.shape',Yr_train.shape)\n                X_train_categorical = df_de_train[['cell_type','sm_name']][mask_tr]\n\n                Yr_valid = reducer.transform(Y[mask_va,:])\n                X_valid_categorical = df_de_train[['cell_type','sm_name']][mask_va]\n\n#                 enc_0 = ce.TargetEncoder(smoothing=smoothing[0])\n#                 enc_1 = ce.TargetEncoder(smoothing=smoothing[1])\n                Yr_valid_pred = np.zeros( Yr_valid.shape )\n                for i_comp in range(Yr_train.shape[1]):\n                    \n#                     if i_comp == i_selected_component:\n#                         #model = Ridge(alpha = alpha_for_selected_component )\n#                         enc_0 = ce.TargetEncoder(smoothing=smoothing[0])\n#                         enc_1 = ce.TargetEncoder(smoothing=smoothing[1])\n#                     else:\n#                         #model = Ridge(alpha = dict_best_alpha_for_component[i_comp] )\n#                         enc_0 = ce.TargetEncoder(smoothing=dict_best_smooth_celltype_for_component[i_comp])\n#                         enc_1 = ce.TargetEncoder(smoothing=dict_best_smooth_compound_for_component[i_comp])\n\n                    enc_0 = ce.TargetEncoder(smoothing=dict_best_smooth_celltype_for_component[i_comp]) # Best smooth0 found in v42\n                    enc_1 = ce.TargetEncoder(smoothing=dict_best_smooth_compound_for_component[i_comp]) # Best smooth1 found in v42\n                    model = Ridge(alpha = dict_best_alpha_for_component[i_comp] ) # Best alpha found in v 30-31 , new v47\n    \n                    y_tr = Yr_train[:,i_comp]\n                    X_train_enc = pd.concat( [enc_0.fit_transform(X_train_categorical[['cell_type']], y_tr),\n                                              enc_1.fit_transform(X_train_categorical[['sm_name']], y_tr) ] , axis = 1)\n                    #print('X_train_enc.shape', X_train_enc.shape)\n                    model.fit(X_train_enc, y_tr)\n\n                    #X_valid_enc = enc.transform(X_valid_categorical)\n                    X_valid_enc = pd.concat( [enc_0.transform(X_valid_categorical[['cell_type']] ),\n                                              enc_1.transform(X_valid_categorical[['sm_name']] ) ] , axis = 1 )\n                    y_pred = model.predict(X_valid_enc)\n                    Yr_valid_pred[:,i_comp] = y_pred\n\n                Y_valid = reducer.inverse_transform( Yr_valid_pred  )\n\n                mrrmse = np.sqrt(np.square(Y[mask_va,:] - Y_valid).mean(axis=1)).mean()\n                r2 = r2_score( Y[mask_va,:] ,Y_valid  )\n                if verbose >= 100:\n                    print(f\"# Fold {fold}: mse: {mrrmse:5.3f} r2: {r2:5.3f} val='{val_cell_type}'\")\n                mrrmse_list.append(mrrmse)   \n\n                IX_stat += 1\n                # df_stat.loc[IX_stat,'Id'] = 'alpha'+str(alpha) + ' ' + 'smooth'+str(smoothing) + ' ' +'tsvd'+str(n_components)\n                df_stat.loc[IX_stat,'Id'] = 'best smooth and alpha'# 'sel_comp'+str(i_selected_component) + ' ' + 'sel_param'+str(smoothing[0]) + ' ' + 'smooth'+str(smoothing)\n                #df_stat.loc[IX_stat,'Id'] = 'sel_comp'+str(i_selected_component) + ' ' + 'sel_param'+str(smoothing[1])  + ' ' + 'smooth'+str(smoothing)\n                df_stat.loc[IX_stat,'Fold'] = fold\n                df_stat.loc[IX_stat,'Score'] = mrrmse\n                df_stat.loc[IX_stat,'r2'] = r2\n                df_stat.loc[IX_stat,'Fold inf'] = val_cell_type\n                df_stat.loc[IX_stat,'Time'] = np.round( time.time() - t0_fold, 1 )\n                \n#                 df_stat.loc[IX_stat,'i_selected_component'] = i_selected_component\n#                 df_stat.loc[IX_stat,'param_for_component'] = alpha_for_selected_component\n#                 df_stat.loc[IX_stat,'i_selected_component'] = i_selected_component\n#                 df_stat.loc[IX_stat,'param_for_component'] = smoothing[1]\n                \n                \n\n            mrrmse = np.array(mrrmse_list).mean()    \n            if verbose >= 10:\n                print(f\"# Overall {mrrmse:5.3f} \", df_stat.loc[IX_stat,'Id'] )    \n            numeric_columns = df_stat.select_dtypes(include='number').columns.tolist()\n            df_stat_aggr = df_stat.groupby('Id')[numeric_columns].mean()# .sort_values('Score') \n            if verbose >= 100:\n                display(  df_stat_aggr.tail(2)  )\n            df_stat_aggr = df_stat_aggr.sort_values('Score') \n            if verbose >= 100:\n                display(  df_stat_aggr.head(2)  )\n            df_stat_aggr.to_csv('df_stat_aggr.csv')\n            if (verbose > 0) and (cc%100 == 1):\n                print(cc, '%.1f secs passed'%(time.time() - t0))\n                display(  df_stat_aggr.head(5)  )\n            cc += 1\n                \n\ndisplay( df_stat )\ndf_stat.to_csv('df_stat.csv')\nnumeric_columns = df_stat.select_dtypes(include='number').columns.tolist()\ndf_stat_aggr = df_stat.groupby('Id')[numeric_columns].mean().sort_values('Score') \ndisplay(  df_stat_aggr  )\ndf_stat_aggr.to_csv('df_stat_aggr.csv')\n","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:03:50.299588Z","iopub.execute_input":"2023-10-12T11:03:50.299944Z","iopub.status.idle":"2023-10-12T11:04:00.816655Z","shell.execute_reply.started":"2023-10-12T11:03:50.299913Z","shell.execute_reply":"2023-10-12T11:04:00.816077Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_stat_aggr.head(10)","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:04:00.817581Z","iopub.execute_input":"2023-10-12T11:04:00.818224Z","iopub.status.idle":"2023-10-12T11:04:00.825576Z","shell.execute_reply.started":"2023-10-12T11:04:00.818201Z","shell.execute_reply":"2023-10-12T11:04:00.825069Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Show best_alpha_for_selected_component ","metadata":{}},{"cell_type":"code","source":"column_best_param = 'param_for_component'#  'smoothing' \n# column_best_param = 'alpha_for_selected_component' # For version 30-32- search of optimal alpha \nif column_best_param in df_stat_aggr.columns:\n    # For version 31,32 - search of optimal alpha \n    sr = pd.Series()\n    for i_selected_component in df_stat_aggr['i_selected_component'].unique():\n        m = df_stat_aggr['i_selected_component'] == i_selected_component\n        best_alpha_for_selected_component = df_stat_aggr[m].sort_values('Score',ascending = True)[column_best_param].iat[0]\n        sr.loc[i_selected_component] = best_alpha_for_selected_component\n\n    sr = sr.sort_index()  \n    sr.to_csv('best_param_for_selected_component.csv')\n    print(  sr  )\n    print()\n    print( dict(sr) )\n\n","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:04:00.826288Z","iopub.execute_input":"2023-10-12T11:04:00.826679Z","iopub.status.idle":"2023-10-12T11:04:00.838478Z","shell.execute_reply.started":"2023-10-12T11:04:00.826658Z","shell.execute_reply":"2023-10-12T11:04:00.837602Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Show top 50 results","metadata":{}},{"cell_type":"code","source":"df_stat_aggr.head(50)","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:04:00.839423Z","iopub.execute_input":"2023-10-12T11:04:00.839684Z","iopub.status.idle":"2023-10-12T11:04:00.861578Z","shell.execute_reply.started":"2023-10-12T11:04:00.839662Z","shell.execute_reply":"2023-10-12T11:04:00.860800Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## r2 and mmse are NOT consistent - strange\n\nnegative correlation and sorting by r2 and mmse gives different results ","metadata":{}},{"cell_type":"code","source":"display(df_stat_aggr.corr(method = 'spearman'))\ndf_stat_aggr.corr()\n","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:04:00.862800Z","iopub.execute_input":"2023-10-12T11:04:00.863071Z","iopub.status.idle":"2023-10-12T11:04:00.892582Z","shell.execute_reply.started":"2023-10-12T11:04:00.863050Z","shell.execute_reply":"2023-10-12T11:04:00.891477Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Sorting by r2 gives different results - strange ! ","metadata":{}},{"cell_type":"code","source":"df_stat.sort_values('r2', ascending = False)","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:04:00.894168Z","iopub.execute_input":"2023-10-12T11:04:00.894537Z","iopub.status.idle":"2023-10-12T11:04:00.906830Z","shell.execute_reply.started":"2023-10-12T11:04:00.894505Z","shell.execute_reply":"2023-10-12T11:04:00.905846Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Inference  Ridge+TE+tsvd - option - blend many folds\n\nWe use AmbrosM CV to choose the best model, we can retrain model on the whole train set with these params,\nbut seems better way to retrain model on subsets and blend the results - compare V1 and V2 here - LB improvement from 3 folds to 50 is 0.623->0.617.\nOn the other hand compare V33 V32 - almost no improvement from that trick 0.624 -> 0.625 \n\nThere is also :\n\nmode_submit = 'train_on_full'\n\nWhere we retrain model on the entire submit part ","metadata":{}},{"cell_type":"code","source":"X_submit_categorical = df_id_map[['cell_type','sm_name']]\ndisplay( X_submit_categorical.head(3) )\n\nY_submit = np.zeros( (255,18211 ) )\nY_submit.shape","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:04:00.908187Z","iopub.execute_input":"2023-10-12T11:04:00.908730Z","iopub.status.idle":"2023-10-12T11:04:00.924162Z","shell.execute_reply.started":"2023-10-12T11:04:00.908698Z","shell.execute_reply":"2023-10-12T11:04:00.923113Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \n\n# mode_submit = 'random_folds'# \nmode_submit = 'train_on_full'\nn_splits = 50 # if mode_submit == 'train_on_full' - that will NOT be used \n\nfrom sklearn.linear_model import Ridge\nfrom sklearn.model_selection import KFold\n\nfrom sklearn.metrics import r2_score\nfrom sklearn.linear_model import Ridge\nimport time\nfrom sklearn.metrics import r2_score\nfrom sklearn.decomposition import TruncatedSVD\n\n# from sklearn.preprocessing import TargetEncoder\nimport category_encoders as ce\n\n\nkf = KFold(n_splits=n_splits, random_state = 42, shuffle = True )\ndf_IX = pd.DataFrame(); df_IX['IX'] = range(len(df_de_train))\n\nverbose = 1\n# smoothing = (1e7,1e7) # smoothing = 100\n# # alpha = 1e6\n# n_components = 35 \n# dict_best_alpha_for_component = {0: 1000000.0, 1: 1000000.0, 2: 1000000.0, 3: 1000000.0, 4: 100000.0, 5: 1000000.0, 6: 1000000.0, 7: 1000000.0, 8: 100000.0, 9: 100000.0, 10: 10000.0, 11: 10000.0, 12: 100000.0, 13: 10000.0, 14: 100000.0, 15: 100000.0, 16: 100000.0, 17: 10000.0, 18: 10000.0, 19: 10000.0, 20: 10000.0, 21: 10000.0, 22: 10000.0, 23: 1000.0, 24: 10000.0}\n# # From versions 31,32 of the notebook:\n# # https://www.kaggle.com/code/alexandervc/op2-models-cv-tuning?scriptVersionId=144944894&cellId=19\n# n_components = 25\n# alpha_default = 1e6\n\ndf_stat_inference = pd.DataFrame(); IX_stat = -1\nY = df_de_train.iloc[:,5:].values\n\ncc = 0\nt0 = time.time()\nmrrmse_list = []\ni_blend = 0\nfold = 0\nfor fold, (train_index, test_index) in enumerate(kf.split(df_IX)):\n    print(f\"Fold {fold}:\")\n    #print(type(train_index), train_index.shape)\n    t0_fold = time.time()\n    reducer = TruncatedSVD(n_components=n_components, n_iter=7, random_state=42)\n\n    mask_va = df_IX['IX'].isin(test_index)# (df_de_train.cell_type == val_cell_type) & ~df_de_train.sm_name.isin(train_sm_names)\n    mask_tr = df_IX['IX'].isin(train_index)# ~mask_va # 485 or 487 training rows\n    if mode_submit == 'train_on_full':\n        mask_tr = pd.Series(index = df_IX.index, data = True)\n        mask_va = pd.Series(index = df_IX.index, data = True)\n\n    if verbose >= 1000:\n        print(mask_va.sum() , mask_tr.sum(), mask_va.sum()+ mask_tr.sum() , 'fold:', fold, val_cell_type )\n\n    Yr_train = reducer.fit_transform(Y[mask_tr,:])\n    if verbose >= 1000:\n        print('Yr_train.shape',Yr_train.shape)\n    X_train_categorical = df_de_train[['cell_type','sm_name']][mask_tr]\n\n    Yr_valid = reducer.transform(Y[mask_va,:])\n    X_valid_categorical = df_de_train[['cell_type','sm_name']][mask_va]\n\n#     enc_0 = ce.TargetEncoder(smoothing=smoothing[0])\n#     enc_1 = ce.TargetEncoder(smoothing=smoothing[1])\n    Yr_valid_pred = np.zeros( Yr_valid.shape )\n    Yr_submit_pred = np.zeros( (255, Yr_valid.shape[1] ) )\n    \n    for i_comp in range(Yr_train.shape[1]):\n\n#                     if i_comp == i_selected_component:\n#                         model = Ridge(alpha = alpha_for_selected_component )\n#                     else:\n#                         model = Ridge(alpha = alpha_default )\n\n        enc_0 = ce.TargetEncoder(smoothing=dict_best_smooth_celltype_for_component[i_comp]) # Best smooth0 found in v42\n        enc_1 = ce.TargetEncoder(smoothing=dict_best_smooth_compound_for_component[i_comp]) # Best smooth1 found in v42\n        model = Ridge(alpha = dict_best_alpha_for_component[i_comp] ) # Best alpha found in v 30-31\n\n        y_tr = Yr_train[:,i_comp]\n        X_train_enc = pd.concat( [enc_0.fit_transform(X_train_categorical[['cell_type']], y_tr),\n                                  enc_1.fit_transform(X_train_categorical[['sm_name']], y_tr) ] , axis = 1)\n        #print('X_train_enc.shape', X_train_enc.shape)\n        model.fit(X_train_enc, y_tr)\n        \n\n        #X_valid_enc = enc.transform(X_valid_categorical)\n        X_valid_enc = pd.concat( [enc_0.transform(X_valid_categorical[['cell_type']] ),\n                                  enc_1.transform(X_valid_categorical[['sm_name']] ) ] , axis = 1 )\n        y_pred = model.predict(X_valid_enc)\n        Yr_valid_pred[:,i_comp] = y_pred\n        \n        ###################################################### Submit part ############################################\n        \n        #X_submit_enc = enc.transform(X_submit_categorical)\n        X_submit_enc = pd.concat( [enc_0.transform(X_submit_categorical[['cell_type']] ),\n                                  enc_1.transform(X_submit_categorical[['sm_name']] ) ] , axis = 1 )\n        y_pred = model.predict(X_submit_enc)\n        Yr_submit_pred[:,i_comp] = y_pred\n        \n\n    Y_submit = (Y_submit*i_blend + reducer.inverse_transform( Yr_submit_pred  )  ) / (i_blend + 1) \n    \n    Y_valid = reducer.inverse_transform( Yr_valid_pred  )\n\n    mrrmse = np.sqrt(np.square(Y[mask_va,:] - Y_valid).mean(axis=1)).mean()\n    r2 = r2_score( Y[mask_va,:] ,Y_valid  )\n    if verbose >= 100:\n        print(f\"# Fold {fold}: mse: {mrrmse:5.3f} r2: {r2:5.3f} val='{val_cell_type}'\")\n    mrrmse_list.append(mrrmse)   \n\n    IX_stat += 1\n    # df_stat_inference.loc[IX_stat,'Id'] = 'alpha'+str(alpha) + ' ' + 'smooth'+str(smoothing) + ' ' +'tsvd'+str(n_components)\n    df_stat_inference.loc[IX_stat,'Id'] = 'best_alpha_for_each_component'#  'sel_comp'+str(i_selected_component) + ' ' + 'sel_alpha'+str(alpha_for_selected_component) + ' ' + 'smooth'+str(smoothing)\n    df_stat_inference.loc[IX_stat,'Fold'] = fold\n    df_stat_inference.loc[IX_stat,'Score'] = mrrmse\n    df_stat_inference.loc[IX_stat,'r2'] = r2\n    df_stat_inference.loc[IX_stat,'Time'] = np.round( time.time() - t0_fold, 1 )\n\n# #                 df_stat.loc[IX_stat,'i_selected_component'] = i_selected_component\n# #                 df_stat.loc[IX_stat,'alpha_for_selected_component'] = alpha_for_selected_component\n\n\n\n    mrrmse = np.array(mrrmse_list).mean()    \n    if verbose >= 10:\n        print(f\"# Overall {mrrmse:5.3f} \", df_stat_inference.loc[IX_stat,'Id'] )    \n    cc += 1\n    if mode_submit == 'train_on_full':\n        break\n    \n                \n\ndisplay( df_stat_inference )\ndf_stat_inference.to_csv('df_stat_inference.csv')\nnumeric_columns = df_stat_inference.select_dtypes(include='number').columns.tolist()\ndf_stat_aggr_inference = df_stat_inference.groupby('Id')[numeric_columns].mean().sort_values('Score') \ndisplay(  df_stat_aggr_inference  )\ndf_stat_aggr_inference.to_csv('df_stat_aggr_inference.csv')\n    ","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:04:00.929203Z","iopub.execute_input":"2023-10-12T11:04:00.929993Z","iopub.status.idle":"2023-10-12T11:04:04.191316Z","shell.execute_reply.started":"2023-10-12T11:04:00.929956Z","shell.execute_reply":"2023-10-12T11:04:04.190696Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## submit CSV prepare/save","metadata":{}},{"cell_type":"code","source":"%%time\nprint(Y_submit.shape)\nY_submit\n\ndf_submit = pd.DataFrame(Y_submit, columns = df_de_train.columns[5:])\ndf_submit.index.name = 'id'\nprint( df_submit.shape )\ndisplay(df_submit)\n\nif 1:\n    df_submit.to_csv('submission_RidgeTEtsvdLB0612.csv')","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:04:04.192124Z","iopub.execute_input":"2023-10-12T11:04:04.192749Z","iopub.status.idle":"2023-10-12T11:04:12.705503Z","shell.execute_reply.started":"2023-10-12T11:04:04.192724Z","shell.execute_reply":"2023-10-12T11:04:12.704375Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Ridge Many TE features for each target \n\nuse only compound - increases CV score (AmbrosM scheme)\n","metadata":{}},{"cell_type":"code","source":"%%time\nverbose = 1000\n\ndf_stat = pd.DataFrame(); IX_stat = -1\nY = df_de_train.iloc[:,5:].values\n\ncc = 0\nt0 = time.time()\nmrrmse_list = []; r2_list = []\nfor fold, val_cell_type in enumerate(validation_cell_types):\n    t0_fold = time.time()\n\n    mask_va = (df_de_train.cell_type == val_cell_type) & ~df_de_train.sm_name.isin(train_sm_names)\n    mask_tr = ~mask_va # 485 or 487 training rows\n    if verbose >= 1000:\n        print(mask_va.sum() , mask_tr.sum(), mask_va.sum()+ mask_tr.sum() , 'fold:', fold, val_cell_type )\n\n    Y_red_train = reducer.fit_transform(Y[mask_tr,:])\n    Y_red_valid = reducer.transform(Y[mask_va,:])\n    Y_red_valid_pred = np.zeros( Y_red_valid.shape )\n    if verbose >= 1000:\n        print('Yr_train.shape',Yr_train.shape)\n    X_train_categorical = df_de_train[['cell_type','sm_name']][mask_tr]\n    X_valid_categorical = df_de_train[['cell_type','sm_name']][mask_va]\n    X_train_enc = np.zeros( (X_train_categorical.shape[0], 2*Y_red_train.shape[1])  )\n    X_valid_enc = np.zeros( (X_valid_categorical.shape[0], 2*Y_red_train.shape[1])  )\n    for i_comp in range(Y_red_train.shape[1]):\n        enc_0 = ce.TargetEncoder(smoothing=dict_best_smooth_celltype_for_component[i_comp]) # Best smooth0 found in v42\n        enc_1 = ce.TargetEncoder(smoothing=dict_best_smooth_compound_for_component[i_comp]) # Best smooth1 found in v42\n        y_tr = Y_red_train[:,i_comp]\n#         X_train_enc[:,2*i_comp] = enc_0.fit_transform(X_train_categorical[['cell_type']], y_tr).values.ravel()\n#         X_valid_enc[:,2*i_comp] = enc_0.transform(X_valid_categorical[['cell_type']] ).values.ravel()\n        X_train_enc[:,2*i_comp+1] = enc_1.fit_transform(X_train_categorical[['sm_name']], y_tr).values.ravel()\n        X_valid_enc[:,2*i_comp+1] = enc_1.transform(X_valid_categorical[['sm_name']] ).values.ravel()\n        \n    vec_alpha = np.array(list(dict_best_alpha_for_component.values() ))\n    model = Ridge(alpha =  vec_alpha ) #\n    model.fit(X_train_enc, Y_red_train)\n\n    Y_red_valid_pred = model.predict(X_valid_enc)\n    Y_valid_pred = reducer.inverse_transform( Y_red_valid_pred  )\n\n    mrrmse = np.sqrt(np.square(Y[mask_va,:] - Y_valid_pred).mean(axis=1)).mean();     mrrmse_list.append(mrrmse)   \n    r2 = r2_score( Y[mask_va,:] ,Y_valid_pred  ); r2_list.append(r2)\n    if verbose >= 100:\n        print(f\"# Fold {fold}: mse: {mrrmse:5.3f} r2: {r2:5.3f} val='{val_cell_type}'\")\n\nmrrmse = np.array(mrrmse_list).mean();  r2 = np.array(r2_list).mean();                   \nif verbose >= 10:\n    print(f\"# Overall {mrrmse:5.3f},  r2 {r2:5.3f} \" )    \n","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:04:12.706762Z","iopub.execute_input":"2023-10-12T11:04:12.707079Z","iopub.status.idle":"2023-10-12T11:04:20.096438Z","shell.execute_reply.started":"2023-10-12T11:04:12.707019Z","shell.execute_reply":"2023-10-12T11:04:20.095481Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Retrain on full data and submit prepare ","metadata":{}},{"cell_type":"code","source":"%%time\nY_red_full = reducer.fit_transform(Y)\nif verbose >= 1000:\n    print('Y_red_full.shape',Y_red_full.shape)\nX_full_categorical = df_de_train[['cell_type','sm_name']]\nX_full_enc = np.zeros( (X_full_categorical.shape[0], 2*Y_red_full.shape[1])  )\nfor i_comp in range(Y_red_full.shape[1]):\n    enc_0 = ce.TargetEncoder(smoothing=dict_best_smooth_celltype_for_component[i_comp]) # Best smooth0 found in v42\n    enc_1 = ce.TargetEncoder(smoothing=dict_best_smooth_compound_for_component[i_comp]) # Best smooth1 found in v42\n    y_tr = Y_red_full[:,i_comp]\n#     X_full_enc[:,2*i_comp] = enc_0.fit_transform(X_full_categorical[['cell_type']], y_tr).values.ravel()\n    X_full_enc[:,2*i_comp+1] = enc_1.fit_transform(X_full_categorical[['sm_name']], y_tr).values.ravel()\n\nvec_alpha = np.array(list(dict_best_alpha_for_component.values() ))\nmodel = Ridge(alpha =  vec_alpha ) #\nmodel.fit(X_full_enc, Y_red_full)\n\n# Y_red_full_pred = np.zeros( Y_red_valid.shape )\nY_red_full_pred = model.predict(X_full_enc)\nY_full_pred = reducer.inverse_transform( Y_red_full_pred  )\n\nmrrmse = np.sqrt(np.square(Y - Y_full_pred).mean(axis=1)).mean();\nr2 = r2_score( Y ,Y_full_pred  ); \nprint(f\"# Fold {fold}: mse: {mrrmse:5.3f} r2: {r2:5.3f} Train Full Data\")\n\nX_submit_categorical = df_id_map[['cell_type','sm_name']]\ndisplay( X_submit_categorical.head(3) )\nX_submit_enc = np.zeros( (X_submit_categorical.shape[0], 2*Y_red_full.shape[1])  )\nfor i_comp in range(Y_red_full.shape[1]):\n    enc_0 = ce.TargetEncoder(smoothing=dict_best_smooth_celltype_for_component[i_comp]) # Best smooth0 found in v42\n    enc_1 = ce.TargetEncoder(smoothing=dict_best_smooth_compound_for_component[i_comp]) # Best smooth1 found in v42\n    y_tr = Y_red_full[:,i_comp]\n#     X_full_enc[:,2*i_comp] = enc_0.fit_transform(X_full_categorical[['cell_type']], y_tr).values.ravel()\n    X_full_enc[:,2*i_comp+1] = enc_1.fit_transform(X_full_categorical[['sm_name']], y_tr).values.ravel()\n#     X_submit_enc[:,2*i_comp] = enc_0.transform(X_submit_categorical[['cell_type']]).values.ravel()\n    X_submit_enc[:,2*i_comp+1] = enc_1.transform(X_submit_categorical[['sm_name']]).values.ravel()\n\n\n# Y_submit.shape\nY_red_submit_pred = model.predict(X_submit_enc)\nY_submit = reducer.inverse_transform( Y_red_submit_pred  )\n\n\ndf_submit = pd.DataFrame(Y_submit, columns = df_de_train.columns[5:])\ndf_submit.index.name = 'id'\nprint( df_submit.shape )\ndisplay(df_submit)\n\nif 1:\n    df_submit.to_csv('submission_RidgeManyTEtsvdCV0993.csv')","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:04:20.098097Z","iopub.execute_input":"2023-10-12T11:04:20.098838Z","iopub.status.idle":"2023-10-12T11:04:31.828477Z","shell.execute_reply.started":"2023-10-12T11:04:20.098803Z","shell.execute_reply":"2023-10-12T11:04:31.827827Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Advanced modeling. Many: models, blends, CV, features.","metadata":{}},{"cell_type":"markdown","source":"## CV-schemes ","metadata":{}},{"cell_type":"code","source":"%%time\ntrain_sm_names = ['Idelalisib', 'Crizotinib', 'Linagliptin', 'Palbociclib',\n       'Dabrafenib', 'Alvocidib', 'LDN 193189', 'R428',\n       'Porcn Inhibitor III', 'Belinostat', 'Foretinib', 'MLN 2238',\n       'Penfluridol', 'Dactolisib', 'O-Demethylated Adapalene',\n       'Oprozomib (ONX 0912)', 'CHIR-99021']\n# ['NK cells', 'T cells CD4+', 'T cells CD8+', 'T regulatory cells']\n# val_cell_type = 'NK cells'\n# mask_va = (df_de_train.cell_type == val_cell_type) & ~df_de_train.sm_name.isin(train_sm_names)\n# mask_tr = ~mask_va # 485 or 487 training rows\n# print('fold:',  mask_va.sum() , mask_tr.sum(), mask_va.sum()+ mask_tr.sum() ,  )\n\nfrom sklearn.metrics import r2_score\n\ndef get_list_fold_ids( CV_scheme, flag_add_full_data=False):\n    if CV_scheme ==   'AmbrosM':\n        list_fold_ids =  ['NK cells', 'T cells CD4+', 'T cells CD8+', 'T regulatory cells']\n    if flag_add_full_data:\n        list_fold_ids += ['full']\n    return list_fold_ids\n\ntrain_sm_names = ['Idelalisib', 'Crizotinib', 'Linagliptin', 'Palbociclib',\n       'Dabrafenib', 'Alvocidib', 'LDN 193189', 'R428',\n       'Porcn Inhibitor III', 'Belinostat', 'Foretinib', 'MLN 2238',\n       'Penfluridol', 'Dactolisib', 'O-Demethylated Adapalene',\n       'Oprozomib (ONX 0912)', 'CHIR-99021']\ndef get_fold_data( fold_id , CV_scheme='AmbrosM'):\n    if CV_scheme ==   'AmbrosM':\n        mask_va = (df_de_train.cell_type == fold_id) & ~df_de_train.sm_name.isin(train_sm_names)\n        mask_tr = ~mask_va # 485 or 487 training rows\n    return mask_tr, mask_va\n\ndef get_cv_scores( Y_true, Y_oof_pred , CV_scheme, verbose = 0 ):\n    list_fold_ids = get_list_fold_ids( CV_scheme, flag_add_full_data=False) \n    mrrmse_list = []; r2_list = []\n    for i_fold,fold_id in enumerate(list_fold_ids):\n        mask_tr,mask_va = get_fold_data(fold_id, CV_scheme = CV_scheme )\n        mrrmse = np.sqrt(np.square(Y_true[mask_va,:] - Y_oof_pred[mask_va,:]).mean(axis=1)).mean()\n        r2 = r2_score( Y_true[mask_va,:] , Y_oof_pred[mask_va,:] ) \n        mrrmse_list.append(mrrmse)   \n        r2_list.append(r2)   \n        if verbose >= 10:\n            print(f\"# Fold {fold_id}: mse: {mrrmse:5.3f} r2: {r2:5.3f} \")   \n\n    return mrrmse_list, r2_list\nmrrmse_list, r2_list =  get_cv_scores( Y, Y , CV_scheme = 'AmbrosM', verbose = 10 ) \nprint(mrrmse_list, r2_list   ) \nmask_tr,mask_va = get_fold_data('NK cells')\nprint(mask_tr.sum() , mask_va.sum() )\n","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:04:31.829754Z","iopub.execute_input":"2023-10-12T11:04:31.830229Z","iopub.status.idle":"2023-10-12T11:04:32.266681Z","shell.execute_reply.started":"2023-10-12T11:04:31.830205Z","shell.execute_reply":"2023-10-12T11:04:32.266059Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Features","metadata":{}},{"cell_type":"code","source":"from sklearn.preprocessing import OneHotEncoder\nX_full_categorical = df_de_train[['cell_type','sm_name']]\nX_full_categorical = df_de_train[['sm_name']]\n\nX_full_categorical.head(1)\nenc = OneHotEncoder(handle_unknown='ignore')\n\nX_full_one_hot = enc.fit_transform(X_full_categorical).toarray()\ntype(X_full_one_hot)","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:04:32.267688Z","iopub.execute_input":"2023-10-12T11:04:32.268075Z","iopub.status.idle":"2023-10-12T11:04:32.278525Z","shell.execute_reply.started":"2023-10-12T11:04:32.268051Z","shell.execute_reply":"2023-10-12T11:04:32.277817Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Models","metadata":{}},{"cell_type":"code","source":"%%time\nfrom sklearn.metrics import r2_score\nfrom sklearn.decomposition import TruncatedSVD\nfrom sklearn.linear_model import Ridge\nfrom sklearn.svm import LinearSVR\nimport lightgbm as lgb\n\nfrom sklearn.multioutput import MultiOutputRegressor\n\ndef get_model(main_config_model_feature_etc , verbose = 0):\n    \n    str_model_id =  str( main_config_model_feature_etc['model'] )\n    if main_config_model_feature_etc['model'] == 'Ridge':\n        alpha = main_config_model_feature_etc.get('alpha',1)\n        model = Ridge( alpha  )\n        str_model_id  = 'Ridge'+str(alpha)\n    elif main_config_model_feature_etc['model'] == 'LGB':\n        params_loc = main_config_model_feature_etc.get('LGB_params',{} )\n        model = lgb.LGBMRegressor( **params_loc ) \n    elif main_config_model_feature_etc['model'] == 'LGBcv1_036':\n        params_best1_cv1_036 = {'random_state': 0, 'n_estimators': 20, 'reg_alpha': 6.764079452929363, 'reg_lambda': 0.41900776876588564, \n        'colsample_bytree': 0.3, 'subsample': 0.7, 'max_depth': 1, 'learning_rate': 0.08456104070184789,  'num_leaves': 682, 'min_child_samples': 102}\n        model = lgb.LGBMRegressor( **params_best1_cv1_036 ) \n        \n    if verbose >= 100:\n        print( str_model_id )\n        print( model )\n        print( main_config_model_feature_etc )\n        \n    return model, str_model_id\n\nmodel, str_model_id = get_model({'model':'Ridge'}, verbose = 100)\nprint(model, str_model_id  ); print()\n\ndef get_reducer( main_config_model_feature_etc , verbose = 0): \n    if 'reducer' not in main_config_model_feature_etc.keys():\n        n_components = 25\n        reducer = TruncatedSVD(n_components=n_components, n_iter=7, random_state=42)\n        str_reducer_id = 'tsvd25'\n        return reducer, str_reducer_id \n    \n    str_reducer_id =  str( main_config_model_feature_etc['reducer'] )\n    if main_config_model_feature_etc['reducer']=='tsvd':\n        n_components = main_config_model_feature_etc.get('n_components',25)\n        reducer = TruncatedSVD(n_components=n_components, n_iter=7, random_state=42)\n        str_reducer_id = 'tsvd'+str( n_components )\n    \n    if verbose >= 100:\n        print( str_reducer_id )\n        print( reducer )\n        print( main_config_model_feature_etc )\n\n    return reducer, str_reducer_id\n\nget_reducer( {'reducer': 'tsvd'} , verbose = 100)","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:04:32.279525Z","iopub.execute_input":"2023-10-12T11:04:32.279744Z","iopub.status.idle":"2023-10-12T11:04:32.295368Z","shell.execute_reply.started":"2023-10-12T11:04:32.279725Z","shell.execute_reply":"2023-10-12T11:04:32.294508Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Modeling. Blend. ","metadata":{}},{"cell_type":"code","source":"%%time\nverbose = 1000\nimport time\n\n\nY = df_de_train.iloc[:,5:].values\nX = X_full_one_hot\n\nY_oof_pred_blend = np.zeros_like( Y ); i_oof_blend = 0 \n# print(); print('Start training models',datetime.datetime.now()) ; print()\n########################## Main modelling  ###########################################\nmain_config_model_feature_etc1 = {'model':'Ridge', 'alpha':1e3} \nmain_config_model_feature_etc2 = {'model':'Ridge', 'alpha':1e2, 'reducer': 'tsvd', 'n_components':10} \nlist_main_config_model_feature_etc = [main_config_model_feature_etc1, main_config_model_feature_etc2 ]\nlist_main_config_model_feature_etc = []\nfor n_comp in range(25,35):\n    main_config_model_feature_etc = {'model':'Ridge', 'alpha':1e2, 'reducer': 'tsvd', 'n_components':n_comp} \n    list_main_config_model_feature_etc += [main_config_model_feature_etc ]\n\nfor main_config_model_feature_etc  in list_main_config_model_feature_etc:\n    Y_oof_pred = np.zeros_like( Y )\n    CV_scheme = 'AmbrosM'\n    list_fold_ids = get_list_fold_ids( CV_scheme, flag_add_full_data=False)\n    reducer, str_reducer_id = get_reducer( main_config_model_feature_etc , verbose = 0) \n\n    t0_all_folds = time.time()\n    for i_fold,fold_id in enumerate(list_fold_ids):\n\n        mask_tr,mask_va = get_fold_data(fold_id, CV_scheme = CV_scheme )\n\n        Y_train = Y[mask_tr,:]\n        Y_valid = Y[mask_va,:]\n        Y_proc_train = reducer.fit_transform(Y_train)\n        Y_proc_valid = reducer.transform(Y_valid)\n\n\n        Y_proc_valid_pred = np.zeros( Y_proc_valid.shape )\n        Y_proc_train_pred = np.zeros( Y_proc_train.shape )\n        for i_target in range(Y_proc_train.shape[1] ):\n\n            X_train = X[mask_tr,:]\n            X_valid = X[mask_va,:]\n\n            model, str_model_id = get_model(main_config_model_feature_etc)\n\n            model.fit(X_train, Y_proc_train[:,i_target] )\n            Y_proc_valid_pred[:, i_target ] =  model.predict(X_valid )\n            Y_proc_train_pred[:, i_target ] =  model.predict(X_train )\n\n        Y_valid_pred = reducer.inverse_transform( Y_proc_valid_pred )\n        Y_train_pred = reducer.inverse_transform( Y_proc_train_pred )\n\n        Y_oof_pred[mask_va,:] = Y_valid_pred\n\n        mrrmse = np.sqrt(np.square(Y_valid - Y_valid_pred).mean(axis=1)).mean()\n        r2 = r2_score( Y_valid , Y_valid_pred ) \n        if verbose >= 100:\n            print(f\"# Fold {fold_id}: mse: {mrrmse:5.3f} r2: {r2:5.3f} Valid\")   \n\n        mrrmse = np.sqrt(np.square(Y_train - Y_train_pred).mean(axis=1)).mean()\n        r2 = r2_score( Y_train , Y_train_pred ) \n        if verbose >= 1000:\n            print(f\"# Fold {fold_id}: mse: {mrrmse:5.3f} r2: {r2:5.3f} Train\")        \n\n\n    mrrmse_list, r2_list =  get_cv_scores( Y, Y_oof_pred , CV_scheme = CV_scheme, verbose = 0 )\n    if verbose >= 10:\n        print(mrrmse_list, r2_list )\n        print(np.mean( mrrmse_list ) , np.mean(r2_list ))\n\n    Y_oof_pred_blend =  ( Y_oof_pred_blend * i_oof_blend + Y_oof_pred)/(i_oof_blend + 1)\n    i_oof_blend += 1 \n    mrrmse_list, r2_list =  get_cv_scores( Y, Y_oof_pred_blend , CV_scheme = CV_scheme, verbose = 0 )\n    if verbose >= 10:\n        print(mrrmse_list, r2_list )\n        print(np.mean( mrrmse_list ) , np.mean(r2_list ))\n\n    if verbose >= 10:\n        print(f'{time.time()-t0_all_folds:.1f} secs passed')\n        ","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:04:32.296488Z","iopub.execute_input":"2023-10-12T11:04:32.296719Z","iopub.status.idle":"2023-10-12T11:06:02.713190Z","shell.execute_reply.started":"2023-10-12T11:04:32.296701Z","shell.execute_reply":"2023-10-12T11:06:02.712200Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Blending with aggregates by solely compound and solely cell-type (following ZXMKCD and LIUDA CHELDIEVA )\n\nIt was observed by ZXMKCD https://www.kaggle.com/code/zmcxjt/streamlined-baseline-approach that blending with simple aggregate by compound improves the score. \n\nAnd then by LIUDA CHELDIEVA https://www.kaggle.com/code/liudacheldieva/streamlined-baseline-approach that adding aggregate by cell-type improves it more. ","metadata":{}},{"cell_type":"markdown","source":"## Prepare simple predictions by aggergating by compounds and cell types separately","metadata":{}},{"cell_type":"code","source":"%%time \n#group by drug and take the mean\ndf_tmp = df_de_train.iloc[:, [1] + list(range(5, df_de_train.shape[1]))] # Take only numeric columns and \"sm_name\"\ndf_aggr = df_tmp.groupby('sm_name').mean().reset_index()\nprint(df_aggr.shape)\ndf_submit_aggr_compound = pd.merge( df_id_map,  df_aggr, on='sm_name', how = 'left' ).sort_values('id').drop(columns = ['cell_type', 'sm_name']).set_index('id')\nprint(df_submit_aggr_compound.shape)\ndisplay(df_submit_aggr_compound.head(3))\n\nprint( )\n\ndf_tmp = df_de_train.iloc[:, [0] + list(range(5, df_de_train.shape[1]))] # Take only numeric columns and \"cell_type\"\ndf_aggr = df_tmp.groupby('cell_type').mean().reset_index()\nprint(df_aggr.shape)\ndf_submit_aggr_cell_type = pd.merge( df_id_map,  df_aggr, on='cell_type', how = 'left' ).sort_values('id').drop(columns = ['cell_type', 'sm_name']).set_index('id')\nprint(df_submit_aggr_cell_type.shape)\ndisplay(df_submit_aggr_cell_type.head(3))\n","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:06:02.714817Z","iopub.execute_input":"2023-10-12T11:06:02.715207Z","iopub.status.idle":"2023-10-12T11:06:03.310285Z","shell.execute_reply.started":"2023-10-12T11:06:02.715173Z","shell.execute_reply":"2023-10-12T11:06:03.309627Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Prepare ensemble","metadata":{}},{"cell_type":"code","source":"%%time\nfn = 'submission_RidgeTEtsvdLB0612.csv'\ndf_submit_model1 = pd.read_csv(fn, index_col = 0)\nfn = 'submission_RidgeManyTEtsvdCV0993.csv'\ndf_submit_model2 = pd.read_csv(fn, index_col = 0)\n\n# LB0.605\ndf_submit = 0.45*(df_submit_model1) + 0.45 * df_submit_aggr_compound + 0.1 * df_submit_aggr_cell_type\ndf_submit.to_csv('submission_ensemble1.csv')\n# LB0.605 also \ndf_submit = 0.45*(0.5*(df_submit_model1+df_submit_model2)) + 0.45 * df_submit_aggr_compound + 0.1 * df_submit_aggr_cell_type\ndf_submit.to_csv('submission_ensemble2.csv')\ndisplay(df_submit)","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:06:03.311519Z","iopub.execute_input":"2023-10-12T11:06:03.312233Z","iopub.status.idle":"2023-10-12T11:06:25.833177Z","shell.execute_reply.started":"2023-10-12T11:06:03.312197Z","shell.execute_reply":"2023-10-12T11:06:25.832292Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Blend some precomputed solutions\n\n    # LB 0.604 blend all 9 submits with priors \n    # LB 0.604 (with priors, but selected 4 ): fn = 'submission_ensemble_'+str(i0)+'_0604-0604selected_'+ subdirectory_path_postfix+'.csv'\n    # LB 0.610 (no priors, 4 selected scores 0.612-0.615 ): fn = 'submission_simple_ensemble_'+str(i0)+'_0612-0615selected_'+ subdirectory_path_postfix+'.csv'\n","metadata":{}},{"cell_type":"code","source":"%%time\nfn = 'submission_Random_20_42_blend_with_priors_0.45_0.1.csv'\nfn = 'submission_simple.csv'\n\nimport os\ni0 = 0\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        #if ('0605' in dirname) or ( '0604' in dirname ):\n        if ('0612' in dirname) or ( '0613' in dirname ) or ( '0615' in dirname ):\n            if filename == fn:\n                full_fn = os.path.join(dirname, filename)\n                print(full_fn)\n                if i0 == 0:\n                    df_submit = pd.read_csv(full_fn, index_col = 0)\n                else:\n                    df_submit += pd.read_csv(full_fn, index_col = 0)\n                i0 += 1\n            \nprint(i0)            \ndf_submit /= i0\ndf_submit","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:32:45.006934Z","iopub.execute_input":"2023-10-12T11:32:45.007309Z","iopub.status.idle":"2023-10-12T11:33:00.034573Z","shell.execute_reply.started":"2023-10-12T11:32:45.007282Z","shell.execute_reply":"2023-10-12T11:33:00.033577Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nimport datetime\ncurrent_datetime = datetime.datetime.now()\nsubdirectory_path_postfix =  str(current_datetime)[5:16].replace(' ','-').replace(':','-')\n    \n# LB 0.604 (with priors): \n#fn = 'submission_ensemble_'+str(i0)+'_0604-0604selected_'+ subdirectory_path_postfix+'.csv'\n# LB 0.610 (no priors ): \nfn = 'submission_simple_ensemble_'+str(i0)+'_0612-0615selected_'+ subdirectory_path_postfix+'.csv'\nprint(fn)\ndf_submit.to_csv( fn)\n","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:39:52.184124Z","iopub.execute_input":"2023-10-12T11:39:52.184858Z","iopub.status.idle":"2023-10-12T11:40:00.500307Z","shell.execute_reply.started":"2023-10-12T11:39:52.184826Z","shell.execute_reply":"2023-10-12T11:40:00.499478Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# LB0.604 (V81) = LB 0.610 (V80) +  0.45 * df_submit_aggr_compound + 0.1 * df_submit_aggr_cell_type\ndf_submit = 0.45*(df_submit) + 0.45 * df_submit_aggr_compound + 0.1 * df_submit_aggr_cell_type\ndf_submit.to_csv('submission_first_blend_4_selected_models_get0610_then_add_priors.csv')\n","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:41:17.377731Z","iopub.execute_input":"2023-10-12T11:41:17.378880Z","iopub.status.idle":"2023-10-12T11:41:25.421930Z","shell.execute_reply.started":"2023-10-12T11:41:17.378837Z","shell.execute_reply":"2023-10-12T11:41:25.421194Z"},"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":"markdown","source":"# Final timing","metadata":{}},{"cell_type":"code","source":"print('%.1f seconds passed total '%(time.time()-t0start) )\nprint('%.1f minutes passed total '%( (time.time()-t0start)/60)  )\nprint('%.2f hours passed total '%( (time.time()-t0start)/3600)  )","metadata":{"execution":{"iopub.status.busy":"2023-10-12T11:06:25.834281Z","iopub.execute_input":"2023-10-12T11:06:25.834520Z","iopub.status.idle":"2023-10-12T11:06:25.839442Z","shell.execute_reply.started":"2023-10-12T11:06:25.834501Z","shell.execute_reply":"2023-10-12T11:06:25.838584Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}