{"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 educational notebook how to use Pytorch Embedding for categorical features and train neural networks on categorical features.\n\nWe use Kaggle \"Open Problems – Single-Cell Perturbations\" as an example . \n\nThe task have only two features - cell type and compound (drug). Both features are categorical. With 6 and 146 unique categiries respectively.\n\nThere are 18211 targets, it is quite a lot, so typical way is: to reduce targets by pca-like transform, predict transformed targets, and then use reducer.inverse_transform to get original targets. It is not only speed up the computation, but often gives better quality since it is kind of denoising. We reduce by tsvd to 35 dimensions. And so neural net is predicting only 35 outputs using 2 categorical features in process learning their embeddings. \n\nSo we use train simple neural network using these two categorical features and 35 targets. \nWe do not fix random seed so results are stochastic, but after some trials we to get test r2-score positive - so model is not so bad locally, and submission with LB 0.635 - confirms that it is quite comparable with other simple baselines.\n\nThe neural net is very simple  - embeddings and one non-linear hidden layer. Sizes for embeddings 2,5, hidden layer 10 - some randomly chosen.  Param tuning have NOT beed done - you are welcome to do it:\n\n\n### Educational task: change params for neural network to get better results.\n\n\nPS\n\nSome EDA and baselines can be found here: https://www.kaggle.com/code/alexandervc/op2-eda-baseline-s\n\nSome other educational notebooks on Pytorch and Keras are here : https://www.kaggle.com/datasets/alexandervc/nns-pytorch-keras-etc,\nand some other examples are: https://www.kaggle.com/code/alexandervc/pytorch-01-basics , more sophisticated examples: https://www.kaggle.com/code/alexandervc/pytorch-2-cv-focal-sophia-etc, https://www.kaggle.com/code/alexandervc/pytorch-keras-etc-3-blend-cafa-metric-etc\n\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\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-09-25T20:09:59.825395Z","iopub.execute_input":"2023-09-25T20:09:59.826251Z","iopub.status.idle":"2023-09-25T20:10:02.582102Z","shell.execute_reply.started":"2023-09-25T20:09:59.826206Z","shell.execute_reply":"2023-09-25T20:10:02.580391Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load Train data ","metadata":{}},{"cell_type":"code","source":"%%time\nfn = '/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet'\ndf_de_train = pd.read_parquet(fn)# , index_col = 0)\nprint(df_de_train.shape)\ndf_de_train","metadata":{"execution":{"iopub.status.busy":"2023-09-25T20:10:02.584543Z","iopub.execute_input":"2023-09-25T20:10:02.585544Z","iopub.status.idle":"2023-09-25T20:10:06.202321Z","shell.execute_reply.started":"2023-09-25T20:10:02.585495Z","shell.execute_reply":"2023-09-25T20:10:06.201009Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Test data and show sample submission example","metadata":{}},{"cell_type":"code","source":"%%time\nfn = '/kaggle/input/open-problems-single-cell-perturbations/id_map.csv'\ndf_id_map = pd.read_csv(fn)\nprint(df_id_map.shape)\ndisplay(df_id_map)\nfn = '/kaggle/input/open-problems-single-cell-perturbations/sample_submission.csv'\ndf = pd.read_csv(fn, index_col = 0)\nprint(df.shape)\ndf","metadata":{"execution":{"iopub.status.busy":"2023-09-25T20:10:06.203987Z","iopub.execute_input":"2023-09-25T20:10:06.204425Z","iopub.status.idle":"2023-09-25T20:10:11.573180Z","shell.execute_reply.started":"2023-09-25T20:10:06.204386Z","shell.execute_reply":"2023-09-25T20:10:11.571913Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import torch \n\ndef r2_score_torch(y_true, y_pred):\n    \"\"\"\n    Compute the R-squared (R2) score between true and predicted values.\n\n    Args:\n    - y_true (torch.Tensor): True target values (ground truth).\n    - y_pred (torch.Tensor): Predicted target values.\n\n    Returns:\n    - r2 (float): R-squared score.\n    \"\"\"\n    y_mean = torch.mean(y_true)\n    ss_tot = torch.sum((y_true - y_mean)**2)\n    ss_res = torch.sum((y_true - y_pred)**2)\n    r2 = 1 - (ss_res / ss_tot)\n    return r2.item()","metadata":{"execution":{"iopub.status.busy":"2023-09-25T20:10:11.576148Z","iopub.execute_input":"2023-09-25T20:10:11.576552Z","iopub.status.idle":"2023-09-25T20:10:15.330212Z","shell.execute_reply.started":"2023-09-25T20:10:11.576520Z","shell.execute_reply":"2023-09-25T20:10:15.328660Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train test split ","metadata":{}},{"cell_type":"code","source":"from sklearn.model_selection import train_test_split\nIX_train, IX_test  = train_test_split(range(len(df_de_train)), test_size=0.1, random_state=42,shuffle = True)\nprint(IX_train[:10], IX_test[:10])","metadata":{"execution":{"iopub.status.busy":"2023-09-25T20:10:15.332605Z","iopub.execute_input":"2023-09-25T20:10:15.333533Z","iopub.status.idle":"2023-09-25T20:10:15.571687Z","shell.execute_reply.started":"2023-09-25T20:10:15.333471Z","shell.execute_reply":"2023-09-25T20:10:15.570619Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Preparation - dimensional reduction of targets \n\nWe have 18211 targets, it quite a lot. Typicaly way in such situations - reduce dimensions by pca/ica/tsvd/... predict these transformed targets and then use reducer.inverse_transform() method to return to the original space. It is not only faster, but often gives better score because such procedure works as a kind of denoising. \n\n\nWe use tsvd and 35 dimensions since some previous experiments suggest that way works better. \nSee https://www.kaggle.com/code/alexandervc/op2-eda-baseline-s ,  and also may look at:   https://www.kaggle.com/code/alexandervc/op2-models-cv-tuning\n","metadata":{}},{"cell_type":"code","source":"%%time \nfrom sklearn.decomposition import TruncatedSVD\n\nn_components = 35\n\nreducer = TruncatedSVD(n_components=n_components, n_iter=7, random_state=42)\nprint(reducer)\n\nY = df_de_train.iloc[:,5:].values\nYr = reducer.fit_transform(Y)\nprint(Y.shape, Yr.shape)","metadata":{"execution":{"iopub.status.busy":"2023-09-25T20:10:15.573677Z","iopub.execute_input":"2023-09-25T20:10:15.574214Z","iopub.status.idle":"2023-09-25T20:10:18.694434Z","shell.execute_reply.started":"2023-09-25T20:10:15.574164Z","shell.execute_reply":"2023-09-25T20:10:18.692999Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Preparation - categorical to ordinal encoding - necessary for pytorch models\n\nPytorch does not accept categorical features directly so we first encode categories by numbers - e.g. sklearn ordginal encoder suits.\n\nThat is technical step\n","metadata":{}},{"cell_type":"code","source":"from sklearn.preprocessing import OrdinalEncoder\nenc = OrdinalEncoder()\nX = enc.fit_transform(df_de_train[ ['cell_type','sm_name']])\nprint( len(np.unique(X[:,0])),  len(np.unique(X[:,1]) )) \nprint( str(enc.categories_)[:30],'   ', str(enc.categories_)[-30:]  )\nprint('X.shape',X.shape)\nprint(X[:3,:2])\nprint()\nX_submit = enc.transform( df_id_map[ ['cell_type','sm_name']]  )\nprint('X_submit.shape', X_submit.shape)\nprint(X_submit[:3,:2])\n\n","metadata":{"execution":{"iopub.status.busy":"2023-09-25T20:10:18.696984Z","iopub.execute_input":"2023-09-25T20:10:18.697416Z","iopub.status.idle":"2023-09-25T20:10:18.721371Z","shell.execute_reply.started":"2023-09-25T20:10:18.697377Z","shell.execute_reply":"2023-09-25T20:10:18.719961Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Prepare/format train and test data\n\nPay attention on types \"long\" (and also better \"float\", sometimes \"double\" may cause errors). ","metadata":{}},{"cell_type":"code","source":"%%time\nimport torch\n\nX_train = torch.tensor( X[IX_train,:], dtype = torch.long  ) \nX_test = torch.tensor( X[IX_test,:], dtype = torch.long )\nX_submit = torch.tensor( X_submit, dtype = torch.long )\nprint(X_train.shape, X_test.shape, X_submit.shape)\n\nYr_train = torch.tensor( Yr[IX_train,:], dtype = torch.float  ) \nYr_test = torch.tensor( Yr[IX_test,:], dtype = torch.float  ) \nprint(Yr_train.shape, Yr_test.shape ) \n","metadata":{"execution":{"iopub.status.busy":"2023-09-25T20:10:18.723781Z","iopub.execute_input":"2023-09-25T20:10:18.724378Z","iopub.status.idle":"2023-09-25T20:10:18.776351Z","shell.execute_reply.started":"2023-09-25T20:10:18.724323Z","shell.execute_reply":"2023-09-25T20:10:18.774923Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# The model","metadata":{}},{"cell_type":"code","source":"%%time\n\nimport torch\nimport torch.nn as nn\nimport torch.optim as optim\n\n\n# Define the neural network model with two embedding layers\nclass EmbNN(nn.Module):\n    def __init__(self, num_categories1=6, num_categories2=146, emb1_size=2, emb2_size=5, hidden_layer1 = 10 ):\n        super(EmbNN, self).__init__()\n        self.embedding1 = nn.Embedding(num_categories1, emb1_size)  # Embedding layer for categorical feature 1\n        self.embedding2 = nn.Embedding(num_categories2, emb2_size)   # Embedding layer for categorical feature 2\n        self.fc1 = nn.Linear(emb1_size + emb2_size , hidden_layer1)  # Input: 10 (from embedding 1) + 5 (from embedding 2) + 1 (continuous), Output: 10 features\n        self.fc2 = nn.Linear(hidden_layer1, Yr.shape[1])  # Input: 10 features, Output: 1 feature\n\n\n    def forward(self, cat_input1, cat_input2):\n        cat_embed1 = self.embedding1(cat_input1.long() )  # Embed categorical feature 1\n        cat_embed2 = self.embedding2(cat_input2.long())  # Embed categorical feature 2\n        x = torch.cat((cat_embed1.view(cat_input1.size(0), -1), cat_embed2.view(cat_input2.size(0), -1)), dim=1)  # Concatenate with continuous input\n        x = torch.relu(self.fc1(x))\n        x = self.fc2(x)\n        return x\n\n    def emb1(self, cat_input):\n        cat_embed1 = self.embedding1(cat_input.long() )  # Embed categorical feature 1\n        return cat_embed1\n    def emb2(self, cat_input):\n        cat_embed2 = self.embedding2(cat_input.long() )  # Embed categorical feature 1\n        return cat_embed2\n    ","metadata":{"execution":{"iopub.status.busy":"2023-09-25T20:10:18.778071Z","iopub.execute_input":"2023-09-25T20:10:18.778483Z","iopub.status.idle":"2023-09-25T20:10:18.793643Z","shell.execute_reply.started":"2023-09-25T20:10:18.778445Z","shell.execute_reply":"2023-09-25T20:10:18.792196Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create the neural network example \nnum_categories1 = len( set(X[:,0]) )\nnum_categories2 = len( set(X[:,1]) )\nemb1_size = 2\nemb2_size = 5\nhidden_layer1 = 5 \nmodel = EmbNN(num_categories1=num_categories1, num_categories2=num_categories2,  emb1_size=emb1_size, emb2_size=emb2_size, hidden_layer1 = hidden_layer1) \n\nprint(model)\n","metadata":{"execution":{"iopub.status.busy":"2023-09-25T20:10:18.798387Z","iopub.execute_input":"2023-09-25T20:10:18.798815Z","iopub.status.idle":"2023-09-25T20:10:18.829471Z","shell.execute_reply.started":"2023-09-25T20:10:18.798771Z","shell.execute_reply":"2023-09-25T20:10:18.828147Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train the model and check your luck :) and welcome to improve\n\nWe do not fix random seed and thus results would always be different. \n\nIt is too simple neural net and sometimes it returns bad scores i.e. r2 can be negative on test.\n\nYou are welcome to improve it or check your luck )\n\n","metadata":{}},{"cell_type":"code","source":"%%time \nimport time \nt0 = time.time()\n\n# Create the neural network\nnum_categories1 = len( set(X[:,0]) )\nnum_categories2 = len( set(X[:,1]) )\nemb1_size = 2\nemb2_size = 5\nhidden_layer1 = 5 \nmodel = EmbNN(num_categories1=num_categories1, num_categories2=num_categories2,  emb1_size=emb1_size, emb2_size=emb2_size, hidden_layer1 = hidden_layer1) \n\nprint(model)\n\n# Define loss function and optimizer\ncriterion = nn.MSELoss()  # Mean Squared Error loss\noptimizer = optim.SGD(model.parameters(), lr=0.01)  # Stochastic Gradient Descent\n\n# Training loop\nnum_epochs = 5000\n\nfor epoch in range(num_epochs):\n    model.train()\n    # Forward pass\n#     print(epoch)\n    categorical_feature1 = X_train[:,0]\n    categorical_feature2 = X_train[:,1]\n    preds = model(categorical_feature1, categorical_feature2)\n#     print(outputs.shape)\n    loss = criterion(preds, Yr_train)# .float() )\n\n    # Backward pass and optimization\n    optimizer.zero_grad()\n    loss.backward()\n    optimizer.step()\n\n    if (epoch + 1) % 500 == 0:\n        print(f'Epoch [{epoch + 1}/{num_epochs}], Loss: {loss.item():.4f}')\n\n\n        model.eval()\n        with torch.no_grad():       \n            categorical_feature1 = X_test[:,0]\n            categorical_feature2 = X_test[:,1]\n            preds = model(categorical_feature1, categorical_feature2)\n            loss_test = criterion(preds, Yr_test)# .float() )\n            print('loss_test:',loss_test.item())\n            s = r2_score_torch(Yr_test, preds )\n            print('r2 test:',s)\n            \n            \n            if s > 0.4: break\n            ","metadata":{"execution":{"iopub.status.busy":"2023-09-25T20:17:17.647351Z","iopub.execute_input":"2023-09-25T20:17:17.648528Z","iopub.status.idle":"2023-09-25T20:17:19.306734Z","shell.execute_reply.started":"2023-09-25T20:17:17.648473Z","shell.execute_reply":"2023-09-25T20:17:19.304451Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fn = 'model_stat_dict'\ntorch.save(model.state_dict(), fn)","metadata":{"execution":{"iopub.status.busy":"2023-09-25T20:17:24.732511Z","iopub.execute_input":"2023-09-25T20:17:24.733227Z","iopub.status.idle":"2023-09-25T20:17:24.739806Z","shell.execute_reply.started":"2023-09-25T20:17:24.733186Z","shell.execute_reply":"2023-09-25T20:17:24.738907Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Visualization of the obtained embeddings ","metadata":{}},{"cell_type":"code","source":"%%time \nX_full = torch.tensor( X, dtype = torch.long  ) \nmodel.eval()\nwith torch.no_grad():       \n    embs1 = model.emb1(X_full[:,0])\n    embs2 = model.emb2(X_full[:,1])\nd1 = pd.DataFrame(embs1)\nd1['cell_type'] = df_de_train['cell_type']\nd1 = d1.drop_duplicates()\nprint(d1.shape)\nsns.scatterplot(x = d1[0], y = d1[1], hue = d1['cell_type'])\nplt.title('Embeddings cell type',fontsize = 20)\nplt.grid()\nplt.show()\nd1.to_csv('embeds_cell_type.csv')\n\nd2 = pd.DataFrame(embs2)\n#d2['compound'] = df_de_train['sm_name']\nd2 = d2.drop_duplicates()\nprint(d2.shape)\nsns.pairplot(d2)\nplt.suptitle('Embeddings compound',fontsize = 20)\n# sns.scatterplot(x = d2[0], y = d2[1])# , hue = d2['compound'])\nplt.grid()\nplt.show()\nd2 = pd.DataFrame(embs2)\nd2['compound'] = df_de_train['sm_name']\nd2 = d2.drop_duplicates()\nprint(d2.shape)\nd2.to_csv('embeds_compound.csv')\n","metadata":{"execution":{"iopub.status.busy":"2023-09-25T20:17:27.201290Z","iopub.execute_input":"2023-09-25T20:17:27.201691Z","iopub.status.idle":"2023-09-25T20:17:36.786379Z","shell.execute_reply.started":"2023-09-25T20:17:27.201659Z","shell.execute_reply":"2023-09-25T20:17:36.785003Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Prepare submission","metadata":{}},{"cell_type":"code","source":"model.eval()\nwith torch.no_grad():       \n    categorical_feature1 = X_submit[:,0]\n    categorical_feature2 = X_submit[:,1]\n    Y_reduced_preds = model(categorical_feature1, categorical_feature2)\n\nprint(Y_reduced_preds.shape)    \nY_submit = reducer.inverse_transform( Y_reduced_preds )\nprint(Y_submit.shape)","metadata":{"execution":{"iopub.status.busy":"2023-09-25T20:17:36.788504Z","iopub.execute_input":"2023-09-25T20:17:36.788921Z","iopub.status.idle":"2023-09-25T20:17:36.844524Z","shell.execute_reply.started":"2023-09-25T20:17:36.788875Z","shell.execute_reply":"2023-09-25T20:17:36.842332Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndf_submit = pd.DataFrame(Y_submit, columns = df_de_train.columns[5:])\ndf_submit.index.name = 'id'\nprint( df_submit.shape )\ndisplay(df_submit)\ndf_submit.to_csv('submission.csv')","metadata":{"execution":{"iopub.status.busy":"2023-09-25T20:17:38.809371Z","iopub.execute_input":"2023-09-25T20:17:38.809849Z","iopub.status.idle":"2023-09-25T20:17:52.609923Z","shell.execute_reply.started":"2023-09-25T20:17:38.809815Z","shell.execute_reply":"2023-09-25T20:17:52.608964Z"},"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":"print('%.1f seconds passed total '%(time.time()-t0start) )\nprint('%.1f minutes passed total '%( (time.time()-t0start)/60)  )\nprint('%.2f hours passed total '%( (time.time()-t0start)/3600)  )","metadata":{"execution":{"iopub.status.busy":"2023-09-25T20:17:55.730159Z","iopub.execute_input":"2023-09-25T20:17:55.730617Z","iopub.status.idle":"2023-09-25T20:17:55.739355Z","shell.execute_reply.started":"2023-09-25T20:17:55.730578Z","shell.execute_reply":"2023-09-25T20:17:55.737702Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}