{"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":"## Here is just a simple demo of how to find the best parameters of KRR by using [Ray Tune](https://docs.ray.io/en/master/tune/index.html). You can use ray turn to look for the best parameter or hyperparameter of almost any model and it can make the best use of the resouces of your computer(both CPU and GPU). You may find more detailed information [here](https://docs.ray.io/en/master/tune/index.html). ","metadata":{}},{"cell_type":"code","source":"from __future__ import print_function\n\nimport argparse\nimport os\nimport torch\nimport torch.optim as optim\n\nimport ray\nfrom ray import air, tune\nfrom ray.tune.schedulers import ASHAScheduler,HyperBandForBOHB\nfrom ray.tune.search.bohb import TuneBOHB\n\nimport time\nfrom tqdm import tqdm\nimport gc\nimport random\nimport glob\nimport copy\nimport pickle\nimport torch\nimport numpy as np\nimport pandas as pd\nfrom torch.utils import tensorboard\nfrom sklearn.model_selection import KFold,GroupKFold\nfrom sklearn.preprocessing import LabelEncoder,StandardScaler,OneHotEncoder\nfrom sklearn.kernel_ridge import KernelRidge\nfrom sklearn.gaussian_process.kernels import RBF","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-11-18T13:43:58.329503Z","iopub.execute_input":"2022-11-18T13:43:58.330008Z","iopub.status.idle":"2022-11-18T13:44:00.142888Z","shell.execute_reply.started":"2022-11-18T13:43:58.329969Z","shell.execute_reply":"2022-11-18T13:44:00.141340Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install -U ray[tune]==2.1.0","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install -U hpbandster ConfigSpace","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Finding the best parameter","metadata":{}},{"cell_type":"code","source":"def correlation_score(y_true, y_pred):\n    \"\"\"Scores the predictions according to the competition rules. \n    \n    It is assumed that the predictions are not constant.\n    \n    Returns the average of each sample's Pearson correlation coefficient\"\"\"\n    if type(y_true) == pd.DataFrame: y_true = y_true.values\n    if type(y_pred) == pd.DataFrame: y_pred = y_pred.values\n    if y_true.shape != y_pred.shape: raise ValueError(\"Shapes are different.\")\n    corrsum = 0\n    for i in range(len(y_true)):\n        corrsum += np.corrcoef(y_true[i], y_pred[i])[1, 0]\n    return corrsum / len(y_true)","metadata":{"execution":{"iopub.status.busy":"2022-11-18T13:44:09.109185Z","iopub.execute_input":"2022-11-18T13:44:09.109694Z","iopub.status.idle":"2022-11-18T13:44:09.118543Z","shell.execute_reply.started":"2022-11-18T13:44:09.109652Z","shell.execute_reply":"2022-11-18T13:44:09.116505Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class TrainKRR(tune.Trainable):\n    def setup(self, config):\n        self.prepare_data()\n        self.config = config      \n\n    def step(self):\n        np.random.seed(42)\n        random.seed(42)\n        kf = GroupKFold(n_splits=3) \n        index = 0\n        mean_score = []\n        for id,(idx_tr, idx_va) in enumerate(kf.split(range(self.train_.shape[0]),groups= self.meta_train.donor)):\n            \n            Xtr, Xva = self.train_[idx_tr], self.train_[idx_va]\n            Ytr, Yva = self.target[idx_tr], self.target[idx_va]\n            random.seed(42)\n            index_to_train = random.choices([i for i in range(Xtr.shape[0])],k = 100) # the great k you choose, the best reults you will get.The training time will also be longer\n            Xtr = Xtr[index_to_train]\n            Ytr = Ytr[index_to_train]\n            # print(f'Fold {id}..')\n        \n            kernel = RBF(self.config[\"gamma\"])\n            model = KernelRidge(kernel=kernel,alpha=self.config[\"alpha\"])\n            model.fit(Xtr,Ytr)\n        \n            del Xtr, Ytr\n            gc.collect()\n            s = correlation_score(Yva, model.predict(Xva))\n            mean_score.append(s)\n            del Xva, Yva\n            gc.collect()\n\n            index += 1\n        mean_score = np.mean(mean_score)\n\n        return {\"mean_score\": mean_score}\n\n    def save_checkpoint(self, checkpoint_dir):\n        checkpoint_path = os.path.join(checkpoint_dir, \"model.pth\")\n        return checkpoint_path\n\n    def load_checkpoint(self, checkpoint_path):\n        pass\n    \n    def prepare_data(self):\n        train = np.load(\"/kaggle/input/cite-final/new_cite_train_final.npz\")[\"arr_0\"]\n        train_index = np.load(\"/kaggle/input/multimodal-single-cell-as-sparse-matrix/train_cite_inputs_idxcol.npz\",allow_pickle=True)\n        \n        meta = pd.read_csv(\"/kaggle/input/open-problems-multimodal/metadata.csv\",index_col = \"cell_id\")\n        meta = meta[meta.technology==\"citeseq\"]\n        lbe = LabelEncoder()\n        meta[\"cell_type\"] = lbe.fit_transform(meta[\"cell_type\"])\n        meta[\"gender\"] = meta.apply(lambda x:0 if x[\"donor\"]==13176 else 1,axis =1)\n        meta_train = meta.reindex(train_index[\"index\"])\n        self.meta_train = meta_train\n        # train_meta = meta_train[\"gender\"].values.reshape(-1, 1)\n        # train = np.concatenate([train,train_meta],axis= -1)\n        train_meta = meta_train[\"cell_type\"].values.reshape(-1, 1)\n        ohe = OneHotEncoder(sparse=False)\n        train_meta = ohe.fit_transform(train_meta)\n        self.train_ = np.concatenate([train,train_meta],axis= -1)\n        \n        target = pd.read_hdf(\"/kaggle/input/open-problems-multimodal/train_cite_targets.h5\").values\n        target -= target.mean(axis=1).reshape(-1, 1)\n        target /= target.std(axis=1).reshape(-1, 1)\n        self.target = target\n    ","metadata":{"execution":{"iopub.status.busy":"2022-11-18T13:45:51.168703Z","iopub.execute_input":"2022-11-18T13:45:51.169182Z","iopub.status.idle":"2022-11-18T13:45:51.191595Z","shell.execute_reply.started":"2022-11-18T13:45:51.169140Z","shell.execute_reply":"2022-11-18T13:45:51.190129Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def main():    \n    ray.init()\n\n    sched = HyperBandForBOHB(metric=\"mean_score\", mode=\"max\")\n    algo = TuneBOHB(metric=\"mean_score\", mode=\"max\")\n\n    tuner = tune.Tuner(\n        tune.with_resources(TrainKRR, resources={\"cpu\": 1, \"gpu\":0 }),\n        run_config=air.RunConfig(\n            stop={\n                \"mean_score\": 0.90,\n            },\n            checkpoint_config=air.CheckpointConfig(\n                checkpoint_at_end=False, checkpoint_frequency=0\n            ),\n            failure_config = air.FailureConfig(\n                max_failures = 0,\n            )\n        ),\n\n        tune_config=tune.TuneConfig(\n            scheduler=sched,\n            num_samples=5,    # the great num_samples you choose, the best reults you will get. The training time will also be longer\n            search_alg = algo\n        ),\n        \n        param_space=dict(\n\n            alpha = tune.quniform(0, 2, 0.1),\n            gamma = tune.randint(9, 21)\n            \n        ),\n    )\n    results = tuner.fit()\n\n#     print(\"Best config is:\", results.get_best_result().config)","metadata":{"execution":{"iopub.status.busy":"2022-11-18T13:47:32.882628Z","iopub.execute_input":"2022-11-18T13:47:32.883136Z","iopub.status.idle":"2022-11-18T13:47:32.895992Z","shell.execute_reply.started":"2022-11-18T13:47:32.883095Z","shell.execute_reply":"2022-11-18T13:47:32.894262Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"main()","metadata":{"execution":{"iopub.status.busy":"2022-11-18T13:47:37.839859Z","iopub.execute_input":"2022-11-18T13:47:37.840283Z","iopub.status.idle":"2022-11-18T14:05:30.018724Z","shell.execute_reply.started":"2022-11-18T13:47:37.840244Z","shell.execute_reply":"2022-11-18T14:05:30.016458Z"},"collapsed":true,"jupyter":{"outputs_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ray.shutdown()","metadata":{"execution":{"iopub.status.busy":"2022-11-18T14:39:16.321185Z","iopub.execute_input":"2022-11-18T14:39:16.322636Z","iopub.status.idle":"2022-11-18T14:39:18.779324Z","shell.execute_reply.started":"2022-11-18T14:39:16.322560Z","shell.execute_reply":"2022-11-18T14:39:18.778034Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Fit with the best parameter","metadata":{}},{"cell_type":"code","source":"train = np.load(\"/kaggle/input/cite-final/new_cite_train_final.npz\")[\"arr_0\"]\ntrain_index = np.load(\"/kaggle/input/multimodal-single-cell-as-sparse-matrix/train_cite_inputs_idxcol.npz\",allow_pickle=True)\n\nmeta = pd.read_csv(\"/kaggle/input/open-problems-multimodal/metadata.csv\",index_col = \"cell_id\")\nmeta = meta[meta.technology==\"citeseq\"]\nlbe = LabelEncoder()\nmeta[\"cell_type\"] = lbe.fit_transform(meta[\"cell_type\"])\nmeta[\"gender\"] = meta.apply(lambda x:0 if x[\"donor\"]==13176 else 1,axis =1)\nmeta_train = meta.reindex(train_index[\"index\"])\n# train_meta = meta_train[\"gender\"].values.reshape(-1, 1)\n# train = np.concatenate([train,train_meta],axis= -1)\ntrain_meta = meta_train[\"cell_type\"].values.reshape(-1, 1)\nohe = OneHotEncoder(sparse=False)\ntrain_meta = ohe.fit_transform(train_meta)\ntrain = np.concatenate([train,train_meta],axis= -1)\n\ntarget = pd.read_hdf(\"/kaggle/input/open-problems-multimodal/train_cite_targets.h5\").values\ntarget -= target.mean(axis=1).reshape(-1, 1)\ntarget /= target.std(axis=1).reshape(-1, 1)\ntrain.shape,target.shape","metadata":{"execution":{"iopub.status.busy":"2022-11-18T14:40:49.297279Z","iopub.execute_input":"2022-11-18T14:40:49.298803Z","iopub.status.idle":"2022-11-18T14:40:53.070412Z","shell.execute_reply.started":"2022-11-18T14:40:49.298750Z","shell.execute_reply":"2022-11-18T14:40:53.068629Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"kf = GroupKFold(n_splits=3) \nindex = 0\nscore = []\n\n# model = Ridge(copy_X=False)\nprint('Train...')\nfor id,(idx_tr, idx_va) in enumerate(kf.split(range(train.shape[0]),groups= meta_train.donor)):\n    if 1:\n        Xtr, Xva = train[idx_tr], train[idx_va]\n        Ytr, Yva = target[idx_tr], target[idx_va]\n        print(f'Fold {id}..')\n        \n        # Since the training process would cost a few hours, so we sample the data here\n        # In order to get the best results, you may fit the whole data\n        # If you fit the whole data, the CV score can reach 0.893259\n        random.seed(42)\n        index_to_train = random.choices([i for i in range(Xtr.shape[0])],k = 1000) \n        Xtr = Xtr[index_to_train]\n        Ytr = Ytr[index_to_train]\n        kernel = RBF(20)\n        model = KernelRidge(kernel=kernel,alpha=1.7)\n        model.fit(Xtr,Ytr)\n    \n        del Xtr, Ytr\n        gc.collect()\n        s = correlation_score(Yva, model.predict(Xva))\n        score.append(s)\n        print(id, s)\n        del Xva, Yva\n        gc.collect()\n        filename = f\"./models/fold{id}.pkl\"\n\n        os.makedirs(\"./models/\",exist_ok=True)\n        index += 1\n        with open(filename,\"wb\") as f:\n            pickle.dump(model,f)\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-11-18T14:42:35.520072Z","iopub.execute_input":"2022-11-18T14:42:35.520648Z","iopub.status.idle":"2022-11-18T14:43:17.418949Z","shell.execute_reply.started":"2022-11-18T14:42:35.520599Z","shell.execute_reply":"2022-11-18T14:43:17.417669Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.mean(score)","metadata":{"execution":{"iopub.status.busy":"2022-11-18T14:43:17.421224Z","iopub.execute_input":"2022-11-18T14:43:17.421651Z","iopub.status.idle":"2022-11-18T14:43:17.432444Z","shell.execute_reply.started":"2022-11-18T14:43:17.421614Z","shell.execute_reply":"2022-11-18T14:43:17.430492Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test = np.load(\"../input/cite-final/new_cite_test_final.npz\")[\"arr_0\"]\ntest_index = np.load(\"../input/multimodal-single-cell-as-sparse-matrix/test_cite_inputs_idxcol.npz\",allow_pickle=True)\nmeta_test = meta.reindex(test_index[\"index\"])\ntest_meta = meta_test[\"cell_type\"].values.reshape(-1, 1)\ntest_meta = ohe.transform(test_meta)\ntest = np.concatenate([test,test_meta],axis= -1)\n\ntest.shape","metadata":{"execution":{"iopub.status.busy":"2022-11-18T14:45:56.791072Z","iopub.execute_input":"2022-11-18T14:45:56.791584Z","iopub.status.idle":"2022-11-18T14:45:57.468688Z","shell.execute_reply.started":"2022-11-18T14:45:56.791544Z","shell.execute_reply":"2022-11-18T14:45:57.467029Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def std(x):\n    return (x - np.mean(x,axis=1).reshape(-1,1)) / np.std(x,axis=1).reshape(-1,1)","metadata":{"execution":{"iopub.status.busy":"2022-11-18T14:46:05.144924Z","iopub.execute_input":"2022-11-18T14:46:05.145486Z","iopub.status.idle":"2022-11-18T14:46:05.153902Z","shell.execute_reply.started":"2022-11-18T14:46:05.145438Z","shell.execute_reply":"2022-11-18T14:46:05.152272Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from tqdm.notebook import tqdm\nimport glob\nmodel_path = \"./models/fold*.pkl\"\nmodel_list = glob.glob(model_path)\npreds = np.zeros((test.shape[0], 140))\nfor id,fn in enumerate(tqdm(model_list)):\n    with open(fn, 'rb') as file:\n        model = pickle.load(file)\n        preds += std(model.predict(test))*score[id]\n        gc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-11-18T14:46:18.408168Z","iopub.execute_input":"2022-11-18T14:46:18.409329Z","iopub.status.idle":"2022-11-18T14:47:20.045580Z","shell.execute_reply.started":"2022-11-18T14:46:18.409264Z","shell.execute_reply":"2022-11-18T14:47:20.044051Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import seaborn as sns\nsns.heatmap(preds)","metadata":{"execution":{"iopub.status.busy":"2022-11-18T14:47:42.589924Z","iopub.execute_input":"2022-11-18T14:47:42.590433Z","iopub.status.idle":"2022-11-18T14:47:52.691689Z","shell.execute_reply.started":"2022-11-18T14:47:42.590392Z","shell.execute_reply":"2022-11-18T14:47:52.690282Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def submit(test_pred,multi_path):\n    submission = pd.read_csv(multi_path,index_col = 0)\n    submission = submission[\"target\"]\n    print(\"data loaded\")\n    submission.iloc[:len(test_pred.ravel())] = test_pred.ravel()\n    assert not submission.isna().any()\n    # submission = submission.round(6) # reduce the size of the csv\n    print(\"start -> submission.zip\")\n    submission.to_csv('submission.zip')\n    print(\"submission.zip saved!\")","metadata":{"execution":{"iopub.status.busy":"2022-11-18T14:48:53.963288Z","iopub.execute_input":"2022-11-18T14:48:53.964806Z","iopub.status.idle":"2022-11-18T14:48:53.974406Z","shell.execute_reply.started":"2022-11-18T14:48:53.964739Z","shell.execute_reply":"2022-11-18T14:48:53.973184Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submit(preds,\"../input/4th-solution-ensemble/submission.zip\")","metadata":{"execution":{"iopub.status.busy":"2022-11-18T14:50:30.989183Z","iopub.execute_input":"2022-11-18T14:50:30.990712Z","iopub.status.idle":"2022-11-18T14:58:22.689743Z","shell.execute_reply.started":"2022-11-18T14:50:30.990635Z","shell.execute_reply":"2022-11-18T14:58:22.688131Z"},"trusted":true},"execution_count":null,"outputs":[]}]}