model_selection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cover Type Prediction: Model selection\n",
    "\n",
    "Sébastien Meyer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datetime import datetime\n",
    "\n",
    "from tqdm import tqdm\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from numpy.random import MT19937, RandomState, SeedSequence\n",
    "\n",
    "import tensorflow as tf\n",
    "\n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn.ensemble import (\n",
    "    RandomForestClassifier, ExtraTreesClassifier, AdaBoostClassifier, ExtraTreesRegressor, StackingClassifier,\n",
    "    VotingClassifier\n",
    ")\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.impute import SimpleImputer\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.metrics import accuracy_score, classification_report\n",
    "from sklearn.model_selection import train_test_split, cross_val_score\n",
    "from sklearn.multiclass import OneVsRestClassifier, OneVsOneClassifier, OutputCodeClassifier\n",
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.preprocessing import PolynomialFeatures, MinMaxScaler, LabelBinarizer, StandardScaler, RobustScaler\n",
    "from sklearn.semi_supervised import LabelPropagation, LabelSpreading\n",
    "from sklearn.svm import SVC\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "from sklearn.utils import shuffle\n",
    "\n",
    "from imblearn.over_sampling import SMOTE\n",
    "from imblearn.under_sampling import ClusterCentroids, RandomUnderSampler\n",
    "\n",
    "from mlxtend.classifier import StackingCVClassifier\n",
    "from mlxtend.feature_selection import SequentialFeatureSelector as SFS\n",
    "\n",
    "import lightgbm\n",
    "from catboost import CatBoostClassifier\n",
    "from lightgbm import LGBMClassifier, LGBMRegressor\n",
    "from xgboost import XGBClassifier\n",
    "\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "font = {\"size\": 22}\n",
    "mpl.rc(\"font\", **font)\n",
    "\n",
    "pd.options.mode.chained_assignment = None  # default='warn'\n",
    "\n",
    "seed = 8005\n",
    "\n",
    "np.random.seed(seed)\n",
    "rs = RandomState(MT19937(SeedSequence(seed)))\n",
    "tf.random.set_seed(seed)\n",
    "\n",
    "merge_kdf = True  # use of knowledge domain features ?\n",
    "merge_pw = True  # use of polynomial features ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Read files\n",
    "test_df = pd.read_csv(\"data/covtype.csv\", index_col=[\"Id\"])\n",
    "\n",
    "training_ids = []\n",
    "\n",
    "with open(\"data/training_ids.txt\", \"r\", encoding=\"utf-8\") as f:\n",
    "\n",
    "    training_ids = f.read().split(\",\")\n",
    "    training_ids = [int(x) for x in training_ids]\n",
    "\n",
    "train_df = test_df.iloc[training_ids, :].copy()\n",
    "\n",
    "# Shuffle\n",
    "train_df = shuffle(train_df, random_state=rs)\n",
    "test_df = shuffle(test_df, random_state=rs)\n",
    "\n",
    "# Eliminate useless features\n",
    "if \"Soil_Type15\" in train_df.columns:\n",
    "    train_df.drop(columns=[\"Soil_Type15\"], inplace=True)\n",
    "    test_df.drop(columns=[\"Soil_Type15\"], inplace=True)\n",
    "\n",
    "# Correct missing values of Hillshade_3pm\n",
    "etr_h3pm = ExtraTreesRegressor(random_state=rs, n_jobs=-1)\n",
    "\n",
    "h3pm_var = \"Hillshade_3pm\"\n",
    "train_h3pm_pos = train_df.index[train_df[h3pm_var] != 0].tolist()\n",
    "train_h3pm_zeros = train_df.index[train_df[h3pm_var] == 0].tolist()\n",
    "test_h3pm_zeros = test_df.index[test_df[h3pm_var] == 0].tolist()\n",
    "\n",
    "etr_h3pm.fit(train_df.drop(columns=[h3pm_var]+[\"Cover_Type\"]).loc[train_h3pm_pos, :],\n",
    "             train_df.loc[train_h3pm_pos, h3pm_var])\n",
    "\n",
    "train_df.loc[train_h3pm_zeros, h3pm_var] = \\\n",
    "    etr_h3pm.predict(train_df.drop(columns=[h3pm_var]+[\"Cover_Type\"]).loc[train_h3pm_zeros, :])\n",
    "test_df.loc[test_h3pm_zeros, h3pm_var] = \\\n",
    "    etr_h3pm.predict(test_df.drop(columns=[h3pm_var]+[\"Cover_Type\"]).loc[test_h3pm_zeros, :])\n",
    "\n",
    "# Binary features and target\n",
    "wild_var = [f\"Wilderness_Area{i}\" for i in range(1, 5)]\n",
    "soil_var = [f\"Soil_Type{i}\" for i in range(1, 41) if i != 15]\n",
    "label_var = [\"Cover_Type\"]\n",
    "\n",
    "# Separate discrete and continuous features\n",
    "all_var = train_df.columns\n",
    "disc_var = wild_var + soil_var + label_var\n",
    "cont_var = [x for x in all_var if x not in disc_var]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Utils"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function allows to drop columns (uniquely) that are too highly correlated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def drop_corr_feat(train_df, test_df=None, corr_threshold=1.):\n",
    "    if corr_threshold >= 1.:\n",
    "        return\n",
    "\n",
    "    print(\"\\nDropping correlated features...\")\n",
    "\n",
    "    # Compute correlation\n",
    "    corr_df = train_df.corr()\n",
    "\n",
    "    # Empty dictionary to hold correlated features\n",
    "    above_threshold_cols = {}\n",
    "\n",
    "    # For each column, record the features that are above the threshold\n",
    "    for col in corr_df:\n",
    "        above_threshold_cols[col] = list(corr_df.index[corr_df[col] > corr_threshold])\n",
    "\n",
    "    # Track columns to remove and columns already examined\n",
    "    cols_to_remove = set()\n",
    "    cols_seen = set()\n",
    "    cols_to_remove_pair = set()\n",
    "\n",
    "    # Iterate through columns and correlated columns\n",
    "    for col, corr_cols in tqdm(above_threshold_cols.items()):\n",
    "\n",
    "        # Keep track of columns already examined\n",
    "        cols_seen.add(col)\n",
    "\n",
    "        for x in corr_cols:\n",
    "\n",
    "            if x != col:  # a variable is totally correlated with itself\n",
    "                # Only want to remove one in a pair\n",
    "                if x not in cols_seen:\n",
    "                    cols_to_remove.add(x)\n",
    "                    cols_to_remove_pair.add(col)\n",
    "\n",
    "    # Remove highly correlated features\n",
    "    list_cols_to_remove = list(cols_to_remove)\n",
    "\n",
    "    train_df.drop(columns=list_cols_to_remove, inplace=True)\n",
    "    if test_df is not None:\n",
    "        test_df.drop(columns=list_cols_to_remove, inplace=True)\n",
    "\n",
    "    print(f\"Number of features after decorrelation: {train_df.shape[1]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Knowledge-domain features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Separate discrete and continuous features\n",
    "all_var = train_df.columns\n",
    "disc_var = wild_var + soil_var + label_var\n",
    "cont_var = [x for x in all_var if x not in disc_var]\n",
    "\n",
    "train_kdf = pd.DataFrame([], index=train_df.index)\n",
    "test_kdf = pd.DataFrame([], index=test_df.index)\n",
    "\n",
    "# # Binary features\n",
    "# soil_var = [\"Soil_Type{}\".format(i) for i in range(1, 41) if i != 15]\n",
    "# wild_var = [\"Wilderness_Area{}\".format(i) for i in range(1, 5)]\n",
    "\n",
    "# Merge Wilderness_Area\n",
    "# s = train_df[wild_var].idxmax(axis=1).str[15:].astype(int) - 1\n",
    "# train_kdf[\"Wilderness_Area\"] = s\n",
    "# train_df = train_df.drop(columns=wild_var)\n",
    "\n",
    "# s = test_df[wild_var].idxmax(axis=1).str[15:].astype(int) - 1\n",
    "# test_kdf[\"Wilderness_Area\"] = s\n",
    "# test_df = test_df.drop(columns=wild_var)\n",
    "\n",
    "# s = train_df[soil_var].idxmax(axis=1).str[9:].astype(int) - 1\n",
    "# train_kdf[\"Soil_Type\"] = s\n",
    "# train_df = train_df.drop(columns=soil_var)\n",
    "\n",
    "# # Add the features on test data\n",
    "\n",
    "# s = test_df[soil_var].idxmax(axis=1).str[9:].astype(int) - 1\n",
    "# test_df[\"Soil_Type\"] = s\n",
    "# test_df = test_df.drop(columns=soil_var)\n",
    "\n",
    "# Add ratio of distances to hydrology\n",
    "# train_kdf[\"Ratio_Distance_To_Hydrology\"] = \\\n",
    "# train_df[\"Vertical_Distance_To_Hydrology\"]/train_df[\"Horizontal_Distance_To_Hydrology\"]\n",
    "# test_kdf[\"Ratio_Distance_To_Hydrology\"] = \\\n",
    "# test_df[\"Vertical_Distance_To_Hydrology\"]/test_df[\"Horizontal_Distance_To_Hydrology\"]\n",
    "\n",
    "# imp = SimpleImputer(strategy=\"median\")  # there might be missing values for the ratio (0 horizontal distance)\n",
    "\n",
    "# train_kdf[[\"Ratio_Distance_To_Hydrology\"]] = imp.fit_transform(train_kdf[[\"Ratio_Distance_To_Hydrology\"]])\n",
    "# test_kdf[[\"Ratio_Distance_To_Hydrology\"]] = imp.transform(test_kdf[[\"Ratio_Distance_To_Hydrology\"]])\n",
    "\n",
    "# Add Log of Elevation\n",
    "# train_kdf[\"Elevation_Log\"] = np.log(1+train_df[\"Elevation\"])\n",
    "# test_kdf[\"Elevation_Log\"] = np.log(1+test_df[\"Elevation\"])\n",
    "\n",
    "# Add Log of distance to hydrology\n",
    "train_kdf[\"Horizontal_Distance_To_Hydrology_Log\"] = np.log(1+train_df[\"Horizontal_Distance_To_Hydrology\"])\n",
    "test_kdf[\"Horizontal_Distance_To_Hydrology_Log\"] = np.log(1+test_df[\"Horizontal_Distance_To_Hydrology\"])\n",
    "\n",
    "# Add Log of distance to roadways\n",
    "train_kdf[\"Horizontal_Distance_To_Roadways_Log\"] = np.log(1+train_df[\"Horizontal_Distance_To_Roadways\"])\n",
    "test_kdf[\"Horizontal_Distance_To_Roadways_Log\"] = np.log(1+test_df[\"Horizontal_Distance_To_Roadways\"])\n",
    "\n",
    "# Add Log of distance to fire points\n",
    "train_kdf[\"Horizontal_Distance_To_Fire_Points_Log\"] = np.log(1+train_df[\"Horizontal_Distance_To_Fire_Points\"])\n",
    "test_kdf[\"Horizontal_Distance_To_Fire_Points_Log\"] = np.log(1+test_df[\"Horizontal_Distance_To_Fire_Points\"])\n",
    "\n",
    "# Add Max of known values (numerical features)\n",
    "train_kdf[\"Max\"] = train_df[cont_var].max(axis=1)\n",
    "test_kdf[\"Max\"] = test_df[cont_var].max(axis=1)\n",
    "\n",
    "# train_kdf[\"Min\"] = train_df[cont_var].min(axis=1)\n",
    "# test_kdf[\"Min\"] = test_df[cont_var].min(axis=1)\n",
    "\n",
    "train_kdf[\"Std\"] = train_df[cont_var].std(axis=1)\n",
    "test_kdf[\"Std\"] = test_df[cont_var].std(axis=1)\n",
    "\n",
    "# train_kdf[\"Mean\"] = train_df[cont_var].mean(axis=1)\n",
    "# test_kdf[\"Mean\"] = test_df[cont_var].mean(axis=1)\n",
    "\n",
    "# Sign of vertical distance\n",
    "# train_kdf[\"Vertical_Distance_To_Hydrology_Sign\"] = (train_df[\"Vertical_Distance_To_Hydrology\"] > 0).astype(int)\n",
    "# test_kdf[\"Vertical_Distance_To_Hydrology_Sign\"] = (test_df[\"Vertical_Distance_To_Hydrology\"] > 0).astype(int)\n",
    "\n",
    "# We have the Aspect variable that is between 0 and 360\n",
    "# train_kdf[\"Shifted_Aspect\"] = train_df[\"Aspect\"] - 180\n",
    "# test_kdf[\"Shifted_Aspect\"] = test_df[\"Aspect\"] - 180\n",
    "\n",
    "# train_df[\"Shifted_Aspect_Sign\"] = (train_kdf[\"Shifted_Aspect\"] > 0).astype(int)\n",
    "# test_df[\"Shifted_Aspect_Sign\"] = (test_kdf[\"Shifted_Aspect\"] > 0).astype(int)\n",
    "\n",
    "# We have the horizontal and vertical distances, let's compute the total distance\n",
    "h_d = \"Horizontal_Distance_To_Hydrology\"\n",
    "v_d = \"Vertical_Distance_To_Hydrology\"\n",
    "train_kdf[\"Distance_To_Hydrology\"] = (train_df[h_d].pow(2) + train_df[v_d].pow(2)).pow(0.5)\n",
    "test_kdf[\"Distance_To_Hydrology\"] = (test_df[h_d].pow(2) + test_df[v_d].pow(2)).pow(0.5)\n",
    "\n",
    "# Make some differences and additions of similar features\n",
    "# hillshades_var = [\"Hillshade_9am\", \"Hillshade_Noon\", \"Hillshade_3pm\"]\n",
    "\n",
    "# for i in range(len(hillshades_var)):\n",
    "#     for j in range(i+1, len(hillshades_var)):\n",
    "\n",
    "#         var1 = hillshades_var[i]\n",
    "#         var2 = hillshades_var[j]\n",
    "\n",
    "#         train_kdf[var1+\"_plus_\"+var2] = train_df[var1] + train_df[var2]\n",
    "#         train_kdf[var1+\"_minus_\"+var2] = np.abs(train_df[var1] - train_df[var2])\n",
    "\n",
    "#         test_kdf[var1+\"_plus_\"+var2] = test_df[var1] + test_df[var2]\n",
    "#         test_kdf[var1+\"_minus_\"+var2] = np.abs(test_df[var1] - test_df[var2])\n",
    "\n",
    "h_dist_var = [\n",
    "    \"Horizontal_Distance_To_Hydrology\", \"Horizontal_Distance_To_Fire_Points\",\n",
    "    \"Horizontal_Distance_To_Roadways\"\n",
    "]\n",
    "\n",
    "for i in range(len(h_dist_var)):\n",
    "    for j in range(i+1, len(h_dist_var)):\n",
    "\n",
    "        var1 = h_dist_var[i]\n",
    "        var2 = h_dist_var[j]\n",
    "\n",
    "        train_kdf[var1+\"_plus_\"+var2] = train_df[var1] + train_df[var2]\n",
    "        train_kdf[var1+\"_minus_\"+var2] = train_df[var1] - train_df[var2]\n",
    "\n",
    "        test_kdf[var1+\"_plus_\"+var2] = test_df[var1] + test_df[var2]\n",
    "        test_kdf[var1+\"_minus_\"+var2] = test_df[var1] - test_df[var2]\n",
    "\n",
    "# Also add the sum of all three distances\n",
    "train_kdf[\"Mean_Distance_To_Points_Of_Interest\"] = \\\n",
    "    train_df[\"Horizontal_Distance_To_Hydrology\"] + train_df[\"Horizontal_Distance_To_Fire_Points\"] + \\\n",
    "    train_df[\"Horizontal_Distance_To_Roadways\"]\n",
    "test_kdf[\"Mean_Distance_To_Points_Of_Interest\"] = \\\n",
    "    test_df[\"Horizontal_Distance_To_Hydrology\"] + test_df[\"Horizontal_Distance_To_Fire_Points\"] + \\\n",
    "    test_df[\"Horizontal_Distance_To_Roadways\"]    \n",
    "    \n",
    "# Our EDA has shown that Elevation and Distances to hydrology share a particular relationship\n",
    "train_kdf[\"Elevation_Shifted_Horizontal_Distance_To_Hydrology\"] = \\\n",
    "    train_df[\"Elevation\"] - 0.2*train_df[\"Horizontal_Distance_To_Hydrology\"]\n",
    "test_kdf[\"Elevation_Shifted_Horizontal_Distance_To_Hydrology\"] = \\\n",
    "    test_df[\"Elevation\"] - 0.2*test_df[\"Horizontal_Distance_To_Hydrology\"]\n",
    "\n",
    "train_kdf[\"Elevation_Shifted_Vertical_Distance_To_Hydrology\"] = \\\n",
    "    train_df[\"Elevation\"] - train_df[\"Vertical_Distance_To_Hydrology\"]\n",
    "test_kdf[\"Elevation_Shifted_Vertical_Distance_To_Hydrology\"] = \\\n",
    "    test_df[\"Elevation\"] - test_df[\"Vertical_Distance_To_Hydrology\"]\n",
    "\n",
    "train_kdf[\"Elevation_Shifted_Horizontal_Distance_To_Roadways\"] = \\\n",
    "    train_df[\"Elevation\"] - 0.02*train_df[\"Horizontal_Distance_To_Roadways\"]\n",
    "test_kdf[\"Elevation_Shifted_Horizontal_Distance_To_Roadways\"] = \\\n",
    "    test_df[\"Elevation\"] - 0.02*test_df[\"Horizontal_Distance_To_Roadways\"]\n",
    "\n",
    "# Binning features\n",
    "# cut_points = [0, 2575, 3100, 8000]\n",
    "# train_kdf[\"Elevation_Plateau\"] = pd.cut(train_df[\"Elevation\"], cut_points, labels=[0, 1, 2]).astype(int)\n",
    "# test_kdf[\"Elevation_Plateau\"] = pd.cut(test_df[\"Elevation\"], cut_points, labels=[0, 1, 2]).astype(int)\n",
    "\n",
    "# Tweaking the binary variables\n",
    "# train_kdf[\"Soil_Type12_32\"] = train_df[\"Soil_Type32\"] + train_df[\"Soil_Type12\"]\n",
    "# test_kdf[\"Soil_Type12_32\"] = test_df[\"Soil_Type32\"] + test_df[\"Soil_Type12\"]\n",
    "\n",
    "# train_kdf[\"Soil_Type23_22_32_33\"] = train_df[\"Soil_Type23\"] + \\\n",
    "#     train_df[\"Soil_Type22\"] + train_df[\"Soil_Type32\"] + train_df[\"Soil_Type33\"]\n",
    "# test_kdf[\"Soil_Type23_22_32_33\"] = test_df[\"Soil_Type23\"] + \\\n",
    "#     test_df[\"Soil_Type22\"] + test_df[\"Soil_Type32\"] + test_df[\"Soil_Type33\"]\n",
    "\n",
    "# Mean hillshade\n",
    "train_kdf[\"Mean_Hillshade\"] = train_df[\"Hillshade_9am\"] + train_df[\"Hillshade_Noon\"] + train_df[\"Hillshade_3pm\"]\n",
    "test_kdf[\"Mean_Hillshade\"] = test_df[\"Hillshade_9am\"] + test_df[\"Hillshade_Noon\"] + test_df[\"Hillshade_3pm\"]\n",
    "\n",
    "# Features drawn from hierarchical clustering\n",
    "train_kdf[\"Aspect Hillshade_3pm\"] = train_df[\"Aspect\"] * train_df[\"Hillshade_3pm\"]\n",
    "test_kdf[\"Aspect Hillshade_3pm\"] = test_df[\"Aspect\"] * test_df[\"Hillshade_3pm\"]\n",
    "\n",
    "# train_kdf[\"Wilderness_Area1_plus_Soil_Type29\"] = train_df[\"Wilderness_Area1\"] + train_df[\"Soil_Type29\"]\n",
    "# test_kdf[\"Wilderness_Area1_plus_Soil_Type29\"] = test_df[\"Wilderness_Area1\"] + test_df[\"Soil_Type29\"]\n",
    "# if \"Wilderness_Area1_plus_Soil_Type29\" not in disc_var:\n",
    "#     disc_var.append(\"Wilderness_Area1_plus_Soil_Type29\")\n",
    "\n",
    "# train_kdf[\"Wilderness_Area1_minus_Soil_Type29\"] = train_df[\"Wilderness_Area1\"] - train_df[\"Soil_Type29\"]\n",
    "# test_kdf[\"Wilderness_Area1_minus_Soil_Type29\"] = test_df[\"Wilderness_Area1\"] - test_df[\"Soil_Type29\"]\n",
    "# if \"Wilderness_Area1_minus_Soil_Type29\" not in disc_var:\n",
    "#     disc_var.append(\"Wilderness_Area1_minus_Soil_Type29\")\n",
    "\n",
    "# train_kdf[\"Wilderness_Area4_plus_Soil_Type3\"] = train_df[\"Wilderness_Area4\"] + train_df[\"Soil_Type3\"]\n",
    "# test_kdf[\"Wilderness_Area4_plus_Soil_Type3\"] = test_df[\"Wilderness_Area4\"] + test_df[\"Soil_Type3\"]\n",
    "# if \"Wilderness_Area4_plus_Soil_Type3\" not in disc_var:\n",
    "#     disc_var.append(\"Wilderness_Area4_plus_Soil_Type3\")\n",
    "\n",
    "# train_kdf[\"Wilderness_Area4_minus_Soil_Type3\"] = train_df[\"Wilderness_Area4\"] - train_df[\"Soil_Type3\"]\n",
    "# test_kdf[\"Wilderness_Area4_minus_Soil_Type3\"] = test_df[\"Wilderness_Area4\"] - test_df[\"Soil_Type3\"]\n",
    "# if \"Wilderness_Area4_minus_Soil_Type3\" not in disc_var:\n",
    "#     disc_var.append(\"Wilderness_Area4_minus_Soil_Type3\")\n",
    "\n",
    "# Associations of soil types\n",
    "ratake = [2, 4]\n",
    "vanet = [2, 5, 6]\n",
    "catamount = [10, 11, 13, 26, 31, 32, 33]\n",
    "leighan = [21, 22, 23, 24, 25, 27, 28, 29, 30, 31, 33, 38]\n",
    "bullwark = [10, 11]\n",
    "como = [29, 30]\n",
    "moran = [38, 39, 40]\n",
    "other = [3, 14, 15, 16, 19, 20, 34, 35, 37]\n",
    "\n",
    "ratake_dict = {i: 1 if i in ratake else 0 for i in range(1, 41)}\n",
    "vanet_dict = {i: 1 if i in vanet else 0 for i in range(1, 41)}\n",
    "catamount_dict = {i: 1 if i in catamount else 0 for i in range(1, 41)}\n",
    "leighan_dict = {i: 1 if i in leighan else 0 for i in range(1, 41)}\n",
    "bullwark_dict = {i: 1 if i in bullwark else 0 for i in range(1, 41)}\n",
    "como_dict = {i: 1 if i in como else 0 for i in range(1, 41)}\n",
    "moran_dict = {i: 1 if i in moran else 0 for i in range(1, 41)}\n",
    "other_dict = {i: 1 if i in other else 0 for i in range(1, 41)}\n",
    "\n",
    "soil_var = [\"Soil_Type{}\".format(i) for i in range(1, 41) if i != 15]\n",
    "\n",
    "train_soil_types = train_df[soil_var].idxmax(axis=1).str[9:].astype(int)\n",
    "\n",
    "train_kdf[\"Ratake_Family_Soil_Type\"] = train_soil_types.map(ratake_dict)\n",
    "train_kdf[\"Vanet_Family_Soil_Type\"] = train_soil_types.map(vanet_dict)\n",
    "train_kdf[\"Catamount_Family_Soil_Type\"] = train_soil_types.map(catamount_dict)\n",
    "train_kdf[\"Leighan_Family_Soil_Type\"] = train_soil_types.map(leighan_dict)\n",
    "train_kdf[\"Bullwark_Family_Soil_Type\"] = train_soil_types.map(bullwark_dict)\n",
    "train_kdf[\"Como_Family_Soil_Type\"] = train_soil_types.map(como_dict)\n",
    "train_kdf[\"Moran_Family_Soil_Type\"] = train_soil_types.map(moran_dict)\n",
    "train_kdf[\"Other_Family_Soil_Type\"] = train_soil_types.map(other_dict)\n",
    "\n",
    "test_soil_types = test_df[soil_var].idxmax(axis=1).str[9:].astype(int)\n",
    "\n",
    "test_kdf[\"Ratake_Family_Soil_Type\"] = test_soil_types.map(ratake_dict)\n",
    "test_kdf[\"Vanet_Family_Soil_Type\"] = test_soil_types.map(vanet_dict)\n",
    "test_kdf[\"Catamount_Family_Soil_Type\"] = test_soil_types.map(catamount_dict)\n",
    "test_kdf[\"Leighan_Family_Soil_Type\"] = test_soil_types.map(leighan_dict)\n",
    "test_kdf[\"Bullwark_Family_Soil_Type\"] = test_soil_types.map(bullwark_dict)\n",
    "test_kdf[\"Como_Family_Soil_Type\"] = test_soil_types.map(como_dict)\n",
    "test_kdf[\"Moran_Family_Soil_Type\"] = test_soil_types.map(moran_dict)\n",
    "test_kdf[\"Other_Family_Soil_Type\"] = test_soil_types.map(other_dict)\n",
    "\n",
    "family_var = [\n",
    "    \"Ratake_Family_Soil_Type\", \"Vanet_Family_Soil_Type\", \"Catamount_Family_Soil_Type\", \"Leighan_Family_Soil_Type\",\n",
    "    \"Bullwark_Family_Soil_Type\", \"Como_Family_Soil_Type\", \"Moran_Family_Soil_Type\", \"Other_Family_Soil_Type\"\n",
    "]\n",
    "# disc_var += family_var\n",
    "\n",
    "# Add rock soil types\n",
    "soil_var = [\"Soil_Type{}\".format(i) for i in range(1, 41) if i != 15]\n",
    "\n",
    "stony_soil_types = [1, 2, 6, 9, 12, 18, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 36, 37, 38, 39, 40]\n",
    "rubly_soil_types = [3, 4, 5, 10, 11, 13]\n",
    "other_soil_types = [7, 8, 14, 15, 16, 17, 19, 20, 21, 22]\n",
    "\n",
    "stony_dict = {i: 1 if i in stony_soil_types else 0 for i in range(1, 41)}\n",
    "rubly_dict = {i: 1 if i in rubly_soil_types else 0 for i in range(1, 41)}\n",
    "other_dict = {i: 1 if i in other_soil_types else 0 for i in range(1, 41)}\n",
    "\n",
    "train_soil_types = train_df[soil_var].idxmax(axis=1).str[9:].astype(int)\n",
    "\n",
    "train_kdf[\"Stony_Soil_Type\"] = train_soil_types.map(stony_dict)\n",
    "train_kdf[\"Rubly_Soil_Type\"] = train_soil_types.map(rubly_dict)\n",
    "# train_kdf[\"Other_Soil_Type\"] = train_soil_types.map(other_dict)\n",
    "\n",
    "test_soil_types = test_df[soil_var].idxmax(axis=1).str[9:].astype(int)\n",
    "\n",
    "test_kdf[\"Stony_Soil_Type\"] = test_soil_types.map(stony_dict)\n",
    "test_kdf[\"Rubly_Soil_Type\"] = test_soil_types.map(rubly_dict)\n",
    "# test_kdf[\"Other_Soil_Type\"] = test_soil_types.map(other_dict)\n",
    "\n",
    "group_var = [\"Stony_Soil_Type\", \"Rubly_Soil_Type\"]\n",
    "# disc_var += group_var\n",
    "\n",
    "# Potential use of engineered variables in polynomial features\n",
    "# kdf_pw = [\"Horizontal_Distance_To_Roadways_Log\",\n",
    "#           \"Elevation_Shifted_Vertical_Distance_To_Hydrology\"]\n",
    "kdf_pw = [\n",
    "    \"Horizontal_Distance_To_Roadways_Log\", \"Horizontal_Distance_To_Fire_Points_Log\",\n",
    "    \"Elevation_Shifted_Vertical_Distance_To_Hydrology\", \"Elevation_Shifted_Horizontal_Distance_To_Hydrology\"\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Polynomial features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the sklearn handler\n",
    "pf = PolynomialFeatures(degree=2, include_bias=False)\n",
    "\n",
    "# Select base features\n",
    "pw_var = [\"Horizontal_Distance_To_Roadways\", \"Elevation\", \"Horizontal_Distance_To_Fire_Points\", \"Wilderness_Area2\"]\n",
    "\n",
    "train_pw = train_df[pw_var].copy()\n",
    "test_pw = test_df[pw_var].copy()\n",
    "\n",
    "# Select knowledge domain features\n",
    "if merge_kdf and len(kdf_pw) > 0:\n",
    "    pw_var += kdf_pw\n",
    "    train_pw.loc[:, kdf_pw] = train_kdf[kdf_pw]\n",
    "    test_pw.loc[:, kdf_pw] = test_kdf[kdf_pw]\n",
    "\n",
    "# Train polynomial features\n",
    "pf.fit(train_pw)\n",
    "\n",
    "# Transform the features\n",
    "train_pw = pf.transform(train_pw)\n",
    "test_pw = pf.transform(test_pw)\n",
    "\n",
    "print(\"PolynomialFeatures shape: \", train_pw.shape)\n",
    "\n",
    "# Re-create the training DataFrame\n",
    "train_pw = pd.DataFrame(train_pw, columns=pf.get_feature_names(pw_var), index=train_df.index)\n",
    "\n",
    "# Re-create the test DataFrame\n",
    "test_pw = pd.DataFrame(test_pw, columns=pf.get_feature_names(pw_var), index=test_df.index)\n",
    "\n",
    "# Create polynomial grown datasets\n",
    "created_pw_var = [x for x in train_pw.columns if x not in kdf_pw+list(train_df.columns)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Merge the features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Merging base and created features\n",
    "\n",
    "if merge_kdf and merge_pw:\n",
    "    \n",
    "    # Add domain knowledge features\n",
    "    train_df_with_kdf = pd.merge(train_df, train_kdf, left_index=True, right_index=True, how=\"left\")\n",
    "    test_df_with_kdf = pd.merge(test_df, test_kdf, left_index=True, right_index=True, how=\"left\")\n",
    "\n",
    "    # Add polynomial features\n",
    "    train_df_final = pd.merge(train_df_with_kdf, train_pw[created_pw_var], left_index=True, right_index=True, how=\"left\")\n",
    "    test_df_final = pd.merge(test_df_with_kdf, test_pw[created_pw_var], left_index=True, right_index=True, how=\"left\")\n",
    "\n",
    "elif not merge_kdf and merge_pw:\n",
    "    \n",
    "    # Add polynomial features\n",
    "    train_df_final = pd.merge(train_df, train_pw[created_pw_var], left_index=True, right_index=True, how=\"left\")\n",
    "    test_df_final = pd.merge(test_df, test_pw[created_pw_var], left_index=True, right_index=True, how=\"left\")\n",
    "\n",
    "elif merge_kdf and not merge_pw:\n",
    "\n",
    "    # Add domain knowledge features\n",
    "    train_df_final = pd.merge(train_df, train_kdf, left_index=True, right_index=True, how=\"left\")\n",
    "    test_df_final = pd.merge(test_df, test_kdf, left_index=True, right_index=True, how=\"left\")\n",
    "    \n",
    "else:\n",
    "    \n",
    "    train_df_final = train_df.copy()\n",
    "    test_df_final = test_df.copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Drop specific features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "drop_var = [\"Soil_Type7\", \"Soil_Type9\", \"Soil_Type25\"]\n",
    "drop_var = []\n",
    "\n",
    "for var in drop_var:\n",
    "    \n",
    "    if var in list(train_df_final.columns):\n",
    "\n",
    "        train_df_final = train_df_final.drop(columns=[var])\n",
    "        test_df_final = test_df_final.drop(columns=[var])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decorrelation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Delete correlated features\n",
    "threshold = 1\n",
    "\n",
    "corr_cols = drop_corr_feat(train_df_final, test_df=test_df_final, corr_threshold=threshold)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the datasets\n",
    "x_train_decorr = train_df_final.drop(columns=label_var)\n",
    "y_train = train_df_final[label_var].to_numpy().flatten()\n",
    "x_test_decorr = test_df_final.drop(columns=label_var)\n",
    "y_test = test_df_final[label_var].to_numpy().flatten()\n",
    "\n",
    "print(\"Number of features: \", train_df_final.shape[1]-1)\n",
    "print(\"Number of features (decorrelated): \", x_train_decorr.shape[1])\n",
    "\n",
    "# Already selected features (previous runs)\n",
    "select_var_idx = [\n",
    "    1, 2, 3, 4, 6, 7, 8, 9, 10, 12, 14, 15, 17, 21, 22, 23, 26, 28, 32, 36, 41, 45, 47, 49, 50, 51, 52, 54, 56, 58, 59, 60,\n",
    "    61, 62, 63, 64, 65, 66, 67, 69, 72, 73, 75, 78, 80, 82, 83, 84, 90, 92, 93, 96, 97, 98, 100, 101, 103, 105, 106, 108\n",
    "]\n",
    "# selected_var = np.array(list(train_df_decorr.columns))[select_var_idx]\n",
    "selected_var = x_train_decorr.columns\n",
    "\n",
    "x_train = x_train_decorr[selected_var]\n",
    "x_test = x_test_decorr[selected_var]\n",
    "\n",
    "# Scale data\n",
    "cont_var = [x for x in list(x_train.columns) if x not in disc_var]\n",
    "sc = StandardScaler()\n",
    "# sc = RobustScaler()\n",
    "# x_train[cont_var] = sc.fit_transform(x_train[cont_var])\n",
    "# x_test[cont_var] = sc.transform(x_test[cont_var])\n",
    "x_train[x_train.columns] = sc.fit_transform(x_train[x_train.columns])\n",
    "x_test[x_test.columns] = sc.transform(x_test[x_test.columns])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca_ratio=-1\n",
    "\n",
    "if pca_ratio > 0:\n",
    "\n",
    "    pca = PCA(n_components=pca_ratio, random_state=rs)\n",
    "\n",
    "    # Fit the PCA\n",
    "    pca.fit(x_train.to_numpy())\n",
    "\n",
    "    # Transform the data (new columns)\n",
    "    pca_col = [\"PCA_{}\".format(i) for i in range(1, pca.n_components_+1)]\n",
    "\n",
    "    x_train = pd.DataFrame(pca.transform(x_train.to_numpy()), columns=pca_col,\n",
    "                            index=x_train.index)\n",
    "    x_test = pd.DataFrame(pca.transform(x_test.to_numpy()), columns=pca_col,\n",
    "                           index=x_test.index)\n",
    "\n",
    "    print(\"Number of features (after PCA): \", x_train.shape[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train test split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Train test split (cross val)\n",
    "x_train_tts, x_eval_tts, y_train_tts, y_eval_tts = train_test_split(\n",
    "    x_train, y_train, test_size=0.2, random_state=seed, shuffle=True\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Try it!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ExtraTreesClassifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Train / validation\n",
    "etc = ExtraTreesClassifier(\n",
    "    n_estimators=300, criterion=\"gini\", max_depth=None, min_samples_split=2, min_samples_leaf=1,\n",
    "    min_weight_fraction_leaf=0., max_features=\"auto\", max_leaf_nodes=None, min_impurity_decrease=1e-7,\n",
    "    bootstrap=False, ccp_alpha=1e-6, random_state=seed, n_jobs=-1, verbose=0\n",
    ")\n",
    "\n",
    "etc.fit(x_train_tts, y_train_tts)\n",
    "y_pred_tts = etc.predict(x_eval_tts)\n",
    "\n",
    "print(\"Accuracy on selected features: \", accuracy_score(y_eval_tts, y_pred_tts))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Submission\n",
    "etc = ExtraTreesClassifier(\n",
    "    n_estimators=300, criterion=\"gini\", max_depth=None, min_samples_split=2, min_samples_leaf=1,\n",
    "    min_weight_fraction_leaf=0., max_features=\"auto\", max_leaf_nodes=None, min_impurity_decrease=1e-7,\n",
    "    bootstrap=False, ccp_alpha=1e-6, random_state=seed, n_jobs=-1, verbose=0\n",
    ")\n",
    "\n",
    "etc.fit(x_train, y_train)\n",
    "y_pred = etc.predict(x_test)\n",
    "\n",
    "print(\"Accuracy on test set: \", accuracy_score(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RandomForestClassifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Train / validation\n",
    "rfc = RandomForestClassifier(n_estimators=100, random_state=seed, n_jobs=-1, verbose=0)\n",
    "\n",
    "rfc.fit(x_train_tts, y_train_tts)\n",
    "y_pred_tts = rfc.predict(x_eval_tts)\n",
    "\n",
    "print(\"Accuracy on train-test-split: \", accuracy_score(y_eval_tts, y_pred_tts))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Submission\n",
    "rfc = RandomForestClassifier(n_estimators=100, random_state=seed, n_jobs=-1, verbose=0)\n",
    "\n",
    "rfc.fit(x_train, y_train)\n",
    "y_pred = rfc.predict(x_test)\n",
    "\n",
    "print(\"Accuracy on test set: \", accuracy_score(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### LightGBM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def evalacc(y_true, y_pred):\n",
    "    preds = y_pred.reshape(len(np.unique(y_true)), -1)\n",
    "    preds = preds.argmax(axis = 0)\n",
    "    return \"Accuracy\", accuracy_score(y_true, preds), True\n",
    "\n",
    "# Train / validation\n",
    "lgbm = LGBMClassifier(\n",
    "    n_estimators=100, num_leaves=73, min_split_gain=1.22e-5, min_child_weight=3.82e-5, \n",
    "    min_child_samples=12, subsample=0.90, subsample_freq=2, reg_alpha=1.14e-6, reg_lambda=5.37e-5,\n",
    "    random_state=rs, n_jobs=-1\n",
    ")\n",
    "\n",
    "lgbm.fit(\n",
    "    x_train_tts, y_train_tts,\n",
    "    eval_set=[(x_eval_tts, y_eval_tts), (x_train_tts, y_train_tts)],\n",
    "    eval_names=[\"Validation\", \"Training\"],\n",
    "    verbose=10,\n",
    "    eval_metric=[\"logloss\", evalacc]\n",
    ")\n",
    "y_pred_tts = lgbm.predict(x_eval_tts)\n",
    "\n",
    "print(\"Accuracy on selected features: \", accuracy_score(y_eval_tts, y_pred_tts))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lightgbm.plot_metric(lgbm, metric=\"Accuracy\", figsize=(16, 8), title=None)\n",
    "# plt.savefig(\"report/figures/lgbmeval.png\", facecolor=\"white\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Submission\n",
    "lgbm = LGBMClassifier(\n",
    "    n_estimators=177, num_leaves=73, min_split_gain=1.22e-5, min_child_weight=3.82e-5,\n",
    "    min_child_samples=12, subsample=0.90, subsample_freq=2, reg_alpha=1.14e-6, reg_lambda=5.37e-5,\n",
    "    random_state=rs, n_jobs=-1\n",
    ")\n",
    "\n",
    "lgbm.fit(x_train, y_train)\n",
    "y_pred = lgbm.predict(x_test)\n",
    "\n",
    "print(\"Accuracy on test set: \", accuracy_score(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stacking"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Training / validation\n",
    "etc = ExtraTreesClassifier(\n",
    "    n_estimators=300, criterion=\"gini\", max_depth=None, min_samples_split=2, min_samples_leaf=1,\n",
    "    min_weight_fraction_leaf=0., max_features=\"auto\", max_leaf_nodes=None, min_impurity_decrease=1e-7,\n",
    "    bootstrap=False, ccp_alpha=1e-6, random_state=rs, n_jobs=-1, verbose=0\n",
    ")\n",
    "rfc = RandomForestClassifier(\n",
    "    n_estimators=80, criterion=\"gini\", min_samples_split=3, ccp_alpha=1e-5,\n",
    "    random_state=rs, n_jobs=-1, verbose=0\n",
    ")\n",
    "lgbm = LGBMClassifier(\n",
    "    n_estimators=100, num_leaves=70, min_split_gain=4e-5, min_child_weight=1.5e-5, min_child_samples=6,\n",
    "    subsample=0.975, subsample_freq=6, reg_alpha=5.5e-4, reg_lambda=4.5e-4, random_state=rs, n_jobs=-1\n",
    ")\n",
    "\n",
    "est = [(\"rfc\", rfc), (\"extra\", etc), (\"lgbm\", lgbm)]\n",
    "\n",
    "final_est = LogisticRegression(multi_class=\"auto\", solver=\"newton-cg\")\n",
    "\n",
    "stacking = StackingClassifier(\n",
    "    est, final_estimator=final_est, cv=5, n_jobs=-1, verbose=1, stack_method=\"predict_proba\", passthrough=True\n",
    ")\n",
    "\n",
    "stacking.fit(x_train_tts, y_train_tts)\n",
    "y_pred_tts = stacking.predict(x_eval_tts)\n",
    "\n",
    "print(\"Accuracy on selected features: \", accuracy_score(y_eval_tts, y_pred_tts))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Submission\n",
    "etc = ExtraTreesClassifier(\n",
    "    n_estimators=300, criterion=\"gini\", max_depth=None, min_samples_split=2, min_samples_leaf=1,\n",
    "    min_weight_fraction_leaf=0., max_features=\"auto\", max_leaf_nodes=None, min_impurity_decrease=1e-7,\n",
    "    bootstrap=False, ccp_alpha=1e-6, random_state=rs, n_jobs=-1, verbose=0\n",
    ")\n",
    "rfc = RandomForestClassifier(\n",
    "    n_estimators=80, criterion=\"gini\", min_samples_split=3, ccp_alpha=1e-5,\n",
    "    random_state=rs, n_jobs=-1, verbose=0\n",
    ")\n",
    "lgbm = LGBMClassifier(\n",
    "    n_estimators=100, num_leaves=70, min_split_gain=4e-5, min_child_weight=1.5e-5, min_child_samples=6,\n",
    "    subsample=0.975, subsample_freq=6, reg_alpha=5.5e-4, reg_lambda=4.5e-4, random_state=rs, n_jobs=-1\n",
    ")\n",
    "\n",
    "est = [(\"rfc\", rfc), (\"extra\", etc), (\"lgbm\", lgbm)]\n",
    "\n",
    "final_est = LogisticRegression(multi_class=\"auto\", solver=\"newton-cg\")\n",
    "\n",
    "stacking = StackingClassifier(\n",
    "    est, final_estimator=final_est, cv=5, n_jobs=-1, verbose=2, passthrough=True\n",
    ")\n",
    "\n",
    "stacking.fit(x_train, y_train)\n",
    "y_pred = stacking.predict(x_test)\n",
    "\n",
    "print(\"Accuracy on test set: \", accuracy_score(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SVC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Training / validation\n",
    "svc = SVC(\n",
    "    kernel=\"rbf\", C=1000, gamma=\"scale\", shrinking=True, decision_function_shape=\"ovr\", break_ties=False,\n",
    "    random_state=seed, verbose=True\n",
    ")\n",
    "\n",
    "svc.fit(x_train_tts, y_train_tts)\n",
    "y_pred_tts = svc.predict(x_eval_tts)\n",
    "\n",
    "print(\"Accuracy on selected features: \", accuracy_score(y_eval_tts, y_pred_tts))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Catboost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Training / validation\n",
    "cat = CatBoostClassifier(iterations=5000, eval_metric=\"Accuracy\", random_state=seed)\n",
    "\n",
    "cat.fit(\n",
    "    x_train_tts, y_train_tts, \n",
    "    eval_set=(x_eval_tts, y_eval_tts), verbose_eval=50\n",
    ")\n",
    "y_pred_tts = cat.predict(x_eval_tts)\n",
    "\n",
    "print(\"Accuracy on selected features: \", accuracy_score(y_eval_tts, y_pred_tts))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Submission\n",
    "cat = CatBoostClassifier(iterations=5000, eval_metric=\"Accuracy\", random_state=seed)\n",
    "\n",
    "cat.fit(\n",
    "    x_train, y_train,\n",
    "    eval_set=(x_test, y_test), verbose_eval=50\n",
    ")\n",
    "y_pred = cat.predict(x_test)\n",
    "\n",
    "print(\"Accuracy on selected features: \", accuracy_score(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sample weights"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having a look at the submission and/or the data set easily proves that the data set is highly imbalanced. We can first try to set weights to the classes during training."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Undersampling / Oversampling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another well-used technique to cope with imbalanced data sets is to use undersampling, oversampling or a mix of both."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SMOTE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_smotes = 3\n",
    "\n",
    "# SMOTE Oversampling\n",
    "x_train_smote = x_train.copy()\n",
    "y_train_smote = y_train.copy()\n",
    "\n",
    "for k in range(n_smotes):\n",
    "\n",
    "    # Fit an ExtraTreesClassifier as usual on whole training set\n",
    "    etc = ExtraTreesClassifier(\n",
    "        n_estimators=300, criterion=\"gini\", max_depth=None, min_samples_split=2, min_samples_leaf=1,\n",
    "        min_weight_fraction_leaf=0., max_features=\"auto\", max_leaf_nodes=None, min_impurity_decrease=1e-7,\n",
    "        bootstrap=False, ccp_alpha=1e-6, random_state=seed, n_jobs=-1, verbose=0\n",
    "    )\n",
    "    etc.fit(x_train_smote, y_train_smote)\n",
    "\n",
    "    # Estimate repartition of the labels in the test set\n",
    "    y_pred = etc.predict(x_test)\n",
    "\n",
    "    _, repart_pred = np.unique(y_pred, return_counts=True)\n",
    "\n",
    "    repart_pred = repart_pred.astype(float)/len(y_test)\n",
    "    min_class = np.argmin(repart_pred)+1\n",
    "\n",
    "    repart_dict = {i: int(2160*repart_pred[i-1]/repart_pred[min_class-1]) for i in range(1, 8)}\n",
    "\n",
    "    print(\"Estimated repartition in test set after {} SMOTE: {}.\".format(k+1, repart_dict))\n",
    "\n",
    "    # SMOTE\n",
    "    smote = SMOTE(sampling_strategy=repart_dict, random_state=rs, n_jobs=-1)\n",
    "\n",
    "    x_train_smote, y_train_smote = smote.fit_resample(x_train, y_train)\n",
    "\n",
    "    # Submission\n",
    "    etc = ExtraTreesClassifier(\n",
    "        n_estimators=300, criterion=\"gini\", max_depth=None, min_samples_split=2, min_samples_leaf=1,\n",
    "        min_weight_fraction_leaf=0., max_features=\"auto\", max_leaf_nodes=None, min_impurity_decrease=1e-7,\n",
    "        bootstrap=False, ccp_alpha=1e-6, random_state=seed, n_jobs=-1, verbose=0\n",
    "    )\n",
    "    etc.fit(x_train_smote, y_train_smote)\n",
    "    y_pred = etc.predict(x_test)\n",
    "    print(\"Accuracy on test set: \", accuracy_score(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clustering + SMOTE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_op = 2\n",
    "\n",
    "# Undersampling/Oversampling\n",
    "x_train_clust = x_train.copy()\n",
    "y_train_clust = y_train.copy()\n",
    "\n",
    "for k in range(n_op):\n",
    "\n",
    "    # Fit an ExtraTreesClassifier as usual on whole training set\n",
    "    etc = ExtraTreesClassifier(\n",
    "        n_estimators=300, criterion=\"gini\", max_depth=None, min_samples_split=2, min_samples_leaf=1,\n",
    "        min_weight_fraction_leaf=0., max_features=\"auto\", max_leaf_nodes=None, min_impurity_decrease=1e-7,\n",
    "        bootstrap=False, ccp_alpha=1e-6, random_state=seed, n_jobs=-1, verbose=0\n",
    "    )\n",
    "    etc.fit(x_train_clust, y_train_clust)\n",
    "\n",
    "    # Estimate repartition of the labels in the test set\n",
    "    y_pred = etc.predict(x_test)\n",
    "\n",
    "    _, repart_pred = np.unique(y_pred, return_counts=True)\n",
    "\n",
    "    repart_pred = repart_pred.astype(float)/len(y_test)\n",
    "    min_class = np.argmin(repart_pred)+1\n",
    "    maj_class = np.argmax(repart_pred)+1\n",
    "    \n",
    "    n_souhaite = 15120\n",
    "    \n",
    "    repart_dict = {i: int(repart_pred[i-1]*n_souhaite) for i in range(1, 8)}\n",
    "    \n",
    "    repart_clust = {i: min(2160, repart_dict[i]) for i in range(1, 8)}\n",
    "    \n",
    "    print(\"Estimated repartition in test set after {} undersampling: {}.\".format(k+1, repart_clust))\n",
    "    print(\"Estimated repartition in test set after {} oversampling: {}.\".format(k+1, repart_dict))\n",
    "    \n",
    "    # ClusteringCentroids + SMOTE\n",
    "    kmeans = KMeans(random_state=rs)\n",
    "    clust = ClusterCentroids(sampling_strategy=repart_clust, random_state=rs, estimator=kmeans)\n",
    "\n",
    "    x_train_clust, y_train_clust = clust.fit_resample(x_train, y_train)\n",
    "    \n",
    "    smote = SMOTE(sampling_strategy=repart_dict, random_state=rs, n_jobs=-1)\n",
    "    \n",
    "    x_train_clust, y_train_clust = smote.fit_resample(x_train_clust, y_train_clust)\n",
    "\n",
    "    # Submission\n",
    "    etc = ExtraTreesClassifier(\n",
    "        n_estimators=300, criterion=\"gini\", max_depth=None, min_samples_split=2, min_samples_leaf=1,\n",
    "        min_weight_fraction_leaf=0., max_features=\"auto\", max_leaf_nodes=None, min_impurity_decrease=1e-7,\n",
    "        bootstrap=False, ccp_alpha=1e-6, random_state=seed, n_jobs=-1, verbose=0\n",
    "    )\n",
    "    etc.fit(x_train_clust, y_train_clust)\n",
    "    y_pred = etc.predict(x_test)\n",
    "    print(\"Accuracy on test set: \", accuracy_score(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Manual sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Manual resampling\n",
    "x_train_hand = x_train.copy()\n",
    "y_train_hand = y_train.copy()\n",
    "\n",
    "repart_dict = {1: 5500, 2: 7000, 3: 1000, 4: 50, 5: 450, 6: 500, 7: 650}\n",
    "repart_clust = {i: min(2160, repart_dict[i]) for i in range(1, 8)}\n",
    "\n",
    "print(\"Estimated repartition in training set after operations: {}.\".format(repart_dict))\n",
    "\n",
    "# ClusteringCentroids + SMOTE\n",
    "kmeans = KMeans(random_state=rs)\n",
    "clust = ClusterCentroids(sampling_strategy=repart_clust, random_state=rs, estimator=kmeans)\n",
    "\n",
    "x_train_hand, y_train_hand = clust.fit_resample(x_train, y_train)\n",
    "\n",
    "smote = SMOTE(sampling_strategy=repart_dict, random_state=rs, n_jobs=-1)\n",
    "\n",
    "x_train_hand, y_train_hand = smote.fit_resample(x_train_hand, y_train_hand)\n",
    "\n",
    "# Submission\n",
    "etc = ExtraTreesClassifier(\n",
    "    n_estimators=300, criterion=\"gini\", max_depth=None, min_samples_split=2, min_samples_leaf=1,\n",
    "    min_weight_fraction_leaf=0., max_features=\"auto\", max_leaf_nodes=None, min_impurity_decrease=1e-7,\n",
    "    bootstrap=False, ccp_alpha=1e-6, random_state=seed, n_jobs=-1, verbose=0\n",
    ")\n",
    "etc.fit(x_train_hand, y_train_hand)\n",
    "y_pred = etc.predict(x_test)\n",
    "print(\"Accuracy on test set: \", accuracy_score(y_test, y_pred))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}