{"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\nRidge models\n\n#### Findings:\n\nFor CD36 alpha Ridge = 100_000 . Calculations is quite slow:  # Wall time: 1h 10min 36s for alpha in [1e-3, 1e-2, 1e-1, 1,1e1,1e2,1e3,1e4,1e5,1e6,1e7]  \n\n\n#### Versions: \n\n    4, 5, 6 - n_trials =  100 \n    7, 8 - n_trials = 2, 10 \n    9 - quick draft save \n    10 -  n_trials = 10 , train_size = 0.1\n\n","metadata":{}},{"cell_type":"markdown","source":"# Key Params","metadata":{}},{"cell_type":"code","source":"target_name = 'CD36'\n\ntrain_size = 0.1\n\nn_trials = 10\n\n#n_top = 100 ","metadata":{"execution":{"iopub.status.busy":"2023-01-09T10:15:48.325079Z","iopub.execute_input":"2023-01-09T10:15:48.325543Z","iopub.status.idle":"2023-01-09T10:15:48.355991Z","shell.execute_reply.started":"2023-01-09T10:15:48.325447Z","shell.execute_reply":"2023-01-09T10:15:48.354665Z"},"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-09T10:15:48.358707Z","iopub.execute_input":"2023-01-09T10:15:48.359149Z","iopub.status.idle":"2023-01-09T10:15:48.392581Z","shell.execute_reply.started":"2023-01-09T10:15:48.359113Z","shell.execute_reply":"2023-01-09T10:15:48.391672Z"},"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-09T10:15:48.393972Z","iopub.execute_input":"2023-01-09T10:15:48.394740Z","iopub.status.idle":"2023-01-09T10:15:49.554156Z","shell.execute_reply.started":"2023-01-09T10:15:48.394704Z","shell.execute_reply":"2023-01-09T10:15:49.552746Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load Data","metadata":{}},{"cell_type":"code","source":"%%time\ndf_rna = pd.read_hdf('/kaggle/input/open-problems-multimodal/train_cite_inputs.h5')\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-09T10:15:49.556413Z","iopub.execute_input":"2023-01-09T10:15:49.557555Z"},"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.preprocessing import StandardScaler\nscaler = StandardScaler()\n\nX = scaler.fit_transform( X )\nprint(X.shape)\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"##  Determine optimal Alpha","metadata":{}},{"cell_type":"code","source":"%%time\nfrom sklearn.linear_model import LassoCV\nfrom sklearn.linear_model import RidgeCV\nfrom sklearn.linear_model import Ridge\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\n\n#alpha_selected = 1e4\n\nn_splits_for_cross_valdition = 2\nrandom_state_cross_validation = 0\nkf = 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)\nmodel = 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\nprint(model)\nalpha_selected = model.alpha_\nprint( alpha_selected )\nprint()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"alpha_selected ","metadata":{"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\nstr_method = 'Ridge '\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.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    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    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))\ns = df_importances.median(axis = 1) \ns.sort_values(key = abs, ascending = False)\n    ","metadata":{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Analysis of the obtained importances ","metadata":{}},{"cell_type":"code","source":"df_importances.describe()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_importances.corr()","metadata":{"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":{"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":{"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":{"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,600])\n    plt.xlim([0,5000])\n    plt.legend(fontsize = 20 )\n    plt.grid()\n    plt.show()","metadata":{"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":{"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":{"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":{"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":{"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":"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"str_method, n_trials","metadata":{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}