{"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\nCreate models on different subsamples and compare features importances.\n\nLasso models\n\n#### Findings:\n\n**alpha = 0.1** for CD36 does not depend on subsample size ( train_size - parameter ) \n\n**Score: 0.79-0.8** correlation (Pearson) with \"y_true\" (CD36) for all subsample sizes \n\n**Lasso seems not always works fine**  For target CD44 alpha = 0.01 is optimal (from the list), but the number of non-zero features is almost 10_000 - that seems quite inappropriately large. Enforcing alpha = 0.1 a bit decreases the quality non-zero features becomes - about 400. Intersection of 100 trials - 147. \n\n\n#### Versions: \n\n    1,2 - n_trials = 10\n    3 - n_trials = 100 \n    4 - draft save  \n    5 - n_trials = 10, train_size = 0.75\n    6 - n_trials = 100, train_size = 0.75\n    7 - n_trials = 20, train_size = 0.75    \n    8 - n_trials = 100, train_size = 0.5 - repeat   \n    9 - n_trials = 10, train_size = 0.5 - repeat   \n    10,11 - n_trials = 10 , train_size = 0.1 \n    12 - n_trials = 100 , train_size = 0.1 \n    13 - n_trials = 100 , train_size = 0.9\n    14,15,16,17 - n_trials = 10 , train_size = 0.9\n    17,18 - crash\n    19 - n_trials = 100 , train_size = 0.9 - 9 hours \n    20,21,22  - n_trials = 10,50,100, train_size = 0.9 - repeat\n    23,24,25 - crash\n    26,27,28  - n_trials = 10 train_size = 0.95 \n    29  quick save - 2 trials, fast mode , train_size = 0.5 - 7 minutes \n    30  - n_trials = 100 train_size = 0.95  # Cancelled after 12 hours\n    31  - n_trials = 100 train_size = 0.9   # Ran in 9 hours and 28 minutes\n    32  - n_trials = 100 train_size = 0.75  # Ran in 8 hours and 34 minutes\n    33  - n_trials = 100 train_size = 0.5   # Ran in 3 hours and 06 minutes\n    34  - n_trials = 100 train_size = 0.1   # Ran in 3 hours and 49 minutes\n    \n    35,36  - n_trials = 10 train_size = 0.1 - 30.4 minutes , 860 non-zero features - for 1 trial, but intersection between two - 126 and later falls down to 18: [823, 126, 59, 30, 21, 20, 19, 19, 19, 18]\n    37 - n_trials = 10 train_size = 0.5\n    38 - n_trials = 10 train_size = 0.75\n    39 - n_trials = 10 train_size = 0.9\n    40 - n_trials = 10 train_size = 0.95 - 120 minutes, 0.802 0.807 158 non-zero feautures  Intersection: [158, 150, 144, 138, 135, 133, 133, 131, 129, 129], Correlations - 0.99 \n    \n    \n    ====================================================================\n    CD44\n    ====================================================================\n    41 n_trials = 10 train_size = 0.1 # Ran in 36 minutes and 27 seconds # 0.675972\t0.744553\t580.500000  alpha:  0.1\n    42 n_trials = 10 train_size = 0.5 # Ran in 34 minutes and 42 seconds # 0.683030\t0.693261\t409.700000  alpha = 0.1\n\n    43 n_trials = 100 train_size = 0.1 # Ran in 2 hours and 56 minutes\n    44 n_trials = 100 train_size = 0.5 # Ran in 6 hours and 22 minutes\n    45 n_trials = 100 train_size = 0.75 # Cancelled after 12 hours\n    46 n_trials = 100 train_size = 0.9 # Cancelled after 12 hours\n    47 n_trials = 100 train_size = 0.95 # Cancelled after 12 hours\n\n    48 n_trials = 10 train_size = 0.75 # Ran in 6 hours and 14 minutes #  0.693954\t0.804693\t9104.300000 alpha: 0.01 (! different alpha ! , standard alpha = 0.1 gives a bit less score: 0.679396 with 391 features )  [9151, 5776, 4465, 3709, 3269, 2928, 2684, 2459, 2286, 2152]\n    49 n_trials = 10 train_size = 0.9  # Ran in 5 hours and 10 minutes\n    50 n_trials = 10 train_size = 0.95 # Run in 4.75 hours # 0.703858\t0.786153\t8613.700000 # alpha: 0.01 \n                All: [8555, 7174, 6502, 6102, 5821, 5619, 5438, 5308, 5176, 5074]\n                Count Non-zero median importances  8762\n                Average count genes non-zero for trials:  8613.7\n                Count genes non-zero for all trials:  5074\n                number of non-zero coeffients common for Pairs of trials:   7330.29\n                number of non-zero coeffients common for TRIPLES of trials:   6531.725\n\n    \n    ====================================================================\n    CD36, CD44, \n    ====================================================================\n    51  - CD36, n_trials = 100 train_size = 0.1 # 3.21 hours # 0.790715\t0.867522\t864.480000 # \n                Count Non-zero median importances  84\n                Average count genes non-zero for trials:  864.48\n                First 30: [897, 119, 49, 30, 27, 26, 22, 20, 18, 16, 15, 15, 14, 13, 11, 9, 9, 9, 9, 8, 8, 8, 7, 7, 7, 7, 7, 7, 7, 7]\n                Last 30:  [5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5]\n                \n    52  - CD44  n_trials = 20 train_size = 0.9  # Expexted is about 10 hours - compare to V49 # 7.52 hours\n                # 0.700470\t0.790130\t8719.250000\n                Count Non-zero median importances  8514\n                Average count genes non-zero for trials:  8719.25\n                Count genes non-zero for all trials:  3246\n                All: [8717, 6662, 5787, 5242, 4876, 4607, 4399, 4229, 4091, 3974, 3870, 3770, 3691, 3601, 3511, 3442, 3384, 3332, 3285, 3246]\n                \n    53  - CD44  n_trials = 15 train_size = 0.75  # Expexted is about 10 hours - compare to V48 # 7.58 hours\n        All: [9018, 5717, 4382, 3657, 3160, 2854, 2619, 2407, 2251, 2120, 2033, 1938, 1858, 1770, 1714]\n        Count Non-zero median importances  7920\n        Average count genes non-zero for trials:  9110\n        Count genes non-zero for all trials:  1714\n\n         Alpha           Train      Test         N_nonzero_feature\n        1.000000e-03\t0.849821\t0.611386\t19942.0\n        1.000000e-02\t0.803358\t0.698067\t9261.0\n        1.000000e-01\t0.688624\t0.683882\t407.0\n        1.000000e+00\t0.498581\t0.499372\t6.0    \n\n    ====================================================================\n    CD44, enforce alpha = 0.1 \n    ====================================================================\n    \n    54 - CD44  enforce alpha = 0.1 , n_trials = 100 train_size = 0.5 # 6.17 hours # \t0.682333\t0.693402\t412.080000\n            Count Non-zero median importances  306\n            Average count genes non-zero for trials:  412.08\n            Count genes non-zero for all trials:  69\n            [390, 251, 202, 174, 163, 154, 148, 139, 131, 129, 127, 124, 122, 121, 118, 117, 115, 111, 107, 105, 104, 104, 102, 100, 98, 96, 96, 96, 94, 94, 92, 89, 88, 87, 86, 86, 86, 85, 84, 84, 84, 84, 84, 84, 84, 84, 84, 83, 82, 82, 82, 81, 81, 81, 80, 79, 79, 79, 78, 78, 78, 78, 77, 77, 77, 77, 76, 76, 76, 76, 76, 76, 76, 75, 75, 75, 75, 75, 74, 74, 73, 72, 72, 72, 72, 72, 72, 72, 71, 71, 71, 71, 71, 71, 71, 70, 70, 70, 70, 69]\n\n    55 - CD44  enforce alpha = 0.1 , n_trials = 100 train_size = 0.75  # 7.28 hours # 0.682817\t0.689567\t400.210000\n        Count Non-zero median importances  363\n        Average count genes non-zero for trials:  400.21\n        Count genes non-zero for all trials:  147\n            [406, 305, 271, 250, 233, 221, 217, 212, 203, 198, 193, 190, 190, 187, 185, 184, 183, 181, 180, 178, 177, 176, 176, 176, 175, 173, 172, 171, 171, 171, 171, 170, 170, 168, 166, 165, 165, 164, 164, 163, 163, 163, 163, 163, 161, 161, 161, 160, 160, 160, 160, 160, 160, 159, 158, 156, 156, 156, 156, 156, 155, 155, 154, 153, 153, 153, 153, 152, 152, 152, 152, 151, 151, 151, 151, 151, 150, 150, 150, 149, 149, 148, 148, 148, 148, 147, 147, 147, 147, 147, 147, 147, 147, 147, 147, 147, 147, 147, 147, 147]\n    \n    56 - CD44  enforce alpha = 0.1 , n_trials = 100 train_size = 0.9 # 8 hours, 24 minutes  \n        # 0.681424\t0.688602\t395.720000\n        Count Non-zero median importances  378\n        Average count genes non-zero for trials:  395.72\n        Count genes non-zero for all trials:  216\n        All: [388, 328, 315, 304, 291, 286, 284, 275, 271, 267, 265, 264, 263, 261, 259, 258, 256, 250, 250, 250, 249, 248, 248, 247, 247, 245, 243, 242, 242, 240, 240, 240, 238, 238, 238, 237, 237, 236, 235, 235, 235, 235, 234, 234, 232, 232, 232, 232, 232, 232, 232, 230, 230, 230, 230, 229, 228, 228, 227, 226, 226, 226, 226, 226, 226, 225, 225, 225, 225, 225, 225, 225, 224, 224, 224, 223, 223, 223, 222, 222, 221, 221, 221, 221, 221, 221, 221, 221, 221, 220, 220, 218, 217, 217, 217, 217, 217, 216, 216, 216]\n        \n    57  - CD44  enforce alpha = 0.1 , n_trials = 100 train_size = 0.1  \n        \n    ====================================================================\n    CD88 \n    ====================================================================\n        \n    58  - CD88  enforce alpha = 0.1 , n_trials = 100 train_size = 0.1  \n    59  - CD88  enforce alpha = 0.1 , n_trials = 100 train_size = 0.9 \n    \n    60  - CD88  NO enforce alpha , n_trials = 100 train_size = 0.1  \n    61  - CD88  NO enforce alpha , n_trials = 100 train_size = 0.5  \n    ","metadata":{}},{"cell_type":"markdown","source":"# Key Params","metadata":{}},{"cell_type":"code","source":"target_name = 'CD88'\n\ntrain_size = 0.5\n\nn_trials = 100\n\nstr_method = 'Lasso '\n\nfast_mode = 0\n\n\nimport time\nt0start = time.time() ","metadata":{"execution":{"iopub.status.busy":"2023-01-10T20:33:14.983616Z","iopub.execute_input":"2023-01-10T20:33:14.984270Z","iopub.status.idle":"2023-01-10T20:33:15.014001Z","shell.execute_reply.started":"2023-01-10T20:33:14.984137Z","shell.execute_reply":"2023-01-10T20:33:15.012913Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Results ","metadata":{}},{"cell_type":"code","source":"# From the previous notebook with Ridge model: \nimport numpy as np\nimport pandas as pd \n\nprint('Unfortunately intersection of lists from 10 trials and 100 both 47 gives 37, so not fully coincident.  ')\n\n# From Version 3: \nl_10trials_top100_gives47intersection = ['ENSG00000135218_CD36', 'ENSG00000229988_HBBP1', 'ENSG00000112077_RHAG', 'ENSG00000137801_THBS1', 'ENSG00000130303_BST2', 'ENSG00000109099_PMP22', 'ENSG00000110092_CCND1', 'ENSG00000029534_ANK1', 'ENSG00000223609_HBD', 'ENSG00000129824_RPS4Y1', 'ENSG00000204103_MAFB', 'ENSG00000101162_TUBB1', 'ENSG00000235169_SMIM1', 'ENSG00000267279_AC090409.1', 'ENSG00000164946_FREM1', 'ENSG00000168685_IL7R', 'ENSG00000244734_HBB', 'ENSG00000073464_CLCN4', 'ENSG00000197993_KEL', 'ENSG00000113924_HGD', 'ENSG00000185198_PRSS57', 'ENSG00000139174_PRICKLE1', 'ENSG00000166091_CMTM5', 'ENSG00000047648_ARHGAP6', 'ENSG00000196565_HBG2', 'ENSG00000065534_MYLK', 'ENSG00000165682_CLEC1B', 'ENSG00000107984_DKK1', 'ENSG00000088053_GP6', 'ENSG00000206172_HBA1', 'ENSG00000250361_GYPB', 'ENSG00000170873_MTSS1', 'ENSG00000187609_EXD3', 'ENSG00000107130_NCS1', 'ENSG00000115461_IGFBP5', 'ENSG00000169704_GP9', 'ENSG00000103522_IL21R', 'ENSG00000138135_CH25H', 'ENSG00000072274_TFRC', 'ENSG00000169403_PTAFR', 'ENSG00000047597_XK', 'ENSG00000073737_DHRS9', 'ENSG00000077984_CST7', 'ENSG00000251002_AC244502.1', 'ENSG00000092621_PHGDH', 'ENSG00000169071_ROR2', 'ENSG00000213719_CLIC1']\n# From Version 2: \nl_100trials_top300_gives47intersection = ['ENSG00000135218_CD36', 'ENSG00000229988_HBBP1', 'ENSG00000130303_BST2', 'ENSG00000112077_RHAG', 'ENSG00000109099_PMP22', 'ENSG00000137801_THBS1', 'ENSG00000110092_CCND1', 'ENSG00000029534_ANK1', 'ENSG00000223609_HBD', 'ENSG00000164946_FREM1', 'ENSG00000073464_CLCN4', 'ENSG00000235169_SMIM1', 'ENSG00000267279_AC090409.1', 'ENSG00000101162_TUBB1', 'ENSG00000129824_RPS4Y1', 'ENSG00000168685_IL7R', 'ENSG00000026025_VIM', 'ENSG00000197993_KEL', 'ENSG00000244734_HBB', 'ENSG00000113924_HGD', 'ENSG00000166091_CMTM5', 'ENSG00000047648_ARHGAP6', 'ENSG00000139174_PRICKLE1', 'ENSG00000185198_PRSS57', 'ENSG00000065534_MYLK', 'ENSG00000115461_IGFBP5', 'ENSG00000169704_GP9', 'ENSG00000088053_GP6', 'ENSG00000187609_EXD3', 'ENSG00000175084_DES', 'ENSG00000133742_CA1', 'ENSG00000107130_NCS1', 'ENSG00000072274_TFRC', 'ENSG00000047597_XK', 'ENSG00000131981_LGALS3', 'ENSG00000205542_TMSB4X', 'ENSG00000170873_MTSS1', 'ENSG00000198673_FAM19A2', 'ENSG00000169403_PTAFR', 'ENSG00000198400_NTRK1', 'ENSG00000251002_AC244502.1', 'ENSG00000092621_PHGDH', 'ENSG00000213719_CLIC1', 'ENSG00000169071_ROR2', 'ENSG00000251562_MALAT1', 'ENSG00000091409_ITGA6', 'ENSG00000174788_PCP2']\nm = pd.Series(index = l_10trials_top100_gives47intersection, dtype = float).index.isin(l_100trials_top300_gives47intersection )\ns = list( pd.Series(index = l_10trials_top100_gives47intersection, dtype = float)[m].index )\nprint(s[:2])\nprint( len(s), len(l_10trials_top100_gives47intersection), len(l_100trials_top300_gives47intersection ) )\nl = [t.split('_')[1] for t in s ]\nprint('Intersection:', len(l), l)\n\n","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:42:58.749897Z","iopub.execute_input":"2023-01-10T16:42:58.750501Z","iopub.status.idle":"2023-01-10T16:42:58.775382Z","shell.execute_reply.started":"2023-01-10T16:42:58.750443Z","shell.execute_reply":"2023-01-10T16:42:58.773814Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Preparations","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\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-01-10T16:42:58.778734Z","iopub.execute_input":"2023-01-10T16:42:58.779780Z","iopub.status.idle":"2023-01-10T16:42:58.801786Z","shell.execute_reply.started":"2023-01-10T16:42:58.779724Z","shell.execute_reply":"2023-01-10T16:42:58.800737Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport seaborn as sns","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:42:58.803809Z","iopub.execute_input":"2023-01-10T16:42:58.804430Z","iopub.status.idle":"2023-01-10T16:42:58.812831Z","shell.execute_reply.started":"2023-01-10T16:42:58.804372Z","shell.execute_reply":"2023-01-10T16:42:58.811650Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load Data","metadata":{}},{"cell_type":"code","source":"%%time\nfilename_rna_data = '/kaggle/input/open-problems-multimodal/train_cite_inputs.h5'\ndf_rna = pd.read_hdf(filename_rna_data)\ndisplay(df_rna) \n\n#%%time\ndf_y = pd.read_hdf('/kaggle/input/open-problems-multimodal/train_cite_targets.h5')\ndisplay(df_y)\n\nfn = '/kaggle/input/open-problems-multimodal/metadata.csv'\ndf_meta = pd.read_csv(fn, index_col = 0 )\ndf_meta\n# Cut only train cite-seq part: \nd = pd.DataFrame(index = df_y.index)\nprint(d.shape)\ndf_meta = d.join(df_meta, how = 'left')\ndf_meta","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:42:58.815734Z","iopub.execute_input":"2023-01-10T16:42:58.816896Z","iopub.status.idle":"2023-01-10T16:43:59.941215Z","shell.execute_reply.started":"2023-01-10T16:42:58.816819Z","shell.execute_reply":"2023-01-10T16:43:59.939783Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Modeling ","metadata":{}},{"cell_type":"code","source":"y = df_y[target_name].values\nprint(y.shape, type(y) )\n\nX = ( df_rna.values )\nprint(X.shape)\n","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:43:59.943071Z","iopub.execute_input":"2023-01-10T16:43:59.943828Z","iopub.status.idle":"2023-01-10T16:43:59.967780Z","shell.execute_reply.started":"2023-01-10T16:43:59.943781Z","shell.execute_reply":"2023-01-10T16:43:59.966443Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom sklearn.preprocessing import StandardScaler\nscaler = StandardScaler()\n\nX = scaler.fit_transform( X )\nprint(X.shape)\n","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:43:59.969513Z","iopub.execute_input":"2023-01-10T16:43:59.969965Z","iopub.status.idle":"2023-01-10T16:44:24.902701Z","shell.execute_reply.started":"2023-01-10T16:43:59.969922Z","shell.execute_reply":"2023-01-10T16:44:24.901509Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"##  Determine optimal Alpha","metadata":{}},{"cell_type":"code","source":"from sklearn.linear_model import LassoCV\nfrom sklearn.linear_model import Lasso\nfrom sklearn.linear_model import RidgeCV\nfrom sklearn.linear_model import Ridge\n\nfrom sklearn.model_selection import cross_val_predict\nfrom sklearn.model_selection import KFold \nfrom sklearn.metrics import r2_score\nfrom sklearn.metrics import mean_squared_error","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:44:24.904625Z","iopub.execute_input":"2023-01-10T16:44:24.905104Z","iopub.status.idle":"2023-01-10T16:44:24.913577Z","shell.execute_reply.started":"2023-01-10T16:44:24.905058Z","shell.execute_reply":"2023-01-10T16:44:24.910916Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\n\nimport time \nprint('X.shape, y.shape:', X.shape, y.shape, 'train_size:', train_size)\n\nif fast_mode and (target_name == 'CD36' ) and (filename_rna_data == '/kaggle/input/open-problems-multimodal/train_cite_inputs.h5') and ( X.shape[1] == 22050) and ( 'Lasso' in str_method ) : \n     alpha_selected =  0.1 # We have already found that 0.1 is optimal alpha for CD36 full rna-data data for KaggleNIPS22 dataset, so we will just use it \nelse:\n    p = np.random.permutation(len(y))\n    N = int(len(y) *  train_size )\n    IX_train = np.arange(len(y))[p][:N]\n    IX_test =  np.arange(len(y))[p][N:]\n    print(\" len(IX_train), len(IX_test):\",  len(IX_train), len(IX_test) )\n\n    df_models_1 = pd.DataFrame()\n\n\n    for i,alpha in enumerate( [1e-3, 1e-2, 1e-1, 1,1e1,1e2 ,1e3,1e4,1e5,1e6,1e7] ):\n        t0 = time.time()\n\n        model = Lasso(alpha = alpha )\n\n        model.fit(X[IX_train,:], y[IX_train])\n\n        col = i\n        df_models_1.loc[col,'alpha'] = alpha\n        y_pred = model.predict(X[IX_train])\n        c = np.corrcoef(y[IX_train], y_pred)[0,1]\n        print(alpha, 'Corr Coef Train', c)\n        df_models_1.loc[col,'Corr Coef Train'] = c\n\n        y_pred = model.predict(X[IX_test])\n        c = np.corrcoef(y[IX_test], y_pred)[0,1]\n        print(alpha, 'Corr Coef Test', c)\n        df_models_1.loc[col,'Corr Coef'] = c\n\n        df_models_1.loc[col,'n_nonzeros'] = (model.coef_ != 0 ).sum()\n\n        print('%.1f secs passed'%(time.time()-t0))\n\n    alpha_selected = df_models_1.sort_values('Corr Coef', ascending = False)['alpha'].iat[0] \n    print('Best alpha: ', alpha_selected )    \n    display( df_models_1 )  \n    \n    \nprint('Best alpha: ', alpha_selected )    \n","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:44:24.916692Z","iopub.execute_input":"2023-01-10T16:44:24.917612Z","iopub.status.idle":"2023-01-10T16:44:24.940493Z","shell.execute_reply.started":"2023-01-10T16:44:24.917564Z","shell.execute_reply":"2023-01-10T16:44:24.939634Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Enforce alpha value \n# alpha_selected =  0.1\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('Best alpha: ',  alpha_selected )    ","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:44:24.941845Z","iopub.execute_input":"2023-01-10T16:44:24.942453Z","iopub.status.idle":"2023-01-10T16:44:24.955959Z","shell.execute_reply.started":"2023-01-10T16:44:24.942420Z","shell.execute_reply":"2023-01-10T16:44:24.954732Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nprint('That standard  way with LassoCV - causes RAM crash. So we do not use it and  search for alpha by direct loop')  \nif 0:\n    from sklearn.linear_model import LassoCV\n    from sklearn.linear_model import RidgeCV\n    from sklearn.linear_model import Ridge\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.metrics import mean_squared_error\n\n    #alpha_selected = 1e4\n\n    n_splits_for_cross_valdition = 2\n    random_state_cross_validation = 0\n    kf = KFold(n_splits=n_splits_for_cross_valdition,  shuffle=True, random_state= random_state_cross_validation )\n\n    #model = RidgeCV(alphas= [ 1e6], cv= kf ).fit(X, y)\n    model = RidgeCV(alphas=[1e-3, 1e-2, 1e-1, 1,1e1,1e2,1e3,1e4,1e5,1e6,1e7], cv= kf ).fit(X, y) # Wall time: 1h 10min 36s\n\n    print(model)\n    alpha_selected = model.alpha_\n    print( alpha_selected )\n    print()","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:44:24.957505Z","iopub.execute_input":"2023-01-10T16:44:24.958296Z","iopub.status.idle":"2023-01-10T16:44:24.969688Z","shell.execute_reply.started":"2023-01-10T16:44:24.958249Z","shell.execute_reply":"2023-01-10T16:44:24.968574Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nif not fast_mode: \n    \n    model = Lasso(alpha = alpha_selected )\n    model.fit(X, y)\n    print('Count non zero coefs for training on the full sample:', (model.coef_ != 0 ).sum()  )\n    df_model_on_full_sample =  pd.DataFrame(index = df_rna.columns)\n    df_model_on_full_sample['Coefs'] = model.coef_\n    df_model_on_full_sample = df_model_on_full_sample.sort_values('Coefs', key = abs, ascending = False )\n    display( df_model_on_full_sample.head(30) )\n    m = df_model_on_full_sample['Coefs'] != 0\n    display( df_model_on_full_sample[m].tail(15) )\n\n\n    fn = target_name+'_model_on_full_sample_'+str_method + '.csv'\n    print(fn); print()\n    df_model_on_full_sample.to_csv(fn)\n\n\n    m = df_model_on_full_sample['Coefs'] != 0\n    print('Total non-zero coefs (model on full data):', m.sum() )\n    fig = plt.figure(figsize = (20,5))\n    plt.plot( df_model_on_full_sample['Coefs'][m].abs().values, '*-', label = 'Abs Coefs')\n    plt.grid()\n    plt.legend()\n    plt.title('Sorted Lasso abs coefficients train on the full sample. Only non-zero coefs', fontsize = 20 )\n    plt.xlabel('Coefs', fontsize = 20)\n    plt.xlabel('Index', fontsize = 20)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:44:24.971356Z","iopub.execute_input":"2023-01-10T16:44:24.972028Z","iopub.status.idle":"2023-01-10T16:44:24.994723Z","shell.execute_reply.started":"2023-01-10T16:44:24.971991Z","shell.execute_reply":"2023-01-10T16:44:24.993512Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Modeling on different random subsamples ","metadata":{}},{"cell_type":"code","source":"%%time\n\nimport time\n\ndf_importances = pd.DataFrame(index = df_rna.columns)\ndf_models = pd.DataFrame()\n\nprint('alpha_selected:', alpha_selected , 'model:', str_method)\n\nfor trial in range( n_trials):\n    t0 = time.time()\n\n    p = np.random.permutation(len(y))\n    N = int(len(y) * train_size )\n    IX_train = np.arange(len(y))[p][:N]\n    IX_test = np.arange(len(y))[p][N:]\n    print('Trial:', trial, \" len(IX_train), len(IX_test):\",  len(IX_train), len(IX_test), '%.1f secs passed'%(time.time() - t0 ) )\n\n    #model = Ridge(alpha = alpha_selected )\n    model = Lasso(alpha = alpha_selected )\n    model.fit(X[IX_train,:], y[IX_train])\n\n    \n    col = str_method + 'Trial'+str(trial)\n    df_importances[col] = model.coef_\n\n    y_pred = model.predict(X[IX_test])\n    c = np.corrcoef(y[IX_test], y_pred)[0,1]\n    print('Corr Coef:',c)\n    df_models.loc[col,'Corr Coef'] = c\n    \n    y_pred = model.predict(X[IX_train])\n    c = np.corrcoef(y[IX_train], y_pred)[0,1]\n    # print(c)\n    df_models.loc[col,'Corr Coef Train'] = c\n    \n    df_models.loc[col,'n_nonzeros'] = (model.coef_ != 0 ).sum()\n    \n    print('Trial', trial, '%.1f secs passed'%(time.time() - t0 ))\n    \ndisplay(df_models )    \ndisplay(df_importances.sort_values( df_importances.columns[0],ascending = False ).head(6))\nmedian_importances = df_importances.median(axis = 1) \nmedian_importances = median_importances.sort_values(key = abs, ascending = False)\nmedian_importances    ","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:44:24.996589Z","iopub.execute_input":"2023-01-10T16:44:24.997214Z","iopub.status.idle":"2023-01-10T16:46:13.372591Z","shell.execute_reply.started":"2023-01-10T16:44:24.997147Z","shell.execute_reply":"2023-01-10T16:46:13.370968Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_importances.to_csv(target_name+'_importances_'+str_method+'_Trials' + str(df_importances.shape[1]) + '.csv')","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:46:13.375583Z","iopub.execute_input":"2023-01-10T16:46:13.377032Z","iopub.status.idle":"2023-01-10T16:46:13.446728Z","shell.execute_reply.started":"2023-01-10T16:46:13.376953Z","shell.execute_reply":"2023-01-10T16:46:13.445494Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_models.to_csv(target_name+'_models_stat_'+str_method+'_Trials' + str(df_importances.shape[1]) + '.csv')","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:46:13.448641Z","iopub.execute_input":"2023-01-10T16:46:13.449090Z","iopub.status.idle":"2023-01-10T16:46:13.456828Z","shell.execute_reply.started":"2023-01-10T16:46:13.449042Z","shell.execute_reply":"2023-01-10T16:46:13.455431Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_models.describe()","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:46:13.459318Z","iopub.execute_input":"2023-01-10T16:46:13.459982Z","iopub.status.idle":"2023-01-10T16:46:13.489665Z","shell.execute_reply.started":"2023-01-10T16:46:13.459942Z","shell.execute_reply":"2023-01-10T16:46:13.488426Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize = (20,5))\nfor col in ['Corr Coef', 'Corr Coef Train'  ]:\n    if col in df_models.columns:\n        plt.plot( df_models[col].values, label = col)\n    \nplt.grid()\nplt.legend()\nplt.title('Correlation of predictions with targets', fontsize = 20 )\nplt.xlabel('Correlation', fontsize = 20)\nplt.xlabel('Trial', fontsize = 20)\nplt.show()\n\nfig = plt.figure(figsize = (20,5))\nplt.plot( df_models['n_nonzeros'].values, label = 'n_nonzeros')\nplt.grid()\nplt.legend()\nplt.title('Count of non-zeros coefficients for '+str_method, fontsize = 20 )\nplt.xlabel('Count', fontsize = 20)\nplt.xlabel('Trial', fontsize = 20)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:46:13.491329Z","iopub.execute_input":"2023-01-10T16:46:13.491777Z","iopub.status.idle":"2023-01-10T16:46:14.028892Z","shell.execute_reply.started":"2023-01-10T16:46:13.491730Z","shell.execute_reply":"2023-01-10T16:46:14.027593Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pd.concat([ median_importances.iloc[:20].to_frame().reset_index(), \nmedian_importances.iloc[20:40].to_frame().reset_index(), median_importances.iloc[40:60].to_frame().reset_index()\n          , median_importances.iloc[60:80].to_frame().reset_index() ], axis = 1)","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:46:14.030421Z","iopub.execute_input":"2023-01-10T16:46:14.030809Z","iopub.status.idle":"2023-01-10T16:46:14.061997Z","shell.execute_reply.started":"2023-01-10T16:46:14.030772Z","shell.execute_reply":"2023-01-10T16:46:14.060560Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"m = (median_importances != 0 )\nprint('Count Non-zero median importances ', m.sum())\n","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:46:14.063279Z","iopub.execute_input":"2023-01-10T16:46:14.063760Z","iopub.status.idle":"2023-01-10T16:46:14.071618Z","shell.execute_reply.started":"2023-01-10T16:46:14.063704Z","shell.execute_reply":"2023-01-10T16:46:14.070234Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"m = (median_importances != 0 )\nprint(m.sum())\nfig = plt.figure(figsize = (20,5))\nplt.plot( median_importances[m].abs().values, '*-',label = 'Abs Coefs')\nplt.grid()\nplt.legend( fontsize = 20)\nplt.title('Sorted '+str_method+' abs importances . Only non-zero.   Median over trials '+str(n_trials), fontsize = 20 )\nplt.xlabel('Coefs', fontsize = 20)\nplt.xlabel('Index', fontsize = 20)\nplt.show()\n\nm = (median_importances != 0 )\nfig = plt.figure(figsize = (20,5))\nplt.plot( median_importances[m].abs().values[1:100],'*-', label = 'Abs Coefs', )\nplt.grid()\nplt.legend( fontsize = 20)\nplt.title('Sorted  '+str_method+'  abs importances 1:100. Only non-zero.   Median over trials '+str(n_trials), fontsize = 20 )\nplt.xlabel('Coefs', fontsize = 20)\nplt.xlabel('Index', fontsize = 20)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:46:14.079402Z","iopub.execute_input":"2023-01-10T16:46:14.079790Z","iopub.status.idle":"2023-01-10T16:46:14.621019Z","shell.execute_reply.started":"2023-01-10T16:46:14.079750Z","shell.execute_reply":"2023-01-10T16:46:14.619667Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Intersection of non-zero importances for all trials ","metadata":{}},{"cell_type":"code","source":"%%time\nprint('df_importances.shape', df_importances.shape)\n\nm = (median_importances != 0 )\nprint('Count Non-zero median importances ', m.sum())\n\ncc  = (df_importances != 0).sum(axis = 0).mean()  ##  == df_importances.shape[1]\nprint('Average count genes non-zero for trials: ', cc )\n\nm  = (df_importances != 0).sum(axis = 1) == df_importances.shape[1]\nprint('Count genes non-zero for all trials: ', m.sum())\n#print()\n\nif df_importances.shape[1] <= 10:\n    l = []\n    for col1 in df_importances.columns:\n        for col2 in df_importances.columns:\n            m  = (df_importances[[col1,col2]] != 0).sum(axis = 1) == df_importances[[col1,col2]].shape[1]\n            l.append(m.sum())\n    print('number of non-zero coeffients common for Pairs of trials:  ' , np.mean(l))        \n    #print()\n\nif df_importances.shape[1] <= 10:\n    import itertools\n\n    l = []\n    for cols in itertools.combinations( df_importances.columns , 3 ) :\n        m  = (df_importances[list(cols) ] != 0).sum(axis = 1) == df_importances[list(cols) ].shape[1]\n        l.append(m.sum())\n    print('number of non-zero coeffients common for TRIPLES of trials:  ' , np.mean(l))        \n    \nif df_importances.shape[1] <= 5:\n    import itertools\n    l = []\n    for cols in itertools.combinations( df_importances.columns , 4 ) :\n        m  = (df_importances[list(cols) ] != 0).sum(axis = 1) == df_importances[list(cols) ].shape[1]\n        l.append(m.sum())\n    print('number of non-zero coeffients common for QUADRUPLETS of trials:  ' , np.mean(l))        \n    for k in range(5, np.min([10, df_importances.shape[1]+1 ]) ):\n        l = []\n        for cols in itertools.combinations( df_importances.columns , k ) :\n            m  = (df_importances[list(cols) ] != 0).sum(axis = 1) == df_importances[list(cols) ].shape[1]\n            l.append(m.sum())\n        print('number of non-zero coeffients common for ', str(k)+'-tuples of trials:  ' , np.mean(l))        \n\n        \nprint()\nprint('Genes with non-zero importance for all trials: ' )\nprint( df_importances[m].sort_values(df_importances.columns[0], ascending = False, key = abs ).index[:200] )\nprint( [t.split('_')[1] for t in df_importances[m].sort_values(df_importances.columns[0], ascending = False, key = abs ).index[:200] ] )\n","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:46:14.622477Z","iopub.execute_input":"2023-01-10T16:46:14.622827Z","iopub.status.idle":"2023-01-10T16:46:14.654386Z","shell.execute_reply.started":"2023-01-10T16:46:14.622794Z","shell.execute_reply":"2023-01-10T16:46:14.653147Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nl = []\nfor k in range(1, df_importances.shape[1]+1):\n    m  = (df_importances.iloc[:,:k] != 0).sum(axis = 1) == k # df_importances.i.shape[1]\n    l.append(m.sum())\n        \nplt.show()\n\nfig = plt.figure(figsize = (20,5))\nplt.plot(l, '*-')\nplt.grid()\nplt.legend( fontsize = 20)\nplt.title('Count of non-zeros elements common in the first Index trials ', fontsize = 20 )\nplt.xlabel('Count non-zero', fontsize = 20)\nplt.xlabel('Index', fontsize = 20)\nplt.show()\n\nprint()\nprint('Counts of intersections of non-zero genes list  for the first k trials for all k:')\nprint('All:', l)\nprint('First 30:', l[:30])\nprint('Last 30:', l[-30:])\nprint()\nprint('Interesecting first k trials gives common non-zero genes (for selected \"k\"):')\nfor k in [0,10, 20, 30, 40, 50, 100]:\n    if len(l)>k: print('for k=',k+2,' common non-zero genes:', l[k-1]  ) \nprint()\nprint(m.sum(), 'genes non-zero for all trials (show not to more than 300):')\ndd = df_importances[m].median(axis = 1).sort_values(ascending = False, key = abs)\nprint(dd.index[:300]  )\nprint()\nprint('Top15:')\ndisplay( dd.head(15).to_frame() ) \nprint()\nprint('Tail15:')\ndisplay( dd.tail(15).to_frame() ) \ndd.to_csv(target_name+'_non_zero_median_sorted_'+str_method+'_Trials' + str(df_importances.shape[1]) + '.csv')\n        \ni_start = 10\nfig = plt.figure(figsize = (20,5))\nplt.plot(range(i_start, len(l)), l[i_start:], '*-')\nplt.grid()\nplt.legend( fontsize = 20)\nplt.title('Count of non-zeros elements common in the first '+str(i_start ) + '+Index trials ', fontsize = 20 )\nplt.xlabel('Count non-zero', fontsize = 20)\nplt.xlabel('Index', fontsize = 20)\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:46:14.655934Z","iopub.execute_input":"2023-01-10T16:46:14.656839Z","iopub.status.idle":"2023-01-10T16:46:15.169943Z","shell.execute_reply.started":"2023-01-10T16:46:14.656800Z","shell.execute_reply":"2023-01-10T16:46:15.168505Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"(df_importances.iloc[:,1] != 0).sum()","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:46:15.171800Z","iopub.execute_input":"2023-01-10T16:46:15.172194Z","iopub.status.idle":"2023-01-10T16:46:15.181130Z","shell.execute_reply.started":"2023-01-10T16:46:15.172135Z","shell.execute_reply":"2023-01-10T16:46:15.179776Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"N0 = df_importances.shape[1]\nN2 = int(N0  /2 )\nN2\nm1 = (df_importances.iloc[:,:N2] != 0 ).sum(axis = 1) == N2\nm2 = (df_importances.iloc[:,N2:] != 0 ).sum(axis = 1) == (N0 - N2)\n\nprint('Count common non-zeros for the first half and second half of trials: ', m1.sum() , m2.sum()  )\ns = set(df_importances.index[m1]) &  set(df_importances.index[m2]) \nprint('Count interesection between the first and the second halves: ', len(s) )\nsr = df_importances[m1&m2].median(axis = 1).sort_values( ascending = False , key = abs )\nprint('Top 20 of intersection:')\ndisplay( sr.head(20).to_frame() )\nprint('Tail 15 of intersection:')\ndisplay( sr.tail(15).to_frame() )","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:46:15.183598Z","iopub.execute_input":"2023-01-10T16:46:15.184196Z","iopub.status.idle":"2023-01-10T16:46:15.219201Z","shell.execute_reply.started":"2023-01-10T16:46:15.184108Z","shell.execute_reply":"2023-01-10T16:46:15.217933Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Analysis of the obtained importances ","metadata":{}},{"cell_type":"code","source":"df_importances.describe()","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:46:15.220612Z","iopub.execute_input":"2023-01-10T16:46:15.220945Z","iopub.status.idle":"2023-01-10T16:46:15.243872Z","shell.execute_reply.started":"2023-01-10T16:46:15.220914Z","shell.execute_reply":"2023-01-10T16:46:15.242561Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_importances.corr()","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:46:15.245549Z","iopub.execute_input":"2023-01-10T16:46:15.245902Z","iopub.status.idle":"2023-01-10T16:46:15.259622Z","shell.execute_reply.started":"2023-01-10T16:46:15.245870Z","shell.execute_reply":"2023-01-10T16:46:15.258273Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nn_x_subplots = 5\nc = 0\n\nfor i in range(5):\n    if i >= df_importances.shape[1]: continue\n    for j in range(5):\n        if i <= j: continue \n        if j >= df_importances.shape[1]: continue\n        col1 = df_importances.columns[i]\n        col2 = df_importances.columns[j]\n\n        if c % n_x_subplots == 0:\n            if c > 0:\n                plt.show()\n            fig = plt.figure(figsize = (20,5) ); c = 0\n            #plt.suptitle(str(k),fontsize = 20 )# str_data_inf + ' n_cells: ' +str(mask.sum()) + ' RED > median expression, BLUE <= median ' )# +' ' + cell_type +' ' + drug )\n\n        c += 1; fig.add_subplot(1,n_x_subplots ,c)\n\n        sns.scatterplot(x = df_importances[col1], y=df_importances[col2] )\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:46:15.262082Z","iopub.execute_input":"2023-01-10T16:46:15.262471Z","iopub.status.idle":"2023-01-10T16:46:15.522557Z","shell.execute_reply.started":"2023-01-10T16:46:15.262438Z","shell.execute_reply":"2023-01-10T16:46:15.521359Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Calculations of intersections between ordered importances ","metadata":{}},{"cell_type":"code","source":"%%time \ndict_counts = {}\n\ndict_sorted = {} # pd.DataFrame()\nfor col in df_importances.columns:\n    dict_sorted[col] = df_importances[col].abs().sort_values(ascending = False )\n\nfor k in range(0, df_importances.shape[0]):\n    if (k%5000 == 1): print(k)\n    col = df_importances.columns[0]\n    s = set(  dict_sorted[col].index[:k] )\n    for col in df_importances.columns:\n        \n        s = s & set(dict_sorted[col].index[:k] )\n        if col not in dict_counts.keys(): \n            dict_counts[col] = []\n        dict_counts[col].append(len(s))\n\nplt.figure(figsize = (20,10))\nfor i,col in enumerate(dict_counts):\n    if len(dict_counts ) >= 50:\n        if (i%10) != 0: continue\n    plt.plot(dict_counts[col], label = 'intersection till: ' + col )\nplt.legend(fontsize = 20 )\nplt.grid()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:46:15.524413Z","iopub.execute_input":"2023-01-10T16:46:15.524767Z","iopub.status.idle":"2023-01-10T16:48:39.591988Z","shell.execute_reply.started":"2023-01-10T16:46:15.524732Z","shell.execute_reply":"2023-01-10T16:48:39.590627Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Plots of parabolic approximation for full list of features","metadata":{}},{"cell_type":"code","source":"for i,col in  enumerate(dict_counts): # = '3 32606 EryP'\n    \n    if len(dict_counts ) >= 50:\n        if (i%10) != 0: continue \n\n    ll = dict_counts[col]\n    plt.figure(figsize = (20,6))\n    plt.plot(ll, label = 'intersection till: ' + col)\n    y = np.array(ll)\n    x = np.arange(len(y))\n    p = np.polyfit(x,y,2)\n    print(p)\n    plt.plot(np.polyval(p,x) )\n    plt.legend(fontsize = 20 )\n    plt.grid()\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:48:39.593630Z","iopub.execute_input":"2023-01-10T16:48:39.594100Z","iopub.status.idle":"2023-01-10T16:48:40.210515Z","shell.execute_reply.started":"2023-01-10T16:48:39.594043Z","shell.execute_reply":"2023-01-10T16:48:40.209236Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Plots of linear and parabolic approximations for ONLY top genes ","metadata":{}},{"cell_type":"code","source":"list_slopes = []\nfor i, col in  enumerate(dict_counts): # = '3 32606 EryP'\n    ll = dict_counts[col]\n    \n    y = np.array(ll)\n    x = np.arange(len(y))\n    \n    M = 100\n    p = np.polyfit(x[:M],y[:M],1)\n    print(p)\n    list_slopes.append(p[0])\n    \n    if len(dict_counts ) >= 50:\n        if (i%10) != 0: continue \n        \n    plt.figure(figsize = (20,6))\n    plt.plot(ll, label = 'intersection till: ' + col)\n    y = np.array(ll)\n    x = np.arange(len(y))\n    p = np.polyfit(x,y,2)\n    print(p)\n    plt.plot(np.polyval(p,x), label = 'quadratic approx'  )\n\n    M = 100\n    p = np.polyfit(x[:M],y[:M],1)\n    print(p)\n    plt.plot(np.polyval(p,x), label = 'linear approx' )\n\n    \n    if i < 3:\n        plt.ylim([0,5000])\n    else:\n        plt.ylim([0,200])\n    plt.xlim([0,5000])\n    plt.legend(fontsize = 20 )\n    plt.grid()\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:48:40.211989Z","iopub.execute_input":"2023-01-10T16:48:40.212365Z","iopub.status.idle":"2023-01-10T16:48:40.878395Z","shell.execute_reply.started":"2023-01-10T16:48:40.212330Z","shell.execute_reply":"2023-01-10T16:48:40.877421Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Analysis of slopes (of linear approximations) changes with trial number changes  - is there stabilization or not ? ","metadata":{}},{"cell_type":"code","source":"    \nplt.figure(figsize = (20,10))\nplt.plot( list_slopes     )\nplt.title('Slopes for linear approximation ', fontsize = 20)\nplt.grid()\nplt.show()\n\nplt.figure(figsize = (20,10))\nx = np.arange( 100,len(list_slopes) )\nplt.plot(x,  np.array(list_slopes)[x]     )\nplt.title('Slopes for linear approximation ' , fontsize = 20)\nplt.grid()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:48:40.879683Z","iopub.execute_input":"2023-01-10T16:48:40.880545Z","iopub.status.idle":"2023-01-10T16:48:41.401125Z","shell.execute_reply.started":"2023-01-10T16:48:40.880506Z","shell.execute_reply":"2023-01-10T16:48:41.399699Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"median_importances = df_importances.median(axis = 1 )\nmedian_importances_sorted = median_importances.sort_values(ascending = False, key = abs)\nmedian_importances_sorted","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:48:41.402991Z","iopub.execute_input":"2023-01-10T16:48:41.404191Z","iopub.status.idle":"2023-01-10T16:48:41.427755Z","shell.execute_reply.started":"2023-01-10T16:48:41.404119Z","shell.execute_reply":"2023-01-10T16:48:41.426563Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \n\ndict_top_stable_features = {}\ndf_top_stable_features_counts = pd.DataFrame(); IX = 0\n\ndict_sorted = {} # pd.DataFrame()\nfor col in df_importances.columns:\n    dict_sorted[col] = df_importances[col].abs().sort_values(ascending = False )\n\nfor k in [1, 5, 10, 20, 30,40, 50 , 100, 200,300,500, 1000, 2000, 5000]: #  range(0, df_importances.shape[0]):\n    #if (k%5000 == 1): \n    #print(k)\n    col = df_importances.columns[0]\n    s = set(  dict_sorted[col].index[:k] )\n    for col in df_importances.columns:\n        \n        s = s & set(dict_sorted[col].index[:k] )\n\n    s2 = s & set(median_importances_sorted.index[:k])    \n    m = median_importances_sorted.index.isin(s)\n    list_top_stable_features_ordered_by_median_importances = list( median_importances_sorted[m].index )\n    dict_top_stable_features[k] = list_top_stable_features_ordered_by_median_importances\n    df_top_stable_features_counts.loc[IX, 'Top K' ] = k\n    df_top_stable_features_counts.loc[IX, 'Intersection' ] = len(s)\n    df_top_stable_features_counts.loc[IX, 'Intersection With Median Importances' ] = len(s2)\n    IX += 1 \n    \n    print('Interesection top ',k, 'features for all trials gives: ', len(s) ,' in common', 'Intersection with median importances', len(s2) )\n    \nplt.figure(figsize = (10,5 ))\nsns.lineplot(data = df_top_stable_features_counts,  x = 'Top K', y = 'Intersection'   )\nsns.lineplot(data = df_top_stable_features_counts,  x = 'Top K', y = 'Intersection With Median Importances'  )\nplt.legend()\nplt.grid()\nplt.show()\ndisplay(df_top_stable_features_counts)\n","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:48:41.429149Z","iopub.execute_input":"2023-01-10T16:48:41.429548Z","iopub.status.idle":"2023-01-10T16:48:41.762899Z","shell.execute_reply.started":"2023-01-10T16:48:41.429515Z","shell.execute_reply":"2023-01-10T16:48:41.761648Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Print top 100 for intersections obtained for different K","metadata":{}},{"cell_type":"code","source":"for k in dict_top_stable_features:\n    l = dict_top_stable_features[k]\n    print(k, len(l))\n    print(l[:100])\n    print()","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:48:41.764325Z","iopub.execute_input":"2023-01-10T16:48:41.764660Z","iopub.status.idle":"2023-01-10T16:48:41.771746Z","shell.execute_reply.started":"2023-01-10T16:48:41.764629Z","shell.execute_reply":"2023-01-10T16:48:41.770398Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pd.set_option('display.max_rows', 500)\npd.set_option('display.max_columns', 500)\npd.set_option('display.width', 1000)","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:48:41.773337Z","iopub.execute_input":"2023-01-10T16:48:41.773674Z","iopub.status.idle":"2023-01-10T16:48:41.785234Z","shell.execute_reply.started":"2023-01-10T16:48:41.773643Z","shell.execute_reply":"2023-01-10T16:48:41.783749Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Dataframe with stable topK features for different K ","metadata":{}},{"cell_type":"code","source":"mx = 0\nfor k in dict_top_stable_features:\n    l = dict_top_stable_features[k]\n    mx = max([len(l), mx])\nprint('max len:', mx)\n\ndf_stable_top_features = pd.DataFrame()\n\nfor k in dict_top_stable_features:\n    df_stable_top_features['Top ' +str(k) + ' Intersection'] = [np.nan]*mx\n    \n    l = dict_top_stable_features[k]\n    df_stable_top_features['Top ' +str(k) + ' Intersection'].iloc[:len(l)] = l\ndf_stable_top_features.head(100)\n","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:48:41.786634Z","iopub.execute_input":"2023-01-10T16:48:41.787618Z","iopub.status.idle":"2023-01-10T16:48:41.885335Z","shell.execute_reply.started":"2023-01-10T16:48:41.787574Z","shell.execute_reply":"2023-01-10T16:48:41.884034Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"str_method, n_trials","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:48:41.886770Z","iopub.execute_input":"2023-01-10T16:48:41.887650Z","iopub.status.idle":"2023-01-10T16:48:41.895267Z","shell.execute_reply.started":"2023-01-10T16:48:41.887613Z","shell.execute_reply":"2023-01-10T16:48:41.893952Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfn = target_name+'_stable_top_features_'+str_method+'_n_trials_'+str(n_trials)+'.csv'\nprint(fn); print()\ndf_stable_top_features.to_csv(fn)","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:48:41.896604Z","iopub.execute_input":"2023-01-10T16:48:41.897103Z","iopub.status.idle":"2023-01-10T16:48:41.927571Z","shell.execute_reply.started":"2023-01-10T16:48:41.897063Z","shell.execute_reply":"2023-01-10T16:48:41.926338Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tt = time.time() - t0start\nprint('%.1f seconds ( = %.1f minutes, = %.2f hours) passed'%( tt, tt/60, tt/3600 ) )","metadata":{"execution":{"iopub.status.busy":"2023-01-10T16:48:41.928868Z","iopub.execute_input":"2023-01-10T16:48:41.929285Z","iopub.status.idle":"2023-01-10T16:48:41.935527Z","shell.execute_reply.started":"2023-01-10T16:48:41.929247Z","shell.execute_reply":"2023-01-10T16:48:41.934501Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}