{"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\nSimple baselines for Ribonanza Kaggle challenge. \n\n\nVersions:\n\n    5 quick save\n    3,4 LB 0.22934 Onehot + Ridge \n    \n    LB 0.28208 Medians for each position of nucleotide in RNA\n    (EDA notebook v17  https://www.kaggle.com/code/alexandervc/ribonanza-1-eda?scriptVersionId=144103296 )\n    \n    2 LB 0.28255 Medians overs good snr subsamples: DMS:0.123, 2A3:0.214\n    1 LB 0.33975 - all zeros, scoring ~ 5 min, save df submission - 13 mins , load it 1.3min, 269_796_671 - n_rows  \n","metadata":{}},{"cell_type":"markdown","source":"# Preliminaries ","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\ncc = 0\nfor dirname, _, filenames in os.walk('/kaggle/input/stanford-ribonanza-rna-folding-converted/'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\nfor dirname, _, filenames in os.walk('/kaggle/input/stanford-ribonanza-rna-folding-data/'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    if 'stanford-ribonanza-rna-folding-converted' in dirname: continue \n    if 'stanford-ribonanza-rna-folding-data' in dirname: continue \n    for filename in filenames:\n        cc += 1\n        if cc < 20:\n            print(os.path.join(dirname, filename))\n        else:\n            break\n    if cc < 20:\n        pass\n    else:\n        break\n\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-09-30T20:18:29.533569Z","iopub.execute_input":"2023-09-30T20:18:29.533932Z","iopub.status.idle":"2023-09-30T20:18:30.946683Z","shell.execute_reply.started":"2023-09-30T20:18:29.533903Z","shell.execute_reply":"2023-09-30T20:18:30.945648Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load Train sequences\n\n\nselected part with high SNR","metadata":{}},{"cell_type":"code","source":"%%time\nfn = '/kaggle/input/stanford-ribonanza-rna-folding-data/cleaned_163428_rna/train_cleaned_DMS_MaP_163428_info_only.csv'\ndf = pd.read_csv(fn)\n# fn = '/kaggle/input/stanford-ribonanza-rna-folding-converted/train_data.parquet'\n# df = pd.read_parquet(fn)\nprint(df.shape)\ndf","metadata":{"execution":{"iopub.status.busy":"2023-09-30T20:18:32.455065Z","iopub.execute_input":"2023-09-30T20:18:32.455878Z","iopub.status.idle":"2023-09-30T20:18:33.891207Z","shell.execute_reply.started":"2023-09-30T20:18:32.455839Z","shell.execute_reply":"2023-09-30T20:18:33.890167Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndf['len'] = [len(t) for t in df['sequence']]\ndf['len'].value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-09-30T20:18:33.893246Z","iopub.execute_input":"2023-09-30T20:18:33.893640Z","iopub.status.idle":"2023-09-30T20:18:33.971119Z","shell.execute_reply.started":"2023-09-30T20:18:33.893600Z","shell.execute_reply":"2023-09-30T20:18:33.969946Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Onehot encoding of sequences ","metadata":{}},{"cell_type":"code","source":"%%time\nimport time\nt0 = time.time()\n\nN = len(df)\nX = np.zeros( (N, 206*4), dtype = np.int8 )\nprint( X.shape, X.nbytes/1e6 )\nrna_onehot_dict = {'A': [1,0,0,0], 'U': [0,1,0,0], 'G':[0,0,1,0], 'C': [0,0,0,1]}\n\nfor k in range(N):\n    seq = df['sequence'].iat[k]\n    L = len(seq)\n    l =  np.array( [rna_onehot_dict[seq[II]] for II in range(L) ] ).ravel()\n    X[k,:len(l)] = l\n    if (k%10_000 == 0): print(k, '%.1f'%(time.time() - t0 ) )","metadata":{"execution":{"iopub.status.busy":"2023-09-30T20:18:35.570548Z","iopub.execute_input":"2023-09-30T20:18:35.571258Z","iopub.status.idle":"2023-09-30T20:18:55.971917Z","shell.execute_reply.started":"2023-09-30T20:18:35.571227Z","shell.execute_reply":"2023-09-30T20:18:55.971206Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load targets","metadata":{}},{"cell_type":"code","source":"%%time \n# import torch\n\n# fn = '/kaggle/input/ribonanza-data-prepare-y-for-selected/Y_163428_165_2_DMS_2A3_clipped01_torchfloat32.pt'\n# Y = torch.load(fn)\nfn = '/kaggle/input/ribonanza-data-prepare-y-for-selected/Y_163428_165_2_DMS_2A3_clipped01_numpyfloat32.npy'\nY = np.load(fn)\nY.shape, Y[0,40,:]","metadata":{"execution":{"iopub.status.busy":"2023-09-30T20:18:55.973245Z","iopub.execute_input":"2023-09-30T20:18:55.973680Z","iopub.status.idle":"2023-09-30T20:18:57.333833Z","shell.execute_reply.started":"2023-09-30T20:18:55.973656Z","shell.execute_reply":"2023-09-30T20:18:57.333201Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Cut part of data\n\nFor simplicity take 50_000 of samples of len 177","metadata":{}},{"cell_type":"code","source":"%%time\nX = X[13627:(13627+50_000),:]\nY = Y[13627:(13627+50_000),:,:]\n\nprint(X.shape, Y.shape)","metadata":{"execution":{"iopub.status.busy":"2023-09-30T20:18:59.840858Z","iopub.execute_input":"2023-09-30T20:18:59.841700Z","iopub.status.idle":"2023-09-30T20:18:59.847518Z","shell.execute_reply.started":"2023-09-30T20:18:59.841669Z","shell.execute_reply":"2023-09-30T20:18:59.846641Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train valid split","metadata":{}},{"cell_type":"code","source":"N = len(X)\nprint(N)\nnp.random.seed(42)\np = np.random.permutation(N)\np","metadata":{"execution":{"iopub.status.busy":"2023-09-30T20:19:04.627539Z","iopub.execute_input":"2023-09-30T20:19:04.628722Z","iopub.status.idle":"2023-09-30T20:19:04.638507Z","shell.execute_reply.started":"2023-09-30T20:19:04.628681Z","shell.execute_reply":"2023-09-30T20:19:04.637335Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"N_train = 30000\nN_valid = 1000\nN_holdout = N - N_train - N_valid\n\nIX_train = p[:N_train]\nIX_val = p[N_train:(N_train+N_valid)]\nIX_holdout = p[(N_train+N_valid):]\nprint(IX_train.shape, IX_val.shape, IX_holdout.shape )","metadata":{"execution":{"iopub.status.busy":"2023-09-30T20:19:06.496174Z","iopub.execute_input":"2023-09-30T20:19:06.496546Z","iopub.status.idle":"2023-09-30T20:19:06.501772Z","shell.execute_reply.started":"2023-09-30T20:19:06.496521Z","shell.execute_reply":"2023-09-30T20:19:06.500840Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Targets for DMS and 2A3,  padded by zero\n\n\nPadding by zeros is made for simplicity of reformating to submission format by just .ravel() for predictions - technocal things\n\nEssential target values are placed at positions 26-126, and very few non-nans till position 165 - that is stored in file we loaded","metadata":{}},{"cell_type":"code","source":"%%time\nY_DMS = np.zeros( (Y.shape[0], 177 )  )\nY_DMS[:,:165] = Y[:,:,0]\n\nY_2A3 = np.zeros( (Y.shape[0], 177 )  )\nY_2A3[:,:165] = Y[:,:,1]\n\nY_DMS[ np.isnan(Y_DMS)] = 0\nY_2A3[ np.isnan(Y_2A3)] = 0\n\n\nprint( Y_DMS.shape , Y_2A3.shape  )\n","metadata":{"execution":{"iopub.status.busy":"2023-09-30T20:19:08.420202Z","iopub.execute_input":"2023-09-30T20:19:08.420883Z","iopub.status.idle":"2023-09-30T20:19:08.605134Z","shell.execute_reply.started":"2023-09-30T20:19:08.420852Z","shell.execute_reply":"2023-09-30T20:19:08.603984Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_train = X[IX_train,:]\nY_train_DMS = Y_DMS[IX_train,:]\nY_train_2A3 = Y_2A3[IX_train,:]\n\nprint( X_train.shape , Y_train_DMS.shape, Y_train_2A3.shape )\n\nX_val = X[IX_val,:]\nY_val_DMS = Y_DMS[IX_val,:]\nY_val_2A3 = Y_2A3[IX_val,:]\n\nprint( X_val.shape , Y_val_DMS.shape, Y_val_2A3.shape )\n\n","metadata":{"execution":{"iopub.status.busy":"2023-09-30T20:20:20.946537Z","iopub.execute_input":"2023-09-30T20:20:20.946895Z","iopub.status.idle":"2023-09-30T20:20:21.010783Z","shell.execute_reply.started":"2023-09-30T20:20:20.946869Z","shell.execute_reply":"2023-09-30T20:20:21.009893Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Fit Ridge model","metadata":{}},{"cell_type":"code","source":"%%time\nfrom sklearn.linear_model import Ridge\n\nalpha = 1e2\n\nmodel_DMS = Ridge(alpha = alpha )\nmodel_DMS.fit(X_train, Y_train_DMS)\n\nmodel_2A3 = Ridge(alpha = alpha)\nmodel_2A3.fit(X_train, Y_train_2A3)","metadata":{"execution":{"iopub.status.busy":"2023-09-30T20:36:10.012075Z","iopub.execute_input":"2023-09-30T20:36:10.012511Z","iopub.status.idle":"2023-09-30T20:36:11.893212Z","shell.execute_reply.started":"2023-09-30T20:36:10.012483Z","shell.execute_reply":"2023-09-30T20:36:11.891918Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nY_pred_DMS = model_DMS.predict( X_val)\nprint( Y_pred_DMS.shape, Y_val_DMS.shape )\n\nY_pred_2A3 = model_2A3.predict( X_val)\nprint( Y_pred_2A3.shape, Y_val_2A3.shape )\n\nfrom sklearn.metrics import mean_squared_error\nprint('DMS:',  mean_squared_error(Y_val_DMS[:,26:126], Y_pred_DMS[:,26:126]) ) \nprint('2A3:',  mean_squared_error(Y_val_2A3[:,26:126], Y_pred_2A3[:,26:126]) ) \n\n","metadata":{"execution":{"iopub.status.busy":"2023-09-30T20:36:11.896294Z","iopub.execute_input":"2023-09-30T20:36:11.899402Z","iopub.status.idle":"2023-09-30T20:36:11.957987Z","shell.execute_reply.started":"2023-09-30T20:36:11.899349Z","shell.execute_reply":"2023-09-30T20:36:11.956625Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nimport torch \nimport torch.nn as nn\n\n\ncriterion = nn.MSELoss()\n\nN1 = Y_val_DMS.shape[0]\nN2 = Y_val_DMS.shape[1]\nN3 = 2\npreds = torch.zeros( (N1,N2,N3 ) , dtype = torch.float32 )\npreds[:,:,0] = torch.tensor(   Y_pred_DMS )\npreds[:,:,1] = torch.tensor(   Y_pred_2A3 )\n\nY_val = torch.zeros( (N1,N2,N3 ) , dtype = torch.float32 )\nY_val[:,:,0] = torch.tensor(   Y_val_DMS )\nY_val[:,:,1] = torch.tensor(   Y_val_2A3 )\n\n\nloss = criterion(preds[:,26:126,:], Y_val[:,26:126,:])\nprint(f'loss on val: {loss.item():12.5f}' )\n","metadata":{"execution":{"iopub.status.busy":"2023-09-30T20:36:11.964204Z","iopub.execute_input":"2023-09-30T20:36:11.968998Z","iopub.status.idle":"2023-09-30T20:36:11.996493Z","shell.execute_reply.started":"2023-09-30T20:36:11.968944Z","shell.execute_reply":"2023-09-30T20:36:11.995670Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load test sequences data","metadata":{}},{"cell_type":"code","source":"%%time\n# fn = '/kaggle/input/stanford-ribonanza-rna-folding/test_sequences.csv'\nfn = '/kaggle/input/stanford-ribonanza-rna-folding-converted/test_sequences.parquet'\ndft = pd.read_parquet(fn)\nprint(dft.shape)\ndft","metadata":{"execution":{"iopub.status.busy":"2023-09-30T20:36:18.440596Z","iopub.execute_input":"2023-09-30T20:36:18.440946Z","iopub.status.idle":"2023-09-30T20:36:21.009850Z","shell.execute_reply.started":"2023-09-30T20:36:18.440921Z","shell.execute_reply":"2023-09-30T20:36:21.009200Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndft['len'] = [len(t) for t in dft['sequence']]\nprint(dft['len'].iat[335822], dft['len'].iat[335823])\ndft['len'].value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-09-30T20:36:23.203230Z","iopub.execute_input":"2023-09-30T20:36:23.203858Z","iopub.status.idle":"2023-09-30T20:36:23.698370Z","shell.execute_reply.started":"2023-09-30T20:36:23.203830Z","shell.execute_reply":"2023-09-30T20:36:23.697309Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Prepare onehot for test data\n\nTo speed up we will only work with public LB data - i.e. first 335823 sequences of length 177","metadata":{}},{"cell_type":"code","source":"%%time\nimport time\nN = 335823\nL = 206\nX_submit = np.zeros( (N,L*4), np.int8 )\nt0 = time.time()\nfor k in range(N):\n    seq = dft['sequence'].iat[k]\n    L = len(seq)\n    l =  np.array( [rna_onehot_dict[seq[II]] for II in range(L) ] ).ravel()\n    X_submit[k,:len(l)] = l\n    if (k%50_000 == 0): print(k, '%.1f'%(time.time() - t0 ) )\nX_submit        ","metadata":{"execution":{"iopub.status.busy":"2023-09-30T20:36:25.938311Z","iopub.execute_input":"2023-09-30T20:36:25.938658Z","iopub.status.idle":"2023-09-30T20:37:08.357277Z","shell.execute_reply.started":"2023-09-30T20:36:25.938634Z","shell.execute_reply":"2023-09-30T20:37:08.356555Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create df_submit","metadata":{}},{"cell_type":"code","source":"%%time\nN = 269796671\ndf_submit = pd.DataFrame(index = range(N) )\ndf_submit.index.name = 'id'\ndf_submit['reactivity_DMS_MaP'] = np.zeros( N, dtype = np.float16 )\ndf_submit['reactivity_2A3_MaP'] = np.zeros( N, dtype = np.float16 )\nprint(df_submit.values.nbytes/1e6 )\n","metadata":{"execution":{"iopub.status.busy":"2023-09-30T20:22:50.732470Z","iopub.execute_input":"2023-09-30T20:22:50.732736Z","iopub.status.idle":"2023-09-30T20:22:52.418708Z","shell.execute_reply.started":"2023-09-30T20:22:50.732715Z","shell.execute_reply":"2023-09-30T20:22:52.417309Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Make prediction and fill submit","metadata":{}},{"cell_type":"code","source":"%%time\nY_DMS_pred = model_DMS.predict(X_submit)\nY_2A3_pred = model_2A3.predict(X_submit)\nprint(Y_DMS_pred.shape, Y_2A3_pred.shape)\nY_DMS_pred[:5,26:30]","metadata":{"execution":{"iopub.status.busy":"2023-09-30T20:23:44.149638Z","iopub.execute_input":"2023-09-30T20:23:44.150021Z","iopub.status.idle":"2023-09-30T20:23:49.744881Z","shell.execute_reply.started":"2023-09-30T20:23:44.149994Z","shell.execute_reply":"2023-09-30T20:23:49.743898Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Y_2A3_pred[:5,26:30]","metadata":{"execution":{"iopub.status.busy":"2023-09-30T20:23:59.173083Z","iopub.execute_input":"2023-09-30T20:23:59.173425Z","iopub.status.idle":"2023-09-30T20:23:59.180390Z","shell.execute_reply.started":"2023-09-30T20:23:59.173401Z","shell.execute_reply":"2023-09-30T20:23:59.179422Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nNN = Y_DMS_pred.shape[0] * Y_DMS_pred.shape[1] \ndf_submit.loc[:(NN-1),'reactivity_DMS_MaP'] = Y_DMS_pred.ravel()\ndf_submit.loc[:(NN-1),'reactivity_2A3_MaP'] = Y_2A3_pred.ravel()\n# df_submit.loc[:NN, 'reactivity_2A3_MaP'] = Y_2A3_pred.ravel()\n\n","metadata":{"execution":{"iopub.status.busy":"2023-09-30T20:28:36.168102Z","iopub.execute_input":"2023-09-30T20:28:36.168704Z","iopub.status.idle":"2023-09-30T20:28:42.555443Z","shell.execute_reply.started":"2023-09-30T20:28:36.168667Z","shell.execute_reply":"2023-09-30T20:28:42.554385Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndisplay( df_submit.iloc[:NN,:].describe()  ) ","metadata":{"execution":{"iopub.status.busy":"2023-09-30T20:28:51.279041Z","iopub.execute_input":"2023-09-30T20:28:51.279441Z","iopub.status.idle":"2023-09-30T20:28:57.043483Z","shell.execute_reply.started":"2023-09-30T20:28:51.279413Z","shell.execute_reply":"2023-09-30T20:28:57.042445Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_submit.iloc[:NN,:].corr()","metadata":{"execution":{"iopub.status.busy":"2023-09-30T20:28:57.045269Z","iopub.execute_input":"2023-09-30T20:28:57.045558Z","iopub.status.idle":"2023-09-30T20:28:58.594514Z","shell.execute_reply.started":"2023-09-30T20:28:57.045534Z","shell.execute_reply":"2023-09-30T20:28:58.593248Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Save submission\n\n 13-30 minutes","metadata":{}},{"cell_type":"code","source":"%%time\ndf_submit.to_csv('submission.csv')\n","metadata":{"execution":{"iopub.status.busy":"2023-09-24T13:30:46.887638Z","iopub.execute_input":"2023-09-24T13:30:46.888142Z","iopub.status.idle":"2023-09-24T13:51:01.057011Z","shell.execute_reply.started":"2023-09-24T13:30:46.888103Z","shell.execute_reply":"2023-09-24T13:51:01.055131Z"},"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_count":null,"outputs":[]}]}