{"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":"### The Fine Art of Fine Tuning\n\nOn this notebook we attempt to further optimize the amazing [notebook](https://www.kaggle.com/code/cabaxiom/tps-jul-22-bgmm-semi-supervised) by [cabaxiom](https://www.kaggle.com/cabaxiom).\n- First, we follow the approach:\n    - Preprocess \n    - Find number of clusters\n    - Select Features\n    - Train multiple GMMs\n    - Train \"semi-supervised\" ensemble of models on the high confidence predictions\n\n- Then, we search (brute-force) the optimal weights for ensembling our 1st stage models.\n- For each weight combination\n    - We further optimize our predictions by fitting a GMM classifier -> Fixing the predictions -> Fitting a GMM classifier -> Fixing the.. (Following the original notebook)","metadata":{"papermill":{"duration":0.005241,"end_time":"2022-07-17T02:43:50.776380","exception":false,"start_time":"2022-07-17T02:43:50.771139","status":"completed"},"tags":[]}},{"cell_type":"code","source":"!pip install sklego\nimport pickle\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nfrom lightgbm import LGBMClassifier\nfrom sklearn.metrics import accuracy_score\nfrom sklego.mixture import BayesianGMMClassifier\nfrom sklearn.ensemble import ExtraTreesClassifier\nfrom sklearn.preprocessing import PowerTransformer\nfrom sklearn.model_selection import StratifiedKFold\nfrom sklearn.mixture import GaussianMixture, BayesianGaussianMixture\nfrom sklearn.metrics import silhouette_score, calinski_harabasz_score, davies_bouldin_score\nfrom sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis, LinearDiscriminantAnalysis","metadata":{"papermill":{"duration":0.031884,"end_time":"2022-07-17T03:10:58.833931","exception":false,"start_time":"2022-07-17T03:10:58.802047","status":"completed"},"tags":[],"_kg_hide-input":true,"jupyter":{"source_hidden":true,"outputs_hidden":true},"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-07-26T09:02:28.577149Z","iopub.execute_input":"2022-07-26T09:02:28.577481Z","iopub.status.idle":"2022-07-26T09:02:39.863733Z","shell.execute_reply.started":"2022-07-26T09:02:28.577453Z","shell.execute_reply":"2022-07-26T09:02:39.862612Z"},"collapsed":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Configuration\n\n- **`FULL_TUNE`:** Run the full brute-force or just load the current optimal params found?\n- **`GEN_PREDICTIONS`:** Generate the *first stage* predictions? \n- **`BEST_COLS`:** Features to use from the original dataset.","metadata":{}},{"cell_type":"code","source":"FULL_TUNE = False\nGEN_PREDICTIONS = False\nBEST_COLS = ['f_07', 'f_08', 'f_09', 'f_10', 'f_11', 'f_12', 'f_13', 'f_22', 'f_23', 'f_24', 'f_25', 'f_26', 'f_27', 'f_28']","metadata":{"papermill":{"duration":0.031884,"end_time":"2022-07-17T03:10:58.833931","exception":false,"start_time":"2022-07-17T03:10:58.802047","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2022-07-26T09:02:39.866099Z","iopub.execute_input":"2022-07-26T09:02:39.866758Z","iopub.status.idle":"2022-07-26T09:02:39.871957Z","shell.execute_reply.started":"2022-07-26T09:02:39.866717Z","shell.execute_reply":"2022-07-26T09:02:39.871005Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def components_number_multiple(max_n, n_seeds):\n    bic_scores = []\n    for n in range(2,max_n):\n        bic_scores_n = []\n        for seed in range(n_seeds):\n            gmm = GaussianMixture(n_components=n, covariance_type = 'full', n_init=3, random_state=seed)\n            gmm.fit(X_scaled)\n            bic_scores_n.append(gmm.bic(X_scaled))\n        bic_scores.append(bic_scores_n)\n    return bic_scores\n\ndef plot_components_number_multiple(max_n, n_seeds):\n    bic_scores = components_number_multiple(max_n + 1, n_seeds)\n    bic_df = pd.DataFrame(data = bic_scores).T\n    bic_df.columns = range(2,max_n+1)\n    f,ax = plt.subplots(figsize=(20,7))\n    for i in range(n_seeds):\n        sns.lineplot(x=bic_df.columns, y=bic_df.loc[i].values)\n    ax.set_xticks(range(2,max_n+1))\n    return bic_df\n\ndef score_clusters(X, predictions, silhouette = True, verbose=False):\n    db_score = davies_bouldin_score(X=X, labels=predictions)\n    ch_score = calinski_harabasz_score(X=X, labels=predictions)\n    s_score = silhouette_score(X=X, labels=predictions, metric='euclidean')\n    if verbose:\n        print(\"David Bouldin score: {0:0.4f}\".format(db_score))\n        print(\"Calinski Harabasz score: {0:0.3f}\".format(ch_score))\n        print(\"Silhouette score: {0:0.4f}\".format(s_score))\n    return db_score, ch_score, s_score\n\ndef soft_voting(predict_number, best_cols = BEST_COLS):\n    predicted_probabilities = pd.DataFrame(np.zeros((len(df),7)), columns=range(1,8))\n    for i in range(predict_number):\n        print(\"=========\", i, \"==========\")\n        X_scaled_sample = X_scaled.sample(40000)\n        gmm = BayesianGaussianMixture(n_components=7, covariance_type = 'full', weight_concentration_prior_type=\"dirichlet_distribution\", max_iter=300, init_params=\"kmeans\", n_init=3, random_state=i)\n        gmm.fit(X_scaled_sample[best_cols])\n        pred_probs = gmm.predict_proba(X_scaled[best_cols])\n        pred_probs = pd.DataFrame(pred_probs, columns=range(1,8))\n        if i == 0:\n            initial_centers = gmm.means_\n        new_classes = []\n        for mean2 in gmm.means_:\n            distances = [np.linalg.norm(mean1-mean2) for mean1 in initial_centers]\n            new_class = np.argmin(distances) + 1\n            new_classes.append(new_class)\n        if len(new_classes) != len(set(new_classes)):\n            print(\"iteration\", i, \"could not determine the cluster label mapping, skipping\")\n            continue\n        pred_probs = pred_probs.rename(columns=dict(zip(range(1,8),new_classes)))\n        predicted_probabilities = predicted_probabilities + pred_probs\n        score_clusters(X_scaled[best_cols], predicted_probabilities.idxmax(axis=1), verbose=True)\n    predicted_probabilities = predicted_probabilities.div(predicted_probabilities.sum(axis=1), axis=0)\n    return predicted_probabilities\n\ndef best_class(df):\n    new_df = df.copy()\n    new_df[\"highest_prob\"] = df.max(axis=1)\n    new_df[\"best_class\"] = df.idxmax(axis=1)\n    new_df[\"second_highest_prob\"] = df.apply(lambda x: x.nlargest(2).values[-1], axis=1)\n    new_df[\"second_best_class\"] = df.apply(lambda x: np.where(x == x.nlargest(2).values[-1])[0][0]+1, axis=1)\n    return new_df\n\ndef k_fold_cv(model,X,y, verbose=True):\n    kfold = StratifiedKFold(n_splits = 5, shuffle=True, random_state = 0)\n    feature_imp, y_pred_list, y_true_list, acc_list  = [],[],[],[]\n    for fold, (train_index, val_index) in enumerate(kfold.split(X, y)):\n        if verbose: print(\"==fold==\", fold)\n        X_train = X.loc[train_index]\n        X_val = X.loc[val_index]\n        y_train = y.loc[train_index]\n        y_val = y.loc[val_index]\n        model.fit(X_train,y_train)\n        y_pred = model.predict(X_val)\n        y_pred_list = np.append(y_pred_list, y_pred)\n        y_true_list = np.append(y_true_list, y_val)\n        acc_list.append(accuracy_score(y_pred, y_val))\n        if verbose: print('Acc', accuracy_score(y_pred, y_val))\n        try: feature_imp.append(model.feature_importances_)\n        except AttributeError: pass\n            \n    return feature_imp, y_pred_list, y_true_list, acc_list, X_val, y_val\n\ndef evaluate_models():\n    for model_name, model in models.items():\n        print(\"===\",model_name,\"===\")\n        feature_imp, y_pred_list, y_true_list, acc_list, X_val, y_val = k_fold_cv(model=model,X=X,y=y, verbose=False)\n        acc_score = accuracy_score(y_pred_list, y_true_list)\n        print(\"{0:0.4f}\".format(acc_score))\n\ndef fit_predict_all():\n    predictions = []\n    model_names = []\n    scores = []\n    for model_name, model in models.items():\n        print(\"===\",model_name,\"===\")\n        model.fit(X[BEST_COLS], y)\n        preds_prob =  model.predict_proba(X_full[BEST_COLS])\n        preds_prob_df = pd.DataFrame(preds_prob, columns=range(1,8), index=X_scaled.index)\n        db, ch, s = score_clusters(X_scaled[BEST_COLS], preds_prob_df.idxmax(axis=1), verbose=True)\n        scores.append((db,ch,s))\n        predictions.append(preds_prob_df)\n        model_names.append(model_name)\n    return predictions, model_names, scores","metadata":{"papermill":{"duration":0.031884,"end_time":"2022-07-17T03:10:58.833931","exception":false,"start_time":"2022-07-17T03:10:58.802047","status":"completed"},"tags":[],"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-07-26T09:02:39.873626Z","iopub.execute_input":"2022-07-26T09:02:39.874333Z","iopub.status.idle":"2022-07-26T09:02:39.909782Z","shell.execute_reply.started":"2022-07-26T09:02:39.874294Z","shell.execute_reply":"2022-07-26T09:02:39.908687Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Dataset Preprocessing\n\n- Following the best practice of this competition: We simply use `PowerTransformer` since it yields the best LB score.\n","metadata":{}},{"cell_type":"code","source":"df = pd.read_csv(\"../input/tabular-playground-series-jul-2022/data.csv\")\ndf = df.drop(columns=\"id\")\n\nint_cols = [i for i in df.columns if df[i].dtype == int]\nfloat_cols = [i for i in df.columns if df[i].dtype == float]\n\ntransformer = PowerTransformer()\nX_scaled = transformer.fit_transform(df)\nX_scaled = pd.DataFrame(X_scaled, columns = df.columns)","metadata":{"papermill":{"duration":0.031884,"end_time":"2022-07-17T03:10:58.833931","exception":false,"start_time":"2022-07-17T03:10:58.802047","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2022-07-26T09:02:39.913872Z","iopub.execute_input":"2022-07-26T09:02:39.914163Z","iopub.status.idle":"2022-07-26T09:02:44.110580Z","shell.execute_reply.started":"2022-07-26T09:02:39.914133Z","shell.execute_reply":"2022-07-26T09:02:44.109628Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 1st Stage Models\n\n- We first attempt to find the optimal number of clusters using a Gaussian Mixture Model. (our goal is to minimise the BIC score)\n- We then plot the features and abserve how distinct the different classes distributions are: The more distinct: the more useful that feature is.\n- We then fit Multiple Bayesian Gausian Mixture Models.\n- We then follow the \"semi-supervised\" approach of fitting a classifier using on the samples we are confident that we got correct. Then use it to predict the points which we are not sure about.\n- **On this satage we use 5 different models:**    \n    - LGBMClassifier\n    - ExtraTreesClassifier\n    - BayesianGMMClassifier\n    - LinearDiscriminantAnalysis\n    - QuadraticDiscriminantAnalysis    ","metadata":{}},{"cell_type":"code","source":"if GEN_PREDICTIONS:\n    bgmm = BayesianGaussianMixture(n_components=7, covariance_type = 'full', n_init=3, random_state=2)\n    predicted_class = bgmm.fit_predict(X_scaled)\n    df[\"class\"] = predicted_class\n\n    pred_probs = soft_voting(10)\n    score_clusters(X_scaled[BEST_COLS], pred_probs.idxmax(axis=1), verbose=True)\n\n    cluster_class_probs = best_class(pred_probs)\n    second_highest_probs_sum = cluster_class_probs.groupby([\"best_class\",\"second_best_class\"])[\"second_highest_prob\"].sum().reset_index()\n    confident_predictions = cluster_class_probs.loc[cluster_class_probs[\"highest_prob\"] >= 0.8]\n    confident_predictions_class = confident_predictions[\"best_class\"]\n    X_scaled[\"class\"] = confident_predictions_class\n\n    train_df = X_scaled.loc[X_scaled[\"class\"] == X_scaled[\"class\"]]\n    test_df = X_scaled.loc[X_scaled[\"class\"] != X_scaled[\"class\"]]\n\n    X = train_df.drop(columns=\"class\").reset_index(drop=True)\n    y = train_df[\"class\"].reset_index(drop=True)\n    X_test = test_df.drop(columns=\"class\").reset_index(drop=True)\n    X_full = X_scaled.drop(columns=\"class\")\n\n    model_et = ExtraTreesClassifier(n_estimators = 2000, n_jobs = -1, random_state = 42)\n    model_lgbm = LGBMClassifier(objective = 'multiclass', n_estimators = 5000, random_state = 42, learning_rate = 0.1, n_jobs = -1)\n    model_qda = QuadraticDiscriminantAnalysis()\n    model_lda = LinearDiscriminantAnalysis()\n    model_bgmm = BayesianGMMClassifier(n_components = 7, random_state = 1, tol = 1e-3, covariance_type = 'full', max_iter = 400, n_init = 4, init_params = 'kmeans')\n\n    models = {\"ET\":model_et, \"LGBM\":model_lgbm, \"QDA\":model_qda, \"LDA\":model_lda, \"BGMM_C\":model_bgmm}\n\n    evaluate_models()\n    feature_imp, y_pred_list, y_true_list, acc_list, X_val, y_val = k_fold_cv(model=model_lgbm,X=X,y=y)\n\n    predictions, model_names, scores = fit_predict_all()\n    cluster_class_probs = cluster_class_probs.loc[:,[1,2,3,4,5,6,7]]\n    predictions.append(cluster_class_probs)\n    model_names.append(\"BGMM\")\n\n    db, ch, s = score_clusters(X_scaled[BEST_COLS], cluster_class_probs.idxmax(axis=1), verbose=True)\n    scores.append((db,ch,s))\n\n    pickle.dump(predictions, open('predictions.pkl', 'wb'))","metadata":{"papermill":{"duration":0.031884,"end_time":"2022-07-17T03:10:58.833931","exception":false,"start_time":"2022-07-17T03:10:58.802047","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2022-07-26T09:02:44.112399Z","iopub.execute_input":"2022-07-26T09:02:44.113042Z","iopub.status.idle":"2022-07-26T09:02:44.128631Z","shell.execute_reply.started":"2022-07-26T09:02:44.113002Z","shell.execute_reply":"2022-07-26T09:02:44.127645Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2nd Stage Models\n- Now that we got our classifiers predictions we combine them using **weighted avg**.\n- Then we in order to improve our performance even further we iteratively use the the predicted labels from the previous iteration's model as our training labels for the current iteration's model.\n\n**On this stage we search for the optimal weights for combining our models**\n> This is a long and on-going search (iteration: +1h) that is currently running on a local machine. The best results found will be updated here.","metadata":{}},{"cell_type":"code","source":"if FULL_TUNE:\n    predictions = pickle.load(open('../input/tps-jul-final-search/predictions.pkl', 'rb'))\n\n    # Yes. I know.\n    for w_1 in [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0]:\n        for w_2 in [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0]:\n            for w_3 in [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0]:\n                for w_4 in [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0]:\n                    for w_5 in [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0]:\n\n                        predictions_df = w_1 * predictions[0] + w_2 * predictions[1] + w_3 * predictions[2] + w_4 * predictions[4] + w_5 * predictions[5]\n                        predictions_df = predictions_df.div(predictions_df.sum(axis = 1), axis = 0)\n                        predictions_df = best_class(predictions_df)\n\n                        db, ch, s = score_clusters(X_scaled[BEST_COLS], predictions_df[\"best_class\"], verbose = True)\n                        scores.append((db,ch,s))\n                        model_names.append(\"combined\")\n                        pd.DataFrame(scores, index=model_names, columns=[\"Davies-Bouldin Index\",\"Calinski-Harabasz Index\",\"Silhouette Coefficient\"])\n                        second_highest_probs_sum = predictions_df.groupby([\"best_class\",\"second_best_class\"])[\"second_highest_prob\"].sum().reset_index()\n\n                        def update_predictions(predict_number, y):\n                            for i in range(predict_number):\n                                print(\"=========\", i, \"==========\")\n                                X_scaled_sample = X_scaled.sample(50000)\n                                y_sample = y.loc[X_scaled_sample.index]\n                                bgmmC = BayesianGMMClassifier(\n                                                                tol = 1e-3,\n                                                                n_init = 3,\n                                                                max_iter = 300,\n                                                                random_state = i,\n                                                                n_components = 7,\n                                                                init_params = 'kmeans',\n                                                                covariance_type = 'full'\n                                                             )\n                                bgmmC.fit(X_scaled_sample[BEST_COLS], y_sample)\n                                pred_probs = bgmmC.predict_proba(X_scaled[BEST_COLS])\n                                pred_probs = pd.DataFrame(pred_probs, columns = range(1, 8))\n                                score_clusters(X_scaled[BEST_COLS], pred_probs.idxmax(axis = 1), verbose = True)\n                                y = pred_probs.idxmax(axis = 1)\n                            return pred_probs\n\n                        predicted_probabilities = update_predictions(predict_number = 75, y = predictions_df[\"best_class\"])\n                        predictions_df = best_class(predicted_probabilities)\n                        second_highest_probs_sum = predictions_df.groupby([\"best_class\",\"second_best_class\"])[\"second_highest_prob\"].sum().reset_index()\n                        submission = pd.read_csv(\"../input/tabular-playground-series-jul-2022/sample_submission.csv\")\n                        submission[\"Predicted\"] = predictions_df[\"best_class\"]\n                        submission.to_csv('submission_%s_%s_%s_%s_%s.csv' % (w_1, w_2, w_3, w_4, w_5), index = False)\n\nelse:\n\n    submission = pd.read_csv(\"../input/tps-jul-final-search/submission_ens.csv\")\n    submission.to_csv('submission.csv', index = False)","metadata":{"papermill":{"duration":0.031884,"end_time":"2022-07-17T03:10:58.833931","exception":false,"start_time":"2022-07-17T03:10:58.802047","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2022-07-26T09:02:44.130357Z","iopub.execute_input":"2022-07-26T09:02:44.130831Z","iopub.status.idle":"2022-07-26T09:02:44.313751Z","shell.execute_reply.started":"2022-07-26T09:02:44.130774Z","shell.execute_reply":"2022-07-26T09:02:44.312879Z"},"trusted":true},"execution_count":null,"outputs":[]}]}