{"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":"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-10-04T18:32:18.258704Z","iopub.execute_input":"2023-10-04T18:32:18.259108Z","iopub.status.idle":"2023-10-04T18:32:18.269608Z","shell.execute_reply.started":"2023-10-04T18:32:18.259077Z","shell.execute_reply":"2023-10-04T18:32:18.268360Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Problems encontered\n\nWe only have two features, cell type and substance. Given that we want a multioutput, that's really complex.\n\n- For this we can use dimensionality reduction, PCA, LDA or anyother technic. There are some technics that use biological aware dimensionality reduction, like the one used in [OP2:🧠 Biologically-aware dimensionality reduction](https://www.kaggle.com/code/pablormier/op2-biologically-aware-dimensionality-reduction/notebook) [FunctionTransformer](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.FunctionTransformer.html)\n\n- Also another problem is that the features is categorical, we need to do some encoding to use these data in our model. We can use the One-hot Encoding(pd.get_dummies) or Target encoding, due to the nominal nature of our variables.\n\n- We need to consider the judge prizes.\n\n## TO DO's\n\n- Point 2 and 3 below\n- Join id Map","metadata":{"execution":{"iopub.status.busy":"2023-09-28T14:05:08.137514Z","iopub.execute_input":"2023-09-28T14:05:08.138031Z","iopub.status.idle":"2023-09-28T14:05:08.177880Z","shell.execute_reply.started":"2023-09-28T14:05:08.137987Z","shell.execute_reply":"2023-09-28T14:05:08.176426Z"}}},{"cell_type":"markdown","source":"## First notebook: Training models\n\n1. We need to make the custom error function:\n$$\nMRRMSE = \\frac{1}{R}\\sum^{R}_{i=1}\\left(\\frac{1}{n}\\sum_{j=1}^{n}(y_{ij}-\\hat{y_{ij}})^2 \\right)^{\\frac{1}{2}}\n$$\n\n2. We need to train a Baseline Model to see our base results\n3. Than we need to test the models","metadata":{}},{"cell_type":"markdown","source":"## Making the error function\n\nWe can adapt the [mean_squared_error](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.mean_squared_error.html) from sklearn, the only change we need to make is that we want the row wise error instead of the column wise.\n\nGIven that the can extend the concept to others metrics","metadata":{}},{"cell_type":"code","source":"from sklearn.metrics import mean_squared_error\ndef mrrmse(y, y_preds):\n    return mean_squared_error(y_preds.T, y.T, squared=False) ","metadata":{"execution":{"iopub.status.busy":"2023-10-04T18:32:19.056293Z","iopub.execute_input":"2023-10-04T18:32:19.056695Z","iopub.status.idle":"2023-10-04T18:32:19.061340Z","shell.execute_reply.started":"2023-10-04T18:32:19.056663Z","shell.execute_reply":"2023-10-04T18:32:19.060593Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Making a Baseline Model\n\nThe Baseline Model wil return the average for each gene and calculate it error in a kfold.\n\nBut first we need to import and transform the data.\n1. Import the training data, ```de_train.parquet```, features are the ```cell_type``` and ```sm_name```, and our targets are from the fifth column foward.\n2. We map the id from each observation, obtained from ```id_map.csv``` **TO DO's**","metadata":{}},{"cell_type":"code","source":"de_train_df = pd.read_parquet(\"/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet\")\nde_train_df","metadata":{"execution":{"iopub.status.busy":"2023-10-04T18:32:19.455769Z","iopub.execute_input":"2023-10-04T18:32:19.456169Z","iopub.status.idle":"2023-10-04T18:32:23.436474Z","shell.execute_reply.started":"2023-10-04T18:32:19.456134Z","shell.execute_reply":"2023-10-04T18:32:23.435347Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"de_train_df.drop([\"cell_type\", \"sm_lincs_id\", \"SMILES\"], axis=1).groupby(\"sm_name\").mean()","metadata":{"execution":{"iopub.status.busy":"2023-10-04T18:32:23.438777Z","iopub.execute_input":"2023-10-04T18:32:23.439194Z","iopub.status.idle":"2023-10-04T18:32:23.843368Z","shell.execute_reply.started":"2023-10-04T18:32:23.439155Z","shell.execute_reply":"2023-10-04T18:32:23.842148Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# TO DO\n##id_map_df = pd.read_csv(\"/kaggle/input/open-problems-single-cell-perturbations/id_map.csv\")\n##id_map_df","metadata":{"execution":{"iopub.status.busy":"2023-10-04T18:32:23.845366Z","iopub.execute_input":"2023-10-04T18:32:23.845820Z","iopub.status.idle":"2023-10-04T18:32:23.850092Z","shell.execute_reply.started":"2023-10-04T18:32:23.845779Z","shell.execute_reply":"2023-10-04T18:32:23.849092Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Baseline Model and getting the base values","metadata":{}},{"cell_type":"code","source":"class Baseline:\n    \n    def __init__(self, type = \"mean\"):\n        self.type = type\n    \n    def fit(self, y: pd.DataFrame, X: pd.DataFrame):\n        if self.type == \"mean\":\n            self.means = y.mean()\n        elif self.type == \"median\":\n            self.means = y.median()\n        elif self.type == \"conditional_mean\":\n            # Mean by sm type\n            df_tmp = pd.concat([X, y],  axis=1)\n            means = df_tmp.iloc[:, 1:].groupby(\"sm_name\").mean()\n            self.means = means\n    def predict(self, X):\n        if self.type == \"conditional_mean\":\n            df = pd.DataFrame()\n            for sm in X[\"sm_name\"].values:\n                try:\n                    df = pd.concat([df, self.means.loc[sm]], axis=1)\n                except:\n                    df = pd.concat([df, self.means.mean()], axis=1)\n            df = df.T\n        else:\n            df = pd.DataFrame([self.means]*X.shape[0])\n        return df","metadata":{"execution":{"iopub.status.busy":"2023-10-04T18:32:23.853434Z","iopub.execute_input":"2023-10-04T18:32:23.853803Z","iopub.status.idle":"2023-10-04T18:32:23.867760Z","shell.execute_reply.started":"2023-10-04T18:32:23.853773Z","shell.execute_reply":"2023-10-04T18:32:23.866892Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.model_selection import KFold\nfrom sklearn.metrics import r2_score","metadata":{"execution":{"iopub.status.busy":"2023-10-04T18:32:23.869190Z","iopub.execute_input":"2023-10-04T18:32:23.869837Z","iopub.status.idle":"2023-10-04T18:32:23.884867Z","shell.execute_reply.started":"2023-10-04T18:32:23.869806Z","shell.execute_reply":"2023-10-04T18:32:23.883199Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"kfold = KFold(shuffle=True)\nfeatures = [\"cell_type\", \"sm_name\"]\n\nfor i, (train, test) in enumerate(kfold.split(de_train_df)):\n    model = Baseline(\"mean\")\n    X_train = de_train_df.loc[train, features]\n    y_train = de_train_df.iloc[train, 5:]\n    X_test = de_train_df.loc[test, features]\n    y_test = de_train_df.iloc[test, 5:]\n    model.fit(y_train, X_train)\n    y_preds = model.predict(X_test)\n    print(f\"fold {i}, r2 score = {r2_score(y_test.T, y_preds.T)}, MRRMSE = {mrrmse(y_test, y_preds)}\")","metadata":{"execution":{"iopub.status.busy":"2023-10-04T18:32:23.886805Z","iopub.execute_input":"2023-10-04T18:32:23.887686Z","iopub.status.idle":"2023-10-04T18:32:25.588346Z","shell.execute_reply.started":"2023-10-04T18:32:23.887645Z","shell.execute_reply":"2023-10-04T18:32:25.587334Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i, (train, test) in enumerate(kfold.split(de_train_df)):\n    model = Baseline(\"median\")\n    X_train = de_train_df.loc[train, features]\n    y_train = de_train_df.iloc[train, 5:]\n    X_test = de_train_df.loc[test, features]\n    y_test = de_train_df.iloc[test, 5:]\n    model.fit(y_train, X_train)\n    y_preds = model.predict(X_test)\n    print(f\"fold {i}, r2 score = {r2_score(y_test.T, y_preds.T)}, MRRMSE = {mrrmse(y_test, y_preds)}\")","metadata":{"execution":{"iopub.status.busy":"2023-10-04T18:32:25.592158Z","iopub.execute_input":"2023-10-04T18:32:25.592566Z","iopub.status.idle":"2023-10-04T18:32:30.367798Z","shell.execute_reply.started":"2023-10-04T18:32:25.592532Z","shell.execute_reply":"2023-10-04T18:32:30.366500Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i, (train, test) in enumerate(kfold.split(de_train_df)):\n    model = Baseline(\"conditional_mean\")\n    X_train = de_train_df.loc[train, features]\n    y_train = de_train_df.iloc[train, 5:]\n    X_test = de_train_df.loc[test, features]\n    y_test = de_train_df.iloc[test, 5:]\n    model.fit(y_train, X_train)\n    y_preds = model.predict(X_test)\n    print(f\"fold {i}, r2 score = {r2_score(y_test.T, y_preds.T)}, MRRMSE = {mrrmse(y_test, y_preds)}\")","metadata":{"execution":{"iopub.status.busy":"2023-10-04T18:32:30.369248Z","iopub.execute_input":"2023-10-04T18:32:30.369941Z","iopub.status.idle":"2023-10-04T18:32:35.363435Z","shell.execute_reply.started":"2023-10-04T18:32:30.369902Z","shell.execute_reply":"2023-10-04T18:32:35.362240Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Training our first model\n\nFor our first model we will do a simple Random Forest, without hyperparameter tunning, and a XGBoost, also without hyperparameter tunnning, only to have some baseline models and in the future we will dig deeper to improve the model performance.","metadata":{}},{"cell_type":"code","source":"from sklearn.ensemble import RandomForestRegressor\nfrom xgboost import XGBRegressor\nfrom sklearn.multioutput import MultiOutputRegressor\nfrom sklearn.preprocessing import OneHotEncoder\nfrom sklearn.decomposition import PCA\nimport warnings","metadata":{"execution":{"iopub.status.busy":"2023-10-04T18:34:44.364089Z","iopub.execute_input":"2023-10-04T18:34:44.364581Z","iopub.status.idle":"2023-10-04T18:34:44.370098Z","shell.execute_reply.started":"2023-10-04T18:34:44.364544Z","shell.execute_reply":"2023-10-04T18:34:44.368965Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# We need to encode our features\nid_map_df = pd.read_csv(\"/kaggle/input/open-problems-single-cell-perturbations/id_map.csv\")\none_hot = OneHotEncoder(sparse_output=False)","metadata":{"execution":{"iopub.status.busy":"2023-10-04T18:32:35.372071Z","iopub.execute_input":"2023-10-04T18:32:35.372617Z","iopub.status.idle":"2023-10-04T18:32:35.454943Z","shell.execute_reply.started":"2023-10-04T18:32:35.372586Z","shell.execute_reply":"2023-10-04T18:32:35.453227Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with warnings.catch_warnings():\n    warnings.simplefilter(\"ignore\")\n    for i, (train, test) in enumerate(kfold.split(de_train_df)):\n        oh = OneHotEncoder(sparse_output=False)\n        df_train = de_train_df.loc[train]\n        y_train = df_train.iloc[:, 5:]\n        pca = PCA()\n        y_train = pca.fit_transform(y_train)\n        X_train_oh = oh.fit_transform(df_train[features])\n        df_test = de_train_df.loc[test]\n        y_test = df_test.iloc[:, 5:]\n        drop_test = [j for j, f in enumerate(df_test[features].values.flatten(\"F\")) if (f not in df_train[features[0]].values) and (f not in df_train[features[1]].values)]\n        drop_test = [j if j < df_test.shape[0] else j - df_test.shape[0] for j in drop_test]\n        df_test = df_test.drop(df_test.index[drop_test])\n        X_test_oh = oh.transform(df_test[features])\n        y_test = df_test.iloc[:, 5:]\n        model = MultiOutputRegressor(RandomForestRegressor(n_estimators = 20, max_depth=3))\n        model.fit(X_train_oh, y_train)\n        y_preds = pca.inverse_transform(model.predict(X_test_oh))\n        print(f\"fold {i}, r2 score = {r2_score(y_test.T, y_preds.T)}, MRRMSE = {mrrmse(y_test, y_preds)}\")   ","metadata":{"execution":{"iopub.status.busy":"2023-10-04T18:36:38.519232Z","iopub.execute_input":"2023-10-04T18:36:38.519622Z","iopub.status.idle":"2023-10-04T18:38:26.021102Z","shell.execute_reply.started":"2023-10-04T18:36:38.519592Z","shell.execute_reply":"2023-10-04T18:38:26.019488Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with warnings.catch_warnings():\n    warnings.simplefilter(\"ignore\")\n    for i, (train, test) in enumerate(kfold.split(de_train_df)):\n        oh = OneHotEncoder(sparse_output=False)\n        df_train = de_train_df.loc[train]\n        y_train = df_train.iloc[:, 5:]\n        pca = PCA()\n        y_train = pca.fit_transform(y_train)\n        X_train_oh = oh.fit_transform(df_train[features])\n        df_test = de_train_df.loc[test]\n        y_test = df_test.iloc[:, 5:]\n        drop_test = [j for j, f in enumerate(df_test[features].values.flatten(\"F\")) if (f not in df_train[features[0]].values) and (f not in df_train[features[1]].values)]\n        drop_test = [j if j < df_test.shape[0] else j - df_test.shape[0] for j in drop_test]\n        df_test = df_test.drop(df_test.index[drop_test])\n        X_test_oh = oh.transform(df_test[features])\n        y_test = df_test.iloc[:, 5:]\n        model = MultiOutputRegressor(XGBRegressor(n_estimators = 20))\n        model.fit(X_train_oh, y_train)\n        y_preds = pca.inverse_transform(model.predict(X_test_oh))\n        print(f\"fold {i}, r2 score = {r2_score(y_test.T, y_preds.T)}, MRRMSE = {mrrmse(y_test, y_preds)}\")   ","metadata":{"execution":{"iopub.status.busy":"2023-10-04T18:41:02.383652Z","iopub.execute_input":"2023-10-04T18:41:02.384015Z","iopub.status.idle":"2023-10-04T18:43:13.772107Z","shell.execute_reply.started":"2023-10-04T18:41:02.383987Z","shell.execute_reply":"2023-10-04T18:43:13.770614Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_train = de_train_df[features]\ny_train = de_train_df.iloc[:, 5:]\nX_test = id_map_df[features]\noh = OneHotEncoder(sparse_output=False)\npca = PCA()\nX_train_oh = oh.fit_transform(X_train)\nX_test_oh = oh.transform(X_test)\ny_train_pca = pca.fit_transform(y_train)\nvalid_sm = X_train[\"sm_name\"].unique()\nvalid_cell_type = X_train[\"cell_type\"].unique()\nvalid_types = X_test[\"sm_name\"].isin(valid_sm) & X_test[\"cell_type\"].isin(valid_cell_type)\nX_test = X_test.loc[valid_types]\nmodel = MultiOutputRegressor(XGBRegressor(n_estimators = 20))\nmodel.fit(X_train_oh, y_train_pca)\ny_preds =pca.inverse_transform(model.predict(X_test_oh))","metadata":{"execution":{"iopub.status.busy":"2023-10-04T19:00:01.626191Z","iopub.execute_input":"2023-10-04T19:00:01.627492Z","iopub.status.idle":"2023-10-04T19:00:41.813522Z","shell.execute_reply.started":"2023-10-04T19:00:01.627436Z","shell.execute_reply":"2023-10-04T19:00:41.811797Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y_preds_df = pd.DataFrame(y_preds, index=id_map_df[\"id\"], columns = y_train.columns)\ny_preds_df.to_csv(\"submission.csv\")","metadata":{"execution":{"iopub.status.busy":"2023-10-04T19:18:03.270803Z","iopub.execute_input":"2023-10-04T19:18:03.271221Z","iopub.status.idle":"2023-10-04T19:18:12.282791Z","shell.execute_reply.started":"2023-10-04T19:18:03.271193Z","shell.execute_reply":"2023-10-04T19:18:12.281863Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}