{"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\nCalculate univariable feature importances and save to csv. \n\n### Versions\n\n####  71,72,73,74,75,76, 77 , 78 Rerurn timilimit crashed  'MutualInfSklearn' Randomly Permuted\n\n####  62,63,64,65,66,67  Rerun timelimit crashed 'MutualInfSklearn'\n\n\n####  56-61  'MutualInfSklearn' Randomly Permuted\n\n#### 55 - quick save\n\n####  48,49,50,51,52,53,54  'MutualInfSklearn' CD: 0-20, 20-40,40-60,60-80, 80-100, 100-120,120-140\n\n    Three version crashed by timelimit 12 hours. \n    \n    \n#### bugs:  42,43,44,45,46,47  'MutualInfSklearn' CD: 0-20, 20-40,40-60,60-80, 80-100, 100-120,120-140\n\n#### BUGGY: do not use:  35,36,37,38,39,40,41 Random Permutation for 'MutualInfSklearn', CD: 0-20, 20-40,40-60, 60-80, 80-100,100-120,120-140\n\n####  BUGGY: do not use:  24,25,26,27,28,29,30,31,32,33,34  'MutualInfSklearn' CD: 0-10, 10-20, 20-30,30-40,40-50,50-60,60-70 , 70-90, 90-110, 110-130,130-140\n\n    About 10 features per second for downsample = 80 \n    down_sample_percent = 80\n    'MutualInfSklearn'  - is very slow \n\n####  BUGGY: do not use:  23  'MutualInfSklearn'\n    It is very slow - run only 4 CD\n    About 4 features per second for downsample = 50\n    down_sample_percent = 50\n    \n#### 22 - bug\n\n#### 21 Spearman \n\n    down_sample_percent = 50\n\n#### 20 code update - introduce Pearson, LGB, etc correlations \n\n    Run Pearson correlation \n    \n    down_sample_percent = 50\n\n#### 18,19 LGB: n_estimators = 5 , RANDOM PERMUTATION to estimate p-valies\n\n    n_estimators = 5\n    \n    down_sample_percent = 50\n\n#### 16,17 LGB: n_estimators = 5 , repeat with other random seed for LGB and permututation changing for each CD \n\n    n_estimators = 5\n    \n    down_sample_percent = 50\n\n#### 14,15 LGB: n_estimators = 1 , repeat with other random seed for LGB and permututation changing for each CD \n\n    n_estimators = 1 \n    \n    down_sample_percent = 50\n\n#### 12,13 LGB: split calcilation into 2 parts\n\n    n_estimators = 1  \n    \n    down_sample_percent = 50\n\n#### 10,11 - bug\n\n#### 6,7,8,9 LGB: split calcilation into 4 parts\n\n    n_estimators = 10  \n    \n    down_sample_percent = 50\n    \n\n#### 4,5 LGB: V4 - targets to 70, V5 - after 70 - expected time about 10 hours - cancelled by time limit 12 hours\n\n    n_estimators = 10  - it should take longer time \n\n    down_sample_percent = 50\n\n#### 2,3 LGB: V2 - targets to 70, V3 - after 70 - expected time about 6 hours - took 10 hours in reality\n\n    down_sample_percent = 50\n    n_estimators = 5\n\n#### 1 LGB test run: with two targets only\n","metadata":{}},{"cell_type":"markdown","source":"# Key params","metadata":{}},{"cell_type":"code","source":"# target_name = 'CD146'#   'CD194' # 'CD44'# 'CD62L'\n\nstr_dataset_name = 'NIPS22'\n\nstr_method_name = 'MutualInfSklearn' #  'Spearman'# 'Pearson' #  'LGB'                \n# 'MutualInfSklearn' - is very slow 4 features per second for downsample = 50% ( n_samples = 70_00)\n# Downsample 75% gives about 9 features per second  - that is better \n\ncompute_for_randomly_permuted_data = 1 # = 1 will make random permutation - to estimate p-values , not for actual feature importances \n\ndown_sample_percent = 80\n\n# For boostings: \nn_estimators = 5\nrandom_state4lgb = 10\n\n\nverbose =  10_000 # 10_000 - Default # Indicating progress of the main calculation each verbose time of features processed \n\n\nimport time\nt0start = time.time() ","metadata":{"execution":{"iopub.status.busy":"2023-02-25T08:41:58.406149Z","iopub.execute_input":"2023-02-25T08:41:58.406920Z","iopub.status.idle":"2023-02-25T08:41:58.439933Z","shell.execute_reply.started":"2023-02-25T08:41:58.406804Z","shell.execute_reply":"2023-02-25T08:41:58.438829Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"n_CD_to_start = 136# Sometimes we need to split calculation by parts otherwise crash by timelimit \nn_CD_to_finish = 140\nstr_technical_postfix = '_Part'+str(n_CD_to_start)+'_'+str(n_CD_to_finish)\nif (n_CD_to_start == 0) and (n_CD_to_finish == 140): str_technical_postfix = ''\n\nfn4save = str_dataset_name+'_'+'importances_univariable_'+str(str_method_name)+'_down_sample_percent'+str(down_sample_percent)\nif str_method_name == 'LGB':\n    fn4save += '_n_estimators'+str(n_estimators) +'_RandSeed'+str(random_state4lgb)\nif compute_for_randomly_permuted_data > 0:\n    fn4save += '_RandomPermuted'\nfn4save += str_technical_postfix+'.csv'\nfn4save","metadata":{"execution":{"iopub.status.busy":"2023-02-25T08:41:58.442030Z","iopub.execute_input":"2023-02-25T08:41:58.442417Z","iopub.status.idle":"2023-02-25T08:41:58.455601Z","shell.execute_reply.started":"2023-02-25T08:41:58.442385Z","shell.execute_reply":"2023-02-25T08:41:58.454509Z"},"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-02-25T08:41:58.456665Z","iopub.execute_input":"2023-02-25T08:41:58.457029Z","iopub.status.idle":"2023-02-25T08:41:58.533378Z","shell.execute_reply.started":"2023-02-25T08:41:58.456995Z","shell.execute_reply":"2023-02-25T08:41:58.531900Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# pd.set_option('display.max_rows', 500)\n# pd.set_option('display.max_columns', 500)\n# pd.set_option('display.width', 1000)","metadata":{"execution":{"iopub.status.busy":"2023-02-25T08:41:58.535774Z","iopub.execute_input":"2023-02-25T08:41:58.536158Z","iopub.status.idle":"2023-02-25T08:41:58.541354Z","shell.execute_reply.started":"2023-02-25T08:41:58.536123Z","shell.execute_reply":"2023-02-25T08:41:58.539834Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pip install git+https://github.com/jundongl/scikit-feature","metadata":{"execution":{"iopub.status.busy":"2023-02-25T08:41:58.542871Z","iopub.execute_input":"2023-02-25T08:41:58.543283Z","iopub.status.idle":"2023-02-25T08:42:19.930620Z","shell.execute_reply.started":"2023-02-25T08:41:58.543249Z","shell.execute_reply":"2023-02-25T08:42:19.929182Z"},"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-02-25T08:42:19.933872Z","iopub.execute_input":"2023-02-25T08:42:19.934423Z","iopub.status.idle":"2023-02-25T08:42:20.626160Z","shell.execute_reply.started":"2023-02-25T08:42:19.934372Z","shell.execute_reply":"2023-02-25T08:42:20.624850Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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\nfrom skfeature.utility.entropy_estimators import midd\nimport numpy as np\nimport time","metadata":{"execution":{"iopub.status.busy":"2023-02-25T08:42:20.627777Z","iopub.execute_input":"2023-02-25T08:42:20.628192Z","iopub.status.idle":"2023-02-25T08:42:20.859081Z","shell.execute_reply.started":"2023-02-25T08:42:20.628157Z","shell.execute_reply":"2023-02-25T08:42:20.857736Z"},"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-02-25T08:42:20.861066Z","iopub.execute_input":"2023-02-25T08:42:20.861629Z","iopub.status.idle":"2023-02-25T08:44:14.338771Z","shell.execute_reply.started":"2023-02-25T08:42:20.861501Z","shell.execute_reply":"2023-02-25T08:44:14.337428Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print( dict( df_meta.value_counts('cell_type') ))\n\n# for t in df_rna.columns:\n#     if 'CD44' in t:\n#         print(t)\n#     if 'ENSG00000197405' in t: # ENSG00000197405_C5AR1 CD88\n#         print(t)\n        \n        \n# 'ENSG00000135218_CD36', 'ENSG00000026508_CD44', 'ENSG00000255443_CD44-AS1'\n       ","metadata":{"execution":{"iopub.status.busy":"2023-02-25T08:44:14.340348Z","iopub.execute_input":"2023-02-25T08:44:14.340704Z","iopub.status.idle":"2023-02-25T08:44:14.354076Z","shell.execute_reply.started":"2023-02-25T08:44:14.340672Z","shell.execute_reply":"2023-02-25T08:44:14.353176Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# y = df_y[target_name].values\n# print(y.shape, type(y) )\n\nX = ( df_rna.values )\nprint(X.shape)\n\n# y = df_y[target_name]\n# print(y.shape)","metadata":{"execution":{"iopub.status.busy":"2023-02-25T08:44:14.357527Z","iopub.execute_input":"2023-02-25T08:44:14.358355Z","iopub.status.idle":"2023-02-25T08:44:14.363923Z","shell.execute_reply.started":"2023-02-25T08:44:14.358319Z","shell.execute_reply":"2023-02-25T08:44:14.363151Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Main calculation","metadata":{}},{"cell_type":"code","source":"%%time\n\nimport time\nimport catboost as cb\nimport lightgbm as lgbm\nfrom sklearn.model_selection import cross_val_predict\nfrom sklearn.model_selection import KFold\nfrom scipy import stats\nfrom sklearn.feature_selection import f_regression, mutual_info_regression\n\nimport time\nt0 = time.time()\n\ndf_importances = pd.DataFrame(index = df_rna.columns ) #  Dataframe to store results and target\n\nc = 0\nfor target_name in df_y.columns[n_CD_to_start:n_CD_to_finish]:\n    c += 1 \n    y = df_y[target_name]\n    print(c, target_name, y.shape, X.shape)\n\n\n    t0 = time.time()\n\n    l = []\n    for i,col in enumerate( df_rna.columns):\n    #     s = np.corrcoef(y.values[p][:N],df_rna[col].values[p][:N])[0,1]\n    #     df_importances.loc[col, target_name +' Corr Pearson DownSample'+str(down_sample_percent)  ]  = s\n\n        p = np.random.permutation(len(df_y))\n        N = int(down_sample_percent*len(df_y) / 100 ) # +str(down_sample_percent) \n\n        y_train = y.values[p][:N]\n        y_test =  y.values[p][N:]\n\n        if compute_for_randomly_permuted_data > 0:\n            # For estimation of the p-values - target and data are different by random permutation - so we expect near zero prediction score \n            X_train = pd.Series(df_rna[col].values[:N]).to_frame()\n            X_test = pd.Series(df_rna[col].values[N:]).to_frame()\n        else:\n            # Normal mode - Y and X are permutated by the same permutation - just a way to do random train-test split\n            X_train = pd.Series(df_rna[col].values[p][:N]).to_frame()\n            X_test = pd.Series(df_rna[col].values[p][N:]).to_frame()\n\n        if str_method_name == 'Pearson':                \n            y_pred = X_test.iloc[:,0].values \n            s = np.corrcoef(y_test,y_pred)[0,1]\n        elif str_method_name == 'Spearman':   \n            y_pred = X_test.iloc[:,0].values\n            res = stats.spearmanr(y_test, y_pred )\n            s = res.correlation\n        elif str_method_name ==  'MutualInfSklearn' :   \n            #y_pred = X_test.iloc[:,0].values\n            mi = mutual_info_regression(X_test, y_test)\n            s = mi[0]\n            #print(s, mi )\n        elif str_method_name ==  'MutualInfScikitFeature' :   \n            #y_pred = X_test.iloc[:,0].values\n            x_loc = df_rna[col].values[p][:N].round(1)\n            y = y.round(1)\n            mi = midd(x_loc , y.values[p][:N] )\n            s = mi\n            #print(s, mi )\n        elif str_method_name == 'LGB':                \n            model = lgbm.LGBMRegressor(random_state = random_state4lgb, n_estimators = n_estimators) # Ridge(alpha = alpha_selected )\n            model.fit(X_train,y_train)\n            y_pred = model.predict(X_test)# ,y_train)\n            s = np.corrcoef(y_test,y_pred)[0,1]\n            \n        if np.isnan(s): s = 0 \n        #df_importances.loc[col, target_name +' Corr LGB DownSample'+str(down_sample_percent) +' Trial'+str(trial)+' n_estimators'+str(n_estimators)  ]  = s\n        l.append(s)\n        if (verbose > 0) and (i%verbose == 0 ):\n            I = np.argmax(l)\n            #print(i, col, '%.1f seconds passed'%(time.time() - t0 ), 'Pearson Score', s )\n            print(i, col,'Mean %.6f, max %.6f'%( np.mean(l), np.max(l) ),  \n                  '%.1f seconds passed'%(time.time() - t0 ), 'Argmax:', I, df_rna.columns[I] )            \n\n    df_importances.loc[df_importances.index[:len(l)],  target_name] = l  ","metadata":{"execution":{"iopub.status.busy":"2023-02-25T08:44:14.365430Z","iopub.execute_input":"2023-02-25T08:44:14.366037Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Save to CSV","metadata":{}},{"cell_type":"code","source":"df_importances.to_csv(fn4save)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Some analysis","metadata":{}},{"cell_type":"code","source":"df_importances.describe()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"d = df_y[df_importances.columns].describe()\nd = d.T\nd.sort_values('max',ascending=False, key = abs ).head(50)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"d = df_importances.describe()\nd = d.T\nd.sort_values('max',ascending=False, key = abs ).head(50)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"d.sort_values('max',ascending=False, key = abs ).tail(20)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nplt.figure(figsize = (20,5))\nplt.plot(d['max'].values)\nplt.grid()\nplt.show()\n\nplt.figure(figsize = (20,5))\nplt.plot(d['max'])\nplt.grid()\nplt.show()\n\nplt.figure(figsize = (20,5))\nplt.plot(d['max'].sort_values(ascending = False, key = abs ).head(20) )\nplt.grid()\nplt.show()\n\nplt.figure(figsize = (20,5))\nplt.plot(d['max'].sort_values(ascending = False, key = abs ).head(30) )\nplt.grid()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_corr = df_importances.corr()\ndf_corr","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import seaborn as sns\n#sns.clustermap(df_corr)\nsns.clustermap(np.abs(df_corr),cmap='vlag', figsize=(25, 25), xticklabels=df_corr.index, yticklabels=df_corr.index); plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_importances.sort_values(df_importances.columns[0], ascending = False).head(50)","metadata":{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_importances_ranks = (-df_importances).rank()\nrk = df_importances_ranks.median(axis=1)\nrk.name = 'Median Rank'\nrk = rk.sort_values().to_frame()\nr = df_importances_ranks.mean(axis=1)\nr.name = 'Mean Rank'\nrk = rk.join(r)\nr = df_importances_ranks.min(axis=1)\nr.name = 'Min Rank'\nrk = rk.join(r)\nr = df_importances_ranks.max(axis=1)\nr.name = 'Max Rank'\nrk = rk.join(r)\n\nrk.head(100)\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tt = time.time() - t0start\nprint('%.1f seconds ( = %.1f minutes, = %.2f hours) passed'%( tt, tt/60, tt/3600 ) )","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}